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

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.


Introduction
In their now classic work on viscous fingering, Saffman & Taylor (1958) consider the situation of a single finger of air-or an otherwise less viscous substancesteadily 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 λ ∈ (0, 1).Experimentally, Saffman and Taylor observed that λ ≈ 1/2.However, their asymptotic analysis, valid at small values of the non-dimensional surface tension parameter ϵ → 0, did not seem to restrict the value of λ.Today, it is known that for a fixed ϵ there exists a countably infinite number of possible values of λ = λ i with the property that 1/2 < λ 1 (ϵ) < λ 2 (ϵ) < • • • < 1. (1.1) In the limit ϵ → 0, every element in the family converges to 1/2 (see e.g.Vanden-Broeck 2010).The resultant plot of the eigenvalue, λ, versus the surface tension parameter, ϵ, is shown in fig.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-allorders asymptotics [see e.g.works by McLean & Saffman 1981;Hong & Langer 1986;Combescot et al. 1986Combescot et al. , 1988;;Tanveer 1987Tanveer , 2000;;Chapman 1999;Pullin & Meiron 2013].
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 wedgescenario 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 within the manner illustrated in fig. 2. Thus, the limited wedge configuration considered in this work serves as a model for each sector of the full source problem.
Consider the wedge problem characterised by a wedge angle, θ 0 , measured at the corner.The key parameter, λ, now describes the angular proportion, λθ 0 , of the cell occupied by the finger.In comparison with the previous fig.1(a) for classic Saffman-Taylor viscous fingering, in fig.1(b) we plot the bifurcation diagram for the case of a wedge of angle θ 0 = 20 • .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 (1991b).The figure shows the existence of distinct solution families, λ = λ i (ϵ, θ 0 ).
In addition to the expected solution families, λ i , there is now an additional phenomenon for the wedge-limited configurations.As seen in fig.1(b), there exist certain critical values of the surface tension parameter, ϵ, where each solution family, λ i , reaches a turning point connected to the adjacent family, i.e. λ i = λ 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 (1991b), the merging of eigenvalue pairs causes a loss of existence of the solutions for sufficiently small ϵ.Physically, in connection with the full geometry in fig.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 section 10.

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 (1981) and Thomé et al. (1989).Then in the early 1990s, a number of works appeared on the beyond-all-orders aspects of the wedge configuration (Brener et al. 1990, Ben Amar 1991a, Ben Amar 1991b, Tu 1991, Combescot 1992, Levine & Tu 1992).Combescot (1992) identified the singularities in the complex plane that are responsible for the selection mechanism.The solvability condition was found exactly by Brener et al. (1990) for the case with a 90 • separation angle, where the solution reduces to a closed analytic form.showing the first 10 branches of selected λ values against the surface tension parameter ϵ 2 .Labels are shown for the first four branches: λ1, λ2, λ3, λ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 θ0 = 20 • .This shows the permitted λ(ϵ) values that are selected by the selection mechanism.Details of the derivation of this selection mechanism follow in rest of the paper.The circles show the numerically calculated values extracted from Ben Amar (1991b).Notice that λ1 and λ2 merge at ϵ ≈ 0.3 and λ3 and λ4 merge at ϵ ≈ 0.07.

Source θ0
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 axi-symmetry, a single finger from this free surface can be considered as arising due to injecting fluid in the corner of Hele-Shaw cell limited by side walls consisting of a wedge of angle θ0.
Both Combescot (1992) and Tu (1991) used WKB methods to make analytic progress on the problem with a general wedge angle whilst numerical results were obtained by Ben Amar (1991b).This paper is most strongly motivated by the work of Tu (1991) and Ben Amar (1991b).Tu (1991) had linearised the free-surface problem with a general wedge angle to obtain a model differential equation.For this new equation, a WKBJ (Liouville-Green) approximation was used to derive a solvability condition that can predict the theoretical zero-surface tension limit for λ, as well as a condition that finds the point in the bifurcation diagram where branch merging occurs.
Ben Amar (1991b) 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 (1991b) 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 that 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 (Chapman et al. 1998) allow us to study the wedge problem in the small surface tension limit without the need to linearise in the same fashion as previous authors.In this paper, we will use these techniques to address the open problem identified by Ben Amar (1991b) and obtain an analytic solvability condition for the selected eigenvalues.
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 (Franco-Gómez et al. 2016;Thompson et al. 2014).Further, there have been many recent experimental (Lawless et al. 2023;Franco-Gómez et al. 2017;Gaillard et al. 2021), numerical (Thompson 2021;Franco-Gómez et al. 2017;Keeler et al. 2019) and analytic (Booth et al. 2023;Keeler et al. 2019) 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.

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, µ.For the case of a circular geometry, like that depicted in fig.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 fig.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 θ 0 .As can be observed experimentally (Thomé et al. 1989), a self-similar shape is reached eventually, where the inviscid fluid occupies an angle λθ 0 , with 0 < λ < 1.This set up is referred to as divergent flow in Ben Amar (1991a).For the case of zero surface tension, a prototypical solution is shown in fig.3(a), and we observe the petal-shaped interface between viscous fluid and inviscid fluid.
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 injection condition can be reposed in the self-similar framework.
The key idea relates to a transformation of the original dimensional lengths, x and ȳ, which are scaled via Figure 3: (a) A numerical plot of the top-down view of the self-similar physical profile in the ẑ-plane is shown for the zero surface tension case, with parameter values θ0 = 20 • and λ = 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 λθ0.The corner of the wedge lies at ẑ = 0 (BF ) and the tip of the finger lies at ẑ = 1 (CE).(b) A sketch of the z-plane, computed via the conformal map z = (2/θ0) log ẑ.This configuration is analogous to the traditional Saffman-Taylor finger (McLean & Saffman 1981).
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 result in effectively the same self-similar problem, given by (A 6), in dimensionless variables with a time-independent effective surface tension parameter σ.This σ, when rescaled, then gives us our small parameter ϵ, defined below in (2.10); see Ben Amar 1991b for further details.
We thus have the governing set of potential-flow equations in (A 6) corresponding to a scaled velocity potential, φ, and self-similar physical-plane coordinates, written in complex form as ẑ = x+iŷ [fig.3(a)].Note that the free-surface is now stationary.We introduce the complex potential, f = φ + i ψ with harmonic conjugate, the streamfunction ψ.Within the f -plane, fluid is located in the infinite strip bounded by ψ ∈ [−1, 1].
Finally, the flow domain is mapped to a channel geometry via so that the walls BA and F G lie on Im z = ±1, respectively, and the tip CE is fixed to the origin z = 0.This is shown in fig.3(b).Following section 3B of Ben Amar (1991b), we review the procedure for developing a set of boundary-integral equations for the potential-flow problem.First, the velocity potential is shifted as 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 ψ * = constant [cf.eqn (3.10) of Ben Amar (1991b)].For the case of the classic Saffman-Taylor viscous fingering problem in a parallel channel, H(z) = z, as shown by McLean & Saffman (1981).
As in Chapman (1999), it is convenient for later analysis to map the fluid region to the upper-half ζ-plane via (2.4) The mapped fluid domain for the configuration in fig. 3 is shown in fig. 4. Thus we see that under the map (2.4), the free surface, denoted BCEF , lies on the real ζ-axis, while the tip of the finger, denoted, CE, is at ζ = 0.
In the governing equations to follow, we shall seek to solve for the unknown free surface location and fluid velocities along the interface, ζ = ξ ∈ R. It is convenient to introduce quantities q and τ via which are analogues of speed and streamline angle, respectively (and reduce to the actual fluid speed and streamline angle in the limit θ 0 → 0).Therefore, we require a set of governing equations for the unknowns (x(ξ), y(ξ), q(ξ), τ (ξ)).
With the various conformal maps now established, at this point, we may follow the same procedures as found in section 3B of Ben Amar (1991b).We find on the free surface, where ζ = ξ is real, that continuity of pressure yields Bernoulli's equation: which can be compared to eqn (3.10) of Ben Amar (1991b).In (2.6), we have defined a number of functions for convenience.Firstly, we have written (2.7) Note that r is a function of x which, in turn, depends on ξ.In (2.6), we have defined Q 0 to be 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 θ 0 , where λ is the proportional finger angle parameter.Finally, we have introduced the key non-dimensional parameter, ϵ, by where σ is the modified surface tension parameter presented eqn (A 9) in Appendix A. We consider ϵ 2 to be a small parameter, corresponding to the small surfacetension regime, and we will therefore study the problem in the asymptotic limit ϵ → 0.
where we have defined the operator, H. Finally, we close the system by integrating the free surface velocity relationships (2.5).This yields (2.12) Thus, the full system consists of equations (2.6), (2.11) and (2.12) for the unknowns (x, y, q, τ ) in addition to the boundary conditions, (2.13b) Above, the first set of boundary conditions correspond to imposing the geometrical constraints while the second set correspond to the velocity and streamline angle constraints.In total, the equations and boundary conditions are equivalent to those of Ben Amar (1991a).In the next section, we shall examine the zero surface-tension solutions of these equations.

Analytic continuation of the free surface
The leading-order profiles, as evaluated on the physical free surface, ζ = ξ ∈ 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 section 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 , τ 0 ), as well as the analytic continuation of the governing equations (2.6), (2.11), (2.12).
4.1.Analytic continuation of the leading-order solutions In the analytic continuation, we allow the previously real-valued ξ to take complex values.Keeping in mind the potential for confusion with the previously introduced ζ, we write ξ = ξ r + iξ c → ζ ∈ C. Note that under this notational choice, q(ζ) is complex-valued within the fluid region, and it is rather the combination Re In theory, one can replace ξ with ζ, and evaluate the special-functions solution (3.7) using standard built-in packages (e.g.Mathematica) to obtain an analyticallycontinued 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 Ricatti equation (3.4) into its real and imaginary parts, where Recall that ℓ is the rescaled angle parameter introduced in (2.9).Here, we prefer to use ζ, as we are now working with the complexified version of the Ricatti equation (3.4).
We may now solve the above system along a chosen parameterised path in the complex ζ-plane by using the exact solution on the free surface as an initial condition.More specifically, we first pre-solve for (x 0 (ζ), y 0 (ζ), q 0 (ζ), τ 0 (ζ)) on the free surface using (3.7) and (3.8) and setting ζ = ξ ∈ R. Then we parameterise a path into the complex plane that starts on the free surface.For example, We can then solve (4.1) along the parameterised path using any standard ODE integrator (we use MATLAB's ode113 with absolute and relative tolerances set to 10 −10 ) and the initial condition ζ 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 , τ 0 )).In fig.6 we show the analytically continued surface for one of these components, 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.
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.
We can track the locations of the singularities as we vary the parameters θ 0 (the wedge angle) and λ (the proportion of the wedge angle occupied by the finger).If θ 0 > 0 and λ > 0.5, the non-central singularities {ζ 1 , ζ1 , ζ 2 , ζ2 } lie off the imaginary axis.In the limit θ 0 → 0 the two singularities ζ 1 and ζ 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 ζ1 and ζ2 in the lower-half ζ-plane.The singularity locations in the limit θ 0 → 0 agree with those found by Chapman (1999) in the classic Saffman-Taylor problem with parallel channel walls.A summary of the locations of these non-central singularities, for varying values of θ 0 and λ, is shown in fig.8.  Finally, the central singularity, ζ C , remains on the imaginary axis between the origin and i for all parameter values.It moves closer to i as either λ decreases towards 0.5, or as θ 0 decreases towards zero.A reflection of this behaviour occurs for ζC in the lower half-plane.Example locations of the singularity ζ C as θ 0 and λ are varied are listed in table 1.
To analytically continue the boundary integral equation (2.11) we must consider the complexification of the Hilbert transform, where Ĥ[τ ] is the complex-valued Hilbert transform, The boundary integral equation (2.11) then becomes and hence the integrated equation for the surface position (2.12) becomes, This results in a set of analytically continued governing equations (4.3) that hold in the upper-half ζ-plane.

Exponential asymptotics
Our procedure for the exponential asymptotic analysis follows similar work by Tanveer (1987) and Chapman (1999), using the methodology established in Chapman et al. (1998).In essence, the goal is to derive the behaviour of the late terms in the asymptotic series.After, in section 6, these late terms are used to study the exponentially-small terms via the Stokes-line switching.
As ϵ → 0, we expand the dependent variables as (5.1) We substitute the above into the analytically continued governing equations (4.3) and this yields, at O (ϵ 2n ) for Bernoulli's equation (4.3a), while for (4.3d) and (4.3e) we have (5.2b) (5.2c) At these later orders we see that the n th terms in the power series for q and τ 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 in the same locations as earlier terms, but with increasing strength.We therefore follow the method of Chapman et al. (1998) and predict a factorial-over-power form for the late-order terms, where Q and Θ are prefactors and χ is a singulant function, which is zero at the singularity.The singulant ensures that each series term has singularities with the correct locations and γ ensures they have the correct strengths.The Gamma function (Abramowitz & Stegun 1972) 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 1973 andCrew &Trinh (2023)].
In the limit n → ∞, the behaviour of the asymptotic expansion will be dictated by the divergence caused by the singularities driving (5.3) (Chapman et al. 1998).
From section 4 we know that the singularities lie in the complex plane away from the free surface.The complex Hilbert transform, Ĥ[τ n ], involves the integrand evaluated along the free surface.Once the ansatzes for q n and τ n via (5.3) are substituted into the Hilbert transform, we may observe that the contribution of Ĥ[τ n ] will be subdominant in the limit that n → ∞, compared to q n and τ n .This follows from the increasing nature of |χ| along Stokes lines, as explained on p. 526 of Chapman (1999).The combination of x n + iy 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 2022).
Using these assumptions gives the dominant behaviour from the boundary integral equation (5.2b) q n q 0 − q 1 q n−1 q 2 0 ∼ iτ n as n → ∞. (5.4) The Bernoulli equation (5.2a) can then be simplified to Here, and for the rest of the paper, we use primes ( ′ ) to denote differentiation with respect to ζ. Substitution of the factorial-over-power ansatz (5.3), and matching terms in the limit that n → ∞ gives at leading-order an equation for the singulant, χ: which can be solved to obtain (5.7) In this expression, ζ * is the singularity location and the integration limits are chosen so that χ(ζ * ) = 0.There will be a singulant function χ 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 ζ axis.Matching terms in (5.5) at the next order as n → ∞ shows that γ is constant, and at the following order we obtain, which can be solved to give an expression for the prefactor, where Λ is a constant that remains to be determined.From the boundary integral equation ( 5.4) we find, (5.10)By definition, ℓ → 0 (see (2.9)) and r(x 0 ) → 1 (see (2.7)) in the limit θ 0 → 0. This means that these results are consistent with those found by Chapman (1999) for the Saffman-Taylor problem in a channel geometry.The late-order series terms in the divergent power series solution for q(ζ) therefore have the form, (5.11)

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, 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 ϵ → 0 that R q q 0 ∼ iR τ .(6.2) Using this in the x + iy equation (5.2c) we see that R x and R y are subdominant compared to the leading-orders of R q and R τ in the limit ϵ → 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 (1999).This ordinary differential equation, known as the remainder equation, reduces to, as ϵ → 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 χ (primes will continue to denote differentiation with respect to ζ) gives, (6.4) Using the definition for the prefactor (5.9), this can be simplified to, To solve (6.5) we pose a Liouville-Green or WKBJ-style ansatz for the form of the remainder given by ϵ .(6.6) Then equating the coefficients at different powers of ϵ for the homogeneous version of the remainder equation ( 6.5), we find that b = ±χ + constant and B ∼ 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 ζ.This quantity is known as the Stokes multiplier, and we denote it by A(χ).The remainder becomes, R N = AQe − χ ϵ , (6.7) where the negative sign in the exponent ensures the remainder is exponentially small.The inhomogeneous equation (6.5) gives 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)This condition gives N ≈ |χ| 2ϵ .As N must be an integer we let N = |χ| 2ϵ + β, where β is bounded as ϵ → 0. This form motivates the transformation to polar coodinates, so we define χ = |χ|e iη .
By the chain rule we have dA dη = dA dχ iχ.(6.10) 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 1972) to approximate N ! in the large N limit.Then using (6.10) we obtain, We substitute in the optimal value of N to give, The right hand side of (6.12) is exponentially small in ϵ unless η = 2πk, k ∈ Z, in which case it will be algebraic in the limit that ϵ → 0. Therefore, the greatest change in the Stokes multiplier A will occur across curves, or Stokes lines, on which η = 2πk, and so Re(χ) > 0, Im(χ) = 0. (6.13)This recovers the classic result of Dingle (1973).We compute the Stokes lines numerically, with the results shown in fig.9.
When Stokes lines intersect the free surface (the real axis in the ζ-plane) then an exponentially small term will be smoothly switched on across this intersection point in the solution.We can see from fig. 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 η and A and consider the behaviour in the neighbourhood of the Stokes line.We apply the rescaling, η = ϵ 2 η and A = Ãϵ −γ−1 .(6.14)We let k = 0 and then (6.12) becomes (6.15)Integrating this gives (6.16) and hence the jump in the Stokes multiplier across the Stokes line is, (6.17) 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, iπ (6.18)where c.c. represents the complex conjugate expression, which arises from the singularity located at the complex conjugate location in the complex plane (fig.9).Similarly, by (5.10), the contribution that is switched on in the τ variable has the form, In the next section we show how these terms switched on across Stokes lines will determine the selection mechanism.

Selection mechanism
In fig. 9 we see that there are three points on the free surface (ζ = ξ ∈ R) that are intersected by Stokes lines.We will continue the notation from the singularities (introduced in section 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 Λ, γ and χ, 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 τ (−∞) = 0 and τ (∞) = −π imply there are no exponentially small contributions present on the free surface as τ → ±∞, 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 ζ = ξ ∈ (S 1 , S 2 ).
To check for existence of a solution we can check that the conditions in both farfields 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, τ (0) = −π/2.A similar symmetry condition is used in Chapman et al. (2023).Enforcing the condition on τ implies, Requiring that the solution satisfies the farfield conditions as ξ → −∞ 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 −∞ to η = 0.The symmetry condition from (7.1) says that the sum of these contributions must be zero at the origin.This means, where We simplify this expression using the conditions at the origin, τ 0 (0) = −π/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.13), i.e.Im(χ C (0)) = 0.Then, (7.3)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 τ variable to derive the selection mechanism.
The constants γ 1 and γ 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 γ 1 = γ C = 1/14 and thus we can simplify the selection condition to Above, the constants Λ 1 and Λ 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 Λ 1 and Λ C depend on both the parameters θ 0 and λ.For example, when θ 0 = 20 • and λ = 0.6 then Λ 1 ≈ −0.48 − 0.15i and Λ C ≈ −0.49− 0.24i.

Results
For each ϵ value we now use a numerical root finding scheme to find the corresponding family of λ values that satisfy the selection condition (7.4).The dependence on λ in the selection condition (7.4) appears through the constants Λ 1 and Λ C as shown in Appendix C and also through χ 1 (0) and χ C (0) as given in (5.7).The selected solution families, for different values of θ 0 , are shown in fig.10(a)-(d).
The figure shows that the bifurcation diagram differs qualitatively between the channel geometry with θ 0 = 0 • and the general wedge geometry with θ 0 > 0. For the channel geometry each branch, λ n (ϵ), continues to exist as ϵ → 0. However, for the wedge geometry we see that the branches merge and disappear in pairs and no solutions will exist when ϵ = 0.The lowest two branches, λ 1 (ϵ) and λ 2 (ϵ), merge at the greatest value of ϵ and then merging continues to happen between higher branches as ϵ → 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 fig. 10 we see that for smaller wedge angles the tip splitting will occur at a smaller value of the effective surface tension ϵ.In the limit θ 0 → 0, the surface tension value at which the branch merging occurs approaches zero.
As noted in e.g.Ben Amar (1991b), it is expected that, as is the case in the channel geometry, only the λ 1 branch is stable.Then the merging of the λ 1 and λ 2 branches as ϵ → 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 ϵ → 0, the classic Saffman-Taylor fingering in a channel involves λ = λ i → λ lim = 0.5 for all indices i.It can be asked whether a similar limiting value of λ is reached as ϵ → 0 and i → ∞ for general θ 0 .
In fig.11 we plot the real parts of χ 1 (0) and χ C (0) for the θ 0 = 20 • case.As indicated in the figure, for λ > λ lim , then 0 < Re χ 1 (0) < Re χ 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).
When Re(χ) is smaller then the corresponding e −Re(χ)/ϵ term will dominate in the selection condition and so fig.11 shows that for smaller λ values the contribution from the 'C' Stokes line dominates and for larger λ values the contribution from the '1' Stokes line dominates, in the small ϵ limit.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 ϵ limit it becomes impossible to satisfy the selection condition and no solutions will exist.The λ value at which Re(χ 1 (0)) = Re(χ C (0)) gives the limiting λ value, λ lim .The λ lim values as a function of θ 0 are shown in fig.11.

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 farfield 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 λ(ϵ) values and associated selected solutions that satisfy the selection condition.The bifurcation diagrams are plotted in fig. 10 and show the merging of branches of λ values as the surface tension parameter, ϵ, 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.

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 correctedness of the asymptotic predictions, we have used some of the early numerical work of Ben Amar (1991b).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.8)]; 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 ζ-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 (2018); Crew & Trinh (2016).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 (fig.10).For the Saffman-Taylor problem in the channel geometry, Tanveer (1987) showed that only the lowest branch, λ 1 , is linearly stable.We conjecture that the same will be true for the wedge geometry and so the merging of branches λ 1 and λ 2 will correspond to the loss of existence of any solutions which could be observed physically.We have written time-dependent code for the circular geometry based on the numerical method described in Dallaston & McCue (2014) 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 (Keeler et al. 2022;Gaillard et al. 2021).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.
This approach has the advantage of keeping the flow-rate, 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 σconceptually, we may consider σ as slowly varying, and even treated as constant within the context of a simulation (see further discussion in Ben Amar 1991b).

Appendix B. Deriving the power γ
In this section we derive the constants γ 1 and γ C , which arise in (7.4), by matching the strengths of the singularities in the inner region.It can be verified, either through dominant balance, or via the numerical continuation of section 4 that all the complex plane singularities of the leadingorder speed, q 0 , are square root singularities.Hence, q 0 ∼ Y − 1 2 as Y → 0. Close to the singularity, the Hilbert transform in the boundary integral equation will be subdominant and so we can use q 0 − iτ 0 = H[τ 0 ] and the behaviour of q 0 to derive the local τ 0 behaviour, and hence e iτ0 ∼ Y −1/2 in this limit.The asymptotic singular behaviour for other key variables is given as Y → 0 by, The local expression for the singulant equation (5.6) is given by, Integrating this expression gives χ = O(Y 7/4 ) as Y → 0. The factorial-over-power ansatz (5.3) then gives the strength of the singularity in the later order terms, Now we take n = 0 and choose γ so that the correct predicted singularity strength is given for q 0 .Although the factorial-over-power ansatz was posed for the limit n → ∞, the singularity strength at lower orders still needs to match for the two expressions to be consistent.The γ that satisfies this condition is, 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 ϵ 2n q n is approximately O(q 0 ) as n → ∞.That is, which gives α = 4/7.The correct scaling for the inner region is thus which is identical to the scaling for the channel geometry (Chapman 1999).
We also require the local behaviour of the dependent variables near the singularities.This is done numerically using the analytic continuation scheme of section 4. Beginning from the free surface, we solve for the analytic continuation along a path that approaches the singularity, ζ * .Once done, the local behaviour of the dependent variables can be determined by numerical fitting to log-log plots.In particular, we find: q 0 ∼ c 1 * ϵ −2/7 ν −1/2 + c 2 * ϵ 2/7 ν 1/2 + O ϵ 6/7 , (C 4b) e iτ0 ∼ d 1 * ϵ −2/7 ν −1/2 + d 2 * ϵ 2/7 ν 1/2 + O ϵ 6/7 .(C 4c) 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 θ 0 and λ 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 θ 0 and λ vary are tabulated in fig.12 and table 2 for the '1' and 'C' singularities respectively.where the constants A n satisfy the recurrence relation, Then the inner expansion gives, which means the n th term is,

Figure 1 :
Figure 1: (a) The bifurcation diagram for the classic Saffman-Taylor problem in a channelshowing the first 10 branches of selected λ values against the surface tension parameter ϵ 2 .Labels are shown for the first four branches: λ1, λ2, λ3, λ4.This is equivalent to the results ofChapman (1999).(b) The solid lines show the bifurcation diagram which we calculate using exponential asymptotics for wedge angle θ0 = 20 • .This shows the permitted λ(ϵ) values that are selected by the selection mechanism.Details of the derivation of this selection mechanism follow in rest of the paper.The circles show the numerically calculated values extracted from BenAmar (1991b).Notice that λ1 and λ2 merge at ϵ ≈ 0.3 and λ3 and λ4 merge at ϵ ≈ 0.07.

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

Figure 6 :
Figure 6: Illustrations of the analytic continuation corresponding to θ0 = 20 • and λ = 0.6 shown via (a) the (Re ζ, Im ζ, Re q0)-plane; and (b) a top-down view of the ζ-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 ζ1 = 0.59 + 1.32i and ζC = 0.96i; these both correspond to square root singularities.There is an additional square root branch point at i.The leading-order solution q0 on the free surface lies on the real ζ axis.

Figure 7 :
Figure 7: Plots of the locations of the three complex conjugate pairs of singularities in the ζ-plane: (a) {ζ1, ζ1}, (b) {ζC, ζC}, (c) {ζ2, ζ2}.These plots are for parameter values θ0 = 20 • and λ = 0.6.The branch cuts at ±i are chosen to show the relevant branches of the Riemann surface which the singularities lie on.The free surface lies on the real ζ axis, ζ = ξ ∈ R.

Figure 8 :
Figure 8: Complex ζ-plane showing how the location of the non-central singularities {ζ1, ζ1}, {ζ2, ζ2} vary with the parameters θ0 and λ.The ζ1 singularity is shown in the top right quadrant with corresponding θ0 and λ values.The locations of the ζ2, ζ2 and ζ1 singularities are shown in the top left, bottom left and bottom right quadrants respectively.

Figure 9 :
Figure 9: Complex ζ-plane showing the Stokes lines emanating from the singularities (circles) and intersecting the free surface (real axis).In this figure, the wedge angle is θ0 = 20 • and λ = 0.6.Branch cuts (shown wavy) lie up and down the imaginary axis from ±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 S1, SC and S2 respectively.

Figure 10 :Figure 11 :
Figure 10: Bifurcation diagrams produced using (7.4).The first ten selected values of λ as functions of ϵ are shown.Subfigure (a) shows the bifurcation diagram for the channel geometry.The other subfigures (b)-(d) show the bifurcation diagram for different values of the wedge angle θ0.
5)This analysis is valid for each singularity, so γ 1 = γ C = 1/14, irrespective of the wedge angle θ 0 .Appendix C. Deriving the pre-factor Λ To find Λ, 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 ζ = ζ * .We define a new inner variable ν such that ζ − ζ * = ϵ α ν.The value of α determines the width of the inner region.From (B 4) with γ = 1/14, the local behaviour of the late terms near the singularity is given by q n = O (ϵ α ν)

Figure 12 :
Figure12: Locations of the '1' singularity for different values of λ and θ0.At each singularity we have included the corresponding triplet of complex constants for the leading-order inner behaviour around the '1' singularity, (b1, c1, d1).

Table 1 :
The locations of the central singularity, ζC (to 4 s.f.) for different values of θ0 and λ. ).

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 θ0 and λ.