# Two-finger selection theory in the Saffman-Taylor problem

## Abstract

We find that solvability theory selects a set of stationary solutions of the Saffman-Taylor problem with coexistence of two unequal fingers advancing with the same velocity but with different relative widths $`\lambda _1`$ and $`\lambda _2`$ and different tip positions. For vanishingly small dimensionless surface tension $`d_0`$, an infinite discrete set of values of the total filling fraction $`\lambda =\lambda _1+\lambda _2`$ and of the relative individual finger width $`p=\lambda _1/\lambda `$ are selected out of a two-parameter continuous degeneracy. They scale as $`\lambda 1/2d_0^{2/3}`$ and $`|p1/2|d_0^{1/3}`$. The selected values of $`\lambda `$ differ from those of the single finger case. Explicit approximate expressions for both spectra are given.

PACS numbers: 47.54.+r, 47.20.Ma, 47.20.Ky, 47.20.Hw

The Saffman-Taylor problem has played a central role in the field of interfacial pattern selection in the last few decades. It deals with the morphological instability of the interface between two inmiscible fluids confined in a quasi twodimensional (Hele-Shaw) cell, when the less viscous fluid is displacing the more viscous one in a channel geometry. In particular, in their seminal work Saffman and Taylor called the attention upon the so-called selection problem, namely the fact that a unique finger-like steady state solution is observed whereas a continuum of solutions is possible if surface tension is neglected. Full analytical understanding of the subtle role of surface tension acting as the relevant selection mechanism was not achieved until much more recently. The resulting scenario of selection has been shown to apply with some genericity to other interfacial pattern forming systems, most remarkably in dendritic growth. On the other hand, despite the relative analytical tractability of the problem, the dynamics of competing fingers is far from being understood even at a qualitative level. Recently, it has been shown that in general surface tension may affect the long time dynamics in an essential way. In the case of the dynamics of finger arrays, the effect of surface tension becomes particularly dramatic, showing that the qualitative picture of finger competition based solely on the concept of screening of the laplacian field or the global instability of a periodic finger array is insufficient.

Existence of multifinger stationary solutions of the zero surface tension problem has been known for a long time. In Ref.  it has been emphasized that multifinger stationary solutions are relevant to the issue of the dynamical role of surface tension. In particular the equal-finger fixed point has been pointed out as the relevant saddle-point to describe competition dynamics. In connection with the phase flow structure around this fixed point, the problem of existence of unequal-finger fixed points with nonzero surface tension has been posed. Here we will extend selection theory to search for such solutions in the case of two fingers. We will follow the formulation of Hong and Langer, which is based on a Fredholm solvability analysis of a non-selfadjoint problem defined through linearization about the zero surface tension solution, together with WKB and steepest descent techniques. Despite the admitted objections to the full quantitative validity of the approximations involved in this method , it has been shown that it leads to the correct qualitative picture of selection and the correct scaling of solutions. It is therefore suitable, for simplicity of calculus and presentation, for a first exploration of novel situations such as the present one.

Our starting point is the dynamical equation for the conformal mapping $`f(w,t)`$ which maps the interior of the unit circle in the complex plane $`w`$ into the viscous fluid region, with the unit circle $`w=e^{i\varphi }`$ being mapped into the interface. Without loss of generality, we will assume a channel width $`W=2\pi `$ in the $`y`$ direction (with periodic boundary conditions) and a velocity $`U_{\mathrm{}}=1`$ of the fluid at infinity. We define the velocity of the stationary solutions of the interface as $`U=\frac{1}{\lambda }`$ where $`\lambda `$ is the total filling fraction of the channel by the invading fluid. The cartesian coordinates in the frame moving with velocity $`U`$ are given by $`z=x+iy=f(w,t)Ut`$. The mapping $`f(w,t)`$ contains a logarithmic singularity which is due to the fact that we are mapping an unbounded domain (the semiinfinite strip) into the unit circle, in such a way that $`f(w,t)+\mathrm{ln}w`$ is always an analytic funtion in the interior of the unit circle. The exact dynamical equation for the mapping can be written in the form

$$Re(i_\varphi f(\varphi ,t)_tf^{}(\varphi ,t))=1Ud_0_\varphi H_\varphi [\kappa ]$$
(1)   
which can be easily derived for instance from Ref.. Here $`d_0`$ is a dimensionless surface tension defined as $`d_0=\frac{\sigma b^2}{12\mu U}`$, where $`\sigma `$ is the surface tension, b is the gap and $`\mu `$ is the viscosity. The curvature $`\kappa `$ can be expressed in terms of $`x(\varphi )`$ and $`y(\varphi )`$ as

$$\kappa (\varphi )=\frac{_\varphi ^2x_\varphi y_\varphi ^2y_\varphi x}{[(_\varphi y)^2+(_\varphi x)^2]^{\frac{3}{2}}}.$$
(2)   
$`H`$ is a linear integral operator (Hilbert transform) which acts on a real $`2\pi `$-periodic function $`g(\varphi )`$ according to the definition

$$H_\varphi [g]=\frac{1}{2\pi }P_0^{2\pi }g(s)\mathrm{cotg}\frac{1}{2}(\varphi s)𝑑s.$$
(3)   
It follows from Eq.(3) that the function $`A(w)`$ defined by $`A(e^{i\varphi })=g(\varphi )+iH_\varphi [g]`$ is analytic in the interior of the unit circle, and $`Im(A(0))=0`$.

In the steady state we will have $`_tf^{}(\varphi ,t)=U`$, and Eq.(1) will read

$$U\frac{dy}{d\varphi }=1Ud_0\frac{dH_\varphi [\kappa ]}{d\varphi }$$
(4)   
We search for solutions of the generic type described in Fig.1. The total filling fraction $`\lambda `$ is splitted into two contributions $`\lambda _1+\lambda _2=\lambda `$, and we define as a new selection parameter the relative finger width $`p=\lambda _1/\lambda `$. For simplicity we will consider fingers which are axisymmetric and for convenience we will fix the tip positions at $`\varphi =\pi /2,3\pi /2`$, for all $`\lambda `$ and $`p`$. The filling fraction $`\lambda `$ ranges from $`0`$ to $`1`$. We take $`\lambda _2\lambda _1`$ so that $`p`$ ranges from $`1/2`$ to $`1`$. In these conditions the two fingers correspond to the intervals $`\varphi _1=\frac{\pi }{2}(12p)`$ to $`\varphi _2=\frac{\pi }{2}(1+2p)`$ and $`\varphi _2`$ to $`\varphi _3=2\pi +\varphi _1`$ respectively.

After integration over $`\varphi `$ of Eq.(4) we obtain

$$Uy(\varphi )=\varphi Ud_0H_\varphi [\kappa ]+c(\varphi )$$
(5)   
where $`c(\varphi )`$ is a piecewise constant function. The values it takes at the intervals $`(\varphi _1,\varphi _2)`$ and $`(\varphi _2,\varphi _3)`$ differ by the finite amount $`\pi (U1)`$ which accounts for the discontinuity of $`y`$ and the finite flux at the points $`\varphi _1,\varphi _2`$.

After Hilbert transform of Eq.(5) and using that $`H_\varphi ^2[g]=g(\varphi )`$, we obtain

$`d_0\kappa (\varphi )+x(\varphi )=x_0(\varphi )`$ (6)   
$`y(\varphi )+\varphi =const+H_\varphi [x].`$ (7)   
where Eq.(7) is just an expression of analyticity of $`f(w,t)+\mathrm{ln}w`$. The function $`x_0(\varphi )`$ is found explicitely as $`x_0(\varphi )=H_\varphi [g]`$ with $`g(\varphi )=(\lambda 1)\varphi +c(\varphi )`$, and by construction it corresponds to the solution of the zero surface tension case. In our case it reads

$`x_0(\varphi )=(1\lambda )\mathrm{ln}(2|\mathrm{sin}\varphi \mathrm{cos}p\pi |).`$ (8)   
Completed with $`y_0(\varphi )=\lambda (\varphi +c(\varphi ))`$, this gives a two parameter class of exact solutions of the type of Fig.1, for $`d_0=0`$. Both $`\lambda `$ and $`p`$ can be varied continuously within their natural range. The difference $`\mathrm{\Delta }_x`$ between the x-coordinate of the two tips is given by

$$\mathrm{\Delta }_x=(1\lambda )\mathrm{ln}\frac{1\mathrm{cos}p\pi }{1+\mathrm{cos}p\pi }.$$
(9)   
These solutions are precisely those studied in Ref..

The present formulation has some interesting advantages over the traditional approach of McLean and Saffman for instance in that the zeroth order solution is obtained naturally as an explicit outcome of the method and that it is more amenable to generalization, for instance to a larger number of fingers.

We now proceed by assuming $`x(\varphi )=x_0(\varphi )+d_0x_1(\varphi )`$ and linearizing on $`x_1(\varphi )`$ but keeping all singular terms necessary for selection. (Nonlinear effects are expected to introduce only a slight quantitative correction to the final spectrum of selection.) Using the relation $`y_1(\varphi )=H_\varphi [x_1]`$ we get

$`d_0{\displaystyle \frac{d^2x_1}{d\varphi ^2}}+d_0p(\varphi ){\displaystyle \frac{d^2H_\varphi [x_1]}{d\varphi ^2}}+r(\varphi )x_1=\mu (\varphi )`$ (10)   
where $`r(\varphi )`$ and $`p(\varphi )`$ are given by

$`r(\varphi )=\lambda ^2{\displaystyle \frac{|q(\varphi )|}{q^4(\varphi )}}[(q(\varphi ))^2+{\displaystyle \frac{1}{\beta ^2}}\mathrm{cos}^2\varphi ]^{\frac{3}{2}}`$ (11)   
$`p(\varphi )={\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\mathrm{cos}\varphi }{q(\varphi )}}`$ (12)   
with $`q(\varphi )=\mathrm{sin}\varphi \mathrm{cos}p\pi `$ and $`\beta =\frac{\lambda }{1\lambda }`$. Explicit knowledge of $`\mu (\varphi )`$ is not necessary for the solvability analysis. First order derivatives are subdominant as $`d_00`$ and have been omitted.

The linear operator on the lhs of Eq.(10) can be seen as a $`2\times 2`$ matrix operator acting on a vector of two components $`x_1^+(\varphi )`$ and $`x_1^{}(\varphi )`$ which are defined respectively on the intervals $`(\varphi _1,\varphi _2)`$ and $`(\varphi _1,\varphi _3)`$. Inserting an ansatz of WKB form with a point of stationary phase of the solution in the upper (or lower) complex plane one can show, using steepest descent techniques, that the off-diagonal terms of Eq.(10) lead to exponentially small contributions. As a consequence, to leading order the problem is decoupled into two separate problems defined in two disjoined intervals. Similarly, neglecting exponentially small terms, the integral part of the diagonal terms takes a purely differential form in the complex plane. The change of variables $`\eta =\frac{1}{\beta }\frac{\mathrm{cos}\varphi }{\mathrm{sin}\varphi \mathrm{cos}p\pi }`$ maps separately each of the two disjoint intervals above into the whole real axis $`\eta ϵ(\mathrm{},\mathrm{})`$. Therefore, to leading order we endup with two (complex) differential equations of the form

$`d_0{\displaystyle \frac{d^2x_1^\pm }{d\eta ^2}}+Q_\pm (\eta )x_1^\pm =R_\pm (\eta ).`$ (13)   
which are mutually independent but linked through the dependence on the parameters $`\lambda `$ and $`p`$. More details of this derivation will be presented elsewhere. We define two solvability functions as

$`\mathrm{\Lambda }_\pm (\lambda ,p;d_0)={\displaystyle _{\mathrm{}}^{\mathrm{}}}\stackrel{~}{x}^\pm (\eta )R_\pm ^{}(\eta )d\eta `$ (14)   
where $`\stackrel{~}{x}^\pm (\eta )`$ are eigenfunctions of the null space of the adjoint operators of the respective homogeneous equations. To enforce solvability we now have to impose the simultaneous vanishing of the two solvability functions $`\mathrm{\Lambda }_\pm (\lambda ,p;d_0)=0`$. These two conditions will fix the discrete spectra of possible values of both $`\lambda `$ and $`p`$.

Within the WKB approximation, the two solvability functions take the form

$`\mathrm{\Lambda }_\pm (\lambda ,p;d_0)={\displaystyle _{\mathrm{}}^{\mathrm{}}}G_\pm (\eta )e^{\frac{1}{\sqrt{d_0}}\mathrm{\Psi }_\pm (\eta )}d\eta `$ (15)   
where

$`\mathrm{\Psi }_\pm (\eta )=i\lambda \beta {\displaystyle _0^\eta }{\displaystyle \frac{(1i\eta ^{})^{\frac{1}{4}}(1+i\eta ^{})^{\frac{3}{4}}}{1+\beta ^2\eta ^2}}`$ (16)   
$`\times \left(1{\displaystyle \frac{\mathrm{cos}p\pi }{\sqrt{1+\eta ^2\beta ^2\mathrm{sin}^2p\pi }}}\right)\mathrm{d}\eta ^{}.`$ (17)   
In order to estimate the solvability functions in the steepest descent approximation, only the form of $`\mathrm{\Psi }_\pm (\eta )`$ is required. The singularity structure of $`\mathrm{\Psi }_\pm (\eta )`$ is such that the cases $`p=1/2`$ and $`p1/2`$ must be treated separately. The first case (two identical fingers) degenerates into the usual single finger problem. For $`p>\frac{1}{2}`$, a more complicated singularity structure is revealed. In the upper half complex plane of $`\eta `$, we find that $`\frac{d\mathrm{\Psi }_+(\eta )}{d\eta }`$ has a new branch point at $`\eta =i/\beta \mathrm{sin}p\pi `$, in addition to the singularities that were present in the single finger problem, namely, a branch point at $`\eta =i`$ and a pole at $`i/\beta `$. On the other hand, $`\frac{d\mathrm{\Psi }_{}(\eta )}{d\eta }`$ has the branch point at $`\eta =i`$ and the new one at $`\eta =i/\beta \mathrm{sin}p\pi `$, whereas the pole at $`i/\beta `$ is suppressed. Since $`1/\beta \mathrm{sin}p\pi >1/\beta `$, we obtain that $`\beta >1`$ is a necessary condition for the first solvability function $`\mathrm{\Lambda }_+(\lambda ,p;d_0)`$ to oscillate, and therefore generate zeroes. We thus recover the condition $`\lambda >1/2`$ of the single finger case, but now for the total filling fraction. The equivalent condition for $`\mathrm{\Lambda }_{}(\lambda ,p;d_0)`$ is $`\beta \mathrm{sin}p\pi >1`$ so that the new singularity at $`\eta =i/\beta \mathrm{sin}p\pi `$ stands below $`\eta =i`$. This condition also implies that in the contour integration for $`\mathrm{\Lambda }_+(\lambda ,p;d_0)`$ we will always pick up a contribution from this new singularity.

By deforming the contour integral as indicated in Fig.2, and following Ref.  in identifying the crossover from oscillating to non oscillating behaviour of the solvability functions, we obtain the scaling of both $`\lambda `$ and $`p`$ with $`d_0`$ to be $`(\lambda \frac{1}{2})d_0^{2/3}`$ and $`|p\frac{1}{2}|d_0^{1/3}`$. According to Eq.(9), the resulting scaling for the tip difference is $`\mathrm{\Delta }_xd_0^{1/3}`$.

An explicit (approximate) discrete spectra of selected values of $`\lambda `$ and $`p`$ for small $`d_0`$ will be given by the condition $`\mathrm{cos}(\frac{\mathrm{\Psi }_\pm (i+0)\mathrm{\Psi }_\pm (i0)}{2i\sqrt{d_0}})=0`$.

From the condition $`\mathrm{\Lambda }_{}(\lambda ,p;d_0)=0`$ we thus obtain

$`{\displaystyle \frac{1}{\sqrt{d_0}}}{\displaystyle \frac{(1\lambda )^2}{\lambda }}I(\beta ,p)=m{\displaystyle \frac{1}{2}}`$ (18)   
with $`m=1,2,\mathrm{}`$ and where

$`I(\beta ,p)={\displaystyle \frac{1}{2\pi }}\mathrm{cotg}p\pi {\displaystyle _0^{u_1}}{\displaystyle \frac{u^{\frac{3}{4}}H(u;\beta ,p)}{(u_3u)(u_1u)^{\frac{1}{2}}}}du`$ (19)   
with the regular part of the integrand $`H(u;\beta ,p)=(2u)^{\frac{1}{4}}(u_4u)^1(u_2u)^{\frac{1}{2}}`$ and $`u_{1,2}=11/\beta \mathrm{sin}p\pi `$, $`u_{3,4}=11/\beta `$.

Finally, from condition $`\mathrm{\Lambda }_+(\lambda ,p;d_0)=0`$, expressing Eq.(18) to leading order and using properties of hypergeometric functions, the two selection conditions can be combined to read

$`{\displaystyle \frac{1}{\sqrt{d_0}}}(2\lambda 1)^{\frac{3}{4}}=n`$ (20)   
$`S(\alpha )={\displaystyle \frac{1}{n}}(m{\displaystyle \frac{1}{2}})`$ (21)   
with $`n=1,2,\mathrm{}`$ and where

$$S(\alpha )=\frac{3\sqrt{2\pi }}{5\mathrm{\Gamma }^2(\frac{1}{4})}_2F_1(\frac{5}{4},\frac{1}{2};\frac{9}{4};1\alpha )(1\alpha )^{\frac{5}{4}}.$$
(22)   
$`{}_{2}{}^{}F_{1}^{}`$ is a hypergeometric function  and $`\alpha =\frac{\pi ^2}{4}\frac{(p\frac{1}{2})^2}{2\lambda 1}`$ is of order $`(d_0)^0`$ and ranges from $`0`$ to $`1`$.

Eq.(20) determines a set of discrete values of $`\lambda `$. Notice that these are given independently of $`p`$ but the set of values are inserted between those of the single finger case ($`p=1/2`$), which in the same approximation are given by $`\frac{1}{\sqrt{d_0}}(2\lambda 1)^{\frac{3}{4}}=n\frac{1}{2}`$ in place of Eq.(20). On the other hand, the lhs of Eq.(21) is a monotonically decreasing function of $`\alpha `$ which varies continuously from $`1/4`$ (at $`\alpha =0`$) to $`0`$ (at $`\alpha =1`$). Solving Eq.(21) for $`\alpha `$ produces solutions with $`p1/2`$. These will exist whenever $`\frac{1}{n}(m\frac{1}{2})<\frac{1}{4}`$. For a given $`n`$, the solutions are labeled by $`m=1,2,\mathrm{}`$ up to the integer part of $`(n+1)/4`$. Therefore, the first solution with $`p1/2`$ will appear at $`n=3`$ and gives $`|p\frac{1}{2}|0.3886d_0^{1/3}+\mathrm{}`$ For fixed $`m`$, $`p`$ is an increasing function of $`n`$ (like $`\lambda `$), but for fixed $`n`$, $`p`$ has its maximum value at $`m=1`$ and then decreases with $`m`$ (See Fig.3). The spectra here derived must be taken with some caution, since they are only approximate. An exact calculation to lowest order in $`d_0`$ should include nonlinear effects and a proper treatment of the turning points in the WKB analysis, but the corrections are expected to be quantitatively small. More details will be presented elsewhere.

Concerning the stability of these solutions, it is reasonable to presume that, in general, they will be globally unstable, such as established numerically for the equal finger array ($`p=1/2`$) in Ref.. This implies that they would only be directly observable as a transient slowing down of the competition dynamics whenever an initial condition is prepared close to any of those solutions. From a dynamical systems point of view, we would like to emphasize that the knowledge of the fixed points, even if unstable, is always relevant to elucidate the topological features of the phase space flow, and therefore to gain insight and qualitative understanding of the dynamics. In particular the location of these new fixed points will definitely affect the path in phase space describing the transient dynamics from a nearly equal finger array towards the single finger attractor. This point of view was developed in Ref. to study the dynamical role of surface tension. In that spirit it was pointed out that a generalized solvability scenario of selection could hold to some extent for the dynamics.

We conclude by remarking that, although the present solvability analysis is not a rigorous proof of existence of solutions it reveals by itself a quite unexpected richness of the problem. It would be interesting to search for these solutions numerically or by other more rigorous means . The sole existence of the predicted solutions and its presumable generalization to a larger number of fingers has important consequences on the physical picture of finger competition, which turns out to be much more complex than common arguments of laplacian screening seem to suggest. The common picture, according to which fingers slightly ahead escape from their neighbors, is not necessarily valid in general because of the existence of growth modes with unequal noncompetiting fingers. For vanishingly small surface tension, however, these modes collapse and only the equal-finger multifinger mode ($`p=1/2`$) survives as a stationary state. Finally, given the genericity of the solvability mechanism of selection, this opens the possibility of finding similar solutions in related problems such as needle crystal growth or viscous fingering in circular geometry.

We acknowledge financial support from the Dirección General de Enseñanza Superior (Spain) under Project PB96-1001-C02-02. F.X.M. also acknowledges financial support from the Comissionat per a Universitats i Recerca (Generalitat de Catalunya).

Fig. 1. Typical configuration of a two-finger stationary solution.

Fig. 2. Deformation of the steepest descent contour of integration in the complex $`\eta `$ plane with $`\beta >1`$ and $`p>\frac{1}{2}`$ (a) for $`\mathrm{\Lambda }_+`$; (b) for $`\mathrm{\Lambda }_{}`$.

Fig. 3. Spectrum of $`p`$ as a function of $`n`$ for different values of $`m`$.