The life and fate of a bubble in a geometrically perturbed Hele-Shaw channel

Abstract Motivated by the desire to understand complex transient behaviour in fluid flows, we study the dynamics of an air bubble driven by the steady motion of a suspending viscous fluid within a Hele-Shaw channel with a centred depth perturbation. Using both experiments and numerical simulations of a depth-averaged model, we investigate the evolution of an initially centred bubble of prescribed volume as a function of flow rate and initial shape. The experiments exhibit a rich variety of organised transient dynamics, involving bubble breakup as well as aggregation and coalescence of interacting neighbouring bubbles. The long-term outcome is either a single bubble or multiple separating bubbles, positioned along the channel in order of increasing velocity. Up to moderate flow rates, the life and fate of the bubble are reproducible and can be categorised by a small number of characteristic behaviours that occur in simply connected regions of the parameter plane. Increasing the flow rate leads to less reproducible time evolutions with increasing sensitivity to initial conditions and perturbations in the channel. Time-dependent numerical simulations that allow for breakup and coalescence are found to reproduce most of the dynamical behaviour observed experimentally, including enhanced sensitivity at high flow rate. An unusual feature of this system is that the set of steady and periodic solutions can change during temporal evolution because both the number of bubbles and their size distribution evolve due to breakup and coalescence events. Calculation of stable and unstable solutions in the single- and two-bubble cases reveals that the transient dynamics is orchestrated by weakly unstable solutions of the system that can appear and disappear as the number of bubbles changes.


Introduction
Fundamental studies of complex fluid mechanical phenomena usually concentrate on the steady or periodic flows that are observed after some time. Experimental observation of steady or periodic flows is used to infer that they are stable solutions of the underlying governing equations. Unstable solutions of the governing equations, on the other hand, are often ignored being thought to be unimportant and unobservable in experiments. In fact, unstable solutions can be glimpsed fleetingly in experiments as they influence the transient evolution of a flow, which is often as important as the final state.
In this paper, we conduct detailed experiments to characterise the transient evolution and ultimate disposition of a simple fluid mechanical system that nonetheless exhibits complex nonlinear behaviour: a finite air bubble propagating within a liquid-filled Hele-Shaw channel containing a relatively small depth perturbation; see figure 1. The stability of finite steady bubbles in unperturbed channels has previously been addressed both experimentally (KopfSill & Homsy 1988) and numerically (Tanveer & Saffman 1987). However, interfacial flows in narrow gaps can also exhibit considerable disorder (Saffman & Taylor 1958;Couder 2000;Vaquero-Stainer et al. 2019), but they are rarely investigated from a dynamical systems perspective. Although our primary interest is to use the system as a model in which to study complex transient dynamics, the flow of bubbles through channels also has a number of industrial, chemical and biological applications, particularly in microfluidic geometries (Anna 2016).
The transient behaviour of a flow is particularly significant when the flow is disordered with no apparent stable steady or periodic states. The connection between disorder in fluid mechanics and the nascent theory of dynamical systems was initiated via lowdimensional mathematical models of turbulence (Hopf 1948) and atmospheric convection (Lorenz 1963). Subsequently, Ruelle & Takens (1971) identified turbulence with deterministic chaos. Characterised by long-term unpredictability due to sensitivity to initial conditions in systems of three or more degrees of freedom, chaos was shown to emerge after only a small number of bifurcations. The first observation of chaos in a fluid system was made by Gollub & Swinney (1975) in Taylor-Couette flow, where geometric confinement reduces the spatial complexity of the flow (Mullin 1993), so that its temporal dynamics can be captured by models of reduced dimensionality tracking a small number of dominant eigenmodes that control the dynamics (Manneville 1990). The relaxation of spatial constraints by increasing the aspect ratio of such closed flows can lead to transitions to spatio-temporal chaos (Bodenschatz et al. 2000;Gollub 1995). However, theoretical insight into the dynamics of these systems tends to be reduced because of the increased number of significant eigenmodes. In recent years, interest in dynamical systems ideas in fluid dynamics has been rekindled through their considerable success in elucidating the subcritical transition to turbulence in shear flows (Kerswell 2005;Schneider et al. 2007;Eckhardt et al. 2008;Duguet et al. 2008). Rather than focusing on the eigenmodes of these weakly-confined, open flows, the approach has been to reduce the dynamics to a small number of invariant objects in phase space.
It is now widely accepted that stable and unstable invariant objects play a key role in organising the dynamics of transition to turbulence (Kawahara et al. 2012;Cvitanović 2013). In this context, an invariant object is a subset of the possible states of the system that does not change under time evolution. For example, a steady solution is invariant because the system does not move away from it unless perturbed. A periodic solution is also invariant because the system evolves through a limited set of states over one period and then the cycle repeats. More complex invariant objects include the stable and unstable manifolds associated with steady and periodic solutions and chaotic attractors. The fact that the system cannot leave an invariant state means that they divide the possible states of the system into regions with different dynamical outcomes; in other words, trajectories of the system cannot cross an invariant object. For the turbulent transition problem, the most important invariant is the so-called 'edge' of turbulence dividing regions that eventually return to laminar flow from those that progress to a fully turbulent state.
A simple conceptual picture of an edge state is to consider the two-dimensional steady stagnation point flow of a fluid. The stagnation point itself is an unstable fixed point of the flow because fluid particles will move away if displaced from the stagnation point. As well as the stagnation point, the important invariant objects are the stable and unstable manifolds, which are simply the dividing streamlines of the flow. These manifolds divide the flow into four different regions and a particle that starts in one region cannot move into another. Hence a complete characterisation of the ultimate fate of a fluid particle follows from knowing its position relative to the dividing streamlines. In other words, the manifolds are the edge states of this simple system.
A difficulty in applying these ideas in fluid dynamics is that the invariant objects of interest can be extremely high dimensional and cannot be easily determined. Nevertheless, recent advances in both computer power and efficient algorithms (Tuckerman et al. 2019;Farano et al. 2019) have allowed ever increasing access to such objects in studies of turbulence (Budanur et al. 2017), convection (Sánchez Umbría & Net 2019) and droplet break up (Gallino et al. 2018). In our chosen system, the bubble exhibits a variety of complex behaviours including symmetry breaking, bistability, steady and periodic states as well as non-trivial transient evolution (Franco-Gómez et al. 2018). The Reynolds number remains small and the nonlinearities arise only through the presence of the air-liquid interface. Hence, direct observation of the interface can be used to distinguish different possible states of the system. We develop a new experimental protocol to generate a variety of different initial bubble configurations and then investigate the bubble's subsequent behaviour under imposition of flow. We find a rich variety of transient evolutions that can, nonetheless, be rationalised and classified because they occur in distinct regions of the parameter space. In addition the system is attractive from a theoretical point of view because its behaviour can be captured both qualitatively and quantitatively by a depth-averaged set of equations (Franco-Gómez et al. 2016), provided that the aspect ratio of the channel's cross-section is sufficiently large. Keeler et al. (2019) presented a purely theoretical study of the depth-averaged model in which unstable periodic orbits were found to be edge states marking the boundary between steady propagation along either the centreline (symmetric) or biased towards one side of the channel (asymmetric). In the present paper, we synthesise experimental data and numerical solutions of the depth-averaged equations to examine the influence of a wide variety of stable and unstable invariant objects on the behaviour of the bubble.
A distinct feature of the system is the propensity of the bubble to break up into two or more smaller bubbles that may or may not recombine in the subsequent dynamics. The transition from single to multiple bubbles allows the evolution of the invariant-object structure for fixed parameter values, which is rather unusual. The resulting complex dynamics are nonetheless orchestrated by a variety of edge states that are invariant solutions of the depth-averaged system.
The paper is organised as follows. The first three experimental sections present the experimental set-up and protocol ( §2), the method used to set the initial bubble shape ( §3) and the experimental results ( §4). This is followed by a presentation of the depth-averaged model and numerical methods ( §5.1) and by an interpretation of the experimental results in terms of dynamical system arguments ( § §5.2-5.4). We draw our conclusions in §6.

Experimental setup
A schematic diagram of the experimental setup is shown in figure 1. The Hele-Shaw channel used in the experiments consisted of two parallel float glass plates of dimensions 170 cm × 10 cm × 2 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 ± 0.01 mm and its width was W * = 40.0 mm with variations up to ±0.1 mm over the length of the channel L * = 170 cm, yielding an aspect ratio α = W * /H * = 40. A thin depth perturbation of rectangular crosssection, width w * = 10.0 ± 0.1 mm and thickness h * = 24 ± 1 µ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 • .
The channel was filled with silicone oil (Basildon Chemicals Ltd) of dynamic viscosity µ = 0.019 Pa.s, density ρ = 951 kg/m 3 and surface tension σ = 21 mN/m at the laboratory temperature of 21 ± 1 • C. 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 1(c)), 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 * = A * /π 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 = 2r * /W * , y c = 2y * c /W * and U b = U * b /U * 0 , respectively. We define a non-dimensional flow rate Q = µU * 0 /σ 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. The capillary number based on the bubble velocity, which measures the ratio of viscous to capillary forces, is Ca = µU * b /σ with values in the range 0.008 Ca 0.48. The ratio of inertial to viscous forces yields a reduced Reynolds number Re = ρU * 0 H * 2 /µW * which takes values less than 0.3 in our experiments, thus indicating negligible inertia. The Bond number Bo = ρgH * 2 /4σ = 0.11, where g is the gravitational acceleration, indicates that bubble buoyancy is also negligible.

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  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 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 about 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 is 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). time scale of the experiment. The thickness of these films increases with the capillary number for a propagating bubble (Bretherton 1961), 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 Ca = 0.48 investigated in this paper, the film thickness averaged over the bubble area is h * f /H * ≈ 0.2 (i.e. h * f ≈ 200 µ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 ± 1 µ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 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 break up 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 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 steady states of the system. 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 × 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. (2017).
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µ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µQ * /(πr * 4 t ) 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 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 to 477 mL/min. This suggests that the solubility of air in oil increases with pressure, as previously reported by Franco-Gómez et al. (2017). 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 parametrized our bubble propagation experiments according to the in-flow equivalent radius of the bubble r Q = (2 A * Q /π)/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.
3. Relaxation of a bubble at Q = 0: selection of the initial bubble shape The relaxation of a centred, elongated bubble in the absence of imposed flow provides a non-invasive method to vary the shape of a bubble while conserving its volume. We selected specific shapes reproducibly as initial conditions for the propagation experiments described in §4 by controlling the time for which the bubble relaxes from its original shape before imposing the flow.
The elongated bubble was prepared by flowing a newly formed bubble through the . The initial slender symmetric shape corresponds to the stable symmetric propagation mode at Q * i = 78 ml/min which was reached before the flow was interrupted at t * = t * i . The time steps between two successive images are ∆t * = 3.2 s (left) and ∆t * = 26.7 s (right). (c) Time evolution of the non-dimensional bubble width d = 2d * /W * during relaxation from di at t * = t * i to a maximum value dm for different bubble sizes. Inset: di and dm plotted as a function of rs. Empirical fits: dm = 1.15rs if rs < 0.87 and dm = 1 if rs > 0.87; di = 0.54 − 1.05 exp(−4.35rs). (d) Variation of dm − d as a function of t * − t * i plotted on a semi-log scale. The initial regime is fitted by an exponential from which we measure a characteristic relaxation time τ * . Inset: τ * /(r * s µ/σ) against confinement ratio r * s /H * fitted with a power law of the form y = 0.63x 2.7 . centring device described in §2.1. A flow rate of Q * i = 78 ml/min (Q i = 0.03) was chosen so that upon exit, the bubble adopted a steadily-propagating elongated shape, which straddled the rail symmetrically about the centreline of the channel (y c = 0); see §4.3.2 for a summary of the steady states of the system. Once this symmetric steady state had been reached a few centimetres downstream of the centring device, the flow was interrupted (at time t * i ) and the bubble was left to relax its shape at Q * = 0. This ensured that the geometry of the centring device did not influence the evolution of the bubble. Note that the interruption of the flow led to the increase of the projected area of the bubble to its static value A * s on a time scale much shorter than that of the bubble relaxation process.
The bubble relaxation process is shown in figure 3. We characterise the shape of the bubble by its maximum length 2l * and maximum width 2d * aligned along x * and y * , respectively. These were measured by identifying the smallest rectangle bounding the bubble contour; see figure 3(a). They are related through the range of shapes adopted by the bubble during relaxation and thus, we choose to parametrize the bubble shape solely by its non-dimensional width d = 2d * /W * .
Typical shape evolutions are shown in figure 3(b), where images taken at a constant time interval are superposed for bubbles of static size r s = (2 A * s /π)/W * = 0.54 (left) and r s = 1.24 (right). Starting from an elongated initial bubble shape at t * i , the length of the bubble decreases and its non-dimensional width d increases from an initial value d i at time t * i to a plateau value d m , as shown in figure 3(c) for the range of bubbles sizes investigated. During relaxation, the bubble area decreased by ∼ 2 % because the fraction of the bubble situated over the rail decreased. Values of d i and d m are plotted against r s in the inset of figure 3(c), which indicates that the bubble relaxes to the width of the channel for r s > 0.87. In the bubble propagation experiments in §4, we are limited to initial bubble widths ranging between d i and d m for a given bubble size and to the shapes explored by the bubble during relaxation. We used the curves shown in figure  3(c) to infer the time required to reach a desired value of d from which to initiate the propagation of the bubble.
The time evolution of d m −d is plotted on a semi-log scale in figure 3(d), revealing that the relaxation of small bubbles (r s 0.46) is approximately exponential. Larger bubbles (e.g. r s = 1.24) exhibit qualitatively different relaxation, which is only exponential at early times, presumably because of the influence of the side walls of the channel. The exponential part of these relaxation curves was fitted to obtain a characteristic relaxation time scale τ * , which saturates to approximately 35 s for bubbles of size r s > 1. The nondimensional decay time τ * σ/(r * s µ) is plotted against the bubble confinement ratio r * s /H * in the inset of figure 3(d). It follows a power law of exponent 2.7 for bubbles smaller than the width of the channel (r s < 1), which is larger than the exponent 2.0 found by Brun et al. (2013) for the exponential relaxation of elliptical bubbles towards a circular shape in a liquid-filled, infinite Hele-Shaw cell of uniform depth, in the absence of external flow. A possible explanation is that the relaxation process is accelerated by the presence of the rail which supports pressure gradients driving the bubble into the deeper parts of the channel.
Due to the presence of the rail, small bubbles systematically migrated towards one of the deeper parts of the channel and reached an asymmetric equilibrium state, consistent with the results of Franco-Gómez et al. (2018) for r 0.87, whereas sufficiently large bubbles (d m ∼ 1) tended to a symmetric equilibrium state. This migration was associated with a decrease in static bubble area of up to 7 % as the fraction of bubble volume residing in the deeper parts of the channel increased. Small bubbles always migrated to the same side of the rail because of a small unavoidable bias in the levelling of channel. When this bias was minimised, the migration time scale was significantly longer than the relaxation time scale (typically one minute compared to 10 s for r s = 0.54), which ensured that the bubble had a centroid offset of |y c | 0.005 when the flow was initiated.

Transient evolution of a centred bubble propagating from rest
Having established how to reproducibly prepare initially centred bubbles of different shapes, we propagate each bubble along the channel by imposing a constant flow rate. We investigate the evolution of these bubbles as a function of initial bubble width and flow rate. Unless otherwise specified, we focus on bubbles sizes r Q = 0.54 measured during propagation.

Influence of the initial bubble shape
The dimensional time t * is indicated in the last snapshot of each time sequence. The in-flow bubble size is rQ = 0.54 (after flow-induced compression) and d is measured at t = 0 (just before flow initiation). 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).
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 4(a-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 4(d-e)). For Q = 0.07, the single bubble straddling the rail at late times in figure 4(a-c) is 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 no 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 break up and aggregation events. The wide range of transient evolutions obtained for relatively modest changes in d (55.5% d 63.7%) is fully reproducible. In some cases, the final state features compound bubbles resulting from the aggregation of two or more bubbles (figure 4(c, 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 will 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 4(a)) and widest (figure 4(f)) 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 bubbles 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. Each point on the bubble interface is displaced along its normal direction and thus, 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 spans a small enough angle from 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 (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 ± 0.1)% b (57.5 ± 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 − 2.6 in figure 4(b-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 the middle tip, resulting in a break up 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 break up. In figure 4(d), this leads to a break up 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).

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. Figure 5 shows the evolution of initially slender bubbles (within a narrow range of widths 53 % d 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 Q 0.18. In order to maintain the in-flow 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.

Initially slender bubbles
For all values of Q, the long-term outcome is a single (simple or compound) bubble which is stable and propagates steadily. 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 5(a, b)), a bubble which straddles the rail symmetrically about the channel centreline, referred to as 'on-rail, symmetric' (final frame in figure 5(c)) and 'on-rail, asymmetric' bubbles (final frames in figure 5(d, 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. 2012;Jisiou et al. 2014;Thompson et al. 2014). 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 5(f)). 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 edge states of the one-bubble system in §5.3.2. For the smallest values of Q ( figure 5(a, 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, whereas for Q = 0.01, 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 break up 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 break up and aggregation events following tip-splitting beyond a threshold value of the flow rate Q ts = 0.085 ± 0.005. This is illustrated in figure 5(e, f), where the bubble breaks into two parts of widely 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. Figure 6 shows the evolution of initially wide bubbles (d = 63.4 %) as a function of the flow rate in the range 0.005 Q 0.18 with a layout similar to figure 5. For low flow rates (figure 6(a, 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.

Initially wide bubbles
However, for higher flow rates (figure 6(c, d)), the long-term outcome is two off-rail bubbles, where the leading bubble is larger 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 6(b)) and separation (figure 6(c)) occurs at a threshold value of the flow rate Q c = 0.017 ± 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 break up of the smaller trailing bubble, as shown in figure 6(d) where tiny bubbles created through break up aggregate with the bubble that spawned them. These secondary break up events occur for Q > 0.135 ± 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.

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.

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.
For Q 0.025±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±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 offrail 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 ± 0.001. For Q < Q c , the single off-rail bubble is reached by migrating off the rail without breaking (•) or by breaking into two parts which subsequently aggregate (•). By contrast, the two bubbles separate after break up for Q c < Q 0.25 ± 0.005 ( ).
For Q > 0.025 ± 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 d d 3 the middle tip grows faster than the side tips whereas the converse is true for d 3 d 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 ( and ) 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 ( ) to six ( )). 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 side walls decreases from six ( ) to four ( )). 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. (2017 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).

Single-bubble invariant modes of propagation
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 P F = 0.0435, obtained from a linear fit of centroid offset values y 2 c (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 are plotted as functions of the non-dimensional flow rate Q for three bubble sizes rQ = 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 yc and U b are averaged over one period.
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).

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 Q ts ; yellow region and hatched region in figure 7).

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 ± 0.001, which marks the transition between aggregation (figure 6(b)) and separation (figure 6(c)). We found that the final outcome of the system is influenced by the value of the size ratio of the two bubbles formed after break up, which depends on y c . For example, at Q = 0.020, the two bubbles resulting from break up 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 increasing 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 that evolutions where break up events had a propensity to generate tiny bubbles featured reduced experimental reproducibility. 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 9(c)) is caused by differences in relative bubble sizes after break up, stemming from small unavoidable differences in the initial bubble shape of the order of 0.5 % in area (t = 0), which is the typical reproducibility limit of the experimental setup.

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 ± 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.
The early-time symmetric evolution, which here involves break up and aggregation, leads to an approximately symmetric compound bubble by t = 4.5. In both figures 10(a) and (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 break up which produces two bubbles of similar sizes in figure 10(b) and widely 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 steadilypropagating 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 break up and aggregation events (hatched region in figure 7). This is illustrated in figure 11 where five successive experiments performed for parameter values corresponding to the yellow 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 figures 11(a, b) result in a single on-rail bubble (t = 17.0) and to two offrail bubbles (t = 13.0) as 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, e) provide examples of more complex time evolutions after centred tip-splitting, involving several break up 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 break up 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), tipsplitting 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 ± 0.005 and becomes increasingly likely the larger the flow rate.
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 et al. (1987) 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.

Comparison between experiment and numerical model
Interpretation of the phenomena described in §4 based on the dynamical systems concepts presented in the introduction requires knowledge of the underlying stable and unstable invariant objects. In this section we present numerical results from a previously validated depth-averaged model that can be used to identify the appropriate invariant objects of the system.

Depth-averaged model for bubble propagation
The depth-averaged model extends McLean & Saffman's (1981) 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 Franco-Gómez et al. 2016) and a single closed air bubble (Franco-Gómez et al. 2017Keeler et al. 2019) 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 break up.
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 b(y) = 1 − 1 2 h [tanh(s(y + w)) − tanh(s(y − w))] , (5.1) 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. ρU * 0 W * /(µα 2 ) 1) and a large aspect ratio (i.e. α = W * /H * 1). Within the fluid domain, the 2D depth-averaged velocity u satisfies where the second equation is a statement of mass conservation. At the channel side walls, 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 at each bubble interface, where R is the 2D position of a point on the bubble boundary, n an outward normal unit vector, and 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 where p i is the internal pressure in bubble i, p the fluid pressure, κ 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: Here Γ i is the interior of bubble i and ∂Γ 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 2006). 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. Timestepping uses BDF1 or BDF2 (backwards difference formulae), with a new triangular mesh generated at least every 5 time steps. The channel length is typically taken as six times its width, to allow for interactions between bubbles of moderate separation.
The experimental results presented in §4 show frequent changes in bubble topology and we must therefore allow for both coalescence and break up in our numerical simulations. To the best of our knowledge, there is no numerical method for solving our model equations in which bubble break up 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.
As the region inside the bubble is treated simply as a region of constant pressure and fixed volume, the need for break up is signalled by self-intersection of the bubble boundary, which here happens in finite time (see figure 12(a, 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 12(c)). In contrast to break up, 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 12(d)), we remove the boundary nodes that are within the tolerance, and reconnect the remainder ( figure 12(e)). Our topological 'surgeries' for both coalescence and break up 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 12(f)). For coalescence, we set the target volume V i for the new bubble to be the sum of the two parent bubbles. For bubble break up, 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. 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 break up between model and experiments in our system indicates that break up 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 is 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  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 πld and size r = √ 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 ( ) 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 P F (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) Q < Q P F (num) , asymmetric on-rail bubbles for Q P F (num) Q < Q H2(num) and oscillatory on-rail bubbles for Q Q H2(num) . simulations starting from a variety of initially elliptical bubble shapes. Specifically, we define the initial bubble shape as 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 = πr 2 = πld with r = 0.54, with the bubble volume calculated from (5.4). 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. Finally steady solutions can be calculated by setting the time derivatives to zero and the effects of the flow rate and volume can be mapped out using continuation and edgetracking techniques (Skufca et al. 2006;Gelfgat 2019).

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 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.
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 twobubble long-term outcomes occurs at a threshold flow rate Q c(num) = 0.025 ± 0.005 (compared to the experimental value Q c = 0.017 ± 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) Q < Q P F (num) ); asymmetric on-rail bubbles (Q P F (num) Q < Q H2(num) ); and finally oscillatory on-rail bubbles (Q 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 and , 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.
A distinction from the experimental phase diagram is the increased range of initial bubble widths 0.40 d 0.80 in figure 13 compared to the experimental range in figure 7. For d 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 ( , 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 ± 0.005 is close to the value Q ts = 0.085 ± 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 twobubble 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 ( ) 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 break up as the flow rate increases. In particular, modest changes in relative bubble volumes after break up 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.

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 plotted as a function of the flow rate Q for a fixed bubble volume V = 0.54 2 π. 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.
The numerical simulations indicate that Q H1(num) = 0.0301 marks a subcritical Hopf bifurcation, Q P F (num) = 0.0714 a supercritical pitchfork bifurcation and Q H2(num) = 0.0789 a supercritical Hopf bifurcation, as previously reported by Keeler et al. (2019). 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.

Single-bubble edge states
The symmetric unstable solution branches associated with double-tipped (S1a and S1b) and triple-tipped (S2) bubbles, calculated numerically for the theoretical model 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. 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. (2018). We believe that broadly equivalent unstable steady 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 the 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. We can test this claim by studying a numerical time simulation at a low flow rate (Q = 0.005, with the limit point at Q F = 0.005085). 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 (b), respectively. At Q = 0.005 the S1a solution has two unstable eigenvalues λ = 1.43560 and λ = 0.27542 (the time dependence is proportional to exp(λt)), whilst the S1b solution has a single unstable eigenvalue λ = 1.09033.
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 . 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 break up into two equal bubbles.
vicinity of S1b via its asymmetric unstable manifold and eventually settling on the 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 twotipped deformation. After leaving the vicinity of the S1a solution the system settles again on an AS1 solution. 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 Q 0.025 the two-tipped shape is seen for all bubble widths and also for Q 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 directly the bubble breaks up. It appears 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. 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.
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 λ = 1.29566, 2.94015 ± 1.48472i, 4.15816, the S2b solution has three unstable eigenvalues λ = 1.44914 ± 2.51394i, 1.85719 and the S2c solution has a complex conjugate pair of unstable eigenvalues λ = 0.19702 ± 3.54153i. For an initially slender bubble (d = 0.4) the system evolves towards the 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 signal 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) are significantly quicker resulting in a symmetric break up 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.

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 break up 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 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 (c) respectively. Steady simulations 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).
In order to demonstrate how the dynamics is related to the identified steady states of the two-bubble system, figures 19(c) and (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 'break up' 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 (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 move 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 ± 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 break up 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 (Keeler et al. 2020).

Conclusion and perspectives
We have studied the dynamics of a bubble evolving from a prescribed initial condition as it propagates through a Hele-Shaw channel driven by the steady motion of a suspending viscous fluid. We have focused on a channel with a depth perturbation in the form of a centred rail, because it supports multiple invariant modes of bubble propagation which lay the foundation for complex dynamics. These invariant modes can be distinguished using only the interfacial morphology because the nonlinearity in this system stems from the presence of the air-oil interface. We observe a rich array of organised transient behaviour, which may involve repeated bubble break up and coalescence events, before a long-term outcome is reached in the form of either a single-bubble invariant mode of propagation or two (or more) steadily separating bubbles (sorted in order of decreasing size and velocity). A highly unusual feature of this system is that the underlying invariant-solution structure evolves dynamically, because changes in bubble topology causes the number of bubbles to vary.
A depth-averaged model of bubble propagation captures the majority of the experimental dynamics and features similar single-bubble, stable invariant states as observed in the experiment. We therefore conjecture that broadly similar unstable invariant states occur in the experiment as in the model and we find that the bubble transiently adopts the shape of unstable solutions during its evolution. Moreover, we show that multiple unstable invariant states play a key role in organising the transient dynamics as edge states, by partitioning the dynamical outcomes of the bubble evolution via their stable and unstable manifolds. The multiplicity of edge states in turn underpins the subtle dependence of the dynamics on the width of the initial bubble. Initial conditions determine the specific paths of exploration between these edge states and thereby the time scales of transient evolution which can vary significantly.
We identify the specific two-bubble edge state that controls the fate of the bubble once it has broken up: in its absence the smaller trailing bubble migrates across the rail to coalesce with the leading bubble, while in its presence the products of the break up remain on opposite sides of the rail and separate indefinitely. The complexity of the two-bubble system is such that we defer its detailed analysis to a forthcoming paper (Keeler et al. 2020). A particularly striking feature to be discussed is that all single-bubble invariant states have equivalents in the two-bubble system. However, additional multiple-bubble states mean that the complexity of the dynamics increases with the number of bubbles.
For moderate flow rates, the evolution of the bubble is fully reproducible in the experiment. However, as the flow rate increases, the experiments exhibit enhanced sensitivity to both initial bubble condition and unavoidable rail perturbations. Hence, multiple repetitions of the same experiment lead to widely different long-term outcomes, e.g., multiple separating bubbles versus a single-bubble invariant mode of propagation. An analogous sensitivity is also found in numerical time simulations for comparably high values of the flow rate. This sensitivity makes the long-term behaviour of the system practically unpredictable.
In both the experiments and simulations the evolution of the bubble at high flow rate follows a complex scenario of transient emulsification and coarsening through repeated cycles of break up and coalescence. The sensitivity arises from differences in the break up of the bubble, rather than direct sensitivity to initial conditions per se, and slight differences in the relative volumes of the daughter bubbles that are created can lead to very different subsequent dynamics. We conjecture that the increased complexity of the dynamics is associated with a growing number of invariant states: the number increases with increases in the flow rate and also with increases in the number of bubbles. An increase in flow rate also promotes break up into an increasing number of bubbles with a broad size distribution which means that the number of potential invariant states can increase very rapidly. The large number of invariant states may entangle into a complex edge that need not persist for all time, but underlies long transients of disordered dynamics for sufficiently large flow rates. Whether the system can eventually exhibit selfsustaining disordered dynamics, equivalent to a fully-developed turbulent state, remains an open question.
In order to estimate the thickness of the wetting films, we performed experiments using dyed silicone oil (APD fluorescent leak detection dye). The rail was removed for simplicity so that the bubble propagated steadily and symmetrically about the channel centreline. The total oil thickness 2h * f was estimated by measuring the intensity of the light transmitted through the channel from the LED light box to the top-view camera moving with the propagating bubble. Calibration of the light intensity was performed by measuring values for air and oil-filled channels, respectively. We also checked that the light intensity was proportional to the thickness of the attenuating medium by performing measurements through an oil-filled channel containing a transparent plastic wedge of fixed angle. In the following, we assume that the top and bottom oil films have the same thickness h * f . A typical distribution of wetting film thickness over the projected area of the bubble is presented in figure 20(a) for a bubble propagating at a flow rate Q * = 100 mL/min. The bubble size during propagation is r Q = 0.59. The film thickness distribution is symmetric about the channel centreline and uniform along x * for fixed values of y * . Wetting films imposition of the flow increased with Q, exhibit behaviour representative of the system that is unlikely to differ from that of a bubble of fixed volume.
We also investigated the effect of flow rate for r Q = 0.8 and r Q = 1.2. At low flow rates (Q < Q c , see figure 7), time evolutions are similar to that of bubbles with r Q = 0.54. For Q > Q c , the main differences with bubbles of size r Q = 0.54 are for time evolutions leading to break up into two bubbles of comparable sizes. For Q 0.025, this happens for d > d 3 or for d < d 3 and Q > Q ts in the case of centred tip-splitting (see figure 7). While at r Q = 0.54 the two bubbles remained on opposite sides of the rail (see figures 4(d-f) and 10(b)), the smaller trailing bubble often broke into two parts after elongation across the rail at r Q = 0.8, leading to a state with potentially three separating bubbles. In contrast, at r Q = 1.2, the smaller trailing bubble crossed the rail, leading to a single off-rail bubble. This is presumably due to the faster retraction of the rear part of the leading bubble when its size increases, which results in lower pressures at its rear, thus intensifying its attraction of the trailing bubble.
Experiments with air fingers (semi-infinite bubbles) yield similar modes of propagation to those shown in figure 8 and the same initial deformation of the finger tip as in figure  7. When initial deformation leads to two or three tips along the finger front, only the fastest ultimately propagates while the other(s) stop growing and retract slowly. This simplifies the evolution because break up and aggregation events do not usually occur.
Appendix C. Onset of oscillations Figure 22(a) shows a snapshot of an oscillating asymmetric bubble of size r Q = 1.2 propagating at flow rate Q = 0.17 (PO branch in figure 8). The dimples periodically generated at the bubble tip (where the bubble front meets the edge of the rail) remain at fixed positions x * in the laboratory frame of reference. Each generated dimple progressively flattens under the effect of capillary pressure gradients and, consequently, the bubble is slightly flatter at the rear than at the front. We define the oscillation amplitude A * H2 as the maximum peak-to-peak amplitude of the interface perturbation, see figure  22(a). The plateau value of the oscillation amplitude 2A * H2 /W * (reached following transient increase) is plotted in figure 22(b) as a function of the non-dimensional flow rate Q for three bubble sizes. The small amplitude values measured at low flow rates indicate the increasing excitability of the system near the transition which leads to damped oscillations of the bubble. The amplitude growth at higher flow rate is accurately captured by a square-root dependence on Q for the two highest bubbles sizes, which yields a critical flow rate Q H2 = 0.118.
The experimental oscillation frequency f * H2 reported in figure 22(b) increases linearly with flow rate and only slightly with bubble size. The associated non-dimensional pulsation 2πf * H2 ((W * /2)/U * 0 ) presented in the inset of figure 22(b) takes an approximate average value of 12.1. Our experimental resolution is not high enough to capture the expected pulsation increase near the bifurcation, which was calculated by Keeler et al. (2019).