Hostname: page-component-848d4c4894-xfwgj Total loading time: 0 Render date: 2024-06-16T22:34:28.101Z Has data issue: false hasContentIssue false

On the selection of Saffman–Taylor viscous fingers for divergent flow in a wedge

Published online by Cambridge University Press:  24 May 2024

Cecilie Andersen*
Affiliation:
Department of Mathematical Sciences, University of Bath, Bath BA2 7AY, UK
Christopher J. Lustri
Affiliation:
School of Mathematics and Statistics, The University of Sydney, Sydney NSW 2006, Australia
Scott W. McCue
Affiliation:
School of Mathematical Sciences, Queensland University of Technology, Brisbane QLD 4001, Australia
Philippe H. Trinh*
Affiliation:
Department of Mathematical Sciences, University of Bath, Bath BA2 7AY, UK Theoretical Sciences Visiting Program (TSVP), Okinawa Institute of Science and Technology Graduate University, Onna 904-0495, Japan
*
Email addresses for correspondence: ca751@bath.ac.uk, p.trinh@bath.ac.uk
Email addresses for correspondence: ca751@bath.ac.uk, p.trinh@bath.ac.uk

Abstract

We study self-similar viscous fingering for the case of divergent flow within a wedge-shaped Hele-Shaw cell. Previous authors have conjectured the existence of a countably infinite number of selected solutions, each distinguished by a different value of the relative finger angle. Interestingly, the associated solution branches have been posited to merge and disappear in pairs as the surface tension decreases. For the first time, we demonstrate how the selection mechanism can be derived based on exponential asymptotics. Asymptotic predictions of the finger-to-wedge angle are additionally given for different sized wedges and surface-tension values. The merging of solution branches is explained; this feature is qualitatively different to the case of classic Saffman–Taylor viscous fingering in a parallel channel configuration. Moreover, because the asymptotic framework does not highly depend on specifics of the wedge geometry, the proposed theory for branch merging in our self-similar problem likely relates much more widely to tip-splitting instabilities in time-dependent flows in circular and other geometries, where the viscous fingers destabilise and divide in two.

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2024. Published by Cambridge University Press.

1. Introduction

In their now classic work on viscous fingering, Saffman & Taylor (Reference Saffman and Taylor1958) consider the situation of a single finger of air – or an otherwise less viscous substance – steadily penetrating a Hele-Shaw cell filled with a viscous fluid. A key quantity of interest is the proportion of the channel occupied by the width of the finger, denoted $\lambda \in (0,1)$. Experimentally, Saffman and Taylor observed that $\lambda \approx 1/2$. However, their asymptotic analysis, valid at small values of the non-dimensional surface-tension parameter $\epsilon \to 0$, did not seem to restrict the value of $\lambda$. Today, it is known that for a fixed $\epsilon$ there exists a countably infinite number of possible values of $\lambda = \lambda _i$ with the property that

(1.1)\begin{equation} 1/2 <\lambda_1(\epsilon)<\lambda_2(\epsilon)< \cdots <1. \end{equation}

In the limit $\epsilon \to 0$, every element in the family converges to 1/2 (see e.g. Vanden-Broeck Reference Vanden-Broeck2010). The resultant plot of the eigenvalue, $\lambda$, vs the surface-tension parameter, $\epsilon$, is shown in figure 1(a). The resolution of the Saffman–Taylor problem, including the asymptotic derivation of (1.1), is obtained using exponential asymptotics; today, it is a prototypical example of such beyond-all-orders asymptotics (see e.g. works by McLean & Saffman Reference McLean and Saffman1981; Combescot et al. Reference Combescot, Dombre, Hakim, Pomeau and Pumir1986, Reference Combescot, Hakim, Dombre, Pomeau and Pumir1988; Hong & Langer Reference Hong and Langer1986; Tanveer Reference Tanveer1987, Reference Tanveer2000; Chapman Reference Chapman1999; Pullin & Meiron Reference Pullin and Meiron2013).

Figure 1. (a) The bifurcation diagram for the classic Saffman–Taylor problem in a channel showing the first 10 branches of selected $\lambda$ values against the surface-tension parameter $\epsilon ^2$. Labels are shown for the first four branches: $\lambda _1$, $\lambda _2$, $\lambda _3$, $\lambda _4.$ This is equivalent to the results of Chapman (Reference Chapman1999). (b) The solid lines show the bifurcation diagram which we calculate using exponential asymptotics for wedge angle $\theta _0 = 20^\circ$. This shows the permitted $\lambda (\epsilon )$ values that are selected by the selection mechanism. Details of the derivation of this selection mechanism follow in the rest of the paper. The circles show the numerically calculated values extracted from Ben Amar (Reference Ben Amar1991b). Notice that $\lambda _1$ and $\lambda _2$ merge at $\epsilon \approx 0.3$ and $\lambda _3$ and $\lambda _4$ merge at $\epsilon \approx 0.07$.

In this paper, we are interested in deriving a similar selection law to (1.1) but for an analogue problem where fluid is injected at the corner of a Hele-Shaw cell limited by side walls consisting of a wedge of specified angle. Although this wedge scenario is interesting in its own right, it gains further importance as a partial model for the fingering seen where fluid is injected into a Hele-Shaw cell from a central source. As the injected fluid moves outwards, the interface destabilises and fingering occurs in the manner illustrated in figure 2. Thus, the limited wedge configuration considered in this work serves as a model for each sector of the full source problem.

Figure 2. Sketch of the free surface for viscous fingering in the full circular geometry, where a central source injects fluids outwards in all directions. Assuming axisymmetry, a single finger from this free surface can be considered as arising due to injecting fluid in the corner of the Hele-Shaw cell limited by sidewalls consisting of a wedge of angle $\theta _0$.

Consider the wedge problem characterised by a wedge angle, $\theta _0$, measured at the corner. The key parameter, $\lambda$, now describes the angular proportion, $\lambda \theta _0$, of the cell occupied by the finger. In comparison with the previous figure 1(a) for classic Saffman–Taylor viscous fingering, in figure 1(b) we plot the bifurcation diagram for the case of a wedge of angle $\theta _0 = 20^\circ$. The solid lines correspond to the new asymptotic predictions developed in this work, while the circles correspond to the previous numerical results of Ben Amar (Reference Ben Amar1991b). The figure shows the existence of distinct solution families, $\lambda = \lambda _i(\epsilon, \theta _0)$.

In addition to the expected solution families, $\lambda _i$, there is now an additional phenomenon for the wedge-limited configurations. As seen in figure 1(b), there exist certain critical values of the surface-tension parameter, $\epsilon$, where each solution family, $\lambda _i$, reaches a turning point connected to the adjacent family, i.e. $\lambda _i = \lambda _{i+1}$. In this work, we will explain how this merging of eigenvalues can be understood from the perspective of the exponential asymptotics. As noted by Ben Amar (Reference Ben Amar1991b), the merging of eigenvalue pairs causes a loss of existence of the solutions for sufficiently small $\epsilon$. Physically, in connection with the full geometry in figure 2, this results in the finger splitting into two through a tip-splitting instability. One should then consider two wedges of half the angle to continue to follow the finger profiles. It is then interesting to consider the consequences of the exponential asymptotic analysis towards the more complicated problem of time-dependent tip-splitting instabilities in the unrestricted planar domain. This will be further discussed in § 10.

1.1. Background and open challenges of the wedge problem

The literature on Hele-Shaw flows and viscous fingering problems is extensive; here, we provide a review of selected papers, primarily focused on wedge configurations or closed bubbles in channels.

Experimental observations and initial analysis for the wedge geometry can be found in the works of Paterson (Reference Paterson1981) and Thomé et al. (Reference Thomé, Rabaud, Hakim and Couder1989). Then, in the early 1990s, a number of works appeared on the beyond-all-orders aspects of the wedge configuration (Brener et al. Reference Brener, Kessler, Levine and Rappel1990; Ben Amar Reference Ben Amar1991a,Reference Ben Amarb; Tu Reference Tu1991; Combescot Reference Combescot1992; Levine & Tu Reference Levine and Tu1992). Combescot (Reference Combescot1992) identified the singularities in the complex plane that are responsible for the selection mechanism. The solvability condition was found exactly by Brener et al. (Reference Brener, Kessler, Levine and Rappel1990) for the case with a $90^{\circ }$ separation angle, where the solution reduces to a closed analytic form. Both Combescot (Reference Combescot1992) and Tu (Reference Tu1991) used WKBJ (Wentzel-Kramers-Brillouin-Jeffreys or Liouville–Green) methods to make analytical progress on the problem with a general wedge angle whilst numerical results were obtained by Ben Amar (Reference Ben Amar1991b). More recently, there has been additional experimental works, showing clear photographs, of the tip-splitting instabilities in the wedge problem by Lajeunesse & Couder (Reference Lajeunesse and Couder2000).

This paper is most strongly motivated by the work of Tu (Reference Tu1991) and Ben Amar (Reference Ben Amar1991b). Tu (Reference Tu1991) had linearised the free-surface problem with a general wedge angle to obtain a model differential equation. For this new equation, a WKBJ approximation was used to derive a solvability condition that can predict the theoretical zero-surface-tension limit for $\lambda$, as well as a condition that finds the point in the bifurcation diagram where branch merging occurs.

Ben Amar (Reference Ben Amar1991b) produced accurate numerical results for parts of the bifurcation diagram. However, on account of the challenges in numerical computations of the eigenvalue problem, Ben Amar (Reference Ben Amar1991b) noted that:

Our predictions concerning levels [branches] higher than the first two require confirmation by a very careful WKB analysis, which is the most suitable treatment at extremely low surface tension. Probably, the results of analytic solutions without surface tension will make this analysis possible.

At this point (and until the present) we do not believe any group has managed to derive the exact selection mechanism (i.e. the missing analysis referenced above).

Modern techniques of exponential asymptotics allow us to study the wedge problem in the small-surface-tension limit without the need to linearise in the same fashion as previous practitioners (Chapman, King & Adams Reference Chapman, King and Adams1998). In this paper, we will use these techniques to address the open problem identified by Ben Amar (Reference Ben Amar1991b) and obtain an analytic solvability condition for the selected eigenvalues.

Because the approach presented here combines a hybrid asymptotic–numerical insight, it is likely that the methods can be extended to much more general flow configurations. More recently, further work has been done on closely related problems in Hele-Shaw channels. There is particular interest in Hele-Shaw channels with a central raised rail, which can change the stability of the finger (Thompson, Juel & Hazel Reference Thompson, Juel and Hazel2014; Franco-Gómez, Thompson & Hazel Reference Franco-Gómez, Thompson and Hazel2016). Further, there have been many recent experimental (Franco-Gómez et al. Reference Franco-Gómez, Thompson, Hazel and Juel2017; Gaillard et al. Reference Gaillard, Keeler, Le Lay, Lemoult, Thompson, Hazel and Juel2021; Lawless et al. Reference Lawless, Keeler, Hazel and Juel2023), numerical (Franco-Gómez et al. Reference Franco-Gómez, Thompson, Hazel and Juel2017; Keeler et al. Reference Keeler, Thompson, Lemoult, Juel and Hazel2019; Thompson Reference Thompson2021) and analytical (Keeler et al. Reference Keeler, Thompson, Lemoult, Juel and Hazel2019; Booth, Griffiths & Howell Reference Booth, Griffiths and Howell2023) developments concerning closed bubbles propagating along Hele-Shaw channels. We will return at the end of the paper to discuss the connection of our work with these newer problems.

2. Mathematical formulation

A traditional Hele-Shaw cell consists of two parallel plates separated by a small distance, $b$, and filled with viscous fluid with viscosity, $\mu$. For the case of a circular geometry, like that depicted in figure 2, an inviscid fluid is injected at a point, and this drives an outwards-expanding interface between viscous and inviscid fluids.

We now consider the related wedge-shaped geometry, as shown in figure 3(a). Here, the inviscid fluid is injected at the wedge corner at some prescribed flow rate. The viscous fluid is constrained to lie between the wedge walls separated by an internal angle $\theta _0$. As can be observed experimentally (Thomé et al. Reference Thomé, Rabaud, Hakim and Couder1989), a self-similar shape is reached eventually, where the inviscid fluid occupies an angle $\lambda \theta _0$, with $0<\lambda <1$. This set-up is referred to as divergent flow in Ben Amar (Reference Ben Amar1991a). For the case of zero-surface-tension, a prototypical solution is shown in figure 3(a), and we observe the petal-shaped interface between viscous fluid and inviscid fluid.

Figure 3. (a) A numerical plot of the top-down view of the self-similar physical profile in the $\hat {z}$-plane is shown for the zero-surface-tension case, with parameter values $\theta _0 = 20^{\circ }$ and $\lambda = 0.6$. The free surface was computed using (3.7). The Hele-Shaw cell is bounded by the thick black lines and is filled with a viscous fluid, shown in grey. An inviscid fluid is injected from the corner of the wedge and forms a finger with angle $\lambda \theta _0$. The corner of the wedge lies at $\hat z = 0$ ($BF$) and the tip of the finger lies at $\hat {z}=1$ ($CE$). (b) A sketch of the $z$-plane, computed via the conformal map $z=(2/\theta _0)\log \hat {z}$. This configuration is analogous to the traditional Saffman–Taylor finger (McLean & Saffman Reference McLean and Saffman1981).

The classic problem of Saffman–Taylor viscous fingering in a channel is typically studied in a travelling frame corresponding to a steady-state finger. In the wedge geometry, the analogue of a travelling wave frame of reference is a self-similar solution. In Appendix A, we demonstrate how the original equations of potential flow for a Hele-Shaw cell, with boundary conditions on the channel walls and free boundary, and an injection condition can be reposed in the self-similar framework. The key idea relates to a transformation of the original dimensional lengths, $\bar {x}$ and $\bar {y}$, which are scaled via

(2.1)\begin{equation} (\hat{x}, \hat{y}) = \frac{(\bar{x}, \bar{y})}{R_0 {A}(t)}, \end{equation}

where $R_0$ is a length scale (chosen to be the distance between the corner of the wedge and the tip of the finger at $t=0$) and ${A}(t)$ is a dimensionless scaling factor that depends on dimensionless time, $t$. As shown in Appendix A, two possible choices of ${A}(t)$ result in effectively the same self-similar problem, given by (A6), in dimensionless variables with a time-independent effective surface-tension parameter $\sigma$. This $\sigma$, when rescaled, then gives us our small parameter $\epsilon$, defined below in (2.10); see Ben Amar (Reference Ben Amar1991b) for further details.

We thus have the governing set of potential-flow equations in (A6) corresponding to a scaled velocity potential, $\hat {\phi }$, and self-similar physical-plane coordinates, written in complex form as $\hat {z} = \hat {x} + \mathrm {i} \hat {y}$ (figure 3a). Note that the free surface is now stationary. We introduce the complex potential, $\hat {f}=\hat {\phi }+\mathrm {i}\hat {\psi }$, with harmonic conjugate given by the streamfunction $\hat{\psi}$. Within the $\hat {f}$-plane, fluid is located in the infinite strip bounded by $\hat {\psi } \in [-1,1]$.

Finally, the flow domain is mapped to a channel geometry via

(2.2)\begin{equation} z= \frac{2}{\theta_0}\log \hat{z}, \end{equation}

so that the walls $BA$ and $FG$ lie on $\operatorname {Im} z=\pm 1$, respectively, and the tip $CE$ is fixed to the origin $z=0$. This is shown in figure 3(b).

Following § 3B of Ben Amar (Reference Ben Amar1991b), we review the procedure for developing a set of boundary-integral equations for the potential-flow problem. First, the velocity potential is shifted as

(2.3)\begin{equation} f^*=\phi^*+\mathrm{i}\psi^*=\frac{\hat{f}-H(z)}{(1-\lambda)Q_0/\lambda}. \end{equation}

Above, $2Q_0$ is the dimensionless flux of fluid across the wedge at infinity in the self-similar frame (cf. later (2.8)). The function $H(z)$ is implicitly defined so that the free surface lies on $\psi ^*=$ constant (cf. (3.10) of Ben Amar Reference Ben Amar1991b). For the case of the classic Saffman–Taylor viscous fingering problem in a parallel channel, $H(z) = z$, as shown by McLean & Saffman (Reference McLean and Saffman1981).

As in Chapman (Reference Chapman1999), it is convenient for later analysis to map the fluid region to the upper-half-$\zeta$-plane via

(2.4)\begin{equation} \zeta^2 = \mathrm{e}^{{\rm \pi} f^*} - 1. \end{equation}

The mapped fluid domain for the configuration in figure 3 is shown in figure 4. Thus we see that, under the map (2.4), the free surface ($BCEF$) lies on the real $\zeta$-axis while the tip of the finger ($CE$) is at $\zeta = 0$.

Figure 4. The fluid region (grey) is mapped to the upper-half-$\zeta$-plane. The key points from figure 3 are labelled here. Within the $\zeta$-plane, there is a branch cut from the point $D$ ($\zeta = \mathrm {i}$). Here, the branch cut is taken vertically up the imaginary axis from $\zeta = \mathrm {i}$.

In the governing equations to follow, we shall seek to solve for the unknown free-surface location and fluid velocities along the interface, $\zeta = \xi \in \mathbb {R}$. It is convenient to introduce quantities $q$ and $\tau$ via

(2.5)\begin{equation} \frac{q}{1-\lambda}\mathrm{e}^{-\mathrm{i}\tau} = \frac{\operatorname{d}{}f^*}{\operatorname{d}{}z} = \frac{\partial\phi^*}{\partial x} + \mathrm{i} \frac{\partial\phi^*}{\partial y}, \end{equation}

which are analogues of speed and streamline angle, respectively (and reduce to the actual fluid speed and streamline angle in the limit $\theta _0\rightarrow 0$). Therefore, we require a set of governing equations for the unknowns $(x(\xi ), y(\xi ), q(\xi ), \tau (\xi ))$.

With the various conformal maps now established, at this point, we may follow the same procedures as found in § 3B of Ben Amar (Reference Ben Amar1991b). We find on the free surface, where $\zeta = \xi$ is real, that continuity of pressure yields Bernoulli's equation

(2.6) \begin{align} &\epsilon^2 \frac{P(\xi)}{4} \frac{\partial}{\partial \xi}\left(r(x) \left[{-}P(\xi) q(\xi)\frac{\partial \tau}{\partial \xi} + \frac{\ell}{2}\sin\tau(\xi)\right]\right) \nonumber\\ &\quad\quad = \frac{Q_0}{1-\lambda} - \frac{2}{\rm \pi}(1+\xi^2){\int\hskip -1,05em -\,}_{-\infty}^0 \frac{K(\tilde{\xi})\operatorname{d}{}{\tilde{\xi}}}{P(\tilde{\xi})(\xi^2- \tilde{\xi}^2)}, \end{align}

which can be compared with (3.10) of Ben Amar (Reference Ben Amar1991b). In (2.6), we have defined a number of functions for convenience. Firstly, we have written

(2.7ac)\begin{equation} P(\xi) = (1+\xi^2)/\xi, \quad r(x) = \exp({-\theta_0 x(\xi)/2}), \quad K(\xi) = \frac{\sin [\tau(\xi)] }{q(\xi) [r(x(\xi))]^2}. \end{equation}

Note that $r$ is a function of $x$ which, in turn, depends on $\xi$. In (2.6), we have defined $Q_0$ to be

(2.8)\begin{equation} Q_0 = \frac{2(1-\lambda)}{\rm \pi} \int_{-\infty}^0 \frac{K(\tilde{\xi})}{P(\tilde{\xi})}\,{\rm d}{\tilde{\xi}}, \end{equation}

which, as mentioned above, represents a dimensionless constant describing the fluid flux. It is also convenient to define a scaled value for the interior wedge angle $\theta _0$

(2.9)\begin{equation} \ell = \ell(\lambda) = \frac{\theta_0}{\rm \pi} (1 - \lambda), \end{equation}

where $\lambda$ is the proportional finger angle parameter. Finally, we have introduced the key non-dimensional parameter, $\epsilon$, by

(2.10)\begin{equation} \epsilon^2 = \frac{4{\rm \pi}^2 \hat{\sigma}}{(1-\lambda)^2}, \end{equation}

where $\hat {\sigma }$ is the modified surface-tension parameter presented in (A9ac) in Appendix A. We consider $\epsilon ^2$ to be a small parameter, corresponding to the small-surface-tension regime, and we will therefore study the problem in the asymptotic limit $\epsilon \rightarrow 0$.

Analyticity of $q\mathrm {e}^{-\mathrm {i} \tau }$ in the upper-half-$\zeta$-plane gives, by the Hilbert transform,

(2.11)\begin{equation} \log q(\xi) = \mathcal{H}[\tau](\xi) \quad \text{where} \ \mathcal{H}[\tau](\xi) = \frac{2}{\rm \pi}{\int\hskip -1,05em -\,}_{-\infty}^0 \frac{\tau(\tilde{\xi}) \, \tilde{\xi}}{\xi^2 - \tilde{\xi}^2} \, {\rm d}{\tilde{\xi}}, \end{equation}

where we have defined the operator, $\mathcal {H}$. Finally, we close the system by integrating the free-surface velocity relationships (2.5). This yields

(2.12)\begin{equation} x(\xi) +\mathrm{i} y(\xi) = \frac{2(1-\lambda)}{\rm \pi}\int_\xi^0 \frac{\mathrm{e}^{\mathrm{i}\tau(\tilde{\xi})}}{q(\tilde{\xi}) P(\tilde{\xi})} \, {\rm d}{\tilde{\xi}}. \end{equation}

Thus, the full system consists of (2.6), (2.11) and (2.12) for the unknowns $(x, y, q, \tau )$ in addition to the boundary conditions

(2.13a,b)$$\begin{gather} x({\pm} \infty) ={-}\infty, \quad y({\pm} \infty) ={\mp} \lambda, \end{gather}$$
(2.13c,d,e)$$\begin{gather}q({\pm} \infty) = 1, \quad \tau(-\infty) = 0, \quad \tau(\infty) ={-}{\rm \pi}. \end{gather}$$

Above, the first set of boundary conditions corresponds to imposing the geometrical constraints while the second set corresponds to the velocity and streamline angle constraints. In total, the equations and boundary conditions are equivalent to those of Ben Amar (Reference Ben Amar1991a). In the next section, we shall examine the zero-surface-tension solutions of these equations.

3. Zero-surface-tension solutions

We first discuss the zero-surface-tension solutions, $(x_0(\xi ), y_0(\xi ), q_0(\xi ), \tau _0(\xi ))$, which correspond to setting $\epsilon = 0$ in the governing equations (2.6), (2.11) and (2.12). We thus approximate $x \sim x_0$, $y \sim y_0$, $q \sim q_0$ and $\tau \sim \tau _0$ in the limit $\epsilon \to 0$. The zero-surface-tension equations are given by

(3.1a)$$\begin{gather} 0 =\frac{Q_0}{1-\lambda} - \frac{2}{\rm \pi} (1+\xi^2){\int\hskip -1,05em -\,}_{-\infty}^0 \frac{K_0(\tilde{\xi}) }{ P(\tilde{\xi})(\xi^2-\tilde{\xi}^2)}\,{\rm d}{\tilde{\xi}}, \end{gather}$$
(3.1b)$$\begin{gather}\log q_0(\xi) = \mathcal{H}[\tau_0](\xi), \end{gather}$$
(3.1c)$$\begin{gather}x_0(\xi) + \mathrm{i} y_0(\xi) ={-} \frac{2(1-\lambda)}{\rm \pi}\int_\zeta^0 \frac{\mathrm{e}^{\mathrm{i}\tau_0(\tilde{\xi})}} {q_0(\tilde{\xi})P(\tilde{\xi})}\,{\rm d}{\tilde{\xi}}, \end{gather}$$

where

(3.2)\begin{equation} K_0(\xi) = \frac{\sin[\tau_0(\xi)]}{q_0(\xi)[r(x_0(\xi))]^2}. \end{equation}

The corresponding boundary conditions are

(3.3ae)\begin{equation} q_0(0) = 0, \quad q_0({\pm} \infty) = 1, \quad \tau_0(-\infty) = 0, \quad \tau_0(0) ={-}\frac{\rm \pi}{2}, \quad \tau_0(\infty) ={-}{\rm \pi}. \end{equation}

As shown by Ben Amar (Reference Ben Amar1991a) and Tu (Reference Tu1991), the leading-order system (3.1) can be re-arranged so as to obtain a Riccati equation

(3.4)\begin{equation} \frac{\mathrm{d}G}{\mathrm{d}\xi} = G^2(\xi) \frac{\xi}{1+\xi^2} + G(\xi) \left(2\ell\frac{\xi}{1+\xi^2} - \frac{1}{\xi}\right) + \ell\left(\frac{1}{(1-\lambda)^2} - 1\right)\frac{1+\xi^2}{\xi}, \end{equation}

with boundary conditions

(3.5)\begin{equation} G({\pm} \infty) = 0, \end{equation}

where the new unknown is defined by

(3.6)\begin{equation} G(\xi) = \frac{\mathrm{e}^{\mathrm{i}\tau_0(\xi)}}{q_0(\xi)} - 1. \end{equation}

It was shown by Ben Amar (Reference Ben Amar1991a) and Tu (Reference Tu1991) that the leading-order (zero-surface-tension) solutions can then be written in terms of the hypergeometric function $F$ (Abramowitz & Stegun Reference Abramowitz and Stegun1972)

(3.7a)$$\begin{gather} x_0 = (1+\xi^2)^{-\ell/2} \times F\left(\frac{\theta_0(2-\lambda)}{2{\rm \pi}}, - \frac{\lambda \theta_0}{2{\rm \pi}}, \frac{1}{2}, \frac{\xi^2}{1+\xi^2}\right), \end{gather}$$
(3.7b)$$\begin{gather}y_0 = \tilde{A}\, \xi\,(1+\xi^2)^{(1-\ell)/2}\times F\left(\frac{1}{2} + \frac{\theta_0 (2-\lambda)}{2{\rm \pi}}, \frac{1}{2} - \frac{\lambda \theta_0}{2{\rm \pi}}, \frac{3}{2}, \frac{\xi^2}{1+\xi^2}\right), \end{gather}$$

where we have also defined the constant

(3.7c)\begin{equation} \tilde{A} = 2 \tan \left(\frac{\lambda \theta_0}{2}\right)\dfrac{\varGamma\left(1-\dfrac{\theta_0(2-\lambda)}{2{\rm \pi}}\right) \varGamma \left(1+ \dfrac{\lambda\theta_0}{2{\rm \pi}}\right)}{\varGamma \left(\dfrac{1}{2} - \dfrac{\theta_0 (2-\lambda)}{2{\rm \pi}}\right)\varGamma \left(\dfrac{1}{2}+\dfrac{\lambda\theta_0}{2{\rm \pi}}\right)}. \end{equation}

By differentiating and rearranging (3.1c) we also obtain expressions for $q_0(\xi )$ and $\tau _0(\xi )$

(3.8a,b)\begin{align} q_0 = \ell\,\frac{\xi}{1+\xi^2}\,\sqrt{\frac{x_0^2+y_0^2}{\left(\dfrac{\mathrm{d}\kern1pt x_0}{\mathrm{d}\xi}\right)^2 + \left(\dfrac{\mathrm{d}y_0}{\mathrm{d}\xi}\right)^2}}, \quad \cos\tau_0 ={-}\frac{q_0}{\ell}\,\frac{1+\xi^2}{\xi}\,\frac{x_0\dfrac{\mathrm{d}\kern1pt x_0}{\mathrm{d}\xi} + y_0 \dfrac{\mathrm{d}y_0}{\mathrm{d}\xi}}{x_0^2 + y_0^2}. \end{align}

In figure 5, we plot example profiles for the leading-order solutions, $(x_0, y_0, q_0, \tau _0)$, along the free surface for parameter values $\theta _0 = 20^\circ$ and $\lambda = 0.6.$ These solutions are generated using (3.7) and (3.8a,b).

Figure 5. Plots of the leading-order (zero-surface-tension) solutions for the four variables $(x_0,y_0,q_0,\tau _0)$ on the free surface, $\zeta = \xi \in \mathbb {R}$. These plots show the solutions for parameter values $\theta _0 = 20^\circ$ and $\lambda = 0.6$ generated using (3.7) and (3.8a,b). The tip of the finger lies at $\xi = 0.$ In figure 3 we plot this solution in the physical plane.

4. Analytic continuation of the free surface

The leading-order profiles, as evaluated on the physical free surface, $\zeta = \xi \in \mathbb {R}$, were shown in the previous section. The analytic continuation of these profiles to the complex plane contains square-root singularities. We shall see that such singularities form one of the key ingredients in the exponential asymptotics procedure of § 6; these points cause the asymptotic expansion to diverge, and will be crucial in determining the eventual selection mechanism. In this section, we discuss the numerical procedure for generating the analytic continuation of the leading-order solutions $(x_0, y_0, q_0, \tau _0)$, as well as the analytic continuation of the governing equations (2.6), (2.11) and (2.12).

4.1. Analytic continuation of the leading-order solutions

In the analytic continuation, we allow the previously real-valued $\xi$ to take complex values. Keeping in mind the potential for confusion with the previously introduced $\zeta$, we write $\xi = \xi _r + \mathrm {i} \xi _c \mapsto \zeta \in \mathbb {C}$. Note that, under this notational choice, $q(\zeta )$ is complex valued within the fluid region, and it is rather the combination $\operatorname {Re}[q\mathrm {e}^{-\mathrm {i}\tau }/(1-\lambda )]$ that is identified with the fluid speed (cf. (2.5)).

In theory, one can replace $\xi$ with $\zeta$, and evaluate the special-function solution (3.7) using standard built-in packages (e.g. Mathematica) to obtain an analytically continued leading-order solution. However, the branch structure of the solutions is complicated and standard software does not easily allow fine-tune control of the branch placement; generation of the full Riemann surface is subsequently difficult. In order to develop the results later in the paper, we must implement a scheme that allows for better control over the generation of the Riemann surface and placement of branches.

For this analytic continuation scheme we first split the Riccati equation (3.4) into its real and imaginary parts

(4.1a)$$\begin{gather} \frac{\mathrm{d}q_0}{\mathrm{d}\zeta} ={-} C(\zeta) \frac{\zeta}{1+\zeta^2}q_0^2 \cos \tau_0 + \frac{q_0}{\zeta} -\ell \frac{\zeta}{1+\zeta^2}\cos\tau_0, \end{gather}$$
(4.1b)$$\begin{gather}\frac{\mathrm{d}\tau_0}{\mathrm{d}\zeta} ={-}C(\zeta)\frac{\zeta}{1+\zeta^2} q_0 \sin\tau_0 + \frac{\ell}{q_0}\frac{\zeta}{1+\zeta^2} \sin\tau_0, \end{gather}$$

where

(4.1c)\begin{equation} C(\zeta) ={-} \ell + \frac{1}{\zeta^2} \left(1+\zeta^2+ \frac{\theta_0\lambda}{\rm \pi}\left(\frac{2-\lambda}{1-\lambda}\right)\right). \end{equation}

Recall that $\ell$ is the rescaled angle parameter introduced in (2.9). Here, we prefer to use $\zeta$, as we are now working with the complexified version of the Riccati equation (3.4).

We may now solve the above system along a chosen parameterised path in the complex $\zeta$-plane by using the exact solution on the free surface as an initial condition. More specifically, we first pre-solve for $(x_0(\zeta ), y_0(\zeta ), q_0(\zeta ), \tau _0(\zeta ))$ on the free surface using (3.7) and (3.8a,b) and setting $\zeta = \xi \in \mathbb {R}$. Then, we parameterise a path into the complex plane that starts on the free surface. For example

(4.2)\begin{equation} \zeta_{path}(s) = \zeta_{IC} + \mathrm{i} s, \quad s\in [0,\infty), \quad \text{where} \ \zeta_{IC}\in \mathbb{R}. \end{equation}

We can then solve (4.1) along the parameterised path using any standard ordinary differential equation integrator (we use MATLAB's ode113 with absolute and relative tolerances set to $10^{-10}$) and the initial condition $\zeta _{IC}$.

The solution consists of eight complexified components (the real and imaginary parts of each of the four variables, $(x_0,y_0,q_0,\tau _0)$). In figure 6, we show the analytically continued surface for one of these components, $\operatorname {Re}(q_0)$. The figure shows one possible path of analytic continuation. By repeatedly solving along a mesh of different paths, the primary Riemann sheet is generated.

Figure 6. Illustrations of the analytic continuation corresponding to $\theta _0 = 20^{\circ }$ and $\lambda = 0.6$ shown via (a) the $(\operatorname {Re} \zeta, \operatorname {Im} \zeta, \operatorname {Re} q_0)$-space; and (b) a top-down view of the $\zeta$-plane. In both, a prototypical path of analytic continuation from the physical free surface is shown with an arrow; branch cuts are shown wavy. Singularities (white circles) lie at $\zeta _1=0.59+1.32\mathrm {i}$ and $\zeta _{C}=0.96\mathrm {i}$; these both correspond to square-root singularities. There is an additional square-root branch point at $\mathrm {i}$. The leading-order solution $q_0$ on the free surface lies on the real $\zeta$ axis.

Next, we find the location of the singularities in the complex plane numerically. Numerical analytic continuation along a closed loop around a branch point demonstrates that the start and end points of the solution differ. Hence, continuation around smaller and smaller loops allows the branch point location to be identified. The singularities in this problem are all square-root singularities. Their locations can be found using the method described above; the singularity strength can be further confirmed by verifying the rate of blow up of the solution as the singularity is approached.

Using this scheme, we find three pairs of complex conjugate singularities, which we denote as $\{\zeta _1, \overline {\zeta _1}\}, \{\zeta _{C}, \overline {\zeta _{C}}\}$ and $\{ \zeta _2,\overline {\zeta _2}\}$. The central singularities, $\{\zeta _{C}, \, \overline {\zeta _{C}}\}$, lie on the imaginary axis and the non-central singularities $\{\zeta _1, \overline {\zeta _1}\}$ and $\{\zeta _2, \overline {\zeta _2}\}$ lie equidistant from the imaginary axis. The conformal map introduces branch points at $\pm \mathrm {i}$ in the $\zeta$-plane. Consequently, the $\{\zeta _1, \overline {\zeta _1}\}$ singularities do not lie on the same Riemann sheet of the complex $\zeta$-plane as the $\{ \zeta _2, \overline {\zeta _2}\}$ singularities. The three groupings of singularity locations are shown in figure 7 panels (ac) for the specific case of $\theta _0 = 20^\circ$ and $\lambda =0.6$.

Figure 7. Plots of the locations of the three complex conjugate pairs of singularities in the $\zeta$-plane: (a) $\{\zeta _1, \overline {\zeta _1}\}$, (b) $\{\zeta _{C}, \overline {\zeta _{C}}\}$, (c) $\{\zeta _2, \overline {\zeta _2}\}$. These plots are for parameter values $\theta _0 = 20^\circ$ and $\lambda = 0.6.$ The branch cuts at $\pm \mathrm {i}$ are chosen to show the relevant branches of the Riemann surface which the singularities lie on. The free surface lies on the real $\zeta$ axis, $\zeta = \xi \in \mathbb {R}.$

We can track the locations of the singularities as we vary the parameters $\theta _0$ (the wedge angle) and $\lambda$ (the proportion of the wedge angle occupied by the finger). If $\theta _0>0$ and $\lambda >0.5$, the non-central singularities $\{\zeta _1, \overline {\zeta _1}, \zeta _2, \overline {\zeta _2}\}$ lie off the imaginary axis. In the limit $\theta _0 \rightarrow 0$ the two singularities $\zeta _1$ and $\zeta _2$ converge to the same point on the imaginary axis, but one will be directly above the other on a separate sheet. A reflection of this behaviour occurs for singularities $\overline {\zeta _1}$ and $\overline {\zeta _2}$ in the lower-half-$\zeta$-plane. The singularity locations in the limit $\theta _0 \to 0$ agree with those found by Chapman (Reference Chapman1999) in the classic Saffman–Taylor problem with parallel channel walls. A summary of the locations of these non-central singularities, for varying values of $\theta _0$ and $\lambda$, is shown in figure 8.

Figure 8. Complex $\zeta$-plane showing how the location of the non-central singularities $\{\zeta _1, \overline {\zeta _1} \}, \{\zeta _2, \overline {\zeta _2} \}$ vary with the parameters $\theta _0$ and $\lambda$. The $\zeta _1$ singularity is shown in the top right quadrant with corresponding $\theta _0$ and $\lambda$ values. The locations of the $\zeta _2,$ $\overline {\zeta _2}$ and $\overline {\zeta _1}$ singularities are shown in the top left, bottom left and bottom right quadrants, respectively.

Finally, the central singularity, $\zeta _{C}$, remains on the imaginary axis between the origin and $\mathrm {i}$ for all parameter values. It moves closer to $\mathrm {i}$ as either $\lambda$ decreases towards 0.5, or as $\theta _0$ decreases towards zero. A reflection of this behaviour occurs for $\overline {\zeta _{C}}$ in the lower-half-plane. Example locations of the singularity $\zeta _{C}$ as $\theta _0$ and $\lambda$ are varied are listed in table 1.

Table 1. The locations of the central singularity, $\zeta _{C}$ (to 4 significant figures) for different values of $\theta _0$ and $\lambda$.

4.2. Analytic continuation for the full problem

If the variable $\zeta$ is analytically continued into the upper-half-$\zeta$-plane, the Bernoulli equation (2.6) becomes

(4.3a)\begin{equation} \epsilon^2 \frac{P(\zeta)}{4} \frac{\partial}{\partial\zeta} \left(r(x) \left[P(\zeta) q(\zeta) \frac{\partial\tau}{\partial\zeta} - \ell \sin\tau(\zeta)\right]\right) = \frac{2}{\rm \pi} \int_{-\infty}^0 \frac{K(\tilde{\zeta}) \,\tilde{\zeta}}{\zeta^2-\tilde{\zeta}^2}\,{\rm d}{\tilde{\zeta}} + \mathrm{i} K(\zeta). \end{equation}

Here, the principal-value integral becomes a normal integral and can then be combined with the $Q_0$ term in (2.6) to simplify the right-hand side. Recall that $P(\zeta ), r(x(\zeta )),\ell, K(\zeta )$ were introduced in (2.7ac) for convenience.

To analytically continue the boundary-integral equation (2.11) we must consider the complexification of the Hilbert transform

(4.3b)\begin{equation} \mathcal{H}[\tau] = \hat{\mathcal{H}}[\tau] + \mathrm{i} \tau, \end{equation}

where $\hat {\mathcal {H}}[\tau ]$ is the complex-valued Hilbert transform

(4.3c)\begin{equation} \hat{\mathcal{H}}[\tau](\xi) = \frac{2}{\rm \pi}\int_{-\infty}^0 \frac{\tau(\tilde{\xi}) \, \tilde{\xi}}{\xi^2 - \tilde{\xi}^2} \, {\rm d}{\tilde{\xi}}. \end{equation}

The boundary-integral equation (2.11) then becomes

(4.3d)\begin{equation} \log q(\zeta) - \mathrm{i} \tau(\zeta) = \hat{\mathcal{H}}[\tau](\zeta), \end{equation}

and hence the integrated equation for the surface position (2.12) becomes

(4.3e)\begin{equation} x(\zeta)+\mathrm{i} y(\zeta) = \frac{2(1-\lambda)}{\rm \pi}\int_\zeta^0 \frac{\exp\left(-\hat{\mathcal{H}}[\tau](\tilde{\zeta})\right)}{P(\tilde{\zeta})}\,{\rm d}{\tilde{\zeta}}. \end{equation}

This results in a set of analytically continued governing equations (4.3) that hold in the upper-half-$\zeta$-plane.

5. Exponential asymptotics

Our procedure for the exponential asymptotic analysis follows similar work by Tanveer (Reference Tanveer1987) and Chapman (Reference Chapman1999), using the methodology established in Chapman et al. (Reference Chapman, King and Adams1998). In essence, the goal is to derive the behaviour of the late terms in the asymptotic series. After, in § 6, these late terms are used to study the exponentially small terms via the Stokes-line switching.

As $\epsilon \rightarrow 0$, we expand the dependent variables as

(5.1) \begin{equation} \begin{gathered} x(\zeta) \sim \sum_{n=0}^\infty \epsilon^{2n}x_n(\zeta), \quad y(\zeta) \sim \sum_{n=0}^\infty \epsilon^{2n}y_n(\zeta), \\ q(\zeta) \sim \sum_{n=0}^\infty \epsilon^{2n}q_n(\zeta), \quad \tau(\zeta) \sim \sum_{n=0}^\infty \epsilon^{2n}\tau_n(\zeta). \end{gathered}\end{equation}

We substitute the above into the analytically continued governing equations (4.3) and this yields, at ${O}(\epsilon ^{2n})$ for Bernoulli's equation (4.3a),

(5.2a) \begin{align} &\frac{P}{4}\frac{\partial}{\partial \zeta}\left(r(x_0)\left[Pq_{n-1} \frac{\partial \tau_0}{\partial \zeta} + Pq_0 \frac{\partial \tau_{n-1}}{\partial \zeta} - \ell \tau_{n-1}\cos\tau_0\right.\right.\nonumber\\ &\qquad \qquad \qquad \ \left.\left.-\frac{\theta_0 x_{n-1}}{2} \left(Pq_0 \frac{\partial \tau_0}{\partial \zeta} - \ell \sin\tau_0 \right) + \ldots \right] \right) \nonumber\\ &\qquad\qquad\qquad\qquad \qquad = \ldots + \frac{\mathrm{i}}{r^2(x_0)\,q_0}\left(\tau_n\cos\tau_0 - \frac{q_n}{q_0}\sin\tau_0 +\ldots\right), \end{align}

while for (4.3d) and (4.3e) we have

(5.2b)$$\begin{gather} \frac{q_n}{q_0} + \ldots - \mathrm{i}\tau_n = \hat{\mathcal{H}}[\tau_n], \end{gather}$$
(5.2c)$$\begin{gather}x_n + \mathrm{i} y_n = \frac{2(1-\lambda)}{\rm \pi}\int_\zeta^0 \frac{\exp(-\hat{\mathcal{H}}[\tau_0]) \hat{\mathcal{H}}[\tau_n] + \ldots}{P}\,{\rm d}{\tilde{\zeta}}. \end{gather}$$

At these later orders we see that the $n$th terms in the power series for $q$ and $\tau$ are obtained by differentiating the $(n-1)$th terms twice. Any singularities in the leading-order solution will grow in strength with each successive differentiation. This means that later terms in the power series will have singularities at the same locations as earlier terms, but with increasing strength. We therefore follow the method of Chapman et al. (Reference Chapman, King and Adams1998) and predict a factorial-over-power form for the late-order terms

(5.3)\begin{equation} q_n(\zeta) \sim \frac{Q(\zeta) \varGamma (2n+\gamma)}{\chi(\zeta)^{2n+\gamma}}, \quad \tau_n(\zeta) \sim \frac{\varTheta(\zeta)\varGamma(2n+\gamma)}{\chi(\zeta)^{2n+\gamma}}, \quad \text{as} \ n \rightarrow \infty, \end{equation}

where $Q$ and $\varTheta$ are prefactors and $\chi$ is a singulant function, which is zero at the singularity. The singulant ensures that each series term has singularities with the correct locations and $\gamma$ ensures they have the correct strengths. The gamma function (Abramowitz & Stegun Reference Abramowitz and Stegun1972) is a consequence of the factorial behaviour caused by repeated differentiation. The late-order terms are a sum of such factorial-over-power terms – one associated with each distinct complex singularity. Note that the prototypical factorial-over-power divergence of singular asymptotic expansions is a consequence of Darboux's theorem (cf. p. 4 of Dingle Reference Dingle1973; Crew & Trinh Reference Crew and Trinh2023).

In the limit $n\to \infty$, the behaviour of the asymptotic expansion will be dictated by the divergence caused by the singularities driving (5.3) (Chapman et al. Reference Chapman, King and Adams1998). From § 4 we know that the singularities lie in the complex plane away from the free surface. The complex Hilbert transform, $\hat {\mathcal {H}}[\tau _n]$, involves the integrand evaluated along the free surface. Once the ansatzes for $q_n$ and $\tau _n$ via (5.3) are substituted into the Hilbert transform, we may observe that the contribution of $\hat {\mathcal {H}}[\tau _n]$ will be subdominant in the limit that $n \rightarrow \infty$, compared with $q_n$ and $\tau _n$. This follows from the increasing nature of $|\chi |$ along Stokes lines, as explained on p. 526 of Chapman (Reference Chapman1999). The combination of $x_n+\mathrm {i} y_n$ will also be subdominant in this limit, although the individual components may still diverge. We will assume in this analysis that the individual components do not diverge, and this assumption will be validated a posteriori, as is done in similar asymptotic analyses of boundary-integral problems (Shelton & Trinh Reference Shelton and Trinh2022).

Using these assumptions gives the dominant behaviour from the boundary-integral equation (5.2b)

(5.4)\begin{equation} \frac{q_n}{q_0} - \frac{q_1q_{n-1}}{q_0^2} \sim \mathrm{i}\tau_n \quad \text{as} \ n \rightarrow \infty. \end{equation}

The Bernoulli equation (5.2a) can then be simplified to

(5.5)\begin{equation} \frac{q_{n-1}''}{q_n} + \left[\frac{\left(r(x_0)P\right)'}{r(x_0)P} - \frac{q_0'}{q_0} + \mathrm{i}\tau_0' - \ell \frac{\cos\tau_0}{Pq_0}\right]\frac{q_{n-1}'}{q_n} = \frac{4(\sin\tau_0 + \mathrm{i} \cos\tau_0)}{P^2q_0^2r(x_0)^3}. \end{equation}

Here, and for the rest of the paper, we use primes ($'$) to denote differentiation with respect to $\zeta$. Substitution of the factorial-over-power ansatz (5.3) and matching terms in the limit that $n \rightarrow \infty$ gives at leading order an equation for the singulant, $\chi$,

(5.6)\begin{equation} (\chi')^2 = \frac{4(\sin\tau_0 + \mathrm{i}\cos\tau_0)}{P^2q_0^2r(x_0)^3}, \end{equation}

which can be solved to obtain

(5.7)\begin{equation} \chi ={-}\int_{\zeta_*}^\zeta \frac{2\sqrt{\sin\tau_0 + \mathrm{i}\cos\tau_0}}{P\,q_0\,r(x_0)^{{3}/{2}}} \,{\rm d}{\tilde{\zeta}}. \end{equation}

In this expression, $\zeta _*$ is the singularity location and the integration limits are chosen so that $\chi (\zeta _*) = 0$. There will be a singulant function $\chi$ for each complex singularity. The negative sign is selected when taking the square root so that the Stokes line intersects the free surface, which lies on the real $\zeta$ axis.

Matching terms in (5.5) at the next order as $n \rightarrow \infty$ shows that $\gamma$ is constant, and at the following order we obtain

(5.8)\begin{equation} \frac{Q'}{Q} ={-} \frac{1}{2}\left[\frac{\chi''}{\chi'} + \frac{(r(x_0)P)'}{r(x_0)P} - \frac{q_0'}{q_0} + \mathrm{i}\tau_0' - \ell \frac{\cos\tau_0}{Pq_0}\right], \end{equation}

which can be solved to give an expression for the prefactor

(5.9)\begin{equation} Q = \varLambda \mathrm{i} \left(\frac{q_0\, \mathrm{e}^{-\mathrm{i}\tau_0} \exp \left[\displaystyle\int\nolimits_0^\zeta \ell \dfrac{\cos\tau_0}{Pq_0} {\rm d}{\tilde{\zeta}}\right]}{\chi' r(x_0)P}\right)^{{1}/{2}}, \end{equation}

where $\varLambda$ is a constant that remains to be determined. From the boundary-integral equation (5.4) we find

(5.10)\begin{equation} \varTheta ={-}\mathrm{i}\frac{Q}{q_0}. \end{equation}

By definition, $\ell \rightarrow 0$ (see (2.9)) and $r(x_0) \rightarrow 1$ (see (2.7ac)) in the limit $\theta _0 \rightarrow 0$. This means that these results are consistent with those found by Chapman (Reference Chapman1999) for the Saffman–Taylor problem in a channel geometry.

The late-order series terms in the divergent power series solution for $q(\zeta )$ therefore have the form

(5.11)\begin{equation} q_n \sim \varLambda \mathrm{i} \left(\frac{q_0 \mathrm{e}^{-\mathrm{i}\tau_0} \exp \left[\displaystyle\int\nolimits_0^\zeta \ell \dfrac{\cos\tau_0}{Pq_0} {\rm d}{\tilde{\zeta}}\right]}{\chi' r(x_0)P}\right)^{{1}/{2}} \frac{\varGamma(2n+ \gamma)}{\chi^{2n+\gamma}}, \quad \text{as} \ n\rightarrow \infty. \end{equation}

6. Optimal truncation and Stokes lines

In the previous section, we noted that complex singularities cause the power series to become divergent and so they will need to be truncated. We will truncate the divergent power series at some to be determined optimal point $N$, and introduce the remainder terms

(6.1a)$$\begin{gather} x = \sum_{n=0}^{N-1} \epsilon^{2n} x_n + R_x, \quad y = \sum_{n=0}^{N-1} \epsilon^{2n}y_n + R_y, \end{gather}$$
(6.1b)$$\begin{gather}q = \sum_{n=0}^{N-1} \epsilon^{2n} q_n + R_q, \quad \tau = \sum_{n=0}^{N-1} \epsilon^{2n} \tau_n + R_\tau. \end{gather}$$

We can substitute these into the governing equations (4.3). Given that the first $N$ orders must exactly satisfy the relationships in (5.2), we derive leading-order relationships between the remainder terms.

From the boundary-integral equation (5.4), we derive at leading order as $\epsilon \rightarrow 0$ that

(6.2)\begin{equation} \frac{R_q}{q_0} \sim \mathrm{i} R_\tau. \end{equation}

Using this in the $x+\mathrm {i} y$ equation (5.2c) we see that $R_x$ and $R_y$ are subdominant compared with the leading orders of $R_q$ and $R_\tau$ in the limit $\epsilon \rightarrow 0$.

We substitute these relationships into the Bernoulli equation (4.3a) and derive a single ordinary differential equation for $R_q$, which we now rename as $R_N$ for consistency with other works, including Chapman (Reference Chapman1999). This ordinary differential equation, known as the remainder equation, reduces to

(6.3)\begin{equation} \epsilon^2R_N'' + \epsilon^2 \left[\frac{(r(x_0)P)'}{r(x_0)P} -\frac{q_0'}{q_0}+ \mathrm{i} \tau_0' - \ell \frac{1}{P q_0}\cos\tau_0\right]R_N' = (\chi')^2(R_N-\epsilon^{2N}q_N), \end{equation}

as $\epsilon \rightarrow 0$, where the terms that do not appear at the leading or second order in this limit have been omitted. Changing the independent variable to $\chi$ (primes will continue to denote differentiation with respect to $\zeta$) gives

(6.4)\begin{equation} \epsilon^2 \frac{\mathrm{d}^2R_N}{\mathrm{d}\chi^2} + \epsilon^2 \frac{1}{\chi'} \left[\frac{\chi''}{\chi'} + \frac{(r(x_0)P)'}{r(x_0)P} -\frac{q_0'}{q_0}+ \mathrm{i} \tau_0' - \ell \frac{1}{P q_0}\cos\tau_0\right]\frac{\mathrm{d} R_N}{\mathrm{d}\chi} - R_N ={-} \epsilon^{2N}q_N. \end{equation}

Using the definition for the prefactor (5.9), this can be simplified to

(6.5)\begin{equation} \epsilon^2 \frac{\mathrm{d}^2R_N}{\mathrm{d}\chi^2} -2\epsilon^2\frac{Q'}{Q\chi'}\frac{\mathrm{d} R_N}{\mathrm{d}\chi} - R_N ={-}\epsilon^{2N}q_N. \end{equation}

To solve (6.5) we pose a Liouville–Green or WKBJ-style ansatz for the form of the remainder given by

(6.6)\begin{equation} R_N = B(\chi){\rm e}^{{b(\chi)}/{\epsilon}}. \end{equation}

Then equating the coefficients at different powers of $\epsilon$ for the homogeneous version of the remainder equation (6.5), we find that $b = \pm \chi + \textrm {const.}$ and $B\sim Q.$ The arbitrary constant in $b$ is equivalent to multiplying the entire expression (6.6) by an arbitrary constant that is not determined by the WKBJ analysis.

To solve the full inhomogeneous remainder equation (6.5) we apply the method of variation of parameters and permit the arbitrary constant to vary in $\zeta$. This quantity is known as the Stokes multiplier, and we denote it by $A(\chi )$. The remainder becomes

(6.7)\begin{equation} R_N = A Q \mathrm{e}^{-{\chi}/{\epsilon}}, \end{equation}

where the negative sign in the exponent ensures the remainder is exponentially small. The inhomogeneous equation (6.5) gives

(6.8)\begin{equation} {-}2\epsilon \frac{\mathrm{d}A}{\mathrm{d}\chi} Q \mathrm{e}^{-{\chi}/{\epsilon}} ={-}\epsilon^{2N} \frac{Q\varGamma(2N+\gamma)}{\chi^{2N+\gamma}}, \end{equation}

where we have substituted in the late-order expression for $q_N$ from (5.3).

The next step involves truncating the series at an optimal point. We define this optimal truncation point, $N$, to be where successive terms in the divergent series are approximately equal in magnitude, so

(6.9)\begin{equation} \left| \frac{\epsilon^{2N+2} q_{N+1}}{\epsilon^{2N}q_N}\right| \approx 1. \end{equation}

This condition gives $N \approx {|\chi |}/{2\epsilon }$. As $N$ must be an integer we let $N = {|\chi |}/{2\epsilon } +\beta$, where $\beta$ is bounded as $\epsilon \rightarrow 0.$ This form motivates the transformation to polar coordinates, so we define $\chi = |\chi |\mathrm {e}^{\mathrm {i}\eta }.$

By the chain rule we have

(6.10)\begin{equation} \frac{\mathrm{d}A}{\mathrm{d}\eta} = \frac{\mathrm{d}A}{\mathrm{d}\chi}\mathrm{i}\chi. \end{equation}

We substitute the optimal value of $N$ into (6.8) and note that $N$ is large, which allows us to use Stirling's formula (Abramowitz & Stegun Reference Abramowitz and Stegun1972) to approximate $N!$ in the large $N$ limit. Then using (6.10) we obtain

(6.11)\begin{equation} \frac{\mathrm{d}A}{\mathrm{d}\eta} \sim \frac{\mathrm{i} \epsilon^{2N-1}}{2}\frac{\exp\left(\dfrac{|\chi|}{\epsilon}\mathrm{e}^{\mathrm{i}\eta}\right) \sqrt{2{\rm \pi}} \left(2N+\gamma\right)^{2N+\gamma-{1}/{2}}}{\mathrm{e}^{2N+\gamma}|\chi|^{2N+\gamma-1}\mathrm{e}^{\mathrm{i}\eta(2N+\gamma-1)}}. \end{equation}

We substitute in the optimal value of $N$ to give

(6.12)\begin{align} \frac{\mathrm{d}A}{\mathrm{d}\eta} \sim \frac{\mathrm{i}\sqrt{\rm \pi}}{\sqrt{2}c}\frac{|\chi|^{{1}/{2}}}{\epsilon^{\gamma+{1}/{2}}}\frac{\mathrm{e}^{2\beta+\gamma}}{ \exp\left(\mathrm{i}\eta\left( \dfrac{|\chi|}{\epsilon}+2\beta+\gamma\right)\right)}\frac{\exp\left(\dfrac{|\chi|}{\epsilon} \mathrm{e}^{\mathrm{i}\eta}\right)}{\exp\left(\dfrac{|\chi|}{\epsilon} + 2\beta+\gamma\right)} \quad \text{as} \ \epsilon \rightarrow 0. \end{align}

The right-hand side of (6.12) is exponentially small in $\epsilon$ unless $\eta = 2{\rm \pi} k,$ $k \in \mathbb {Z},$ in which case it will be algebraic in the limit that $\epsilon \rightarrow 0.$ Therefore, the greatest change in the Stokes multiplier $A$ will occur across curves, or Stokes lines, on which $\eta = 2{\rm \pi} k$, and so

(6.13a,b)\begin{equation} \mathrm{Re}(\chi)>0, \quad \mathrm{Im}(\chi) = 0. \end{equation}

This recovers the classic result of Dingle (Reference Dingle1973). We compute the Stokes lines numerically, with the results shown in figure 9.

Figure 9. Complex $\zeta$-plane showing the Stokes lines emanating from the singularities (circles) and intersecting the free surface (real axis). In this figure, the wedge angle is $\theta _0 = 20^\circ$ and $\lambda = 0.6$. Branch cuts (shown wavy) lie up and down the imaginary axis from $\pm \mathrm {i}$. Stokes lines are shown with dashed when they lie on a different Riemann sheet to the free surface. The three points where the Stokes lines intersect the real axis are labelled $S_1,$ $S_{C}$ and $S_2$.

When Stokes lines intersect the free surface (the real axis in the $\zeta$-plane) then an exponentially small term will be smoothly switched on across this intersection point in the solution. We can see from figure 9 that there are three points on the free surface where such exponentially small terms will be switched on.

To find the jump in solution behaviour across a Stokes line, we rescale $\eta$ and $A$ and consider the behaviour in the neighbourhood of the Stokes line. We apply the rescaling

(6.14a,b)\begin{equation} \eta = \epsilon^2 \tilde{\eta} \quad \text{and} \quad A = \tilde{A}\epsilon^{-\gamma-1}. \end{equation}

We let $k = 0$ and then (6.12) becomes

(6.15)\begin{equation} \frac{\mathrm{d}\tilde{A}}{\mathrm{d}\tilde{\eta}} \sim \frac{\mathrm{i}\sqrt{\rm \pi}}{\sqrt{2}} |\chi|^{{1}/{2}}\exp\left(-\frac{1}{2}|\chi| \tilde{\eta}^2\right) \quad \text{as}\ \epsilon\rightarrow 0. \end{equation}

Integrating this gives

(6.16)\begin{equation} \tilde{A} = \text{const.} + \mathrm{i}\sqrt{\rm \pi} \int_{-\infty}^{\tilde{\eta}\sqrt{{|\chi|}/{2}}} \mathrm{e}^{-\tilde{\xi}^2}\,{\rm d}{\tilde{\xi}}, \end{equation}

and hence the jump in the Stokes multiplier across the Stokes line is

(6.17)\begin{equation} \lim_{\eta \rightarrow 0+}A(\eta) - \lim_{\eta \rightarrow 0-}A(\eta) = \frac{\mathrm{i}\sqrt{\rm \pi}}{\epsilon^{\gamma+1}} \int_{-\infty}^\infty \mathrm{e}^{-\tilde{\xi}^2}{\rm d}{\tilde{\xi}} = \frac{\mathrm{i}{\rm \pi}}{ \epsilon^{\gamma+1}}. \end{equation}

That is, upon crossing a Stokes line at its intersection with the real axis, a contribution in the $q$ variable is switched on. This has the form

(6.18)\begin{equation} \frac{\mathrm{i} {\rm \pi}}{ \epsilon^{\gamma+1}}Q\exp\left(-\frac{\chi}{\epsilon}\right) + \text{c.c.}, \end{equation}

where c.c. represents the complex conjugate expression, which arises from the singularity located at the complex conjugate location in the complex plane (figure 9). Similarly, by (5.10), the contribution that is switched on in the $\tau$ variable has the form

(6.19)\begin{equation} \frac{\rm \pi}{ \epsilon^{\gamma+1}} \frac{Q}{q_0}\exp\left(-\frac{\chi}{\epsilon}\right) + \text{c.c.}. \end{equation}

In the next section, we show how these terms switched on across Stokes lines will determine the selection mechanism.

7. Selection mechanism

In figure 9 we see that there are three points on the free surface ($\zeta = \xi \in \mathbb {R}$) that are intersected by Stokes lines. We will continue the notation from the singularities (introduced in § 4.1) and use the subscripts ‘$1$’, ‘${C}$’ and ‘$2$’ to label the Stokes lines, i.e. each subscript matches the respective label of the associated Stokes-line singularity. We denote the three intersection points as $S_1$, $S_{C}$ and $S_2$. Each Stokes line will have different corresponding values for $\varLambda,$ $\gamma$ and $\chi$, which will also be labelled with the same subscripts.

At each intersection point an exponentially small asymptotic contribution of the form (6.19) will be switched on or off. The far-field conditions $\tau (-\infty ) = 0$ and $\tau ( \infty ) = -{\rm \pi}$ imply there are no exponentially small contributions present on the free surface as $\tau \to \pm \infty$, so the exponential terms must only be present in the region of the free surface between the Stokes-line intersection points; that is, in the range $\zeta =\xi \in (S_1, S_2)$.

To check for existence of a solution, we can check that the conditions in both far fields are satisfied. However, due to the symmetry of the problem, it is equivalent to consider half the domain and impose symmetry conditions at the origin: $q(0) = 0,$ $\tau (0) = -{\rm \pi} /2$. A similar symmetry condition is used in Chapman et al. (Reference Chapman, Dallaston, Kalliadasis, Trinh and Witelski2023). Enforcing the condition on $\tau$ implies

(7.1)\begin{equation} \tau_0(0) ={-}\frac{\rm \pi}{2},\quad \tau_n(0) = 0 \quad \mathrm{for} \ n\geqslant 1, \quad R_\tau(0) = 0. \end{equation}

Requiring that the solution satisfies the far field conditions as $\xi \to -\infty$ and the conditions at the origin from (7.1) will result in a selection condition that must be satisfied for solutions to exist.

At the point $S_1$, a contribution of the form (6.19) is switched on as the Stokes line is crossed from left to right. And then to reach the origin we switch on half of another contribution of the form (6.19) (as we have only crossed half the boundary layer about the ‘${C}$’ Stokes line). That is, the integral in (6.16) will only range from $-\infty$ to $\tilde {\eta } = 0.$ The symmetry condition from (7.1) says that the sum of these contributions must be zero at the origin. This means

(7.2) \begin{align} &{\rm \pi} \epsilon^{-\gamma_1-1}\frac{\varLambda_1}{\sqrt{2}} \mathrm{e}^{-\mathrm{i}{\tau_0}/{4}}I(\zeta)r(x_0)^{{1}/{4}}\mathrm{e}^{3\mathrm{i}{\rm \pi}/{8}} \exp\left(-\frac{\chi_1}{\epsilon}\right) + \mathrm{c. c.}\nonumber\\ &\qquad +\left.\frac{1}{2}{\rm \pi} \epsilon^{-\gamma_{\text{C}}-1}\frac{\varLambda_C}{\sqrt{2}} \mathrm{e}^{-\mathrm{i}{\tau_0}/{4}}I(\zeta)r(x_0)^{{1}/{4}}\mathrm{e}^{3\mathrm{i}{\rm \pi}/{8}} \exp\left(-\frac{\chi_{\text{C}}}{\epsilon}\right) + \mathrm{c.c.}\right|_{\zeta=0} = 0, \end{align}

where $I(\zeta ) := \exp [\int _0^\zeta \ell ({\cos \tau _0}/{2Pq_0})\,\textrm {d}{\tilde {\zeta }}].$ We simplify this expression using the conditions at the origin, $\tau _0(0) = -{{\rm \pi} }/{2}$ and $x_0(0) = 0$. Noting that the origin lies on the ‘${C}$’ Stokes line, we use the definition for a Stokes line obtained in (6.13a,b), i.e. $\mathrm {Im}(\chi _{C}(0)) = 0$. Then

(7.3) \begin{align} &\frac{|\varLambda_1|}{|\varLambda_{C}|} \epsilon^{-\gamma_1+\gamma_{C}}e^{{-\mathrm{Re}(\chi_1(0))/\epsilon+\mathrm{Re}(\chi_{C}(0))/\epsilon}} \cos\left(\arg(\varLambda_1) + \frac{\rm \pi}{2} - \frac{\mathrm{Im}(\chi_1(0))}{\epsilon}\right)\nonumber\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad + \frac{1}{2}\cos\left(\arg(\varLambda_{C}) + \frac{\rm \pi}{2}\right) = 0. \end{align}

By the symmetry of the problem, this condition could also be obtained by calculating the values that appear across the ‘$2$’ and ‘${C}$’ Stokes lines. Note that, for this problem, we always have symmetry of $q_0$ at the origin and so $q_0(0) = q_0'(0) = 0$. This means that the possible selection condition, $q'(0) = 0$, is automatically satisfied. We therefore use the $\tau$ variable to derive the selection mechanism.

The constants $\gamma _1$ and $\gamma _{C}$ are derived by an inner matching of the local singularity strengths at lower orders in the asymptotic series. The details of the derivation can be found in Appendix B. We find that $\gamma _1 = \gamma _{C} = {1}/{14}$ and thus we can simplify the selection condition to

(7.4) \begin{align} &\frac{|\varLambda_1|}{|\varLambda_{C}|} {\rm e}^{-\mathrm{Re}(\chi_1(0))/\epsilon+\mathrm{Re}(\chi_{C}(0))/\epsilon} \cos\left(\arg(\varLambda_1) + \frac{\rm \pi}{2} - \frac{\mathrm{Im}(\chi_1(0))}{\epsilon}\right)\nonumber\\ &\qquad\qquad\qquad\qquad\qquad\quad\qquad\qquad\qquad\qquad \ + \frac{1}{2}\cos\left(\arg(\varLambda_{C}) + \frac{\rm \pi}{2}\right) = 0. \end{align}

Above, the constants $\varLambda _1$ and $\varLambda _{C}$ are derived in Appendix C by using an asymptotic matching process to connect the late-order term behaviour (5.3) with a local expansion of the solution in a neighbourhood of the relevant singularity. We perform this matching numerically, and find that the values of $\varLambda _1$ and $\varLambda _{C}$ depend on both the parameters $\theta _0$ and $\lambda$. For example, when $\theta _0 = 20^\circ$ and $\lambda = 0.6$ then $\varLambda _1 \approx -0.48-0.15\mathrm {i}$ and $\varLambda _{C} \approx -0.49-0.24\mathrm {i}$.

8. Results

For each $\epsilon$ value we now use a numerical root finding scheme to find the corresponding family of $\lambda$ values that satisfy the selection condition (7.4). The dependence on $\lambda$ in the selection condition (7.4) appears through the constants $\varLambda _1$ and $\varLambda _{C}$, as shown in Appendix C, and also through $\chi _1(0)$ and $\chi _{C}(0)$, as given in (5.7). The selected solution families, for different values of $\theta _0$, are shown in figures 10(a)–10(d).

Figure 10. Bifurcation diagrams produced using (7.4). The first ten selected values of $\lambda$ as functions of $\epsilon$ are shown. Panel (a) shows the bifurcation diagram for the channel geometry. The other panels (bd) show the bifurcation diagram for different values of the wedge angle $\theta _0$.

The figure shows that the bifurcation diagram differs qualitatively between the channel geometry with $\theta _0 = 0^\circ$ and the general wedge geometry with $\theta _0>0$. For the channel geometry, each branch, $\lambda _n(\epsilon )$, continues to exist as $\epsilon \rightarrow 0$. However, for the wedge geometry, we see that the branches merge and disappear in pairs and no solutions will exist when $\epsilon =0$. The lowest two branches, $\lambda _1(\epsilon )$ and $\lambda _2(\epsilon )$, merge at the greatest value of $\epsilon$ and then merging continues to happen between higher branches as $\epsilon \rightarrow 0.$ Each time such a merge happens, the solutions corresponding to those two branches will cease to exist. In practice, we hypothesise that the existence of the self-similar solution is lost due to a tip-splitting instability which occurs when the flow rate or surface tension reaches a critical value. In figure 10 we see that, for smaller wedge angles, the tip-splitting will occur at a smaller value of the effective surface tension $\epsilon$. In the limit $\theta _0 \rightarrow 0$, the surface-tension value at which the branch merging occurs approaches zero.

As noted in e.g. Ben Amar (Reference Ben Amar1991b), it is expected that, as is the case in the channel geometry, only the $\lambda _1$ branch is stable. Then the merging of the $\lambda _1$ and $\lambda _2$ branches as $\epsilon \to 0$ is associated with a loss of stable solutions, resulting in a tip-splitting instability.

The key observation is that the relative size of the terms in the selection condition (7.4) provides an approximate criterion for the existence of solutions. We recall that, as $\epsilon \to 0$, the classic Saffman–Taylor fingering in a channel involves $\lambda = \lambda _i \to \lambda _{lim} = 0.5$ for all indices $i$. It can be asked whether a similar limiting value of $\lambda$ is reached as $\epsilon \to 0$ and $i \to \infty$ for general $\theta _0$.

In figure 11 we plot the real parts of $\chi _1(0)$ and $\chi _{C}(0)$ for the $\theta _0 = 20^\circ$ case. As indicated in the figure, for $\lambda > \lambda _{lim}$, then $0 < \operatorname {Re}\chi _1(0) < \operatorname {Re}\chi _{C}(0)$. Therefore, from the selection condition, the contribution from the central Stokes line is exponentially dominated by the contribution from the ‘1’ Stokes line(s).

Figure 11. (a) Relative sizes of the real parts of $\chi _1(0)$ and $\chi _{C}(0)$ for different values of $\lambda$ in the $\theta _0 = 20^\circ$ case. The $\lambda$ value where $\mathrm {Re}(\chi _{C}(0)) =\mathrm {Re}(\chi _1(0))$ is denoted $\lambda _{lim}$. Above this value the selection condition (7.4) can be satisfied. (b) The limiting $\lambda$ value, $\lambda _{lim}$, in the small-surface-tension limit as a function of the wedge angle.

When $\mathrm {Re}(\chi )$ is smaller, the corresponding $\exp ({-\mathrm {Re}(\chi )/\epsilon })$ term dominates in the selection condition. As figure 11 shows, in the small $\epsilon$ limit, for smaller $\lambda$ values, the contribution from the ‘$C$’ Stokes line dominates; while for larger $\lambda$ values, the contribution from the ‘$1$’ Stokes line dominates. When the contribution from the ‘1’ Stokes line dominates, then the exponentially small oscillations can be switched off by the symmetric ‘2’ Stokes line and solutions will exist, similarly to the channel case. However, when the ‘C’ Stokes line dominates, then in the small $\epsilon$ limit it becomes impossible to satisfy the selection condition and no solutions will exist. The $\lambda$ value at which $\mathrm {Re}(\chi _1(0)) = \mathrm {Re}(\chi _{C}(0))$ gives the limiting $\lambda$ value, $\lambda _{lim}$. The $\lambda _{lim}$ values as a function of $\theta _0$ are shown in figure 11.

9. Conclusion

We have used exponential asymptotics to derive the selection mechanism in the small-surface-tension limit for the divergent Saffman–Taylor fingers in a wedge of an arbitrary angle. The selection mechanism is based on the requirement for the exponentially small contributions, which are switched on at Stokes lines, to satisfy the far-field boundary conditions. This Stokes-line structure in the complex plane relies on an understanding of the singularities in the analytic continuation of the leading-order solution. Here, we must obtain this using a numerical scheme for the analytic continuation. We find the countable family of $\lambda (\epsilon )$ values and associated selected solutions that satisfy the selection condition. The bifurcation diagrams are plotted in figure 10 and show the merging of branches of $\lambda$ values as the surface-tension parameter, $\epsilon$, is decreased. We hypothesise that the loss of existence of solutions through this branch merging relates to the tip-splitting instability observed in the time-dependent flows in a circular geometry.

10. Discussion

There are a number of interesting extensions to the work in this paper in the fields of exponential asymptotics and Hele-Shaw problems.

In order to verify the correctness of the asymptotic predictions, we have used some of the early numerical work of Ben Amar (Reference Ben Amar1991b). As it turns out, there are a few aspects of numerical computations of the governing equations (2.6), (2.11) and (2.12) that render it more challenging. In a forthcoming work, we will present a specialised scheme that solves for the free-surface profile of the fingers.

One of the main difficulties with the analysis in this paper arises from the form of the leading-order solutions. Traditionally, in problems considered with exponential asymptotics, the leading-order solution will have a simple closed form; then, locations and strengths of the singularities in the complex plane are easily identified. For this problem, the leading-order solution is written in terms of special functions (cf. (3.7) and (3.8a,b)); however, the singularities of these expressions are not so easily obtained. In practice, we use a numerical scheme to analytically continue the leading-order solutions into the complex $\zeta$-plane and search for singularities. This numerical scheme could equally be used to analytically continue a fully numerical leading-order solution into the complex plane. Similar numerical analytic continuation schemes can be found in Chandler & Trinh (Reference Chandler and Trinh2018) and Crew & Trinh (Reference Crew and Trinh2016). With the advancement of such schemes, in the field of exponential asymptotics we are no longer restricted only to problems with a simple analytic leading-order solution.

To complete the understanding of the tip-splitting behaviour in the circular geometry it will be necessary to check the stability of the branches plotted in the bifurcation diagram (figure 10). For the Saffman–Taylor problem in the channel geometry, Tanveer (Reference Tanveer1987) showed that only the lowest branch, $\lambda _1$, is linearly stable. We conjecture that the same will be true for the wedge geometry and so the merging of branches $\lambda _1$ and $\lambda _2$ will correspond to the loss of existence of any solutions which could be observed physically. We have written a time-dependent code for the circular geometry based on the numerical method described in Dallaston & McCue (Reference Dallaston and McCue2014) and we hope that this will provide an insight into the stability of the branches observed in the bifurcation diagram for the wedge problem.

Finally, there are many recent experimental and numerical results for different Hele-Shaw problems, some of which appear to have similar selected branches of permitted solutions (Gaillard et al. Reference Gaillard, Keeler, Le Lay, Lemoult, Thompson, Hazel and Juel2021; Keeler et al. Reference Keeler, Gaillard, Lawless, Thompson, Juel and Hazel2022). We expect that the techniques used here can also be applied to these problems to derive selection conditions, governed by exponentially small components of the solutions. In particular, the ability to do exponential asymptotics without a simple analytic leading-order solution will be necessary to make progress in these problems.

Acknowledgements

We especially thank J. Chapman (Oxford) for helping us navigate the formulation of this problem and getting us started on the analysis. His time and effort is very much appreciated. Further, we thank J. Shelton (Bath) and J. King (Nottingham) for many stimulating and helpful discussions. We would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programmes ‘Complex Analysis: Techniques, Applications and Computations’ and ‘Applicable Resurgent Asymptotics: Towards a Universal Theory’, where some work on this paper was undertaken (EPSRC grant no. [EP/R014604/1]). Some work on this paper was also conducted while visiting the Okinawa Institute of Science and Technology (OIST) through the Theoretical Sciences Visiting Program (TSVP). C.A. and P.H.T. gratefully acknowledge support by the Engineering and Physical Sciences Research Council (EPSRC) [EP/V012479/1]. C.J.L. gratefully acknowledges support by Australian Research Council Discovery Project DP190101190.

Declaration of interests

The authors report no conflict of interest.

Appendix A. Further details on governing equations

In this section, we briefly present the details for deriving the governing equations. The fluid within the Hele-Shaw cell is governed by Stokes flow, and we consider the setup as shown in figure 3(a). We shall use bars for dimensional quantities, including Cartesian coordinates $(\bar {x},\bar {y})$ and polar coordinates $(\bar {r},\theta )$. Let ${\bar {\boldsymbol {q}}}$ be the dimensional velocity averaged over the gap. The fluid pressure, $\bar {p}$, is given by Darcy flow, ${\bar {\boldsymbol {q}}}=-(b^2/12\mu )\bar {\boldsymbol {\nabla }} \bar {p}$, where $b$ is the Hele-Shaw cell gap thickness and $\mu$ is the viscosity. We introduce a potential via $\bar {\phi }=-(b^2/12\mu )p$. Then, following the standard derivation via Stokes flow (Ockendon & Ockendon Reference Ockendon and Ockendon1995), we see the viscous fluid within the Hele-Shaw cell is governed by Laplace's equation, $\nabla ^2\bar \phi =0$, for a velocity potential with ${\bar {\boldsymbol {q}}}=\bar {\boldsymbol {\nabla }}\bar {\phi }$.

Inviscid fluid is injected at the origin and the interface between viscous and inviscid fluids is given by $\bar {r}=\bar {R}(\theta,\bar {t})$. The full set of governing equations, including interfacial and flux conditions are

(A1a)$$\begin{gather} \bar{\nabla}^2\bar{\phi}=0 \quad \text{in the fluid}, \end{gather}$$
(A1b)$$\begin{gather}\bar{v}_n=\frac{\partial\bar{\phi}}{\partial \bar{n}} \quad \text{on $\bar{r}=\bar{R}(\theta,\bar{t})$}, \end{gather}$$
(A1c)$$\begin{gather}\bar{\phi} - \bar{\phi}_0={-}\frac{b^2\sigma}{12\mu}\,\bar{\kappa} \quad \text{on $\bar{r}=\bar{R}(\theta,\bar{t})$}, \end{gather}$$
(A1d)$$\begin{gather}\frac{\partial\bar{\phi}}{\partial\bar{n}} =0 \quad \text{on $\bar{y}={\pm} \bar{x}\tan (\theta_0/2)$,} \end{gather}$$
(A1e)$$\begin{gather}\frac{\partial\bar{\phi}}{\partial\bar{r}} \sim \frac{\bar{Q}}{2{\rm \pi} b\bar{r}} \quad \text{as $\bar{r}\rightarrow\infty$} , \end{gather}$$

where $\sigma$ is the dimensional surface tension, $\bar {v}_n$ is the velocity normal to the interface, $\bar {\kappa }$ is the curvature of the interface and $\bar {\phi }_0$ an arbitrary shift of the potential. Equations (A1b) and (A1d) correspond to kinematic conditions on the interface and wall, respectively, (A1c) to the dynamic boundary condition on the interface and (A1e) the fluid injection condition for a dimensional source strength $\bar {Q}$.

A.1. Temporal scaling and non-dimensionalisation

We scale time by a yet-to-be-specified time scale, $T$, so that dimensionless time is $t=\bar {t}/T$ and introduce a length scale, $R_0$, which is the initial distance ($t = 0$) between the tip of the bubble and the apex. We shall apply a time-dependent stretching of the domain via

(A2)\begin{equation} \bar{R}(\theta,\bar{t})={A}(t)\bar{R}(\theta,0) = {A}(t)R_0\hat{R}(\theta), \end{equation}

where $A(t)$ is to be specified. Our dimensionless variables are now denoted without the overbar and we write

(A3ac)\begin{equation} (\hat{x},\hat{y})=\frac{(\bar{x},\bar{y})}{R_0{A}(t)}, \quad \hat{r}=\frac{\bar{r}}{R_0{A}(t)}, \quad \phi= \frac{T\bar{\phi}}{R_0^2}, \end{equation}

and so forth. The governing equations are now

(A4a)$$\begin{gather} \hat{\nabla}^2\phi =0 \quad \text{in the fluid}, \end{gather}$$
(A4b)$$\begin{gather}\hat{\boldsymbol{\nabla}}\phi\boldsymbol{\cdot} {\hat{\boldsymbol{n}}} = {A}(t){A}'(t){\hat{\boldsymbol{R}}}\boldsymbol{\cdot}{\hat{\boldsymbol{n}}}\quad \text{on $\hat{r} = \hat{R}(\theta)$}, \end{gather}$$
(A4c)$$\begin{gather}\phi - \phi_0={-}\frac{b^2T\sigma}{12\mu R_0^3{A}(t)}\,\hat{\kappa} \quad \text{on $\hat{r} = \hat{R}(\theta)$}, \end{gather}$$
(A4d)$$\begin{gather}\frac{\partial{\phi}}{\partial \hat{n}} =0 \quad \text{on $\hat{y}={\pm} \hat{x}\tan (\theta_0/2)$,} \end{gather}$$
(A4e)$$\begin{gather}\frac{\partial\phi}{\partial \hat{r}} \sim \frac{\bar{Q}\,T}{2{\rm \pi} b R_0^2 \hat{r}} \quad \text{as $\hat{r}\rightarrow\infty$.} \end{gather}$$

We review two ways to further reduce the above set of equations via the choice of $\textit {A}$. The first is to define $\hat {\phi }={A}(t)(\phi -\phi _0)$ and

(A5ac)\begin{equation} {A}(t)^2{A}'(t)=1, \quad \hat{\sigma}=\frac{b^2T\sigma}{12\mu R_0^3}, \quad \frac{\bar{Q}T{A}(t)}{b R_0^2}=1. \end{equation}

Then, the governing equations (A4) become

(A6a)$$\begin{gather} \hat{\nabla}^2\hat{\phi} =0 \quad \text{in the fluid}, \end{gather}$$
(A6b)$$\begin{gather}\hat{\boldsymbol{\nabla}}\hat{\phi}\boldsymbol{\cdot} {\hat{\boldsymbol{n}}} = {\hat{\boldsymbol{R}}}\boldsymbol{\cdot}{\hat{\boldsymbol{n}}} \quad \text{on the interface}, \end{gather}$$
(A6c)$$\begin{gather}\hat{\phi} ={-}\hat{\sigma}\hat{\kappa} \quad \text{on the interface}, \end{gather}$$
(A6d)$$\begin{gather}\frac{\partial \hat{\phi}}{\partial \hat{n}} =0 \quad \text{on $\hat{y}={\pm} \hat{x}\tan(\theta_0/2)$,} \end{gather}$$
(A6e)$$\begin{gather}\frac{\partial \hat{\phi}}{\partial \hat{r}} \sim \frac{1}{2{\rm \pi} \hat{r}} \quad \text{as $\hat{r}\rightarrow\infty$}. \end{gather}$$

In this way, our domain stretching function and time scale become

(A7a,b)\begin{equation} {A}(t)=(3t+1)^{1/3} \quad \text{and} \quad T=\frac{b R_0^2}{\bar{Q}(0)}, \end{equation}

our surface-tension parameter $\hat {\sigma }$ is a constant and the flow rate must be

(A8)\begin{equation} \bar{Q}(t)=\frac{\bar{Q}(0)}{(3t+1)^{1/3}}. \end{equation}

This way of reducing the equations has the attraction of keeping the surface tension constant (which is what happens in reality), while also requiring a non-constant flow rate with a $t^{1/3}$ scaling. Such a flow rate is perfectly realisable in practice (Li et al. Reference Li, Lowengrub, Fontana and Palffy-Muhoray2009).

The other way of reducing the equations is to set $\hat {\phi }=\phi -\phi _0$ and define

(A9ac)\begin{equation} {A}(t){A}'(t)=1, \quad \hat{\sigma}=\frac{b^2T\sigma}{12\mu R_0^3{A}(t)}, \quad \frac{\bar{Q}T}{b R_0^2}=1, \end{equation}

so the governing equations (A4) are the same as above (A6). This time, we have

(A10a,b)\begin{equation} {A}(t)=(2t+1)^{1/2} \quad \text{and} \quad T=\frac{b R_0^2}{\bar{Q}}. \end{equation}

This approach has the advantage of keeping the flow rate, $\bar {Q}$, constant; this is the most common realisation in an experimental set-up, but has the consequence of requiring the consideration of a time-dependent surface-tension parameter $\hat {\sigma }$ – conceptually, we may consider $\sigma$ as slowly varying, and even treated as constant within the context of a simulation (see further discussion in Ben Amar Reference Ben Amar1991b).

Appendix B. Deriving the power $\gamma$

In this section we derive the constants $\gamma _1$ and $\gamma _{C}$, which arise in (7.4), by matching the strengths of the singularities in the inner region. Firstly, we introduce a new variable, $Y = \zeta - \zeta _*$, where $\zeta _*$ is one of the singularities, $\zeta _1$ or $\zeta _{C}$.

It can be verified, either through dominant balance or via the numerical continuation of § 4, that all the complex plane singularities of the leading-order speed, $q_0$, are square-root singularities. Hence, $q_0 \sim Y^{-{1}/{2}}$ as $Y \rightarrow 0$. Close to the singularity, the Hilbert transform in the boundary-integral equation will be subdominant and so we can use $q_0 - \textrm {i}\tau _0 = \mathcal {H}[\tau _0]$ and the behaviour of $q_0$ to derive the local $\tau _0$ behaviour

(B1)\begin{equation} \tau_0 \sim \frac{\mathrm{i}}{2}\log Y \quad \text{as} \ Y \rightarrow 0, \end{equation}

and hence $\mathrm {e}^{\mathrm {i}\tau _0} \sim Y^{-1/2}$ in this limit. The asymptotic singular behaviour for other key variables is given as $Y \to 0$ by

(B2a)$$\begin{gather} x_0 + \mathrm{i} y_0 \sim \int_Y^{-\zeta_*} 1 \, {\rm d}{\tilde{Y}} = O(Y), \quad P = O(1), \end{gather}$$
(B2b)$$\begin{gather}I(\zeta) = O(1), \quad Q = O\left(Y^{{-}3/8}\right). \end{gather}$$

The local expression for the singulant equation (5.6) is given by

(B3)\begin{equation} \chi' = \frac{2\sqrt{\mathrm{i}}\exp({-\mathrm{i}\tau_0/2})}{q_0P}\exp({{3\theta_0}x_0/4}) = O\left(Y^{3/4}\right) \quad \mathrm{as} \ Y \to 0. \end{equation}

Integrating this expression gives $\chi = {O}(Y^{7/4})$ as $Y \to 0$. The factorial-over-power ansatz (5.3) then gives the strength of the singularity in the later-order terms

(B4)\begin{equation} q_n \sim \frac{Q\varGamma(2n+\gamma)}{\chi^{2n+\gamma}} = O\left(Y^{{-}3/8}Y^{-(2n+\gamma)7/4}\right) \quad \text{as} \ Y\rightarrow 0. \end{equation}

Now we take $n=0$ and choose $\gamma$ so that the correct predicted singularity strength is given for $q_0.$ Although the factorial-over-power ansatz was posed for the limit $n \rightarrow \infty,$ the singularity strength at lower orders still needs to match for the two expressions to be consistent. The $\gamma$ that satisfies this condition is

(B5)\begin{equation} \gamma = \tfrac{1}{14}. \end{equation}

This analysis is valid for each singularity, so $\gamma _1 = \gamma _{C} = 1/14$, irrespective of the wedge angle $\theta _0$.

Appendix C. Deriving the pre-factor $\varLambda$

To find $\varLambda$, which appears in the late-order expression (5.11) and also the exponential switching (7.4) we consider the solution in an inner region near the singularity at $\zeta = \zeta ^*$. We define a new inner variable $\nu$ such that $\zeta - \zeta _* = \epsilon ^\alpha \nu$. The value of $\alpha$ determines the width of the inner region. From (B4) with $\gamma = 1/14$, the local behaviour of the late terms near the singularity is given by

(C1)\begin{equation} q_n = O\left(\left(\epsilon^\alpha \nu\right)^{{-}3/8 - (2n+{1}/{14})7/4}\right) \quad \text{as} \ \epsilon \rightarrow 0. \end{equation}

We locate the inner region by identifying where terms in the power series (5.1) reorder, and the power series expansion is therefore no longer valid. This occurs when $\epsilon ^{2n}q_n$ is approximately ${O}(q_0)$ as $n \rightarrow \infty$. That is

(C2)\begin{equation} 2n + \alpha\left(-\frac{3}{8} - \frac{7n}{2} - \frac{1}{8}\right) ={-}\frac{\alpha}{2}, \end{equation}

which gives $\alpha = {4}/{7}$. The correct scaling for the inner region is thus

(C3)\begin{equation} \zeta = \zeta_* + \epsilon^{4/7} \nu, \end{equation}

which is identical to the scaling for the channel geometry (Chapman Reference Chapman1999).

We also require the local behaviour of the dependent variables near the singularities. This is done numerically using the analytic continuation scheme of § 4. Beginning from the free surface, we solve for the analytic continuation along a path that approaches the singularity, $\zeta _*$. Once done, the local behaviour of the dependent variables can be determined by numerical fitting to log–log plots. In particular, we find

(C4a)$$\begin{gather} x_0 \sim b_{1*} + b_{2*} \epsilon^{2/7}\nu^{1/2} + O\left(\epsilon^{4/7}\right), \end{gather}$$
(C4b)$$\begin{gather}q_0 \sim c_{1*}\epsilon^{{-}2/7}\nu^{{-}1/2} + c_{2*}\epsilon^{2/7}\nu^{1/2} + O\left(\epsilon^{6/7}\right), \end{gather}$$
(C4c)$$\begin{gather}\mathrm{e}^{i\tau_0} \sim d_{1*}\epsilon^{{-}2/7}\nu^{{-}1/2} + d_{2*} \epsilon^{2/7}\nu^{1/2} + O\left(\epsilon^{6/7}\right). \end{gather}$$

Here, the coefficients $b_{1*}, b_{2*}, c_{1*}, c_{2*}, d_{1*}, d_{2*}$ are constant with respect to the independent variables, but they depend on the parameters $\theta _0$ and $\lambda$ and also vary between the different singularities. We indicate the choice of singularity using the $*$ subscript. The constants are typically complex. Selected values of the constants as $\theta _0$ and $\lambda$ vary are tabulated in figure 12 and table 2 for the ‘1’ and ‘${C}$’ singularities, respectively.

Figure 12. Locations of the ‘1’ singularity for different values of $\lambda$ and $\theta _0$. At each singularity we have included the corresponding triplet of complex constants for the leading-order inner behaviour around the ‘1’ singularity, $(b_1, c_1, d_1).$

Table 2. The triplet of complex constants for the leading-order inner behaviours in the inner region about the ‘C’ singularity for different values of $\theta _0$ and $\lambda$.

From the expansions given above it follows that $\tau$ can be written in terms of $q$

(C5)\begin{align} \log q - \mathrm{i}\tau &= \hat{\mathcal{H}}[\tau_0] + O(\epsilon^2) = \log q_0 - \mathrm{i}\tau_0 + O(\epsilon^2) \nonumber\\ &= \log \left(\frac{c_{1*}}{d_{1*}}\right) + \left(\frac{c_{2*}}{c_{1*}} - \frac{d_{2*}}{d_{1*}}\right)\epsilon^{4/7}\nu + O\left(\epsilon^{8/7}\right), \end{align}

and $r$ can be expanded

(C6)\begin{equation} r(x) = r(x_0) + O(\epsilon^2) = r(b_{1*}) - r(b_{1*})\frac{\theta_0}{2}\epsilon^{2/7}\nu^{1/2} + O(\epsilon^{4/7}). \end{equation}

Next, we will derive the inner equation from the Bernoulli equation in terms of the inner variable for $q$, which is defined as

(C7)\begin{equation} q = c_{1*} \epsilon^{{-}2/7} \nu^{{-}1/2} q_{in}. \end{equation}

The analytically continued Bernoulli equation (4.3a) is

(C8) \begin{equation} \epsilon^2 \frac{P}{4}\frac{\partial}{\partial \zeta}\left(r(x)\left[Pq\frac{\partial \tau}{\partial \zeta} - \ell \sin\tau\right]\right) = \frac{2}{\rm \pi}\int_{-\infty}^0 \frac{K(\tilde{\zeta}) \tilde{\zeta}}{\zeta^2 - \tilde{\zeta}^2}\,{\rm d}{\tilde{\zeta}} + \mathrm{i} K(\zeta), \end{equation}

where $P,r(x), \ell, K(\zeta )$ are as defined in (2.7ac) and (2.9). The left-hand side and second term on the right-hand side can be straightforwardly expanded using the expressions above. The integral term on the right-hand side needs to be treated more carefully. Firstly, we expand using the asymptotic series

(C9)\begin{equation} \frac{2}{\rm \pi}\int_{-\infty}^0 \frac{K(\tilde{\zeta}) \tilde{\zeta}}{\zeta^2 - \tilde{\zeta}^2}\,{\rm d}\tilde{\zeta} = \frac{2}{\rm \pi}\int_{-\infty}^0 \frac{K_0(\tilde{\zeta}) \tilde{\zeta}}{\zeta^2 - \tilde{\zeta}^2}\,{\rm d}{\tilde{\zeta}}+ O(\epsilon^2), \end{equation}

where $K_0(\zeta )$ is the leading-order part of $K$ as defined in (3.2). Then we apply the residue theorem to the leading-order contribution by deforming the contour to a semicircular contour in the upper-half-$\zeta$-plane and noting the symmetry between the negative and positive real axes. At the first two orders the contribution from the integral cancels with the contribution from the second term on the right-hand side of Bernoulli's equation (4.3a). The leading-order inner equation therefore becomes

(C10)\begin{equation} -\mathrm{i}\frac{1}{4}P^2(\zeta_*)r(b_{1*})c_{1*}\frac{\mathrm{d}^2}{\mathrm{d}z^2}\left(\nu^{{-}1/2}q_{in}\right) = \frac{\nu}{2c_{1*}d_{1*}r^2(b_{1*})}\left(1-\frac{1}{q_{in}^2}\right). \end{equation}

It is natural to change variables

(C11)\begin{equation} \nu = A_*^{2/7} \nu_{in} = \left(-\frac{\mathrm{i} P^2(\zeta_*)r^3(b_{1*})c_{1*}^2d_{1*}}{2}\right)^{2/7} \nu_{in}, \end{equation}

so then the inner equation becomes

(C12)\begin{equation} - \frac{1}{\nu_{in}}\frac{\mathrm{d}^2}{\mathrm{d}\nu_{in}^2} \left(\nu_{in}^{{-}1/2}q_{in}\right) = \frac{1}{q_{in}} - 1, \end{equation}

which is the same as the inner equation for the classic Saffman–Taylor problem in a channel (Chapman Reference Chapman1999). In the outer limit (as $\tilde {\nu } \rightarrow \infty$) we thus derive the same relation as in the channel

(C13)\begin{equation} q_{in} = \sum_{n=0}^\infty A_n \nu_{in}^{{-}7n/2}, \end{equation}

where the constants $A_n$ satisfy the recurrence relation

(C14a)$$\begin{gather} \sum_{m=1}^n \left(\frac{7m}{2}-3\right)\left(\frac{7m}{2}-2\right)A_{m-1} \sum_{k=0}^{n-m}A_{n-m-k}A_k = \sum_{m=0}^n A_{n-m}A_m, \quad n\geqslant 1, \end{gather}$$
(C14b)$$\begin{gather}A_0 = 1. \end{gather}$$

Then the inner expansion gives

(C15)\begin{equation} q = c_{1*}\epsilon^{{-}2/7}\nu^{{-}1/2}\sum_{n=0}^\infty A_n \nu^{{-}7n/2}A_*^{n}, \end{equation}

which means the $n{\textrm {th}}$ term is

(C16)\begin{equation} \epsilon^{2n}q_n \sim c_{1*}\epsilon^{{-}2/7}\nu^{{-}1/2}A_n \nu^{{-}7n/2}\left(-\frac{\mathrm{i} P^2(\zeta_*)r^3(b_{1*})c_{1*}^2d_{1*}}{2}\right)^n, \quad \text{as} \quad n \rightarrow \infty. \end{equation}

Now, writing the outer expansion (5.7) in terms of the inner variable, we find

(C17)\begin{align} \chi &={-}\int_{\zeta_*}^\zeta \frac{2\mathrm{e}^{\mathrm{i}{\rm \pi}/4}\mathrm{e}^{-\mathrm{i}\tau_0/2}}{Pq_0r(x_0)^{3/2}}\,{\rm d}{\tilde{\zeta}} \nonumber\\ &\sim{-}\int_0^\nu \frac{2\mathrm{e}^{\mathrm{i}{\rm \pi}/4}d_{1*}^{{-}1/2}\epsilon^{1/7}\tilde{\nu}^{1/4}}{P(\zeta_*)r(b_{1*})^{3/2} c_{1*}\epsilon^{{-}2/7}\tilde{\nu}^{{-}1/2}}\epsilon^{4/7}\,{\rm d}{\tilde{\nu}} \nonumber\\ &= \frac{8}{7} \frac{\mathrm{e}^{{-}3\mathrm{i}{\rm \pi}/4}d_{1*}^{{-}1/2}}{c_{1*}P(\zeta_*)r(b_{1*})^{3/2}}\epsilon \nu^{7/4}, \end{align}

and so from the factorial over power ansatz (5.3)

(C18)\begin{align} \epsilon^{2n}q_n &\sim \varLambda_* \left(\frac{64}{49} \frac{\mathrm{i}}{d_{1*}c_{1*}^2 P^2(\zeta_*) r^3(b_{1*})}\right)^{{-}n} \nu^{{-}7n/2} \left(\frac{8\mathrm{e}^{{-}3\mathrm{i}{\rm \pi}/4}d_{1*}^{{-}1/2}\epsilon}{7c_{1*}P(\zeta_*)r(b_{1*})^{3/2}}\right)^{{-}1/14} \nonumber\\ &\quad \times \nu^{{-}1/8} \mathrm{e}^{7\mathrm{i}{\rm \pi}/8}2^{{-}1/2}d_{1*}^{{-}1/4}\epsilon^{{-}3/14} \nu^{{-}3/8} c_{1*}r(b_{1*})^{1/4}\varGamma(2n+\gamma), \end{align}

where the $*$ subscript on the $\varLambda$ labels the corresponding singularity. Then equating the inner and outer expansions gives

(C19)\begin{align} \varLambda_* = \frac{2^{5/7}}{7^{1/14}}\exp({-13\mathrm{i}{\rm \pi}/14})c_{1*}^{{-}1/14}d_{1*}^{3/14}P(\zeta_*)^{{-}1/14}r(b_{1*})^{{-}5/14} \lim_{n\rightarrow \infty} \frac{A_n}{\varGamma(2n+\gamma)}\left(\frac{32}{49}\right)^n. \end{align}

Using the large $n$ solution to the recurrence relation from Chapman (Reference Chapman1999) we obtain a final expression for the constant

(C20)\begin{equation} \varLambda_* \approx 0.46048 \times 2^{{-}1/7}\exp({-13\mathrm{i}{\rm \pi}/14})\left(\frac{d_{1*}^3}{c_{1*}P(\zeta_*)r(b_{1*})^5}\right)^{1/14} . \end{equation}

References

Abramowitz, M. & Stegun, I. 1972 Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover.Google Scholar
Ben Amar, M. 1991 a Exact self-similar shapes in viscous fingering. Phys. Rev. A 43 (10), 5724.CrossRefGoogle Scholar
Ben Amar, M. 1991 b Viscous fingering in a wedge. Phys. Rev. A 44 (6), 36733685.CrossRefGoogle Scholar
Booth, D.J., Griffiths, I.M. & Howell, P.D. 2023 Circular bubbles in a Hele-Shaw channel: a Hele-Shaw Newton's cradle. J. Fluid Mech. 954, 120.CrossRefGoogle Scholar
Brener, E.A., Kessler, D.A., Levine, H. & Rappel, W.-J. 1990 Selection of the viscous finger in the $90^\circ$ geometry. Europhys. Lett. 13 (2), 161166.CrossRefGoogle Scholar
Chandler, T.G.J. & Trinh, P.H. 2018 Complex singularities near the intersection of a free surface and wall. Part 1. Vertical jets and rising bubbles. J. Fluid Mech. 856, 323350.CrossRefGoogle Scholar
Chapman, S.J. 1999 On the role of Stokes lines in the selection of Saffman–Taylor fingers with small surface tension. Eur. J. Appl. Maths 10, 513534.CrossRefGoogle Scholar
Chapman, S.J., Dallaston, M.C., Kalliadasis, S., Trinh, P.H. & Witelski, T.P. 2023 The role of exponential asymptotics and complex singularities in self-similarity, transitions, and branch merging of nonlinear dynamics. Physica D 453, 114.CrossRefGoogle Scholar
Chapman, S.J., King, J.R. & Adams, K.L. 1998 Exponential asymptotics and Stokes lines in nonlinear ordinary differential equations. Proc. R. Soc. Lond. A 454, 27332755.CrossRefGoogle Scholar
Combescot, R. 1992 Saffman–Taylor fingers in the sector geometry. Phys. Rev. A 45 (2), 873884.CrossRefGoogle ScholarPubMed
Combescot, R., Dombre, T., Hakim, V., Pomeau, Y. & Pumir, A. 1986 Shape selection of Saffman–Taylor fingers. Phys. Rev. Lett. 56 (19), 20362039.CrossRefGoogle ScholarPubMed
Combescot, R., Hakim, V., Dombre, T., Pomeau, Y. & Pumir, A. 1988 Analytic theory of the Saffman–Taylor fingers. Phys. Rev. A 37 (4), 12701283.CrossRefGoogle ScholarPubMed
Crew, S. & Trinh, P.H. 2023 Resurgent aspects of applied exponential asymptotics. Stud. Appl. Maths 153 (3), 9741025.Google Scholar
Crew, S.C. & Trinh, P.H. 2016 New singularities for Stokes waves. J. Fluid Mech. 798, 256283.CrossRefGoogle Scholar
Dallaston, M.C. & McCue, S.W. 2014 Corner and finger formation in Hele-Shaw flow with kinetic undercooling regularisation. Eur. J. Appl. Maths 25 (6), 707727.CrossRefGoogle Scholar
Dingle, R.B. 1973 Asymptotic Expansions: their Derivation and Interpretation. Academic.Google Scholar
Franco-Gómez, A., Thompson, A.B. & Hazel, A.L. 2016 Sensitivity of Saffman–Taylor fingers to channel-depth perturbations. J. Fluid Mech. 794, 343368.CrossRefGoogle Scholar
Franco-Gómez, A., Thompson, A.B., Hazel, A.L. & Juel, A. 2017 Bubble propagation on a rail: a concept for sorting bubbles by size. Soft Matt. 13, 86848697.CrossRefGoogle Scholar
Gaillard, A., Keeler, J.S., Le Lay, G., Lemoult, G., Thompson, A.B., Hazel, A.L. & Juel, A. 2021 The life and fate of a bubble in a geometrically perturbed Hele-Shaw channel. J. Fluid Mech. 914, 140.CrossRefGoogle Scholar
Hong, D. & Langer, J. 1986 Analytic theory of the selection mechanism in the Saffman–Taylor problem. Phys. Rev. Lett. 56 (19), 20322035.CrossRefGoogle ScholarPubMed
Keeler, J.S., Gaillard, A., Lawless, J., Thompson, A.B., Juel, A. & Hazel, A.L. 2022 The interaction of multiple bubbles in a Hele-Shaw channel. J. Fluid Mech. 946, 126.CrossRefGoogle Scholar
Keeler, J.S., Thompson, A.B., Lemoult, G., Juel, A. & Hazel, A.L. 2019 The influence of invariant solutions on the transient behaviour of an air bubble in Hele-Shaw channel. Proc. R. Soc. Lond. A 475, 121.Google ScholarPubMed
Lajeunesse, E. & Couder, Y. 2000 On the tip-splitting instability of viscous fingers. J. Fluid Mech. 419, 125149.CrossRefGoogle Scholar
Lawless, J., Keeler, J.S., Hazel, A.L. & Juel, A. 2023 Stable bubble formations in a Hele-Shaw channel. J Fluid Mech. (submitted). https://arxiv.org/abs/2302.01038v3.Google Scholar
Levine, H. & Tu, Y. 1992 Mean-field diffusion-limited aggregation in radial geometries. Phys. Rev. A 45 (2), 10531059.CrossRefGoogle ScholarPubMed
Li, S., Lowengrub, J.S., Fontana, J. & Palffy-Muhoray, P. 2009 Control of viscous fingering patterns in a radial Hele-Shaw cell. Phys. Rev. Lett. 102, 14.CrossRefGoogle Scholar
McLean, J. & Saffman, P. 1981 The effect of surface tension on the shape of fingers in a Hele-Shaw cell. J. Fluid Mech. 102, 455469.CrossRefGoogle Scholar
Ockendon, H. & Ockendon, J.R. 1995 Viscous Flow, vol. 13. Cambridge University Press.CrossRefGoogle Scholar
Paterson, L. 1981 Radial fingering in a Hele-Shaw cell. J. Fluid Mech. 113, 513529.CrossRefGoogle Scholar
Pullin, D.I. & Meiron, D.I. 2013 Philip G. Saffman. Annu. Rev. Fluid Mech. 45, 1934.CrossRefGoogle Scholar
Saffman, P.G. & Taylor, G. 1958 The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Proc. R. Soc. Lond. A 245 (1242), 312329.Google Scholar
Shelton, J. & Trinh, P.H. 2022 Exponential asymptotics for steady parasitic capillary ripples on steep gravity waves. J. Fluid Mech. 939 (A17), 135.CrossRefGoogle Scholar
Tanveer, S. 1987 Analytic theory for the selection of a symmetric Saffman–Taylor finger in a Hele-Shaw cell. Phys. Fluids 30 (6), 15891605.CrossRefGoogle Scholar
Tanveer, S. 2000 Surprises in viscous fingering. J. Fluid Mech. 409, 273308.CrossRefGoogle Scholar
Thomé, H., Rabaud, M., Hakim, V. & Couder, Y. 1989 The Saffman–Taylor instability: from the linear to the circular geometry. Phys. Fluids A: Fluid 1 (2), 224240.CrossRefGoogle Scholar
Thompson, A.B. 2021 Bifurcations of drops and bubbles propagating in variable-depth Hele-Shaw channels. J. Engng Maths 129 (12), 117.CrossRefGoogle Scholar
Thompson, A.B., Juel, A. & Hazel, A.L. 2014 Multiple finger propagation modes in Hele-Shaw channels of variable depth. J. Fluid Mech. 746, 123164.CrossRefGoogle Scholar
Tu, Y. 1991 Saffman–Taylor problem in sector geometry: solution and selection. Phys. Rev. A 44 (2), 12031210.CrossRefGoogle ScholarPubMed
Vanden-Broeck, J.-M. 2010 Gravity-Capillary Free-Surface Flows. Cambridge University Press.CrossRefGoogle Scholar
Figure 0

Figure 1. (a) The bifurcation diagram for the classic Saffman–Taylor problem in a channel showing the first 10 branches of selected $\lambda$ values against the surface-tension parameter $\epsilon ^2$. Labels are shown for the first four branches: $\lambda _1$, $\lambda _2$, $\lambda _3$, $\lambda _4.$ This is equivalent to the results of Chapman (1999). (b) The solid lines show the bifurcation diagram which we calculate using exponential asymptotics for wedge angle $\theta _0 = 20^\circ$. This shows the permitted $\lambda (\epsilon )$ values that are selected by the selection mechanism. Details of the derivation of this selection mechanism follow in the rest of the paper. The circles show the numerically calculated values extracted from Ben Amar (1991b). Notice that $\lambda _1$ and $\lambda _2$ merge at $\epsilon \approx 0.3$ and $\lambda _3$ and $\lambda _4$ merge at $\epsilon \approx 0.07$.

Figure 1

Figure 2. Sketch of the free surface for viscous fingering in the full circular geometry, where a central source injects fluids outwards in all directions. Assuming axisymmetry, a single finger from this free surface can be considered as arising due to injecting fluid in the corner of the Hele-Shaw cell limited by sidewalls consisting of a wedge of angle $\theta _0$.

Figure 2

Figure 3. (a) A numerical plot of the top-down view of the self-similar physical profile in the $\hat {z}$-plane is shown for the zero-surface-tension case, with parameter values $\theta _0 = 20^{\circ }$ and $\lambda = 0.6$. The free surface was computed using (3.7). The Hele-Shaw cell is bounded by the thick black lines and is filled with a viscous fluid, shown in grey. An inviscid fluid is injected from the corner of the wedge and forms a finger with angle $\lambda \theta _0$. The corner of the wedge lies at $\hat z = 0$ ($BF$) and the tip of the finger lies at $\hat {z}=1$ ($CE$). (b) A sketch of the $z$-plane, computed via the conformal map $z=(2/\theta _0)\log \hat {z}$. This configuration is analogous to the traditional Saffman–Taylor finger (McLean & Saffman 1981).

Figure 3

Figure 4. The fluid region (grey) is mapped to the upper-half-$\zeta$-plane. The key points from figure 3 are labelled here. Within the $\zeta$-plane, there is a branch cut from the point $D$ ($\zeta = \mathrm {i}$). Here, the branch cut is taken vertically up the imaginary axis from $\zeta = \mathrm {i}$.

Figure 4

Figure 5. Plots of the leading-order (zero-surface-tension) solutions for the four variables $(x_0,y_0,q_0,\tau _0)$ on the free surface, $\zeta = \xi \in \mathbb {R}$. These plots show the solutions for parameter values $\theta _0 = 20^\circ$ and $\lambda = 0.6$ generated using (3.7) and (3.8a,b). The tip of the finger lies at $\xi = 0.$ In figure 3 we plot this solution in the physical plane.

Figure 5

Figure 6. Illustrations of the analytic continuation corresponding to $\theta _0 = 20^{\circ }$ and $\lambda = 0.6$ shown via (a) the $(\operatorname {Re} \zeta, \operatorname {Im} \zeta, \operatorname {Re} q_0)$-space; and (b) a top-down view of the $\zeta$-plane. In both, a prototypical path of analytic continuation from the physical free surface is shown with an arrow; branch cuts are shown wavy. Singularities (white circles) lie at $\zeta _1=0.59+1.32\mathrm {i}$ and $\zeta _{C}=0.96\mathrm {i}$; these both correspond to square-root singularities. There is an additional square-root branch point at $\mathrm {i}$. The leading-order solution $q_0$ on the free surface lies on the real $\zeta$ axis.

Figure 6

Figure 7. Plots of the locations of the three complex conjugate pairs of singularities in the $\zeta$-plane: (a) $\{\zeta _1, \overline {\zeta _1}\}$, (b) $\{\zeta _{C}, \overline {\zeta _{C}}\}$, (c) $\{\zeta _2, \overline {\zeta _2}\}$. These plots are for parameter values $\theta _0 = 20^\circ$ and $\lambda = 0.6.$ The branch cuts at $\pm \mathrm {i}$ are chosen to show the relevant branches of the Riemann surface which the singularities lie on. The free surface lies on the real $\zeta$ axis, $\zeta = \xi \in \mathbb {R}.$

Figure 7

Figure 8. Complex $\zeta$-plane showing how the location of the non-central singularities $\{\zeta _1, \overline {\zeta _1} \}, \{\zeta _2, \overline {\zeta _2} \}$ vary with the parameters $\theta _0$ and $\lambda$. The $\zeta _1$ singularity is shown in the top right quadrant with corresponding $\theta _0$ and $\lambda$ values. The locations of the $\zeta _2,$ $\overline {\zeta _2}$ and $\overline {\zeta _1}$ singularities are shown in the top left, bottom left and bottom right quadrants, respectively.

Figure 8

Table 1. The locations of the central singularity, $\zeta _{C}$ (to 4 significant figures) for different values of $\theta _0$ and $\lambda$.

Figure 9

Figure 9. Complex $\zeta$-plane showing the Stokes lines emanating from the singularities (circles) and intersecting the free surface (real axis). In this figure, the wedge angle is $\theta _0 = 20^\circ$ and $\lambda = 0.6$. Branch cuts (shown wavy) lie up and down the imaginary axis from $\pm \mathrm {i}$. Stokes lines are shown with dashed when they lie on a different Riemann sheet to the free surface. The three points where the Stokes lines intersect the real axis are labelled $S_1,$ $S_{C}$ and $S_2$.

Figure 10

Figure 10. Bifurcation diagrams produced using (7.4). The first ten selected values of $\lambda$ as functions of $\epsilon$ are shown. Panel (a) shows the bifurcation diagram for the channel geometry. The other panels (bd) show the bifurcation diagram for different values of the wedge angle $\theta _0$.

Figure 11

Figure 11. (a) Relative sizes of the real parts of $\chi _1(0)$ and $\chi _{C}(0)$ for different values of $\lambda$ in the $\theta _0 = 20^\circ$ case. The $\lambda$ value where $\mathrm {Re}(\chi _{C}(0)) =\mathrm {Re}(\chi _1(0))$ is denoted $\lambda _{lim}$. Above this value the selection condition (7.4) can be satisfied. (b) The limiting $\lambda$ value, $\lambda _{lim}$, in the small-surface-tension limit as a function of the wedge angle.

Figure 12

Figure 12. Locations of the ‘1’ singularity for different values of $\lambda$ and $\theta _0$. At each singularity we have included the corresponding triplet of complex constants for the leading-order inner behaviour around the ‘1’ singularity, $(b_1, c_1, d_1).$

Figure 13

Table 2. The triplet of complex constants for the leading-order inner behaviours in the inner region about the ‘C’ singularity for different values of $\theta _0$ and $\lambda$.