## Abstract

Analytical solutions are constructed for an assembly of any finite number of bubbles in steady motion in a Hele-Shaw channel. The solutions are given in the form of a conformal mapping from a bounded multiply connected circular domain to the flow region exterior to the bubbles. The mapping is written as the sum of two analytic functions—corresponding to the complex potentials in the laboratory and co-moving frames—that map the circular domain onto respective degenerate polygonal domains. These functions are obtained using the generalized Schwarz–Christoffel formula for multiply connected domains in terms of the Schottky–Klein prime function. Our solutions are very general in that no symmetry assumption concerning the geometrical disposition of the bubbles is made. Several examples for various bubble configurations are discussed.

## 1. Introduction

Many free boundary problems naturally arise from the consideration of different types of Hele-Shaw systems. A Hele-Shaw system is one where two viscous fluids (typically with one fluid much less viscous than the other) are sandwiched between two closely spaced parallel plates so as to produce a flow that is essentially two dimensional. Hele-Shaw flows in various geometries have been extensively studied over the years, and many processes in physics involving the evolution of interfacial boundaries, such as dendritic crystal growth, direct solidification and fluid displacement in porous media, can be modelled mathematically (under certain assumptions) as a free boundary problem of the Hele-Shaw type. This diverse array of free boundary problems has a plethora of analytical solutions and a wide range of mathematical methods can be used to solve them [1,2]. The models defining these free boundary problems also go by the name of Laplacian growth processes [3] because the governing field equation in the viscous fluid region is Laplace's equation and the evolution of the fluid interfaces is governed through surface derivatives of this field. In the case of Hele-Shaw bubbles, the flow is governed by Darcy's law and the bubble interfaces evolve with a velocity proportional to minus the local gradient of the fluid pressure.

In this paper, we present a general analytical solution for the problem of multiple bubbles steadily translating along a Hele-Shaw channel when surface tension effects on the bubble boundaries are neglected. The mathematical problem to be solved belongs to a class of free boundary problem defined in a multiply connected domain which is related to a special class of scalar Riemann–Hilbert problem [4] recently considered by Crowdy [5]. Here, however, we shall use a more direct approach which allows us to reduce our original free boundary problem to two fixed boundary problems that can be more easily solved. More specifically, we introduce a conformal mapping from a circular domain in an auxiliary complex *ζ*-plane to the physical flow domain (i.e. the exterior of the bubbles bounded by the channel walls) in the complex *z*-plane and show that this mapping can be written as the sum of two analytic functions, corresponding to the complex potentials in the laboratory frame and in the co-moving frame, respectively. Each one of these functions maps the circular domain onto a multiply connected degenerate polygonal domain and hence can be constructed using a variation of the generalized Schwarz–Christoffel mapping for multiply connected domains recently obtained by Crowdy [6,7] in terms of the Schottky–Klein prime function. Our final expression for the conformal map revealing the bubble shapes will be given as an explicit indefinite integral whose integrand consists of products of Schottky–Klein prime functions and their derivatives.

There are several prior results pertaining to steady multiple bubbles in Hele-Shaw systems which we wish to survey to motivate the free boundary problem considered in this paper. An important assumption that is made in each of these works is the exclusion of surface tension effects. From a theoretical standpoint, this makes the problem analytically tractable and allows for exact solutions to be found. Taylor & Saffman [8] found an exact solution for a single bubble in a channel with reflectional symmetry about the channel centreline. Tanveer [9] was able to generalize this solution using elliptic function theory to describe a single asymmetric bubble in the channel. Vasconcelos [10] reported exact solutions for a finite number of steadily translating bubbles in a Hele-Shaw channel. He considered symmetrical bubble shapes and, by reducing the problem to a simply connected flow domain, was able to derive Schwarz–Christoffel-type formulae for the conformal mappings determining the bubble interfaces. Adopting a similar approach, Silva & Vasconcelos [11] have found exact solutions for a doubly periodic array of multiple symmetrical bubbles, with Schwarz–Christoffel methods again proving to be fruitful. Exploiting symmetry arguments, Vasconcelos [12,13] has also found families of exact solutions for various infinite streams of bubbles in the Hele-Shaw system. If symmetry is not enforced in configurations with multiple bubbles, the flow domain is inherently multiply connected, making the problem more mathematically challenging. Indeed, only a few exact solutions are known for Hele-Shaw flows in multiply connected regions, mostly for doubly connected cases. For instance, the Tanveer solution [9] for one asymmetric bubble involves conformal mappings between doubly connected domains. More recently, Silva & Vasconcelos [14] obtained a solution for a stream of asymmetric bubbles in a Hele-Shaw channel by deploying the Schwarz–Christoffel formula for doubly connected regions. A few exact solutions for doubly connected time-dependent Hele-Shaw flows have also been obtained [15–18].

Most relevant to our present free boundary problem is the work of Crowdy [19] who found analytical solutions determining the shapes of any finite number of steadily translating bubbles in an unbounded Hele-Shaw cell. He derived analytical expressions for both the complex potential and the conformal map from a bounded multiply connected circular domain to the exterior of the bubble assembly by using the Schottky–Klein prime function. The solutions presented herein can be viewed as the generalization to the channel geometry of the solutions found by Crowdy [19] for multiple bubbles in an unbounded Hele-Shaw cell. Our solutions thus account for the effect of the two channel walls, which greatly influence the nature of the free boundary problem. An earlier attempt at solving this problem was made by Crowdy [20] by introducing two different Schottky groups and their associated Schottky–Klein prime functions, but it was subsequently found that the construction of the second Schottky group was not entirely correct (DC Crowdy 2011, personal communication). Herein, the problem is solved by employing the generalized Schwarz–Christoffel formula. The results reported here give the complete set of solutions for a finite number of bubbles steadily moving in a Hele-Shaw channel, in the sense that for any given set of physical parameters, namely the velocity, areas and centroids of the bubbles, specifying a bubble configuration, the corresponding solution can in principle be obtained from our general formulae.

## 2. Problem formulation

### (a) The complex potentials

We consider the problem of *M* finite-area bubbles translating uniformly with speed *U* parallel to the *x*-axis in a Hele-Shaw channel filled with an incompressible viscous fluid. Without the loss of generality, we will assume that the Hele-Shaw channel has a width equal to 2 and that the viscous fluid (outside the bubbles) has a uniform speed *V* =1 in the far field; see figure 1 for a schematic. Our model will be centred around some simplifying assumptions in order to render the problem analytically tractable. We shall assume that the fluid inside the bubbles (say, air) has negligible viscosity so that the pressure inside each bubble is constant. We also neglect surface tension effects so that the viscous fluid pressure will have a constant value on each bubble boundary. We assume furthermore that the Hele-Shaw channel is horizontally placed so that the effects of gravity can be neglected. Finally, we shall neglect three-dimensional thin film effects.

Under the assumptions above, the motion of the viscous fluid in our Hele-Shaw channel is governed by Darcy's law
2.1where the velocity potential *ϕ*(*x*,*y*) is given by
2.2Here, **v** is the averaged fluid velocity across the channel, *p* is the viscous fluid pressure, *b* is the gap between the plates and *μ* is the fluid viscosity. From the incompressibility condition, ∇⋅**v**=0, it follows that *ϕ* satisfies Laplace's equation, ∇^{2}*ϕ*=0. It is therefore natural to formulate the free boundary problem to be solved in the complex *z*-plane, where *z*=*x*+i*y*.

Let us thus introduce the complex potential, *w*(*z*), defined by
2.3where *ψ* is the stream function associated with the velocity potential *ϕ*. As the viscous fluid is assumed to have unit speed in the far field, it follows that
2.4We shall label the viscous fluid region by *D*_{z} and the boundary of the *j*th bubble by ∂*D*_{j}, *j*=1,…,*M* (figure 1). Apart from a simple pole at infinity, the complex potential *w*(*z*) is analytic everywhere in the viscous fluid region *D*_{z} and must satisfy the following boundary conditions:
2.5and
2.6Condition (2.5) simply states that the channel walls, *y*=±1, are streamlines of the flow, whereas (2.6) follows from the fact that pressure *p* is constant on each bubble boundary. From conditions (2.5) and (2.6), one sees that the flow domain in the complex potential *w*-plane consists of a horizontal strip of width 2 with *M* vertical slits in its interior, where each slit corresponds to a bubble in the *z*-plane (figure 2*a*).

Let us also introduce the function *τ*(*z*) corresponding to the complex potential in a frame of reference co-travelling with the bubbles
2.7From (2.4) and (2.7), one immediately sees that the far-field behaviour of *τ*(*z*) is
2.8which implies the following boundary conditions for *τ*(*z*) on the channel walls
2.9Furthermore, in the co-travelling frame the bubble boundaries are necessarily streamlines of the flow, thus
2.10It then follows from (2.9) and (2.10) that the flow domain in the *τ*-plane is a horizontal strip of width 2(*U*−1) with horizontal slits in its interior corresponding to the bubbles (figure 2*b*).

### (b) The conformal mapping

In order to determine the shapes of the bubbles, we shall construct the conformal map *z*(*ζ*) from a bounded *M*+1 connected circular domain *D*_{ζ} to the *M*+1 connected fluid region *D*_{z} exterior to *M* bubbles in the *z*-plane. We choose *D*_{ζ} to be the unit circle *ζ*-disc with *M* smaller non-overlapping discs excised from it. Label the unit circle by *C*_{0} and label the *M* inner circular boundaries as *C*_{1},…,*C*_{M}, and let the centre and radius of *C*_{j} be *δ*_{j} and *q*_{j}, respectively. A schematic of *D*_{ζ} is shown in figure 3 in the case where *M*=2 (triply connected). Let the unit circle *C*_{0} map to the channel walls and the interior circles *C*_{1},…,*C*_{M}, to the bubble boundaries ∂*D*_{1},…,∂*D*_{M}, respectively. This implies that the mapping function *z*(*ζ*) will necessarily have two logarithmic singularities on *C*_{0}. By the degrees of freedom afforded by the Riemann–Koebe mapping theorem [21], we can place these logarithmic singularities at *ζ*=±1.

Let us define the functions *W*(*ζ*) and *T*(*ζ*) through the following compositions:
and
These functions must be analytic in the circular domain *D*_{ζ} and satisfy appropriate boundary conditions on the circles *C*_{j}, as discussed next. From conditions (2.5) and (2.6), it follows that *W*(*ζ*) must be such that
2.11and
2.12Similarly, for *T*(*ζ*) one has from (2.9) and (2.10) that
2.13and
2.14In (2.11) and (2.13), the upper and lower signs correspond to the upper and lower semicircles on the unit circle, respectively.

Now, from (2.7) it follows that the mapping function *z*(*ζ*) can be written as
2.15We have thus reduced our original free boundary problem to the more manageable task of computing two analytic functions, *W*(*ζ*) and *T*(*ζ*), that map the circular domain *D*_{ζ} to multiply connected slit domains which can be viewed as degenerate polygonal domains. These functions can then be readily obtained using a generalized Schwarz–Christoffel mapping for multiply connected polygonal domains, which is briefly reviewed in §3.

## 3. Generalized Schwarz–Christoffel mapping

Consider the function *z*(*ζ*) defined on the circular domain *D*_{ζ} by the following expression:
3.1where is a complex constant and *ω*(*ζ*,*γ*) is the Schottky–Klein prime function associated with the domain *D*_{ζ}. For a definition of the Schottky–Klein prime function and a discussion of some of its properties, see, e.g. Crowdy [6]. The Schottky–Klein prime function has deep connections with Riemann surface theory [22], but for the purposes of this paper it suffices to think of it as a special computable function [23].

It is shown elsewhere [24] that the function defined in (3.1) conformally maps the circular domain *D*_{ζ} onto a bounded (*M*+1)-connected polygonal domain in the *z*-plane, where the unit circle *C*_{0} is mapped to the outer boundary polygon *P*_{0} and the interior circles *C*_{1},…,*C*_{M} are mapped to the inner polygonal boundaries *P*_{j}, *j*=1,…,*M*, respectively. The points , *k*=0,1,…,*n*_{j}, are the preimages of the vertices on polygon *P*_{j}, with being the corresponding turning angles so that the internal angle at the respective vertex is , with . The outer polygon can have one or more vertices at infinity, i.e. , for some *k*, in which case because the angle is measured with respect to the variable 1/*z* for . Inner polygons are also allowed to degenerate to line segments (i.e. slits).

The set of points appearing in formula (3.1) are obtained by computing the zeros (on each inner circle) of the following equation:
3.2This guarantees that the derivative *z*_{ζ}(*ζ*) has no poles at the points and , so that the only singularities of *z*(*ζ*) in *D*_{ζ} are the appropriate branch points at . Formula (3.1) is an alternative version of the generalized Schwarz–Christoffel mapping derived by Crowdy [6] that is more convenient to treat multiply connected strip domains.

For later use, we quote here an important property of the Schottky–Klein prime function [6]: the functions defined by
3.3have constant argument on all circles *C*_{k}. More specifically, one has [24]
3.4where
3.5Here, the functions *v*_{j}(*ζ*), *j*=1,…,*M*, are the *M* integrals of the first kind associated with the domain *D*_{ζ}. An algorithm for computing *v*_{j}(*ζ*) can be found in [23]. For convenience, we have defined *v*_{0}(*ζ*)≡0 and introduced the notation , with similar definition for . In §4, we shall use generalized Schwarz–Christoffel formula (3.1) and property (3.4) to give an explicit construction of the complex potentials *W*(*ζ*) and *T*(*ζ*) defined above.

## 4. The general solution

### (a) The function *T*(*ζ*)

Recall that the mapping *τ*=*T*(*ζ*) conformally maps the circular domain *D*_{ζ} onto a strip domain in the *τ*-plane where the unit circle *C*_{0} is mapped to the strip boundaries, Im[*τ*]=±(*U*−1), and the inner circles *C*_{j} are mapped to horizontal slits (figure 2*b*). This domain can alternatively be viewed as a multiply connected degenerate polygonal domain, where the outer polygonal boundary has only two edges that meet at and the inner polygons degenerate to slits. The corresponding turning angle parameters are therefore given by at the two vertices at infinity and , *j*=1,…,*M*, at the endpoints of the slits. Following previous notation, we label by the preimages in the *ζ*-plane of the endpoints of the *τ*-strip and by , *j*=1,…,*M*, the preimages of the slit endpoints. Now recall that the preimages in the *ζ*-plane of the channel left and right ends (which correspond to ) have been chosen to be *ζ*=∓1, and so we have and . Applying the Schwarz–Christoffel formula (3.1) to this case then yields
4.1

As discussed in §3, the points , are obtained by computing the solutions of (3.2). In this case, because the slits are all horizontal, one can select the points and to be the preimages of the slit endpoints [24], i.e.
4.2so that the terms in the numerator and in the denominator of the product appearing in (4.1) will precisely cancel out. Hence, (4.1) simplifies to
4.3which upon integration becomes
4.4where we have set to account for the fact that width of the strip in the *τ*-plane is 2(*U*−1) (figure 2*b*). It is worth noting that mapping (4.4) can be obtained directly from the properties of the Schottky–Klein prime function (e.g. Crowdy [20]). It is instructive however to demonstrate that it can be recovered from Schwarz–Christoffel mapping (3.1), as shown above.

### (b) The function *W*(*ζ*)

As discussed in §2*a*, the complex potential *W*(*ζ*) maps the circular domain *D*_{ζ} onto a multiply connected domain in the *w*-plane consisting of a horizontal strip of width 2 with *M* vertical slits in its interior, as shown in figure 2*a*. This slit strip domain is therefore of the same type as the flow domain in the *τ*-plane discussed above, the only difference being that the slits are now vertical. Let us then denote by the set of preimages in the *ζ*-plane of the endpoints of the slits in the *w*-plane. It should be clear from the preceding discussion that the function *W*_{ζ} is given by
4.5where is a complex constant. The modulus of has to be chosen such that the width of the strip in the *w*-plane equals 2; see below. On the other hand, the argument of and the parameters are determined from the requirement that *C*_{0} is mapped by *W*(*ζ*) to horizontal straight lines and that the inner circles map to vertical slits, as shown next.

First, it is convenient to rewrite (4.5) in terms of the functions *F*_{j}(*ζ*;*ζ*_{1},*ζ*_{2}) defined in (3.3)
4.6The requirement that *C*_{0} is mapped by *W*(*ζ*) to horizontal walls means that
4.7which implies that , for some integer *n*, as one traverses an infinitesimal angle d*θ* on the unit circle. This in turn implies from (4.6) and (3.4) that
4.8where *Q*_{jk}(*ζ*_{1},*ζ*_{2}) is given in (3.5). Similarly, the condition that the inner circles *C*_{k} are mapped to vertical slits can be written as
4.9which implies that , for some integer *m*, as one traverses an infinitesimal angle d*θ* on the inner circle *C*_{k}. From this condition and (4.6), it follows that
4.10

Equations (4.8) and (4.10) give *M*+1 conditions on the 2*M*+1 parameters corresponding to the argument of and the points . The other set of *M* conditions necessary to determine these parameters comes from the requirement that *W*(*ζ*) be single-valued in *D*_{ζ}. This implies that a 2*π*-traversal around *C*_{j} should correspond to returning to the same point on the *j*th vertical slit, i.e.
4.11Thus, once the conformal moduli *q*_{j} and *δ*_{j} are specified, the conditions (4.8), (4.10) and (4.11) give 2*M*+1 real equations which can be solved numerically to determine and the points . From a numerical viewpoint, however, working with equations (4.8) and (4.10) is more cumbersome because one does not know *a priori* the values of the integers *n* and *m* (which have to be found by trial and error). To circumvent this problem, we enforce boundary conditions (4.7) and (4.9) on a particular point, say, *θ*=*π*/2, and then solve these equations, together with (4.11), via a multivariate Newton's method for root finding. (In the examples shown in the next section, the corresponding accessory parameters were obtained using this method.)

The last parameter that needs to be determined is the modulus of the premultiplier . This is obtained from the requirement that the logarithmic singularities of *W*(*ζ*) at *ζ*=±1 have the appropriate strength, so that the jump in *W*(*ζ*) when crossing either of these singularities equals 2i, which corresponds to the channel width in the *w*-plane (figure 2*a*). The modulus is then found to be
4.12where we have used that *ω*(*ζ*,*γ*)=*ζ*−*γ*, as *ζ*→*γ*. This completes the construction of the function *W*_{ζ}(*ζ*).

### (c) Conformal map *z*(*ζ*)

In light of definition (2.15), the conformal map *z*(*ζ*) that we seek will have the following integral form:
4.13where is a constant, *ζ*_{0}∈C is an arbitrary point inside *D*_{ζ} and expressions for *T*_{ζ}(*ζ*) and *W*_{ζ}(*ζ*) are given in (4.3) and (4.5), respectively. Without the loss of generality, we can set the bubble velocity to *U*=2; it is demonstrated by Vasconcelos [10] that all other bubble assemblies corresponding to different values of *U* can be obtained from the *U*=2 solutions by a simple rescaling.

Let us recall that we have used up two of the three degrees of freedom associated with the Riemann–Koebe mapping theorem, namely, choosing *ζ*=±1 to map to the channel ends. The remaining degree of freedom can now be used to fix the value of the constant . We are then left with 3*M* free real parameters corresponding to the 3*M* conformal moduli of our circular domain *D*_{ζ}. Physically, these parameters correspond to the area and centroids for each of the *M* bubbles. Once the 3*M* conformal moduli are prescribed, we can determine all other parameters entering (4.13), namely, the set of points on the inner circles and the premultiplier , and thus a specific solution corresponding to a particular bubble assembly is obtained. In §5, we illustrate the foregoing theory by considering some specific examples of various bubble configurations.

## 5. Examples

Solutions for multiple steady bubbles that either are reflectionally symmetric about the channel centreline or have fore-and-aft symmetry have been obtained by Vasconcelos [10] by reducing the problem—on account of the symmetry—to a simply connected domain, and then applying the standard Schwarz–Christoffel formula. These symmetric solutions are readily recovered in our formalism by choosing circular domains *D*_{ζ} with the appropriate symmetry: for bubbles with centreline (fore-and-aft) symmetry, one must choose a domain *D*_{ζ} that is symmetric with respect to reflections about real (imaginary) axis. Figure 4 shows an example of two bubbles whose centroids are aligned along the channel centreline and which are reflectionally symmetric about it. In fact, we were able to successfully recover similar bubble shapes as in figure 4 using the analytical solutions of Vasconcelos [10]. The formalism presented here is however much more general is that it applies to arbitrary bubble configurations, i.e. with no imposed symmetry *a priori*. Figure 5 shows an example of two asymmetric bubbles. Several streamlines of flow field in the co-travelling frame have also been plotted—these streamlines provide a qualitative check on the solutions. Streamlines are shown only for this case because they demand considerable numerical work.

The versatility and generality of our method can be demonstrated through solving for parameters yielding a larger number of bubbles in some asymmetric configuration. Figure 6 shows an example of the bubble shapes for a particular asymmetric assembly of three bubbles, whereas figure 7 reveals the shape of the bubble boundaries in a particular four-bubble configuration. These assemblies are not symmetric about any axis and the bubbles in each assembly all have different areas (corresponding to the different choices of *q*_{j} values). Solutions for a higher number of bubbles can be treated in a similar manner but the numerical computation of the parameters becomes increasingly more expensive.

In each of the figures 4–7, the interaction between the bubbles themselves, and between the bubbles and the channel walls, is clearly visible from the shapes of their boundaries. As specified in the respective figure captions, we have made particular choices of the conformal moduli defining *D*_{ζ} in order to find the finite set of accessory parameters in our conformal map determining the bubble shapes. Alternatively, we could have specified the areas and the centroids of each of the individual bubbles and solved for the conformal moduli of *D*_{ζ}. However, this would have been a rather challenging numerical undertaking, owing to the fact that the parameters (the set of preimages of the endpoints of the *M* horizontal slits in the *τ*-plane) would need to be calculated on each iterative step.

## 6. Discussion

We have presented analytical solutions to the free boundary problem of determining the interface shapes of a finite number *M* of bubbles steadily translating along a Hele-Shaw channel. To do this, we found a concise formula in the form of an explicit indefinite integral for the conformal map from a bounded (*M*+1)-connected circular domain to the fluid region in the channel exterior to the *M* bubbles. The integrand of this indefinite integral is neatly expressed in terms of products of Schottky–Klein prime functions and their derivatives, and is known up to a finite set of accessory parameters to be found as part of the solution.

Our formula for the conformal map determining the bubble shapes is very general and makes no *a priori* symmetry assumptions concerning the geometrical arrangement of the bubbles. Indeed, all the geometrical information about the physical domain is encapsulated in the prescription of the preimage domain *D*_{ζ} over which each of the Schottky–Klein prime functions appearing in (4.13) is defined. In previous works on steady bubbles in a Hele-Shaw channel, symmetry had to be enforced in order to make progress. Solutions for a single bubble with centreline symmetry were first found by Taylor & Saffman [8]. These solutions were later extended by Tanveer [9] to include asymmetric bubbles; an alternative derivation of the Tanveer solutions was subsequently given by Combescot & Dombre [25] using Riemann–Hilbert methods. A general class of exact solutions for multiple steady bubbles with either centreline or fore-and-aft symmetry were obtained by Vasconcelos [10] by reducing the problem—on account of symmetry—to a simply connected flow domain. Our solution scheme readily incorporates the solutions found in these works, and so these solutions can be viewed as special cases of ours. It also generalizes the solutions for multiple bubbles in an unbounded Hele-Shaw cell obtained by Crowdy [19] by including the effect of the channel walls. Indeed, it is possible to retrieve the solutions of Crowdy [19] from those presented here by taking the limit as the channel width becomes infinite. This can be accomplished in our formalism by choosing a circular domain *D*_{ζ} where an inner circle is mapped to the channel walls (with the unit circle mapped to one of the bubbles) and then taking the radius of this inner circle to zero, but we have not pursued this detail here.

In conclusion, we have devised a constructive method for finding solutions to the free boundary problem of multiple bubbles (with no assumed symmetry) in a Hele-Shaw channel. The crucial step in our scheme was to represent the conformal mapping from a circular domain to the physical flow region as a sum of two analytic functions that map the circular domain to multiply connected degenerate polygonal domains. These analytic functions, which correspond to the complex potentials in the laboratory and moving frames, were then obtained from the generalized Schwarz–Christoffel formula for multiply connected domains in terms of Schottky–Klein prime functions. We considered examples of various bubble assemblies, demonstrating that our solution scheme is capable of modelling any finite number of bubbles in a particular assembly.

An interesting extension of this work would be to consider a doubly periodic array of asymmetric bubble assemblies. This, in turn, would generalize the solutions of Silva & Vasconcelos [11] for a doubly periodic array of symmetric bubbles. It is worth noting that an exact solution for an asymmetric stream of bubbles (with one bubble per unit cell) in a Hele-Shaw channel was recently found by Silva & Vasconcelos [14] using conformal mapping between doubly connected domains. It is possible to extend these solutions using the ideas presented in this paper to include a periodic assembly of bubbles with multiple asymmetric bubbles per unit cell; work is currently in progress in this direction. Another interesting line of enquiry would be to investigate the effects of surface tension on the bubble boundaries. This gives rise to the so-called selection problems: for non-zero surface tension, there is no longer a continuum of bubble velocities for which solutions can exist, and so only a discrete set of velocities are allowed. It is also important to point out that obtaining exact analytical solutions for steady Hele-Shaw systems naturally paves the way to finding time-dependent solutions which are also of great physical and mathematical interest. An example of the extension of the steady theory to the time-dependent case is the recent study by Vasconcelos & Mineev-Weinstein [18], who have been able to determine an exact solution without surface tension for the time evolution of a bubble of arbitrary initial shape using conformal mappings between doubly connected domains. These authors demonstrated that the solution with *U*=2 is the only attractor for the non-singular solutions, thus showing that selection is inherently determined by the zero surface tension dynamics. They also conjectured [18] that a similar scenario holds in domains of arbitrary connectivity. Our new solutions—in terms of Schottky–Klein prime functions for multiply connected circular domains—should therefore serve as a starting point for future endeavours looking into constructing time-dependent solutions for multiple bubbles in a Hele-Shaw channel, from which the general selection problem can hopefully be addressed.

## Funding statement

C.C.G. acknowledges financial support from a Doctoral Prize Fellowship of the Engineering and Physical Sciences Research Council (UK). G.L.V. acknowledges financial support from a scholarship from the Conselho Nacional de Desenvolvimento Cientifico e Tecnologico (Brazil) for a sabbatical stay at Imperial College London.

## Acknowledgements

The authors thank D. G. Crowdy for helpful discussions. C.C.G. is appreciative of the hospitality of the Department of Physics at the Federal University of Pernambuco where part of this work was carried out. G.L.V thanks the hospitality of the Department of Mathematics at Imperial College London where this work was completed.

- Received October 17, 2013.
- Accepted December 2, 2013.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.