2.1. Experimental set-up

A schematic diagram of the experimental set-up is shown in figure 1. The Hele-Shaw channel used in the experiments consisted of two parallel float glass plates of dimensions $170\ \textrm {cm} \times 10\ \textrm {cm} \times 2\ \textrm {cm}$ separated by two strips of stainless steel shim. The shims were sprayed with a thin layer of adhesive (3M spray mount^TM) and bonded to the bottom glass plate with a distance $W^*$ between them, which was imposed using gauge blocks. The channel height was $H^* = 1.00 \pm 0.01$ mm and its width was $W^* = 40.0$ mm with variations up to $\pm 0.1$ mm over the length of the channel $L^* = 170$ cm, yielding an aspect ratio $\alpha = W^*/H^* = 40$. A small depth perturbation of rectangular cross-section, width $w^* = 10.0 \pm 0.1$ mm and thickness $h^* = 24 \pm 1\ \mathrm {\mu }$m, was positioned on the bottom glass plate, symmetrically about the channel centreline; see figure 1(b). This depth perturbation, henceforth referred to as a rail, was made from a translucent adhesive tape strip bonded to the bottom glass plate and its surface gently smoothed with sandpaper of grit size P1500. The strip width was larger than $w^*$ and the excess width on each side was trimmed using a sharp blade placed against a 15 mm gauge block sliding along the fixed steel shims. The top glass plate was then placed onto the shims and the channel was sealed with clamps and levelled horizontally to within 0.03$^{\circ }$.

The channel was filled with silicone oil (Basildon Chemicals Ltd) of dynamic viscosity $\mu = 0.019$ Pa s, density $\rho = 951$ kg m$^{-3}$ and surface tension $\sigma = 21$ mN m$^{-1}$ at the laboratory temperature of $21\pm 1\,^{\circ }$C. Three cylindrical brass pieces embedded into the top glass plate provided oil inlet and outlet ports as well as an air injection port (see figure 1c), and all connections were made with rigid plastic tubes. The output from three syringe pumps in parallel (Legato 200 series) was fed to a three-way solenoid valve connected to the channel inlet and to an oil reservoir at atmospheric pressure, which was also linked to the outlet of the channel via a two-way solenoid valve. Hence, depending on the setting of the three-way inlet valve, oil could be injected into the channel at a constant flow rate $Q^*$ (open outlet valve) or withdrawn from the reservoir to refill the syringes.

The air injection port was connected to ambient air through a two-way solenoid valve. An air bubble was injected into the channel through the air port by closing the outlet valve, opening the air valve and slowly withdrawing a set volume of oil from the channel via a syringe pump connected to the outlet. Once the bubble was formed, the air valve was closed, the outlet valve was opened, and a small amount of oil was slowly injected into the channel to nudge the bubble forward. This prevented the compression of the bubble into the air port when flow was initiated. A bubble centring device was positioned downstream of the bubble injection port, where the width of the channel was reduced to $W^*_i = 17$ mm over a length of 110 mm which expanded linearly over 40 mm into the channel of width $W^* = 40$ mm, thus featuring an expansion ratio of $W^*/W_i^* = 2.4$. The rail started within this device.

A centred bubble of controlled shape was obtained in the channel using a protocol described in § 3. This bubble was then set in motion from rest by imposing a constant oil flow rate $Q^*$. The propagating bubble was filmed in top view using a CMOS camera mounted on a motorised translation stage, which was placed above the channel and spanned its entire length $L^*$. An empirical relation between the bubble speed and oil flow rate was used to set a constant translation speed for the camera so that the propagating bubble remained within the field of view of the camera for the duration of the experiment. The channel was uniformly back lit with a custom made LED light box, so that the oil–air interface appeared as a black line through light refraction. The syringe pumps, valves, translation stage and camera were computer controlled via a LabVIEW code.

We measured the projected area $A^*$ of the bubble and its centroid position from the bubble contour, which was identified by an edge detection algorithm in terms of $x^*$ and $y^*$ coordinates spanning the length and width of the channel, respectively (see figure 1). We will examine the effects of flow on $A^*$ in § 2.2. We also determined the bubble velocity $U_b^*$ along the $x^*$ direction from the camera speed and the centroid position on successive frames. Henceforth, we will refer to the bubble size in terms of an equivalent radius $r^* = \sqrt {A^*/{\rm \pi} }$ based on the projected area of the bubble. We only considered bubbles larger than the rail, i.e. $2r^* > w^*$, for which the bubble confinement ratio $r^*/H^*$ was therefore larger than $5$. We denote as $y_c^*$ the distance from the bubble centroid to the centreline of the channel. We choose the channel half-width $W^*/2$ and the average oil velocity in a channel without a rail $U_0^* = Q^*/(W^* H^*)$ as the characteristic length and velocity scales in the ($x^*,y^*$) plane, respectively. The non-dimensional bubble size, centroid offset and velocity are therefore $r = 2 r^*/W^*$, $y_c = 2y^*_c/W^*$ and $U_b = U_b^*/U_0^*$, respectively. We define a non-dimensional flow rate $Q = \mu U_0^* / \sigma$ which is a capillary number based on the mean oil velocity (in the channel without a rail) and which ranges between $0.005$ and $0.19$ in the experiment. Note that $U_b$ measures the bubble velocity relative to that of the fluid far from the bubble, and need not vary monotonically with $Q$.

The capillary number based on the bubble velocity, which measures the ratio of viscous to capillary forces, is ${\textit {Ca}} = \mu U_b^* / \sigma$ with values in the range $0.008\le {\textit {Ca}} \le 0.48$. The Bond number ${\textit {Bo}} = \rho g H^{*2} / 4 \sigma =0.11$, where $g$ is the gravitational acceleration, indicates that bubble buoyancy is negligible because Jensen et al. (Reference Jensen, Libchaber, Pelcé and Zocchi1987) showed that for ${\textit {Bo}} < 1$, there is no change to the pressure jump across the fluid interface in the small-${\textit {Ca}}$ asymptotic limit of the Saffman–Taylor model. The ratio of inertial to viscous forces yields a reduced Reynolds number $\textit {Re} = \rho U_0^* H^{*2} / \mu W^*$ which takes values less than $0.3$ in our experiments, thus indicating broadly negligible inertia.

2.2. Effect of flow on the projected area of the bubble

The projected area of a propagating bubble is affected by two competing effects: films of oil separating the bubble from the top and bottom glass plates tend to increase the projected area of the bubble, while the pressure increase associated with the imposition of flow in the channel reduces the projected area through compression of the bubble.

Films of oil on the top and bottom glass plates encapsulate a newly created bubble because silicone oil wets the channel walls. When the bubble is at rest ($Q^*=0$), these films tend to be very thin, but they drain so slowly that they are always present on the time scale of the experiment. The thickness of these films increases with the capillary number for a propagating bubble (Bretherton Reference Bretherton1961), so that for a fixed bubble volume, the projected area of the bubble increases accordingly.

We estimated the effect of wetting films in our experiments by measuring their thickness for a bubble propagating in a channel where the rail was absent; see appendix A. We found that for the maximum value of ${\textit {Ca}}=0.48$ investigated in this paper, the film thickness averaged over the bubble area is $\langle h_f^* \rangle / H^* \approx 0.2$ (i.e. $\langle h_f^* \rangle \approx 200\ \mathrm {\mu }$m). This means that up to $40\,\%$ of the channel depth is filled with oil, and that the bottom film could be up to $8$ times thicker than the rail ($h^* = 24 \pm 1\ \mathrm {\mu }$m). We expect that the increasing film thickness will lead to a reduced influence of the rail on the bubble propagation as the flow rate increases.

Upon flow initiation, we systematically observed a net reduction of the bubble projected area despite thickening wetting films, which suggests that compression of the bubble due to the imposed flow plays a dominant role. This is illustrated in figure 2(a) in which a bubble initially at rest ($Q^*=0$) with a static projected area $A_s^* = 4.9$ cm$^2$ compresses rapidly when a flow rate $Q^* = 150$ ml min$^{-1}$ is imposed to reach a minimum area of $A_Q^*=4.3$ cm$^2$ (13 % reduction). This rapid decrease in bubble area is followed by a slow increase of approximately 2 % as the bubble propagates towards the channel outlet. We chose to propagate strongly asymmetric bubbles to measure compression effects because they persist over the entire range of flow rates investigated, whereas symmetric bubbles are prone to breakup when propagating from rest on the rail; see § 4. Such a bubble was prepared by initially propagating it at a low flow rate $Q^*=10$ ml min$^{-1}$ for which it systematically migrated sideways towards one of the deeper sides of the channel before interrupting the flow; see § 4.3.2 for a summary of the states of steady propagation.

Figure 2. Bubble compression in flow. (a) Variation of the projected area $A^*$ of a bubble initially at rest ($Q^*=0$), asymmetric about the channel centreline, when a flow rate $Q^* = 150$ ml min$^{-1}$ is imposed at $t^*=0$ as indicated by the arrow. The bubble area rapidly decreases from its static value $A_s^* = 4.9$ cm$^2$ to a minimum value $A_Q^*=4.3$ cm$^2$ (13 % reduction) due to the sudden pressure increase in the channel, and then slowly increases by approximately 2 % as it propagates towards the channel outlet where pressure is lower. Inset snapshots indicate the shape of the bubble before and after flow imposition and the grey band is the rail. (b) Compression ratio $A_Q^*/A_s^*$ as a function of the dimensional flow rate $Q^*$ and non-dimensional flow rate $Q$ for three static bubble sizes. The data are accurately captured by the linear relation $A_Q^*/A_s^* = 0.98 - K^* Q^*$ where $K^* = 4.1 \, 10^{4}$ s m$^{-3}$ (dash-dotted line). The dotted line corresponds to the prediction of the compression ratio when assuming ideal gas behaviour and neglecting wetting films (see main text).

Figure 2(b) shows the compression ratio $A_Q^*/A_s^*$ for different flow rates $Q^*$ and three representative static bubble sizes. $A_Q^*/A_s^*$ decreases monotonically and approximately linearly with increasing flow rate and does not depend significantly on the bubble size. A least-square linear fit $A_Q^*/A_s^* = K_0 - K^* Q^*$ yields $K_0 = 0.98$ and $K^* = 4.1 \!\times\! 10^{4}$ s m$^{-3}$. The intercept $A_Q^*/A_s^* =0.98$ for $Q^*=0$ suggests a steep increase towards $A_Q^*/A_s^* =1.00$ for vanishing flow rates, which was also reported by Franco-Gómez et al. (Reference Franco-Gómez, Thompson, Hazel and Juel2017).

The compression ratio $A_Q^*/A_s^*$ can be estimated by considering that the fluid pressure acting on the bubble increases rapidly from atmospheric pressure $P_A^*$ to $P_A^* + G^* l^* + G_t^* l_t^*$ when the flow is imposed. $G^* = 12 \mu Q^* / (W^* H^{*3})$ is the pressure gradient in the channel, $l^*$ is the distance initially separating the bubble from the channel outlet, typically $110$ cm in our experiments and $G_t^* = 8 \mu Q^* / ({\rm \pi} r_t^{*4})$ is the pressure gradient in the tube of radius $r_t^* = 2$ mm and length $l_t^* = 120$ cm connecting the channel outlet to the reservoir at atmospheric pressure; see figure 1. Assuming ideal gas behaviour and neglecting wetting films for simplicity, we get a compression ratio $A_Q^* / A_s^* = (1+ (G^*l^* +G_t^*l_t^*)/P^*_A)^{-1}$. This expression is plotted as a function of flow rate in figure 2(b). The discrepancy with the experimental data is due to the increase in projected area resulting from the increase in wetting film thickness when the bubble propagates, which partly compensates for the reduction in bubble volume due to compression.

The bubble does not recover its initial projected area $A_s^*$ when approaching the channel outlet partly because the pressure is still larger than $P_A^*$. However, the remaining pressure head $G_t^* l_t^*$ is too small to account for the modest $2\,\%$ increase observed in figure 2(a). The missing area is consistent with air diffusion across the oil–air interface during bubble propagation. Experiments for $A_s^*=3.9$ cm$^2$, where the flow rate was interrupted when the bubble was close to the channel outlet, showed that the projected area of the bubble at rest after propagation was less than the initial value $A_s^*$ by $0.35\,\%$ to $11\,\%$ for flow rates ranging from $10$ ml min$^{-1}$ to $477$ ml min$^{-1}$. This suggests that the solubility of air in oil increases with pressure, as previously reported by Franco-Gómez et al. (Reference Franco-Gómez, Thompson, Hazel and Juel2017). Hence, in our bubble propagation experiments, air diffusion into the silicone oil helps to retain an approximately constant bubble area after initial compression. The maximum increase in bubble area was $7\,\%$ for the highest flow rates, which was small enough to avoid a measurable change in bubble propagation; see also § 4.3.2 and appendix B for a discussion of the influence of the bubble size.

In order to account for these effects, we parametrised our bubble propagation experiments according to the inflow equivalent radius of the bubble $r_Q = (2 \sqrt {A_Q^*/ {\rm \pi}}) / W^*$ and we used the relation presented in figure 2(b) to obtain bubbles with the required value of projected area $A_Q^*$. Note that compression ratios similar to figure 2(b) were found for the initially symmetric bubbles prepared with the protocol described in § 3.

4.1. Influence of the initial bubble shape

Figure 4 shows the evolution of bubbles propagating from rest at a flow rate $Q^* = 186$ ml min$^{-1}$ ($Q = 0.07$), with non-dimensional time elapsed since flow initiation, $t = 2 U_0^* t^*/ W^*$. Each row of images corresponds to a time sequence of bubble evolution of an initial bubble shape quantified by its non-dimensional width $d$ measured at $t = 0$ (just before flow initiation) which ranges from $55.5\,\%$ (a) to $63.7\,\%$ (f).

Figure 4. Time evolution of bubbles of different initial shapes propagating from rest at a flow rate $Q^* = 186$ ml min$^{-1}$ ($Q = 0.07$) imposed at $t^*=0$. Each row of top-view images corresponds to a time sequence with a different value of non-dimensional initial bubble width $d = 2d^*/W^*$ in the range $55.5\,\% \le d \le 63.7\,\%$. The grey band in each image indicates the rail. The non-dimensional time $t = 2 U_0^* t^*/ W^*$ elapsed since flow initiation is indicated in each snapshot, where $U_0^* = Q^*/(W^* H^*)$. The dimensional time $t^*$ is indicated in the last snapshot of each time sequence. The inflow bubble size is $r_Q = 0.54$ (after flow-induced compression) and $d$ is measured at $t = 0$ (just before flow initiation). The final states are either a single on-rail bubble (a–c), which is slightly asymmetric at this value of $Q$, or two off-rail bubbles (d–f). See also supplementary movie 1 available at https://doi.org/10.1017/jfm.2020.844.

For all initial bubble shapes, the long-term outcome falls into only two categories: (i) a single, steadily propagating bubble straddling the rail (final frames in figure 4a–c); or (ii) two bubbles on opposite sides of the rail each moving at a constant velocity, where the larger, leading bubble propagates faster than the smaller, trailing one so that their relative distance increases in time (final frames in figure 4d,e). For $Q=0.07$, the single bubble straddling the rail at late times in figure 4(a–c) is slightly asymmetric about the centreline of the channel, and this state will henceforth be referred to as ‘on-rail, asymmetric’. The two bubbles propagating on opposite sides of the rail at late times in figure 4(d,e) are termed ‘off-rail’ bubbles because they do not straddle the rail, only partially overlapping it from one side.

These long-term outcomes are associated with an organised set of transient evolutions that are determined by the initial deformation of the centred bubble into single, triple or double-tipped bubbles as $d$ increases (see second column of snapshots in figure 4). Following this initial deformation, the bubble evolution may feature complex bubble breakup and aggregation events. The wide range of transient evolutions obtained for relatively modest changes in $d$ ($55.5\,\% \le d \le 63.7\,\%$) is fully reproducible. In some cases, the final state features compound bubbles resulting from the aggregation of two or more bubbles (figure 4c,d) which persist until reaching the end of the channel at $Q = 0.07$. However, slow drainage of the thin oil film separating the different bubbles would eventually result in their coalescence into simple bubbles.

The simplest evolution scenarios towards single and two-bubble final states occur for the most slender (figure 4a) and widest (figure 4f) initial bubble shapes, respectively. In figure 4(a), the centred bubble deforms symmetrically about the channel centreline until $t = 3.5$, after which it migrates sideways to reach the asymmetric, on-rail state at late times ($t=7.2$). When the flow is initiated, the front of the bubble located over the rail is pulled forward. By $t=2.1$, dimples on the interface (i.e. regions of negative curvature) have formed at the edges of the rail. These dimples are subsequently advected towards the rear part of the bubble (in the frame of reference of the moving bubble) where they vanish shortly after $t=2.6$, so that an approximately symmetric, elongated bubble with entirely positive curvature emerges at $t=3.5$.

In figure 4(f), the bubble initially features a wide front that encroaches into the deeper regions of the channel on either side of the rail. When flow is initiated, two fingers form in the deeper side channels separated by a dimple of negative curvature over the rail ($t=1.8$), which is essentially the reverse of the initial deformation in figure 4(a). The dimple deepens into a cleft ($t=2.4$) which continues to grow as the fingers lengthen, until the bubble breaks in two ($t=3.0$). Because of the small, unavoidable asymmetry in the initial bubble shape discussed in § 3 (here, centroid offset $y_c = -0.005$ at $t=0$), the two fingers do not lengthen at precisely the same rate and thus, the newly formed bubbles propagating on opposite sides of the rail have slightly different sizes (the projected area of the bottom bubble is $5$ % larger than that of the top one). We find that larger bubbles always propagate faster than smaller ones within the range of flow rates and bubble sizes investigated experimentally and we quantify this effect in § 4.3.2. In figure 4(f), the larger bubble has visibly moved ahead of the smaller one by $t=6.0$ and the two bubbles continue to separate, which indicates that the change in topology is permanent. Both bubbles propagate steadily by $t=9.4$, with the leading bubble propagating $2\,\%$ faster than the trailing one.

Prior to flow initiation, the bubble is characterised by a capillary-static pressure distribution which depends on the local in-plane curvature of the bubble and depth of the channel through the Young–Laplace relation. If the bubble is left to relax, local pressure gradients generate flows in the viscous fluid which act to broaden slender bubbles as shown in § 3. However, the imposition of a flow rate driving the bubble forward tends to override this initial capillary-static distribution. The bubble interface evolves via movement in its normal direction and its normal velocity is largest when the interface is orthogonal to the direction of imposed flow along $x$. In addition, the presence of the rail means that the mobility of the fluid is reduced over the rail compared with the deeper side channels. Hence, if a portion of the bubble front encroaches into the deeper side channels, these regions of interface can advance more rapidly than the part overlapping the rail, provided that these portions are oriented so that their normal is sufficiently aligned with the $x$ direction. In order for bulges to form, these portions of interface also need to be sufficiently wide to overcome capillary forces, as is the case in figure 4(f). If these criteria are not fulfilled, the portion of interface over the rail will advance most rapidly to form a central bulge as shown in figure 4(a).

For initial bubble shapes of widths intermediate between the slender and wide bubbles shown in figures 4(a) and 4(f), the front of the bubble initially deforms into a shape that is hybrid between those described above, i.e. by simultaneously developing a central bulge as well as bulges in the deeper side channels. This intermediate initial bubble deformation is associated with the most complex bubble evolutions, which occur within a narrow range of initial shapes, $(56.8 \pm 0.1)\% \le d \le (57.5 \pm 0.1)\%$. We refer to the developing bulges as ‘side’ and ‘middle’ tips as labelled in figure 4(b) ($t=2.4$). As $d$ is increased, the growth rate of the side tips increases relative to that of the middle tip, which is reduced ($t=2.1\text {--}2.6$ in figure 4b–e). The competition between side and middle tips is a key factor in the evolution of the bubble. In figure 4(b), the middle tip grows and the side tips retract, so that by $t=4.0$, the bubble evolves similarly to that in figure 4(a). By figure 4(c), the rate of growth of the side tips has increased and is only marginally smaller than that of the middle tip, resulting in a breakup into three bubbles ($t=3.2$). The larger middle bubble, which is approximately centred on the rail, pulls the smaller side bubbles in behind it where pressure is lower ($t=3.9$) to form a centred compound bubble ($t=5.1$), which finally drifts sideways to reach an asymmetric, on-rail state ($t=7.6$).

The transition in long-term behaviour, from a single on-rail bubble to two off-rail bubbles, occurs between figures 4(c) and 4(d) and coincides with the side tips overtaking the middle tip prior to bubble breakup. In figure 4(d), this leads to a breakup into two side bubbles propagating ahead of a smaller centred bubble ($t=3.2$). This smaller bubble is stretched across the rail as it is attracted similarly towards both side bubbles ($t=3.9$), so that in turn, it breaks into two small bubbles which are pulled in behind the side bubbles ($t=4.8$). In figure 4(e) the side tips have become dominant so that the bubble breaks into two parts. In both cases, the remaining two-bubble evolution is similar to that in figure 4(f).

4.2. Influence of the flow rate

We now turn to the influence of the flow rate on bubble evolution from a prescribed initial condition and successively examine the limits of slender and wide initial bubbles.

4.2.1. Initially slender bubbles

Figure 5 shows the evolution of initially slender bubbles (within a narrow range of widths $53\,\% \le d \le 55\,\%$) as a function of flow rate, where each row of images corresponds to a time sequence from the instant $t = 0$ when a fixed value of the flow rate was imposed in the range $0.005 \le Q \le 0.18$. In order to maintain the inflow bubble size at a fixed value $r_Q = 0.54$, the size of the bubble before flow-induced compression needs to increase with $Q$; see the initial shapes at $t=0$ in the column of magenta frames in figure 5 where the bubble length increases while $d$ remains approximately constant.

Figure 5. Time evolution of initially slender bubbles of non-dimensional width $53\,\% \le d \le 55\,\%$ with $r_Q=0.54$, where each row of top-view images corresponds to a time sequence at a different value of the flow rate imposed at $t^* = 0$ from $Q^* = 13 (Q = 0.005$) to $Q^* = 477$ ml min$^{-1}$ ($Q = 0.18$). The time labelling is the same as in figure 4. In (f), the transient evolution is shown on a time scale similar to previous sequences on the first line ($t\le 7.7$) and the sequence of images on the second line shows one period of oscillation of the final, time-periodic state ($23.82 \le t \le 24.33$). Between $t=7.7$ and $t=23.82$, the compound bubble coalesces into a single bubble due to the drainage of the oil film separating the aggregated bubbles, which is enhanced by the bubble oscillations.

For all values of $Q$, the long-term outcome is a single (simple or compound) bubble. We observe three types of steadily propagating, invariant bubbles: strongly asymmetric bubbles, which partially overlap the rail from one side, referred to as ‘off rail’ (final frames in figure 5a,b), a bubble which straddles the rail symmetrically about the channel centreline, referred to as ‘on-rail, symmetric’ (final frame in figure 5c) and ‘on-rail, asymmetric’ bubbles (final frames in figure 5d,e). In figure 5(f), the bubble exhibits periodic oscillations with the green frames showing one period. This state is referred to as ‘on-rail, asymmetric, oscillatory’ and is related to previously studied bubble and finger oscillations induced by channel depth perturbations (Pailha et al. Reference Pailha, Hazel, Glendinning and Juel2012; Jisiou et al. Reference Jisiou, Dawson, Thompson, Mohr, Fielden, Hazel and Juel2014; Thompson, Juel & Hazel Reference Thompson, Juel and Hazel2014). It features the advection of interface dimples along the side of the bubble that is furthest from the channel centreline (top region in the final frames of figure 5f). These dimples are periodically generated at the point where the front of the bubble meets the edge of the rail.

All bubbles in figure 5 initially deform symmetrically about the centreline of the channel. This key feature of the bubble evolution will be interpreted in terms of the unstable steadily propagating states of the one-bubble system in § 5.3.2. For the smallest values of $Q$ (figure 5a,b), the bubble initially broadens in response to the imposed flow because surface tension forces are dominant and the bubble evolves as though relaxing to static equilibrium. This results in the bubble migrating off the rail without breaking for $Q=0.005$ (figure 5a), whereas for $Q=0.01$ (figure 5b), the decreasing influence of surface tension is associated with the growth of a dimple of negative curvature on the front part of the bubble ($t=1.8$), which leads to breakup into two bubbles on opposite sides of the rail ($t=3.1$). These have different sizes because of the small, unavoidable asymmetry in the original bubble shape (centroid offset $y_c = 0.002$ at $t=0$). The smaller trailing bubble migrates across the rail ($t=6.4$) and is pulled in behind the larger leading bubble resulting in an off-rail compound bubble ($t=11.5$). In contrast with figure 4(c,d), we observe coalescence into a simple bubble before the bubble reaches the end of the channel ($t=13.8$) because the oil film separating the bubbles has sufficient time to drain due to the smaller flow rate. The bubble velocity is the same before and after coalescence.

As $Q$ increases, the early-time broadening of the bubble promoted at low $Q$ is progressively replaced by the growth of a bulge over the rail. This is illustrated in figure 5(c–f) which shows transient bubble evolutions similar to figure 4(a) until the bubble reaches an elongated, symmetric shape, see e.g. final frame in figure 5(c) and yellow frames in figure 5(d–f). We find that the next stage of evolution may involve breakup and aggregation events following tip splitting beyond a threshold value of the flow rate $Q_{ts} = 0.085 \pm 0.005$. This is illustrated in figure 5(e,f), where the bubble breaks into two parts of very different sizes which subsequently aggregate. These events are increasingly likely the larger the flow rate and we refer to § 4.4.2 for a discussion of this late stage of evolution.

Note that the compression of the bubble occurs gradually over increasing time intervals as flow rate increases, so that in figure 5(f) ($Q=0.18$) for example, an approximately constant projected area is only reached for $t>4.0$. We provide evidence in appendix B that increasing the size of the bubble beyond $r_Q=0.54$ does not influence its evolution significantly and thus, the bubble evolutions presented in figure 5 are not expected to differ from the idealised case of a bubble of fixed volume.

4.2.2. Initially wide bubbles

Figure 6 shows the evolution of initially wide bubbles ($d = 63.4\,\%$) as a function of the flow rate in the range $0.005 \le Q \le 0.18$ with a layout similar to figure 5. For low flow rates (figure 6a,b), the bubble evolution is similar to that of the slender bubbles in figure 5(a,b), resulting in a steadily propagating off-rail bubble. This suggests that the bubble evolution is insensitive to the initial shape of the bubble within this range of flow rates.

Figure 6. Time evolution of initially wide bubbles with non-dimensional width $d = 63.4\,\%$ and $r_Q=0.54$. Each row of top-view images corresponds to a time sequence at a different value of the flow rate imposed at $t^* = 0$ between $Q^* = 13 (Q = 0.005$) and $Q^* = 477$ m min$^{-1}$ ($Q = 0.18$). The time labelling is the same as in figure 4.

However, for higher flow rates (figure 6c,d), the long-term outcome is two steadily propagating off-rail bubbles, where the leading bubble is faster than the trailing one so that their relative distance increases with time. The early-time evolution is similar to that of figure 6(b), in that the bubble breaks into two parts of similar but not identical sizes. Whereas in figure 6(b), the smaller trailing bubble migrates across the rail and aggregates with the larger leading bubble, in figure 6(c), the two bubbles remain on opposite sides of the rail and separate, like in figure 4(f). The transition between the time evolutions leading to aggregation (figure 6b) and separation (figure 6c) occurs at a threshold value of the flow rate $Q_c = 0.017 \pm 0.001$. We will show in § 5.4 how these distinct evolutions are linked to the presence of two-bubble edge states.

In figure 6(c), the smaller bubble exhibits transient oscillations (at around $t=10.1$) before reaching a state of steady propagation. As the flow rate increases, these oscillations increase in amplitude and promote the repeated breakup of the smaller trailing bubble, as shown in figure 6(d) where tiny bubbles created through breakup aggregate with the bubble that spawned them. These secondary breakup events occur for $Q>0.135 \pm 0.005$ and increase in complexity and unpredictability as $Q$ increases, leading to complex final states involving more than two bubbles for the highest flow rates investigated.

4.3. Summary of bubble evolution

We now proceed to summarise the transient bubble evolutions reported in §§ 4.1 and 4.2 in a phase diagram and map the different long-term behaviours of single bubbles onto a bifurcation diagram.

4.3.1. Phase diagram of bubble evolution

The influence of the initial bubble width and the flow rate is summarised in a phase diagram in figure 7(a), where we classify more than two hundred experimental time evolutions into eight categories, each denoted by a different symbol. Each type of transient evolution occurs in a simply connected region of the diagram. The diagram is reproduced in figure 7(b) without the data points to allow space for annotations.

Figure 7. (a) Phase diagram of the different time evolutions of a bubble with $r_Q = 0.54$ in terms of its initial non-dimensional width $d$ and the non-dimensional flow rate $Q$. (b) Reproduction of (a) without the data points to allow for annotations. We distinguish eight different types of transient evolution, each denoted by a symbol illustrated by a snapshot of the associated initial deformation on the right-hand side of the diagram in (a). The piecewise linear solid lines labelled $d_i$ ($i=1-5$) in (b) indicate approximate boundaries between different regions. The background colour indicates the long-term outcome of the bubble evolution as illustrated in the legend on the right-hand side of the diagram in (b). The hatching at high flow rate highlights the increasing diversity of final outcomes as $Q$ increases (see figure 11). Figure 4 is representative of the bubble evolution along a horizontal line of constant $Q$ for $Q > 0.025 \pm 0.005$, while figures 5 and 6 correspond to the behaviour along vertical lines of constant $d$ for slender ($d<d_3$) and wide ($d>d_3$) initial bubbles, respectively.

For $Q \le 0.025 \pm 0.005$, the evolution of the bubble does not depend on initial conditions, and three different types of evolution are observed. When $Q > 0.025 \pm 0.005$, the variation of the initial bubble width yields six distinct types of early-time evolution, previously introduced in figure 4. The rich variety of transient evolutions is concentrated over a small range of initial bubble widths for moderate flow rates above this threshold, and the ranges of $d$ for which these initial evolutions are observed widen as $Q$ increases.

The long-term outcomes are indicated by the background colour: blue for single off-rail bubbles, green for single on-rail bubbles (symmetric, asymmetric or oscillating) and red for two off-rail bubbles (or occasionally three as will be shown in § 4.4.1). The yellow region corresponds to bubbles which may undergo tip splitting following their early-time symmetric evolution for $Q>Q_{ts}$, resulting in either a single on-rail bubble or two off-rail bubbles; see § 4.4.2. The hatching highlights the increasing diversity of final outcomes (with possibly more than two bubbles) as the flow rate increases; see § 4.4.2.

For all initial bubble widths investigated, the long-term behaviour of the bubble transitions from one off-rail bubble to two off-rail bubbles at a flow rate $Q_c = 0.017 \pm 0.001$. For $Q < Q_c$, the single off-rail bubble is reached by migrating off the rail without breaking ($\bullet$) or by breaking into two parts which subsequently aggregate ($\circ$). By contrast, the two bubbles separate after breakup for $Q_c < Q \le 0.25 \pm 0.005$ ($\square$).

For $Q > 0.025 \pm 0.005$, the boundary $d_3$ separates the red region, where the long-term outcome is always two (or three) off-rail separating bubbles, from the green and yellow regions where a single on-rail bubble is a possible long-term outcome. In the green and yellow regions, dash-dotted lines in figure 7(b) indicate the critical values of flow rate at which different single on-rail bubble invariant states emerge. These are detailed in an experimental bifurcation diagram in § 4.3.2.

The boundary $d_3$ corresponds to the initial bubble width for which side and middle tips of the three-tipped transient bubbles introduced in figure 4 grow at the same rate during early-time evolution. For $d_1 \le d \le d_3$ the middle tip grows faster than the side tips whereas the converse is true for $d_3 \le d \le d_5$. For $d<d_2$, the bubble does not break following initial deformation, whereas for $d_2<d<d_4$ it breaks in three ($\bigtriangleup$ and $\bigtriangledown$) and for $d>d_4$ it breaks in two. The value $d_1$ indicates the development of side tips (the maximum number of intersections of the bubble contour with a line perpendicular to the side walls increases from two ($\blacktriangleleft$) to six ($\blacktriangledown$)). Similarly, the value $d_5$ indicates the disappearance of the middle tip (the maximum number of intersections of the bubble contour with a line perpendicular to the sidewalls decreases from six ($\vartriangleright$) to four ($\square$)).

4.3.2. Single-bubble invariant modes of propagation

Figure 8 shows the single-bubble invariant modes of propagation observed experimentally in terms of their non-dimensional centroid offset $y_c$ (a) and velocity $U_b$ (b), both plotted as functions of $Q$ for three values of the bubble size $r_Q$. These bifurcation diagrams extend those reported by Franco-Gómez et al. (Reference Franco-Gómez, Thompson, Hazel and Juel2017, Reference Franco-Gómez, Thompson, Hazel and Juel2018) for $Q<0.014$ to larger values of the flow rate up to $Q=0.19$. As shown in figure 8(a), the branch of off-rail asymmetric bubbles (AS1) extends from the capillary-static limit to the largest flow rate investigated. The AS1 mode of propagation was reached by initially positioning the bubble off-rail prior to imposing the flow rate $Q$ (following the same procedure as in § 2.2).

Figure 8. Experimental bifurcation diagram where the non-dimensional bubble $y$-offset $y_c = y_c^*/(W^*/2)$ (a) and bubble velocity relative to the flow far from the bubble $U_b = U_b^*/U_0^* = U_b^* W^* H^*/Q^*$ (b) are plotted as functions of the non-dimensional flow rate $Q$ for three bubble sizes $r_Q = 0.54$, $0.8$ and $1.2$. The off-rail branch AS1 and the three on-rail branches, which are symmetric (S3) asymmetric (AS1) and oscillating asymmetrically (PO) are illustrated by the final-state snapshots of figure 5(b–f), respectively. H1 indicates the first appearance of on-rail bubbles, PF and H2 refer to the bifurcations between the on-rail branches. For the periodic orbit PO, the reported values of $y_c$ and $U_b$ are averaged over one period.

The other invariant modes of propagation are those classified as ‘on rail’ in figure 7. As shown in figure 8(a), the symmetric branch S3 emerges beyond a threshold flow rate $Q_{H1}$. We obtained the smallest value $Q_{H1}=0.016$ by imposing the flow while the bubble was at rest inside the centring device depicted in figure 1(c). However, the true value could be smaller because the long-term outcome is influenced by the expansion ratio of the centring device which determines the bubble shape at its exit and thereby the bubble's dynamical evolution. The subsequent transition to an on-rail, asymmetric state (AS2) at $Q_{PF}=0.0435$, obtained from a linear fit of $y_c^2(Q)$ near the transition, is consistent with an imperfect supercritical pitchfork bifurcation. Finally, the transition to an on-rail, asymmetric and oscillatory state (PO) at $Q_{H2}= 0.118$ is consistent with a supercritical Hopf bifurcation; see appendix C for measurements of the amplitude and frequency of these oscillations. The variation of these critical flow rates with bubble size was less than $3\,\%$ within the range investigated.

As shown in figure 8(b), off-rail bubbles propagate faster than on-rail bubbles at a given flow rate due partly to lower viscous resistance in the deeper parts of the channel. Also, bubble velocity systematically increases with bubble size over the range of bubble sizes and flow rates investigated. This is consistent with the bubble evolutions where the larger leading bubble propagates faster than the smaller trailing one; see e.g. figure 4(f).

4.4. Sensitivity of the transient evolution

We now turn our attention to dynamical evolutions of the bubble that were particularly sensitive to either the initial condition of the system or imperfections of the rail. We distinguish situations where sensitivity occurs at moderate flow rate ($Q<Q_{ts}$; see figure 7) from those where sensitivity occurs at high flow rate ($Q \ge Q_{ts}$; yellow region and hatched region in figure 7).

4.4.1. Moderate flow rates

We observed that the transient evolution of the bubble exhibited enhanced sensitivity to the value of the initial offset of the bubble from the channel centreline $y_c$, for flow rates slightly larger than the value $Q_c=0.017 \pm 0.001$ which marks the transition between aggregation (figure 6b) and separation (figure 6c). We found that the final outcome of the system is influenced by the value of the size ratio of the two bubbles formed after breakup, which depends on $y_c$. For example, at $Q=0.020$, the two bubbles resulting from breakup only separated when the leading bubble was at most $10\,\%$ larger than the trailing bubble ($|y_c| < 0.002$), while they aggregated when its size exceeded that of the trailing bubble by $20\,\%$ to $80\,\%$ ($0.003 < |y_c| < 0.008$). For intermediate values, the smaller trailing bubble could stretch sufficiently across the rail to break in two, leading to a final state where three bubbles, organised in order of decreasing size and velocity, separate indefinitely. For this reason, extra care was taken when estimating $Q_c$ to minimise the asymmetry of the initial bubble to offsets $|y_c|< 0.002$, which was achieved through accurate levelling of the channel in order to minimise migration of the bubble during the relaxation process described in § 3.

In addition, we observed reduced experimental reproducibility of evolutions where breakup events had a propensity to generate tiny bubbles. This is illustrated in figure 9 where successive experiments, performed at the same value of $d=57.4\,\%$ for $Q=0.07$, evolve towards different final outcomes involving either two or three separating bubbles. The difference between the time evolutions of figures 9(a) and 9(b) (superposed in figure 9c) is caused by differences in relative bubble sizes after breakup, stemming from small unavoidable differences in the initial bubble area of the order of $0.5\,\%$ ($t=0$), which is the limit of reproducibility in the experiments.

Figure 9. (a,b) Time evolution of a bubble of initial width $d = 57.4\,\%$ propagating from rest at a flow rate $Q^* = 186$ ml min$^{-1}$ ($Q = 0.07$) with $r_Q=0.54$. (c) Superposition of (a) in blue and (b) in red. Both experiments performed for the same experimental parameters feature breakup into three bubbles. However, the size of the middle bubble differs visibly, which leads to different long term outcomes. The time labelling is the same as in figure 4.

4.4.2. High flow rates

We mentioned in § 4.2.1 that initially slender bubbles may break up via tip splitting for flow rates beyond the threshold value $Q_{ts}= 0.085 \pm 0.005$, following a reproducible symmetric early-time evolution. This second stage of bubble evolution is not reproducible in the experiment, as illustrated in figure 10 where two successive experiments at the same parameters evolve towards radically different final outcomes.

Figure 10. (a,b) Time evolution of a bubble of initial width $d = 57.7\,\%$ propagating from rest at a flow rate $Q^* = 265$ ml min$^{-1}$ ($Q = 0.10$) with $r_Q=0.54$. (c) Superposition of (a) in blue and (b) in red. Both experiments performed for the same experimental parameters lead to the formation of an approximately symmetric transient compound bubble ($t=4.5$), which either undergoes off-centre (a) or centred (b) tip splitting leading to different final outcomes. The time labelling is the same as in figure 4.

The early-time symmetric evolution, which here involves breakup and aggregation, leads to an approximately symmetric compound bubble by $t=4.5$. In both figures 10(a) and 10(b), a dimple forms at the bubble tip ($t=5.2$), but at different locations (more centrally in (b) than in (a)). This results in a second breakup which produces two bubbles of similar sizes in figure 10(b) and very different sizes in figure 10(a) ($t=5.9$). In the latter case, the system then evolves as in figure 5(e), resulting in a steadily propagating on-rail, asymmetric bubble ($t=9.8$). However, the two bubbles in figure 10(b) ultimately propagate on opposite sides of the rail ($t=9.8$), with the larger leading bubble propagating faster so that they separate with increasing time. This indicates that the position of the dimple relative to the bubble tip is a key factor in predicting the long-term evolution of the system.

The time evolution becomes increasingly complex and unpredictable as the flow rate is increased, typically involving an increasing number of breakup and aggregation events (hatched region in figure 7). This is illustrated in figure 11 where five successive experiments performed at identical parameter values corresponding to the yellow hatched region in figure 7 evolve towards five different final outcomes. At $t=0$, the relative difference in bubble area is within $0.5\,\%$ and the centroid positions along $x^*$ are within $1$ mm of each other. The time sequences in figure 11(a,b) result in a single on-rail bubble ($t=17.0$) and in two off-rail bubbles ($t=13.0$) because tip splitting occurs off centre and on the centreline, respectively, similarly to figure 10. The flow rate is large enough for the on-rail, asymmetric bubble obtained in figure 11(a) to oscillate. Figures 11(d) and 11(e) provide examples of more complex time evolutions after centred tip splitting, involving several breakup and aggregation events. In figure 11(d), an on-rail compound bubble (asymmetric and oscillating) is recovered followed by a tiny slowly propagating bubble, while in figure 11(e), a large off-rail bubble is followed by two much smaller bubbles on either side of the rail, which will ultimately be sorted by size due to their differing velocities. A recurrent event during these complex evolutions is the rapid stretching and breakup of bubbles trailing behind a leading bubble, whose products may in turn recombine with the bubbles ahead depending on their sizes and positions (figure 11(d), $t=9.4 - 10.8$). In contrast, in figure 11(c), tip splitting does not occur following initial symmetric deformation and the bubble remains on rail, wobbling in a non-periodic manner with its centroid remaining on average at $y_c = 0$ until the end of the channel. This scenario only occurs for $Q > 0.155 \pm 0.005$ and becomes increasingly likely the larger the flow rate.

Figure 11. Time evolution of a bubble of initial width $d = 57.3\,\%$ propagating from rest at a flow rate $Q^* = 477$ ml min$^{-1}$ ($Q = 0.18$) with $r_Q=0.54$. All experiments are performed for the same experimental parameters and lead to the formation of a transient approximately symmetric bubble ($t=4.6$) which, except in (c), undergoes tip splitting and evolves towards different final outcomes. The time labelling is the same as in figure 4. See also supplementary movie 2.

The lack of reproducibility at high flow rate (figures 10 and 11) appears to stem from the increased sensitivity to both initial conditions and perturbations in the channel such as defects on the rail (of the order of the tape roughness). In order to test the influence of rail defects, we varied the portion of rail travelled by the bubble by initiating their propagation from different initial positions along the rail at high flow rate ($Q>Q_{ts}$) while keeping all other experimental parameters constant. By doing so, the bubble was subjected to different patterns of defects while propagating. We found that different initial bubble positions along the rail usually resulted in different locations of the dimple formed shortly after the bubble finishes its symmetric evolution. However, the associated tip splitting always occurred at the same time after flow initiation for a given flow rate regardless of the initial bubble position. In contrast, when repeating experiments from the same initial position, the bubble generally followed the same late-time evolution, with rare exceptions. This suggests that the location of the dimple leading to tip splitting coincides with the location of the largest rail defect ahead of the bubble once it reaches the end of its symmetric evolution. This is reminiscent of the findings of Tabeling, Zocchi & Libchaber (Reference Tabeling, Zocchi and Libchaber1987) in a rectangular Hele-Shaw channel (in the absence of a rail), where symmetric fingers propagating at sufficiently large flow rates developed dimples at fixed positions in the channel due to localised perturbations on the top and bottom glass plates.

However, in our experiments, the same initial bubble position (within $<$1 mm) could also result in different dimple locations and therefore lead to different long-term outcomes. These relatively rare cases are those previously shown in figures 10 and 11. Ten experiments performed under the experimental conditions of figure 11 led to two off-centred tip splittings, five centred tip splittings and three cases without tip splitting like in figure 11(c). We suspect that tip splitting does not occur if the bubble does not encounter a sufficiently large rail defect at the end of its symmetric evolution phase. In all cases, the bubble fate is determined by its exact shape at that time, which exhibits small variations between experiments stemming from small differences in initial bubble shape at $t=0$.

5.1. Depth-averaged model for bubble propagation

The depth-averaged model extends McLean & Saffman's (Reference McLean and Saffman1981) approach to account for a non-uniform channel height. We have previously used the same depth-averaged model to study the propagation of a semi-infinite air finger (Thompson et al. Reference Thompson, Juel and Hazel2014; Franco-Gómez et al. Reference Franco-Gómez, Thompson, Hazel and Juel2016) and a single closed air bubble (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) within the constricted channel. Two extensions are required here: (i) the simulation of multiple air bubbles; and (ii) the inclusion of topological changes in the form of coalescence and breakup.

We work in a frame moving with the centroid position of the whole collection of bubbles and non-dimensionalise the physical system as shown in figure 1: scaling the $x^*$ and $y^*$ directions by $W^*/2$, the $z^*$ direction by $H^*$ and the velocity by $U_0^* = Q^*/(W^*H^*)$. The non-dimensional variable channel height is approximated by the smoothed tanh profile

(5.1)\begin{equation} b(y) = 1 - \tfrac{1}{2}h[\tanh(s(y + w)) - \tanh(s(y - w))], \end{equation}

where $h = h^*/H^*$ and $w = w^*/W^*$ are the non-dimensional height and width of the rail, respectively, and where $s$ sets the sharpness of the sides of the rail. We use the experimental values $h=0.024$ and $w = 0.25$, and choose $s=40$.

We describe the flow of the viscous fluid in terms of a depth-averaged velocity, assuming Stokes flow (i.e. $\rho U_0^*W^*/(\mu \alpha ^2) \ll 1$) and a large aspect ratio (i.e. $\alpha = W^*/H^* \gg 1$). Within the fluid domain, the two-dimensional depth-averaged velocity $\boldsymbol {u}$ satisfies

(5.2a,b)\begin{equation} \boldsymbol{u} = -b^2(y) \boldsymbol{\nabla} p, \quad \boldsymbol{\nabla} \boldsymbol{\cdot} (b(y)\boldsymbol{u})=0, \end{equation}

where the second equation is a statement of mass conservation. At the channel sidewalls, we enforce no-penetration conditions. The pressure is fixed to zero at the inflow, and a non-zero constant at the outflow to ensure the dimensionless volume flux is 2.

For both kinematic and dynamic conditions at the interface of each bubble, we assume that the air bubble fills the full height of the channel $b(y)$. Mass conservation then leads to the kinematic boundary condition

(5.3)\begin{equation} \frac{\partial \boldsymbol{R}}{\partial t}\boldsymbol{\cdot} \boldsymbol{n} + \boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{n} + b^2(y)\boldsymbol{\nabla} p\boldsymbol{\cdot} \boldsymbol{n} = 0,\end{equation}

at each bubble interface, where $\boldsymbol {R}$ is the two-dimensional position of a point on the bubble boundary, $\boldsymbol {n}$ an outward normal unit vector and $\boldsymbol {U} = (U_b(t),0)$ with $U_b(t)$ the dynamically varying frame speed. We allow for multiple incompressible bubbles, each of different internal pressures. Hence the dynamic condition at the boundary of bubble $i$ is

(5.4)\begin{equation} p_i - p = \frac{1}{3\alpha Q}\left(\frac{1}{b(y)}+\frac{\kappa}{\alpha} \right), \end{equation}

where $p_i$ is the internal pressure in bubble $i$, $p$ the fluid pressure, $\kappa$ is the lateral curvature of the bubble in $(x,y)$ plane and the transverse curvature contribution gives rise to the $1/b(y)$ term.

The pressure in each bubble, $p_i(t)$ is not known a priori, but instead is determined so that the dimensionless volume $V_i$ of each bubble remains constant, where $V_i$ is calculated assuming that the bubble occupies the full height of the channel

(5.5a,b)\begin{equation} V_i = \int_{\varGamma_i} b(y)\, \textrm{d}\kern0.7pt x\, \textrm{d} y= -\int_{\partial \varGamma_i} x b(y) \,\textrm{d} y, \quad A_i = \int_{\varGamma_i} \, \textrm{d}\kern0.7pt x\, \textrm{d} y. \end{equation}

Here, $\varGamma _i$ is the interior of bubble $i$ and $\partial \varGamma _i$ its bounding curve. In the model, the dimensionless area and volume of the bubble are almost identical because we assume that the bubble fills the available height, which differs from 1 by at most $2.6\,\%$. As discussed in § 2.2, the projected area of the bubble $A_{i}$ varies with flow rate in the experiments, but it remains approximately constant at a given flow rate.

We discretise the system in space using a finite-element method, implemented in the open source oomph-lib package (Heil & Hazel Reference Heil and Hazel2006). The spatial discretisation uses a piecewise quadratic, triangular mesh, fitted to the bubble boundary. Mesh deformation is handled by treating the nodes of the mesh as embedded in a fictitious elastic solid. Time stepping uses BDF1 or BDF2 (backwards difference formulae), with a new triangular mesh generated at least every 5 time steps, and a typical time step is $\delta t \approx 0.01$. The element distribution in the new mesh is determined by using a ZZ error estimator (Zienkiewicz & Zhu Reference Zienkiewicz and Zhu1992) based on the continuity of pressure gradients between elements. Typically, the minimum element area is approximately $10^{-7}$ and the maximum area is approximately $10^{-4}$. For a single-bubble simulation, the standard mesh resolution is around 4000 quadratic elements, whereas it is approximately 6000 elements in cases with two bubbles. The channel length is typically taken as six times its width, to allow for interactions between bubbles of moderate separation. The resulting set of discrete nonlinear algebraic equations is solved via a Newton method and the associated linear systems of equations are solved using SuperLU (Demmel, Gilbert & Li Reference Demmel, Gilbert and Li1999).

The experimental results presented in § 4 show frequent changes in bubble topology and we must therefore allow for both coalescence and breakup in our numerical simulations. To the best of our knowledge, there is no numerical method for solving our model equations in which bubble breakup and coalescence do not require some intervention at the moment of topological change. In our implementation, the intervention involves detection of a topological event, reconnection of the lines defining the bubble interface, creating a new bulk mesh and restarting the simulation with a new set of volume constraints. It is important to note that the only time derivatives in our model are those associated with motion of the interface. Hence, after a change in topology, we need only specify the new configuration of the interface, and do not require an initial condition for the velocity of the fluid. In our code, the intervention is completely automated and based on the criteria described in the following paragraphs.

As the region inside the bubble is treated simply as a region of constant pressure and fixed volume, the need for breakup is signalled by self-intersection of the bubble boundary, which here happens in finite time (see figure 12a,b). After each timestep, we check for self-intersection. If intersection has occurred, we remove the intersected segments of the interface, and reconnect the two other segments with straight lines to create two simply connected daughter bubbles (figure 12c). In contrast to breakup, coalescence does not occur inevitably in finite time. Instead we must set some criteria that triggers coalescence, and we choose to base this on a simple minimum distance between the interfaces of two different bubbles. If bubbles approach below this minimum (figure 12d), we remove the boundary nodes that are within the tolerance, and reconnect the remainder (figure 12e). Our topological ‘surgeries’ for both coalescence and breakup produce sharp corners on the bubble interface and are not volume conserving. To recover from the sharp corners, we take smaller time steps after a topological change, and remesh frequently. The interface typically regains a smooth shape very quickly (figure 12f). For coalescence, we set the target volume $V_i$ for the new bubble to be the sum of the two parent bubbles. For bubble breakup, we make a small adjustment to the target volumes for the two daughter bubbles (allocated in proportion to their measured volumes) in order to fully recover the volume of the parent bubble.

Figure 12. Mesh interventions for breakup and coalescence. (a) Bubble is simply connected before coalescence, $t=1.24$. (b) Self-intersection is detected, $t=1.28$. (c) Nodes in the intersected segment are removed, with the bubble reconnected via straight lines into two closed daughter bubbles, and the bulk is completely remeshed. (d) Two bubbles approach within the tolerance distance for coalescence, $t=11.022500$. (e) Nodes within the tolerance are removed, and the two bubbles joined via straight lines. The bulk region is remeshed. (f) After several small time steps, the sharp corners from (e) have smoothed into recoiling tendrils, $t=11.022515$.

We choose a dimensionless minimum distance for coalescence of $0.01$ which, although relatively large and significantly above the distance at which van der Waals forces become relevant, avoids the high computational cost of meshing very thin drainage regions. Indeed, once the drainage regions become thinner than the height of the channel our model is not expected to describe their behaviour accurately.

Comparison of bubble breakup between model and experiments in our system indicates that breakup is reasonably well approximated by our approach. In contrast, the approach to coalescence is different between model and experiments. As described in § 4, compound bubbles that can travel large distances before eventual coalescence are observed in the experiments. There is typically very little change in the propagation speed of the bubble and in the outline shape of the bubble front before and after coalescence. Therefore, we believe that the final state of the system is reasonably insensitive to details of the coalescence process, except perhaps in cases where bubble shapes are changing rapidly. The same general phenomena are observed in the simulations, but coalescence happens much more rapidly. We have confirmed that reducing the minimum distance for coalescence in the model to $0.001$ does not significantly affect the behaviour and merely increases the time taken for a compound bubble to coalesce.

In order to explore the time-dependent evolution of a single bubble, we perform time simulations starting from a variety of initially elliptical bubble shapes. Specifically, we define the initial bubble shape as

(5.6a,b)\begin{equation} x = l\cos(\theta),\quad y = d\sin(\theta),\quad \theta \in [0,2{\rm \pi}]. \end{equation}

We note that this initial bubble shape is symmetric about the line $y=0$, but the triangular mesh generated in the fluid bulk is not perfectly symmetric. All the initial value calculations presented here start from a single bubble with fixed initial area $A={\rm \pi} r^2 = {\rm \pi}l d$ with $r=0.54$, with the bubble volume calculated from (5.5a,b). We vary the flow rate $Q$ and the bubble width $d$ so that the phase diagram of the transient behaviour in the $(d,Q)$ plane can be compared with the experiment.

Steady solutions are calculated by setting the time derivatives to zero and again solving the discrete nonlinear system via a Newton method, using SuperLU (Demmel et al. Reference Demmel, Gilbert and Li1999) to solve the linear system of equations. Once a steady solution has been obtained the effects of the flow rate and volume can be mapped out using continuation and edge-tracking techniques (Skufca, Yorke & Eckhardt Reference Skufca, Yorke and Eckhardt2006; Gelfgat Reference Gelfgat2019). We assess the linear stability of the steady solutions by solving a discrete generalised eigenproblem using the Anasazi eigensolver from the Trilinos project (Heroux et al. Reference Heroux2005). Our sign convention is chosen so that eigenvalues with a positive real part correspond to unstable modes.

5.2. Transient evolution and phase diagram

The model described in § 5.1 captures qualitatively all the bubble evolutions observed experimentally. A total of 273 simulations performed for $Q \le 0.12$ are organised in a numerical phase diagram in figure 13. The phase diagram uses the same structure, symbols and colours as the experimental phase diagram shown in figure 7.

Figure 13. Numerical phase diagram of the different time evolutions as a function of the non-dimensional initial bubble width $d$ and non-dimensional flow rate $Q$. The initial bubble shape is an ellipse of non-dimensional projected area ${\rm \pi} ld$ and size $r = \sqrt {ld} = 0.54$. The layout, symbols, colours and notations are the same as in figure 7 including the threshold values $Q_{c(num)}$ and $Q_{ts(num)}$. The symbols distinguish different transient evolutions shown with inset images of the associated initial deformations in (b). A ninth type of initial deformation is defined ($\diamond$) which involves a tip-split finger and occurs in the upper left corner of the diagram, for small values of $d$ which were not probed in the experiment. The values $Q_{F(num)}$, $Q_{H1(num)}$, $Q_{PF(num)}$ and $Q_{H2(num)}$ indicated by arrows in (b) are values calculated from a linear stability analysis indicating bifurcations in the solution structure of the single-bubble system; see figure 15. The last three values mark transitions between different stable one-bubble invariant solutions. Specifically, the long term outcomes observed are off-rail bubbles for $Q< Q_{H1(num)}$, symmetric on-rail bubbles for $Q_{H1(num)} \le Q < Q_{PF(num)}$, asymmetric on-rail bubbles for $Q_{PF(num)} \le Q < Q_{H2(num)}$ and oscillatory on-rail bubbles for $Q \ge Q_{H2(num)}$.

The initial bubble deformations and long-term outcomes are qualitatively similar in both experimental and numerical phase diagrams, although the threshold flow rates do not coincide exactly. For example, the blue–red border separating one-bubble from two-bubble long-term outcomes occurs at a threshold flow rate $Q_{c(num)}=0.025 \pm 0.005$ (compared to the experimental value $Q_c=0.017 \pm 0.001$). The single-bubble invariant states arise in the same order as in the experiment as $Q$ is increased: off-rail bubbles ($Q < Q_{H1(num)}$); symmetric on-rail bubbles ($Q_{H1(num)} \le Q < Q_{PF(num)}$); asymmetric on-rail bubbles ($Q_{PF(num)} \le Q < Q_{H2(num)}$); and finally oscillatory on-rail bubbles ($Q \ge Q_{H2(num)}$); see § 4.3.2 for a discussion of the associated bifurcation structure.

Two examples of numerical time evolutions are shown in figure 14(a,b), which are representative of regions in figure 13 indicated by symbols $\bigtriangleup$ and $\bigtriangledown$, respectively. These time sequences illustrate the key role of the relative size of side and middle tips in determining the long-term outcome when the front of the bubble initially deforms into three tips, consistent with experimental findings.

Figure 14. Numerical time evolutions of an initially centred elliptical bubble with parameters $Q=0.07$ and $d=0.60$ ($\bigtriangleup$) (a) and $Q=0.10$ and $d=0.68$ ($\bigtriangledown$) (b); see figure 13. The dimensionless time is given at the bottom of each snapshot. The bubble size is $r=0.54$.

A distinction from the experimental phase diagram is the increased range of initial bubble widths $0.40 \le d \le 0.80$ in figure 13 compared to the experimental range in figure 7. For $d \le 0.48$ (not accessed experimentally), numerical simulations reveal only single-bubble long-term outcomes at moderate flow rates, as indicated by the blue-green border in figure 13. At high flow rates within this range of initial bubble widths, we identify a ninth type of bubble evolution where the bubble breaks in two parts of unequal sizes following asymmetric tip splitting ($\diamond$, red region in the top left corner of the phase diagram in figure 13).

At flow rates $Q>Q_{ts(num)}$, where $Q_{ts(num)} = 0.095 \pm 0.005$ is close to the value $Q_{ts}=0.085 \pm 0.005$ measured experimentally, the late-time evolution of the bubble typically involves tip splitting, like in the experimental evolutions of figures 10 and 11. In this yellow region of the phase diagram of figure 13, small differences in the details of tip splitting between two neighbouring data points may lead to either one or two-bubble final outcomes. More complex final outcomes occur in the hatched region as the complexity of the bubble evolution increases with flow rate. We hypothesise that the initial asymmetric tip splitting observed in the ninth region of bubble evolution ($\diamond$) in figure 13 is analogous to the tip splitting observed at late times in the yellow region but that it occurs sooner, superseding the initial symmetric deformation observed in the yellow region at early times.

The increasing sensitivity of the numerical model as $Q$ increases is reminiscent of the sensitivity to initial conditions and rail defects reported experimentally in § 4.4.2; and represented as a yellow and a hatched region in the phase diagram of figure 7. In the experiments, variations occurred between repetitions at fixed parameter values due to unavoidable differences in the initial conditions. In the numerical simulations the results are, of course, identical for the same initial conditions, but we find that smaller differences in initial conditions are required to provoke dynamically significant differences in the breakup as the flow rate increases. In particular, modest changes in relative bubble volumes after breakup can lead to quite different final outcomes, which may be a consequence of alterations in the associated invariant solution structure.

Overall, the qualitative agreement between numerical simulations and experiments is remarkable given the simplicity of the mathematical model and approximation of the initial experimental conditions by ellipses.

5.3. Single-bubble invariant modes of propagation

5.3.1. Bifurcation diagram

Steady solutions of the dynamical system are summarised in the numerical bifurcation diagram of figure 15 where the bubble centroid offset $y_c$ (a) and velocity $U_b$ (b) are each plotted as a function of the flow rate $Q$ for a fixed bubble volume $V=0.54^2 {\rm \pi}$. Stable and unstable solution branches are plotted with solid blue lines and dotted red lines, respectively. The stable solution branches are in agreement with the stable modes of propagation reported in the experimental bifurcation diagram of figure 8.

Figure 15. Numerical bifurcation diagram showing steadily propagating solutions for a single bubble. The non-dimensional bubble offset $y_c$ (a) and velocity $U_b$ (b) are plotted as a function of the non-dimensional flow rate $Q$ for a fixed bubble volume $V=0.54^2 {\rm \pi}$ rather than a fixed projected area. Solid blue and dashed red lines indicate stable and unstable solutions, respectively. Each solution branch is illustrated by a snapshot corresponding to a flow rate indicated by a black cross, except the S2b and S2c branches where snapshots are shown in figure 18. Solid circles indicate locations of bifurcation points. More details about the solution structure can be found in Keeler et al. (Reference Keeler, Thompson, Lemoult, Juel and Hazel2019).

The numerical simulations indicate that $Q_{H1(num)}=0.0301$ marks a subcritical Hopf bifurcation, $Q_{PF(num)}=0.0714$ a supercritical pitchfork bifurcation and $Q_{H2(num)}=0.0789$ a supercritical Hopf bifurcation, as previously reported by Keeler et al. (Reference Keeler, Thompson, Lemoult, Juel and Hazel2019). In the experiments, the first two bifurcations occur at values of $Q$ lower by factors of $1.9$ and $1.6$, while $Q_{H2}$ is larger than its numerical counterpart by a factor of $1.5$.

In addition, the numerical bifurcation diagram reveals the organisation of the unstable solution branches. We will discuss their influence on the bubble evolution in § 5.3.2, but we note the absence of any symmetric solution for $Q > 0.1345$, where the unstable symmetric branch S4 reaches a limit point.

5.3.2. Single-bubble edge states

The symmetric unstable solution branches associated with double-tipped (S1a and S1b) and triple-tipped (S2) bubbles, calculated numerically from the theoretical model and reported in figure 15, are strikingly reminiscent of the early-time bubble shapes transiently observed in experiments; see figure 7. We also observe that the order in which the double and triple-tipped shapes appear in the numerical solution space is consistent with the order in which the corresponding shapes are observed experimentally as a function of the flow rate, as also reported by Franco-Gómez et al. (Reference Franco-Gómez, Thompson, Hazel and Juel2018). We believe that broadly equivalent unstable steadily propagating states exist in the model and the experiment. We now show that these unstable states have a dramatic impact on the transient dynamics of the system as a whole, with the richness of the unstable solution structure reflected in the transient dynamics.

In the model, the first pair of solution branches exists only for very low flow rates, merging at a limit point at $Q_F=0.005085$. The S1b and S1a branches have respectively one and two unstable eigenvalues and only the shape associated with the S1a branch features two pronounced tips. In the experiments, shapes similar to the S1b and S1a solutions are transiently explored at $Q=0.005$ in figures 5(a) ($t=2.4$) and 6(a) ($t=2.3$), respectively, before the bubble eventually migrates towards a stable AS1 solution. Hence it appears that slender bubbles are initially attracted towards an equivalent of the S1b solution via its symmetric stable manifold and are finally repelled via its asymmetric unstable manifold, while wider bubbles are attracted to the symmetric S1a solution before being repelled asymmetrically. In this scenario, the S1a and S1b solutions are edge states within a basin boundary that separates the two possible asymmetric solutions corresponding to the bubble being biased towards one side or the other of the channel.

We can test this claim by studying a numerical time simulation at a low flow rate $Q=0.005 < Q_F$. Figure 16 shows the evolution of a single bubble from initial elliptical shapes of width $d = 0.4$ (a) and $0.8$ (b). The bubble contours corresponding to the S1b and S1a unstable solutions are shown with coloured dashed lines in figures 16(a) and 16(b), respectively. At $Q=0.005$ the S1a solution has two unstable eigenvalues $\lambda = 1.43560$ and $\lambda =0.27542$ (the time dependence is proportional to $\exp (\lambda t)$), whilst the S1b solution has a single unstable eigenvalue $\lambda =1.09033$.

Figure 16. Numerical time evolution of initially elliptical bubbles of width $d=0.4$ (a) and $d=0.8$ (b) and size $r = 0.54$. The flow rate $Q=0.005$ is just below the value at which the S1a and S1b branches disappear through a limit point. The dashed bubble contours correspond to the S1b (a) and S1a (b) unstable steady states; see figure 15. In both cases, the bubble eventually reaches the stable AS1 solution, although in (a) it converges to the ‘upper’ asymmetric branch and in (b) to the ‘lower’ asymmetric branch. Trajectories of the system are shown in (c) and (d) in a projection of the phase space $(p_b, U_b)$. In (d), the dashed line reproduces the trajectory shown in (c). Some of the times indicated in the snapshots (a) and (b) are reported on the trajectories in (c) and (d), respectively. Coloured crosses correspond to the positions of the S1a and S1b steady states. Note that in (c) the trajectory appears to explore AS1 twice because the phase-space projection does not capture offset variations, but in fact the bubble only reaches AS1 at late times.

In figure 16(c,d) we visualise the time evolution of the bubbles shown in figure 16(a,b) as trajectories projected into a two-dimensional phase plane spanned by bubble pressure $p_b$ and bubble velocity $U_b$, where the steady solutions S1a and S1b are indicated with coloured crosses. Figure 16(c) reveals that initially slender bubbles ($d=0.4$) visit the symmetric stable manifold of the least unstable (S1b) solution before leaving the vicinity of S1b via its asymmetric unstable manifold and eventually settling on a stable asymmetric AS1 solution (see inset in (c)). Initially wider bubbles ($d=0.8$) in figure 16(d) are instead attracted to the more unstable (S1a) solution resulting in a more visible two-tipped deformation. After leaving the vicinity of the S1a solution, the system again settles on an AS1 solution but in this case the bubble is nearer the ‘lower’ rather than the ‘upper’ wall. As shown in figure 16(a,b), the shape of the evolving bubble and that of the associated unstable steady state are virtually indistinguishable at times $t = 5.0$ and $t = 1.0$, respectively. Note, however, that the chosen phase plane projection does not capture a complete picture because it does not, for example, convey information about the bubble offset.

The highly deformed double-tipped shape is also transiently observed in the numerical calculations for larger values of $Q$, for which the S1 branch does not exist. Figure 13 shows that for $0.010 \le Q \le 0.025$ the two-tipped shape is seen for all bubble widths and also for $Q \ge 0.025$ when the bubble is sufficiently wide. Although the S1 branch does not exist in this range, the trajectories in phase space are expected to vary smoothly beyond the limit point at $Q_F$. Hence it is plausible that the system will trace a trajectory through phase space that results in two-tipped bubble shapes even if the steady solution is not present. Figure 17 shows the trajectories of slender and wide bubbles at $Q=0.006$ where the S1 branch does not exist. For both cases the trajectory initially follows a qualitatively similar path through phase space compared to when $Q=0.005$. However, the late evolution is different in that rather than reaching the stable AS1 state, the bubble breaks up. It appears that the unstable S1 branch plays a vital role in enabling a single bubble to approach the AS1 state as its absence means that the bubble will break up. In other words, the edge states on the S1 branch organise the dynamics so that breakup is avoided. Therefore the value of $Q_F$, the highest flow rate for which the unstable branch exists, marks the threshold between the solid and hollow circular markers in figure 13.

Figure 17. Numerical time evolution of initially elliptical bubbles of width $d=0.4$ (a) and $d=0.8$ (b) and size $r = 0.54$. The flow rate $Q=0.006$ is just above the value at which the limit point $F$ occurs so that the S1 branch does not exist; see figure 15. (c,d) Trajectories are shown in a ($p_b$,$U_b$) projection of the phase space with a layout similar to that of figure 16.

The same approach can be used to examine the influence of the unstable S2 solution on the dynamics of the system. Figure 18 shows trajectories for a slender and wide bubble when $Q=0.029$. At this flow rate, in addition to the stable AS1 solution, there are three unstable steady states on the S2 branch. The most unstable S2a solution has four unstable eigenvalues $\lambda =1.29566$, $2.94015\pm 1.48472i$, 4.15816, the S2b solution has three unstable eigenvalues $\lambda =1.44914\pm 2.51394i$, 1.85719 and the S2c solution has a complex conjugate pair of unstable eigenvalues $\lambda =0.19702\pm 3.54153i$. For an initially slender bubble ($d=0.4$) the system evolves towards the least unstable S2c solution where it stays in the stable manifold for a long time; see figure 18(a,c). Due to the complex eigenvalues the trajectory oscillates in the vicinity of S2c (see the inset time evolution of $U_b$ in (c)) before eventually breaking up in a highly asymmetric manner that resembles a thwarted migration towards the stable AS1 solution. In contrast the dynamics for an initially wide bubble ($d=0.8$) is significantly quicker resulting in a symmetric breakup event. These results demonstrate not only that the unstable states are transiently observed during time evolution, as claimed, but also that their presence influences the dynamics of single bubbles.

Figure 18. Numerical time evolution of initially elliptical bubbles of width $d=0.4$ (a) and $d=0.8$ (b) and size $r = 0.54$ for $Q=0.029$. The dashed contours in (a) indicate the S1c unstable steady state; see figure 15. (c,d) Trajectories are shown in a ($p_b$,$U_b$) projection of the phase space with a layout similar to that of figure 16. Coloured crosses in (c) and (d) correspond to the positions of the S2a, S2b and S2c steady states; see figure 15. In both cases, the bubble eventually breaks up, but the trajectories are very different. The initially slender bubble lingers for a long time near the S2c solution, which is weakly unstable with oscillatory eigenmodes, before eventually breaking up asymmetrically. The initially wide bubble experiences a much more straightforward evolution towards breakup into two equal bubbles.

5.4. Two-bubble edge states

In both the experiments and numerical simulations at relatively low flow rates, an initially wide bubble splits in two and a critical flow rate $Q_{c}$ determines whether the two bubbles coalesce or separate indefinitely; see figures 7 and 13. We now provide evidence that after bubble breakup, unstable two-bubble steady states are so-called edge states responsible for the organisation of these dynamical outcomes of the system.

Numerical simulations of the time evolution of an initially elliptical bubble of width $d=0.72$ at $Q=0.02$ are shown in figure 19(a), demonstrating that the bubble splits into two bubbles which later coalesce. Simulations of the evolution of the same initial bubble at a higher flow rate, $Q=0.04$, are shown in figure 19(b) and in this case the two bubbles do not coalesce. The computed time evolutions are qualitatively similar to the experimental evolutions shown in figures 6(b) and 6(c), respectively. Steady calculations reveal that there is an unstable steady ‘up–down’ (UD) two-bubble solution at $Q=0.02$, shown as blue dashed lines in figure 19(a). At $Q=0.04$, in addition to the UD solution there is another ‘barrier’ (B) unstable steady solution, shown as red dashed lines in figure 19(b).

Figure 19. Numerical time evolutions at flow rates $Q=0.02$ (a) and $Q=0.04$ (b) for an initially elliptical bubble of width $d=0.72$ and size $r=0.54$. The dimensionless time is given at the bottom of each snapshot. The two-bubble coloured dashed shapes correspond to unstable steady states of the two-bubble system, labelled ‘UD’ for ‘up–down’ (blue) and ‘B’ for ‘barrier’ (red). According to our depth-averaged model, the B solution emerges for flow rates greater than $0.029 \pm 0.001$. The position $(x_1,y_1)$ of the centroid of the trailing bubble is identified by a point. The associated trajectories of the two-bubble system in the $(x_1,y_1)$ phase-space projection are represented by the solid lines in (c) ($Q=0.02$) and (d) ($Q=0.04$). In (d), the dashed line is a reproduction of the trajectory in (c). Times corresponding to the snapshots of (a) and (b) are reported on the trajectories (c) and (d), respectively. Coloured crosses in (c) and (d) correspond to the values of $(x_1,y_1)$ associated with the ‘UD’ and ‘B’ unstable steady states. See also supplementary movie 3.

In order to demonstrate how the dynamics is related to the identified steady states of the two-bubble system, figures 19(c) and 19(d) show the time evolutions in a projection of the phase space of the system created by plotting the position $(x_1,y_1)$ of the centroid of the trailing bubble when two bubbles are present. Note that the projected phase space is not the same as that in § 5.3.2. The steady two-bubble solutions are also marked on the same projection of the phase space. The two bubbles are formed at the point labelled ‘breakup’ and the subsequent evolution is shown by the solid black line. In both cases, the two bubbles initially propagate side-by-side as the system is attracted towards the UD steady state, which has only a single unstable eigenvalue and is therefore strongly attracting. As seen in figures 19(a) and 19(b) the bubble shapes are virtually indistinguishable from those of the UD steady state at $t=1.8$ and $t=1.7$, respectively. The system then moves away from the unstable steady state in the direction associated with the single unstable eigenvalue: in the projection of the phase space, the unstable steady state is analogous to a two-dimensional saddle point.

At $Q=0.02$, the lower bubble moves up and over the rail and two bubbles aggregate and coalesce. In terms of movement through the phase space, the system is attracted to a compound bubble that is analogous to a stable single-bubble state. Examples of such compound bubbles are found throughout § 4. At $Q=0.04$ the presence of the second unstable steady state and its associated unstable manifold prevents the system from reaching a coalescence point by deflecting the trajectory. Although the lower bubble starts to move over the rail, it remains on the lower side of the rail because the B state is sufficiently repulsive. Unable to reach the coalescence point, the two bubbles continue to separate because they propagate at different speeds. The B steady state arises through a saddle-node bifurcation and only exists for Q greater than a threshold flow rate $0.029 \pm 0.001$ which can hence be identified with the value $Q_{c(num)}$ roughly estimated in figure 13.

These calculations provide compelling evidence that the unstable two-bubble steady states are edge states whose stable and unstable manifolds partition the dynamical outcomes of the system. The two-bubble system is more complicated than the single-bubble case, in part because the volume ratio between the two bubbles is an additional parameter. For example, if the bubble is not initially symmetric, then breakup may result in two bubbles of significantly different sizes for which the impact of the edge states described here may differ from figure 19. For all size ratios we have examined, we find that every single-bubble steady state has an equivalent two-bubble steady state with near identical flow rate–propagation speed characteristics and hence these solution branches are qualitatively unchanged. There are, however, additional two-bubble steady states that do not have single-bubble equivalents and a thorough investigation and discussion of two-bubble steady states is provided in the companion paper (authors' unpublished observations).