# Stable bubble formations in a depth-perturbed Hele-Shaw channel

J. Lawless  Department of Physics and Astronomy, University of Manchester,
Oxford Road, Manchester, M13 9PL, UK  J. S. Keeler  School of Mathematics,
University of East Anglia, Norwich Research Park, Norwich, NR4 7TJ, UK  A. L.
Hazel  Department of Mathematics, University of Manchester, Oxford Road,
Manchester, M13 9PL, UK  A. Juel  Department of Physics and Astronomy,
University of Manchester, Oxford Road, Manchester, M13 9PL, UK
<EMAIL_ADDRESS>
###### Abstract

Deformable bubbles propagated by the flow of a viscous liquid in a Hele-Shaw
channel of uniform depth tend to travel along the channel’s axis. Pairs of
neighbouring bubbles will either separate or coalesce because an individual
bubble’s propagation speed increases monotonically with its size; this means
that any group of bubbles will eventually rearrange itself in order of
decreasing bubble size and all bubbles will separate. We show that, by
introducing a small geometric perturbation to the channel in the form of an
axially-uniform depth reduction along its centreline, the system supports a
multitude of stable bubble formations and this can disrupt the necessary
ordering by bubble size. The constituent bubbles of a stable formation lie in
alternation on opposite sides of the depth-perturbation, retain fixed shapes
and propagate steadily at the same speed. A stable formation is always led by
the smallest of its constituent bubbles, which would be the slowest in
isolation. The leading bubble propagates as if it were isolated and the
trailing bubbles reduce their speeds by adjusting their shapes and overlap of
the depth-perturbation in the perturbation fields of their preceding nearest-
neighbours in order to match that of the leading bubble. The trailing bubbles
can be arranged in any order and, hence, the number of stable formations
increases factorially as the number of bubbles is increased.

###### keywords:

##  1 Introduction

Dispersed two-phase flows are found in a variety of natural and industrial
processes. The dispersed phase may consist of gas bubbles or liquid droplets
that are suspended within a continuous liquid phase, examples of which include
oceanic flows (Deike, 2022), volcanic flows (Namiki et al., 2022), bubble
column reactors (Kantarci et al., 2005) and oil recovery processes (Xiao-Xuan,
2007). Due to their complexities, it is common for such flows to appear
disordered. However, dispersed flows have been observed to self-organise and
form coherent spatial structures in confined geometries. This behaviour has
been reported in a host of recent studies in microfluidics, in which
techniques such as flow-focusing are used to generate large quantities of
droplets or bubbles (Anna, 2016). For example, one-dimensional trains of
monodisperse microdroplets in a quasi-two-dimensional channel were observed to
exhibit unsteady collective dynamics in the form of longitudinal and
transverse vibrational modes when subjected to a finite amplitude perturbation
(Beatus et al., 2006; Schiller et al., 2015). This behaviour was found to be
driven by hydrodynamic interactions between neighbouring droplets, arising
from the dipolar flow field that each droplet imposes on its surroundings
(Beatus et al., 2012). A variety of geometrical features have also been found
to promote self-organisation in the flow of microdroplets. Upon encountering a
channel expansion, one-dimensional double emulsion droplet trains formed
regular droplet clusters through coalescence of their coating films
(Schirrmann et al., 2021). In a similar study, clusters of microdroplets were
found to spontaneously rearrange and form a variety of complex ‘building
block’-like shapes due to a coupling between inter-droplet hydrodynamic
interactions and adhesive depletion forces (Shen et al., 2016). In both cases,
variation of the flow parameters allowed for robust control of the generated
structures. Furthermore, an axial gradient in the channel depth profile was
found to promote the self-organisation of microdroplets into ordered two and
three-dimensional arrays (Parthiban et al., 2019). External actuation has also
been used to promote self-organisation. For example, groups of microbubbles
were observed to form periodic crystal-like lattices when excited by an
acoustic field (Rabaud et al., 2011). In this case, inter-bubble interactions
are mediated by acoustic surface waves emitted by the bubbles. We note that in
each of these examples, the bubbles or droplets exhibit little to no
deformation because they are sufficiently small such that surface tension has
the dominant influence on the interface shapes. However, in this paper, we
focus on the self-organisation of larger bubbles, for which the restoring
effects of surface tension cannot entirely counteract the deformation that is
induced by the surrounding fluid flow.

An idealised confined geometry to study multi-bubble dynamics is that of a
Hele-Shaw channel. A Hele-Shaw channel consists of two parallel plates that
are separated by a small gap relative to their width. In this particular
geometry, the interface that separates the two phases has a propensity to
deform if the bubbles are sufficiently large. Interfacial deformation is
promoted by viscous forces and resisted by surface tension; the capillary
number $\mathrm{Ca}$ quantifies the ratio between these two competing effects,
and is, therefore, the main dynamic parameter of interest. In general, pairs
of bubbles in a Hele-Shaw channel of uniform depth either separate or
coalesce. The terminal rising speed of a free-rising bubble in an inclined
Hele–Shaw channel increases monotonically with its size at low values of the
Reynolds number $\mathrm{Re}$ and, hence, a pair of co-axial bubbles that is
led by the smallest bubble will typically aggregate. However, in the reverse
configuration, the long-term outcome of the bubbles is determined by their
initial separation; well-separated bubbles tend to separate indefinitely
whereas aggregation may occur if their initial separation is sufficiently
small. The combination of these two dynamical behaviours leads to the
formation of steadily propagating compound bubble clusters, similar to the one
shown in the final panel of Figure 1(b), that separate indefinitely
(Maxworthy, 1986). However, the constituent bubbles of a cluster ultimately
coalesce due to the slow drainage of the lubrication films that separate them
(Chan et al., 2011). Although similar long-term dynamics are observed in the
high-$\mathrm{Re}$ regime, the transient dynamics are considerably richer
because free-rising bubbles generate unsteady wake flows. The core of the wake
generated by the leading bubble vertically entrains the trailing bubble and
regularly sheds counter-rotating vortices that lead to the periodic ejection
and subsequent realignment of the trailing bubble with the core of the wake as
it rises (Huisman et al., 2012; Filella et al., 2015, 2020). Broadly similar
dynamics were reported in generalised studies of homogeneous swarms of free-
rising bubbles dispersed throughout a channel (Bouche et al., 2012; Ruiz-Rus
et al., 2022).

Figure 1: (a, b) The dynamics of a pair of air bubbles in a horizontal Hele-
Shaw channel of uniform depth when subjected to a constant volume flux flow
$Q$ of a surrounding viscous liquid. Time increases from left-to-right. (a)
The larger bubble is initially situated in the leading position and the two
bubbles separate indefinitely. (b) The smaller bubble is initially situated in
the leading position and the two bubbles aggregate. These experiments were
performed in the setup that is described in section 2.1.1 without a rail.

When driven by the constant flux flow of a surrounding viscous fluid in a
horizontal Hele-Shaw channel, pairs of bubbles exhibit qualitatively similar
long-term behaviours to those of free-rising bubbles. The propagation speed of
an isolated bubble increases monotonically with its size and, hence, pairs of
bubbles of different sizes either separate indefinitely or aggregate, see
representative experimental time-sequences in figure 1. However, alternative
forms of collective long-term behaviours have been uncovered in analytical
studies of multiple bubbles in a Hele-Shaw channel. In the limit of infinite
surface tension, i.e. in the limit of vanishing $\mathrm{Ca}$, one-dimensional
trains of small, circular bubbles were found to exhibit longitudinal
vibrational modes when driven by the flow of a surrounding viscous liquid
(Booth et al., 2023). Interestingly, such modes were found to give rise to
collective dynamics likened to that of a Newton’s cradle through the
aggregation and subsequent breakup of neighbouring bubble pairs. In the
absence of surface tension, i.e. in the limit of infinite $\mathrm{Ca}$,
groups of deformable bubbles were found to propagate steadily with fixed
shapes in a variety of spatial formations when driven by the flow of a viscous
liquid (Vasconcelos, 2015). The bubbles in a formation neither separate nor
aggregate but, instead, maintain a constant separation from one another.
However, because this study did not include a stability analysis, it is
unknown if these formations could be observed in a real-world channel.

Figure 2: (a) A schematic diagram of the Hele-Shaw channel that contains a
small depth reduction or ‘rail’ along its centreline. Bubbles are generated at
the air port and are then guided to an asymmetric initial position by
constricting the channel with a steel washer. Oil is injected into the channel
at a constant volumetric flux $Q^{*}$ at the injection port. (b) Cross-
sectional schematic of the channel, where $W^{*}\gg H^{*}$. (C) A schematic
diagram of the procedure that is used to initialise a bubble on a particular
side of the rail.

In this paper, we perform experiments at finite $\mathrm{Ca}$ and demonstrate
how groups of bubbles of different sizes are able to self-organise into a
multitude of stable formations in a Hele-Shaw channel that possesses a small
non-uniformity in its geometry. The result is particularly striking because,
to the best of our knowledge, the only such states that exist in uniform
channels are the trivial states consisting of bubbles that travel at the same
speeds in isolation spaced far enough apart such that there are no
interactions between them. Throughout this paper, the geometric non-uniformity
is a small axially-uniform depth reduction along the centreline of the
channel’s cross-section, which is henceforth referred to as a rail. The
Reynolds number remains small throughout our experiments and, hence, inertial
forces are negligible, but the system is nonlinear as a consequence of the
presence of air-fluid interfaces. The behaviour of isolated bubbles in the
system with a rail has been studied extensively both experimentally (Gaillard
et al., 2021; Lawless et al., 2022) and numerically (Franco-Gómez et al.,
2018; Keeler et al., 2019), revealing a rich dynamical landscape that arises
as a direct consequence of the variation in channel depth. For sufficiently
large bubbles in an equivalent uniform channel, only a single solution branch
is found to be stable at all values of the driving flow rate (Tanveer &
Saffman, 1987). This solution branch corresponds to the family of steady,
symmetric, single-tipped bubbles that are observed in experiments. However, in
the system with a rail, additional solution branches are stabilised due to the
prescribed variation in the channel’s depth. Notably, a steady asymmetric
branch is found to be stable for all $\mathrm{Ca}$ (Keeler et al., 2019;
Gaillard et al., 2021). Furthermore, in contrast to a Hele-Shaw channel of
uniform depth, bubble breakup occurs more readily in the system with a rail,
which motivated Keeler et al. (2022) to study the dynamics of two interacting
bubbles. A family of stable two-bubble formations, characterised by the
bubbles propagating steadily on opposite sides of the rail, were uncovered in
their study through a combination of experiments and numerical simulations of
a two-dimensional depth-averaged lubrication model. However, the detailed
behaviour of these two-bubble states and the implications for steadily
propagating states with more than two bubbles were not investigated.

In this paper, we will establish the conditions for which stable two-bubble
formations occur by identifying the key physical mechanisms at play. We will
proceed to develop our understanding of stable two-bubble formations in order
to establish a set of “design rules” for three-bubble formations and examine
their behaviours in a response to changes in parameters. We then generalise
our results for formations with larger numbers of bubbles. The paper is
organised as follows. The experimental set-up and protocols implemented in
order to reproducibly generate bubbles of prescribed sizes are described in
section 2.1.1. Explaining the origin of the multi-bubble states relies on a
thorough understanding of the behaviour of individual bubbles and, so, we
benchmark the propagation of isolated bubbles in section 2.1.2. An outline of
the two-dimensional depth-averaged lubrication model and its numerical
implementation is provided in section 2.2. We investigate the range of
existence of stable two-bubble formations and their behaviours in response to
varying the bubble sizes, flow rate and rail height in section 3.1. The design
rules governing the existence of larger stable bubble formations are
established in section 3.2. We proceed to discuss the long-term behaviour of
an arbitrary train of bubbles in section 3.4 and conclusions are drawn in
section 4.

##  2 Methods

###  2.1 Experimental methods

####  2.1.1 Experimental set-up and protocols

Figure 3: (a) Variation of $U_{b}$ with $Q$ for isolated bubbles of various
sizes. Inset: experimental snapshots of the bubbles at $Q=0.01$. (b) Variation
of $U_{b}$ with $r$ at three values of $Q$ that span the range of flow rates
investigated. Inset: superimposed experimental snapshots of saturated bubbles
($r\geq r_{s}$) at $Q=0.01$, each of which are aligned in the streamwise
direction at their tips. The shape of the bubble’s front and rear are
unchanged for $r\geq r_{s}$. The rail height is $h=0.024$.

The experiments in this study were performed in the same Hele-Shaw channel
described by Gaillard et al. (2021), who provide detailed descriptions of both
the experimental set-up and protocols implemented to reproducibly generate and
propagate bubbles of prescribed sizes. Here, we recall the salient details.
Schematic diagrams of the set-up and channel cross-section are shown in figure
2. The channel consisted of two horizontally levelled, rectangular float glass
plates that were separated by two parallel strips of stainless steel shim. The
shim’s thickness $H^{*}=$ $1.00\pm 0.01\text{\,}\mathrm{mm}$ was measured at
uniformly spaced points along the channel’s length with a micrometer screw
gauge and they were separated by a distance $W^{*}=$ $40.0\pm
0.1\text{\,}\mathrm{mm}$ with a gauge block. The rail was formed by applying a
thin strip of translucent PET tape along the centreline of the lower glass
plate. The width of the tape was $w^{*}=$ $10.0\pm 0.1\text{\,}\mathrm{mm}$.
The experiments that are described in this paper were performed with two
different rail thicknesses $h^{*}$. For the majority of the experiments,
$h^{*}=$ $24\pm 1.\text{\,}\mathrm{\SIUnitSymbolMicro m}$, whilst a small
number of experiments were carried out at a reduced rail thickness $h^{*}=$
$10\pm 1.\text{\,}\mathrm{\SIUnitSymbolMicro m}$. Hence, the effective height
of the channel in the region above the rail was reduced either by $2.4\%$ or
$1\%$. The channel was filled with silicone oil (Basildon Chemicals Ltd) of
dynamic viscosity $\mu=$ $0.019\text{\,}\mathrm{Pa}\text{\,}\mathrm{s}$,
density $\rho=$ $951\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ and
surface tension $\sigma=$ $21\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}$
at the ambient laboratory temperature of $21\pm
1.\text{\,}\mathrm{\SIUnitSymbolCelsius}$. The flow of silicone oil was
controlled by a network of syringe pumps and solenoid valves. A three-way
solenoid valve was used to connect the injection syringes, channel inlet and
oil reservoir, whereas a two-way solenoid valve was used to connect the
channel outlet and oil reservoir. By changing the configuration of the valves,
oil could either be driven through the channel at a constant flow rate $Q^{*}$
or withdrawn from the reservoir in order to refill the injection syringes. The
syringes and valves were controlled by a custom LabVIEW script.

Air bubbles were generated at the air injection port situated a short distance
downstream of the channel inlet. In order to produce a single bubble, a
prescribed volume of oil was withdrawn slowly through a syringe connected to
the channel outlet whilst the air valve was open to the atmosphere. Once
generated, the bubble was detached and propagated away from the air valve by
imposing a small flow of oil. This process was repeated several times in
succession in order to generate multiple bubbles, each with a prescribed
volume. The initial separation between the bubbles was controlled by varying
the amount of oil that was injected following detachment. The bubbles were
then propagated through a ‘centring device’, which consisted of a symmetric
channel constriction and a linear expansion region. Upon exiting the centring
device, each of the bubbles were guided to an asymmetric initial position by
applying an asymmetric constriction to the channel with a magnetic washer. The
position of the washer was controlled non-invasively by an N52-grade neodymium
magnet situated beneath the lower glass plate and, by varying the side on
which the constriction was applied, bubbles could be initialised
asymmetrically on either side of the rail. A schematic of this procedure is
provided in figure 2(c).

Figure 4: (a) A schematic diagram of the truncated dimensionless computational
domain. The streamwise component of the bubbles’ centre of mass is constrained
to remain fixed at zero. (b) A cross-sectional schematic diagram of the
computational domain. The variable channel depth is modelled by a smoothed
depth profile $b(y)$, whose form is provided in equation 1.

The bubbles were set into motion by imposing a constant flow rate $Q^{*}$ and
filmed in top-view by an overhead CMOS camera, which was mounted onto a
motorised translation stage. The camera was programmed to translate at a
constant speed, which was selected via an empirical relationship between
$U_{b}^{*}$ and $Q^{*}$ to ensure that the bubbles remained within its field
of view. The observation window of the channel was $1760\crossproduct 330$
pixels and, depending on the value of $Q^{*}$, the camera’s frame rate was
varied between 20 and 60. The channel was uniformly illuminated from below by
a custom-built LED light box, which led to the air-fluid interface appearing
darkened due to light refraction. We used a Canny edge-detection algorithm in
order to identify the bubble contours.

We adopt the channel’s half-width $W^{*}/2$ and the mean speed of oil in an
equivalent uniform channel $U_{0}^{*}=Q^{*}/W^{*}H^{*}$ as our characteristic
length and velocity scales, respectively. We define a dimensionless flow rate
$Q=\mu U_{0}^{*}/\sigma$, which is a capillary number based on a
representative speed of the oil. The dimensionless bubble speed is
$U_{b}=U_{b}^{*}/U_{0}^{*}$, where the dimensional speed of a bubble
$U_{b}^{*}$ was determined by calculating the streamwise displacement of its
centre of mass across a series of consecutive frames. We note that in related
literature, the capillary number $\mathrm{Ca}$ is usually based on the speed
of the bubble, $\mathrm{Ca}=QU_{b}$, and, in an attempt to avoid confusion,
rather than redefining $\mathrm{Ca}$ we will use $Q$ as the principal dynamic
parameter throughout this paper. We parameterise bubble sizes in terms of a
dimensionless radius $r=2r^{*}/W^{*}$, where $r^{*}=\sqrt{A^{*}/\pi}$ and
$A^{*}$ is the dimensional in-plane area of the bubble.

####  2.1.2 Isolated bubbles

Throughout this paper, we will primarily work with four bubbles of size
$r=0.33$, $r=0.40$, $r=0.47$ and $r=0.54$ and they will be represented by
blue, red, green and magenta-coloured contours, respectively. We will first
benchmark the behaviour of each bubble when propagated in isolation. Here,
results are shown for the specific asymmetric states that combine to form the
multi-bubble states of interest. Although general information about these and
other isolated-bubble states is given in (Gaillard et al., 2021; Keeler et
al., 2022), the data presented here are new. Figure 3(a) shows the speed of
each isolated bubble as a function of flow rate and, in each case, the
bubble’s speed increases sharply from the quasistatic limit and saturates for
larger flow rates. We note that this behaviour is broadly similar to that of a
semi-infinite air finger, which is the limiting case as the bubble radius
tends to infinity. Figure 3(b) shows the speed of an isolated bubble as a
function of its size for several values of $Q$. The bubble’s speed increases
monotonically with size from the tracer-like limit, in which it moves at the
fluid speed, and saturates at $r=r_{s}\approx 0.90$. As highlighted in the
inset of figure 3(b), the bubble elongates in the streamwise direction whilst
the in-plane curvature of its front and rear remain unchanged for $r\geq
r_{s}$. The bubble speed is determined by the local pressure gradient, which
is affected by changes in the viscous dissipation that arise as the
propagating bubble displaces the viscous fluid, as well as the surface-
tension-induced pressure jump at the interface. The details depend on both the
shape and location of the bubble within the non-uniform channel and we are not
aware of any simpler predictive models than numerical solutions of the depth-
averaged system described in the next subsection.

###  2.2 Mathematical model and numerical methods

Figure 5: (a) Variation of $U$ with $y$ in a channel without a bubble at
$Q=0.04$ in the depth-averaged lubrication model. (b) The flow field in the
vicinity of an isolated bubble of size $r=0.33$ in the depth-averaged
lubrication model at $Q=0.04$. The dimensionless fluid speed $U$ is
represented by the heat map. (c) Variation of $U$ with the dimensionless
streamwise distance $D$ from the bubble’s tip.

We complement experiments with numerical simulations of a two-dimensional
depth-averaged lubrication model. The validity of this model has been
confirmed in studies of isolated bubbles (Gaillard et al., 2021) and multiple
bubbles (Keeler et al., 2022). We opt to describe only the model’s salient
details because its general implementation is unchanged.

####  2.2.1 Model outline

The truncated dimensionless computational domain and its cross-section are
shown in figure 4. The rail is modelled by a piecewise-smooth depth profile
$b(y)$ of the form

$b(y)=1-\frac{h}{2}\bigg{[}\tanh{\big{(}s(y+w)\big{)}}-\tanh{\big{(}s(y-w)\big{)}}\bigg{]},\vskip
6.0pt plus 2.0pt minus 2.0pt$ (1)

where $h=0.024$ is the dimensionless height and $w=0.25$ is the dimensionless
width of the rail. The parameter $s=40$ determines the sharpness of the rail’s
edges (Thompson et al., 2014).

We perform simulations in a co-moving frame of reference that translates with
velocity $\textit{{U}}(t)=(U_{b}(t),\>0)$, where $U_{b}(t)$ is an unknown of
the problem that is obtained by requiring that the streamwise component of the
centre of mass of the group of bubbles is fixed at zero. Upon applying
lubrication theory (Reynolds, 1886) and depth-averaging, the three-dimensional
Navier–Stokes equations in the fluid domain reduce to a single two-dimensional
lubrication equation in terms of the dimensionless fluid pressure $p$,

$\gradient\cdot\big{[}(b(y)^{3}\>\gradient{p}\big{]}=0.\vskip 6.0pt plus 2.0pt
minus 2.0pt$ (2)

We impose no-penetration conditions on the channel side walls, i.e. $p_{y}=0$
at $y=\pm 1$. We denote the coordinates of a point on the air-fluid interface
by R and impose a no-penetration condition:

$\frac{\partial\textbf{R}}{\partial
t}\cdot\hat{\textbf{n}}=\big{[}-b(y)^{2}\>\gradient{p}-\textit{{U}}(t)\big{]}\cdot\hat{\textbf{n}},\vskip
6.0pt plus 2.0pt minus 2.0pt$ (3)

where $\hat{\textbf{n}}$ is the outward-pointing unit normal vector to the
interface at R. The pressure jump across the interface, arising exclusively
due to surface tension, is described by the Young–Laplace equation:

$\llbracket
p\rrbracket^{\textnormal{bubble}}_{\textnormal{fluid}}=\frac{1}{3\alpha
Q}\bigg{[}\frac{\kappa}{\alpha}+\frac{1}{b(y)}\bigg{]},\vskip 6.0pt plus 2.0pt
minus 2.0pt$ (4)

where the dimensionless in-plane curvature of the interface at R is denoted by
$\kappa$. The transverse curvature term $1/b(y)$ arises under two simplifying
assumptions: (i) the bubble occupies the entirety of the channel depth and
(ii) the interface is semi-circular in the transverse direction, with radius
$b(y)/2$. The (constant) internal pressure of each bubble is determined by
enforcing that their dimensionless volumes remain constant. We note that thin
films, which are deposited on the upper and lower channel boundaries in
experiments, are not incorporated into the mathematical model. However, as in
previous studies (Gaillard et al., 2021; Lawless et al., 2022) we find that
the model captures all of the dynamics that are observed in experiments. There
is a quantitative discrepancy in the values of $Q$ at which particular
phenomena occur between experiments and numerical simulations. This
discrepancy increases as $Q$ is increased, which we believe is due to the
corresponding increases in the thicknesses of the neglected fluid films. We
will discuss this issue in section 3.1.

The system is spatially discretised by using a finite element method in the
open-source library oomph-lib (Heil & Hazel, 2006). We propagate elliptical
bubbles from specified initial placements in the computational domain by
imposing a constant pressure gradient in the streamwise direction. The value
of the pressure gradient is chosen in order to obtain a specified
dimensionless flow rate $Q$ of the surrounding fluid. The bubbles are
initialised sufficiently far apart in order to prevent their immediate
aggregation following the imposition of flow. The temporal evolution of the
bubbles is simulated by employing a backwards difference time-stepping
algorithm (BDF2) with a time-step $\Delta t=0.01$. We choose to omit details
relating to the handling of the finite element mesh and convergence tests, but
confirm that all the presented numerical results are resolved; details of
convergence tests are given in Gaillard et al. (2021).

####  2.2.2 Isolated bubbles

The depth-averaged lubrication model is used to quantify the flow field
perturbations that are imposed by the asymmetric propagating bubbles. The
single-phase flow field $U=U(y)$ is shown in figure 5(a); $U=0.97$ in the
occluded region and $U=1.02$ in the unoccluded region sufficiently far away
from the rail’s edges. The decrease in the fluid’s speed in the occluded
region is a consequence of the local increase in the channel’s viscous
resistance.

The fluid’s speed in the vicinity of an isolated bubble of size $r=0.33$ is
represented by the heat map in figure 5(b), where blue-coloured regions
correspond to $U<1$ and red-coloured regions correspond to $U>1$. The bubble’s
speed is more than twice that of the mean flow and a localised perturbation is
imposed onto the surrounding flow field. The fluid’s speed increases in the
vicinity of its front and rear, whereas it decreases in the vicinity of its
sides. The sharp decay of the flow field perturbation with increasing distance
is evident in figure 5(b), where $U$ is plotted as a function of the
dimensionless streamwise distance $D$ from the bubble’s tip. The fluid’s speed
at the bubble’s tip is equal to $U_{b}$ as a consequence of equation 3 and the
ambient fluid speed is reattained (within $0.5\%$) at $D=2$. The decay rate of
the perturbation is slower than that predicted by simple dipole models used to
model small bubbles in Hele-Shaw channels (Green & Vasconcelos, 2014),
indicating that these bubbles are in a regime in which finite-size effects are
important.

##  3 Results

Figure 6: (a, b) Experimental time-sequences of two bubbles of size $r=0.33$
and $r=0.40$ when propagated from rest at $Q=0.01$. The smaller bubble is
initially situated in the leading position in (a), whereas the larger bubble
is initially situated in the leading position in (b). (c) The corresponding
time-evolutions of $D$. (d - e) Superposition of the final bubble shapes with
their corresponding isolated bubbles (dashed black-coloured contours).

###  3.1 Stable two-bubble formations

We will begin by characterising the fundamental features of the stable two-
bubble formations that were identified by Keeler et al. (2022). The
experimental time-sequence in figure 6(a) is an example of two bubbles
organising into a stable formation when propagated from rest on opposite sides
of the rail at $Q=0.01$; the leading bubble’s size is $r_{1}=0.33$ and the
trailing bubble’s size is $r_{2}=0.40$. Thus, in the absence of any
interaction, the trailing bubble would propagate faster than the leading
bubble — see §2.1.2. If the two bubbles were located on the same side of the
rail, or if there were no rail, the trailing bubble would catch up to the
leading bubble and form an aggregate. However, the barrier between bubbles
presented by the rail provides enough additional resistance such that they do
not aggregate. Instead, the bubbles preserve their initial order following the
imposition of flow and evolve transiently until their shapes and separation
remain unchanged with time. This behaviour is reflected in the temporal
evolution of the dimensionless streamwise separation $D$ between the trailing
bubble’s front and the leading bubble’s rear in figure 6(c). The separation
between the two bubbles in their stable formation is $D=1.00$ and their speed
is $U_{b}=1.89$, which is equal to the leading bubble’s isolated speed. Upon
reversing the order of the two bubbles in figure 6(b), because the trailing
bubble propagates slower than the leading bubble, they proceed to separate
indefinitely. We note that this behaviour is reminiscent of figure 1(a). The
gradient of $D$ decreases monotonically and tends towards a constant value
$\textnormal{d}D/\textnormal{d}t=0.05$ for $D>2$, which is equal to the
difference between the speeds of the corresponding steadily propagating
isolated bubbles. We superimpose the shapes of these isolated bubbles (dashed
black-coloured contours) onto the final panels in figure 6(d, e). The two
separating bubbles in figure 6(d) adopt identical shapes to their isolated
counterparts for $D>2$ because they are effectively non-interacting. The two
bubbles adopt different shapes in their stable two-bubble formation, as
highlighted in figure 6(e). The leading bubble’s front is visibly unchanged,
whereas its rear inclines modestly away from the channel’s centreline.
However, the trailing bubble changes its shape considerably as it broadens and
increases its overlap of the rail. The changes in shape and position of the
trailing bubble are accompanied by a reduction in speed to that of the leading
bubble, which is what allows the stable formation to exist.

Figure 7: (a) A stable two-bubble formation obtained from a numerical
simulation of the depth-averaged lubrication model at $Q=0.04$. The two
bubbles have sizes $r_{1}=0.33$ and $r_{2}=0.40$. (b) The formation’s twelve
least-stable eigenvalues. (c) A heat map of the dimensionless fluid speed $U$
in the vicinity of the bubbles.

We find qualitatively similar behaviour in numerical simulations of the depth-
averaged lubrication model. A numerically computed stable two-bubble formation
is shown in figure 7(a). We used the same bubble sizes as in the experiment
but had to increase the flow rate to $Q=0.04$ in order to obtain the same
separation between the bubbles as in the experiments. The increased flow rate
is required because bubbles propagate at slower speeds in the numerical
simulations compared to the experiments for fixed $Q$, which is a consequence
of the neglect of the fluid films. The formation’s twelve least-stable
eigenvalues $\phi_{i}$ that were obtained from a linear stability analysis are
plotted in figure 7(b); a combination of real, negative eigenvalues and
complex conjugate eigenvalues with negative real components indicate that the
formation corresponds to a stable spiral in the system’s phase space. The
numerical model captures all of the qualitative characteristics that were
identified in experiments. For example, the formation’s speed $U_{b}=2.16$ is
marginally $(<1\%)$ slower than the leading bubble’s isolated speed
$U_{b}=2.17$. We superimpose the corresponding isolated bubble shapes (dashed
black-coloured contours) onto the stable two-bubble formation in figure 7(a).
The leading bubble’s shape is visibly unchanged from that of its isolated
counterpart, whereas the trailing bubble broadens and increases its overlap of
the rail (albeit by a smaller amount than in experiments). The flow field
ahead of the leading bubble is unaffected by the flow field perturbation that
is imposed by the trailing bubble, as highlighted in figure 7(c). However, the
two bubbles form an adjoined flow field perturbation in the region between the
leading bubble’s rear and the trailing bubble’s front because the bubbles are
within the interaction range that was established in figure 5(c). We obtain a
similar qualitative agreement for a variety of different flow rates and bubble
sizes and, having thus established a broad agreement between experiments and
the mathematical model, we will now proceed to explore the effect of key
parameter variations on the behaviour of a stable two-bubble formation.

Figure 8: (a - c) Experimental snapshots of stable two-bubble formations as
the middle bubble’s size $r_{2}$ is increased. The leading bubble’s size
$r_{1}=0.33$ and the flow rate $Q=0.04$ are fixed. (d) Variation of $D$, the
streamwise separation between the leading bubble’s rear and the trailing
bubble’s front, with $r_{2}$. (e) Variation of $\Delta U_{b}$, the difference
between the two isolated bubble speeds, with $r_{2}$. Inset: log-log plot of
$D$ versus $\Delta U_{b}$. (f) Variation of $\lambda$, the maximum width of
the trailing bubble’s cross-section, and $\lambda_{\textnormal{rail}}$, the
trailing bubble’s maximum overlap of the rail, with $r_{2}$. Error bars have
been omitted in (e) and (f) because they are smaller than the data markers.

####  3.1.1 Variation of $D$ with $r_{2}$

The experimental snapshots in figure 8(a - c) show the stable two-bubble
formation as the size of the trailing bubble is increased whilst the leading
bubble’s size $r_{1}=0.33$ and flow rate $Q=0.04$ are fixed. Remarkably, the
leading bubble’s shape is hardly affected by the change in the trailing
bubble’s size; the shape of its front is retained, whereas its rear
progressively inclines away from the channel’s centreline as $r_{2}$ is
increased. Moreover, the leading bubble’s speed and, therefore, the speed of
the entire formation, $U_{b}=1.89$, is unchanged for all values of $r_{2}$.
However, as evidenced in figure 8(d), the separation between the bubbles
decreases as $r_{2}$ is increased. This is because as $r_{2}$ increases, the
trailing bubble’s isolated speed increases monotonically, as per figure 3(b).
Hence, the reduction in the trailing bubble’s speed $\Delta U_{b}$ required in
order for it to match the leading bubble’s speed also increases monotonically
with $r_{2}$. We plot $\Delta U_{b}$ as a function of $r_{2}$ in figure 8(e).
An increased speed reduction requires an increased perturbation to the
isolated bubble shape, which is achieved by the trailing bubble moving closer
to the leading bubble so that it experiences a local pressure perturbation of
increased magnitude. We establish that $D$ and $\Delta U_{b}$ are
approximately related by a power law, as indicated by the linear variation of
$D$ with $\Delta U_{b}$ on a log-log scale in the inset of figure 8(e).
Eventually the bubble separation, $D$, saturates at a minimum value
$D=D_{\textnormal{min}}\equiv 0.72$ for $r_{2}\geq r_{s}$ because the speed of
the trailing bubble reaches a maximum value (the limit of a semi-infinite air
finger) and so $\Delta U_{b}$ saturates at a maximum value $\Delta
U_{b}=\Delta{U_{b}}_{\textnormal{max}}=0.15$. The bubble separation increases
sharply in the limit $r_{2}\rightarrow r_{1}$ because $\Delta U_{b}\rightarrow
0$: when the bubbles are the same size, they will travel at the same speed.

Figure 9: (a - c) Experimental snapshots of a stable two-bubble formation as
the flow rate $Q=\mu U_{0}^{*}/\sigma$ is increased. The leading bubble’s size
$r_{1}=0.33$ and the trailing bubble’s size $r_{2}=0.40$ are fixed. (d)
Variation of $D$, the streamwise separation between the leading bubble’s rear
and the trailing bubble’s front, with $Q$. (e) Variation of $U_{b}$, the
bubble’s speed relative to the mean fluid speed, with $Q$. (f) Variation of
$\lambda$, the maximum width of the trailing bubble’s cross-section and
$\lambda_{\textnormal{rail}}$, the trailing bubble’s maximum overlap of the
rail, with $Q$. Error bars have been ommitted in (e) because they are smaller
than the data markers.

As mentioned above, the trailing bubble’s shape becomes increasingly perturbed
from the isolated state as its size increases because a relatively greater
reduction in speed is required to match the speed of the leading bubble. There
are two different mechanisms by which reconfiguration of the trailing bubble
can lead to a change in speed: the width of the bubble can change, or the
position of the bubble within the channel can change. Figure 8(f) shows that
the maximum width of the trailing bubble’s cross-section $\lambda$ and the
maximum width of its cross-section that overlaps the rail
$\lambda_{\textnormal{rail}}$ (red-coloured data points) deviate increasingly
from their isolated counterparts (blue-coloured data points) as $r_{2}$ is
increased, see the inset in figure 8(f) for definitions of $\lambda$ and
$\lambda_{\textnormal{rail}}$. The trailing bubble fully overlaps the rail at
$r_{2}=0.62$, i.e. when $\lambda_{\textnormal{rail}}=0.50$, and continues to
broaden thereafter; this behaviour indicates that the trailing bubble’s
reduction in speed is achieved by a complex interplay between the variation in
its width (shape) and overlap of the rail (position).

####  3.1.2 Variation of $D$ with $Q$

The experimental snapshots in figure 9(a - c) show the stable two-bubble
formation as the flow rate is increased whilst the bubble sizes $r_{1}=0.33$
and $r_{2}=0.40$ are fixed. Figure 9(d) shows that the bubble separation
decreases as $Q$ is increased from the quasistatic limit. The behaviour is
more complex than section 3.1.1 because the leading bubble’s speed and shape
vary significantly as $Q$ is increased. The formation’s speed is equal to the
leading bubble’s isolated speed at all values of $Q$, as indicated in figure
9(e). However, there is a negligible increase in $\Delta U_{b}$ as $Q$ is
increased because the speed of both bubbles increases by approximately the
same amount. Figure 9(f) shows that the trailing bubble’s width and overlap of
the rail both increase by an approximately constant value as $Q$ is increased
because the same reduction in speed requires a similar degree of geometric
remodelling. We, therefore, attribute the decrease in the bubble separation
with flow rate to the decreasing significance of the pressure perturbation
that is induced by the leading bubble relative to the background pressure
gradient. The trailing bubble needs to move closer to the leading bubble in
order to experience the required perturbation that results in the necessary
geometric change.

####  3.1.3 Two-bubble phase diagram

The phase diagram in figure 10(a) classifies the long-term behaviour of all
possible bubble pairs $(r_{1},\>r_{2})$ when propagated from rest on opposite
sides of the rail at a fixed value of the flow rate $Q=0.01$. We note that the
bubbles were initialised sufficiently far apart in order to prevent their
immediate aggregation following the imposition of flow. We find that their
dynamics can be grouped into three simply connected regions, denoted by I, II
and III, which are described as follows.

Figure 10: (a) A phase diagram that classifies the long-term behaviour of two
bubbles initialised on opposite sides of the rail when propagated from rest at
$Q=0.01$. The leading bubble has size $r_{1}$ and the trailing bubble has size
$r_{2}$. The solid blue markers correspond to the bubbles organising into a
stable formation, whereas the solid white markers correspond to an overtaking
of the leading bubble by the trailing bubble. (b - c) Experimental time-
sequences of two bubbles when propagated from rest at $Q=0.01$. The trailing
bubble’s size is fixed at $r_{2}=1.00$. The leading bubble’s size is
$r_{1}=0.30$ in (b) and $r_{1}=0.29$ in (c).

  1. I.

Within this region, where $r_{1}>r_{2}$, the trailing bubble’s isolated speed
is less than that of the leading bubble and, hence, the two bubbles separate
indefinitely. The two bubbles ultimately propagate independently because their
individual flow field perturbations become increasingly separated.

  2. II.

Within this region, where $r_{1}<r_{2}$, the two bubbles do not organise into
a stable two-bubble formation and, instead, the trailing bubble overtakes the
leading bubble. In this case, the flow field perturbation imposed by the
leading bubble is too small to provoke the geometric changes required in the
trailing bubble to reduce its speed to that of the leading bubble.

  3. III.

Within this region, where $r_{1}<r_{2}$, the two bubbles organise into a
stable two-bubble formation because the flow field perturbation imposed by the
leading bubble is able to remodel the trailing bubble sufficiently and reduce
its speed to that of the leading bubble. A stable two-bubble formation exists
for all $0.30\leq r_{1}<r_{2}$. However, a stable two-bubble formation only
exists within a limited interval of $r_{2}$ for $0.24\leq r_{1}\leq 0.29$ and
the size of the interval decreases as $r_{1}$ decreases. We were unable to
resolve any stable two-bubble formations for $r_{1}<0.24$.

The experimental time sequences in figure 10(b - c) show the transition
between regions II and III for $r_{2}=1.00$. The two bubbles organise into a
stable two-bubble formation in figure 10(b) for $r_{1}=0.30$, whereas the
trailing bubble overtakes the leading bubble in figure 10(c) for $r_{1}=0.29$.
The two bubbles aggregate and form a compound bubble, which propagates at a
speed equal to that of an isolated bubble of size $\sqrt{r_{1}^{2}+r_{2}^{2}}$
(Gaillard et al., 2021). We note that the compound bubble only exists
transiently because the lubrication film that separates the two bubbles drains
over time and this ultimately results in their coalescence, although a longer
channel would be required to observe this in the experiments.

In summary, we find that stable steadily propagating two-bubble states are
possible if the bubbles are located on opposite sides of the rail and the
leading bubble is smaller than the trailing bubble. However, the leading
bubble must also be sufficiently large such that its local perturbation to the
fluid pressure is sufficient to cause the required deformation of the trailing
bubble that reduces its speed to match that of the leading bubble.

###  3.2 Stable three-bubble formations

####  3.2.1 Variation of bubble order

Figure 11: (a - f) Experimental time-sequences of the six distinct
arrangements of three bubbles initialised in alternation on opposite sides of
the rail when propagated from rest at $Q=0.01$. The bubble sizes are $r=0.33$
(blue-coloured), $r=0.40$ (red-coloured) and $r=0.47$ (green-coloured). An
arrow between a pair of bubbles in the final panels indicates that they
separate indefinitely, whereas a dashed box around a pair of bubbles indicates
that they organise into a stable two-bubble formation.

Building on our understanding of stable two-bubble formations, we will proceed
to explore the stability of larger groups of bubbles. There are $3!=6$
distinct initial arrangements of three bubbles lying in alternation on
opposite sides of the rail and we index the bubbles in ascending order in the
direction of outlet-to-inlet. We retain the two bubble sizes that were used in
figure 6 and add a third bubble of size $r=0.47$. The experimental time-
sequences in figure 11(a - f) show the evolution of each arrangement when
propagated from rest at $Q=0.01$. The two arrangements in figure 11(a, b)
organise into stable three-bubble formations and, as for a stable two-bubble
formation, they are both led by the smallest bubble and their speeds are equal
to the leading bubble’s isolated speed. The four arrangements in figure 11(c -
f) do not organise into stable three-bubble formations and, instead, at least
one pair of neighbouring bubbles separate. For example, the three bubbles in
figure 11(c) are arranged in order of decreasing size and all separate in
accordance with the differences between their isolated speeds. However, two of
the bubbles in figure 11(d - f) organise into a stable two-bubble formation
and proceed to separate from the remaining bubble.

We have superimposed the shapes of the two smallest bubbles in their stable
two-bubble formation (dashed black-coloured contours) onto the stable three-
bubble formations in figure 12(a, b). The leading bubble retains its shape in
both of the formations. The order of the two smallest bubbles is preserved in
figure 12(a) and the shape of the middle bubble’s front and its separation
from the leading bubble’s rear are retained. However, its rear inclines away
from the channel’s centreline because of the trailing bubble’s imposed flow
field perturbation. We note that this behaviour is consistent with that of the
leading bubble in a stable two-bubble formation. The middle and trailing
bubbles are interchanged in figure 12(b) and the value of $D_{1-2}$ is equal
to that which is obtained for a stable two-bubble formation in figure 8(b).
The trailing bubble retains its shape from the two-bubble formation and this
confirms that a bubble’s shape and position in the channel are the predominant
factors that influence its speed. Hence, because each individual bubble’s
speed must remain the same in order for it to remain part of the formation,
the bubble shapes are only weakly influenced by their relative axial position
within a stable formation.

Figure 12: (a, b) Stable three-bubble formations in experiments at $Q=0.01$.
(c, d) Stable three-bubble formation in the depth-averaged lubrication model
at $Q=0.04$. The bubble sizes are $r=0.33$ (blue-coloured), $r=0.40$ (red-
coloured) and $r=0.47$ (green-coloured). The dashed black-coloured contours
are the stable two-bubble formation bubble shapes for the blue and red-
coloured bubbles. (e) Heat map of the fluid’s speed in the vicinity of the
bubbles in (c). Figure 13: (a - c) Experimental snapshots of various stable
three-bubble formations at $Q=0.01$. The bubble sizes are $r_{2}=0.40$,
$r_{3}=0.47$ and $r_{1}$ is varied. (d) Variation of $D_{1-2}$ and $D_{2-3}$
with $r_{1}$. (e) Log-log plot of $D$ versus $\Delta U_{b}$. (f) Variation of
the three bubble widths $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ with
$r_{1}$. (g) Experimental time-sequence showing the evolution of three bubbles
of size $r_{1}=0.26$, $r_{2}=0.40$ and $r_{3}=0.47$ when propagated from rest
at $Q=0.01$. The leading bubble is too small for a stable three-bubble
formation to exist. Error bars have been omitted in (d) and (e) because they
are smaller than the data markers.

We identified qualitatively similar behaviour in numerical simulations of the
depth-averaged lubrication model. The analogous plots to figure 12(a, b) are
shown in figure 12(c, d) at $Q=0.04$ for the same bubble sizes. The speed of
each formation is marginally ($<1\%$) slower than the leading bubble’s
isolated speed. We note that $D_{1-2}$ increases by $0.07$ in figure 12(c),
which contrasts with experiments because this separation remained unchanged.
However, this is the only discrepancy that was identified between experiments
and numerical simulations. The heat map in figure 12(e) confirms that the
fluid’s speed in the vicinity of the leading and middle bubbles is largely
unchanged relative to that in their stable two-bubble formation shown in
figure 7(c).

####  3.2.2 Variation of $D$ with $r_{1}$

We have determined that the leading bubble sets the speed of a stable three-
bubble formation and the two trailing bubbles adjust their shapes and overlap
of the rail in the perturbation fields of their preceding neighbours in order
to reduce their speeds to that of the leading bubble. We will now explore how
a stable three-bubble formation responds to a change in its speed, which is
imposed by varying the size of the leading bubble whilst keeping the size of
the other two bubbles and the flow rate fixed. The experimental snapshots in
figure 13(a - c) show the stable three-bubble formation as the leading
bubble’s size is decreased whilst $r_{2}=0.40$, $r_{3}=0.47$ and $Q=0.01$ are
fixed. The bubbles become increasingly deformed and, as evidenced in figure
13(d), their separations decrease. Figure 13(e) shows that $D_{2-3}$ decreases
linearly with $\Delta U_{b}$ on a log-log scale and this is consistent with
figure 8(e). However, this is not the case for $D_{1-2}$ and we attribute this
to the change in the leading bubble’s shape. The blue-coloured data points in
figure 13(f) show that the leading bubble’s width decreases significantly as
$r_{1}$ is decreased, which drives an additional decrease in $D_{1-2}$.

The experimental time-sequence in figure 13(g) shows that a stable three-
bubble formation does not exist if the leading bubble is too small. For
$r_{1}\leq 0.26$, the leading bubble is overtaken by the middle bubble
following the imposition of flow. The initial overtaking dynamics are
qualitatively identical to the experimental time-sequence in figure 10(c) and,
hence, we infer that they are uninfluenced by the trailing bubble.

Figure 14: (a) Variation of $\lambda$ and $\lambda_{\textnormal{rail}}$ with
$r$ for an isolated bubble at $h=0.024$ and $h=0.010$. (b) Variation of
$U_{b}$ with $r$ for an isolated bubble at $h=0.024$ and $h=0.010$. (c, d) A
stable three-bubble formation at (c) $h=0.024$ and (d) $h=0.010$ in
experiments. The flow rate is $Q=0.04$ and the bubble sizes are $r_{1}=0.33$,
$r_{2}=0.40$ and $r_{3}=0.47$. (e) Variation of $\mathrm{d}U_{b}/\mathrm{d}r$
with $r$ for an isolated bubble at $Q=0.04$; these curves were obtained by
differentiating numerical power-law fits to the data points in (b). Error bars
have been omitted in all figures because they are smaller than the data
markers.

####  3.2.3 Variation of $D$ with $h$

We will now investigate the role of the rail height on the behaviour of a
stable three-bubble formation. The rail height was reduced to $h=0.010$ in the
experimental channel by applying a thinner strip of PET tape. Figure 14(a)
shows that an isolated bubble broadens and increases its overlap of the rail
as the rail height is reduced. The remodelling of an isolated bubble’s shape
reduces its speed towards that of an equivalent-sized bubble in a channel of
uniform depth ($U_{b}\approx 2$), as evidenced in figure 14(b). We note that
the difference between an isolated bubble’s speed at the two rail heights
increases with its size because larger bubbles occupy a greater proportion of
the rail and, hence, they are increasingly affected by changes in its
geometry.

Figure 15: (a - c) Snapshots of stable three-bubble formations at various
values of $h$ in the depth-averaged lubrication model at $Q=0.04$. The bubble
sizes are $r_{1}=0.33$, $r_{2}=0.40$ and $r_{3}=0.47$. (d) Variation of
$D_{1-2}$ and $D_{2-3}$ with $h$. (e) Variation of $\Delta U_{b}/\Delta r$
with $r$ at various values of $h$. Inset: Variation of $U_{b}$ with $r$ at
various values of $h$. (f) Log-log plot of $D$ versus $\Delta U_{b}$.

The experimental snapshots in figure 14(c, d) show the stable three-bubble
formation that is formed by three bubbles of size $r_{1}=0.33$, $r_{2}=0.40$
and $r_{3}=0.47$ at the two rail heights. The bubbles increase in separation
as the rail height is reduced and this is attributed to the changes in their
isolated speeds. The two solid curves in figure 14(e) show that the rate of
change of $U_{b}$ with $r$ decreases as the rail height is reduced and, hence,
the isolated bubble speeds become more similar. We find that the two curves
are directly proportional and the dashed black-coloured curve corresponds to a
rescaling of the blue-coloured curve by the scale factor $K=0.62$. The value
of $\Delta U_{b}$ between any two bubbles is, hence, multiplied by $K=0.62$ as
the rail height is reduced and, because a smaller change in relative speed is
required, this leads to an increase in separation between the bubbles. We note
that $D_{1-2}$ and $D_{2-3}$ increase by approximately the same multiplicative
factor and this is presumably because the bubble sizes are similar.

We find qualitatively similar behaviour in numerical simulations of the depth-
averaged lubrication model. We retain the three bubble sizes that were used in
experiments and, as evidenced in figure 15(a - d), the bubbles become
increasingly separated as the rail height is reduced. We plot
$\textnormal{d}U_{b}/\textnormal{d}r$ as a function of $r$ for various values
of $h$ in figure 15(e) and, in qualitative agreement with experiments, the
curves are found to be directly proportional and the dashed black-coloured
curves correspond to a rescaling of the $h=0.024$ curve, each by a determined
scale factor $K<1$. The scale factor $K$ decreases as $h$ is reduced and,
hence, the isolated bubble speeds become increasingly similar. As highlighted
in figure 15(f), $D_{1-2}$ and $D_{2-3}$ decrease near-linearly with $\Delta
U_{b}$ on a log-log scale and this indicates that the reduction in their
separations is predominantly driven by their isolated speeds becoming
increasingly similar. The small deviation from a perfect linear relationship
is attributed to the accompanying changes in the bubble shapes.

Figure 16: (a, b) The six stable four-bubble formations that were identified
(a) in experiments at $Q=0.01$ and (b) in numerical simulations at $Q=0.04$.
The bubble sizes are $r=0.33$ (blue-coloured), $r=0.40$ (red-coloured),
$r=0.47$ (green-coloured) and $r=0.54$ (purple-coloured). The shapes of the
blue, red and green-coloured bubbles in their stable three-bubble formation
have been superimposed (dashed black-coloured contours) onto each panel. (c)
Heat map of the fluid’s speed in the vicinity of the bubbles in the top-left
panel of (b).

###  3.3 Larger bubble formations

####  3.3.1 Four-bubble formations

The main design rule for the stable two- and three-bubble formations is that
the leading bubble must be the smallest, and, therefore, the slowest. We find
the same rule for a group of four bubbles. There are $4!=24$ distinct initial
arrangements of four bubbles lying in alternation on opposite sides of the
rail. We retain the three bubble sizes that were used in figure 11(a - f) and
add a fourth bubble of size $r=0.54$. The six stable four-bubble formations
that were identified in experiments at $Q=0.01$ and in numerical simulations
at $Q=0.04$ are shown in figure 16(a, b). The six formations are all led by
the smallest bubble and their speeds are equal to the leading bubble’s
isolated speed in experiments and marginally ($<1\%$) slower than the leading
bubble’s isolated speed in numerical simulations. The bubble shapes are
approximately constant in all six formations, again confirming that a bubble’s
speed is predominantly determined by its shape and location within the cross-
section. The heat map in figure 16(c) confirms that the fluid’s speed in the
vicinity of the first three bubbles is largely unchanged relative to that in
their stable three-bubble formation in figure 12(e). The remaining eighteen
arrangements are not stable and partition into subgroups of independently
propagating bubbles which, depending on the order of the bubbles, are
comprised of stable three-bubble formations, stable two-bubble formations and
single bubbles.

Figure 17: (a, b) The process of building up to a stable six-bubble formation
in (a) experiments at $Q=0.01$ and (b) numerical simulations at $Q=0.04$. The
previous stable formation (dashed black-coloured contours) has been
superimposed onto each panel. The bubble sizes are $r=0.33$ (blue-coloured),
$r=0.40$ (red-coloured), $r=0.47$ (green-coloured), $r=0.54$ (purple-
coloured), $r=0.61$ (yellow-coloured) and $r=0.68$ (magenta-coloured).

####  3.3.2 Larger bubble formations

The observed behaviours in both experiments and numerical simulations suggest
that any particular bubble is primarily influenced by the flow field
perturbation that is imposed by the bubble situated immediately ahead of it.
The leading bubble propagates as if it were isolated because the flow field in
the vicinity of its front is undisturbed by the flow field perturbations
arising from the trailing bubbles. In order for there to be a stable steadily
propagating state, the trailing bubbles must be successively slowed down by a
‘domino effect’ of the mechanism that was described in section 3.1: the
trailing bubbles adjust their shapes and overlap of the rail in the
perturbation fields of their preceding neighbours in order to reduce their
speeds to that of the leading bubble. The local influence of each bubble’s
perturbation means that the established design rules will remain unchanged
when we consider the addition of more bubbles to the system and, hence, we
infer that such formations can be extended indefinitely. We expect that any
arrangement of bubbles lying in alternation on opposite sides of the rail will
organise into a stable formation, provided that they are led by the smallest
bubble. The smallest bubble’s size must also exceed the minimum bubble size
that can slow down the fastest ($r\geq r_{s}$) isolated bubble. We have
confirmed this up to six bubbles in both experiments and numerical simulations
and figure 17(a, b) are examples of the sequential process of building up to a
stable six-bubble formation from an isolated bubble.

###  3.4 Bubble train dynamics

Having established the design rules of stable bubble formations, we will
proceed to discuss the expected long-term behaviour of an arbitrary train of
initially well-separated bubbles. We first consider a train of $N$ bubbles
lying entirely in alternation on opposite sides of the rail, for which there
are two distinct behaviours.

  1. 1.

The bubbles will organise into a stable $N$-bubble formation if they are led
by the smallest bubble.

  2. 2.

The bubbles will partition into $k$ stable formations if they are not led by
the smallest bubble, where $2\leq k\leq N$. The stable formations are arranged
in order of decreasing speed in the direction of outlet-to-inlet and any two
consecutive stable formations will separate at a rate equal to the difference
between their leading bubbles’ isolated speeds.

The dynamics will be considerably more complex if the bubbles do not lie
entirely in alternation on opposite sides of the rail because neighbouring
bubbles on the same side of the rail are susceptible to aggregation events.
Namely, aggregation will occur if the leading bubble of a formation is faster
than the bubble ahead of it, or if the outermost-trailing bubble of a
formation is slower than the bubble behind it. Following aggregation, the
bubbles will reorganise in accordance with the established design rules. There
are several practical difficulties (e.g. an inadequate channel length) that
prevented a detailed experimental investigation of this behaviour and,
similarly, numerical simulations are constrained by the finite size of the
computational domain. However, based on previous results, we will proceed to
infer the dynamics that would occur post-aggregation. The schematics provided
in figure 18(a - d) act as visual aids in order to support the following
descriptions of this more general case.

Figure 18: (a - d) Schematic diagrams of the reorganisation dynamics that
could occur following an aggregation event involving either the leading or
outermost-trailing bubbles of a stable formation. Time increases from left-to-
right. (a) The size of the outermost-trailing bubble increases post-
aggregation. (b - d) The size of the leading bubble increases post-
aggregation. Dashed boxes indicate independently propagating formations.

  1. 1.

If the outermost-trailing bubble of a formation aggregates with the bubble
behind it, the formation will be retained. However, the separation between the
formation’s two outermost-trailing bubbles will decrease because the
outermost-trailing bubble’s value of $\Delta U_{b}$ will increase — see figure
18(a).

  2. 2.

If the leading bubble of a formation aggregates with the bubble ahead of it
and it remains the smallest of the formation’s bubbles post-aggregation, the
formation will be retained. However, the separation between each of the
bubbles will increase because each bubble’s value of $\Delta U_{b}$ will
decrease — see figure 18(b).

  3. 3.

If the leading bubble of a formation aggregates with the bubble ahead of it
and it does not remain the smallest of the formation’s bubbles post-
aggregation, the formation will be destroyed. The remaining bubbles will then
reorganise in accordance with the established design rules. The two most
extreme examples are those in which the remaining bubbles are arranged in
order of decreasing or increasing size. In the former case, they will
reorganise into a smaller formation and separate from the aggregated bubble —
see figure 18(c). In the latter case, all of the bubbles will separate — see
figure 18(d). In all other cases, the remaining bubbles will reorganise into a
number of smaller formations.

An arbitrary bubble train will exhibit a combination of all such behaviours
and a variety of stable bubble formations will transiently appear and
disappear during the system’s temporal evolution before it ultimately settles
in a particular long-term state comprised of various steadily propagating,
indefinitely separating subgroups.

##  4 Conclusion

We have investigated the behaviour of stable multi-bubble formations in a
Hele-Shaw channel that contains an axially uniform depth reduction along its
centreline. The presence of the depth reduction, or rail, allows individual
bubbles to propagate at different locations within the channel’s cross-section
(Franco-Gómez et al., 2018) via the existence of stable asymmetric states,
which can prevent the bubble aggregation that typically occurs in Hele-Shaw
channels of uniform depth. The emergence of the multi-bubble formations is an
example of complex self-organising dynamics in a confined two-phase flow and
we find that there can be multiple stable, steadily propagating formations
that exist for fixed system parameters. However, in contrast, bubbles in Hele-
Shaw channels of uniform depth will only form multi-bubble states in the
trivial, degenerate, case when all bubbles are the same size.

The constituent bubbles of a stable formation must lie in alternation on
opposite sides of the rail in order to prevent aggregation and they each
propagate at the same speed, advancing steadily with fixed shapes. We note
that Keeler et al. (2022) found two-bubble steady states in which the bubbles
were placed on the same side of the rail or aligned along the centre of the
channel in the depth-averaged model, but these were all unstable. A stable
bubble formation is always led by the smallest, which is the slowest, of its
constituent bubbles. The flow field ahead of the leading bubble is undisturbed
by the trailing bubbles and, hence, it propagates as if it were isolated. As a
consequence, the speed of a stable formation is set by the leading bubble’s
isolated speed. The presence of the rail allows the trailing bubbles to reduce
their speeds to that of the leading bubble by adjusting their shapes and
degrees of overlap of the rail within the perturbation fields of their
preceding neighbours. The extent of a bubble’s shape perturbation generally
increases as the difference between its isolated speed and that of the leading
bubble increases. However, if the leading bubble is too small ($r_{1}<0.24$),
then its perturbation field cannot change the shape of the trailing bubble
enough to match the speed of the leading bubble and a stable state is not
possible.

The speed of each bubble is set predominantly by its shape and relative
position within the channel’s cross-section, and each bubble’s local
environment is only influenced by its nearest downstream neighbour because the
spatial decay of the perturbation caused by each bubble is negligible within
one channel width. The required speed of each bubble is dictated by the speed
of the formation, which means that the trailing bubbles can be arranged in any
order and their shapes are largely unaffected when they are interchanged. As a
consequence, the number of possible stable arrangements for a group of bubbles
increases rapidly as the number of bubbles increases. The largest number of
bubbles that we investigated was six because it becomes increasingly difficult
to initialise larger groups of bubbles in experiments and numerical
simulations are constrained by the finite size of the computational domain.
However, our results suggest that stable formations will exist for arbitrarily
large groups of bubbles. We established an excellent qualitative agreement
between experiments and numerical simulations of a depth-averaged lubrication
model. However, a small degree of quantitative disagreement arises because the
model does not incorporate the fluid films that are deposited on the upper and
lower channel boundaries, which means that, in general, bubbles in the
simulations propagate more slowly than those in the experiments for the same
flow rate. Nevertheless, the broad agreement between our two approaches
indicates that the fluid films do not play a fundamental role in the
stabilisation of a formation.

The fact that each bubble is influenced only by its nearest downstream
neighbour means that we can understand the general behaviour of an arbitrary
train of bubbles of different sizes. The bubbles will undergo a variety of
complex reorganisations, including aggregation and merging events when
neighbouring bubbles are on the same side of the rail, until the train settles
into a final state consisting of steadily propagating formations of bubbles
that either propagate at the same speeds, maintaining fixed separations, or
separate indefinitely when the trailing formation propagates more slowly than
the leading one. The evolution towards the final state will contain
considerable structure, however, because the individual bubbles will organise
into local formations that then interact. Although the bubbles might be
expected to end up in approximate size order, the fact that the trailing
bubbles in an individual formation can be arranged in any order means that the
only definitive conclusion is that the smallest bubble will lead the most
downstream formation. In contrast, in a uniform channel, the smallest bubble
will always merge with any larger trailing bubbles until the final bubble in
the train is the smallest remaining. In other words, the presence of the rail
allows small bubbles to propagate as part of the overall train.

The geometric variation provided by the rail plays two crucial roles in this
system: (i) it enables simultaneous stable, steadily propagating states for
individual bubbles whose centres are not aligned along the channel’s axis;
(ii) it enables an extra mechanism by which bubble speeds can be changed
because the volume of fluid displaced, and therefore the local viscous
dissipation, can be varied without changing the bubble shape by altering the
proportion of the bubble that overlaps the rail. In the absence of the rail,
the speeds of the faster bubbles cannot be altered sufficiently to form
stable, steady multi-bubble states. As the height of the rail decreases, the
differences between speeds of the bubbles of different sizes decreases, the
bubble centres move closer to the centreline of the channel and the distance
between the bubbles increases in the steady formations. We expect that the
stability of the multi-bubble formations will be lost in the limit of a
vanishing rail height because the asymmetric, steadily propagating state is
lost for individual bubbles. Furthermore, we also expect that their stability
will be lost in the limit of the rail’s height approaching that of the channel
because, for sufficiently large rail heights, asymmetric bubbles no longer
overlap the rail and, instead, they propagate in disconnected side-channels on
opposite sides of the rail. The rail’s width is also expected to have a
similar influence on the stability of multi-bubble formations. As the rail’s
width is decreased, it will approach the limit of a wire and impose a
localised perturbation at the tip of a bubble similar to the experiments
described by Zocch et al. (1987). Conversely, in the limit that the rail’s
width approaches that of the channel, the bubbles will fully overlap it
propagate along the channel’s centreline. However, we have not pursed a
detailed investigation of these limits.

The transition from simple aggregation and separation events in a channel of
uniform depth to complex dynamics involving large numbers of steadily
propagating states containing multiple bubbles, simply by introducing a small
perturbation to the channel’s geometry, is remarkable. The presence of the
additional steady states will slow down or even remove a number of bubble
aggregation events and it may be, therefore, that such a small geometric
variation could help to stabilise bubble trains propagating over long
distances to slight variations in bubble size.

[Supplementary data]Supplementary material and movies are available at…

[Acknowledgements]Acknowledgements may be included at the end of the paper,
before the References section or any appendices. Several anonymous individuals
are thanked for contributions to these instructions.

[Funding]The authors acknowledge funding from the Engineering and Physical
Sciences Research Council through grant no. EP/P026044/1 and from the
Leverhulme Trust through grant no. ECF-2012-17.

[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…

[Author ORCIDs]
J. Lawless https://orcid.org/0000-0003-4074-7494;
J. S. Keeler https://orcid.org/0000-0002-8653-7970;
A. L. Hazel https://orcid.org/0000-0003-1130-0554;
A. Juel https://orcid.org/0000-0003-3342-7388.

[Author contributions]J. L. performed all of the experiments. J. S. K.
performed all of the simulations. J. L. processed and analysed all of the
data. A. J and A. L. H. conceived the scientific study. All of the authors
discussed the results and contributed to the writing and review of the final
manuscript.

## References

  * Anna (2016) Anna, S. L. 2016 Droplets and bubbles in microfluidic devices. Annu. Rev. Fluid. Mech. 48 (1), 285–309.
  * Beatus et al. (2012) Beatus, T., Bar-Ziv, R. & Tlusty, T. 2012 The physics of 2d microfluidic droplet ensembles. Phys. Rep. 516, 103–145.
  * Beatus et al. (2006) Beatus, T., Tlusty, T. & Bar-Ziv, R. 2006 Phonons in a one-dimensional microfluidic crystal. Nat. Phys. 2, 743–748.
  * Booth et al. (2023) Booth, D. J., Griffiths, I. M. & Howell, P. D. 2023 Circular bubbles in a Hele-Shaw channel: a Hele-Shaw Newton’s cradle. J. Fluid Mech. 954, A21.
  * Bouche et al. (2012) Bouche, E., Roig, V., Risso, F. & Billet, A-M. 2012 Homogeneous swarm of high-Reynolds-number bubbles rising within a thin gap. Part 1. Bubble dynamics. J. Fluid Mech. 704, 211–231.
  * Chan et al. (2011) Chan, D., Klaseboer, E. & Manica, R. 2011 Film drainage and coalescence between deformable drops and bubbles. Soft Matter 7, 2235–2264.
  * Deike (2022) Deike, L. 2022 Mass transfer at the ocean–atmosphere interface: The role of wave breaking, droplets, and bubbles. Annu. Rev. Fluid. Mech. 54 (1), 191–224.
  * Filella et al. (2015) Filella, A., Ern, P. & Roig, V. 2015 Oscillatory motion and wake of a bubble rising in a thin-gap cell. J. Fluid Mech. 778, 60–88.
  * Filella et al. (2020) Filella, A., Ern, P. & Roig, V. 2020 Interaction of two oscillating bubbles rising in a thin-gap cell: vertical entrainment and interaction with vortices. J. Fluid Mech. 888, A13.
  * Franco-Gómez et al. (2018) Franco-Gómez, A., Thompson, A. B., Hazel, A. L. & Juel, A. 2018 Bubble propagation in Hele-Shaw channels with centred constrictions. Fluid Dyn. Res. 50 (2), 021403.
  * Gaillard et al. (2021) Gaillard, A., Keeler, J., Le Lay, G., Lemoult, G., Thompson, A., Hazel, A. & Juel, A. 2021 The life and fate of a bubble in a geometrically perturbed Hele-Shaw channel. J. Fluid Mech. 914, A34.
  * Green & Vasconcelos (2014) Green, C. & Vasconcelos, G. L. 2014 Multiple steadily translating bubbles in a hele-shaw channel. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 470.
  * Heil & Hazel (2006) Heil, M. & Hazel, A. L. 2006 oomph-lib – an object-oriented multi-physics finite-element library. In Fluid-Structure Interaction (ed. H-J Bungartz & M. Schäfer), pp. 19–49. Springer Berlin Heidelberg.
  * Huisman et al. (2012) Huisman, S., Ern, P. & Roig, V. 2012 Interaction and coalescence of large bubbles rising in a thin gap. Phys. Rev. E 85 2 Pt 2, 027302.
  * Kantarci et al. (2005) Kantarci, N., Borak, F. & Ulgen, K. 2005 Bubble column reactors. Process Biochem. 40, 2263–2283.
  * Keeler et al. (2019) Keeler, J., Thompson, A., Lemoult, G., Hazel, A. & Juel, A. 2019 The influence of invariant solutions on the transient behaviour of an air bubble in a Hele-Shaw channel. Proc. R. Soc. 475.
  * Keeler et al. (2022) Keeler, J. S., Gaillard, A., Lawless, J., Thompson, A. B., Juel, A. & Hazel, A. L. 2022 The interaction of multiple bubbles in a Hele-Shaw channel. J. Fluid Mech. 946, A40.
  * Lawless et al. (2022) Lawless, J., Keeler, J., Gaillard, A., Hazel, A. L. & Juel, A. 2022 The unpredictable nature of bubble evolution. Sci. Rep. 12 (20752).
  * Maxworthy (1986) Maxworthy, T. 1986 Bubble formation, motion and interaction in a Hele-Shaw cell. J. Fluid Mech. 173, 95–114.
  * Namiki et al. (2022) Namiki, A., Lev, E., Birnbaum, J. & Baur, J. 2022 An experimental model of unconfined bubbly lava flows: Importance of localized bubble distribution. J. Geophys. Res. 127.
  * Parthiban et al. (2019) Parthiban, P., Doyle, P. S. & Hashimoto, M. 2019 Self-assembly of droplets in three-dimensional microchannels. Soft Matter 15, 4244–4254.
  * Rabaud et al. (2011) Rabaud, D., Thibault, P., Mathieu, M. & Marmottant, P. 2011 Acoustically bound microfluidic bubble crystals. Phys. Rev. Lett. 106, 134501.
  * Reynolds (1886) Reynolds, O. 1886 On the theory of lubrication and its application to Mr. Beauchamp tower’s experiments, including an experimental determination of the viscosity of olive oil. Philos. Trans. R. Soc. Lond. 177, 157–234.
  * Ruiz-Rus et al. (2022) Ruiz-Rus, J., Ern, P., Roig, V. & Martínez-Bazán, C. 2022 Coalescence of bubbles in a high Reynolds number confined swarm. J. Fluid Mech. 944, A13.
  * Schiller et al. (2015) Schiller, U. D., Fleury, J-B., Seemann, R. & Gompper, G. 2015 Collective waves in dense and confined microfluidic droplet arrays. Soft Matter 11, 5850–5861.
  * Schirrmann et al. (2021) Schirrmann, K., Caceres, G. & Juel, A. 2021 Self-assembly of coated microdroplets at the sudden expansion of a microchannel. Microfluid. Nanofluidics 25.
  * Shen et al. (2016) Shen, B., Ricouvier, J., Malloggi, F. & Tabeling, P. 2016 Designing colloidal molecules with microfluidics. Adv. Sci. 3 (6), 1600012.
  * Tanveer & Saffman (1987) Tanveer, S. & Saffman, P. G. 1987 Stability of bubbles in a Hele–Shaw cell. Phys. Fluids 30 (9), 2624–2635.
  * Thompson et al. (2014) Thompson, A. B., Hazel, A. L. & Juel, A. 2014 Multiple finger propagation modes in Hele-Shaw channels of variable depth. J. Fluid Mech. 746, 123–164.
  * Vasconcelos (2015) Vasconcelos, G. L. 2015 Multiple bubbles and fingers in a Hele-Shaw channel: complete set of steady solutions. J. Fluid Mech. 780, 299–326.
  * Xiao-Xuan (2007) Xiao-Xuan, X. 2007 Study on oil–water two-phase flow in horizontal pipelines. J. Pet. Sci. Eng. 59 (1), 43–58.
  * Zocch et al. (1987) Zocch, G., Shaw, B. E., Libchaber, A. & Kadanoff, L. P. 1987 Finger narrowing under local perturbations in the Saffman-Taylor problem. Phys. Rev. A 36, 1894–1900.