Effects of small surface tension in Hele-Shaw multifinger dynamics: an analytical and numerical study

We study the singular effects of vanishingly small surface tension on the dynamics of finger competition in the Saffman-Taylor problem, using the asymptotic techniques described in [S. Tanveer, Phil. Trans. R. Soc. Lond. A 343, 155 (1993)]and [M. Siegel, and S. Tanveer, Phys. Rev. Lett. 76, 419 (1996)] as well as direct numerical computation, following the numerical scheme of [T. Hou, J. Lowengrub, and M. Shelley,J. Comp. Phys. 114, 312 (1994)]. We demonstrate the dramatic effects of small surface tension on the late time evolution of two-finger configurations with respect to exact (non-singular) zero surface tension solutions. The effect is present even when the relevant zero surface tension solution has asymptotic behavior consistent with selection theory.Such singular effects therefore cannot be traced back to steady state selection theory, and imply a drastic global change in the structure of phase-space flow. They can be interpreted in the framework of a recently introduced dynamical solvability scenario according to which surface tension unfolds the structually unstable flow, restoring the hyperbolicity of multifinger fixed points.


I. INTRODUCTION
The displacement of a viscous fluid by a less-viscous one in a Hele-Shaw cell, the so-called Saffman-Taylor problem [1,2,3,4,5], is a prototypical pattern formation problem. Since the seminal work of Saffman and Taylor [1] a considerable effort has been aimed at understanding both steady and unsteady interfacial patterns formed during this flow. The Saffman-Taylor problem is the simplest member of a wide class of interfacial pattern formation problems such as free dendritic growth, directional solidification, or chemical electro-deposition [6,7,8]. As such, a theoretical understanding of Hele-Shaw flow may help elucidate generic behavior common to many pattern forming systems. Despite its relatively simple formulation and the large amount of work devoted to it, however, several aspects of interfacial dynamics in Hele-Shaw flow are still poorly understood, in particular concerning the highly nonlinear and nonlocal dynamics of finger competition.
One of the main reasons for the wide interest in Hele-Shaw flow, at least from a mathematical point of view, is that explicit time-dependent solutions can be found in the case of zero surface tension [9,10,11,12]. However, it is also known that the zero surface tension Saffman-Taylor (ST) problem is ill-posed as an initial value problem [13] and finite-time singularities appear frequently [14]. Nevertheless, rather large classes of zero surface tension solutions have been found which exhibit the variety of morphologies observed both in experiments and numerical simulations. Then, the question that naturally arises is to what extent smooth (nonsingular) zero surface tension solutions reproduce the dynamics of the physical problem with finite surface tension, in particular in the limit of vanishing dimensionless surface tension, B → 0.
It is well known that surface tension is a singular per-turbation to the zero surface tension problem [13]. This singular character manifests dramatically in the classical selection problem posed by Saffman and Taylor [1] and only solved three decades later [15,16,17,18], where an arbitrarily small surface tension selects out a single, stable solution from the continuum of steady single-finger B = 0 solutions. Another manifestation of the singular nature of surface tension which is directly relevant to the present work is its effect on the dynamics. Siegel, Tanveer and Dai [19,20] showed that interfacial evolution for the regularized problem (i.e., vanishingly small B) may differ signifficantly from that for the B = 0 problem in order one time. Then, smooth time dependent solutions of the B = 0 case do not coincide, in general, with the limiting solutions for B → 0. Accordingly, the study of finite B dynamics encounters considerable difficulties since B = 0 solutions cannot be naively used as a starting point for the study of the problem with finite B.
The physical content of exact zero surface tension solutions with pole-like singularities has been recently addressed in Refs. [5,21,22] using a dynamical systems approach. Through a detailed study it has been shown that the exact zero-surface tension phase flow, considered in a global sense, is structurally unstable. In other words, the zero surface tension phase dynamics are not topologically equivalent to the phase space flow of the physical problem, regularized by surface tension. Indeed, the zero surface tension phase flow omits the necessary saddle-point structure of multifinger fixed points, which is crucial to the physical finger competition process [22]. A natural extension of the well known solvability mechanism (first applied to 'select' a finger of width 1/2 out of a continuum of solutions in the single finger case) was proposed for multi-finger solutions in [22]; this helps clarify how the introduction of surface tension modifies the global phase space structure of the flow and restores the hyperbolicity of multi-finger fixed points.
The approach of Ref. [22], however, was qualitative in nature and could not quantify the extent to which zero surface tension trajectories might resemble the evolution with small surface tension. In particular it was recognized that, while some trajectories appear to be qualitatively correct for infinite time, others may have a dramatically different evolution. In particular, adding an infinitessimal surface tension could give the opposite outcome in the finger competition, that is, make the 'losing' finger for B = 0 become the 'winning' finger when B > 0, for sufficiently generic sets of initial conditions.
A satisfactory analytical understanding of the problem with B = 0 has been achieved in two regimes: the initial linear instability of the flat interface followed by the weakly non-linear regime [23], and the asymptotic regime where surface tension selects the width of the single finger [15,16,17,18]. The highly non-linear intermediate regime that connects the quasi-planar interface with the asymptotic single-finger regime has mostly been studied through numerical computation. The first exhaustive numerical studies were reported by Tryggvason and Aref [24,25], who paid considerable attention to the influence of viscosity contrast on the problem and studied both single-finger and multi-finger configurations. Later, Casademunt and Jasnow [26,27] showed that the basin of attraction of the single-finger solution depends strongly on viscosity contrast and that only when one of the two fluid viscosities is negligible it can be claimed that the single finger is the universal attractor of the problem. In the present work we will restrict to this limiting case. DeGregoria and Schwartz [28,29] observed that welldeveloped fingers split when surface tension is sufficiently decreased. This tip-splitting instability is related to the fact that the Saffman-Taylor finger is linearly stable but non-linearly unstable, and the size of the perturbation that triggers the tip-splitting decreases quickly with surface tension [30]. Dai and Shelley [31] showed that for small B numerical computations are extremely sensitive to the precision used in the computations. As a consequence noise level has to be controlled with care in order to ensure that the computation is sufficiently accurate. Hou et al. [32] developed a numerical method that deals with the numerical stiffness of the problem in an efficient manner, aiding the ability to perform long time computations. More recently Ceniceros et al. [33,34], using very high precision arithmetic have been able to study the effect of extremely small surface tension in the circular geometry with suction, and they have observed that surface tension can produce complex ramified patters even without the presence of noise.
An analytical treatment of this highly nonlinear and nonlocal free-boundary problem faces challenging difficulties. In particular, a perturbative study for small B is complicated by the ill-posedness of the zero surface tension problem. Tanveer [13] was able to overcome this obstacle by embedding the zero surface tension problem in a well-posed one. In addition, this well-posed extension of the B = 0 problem allowed Baker et al. [35] to develop a numerical method to compute the time evolution of zero surface tension dynamics in a well-posed manner. Once the B = 0 problem is formulated in a well-posed way the B = 0 case can be studied using a perturbative approach. The main result of the asymptotic perturbative theory developed by Tanveer [13] is that the effect of surface tension may be manifest in a O(1) time: the evolution of the same initial interface for B = 0 and B = 0 will in general differ after a time of order one, even if the B = 0 solution is smooth for all time. Siegel et al. [20] have extended the work of Ref. [13] to later stages of the evolution, and through numerical computation for very small values of B they showed that smooth B = 0 solutions are indeed significantly affected by the presence of arbitrarily small B in order-one time, thus confirming the predictions of the perturbative theory. The zero surface tension solutions studied by Siegel et al. [19,20] in the channel geometry were single-finger solutions with an asymptotic width λ, specifically chosen to be incompatible with selection theory for vanishing surface tension. They found that the singular effect of surface tension was to widen the finger in order to reach the selected width. The surprising feature here is that the effect of surface tension is felt in order-one time, i.e., that the time lapse for which the regularized solution approaches the unperturbed one as B → 0 is bounded.
The present paper expands the work of Refs. [19,20] in the spirit of Ref. [22], towards the study of multi-finger solutions. However, unlike the studies of [19,20] we chose zero surface tension multi-finger soltuions which are compatible with asymptotic selection theory, that is, with an asymptotic finger width λ = 1/2-the selected value in the limit B → 0. In this way we isolate the intrinsic finger competition dynamics from the selection effects responsible of restoring the asymptotic width. Two different kinds of two-finger zero surface tension solutions are studied, and in both cases it is shown that surface tension acts as a singular perturbation to the dynamics in orderone time, modifying dramatically the late time configuration of the interface not only quantitatively but also qualitatively. Specifically we show that paths in phase space associated with zero and nonzero surface tension evolution, and indeed the global topological structure of the the phase spaces, may differ appreciably, even for arbitrarily small B. In physical terms, our evidence suggests that the presence of arbitrarily small surface tension can completely alter outcome of finger competition when compared with the zero surface tension evolution.
The paper is organized as follows: in Sec. II the equations describing Hele-Shaw flow are introduced, and a class of two-finger zero surface tension solutions relevant to two-finger competition is presented and briefly discussed. This class of solutions will be used as initial condition for numerical computation with B > 0. In Sec. III the basic features of the asymptotic theory are recalled, and the theory is applied to the zero surface tension solutions introduced in the previous section. The numerical computations with finite (but small) B are presented in Sec. IV. Sec. V discusses and summarizes the results obtained in previous sections.

II. ZERO SURFACE TENSION
In this section we present the equations which govern the interfacial dynamics in a rectilinear Hele-Shaw cell, following the formalism of [13]. We consider a class of exact, time-dependent zero surface tension solutions that are relevant to the finger competition problem, and briefly describe the solutions within this class.
Consider Hele-Shaw flow in the channel geometry, in which a fluid of negligible viscosity displaces a viscous liquid. The equations governing the interfacial evolution can be conveniently formulated by first introducing a conformal map z(ζ, t) which takes the interior of the unit semicircle in the ζ plane into the region occupied by the viscous fluid in the complex plane z = x + iy, in such a way that the arc ζ = e is for s ∈ [0, π] is mapped to the interface and the diameter of the semi-circle is mapped to the channel walls [38]. The mapping function z(ζ, t) has the form z(ζ, t) = −(2/π) ln ζ + i + f (ζ, t), and inside and on the unit semi-circle we require f (ζ, t) to be analytic and z ζ (ζ, t) = 0. In addition, we require that on the real diameter of the semi-circle. This latter condition ensures that z maps the diameter to the channel walls. Under suitable assumptions (see [13]) the Schwartz reflection principle may be applied to show that f is analytic and z ζ = 0 for |ζ| ≤ 1. The effective velocity field, averaged across the plate gap, is a two-dimensional potential flow satisfying Darcy's law Here ϕ is a velocity potential defined by ϕ = −(b 2 /12µ) p, where p is the pressure, µ is the viscosity and b is the gap width. Under the assumption of incompressibility (∇ · u = 0) the potential satisfies Laplace's equation Incompressibility also implies the existence of a stream function ψ. Therefore, one can define a complex velocity potential W (z, t) = ϕ + iψ which is analytic for z in the fluid region of the channel. Its form as a function of ζ reads where ω(ζ, t) is an analytic function inside the unit circle. The condition that there is no flow through the walls implies that Im ω = 0 must hold on the real diameter of the unit semi-circle. In the absence of surface tension, ω = 0 (see Eq. (6)).
At the interface we impose the usual boundary conditions. The kinematic boundary condition states that the normal component of fluid velocity at a point on the interface equals the normal velocity of the interface at that point, and takes the form The dynamic boundary condition specifies that the pressure jump across the interface is balanced by surface tension, and is given by The parameter B is the nondimensional surface tension and is defined by where T is the surface tension, V is the fluid velocity at infinity and a is half the cell width. The equations given in (1)-(5) are in nondimensional form, with lengths and velocities nondimensionalized by a and V , respectively. When B = 0 it is well known that pole singularities in z ζ (i.e. in f ζ ) present in the exterior of the unit disk are preserved under the dynamics, i.e., such singularities are neither created nor destroyed, although the location of those which are initially present will evolve with time. Exact B = 0 solutions consisting of a collection of pole singularities with constant amplitude have been the focus of extensive studies (see e.g. [9] ). The simplest such solution leading to nontrivial finger competition consists of a pair of singularites in the upper halfplane of |ζ| > 1, located at positions that are symmetric with respect to the y-axis. A second pair of poles conjugate to the first pair is required to satisfy the symmetry restriction (1). This exact solution takes the form [5,21,22] where λ and ǫ are real constants with 0 < λ < 1 and ǫ ≥ 0, and d(t) is real. The singularity locations are given by the complex parameter ζ s (t), which satisfies a simple differential equation given in [21]. Analyticity of f (ζ, t) in the unit circle implies that |ζ s (t)| > 1. We employ the convention that ζ s (t) is a complex number in the first quadrant. The amplitudes of the singularities, given here by the numbers 1 − λ + iǫ and its conjugate, are chosen so that the asymptotic form of the solution consists of one or two steadily propagating fingers of total width λ. The parameter ǫ determines the nature of the finger competition for B = 0.
The solution (7) describes generically two different fingers of the nonviscous phase penetrating the viscous one. In the linear regime |ζ s (t)| >> 1 the interface consists of a single bump or finger, and as time increases a second finger may develop and grow, depending on the value of Arg ζ s (0).
We summarize the features of the solution (7) that are most relevant to the study of finger competition. Consider first ǫ = 0. In this case the asymptotic configuration consists of one or two fingers of total width λ, depending on the initial condition. The singularities move toward the unit disk, with the limit as t → ∞ denoted by ζ s (t) → e iθ . When θ = 0 the asymptotic configuration is a single Saffman-Taylor finger growing in the center of the channel (this asymptotic configuration is denoted ST(R)), for θ = π/2 it is a 'side' Saffman-Taylor finger i.e. a pair of half fingers of total width λ with tips located at the cell walls (denoted ST(L)), and for θ = π/4 it is a 'double' Saffman-Taylor finger, namely two identical fingers of width λ/2 with tips at x = 0, ±1 (denoted 2ST). For any other value of θ the asymptotic configuration consists of two unequal steadily growing fingers. The two-finger asymptotic configuration is a consequence of the continuum of fixed points that is present in the phase portrait of the dynamical variables which govern the shape of the interface, namely (Reζ s (t), Imζ s (t)). In order to correspond to the notation of [21] . Then the planar interface corresponds to α = 0, the center Saffman-Taylor finger to α = −i, the side Saffman-Taylor finger to α = i, and the double Saffman-Taylor finger to α = 1. Figure 1 shows the phase portrait of the dynamical system obtained from the substitution of the mapping defined by Eq. (7) into the evolution Eqs. (5,6) for B = 0, using the dynamical variables (α ′ , α ′′ ) . The asymptotic states with two advancing fingers corresponds in the dynamical system (α ′ , α ′′ ) to a continuum of fixed points given by |α| = 1. Therefore, for ǫ = 0 the solution (7) does not exhibit finger competition. In addition, it is important to note that the evolution of (7) with ǫ = 0 is free of finite time singularities, i.e., z ζ = 0 in the domain |ζ| ≤ 1 for all time.
For ǫ = 0 the continuum of fixed points is removed (see Fig. 1), as is the double Saffman-Taylor finger fixed point. Consequently, the solution to Eq. (7) exhibits 'successful' competition, in the sense that the asymptotic interface shape consists of a single Saffman-Taylor finger or side Saffman-Taylor finger. The price to pay is the appearance of finite time singularities for a certain subset of initial conditions, in the form of a zero of z ζ impacting the unit disk (this is a generic feature of conformal map solutions z ζ composed of a finite number of pole singularities-see [9]). Then, only the subset of initial conditions free of finite time singularities is capable of sustaining finger competition all the way to the t → ∞ outcome. Nevertheless, one may ask whether the class of B = 0 solutions that are free of finite time singularities may describe, at least qualitatively, the physical finger competition for positive surface tension in the limit B → 0.
In the following sections we focus on the class of initial data for which the B = 0 solutions are devoid of finite time singularities, and examine the B > 0 dynamics. We develop a general theory for how the presence of positive surface tension affects the outcome of finger competition. This will enable us to predict the winner of finger competition, i.e., the eventual asymptotic state. Most interestingly, we find instances in which the presence of arbitrarily small surface tension leads to dramatically different outcomes in finger competition when compared with the zero surface tension evolution.

III. ASYMPTOTIC THEORY
Little is known about the effect of finite (but small) surface tension B on the dynamics of zero surface tension multifinger solutions, and in particular on the class of exact solutions (7). For single finger configurations, steady state selection theory predicts that the finger cannot have an arbitrary width. Indeed, for vanishing surface tension B → 0 the width λ = 1/2 is selected, asymptotically in time. Thus, it is clear that surface tension has a critical influence on single finger solutions with λ = 1/2. The nature of this influence in the limit B → 0 has been investigated by Siegel, Tanveer and Dai [19,20], who present evidence that zero surface tension single finger solutions with λ < 1/2 are significantly perturbed by the inclusion of an arbitrarily small amount of surface tension in order one time. The effect of surface tension is an increase of the finger width to reach the width predicted by selection theory.
Consider now the effect of small surface tension on the exact (B = 0) two finger solution (7). When 0 < B ≪ 1 the asymptotic perturbation theory developed in Refs. [13,19,20] can be applied. This perturbation theory describes the effects of the introduction of a small amount of surface tension on initial data z(ζ, 0) specified in the extended complex plane, i.e., in a domain including the 'unphysical' region |ζ| > 1 (the extended domain is required to make the B = 0 problem well-posed). The effect of finite B is most important near isolated zeros and singularities of z ζ (ζ, 0), where a regular perturbation expansion in B breaks down. (Away from these points the perturbation expansion is regular, at least initially.) For the class of solutions (7) we are discussing, the isolated singularities of z ζ (ζ, 0) are simple poles. The theory suggests that the introduction of finite surface tension modifies the poles (ζ s ) by transforming them into localized clusters of −4/3 singularities, but these clusters move at leading order according to the B = 0 dynamics. Thus the effect of one of these clusters on the interface is approximately equivalent to that of the unperturbed (B = 0) pole-like singularity that has given birth to it.
The influence of surface tension on the zeros of z ζ (ζ, 0) is more complex. Each initial zero instantly gives birth to two localized inner regions, i.e., regions where the B = 0 and B > 0 solutions differ by O(1). One of the two inner regions moves, at least initially, according to the B = 0 dynamics of the original zero ζ 0 [39]. Since the particular zero surface tension solutions considered here have zeros that are either bounded away from the unit disk for all time or impact the unit disk only after long times, the inner region around ζ 0 (t) has a negligible influence on the interface. The second inner region created around ζ 0 (0) moves differently. The theory suggests that this inner region consists of a cluster of singularities, whose size scales like B 1/3 . Unlike the case discussed above this second inner region moves away from the B = 0 zero since, to leading order in B, it moves like a singularity of the zero surface tension problem and this speed is different form the speed of the zero ζ 0 (t) which spawned the cluster. As this singularity cluster approaches the physical domain it may perturb the flow and the interface shape may differ significantly from that at B = 0 shape. The location of this singularity cluster will be denoted by ζ d (t), and following [13] we shall call it the daughter singularity. We emphasize that the dynamics of the daughter singularity cluster is determined at lowest order solely by the B = 0 solution z 0 (ζ, t), at least until it arrives at the surroundings of the unit circle, and therefore can be simply computed once the initial locations of the zeros of z ζ (ζ, 0) are determined.
The daughter singularity evolution equation is given by (see [13]) where q 0 1 is defined by and the superscript 0 denotes that the function evaluations are done using the corresponding B = 0 solution. The function −q 0 1 (ζ, t) also gives the characteristic velocity of a pole or branch point singularity of z ζ (ζ, t) located at position ζ in the region |ζ| > 1. The initial position ζ d (0) is a consequence of the fact that each zero ζ 0 (0) of the zero surface tension solution will give birth to a daughter singularity. From Eq. (8) it can be shown [13] that d|ζ d |/dt < 0, so that the daughter singularity approaches the unit circle and it can impact it in a finite time t d , the daughter singularity impact time, satisfyng |ζ d (t d )| = 1. In the limit B → 0, the daughter singularity impact time t d signals the time when the effects of the surface tension are felt on the physical interface. For times larger than t d the B = 0 interface and the B → 0 are expected to differ significantly.
For the family of exact B = 0 solutions the mapping function (7) has four pole-like singularities: ±ζ s and ±ζ s , and four zeros ±ζ 0+ and ±ζ 0− of z ζ located at For the particular case λ = 1/2 this solution presents only one pair of zeros ±ζ 0 located at In the following it will be useful to define the real quantity s which appears in (10) and (11).
Depending on the value of λ the initial data may have zeros on both the real and imaginary axes, or all the zeros may lie on a single axis. This difference has significant consequences in the finite surface tension dynamics. More specifically, when λ < 1/2 the zeros described in (10a) and (10b) are located on both the real and imaginary axes of |ζ| > 1, namely at ±|ζ 0+ | and ±i|ζ 0− |. The situation is different for λ > 1/2, which is further divided into two cases, depending on whether β 2 + 4(1 − 2λ)|ζ s | 2 > 0 or < 0. In the former case all four singularities lie on the real axis (for β > 0) or on the imaginary axis (for β < 0). In the latter case the four zeros are located off the axes in conjugate pairs, i.e. at ±ζ 0 and ±ζ 0 . Finally, when λ = 1/2 the solution (7) has only two zeros, located on the real axis at ±|ζ s | 2 / √ −2β when β < 0 and on the imaginary axis at ±|ζ s | 2 / √ 2β when β > 0. Note that for λ = 1/2 the B = 0 solution has two less zeros than for λ = 1/2. The initial zero locations described above have a critical bearing on whether the daughter singularity will impact the unit disk [40]. Although all daughter singularities approach the unit disk, their impact may be shielded by the presence of an inner region corresponding to a pole singularity. More precisely, since ζ d and ζ s obey the same dynamical equation, they will move together if they get close enough to each other. However, the inner region around a pole moves to leading order like the B = 0 pole, i.e., it moves exponentially slowly toward |ζ| = 1 when |ζ s | − 1 << 1, and does not impinge upon the unit disk in finite time [9]. In this case the O(B 1/3 ) inner region around the daughter singularity will not affect the dynamics on |ζ| = 1, at least until t = O(− ln B). Before this time, we expect the interface to be uninfluenced by the presence of the daughter singularity. This shielding mechanism is discussed in the context of single fingers in [20].
Knowledge of the t → ∞ asymptotic state and the initial locations of zeros can be used to ascertain whether shielding can occur. The B = 0 asymptotic state corresponds to ζ 2 s (t → ∞) → ±1. Thus, for λ < 1/2, only one pair of daughter singularities may be shieldednever both-so at least one pair of daughter singularities will impinge on the unit disk. The daughter singularities will also not be shielded when λ > 1/2 and β 2 + 4(1 − 2λ)|ζ s | 2 < 0. However, for λ > 1/2 and β 2 + 4(1 − 2λ)|ζ s | 2 > 0 it is possible for all the daughter singularities to be shielded, since they lie on a single axis. The daughter singularities can also be completely shielded when λ = 1/2. The different possibilities are schematically depicted in Fig. 2.
We have numerically computed the daughter singularity impact time t d for various values of λ and ǫ, using initial conditions close to the planar interface, |ζ s | 2 = 20 and various values of Arg[ζ 2 s ]. Figure 3 shows the phase portrait for different values of λ and ǫ with the daughter singularity impact indicated. From the plots it is immediately seen that for λ < 1/2 at least one daughter singularity always hits the unit circle, and for λ ≥ 1/2 some trajectories are free from daughter singularity impact. In addition, it is observed that for fixed λ a larger value of ǫ causes the daughter singularities to hit in shorter times (or less developed fingers) than a smaller value of ǫ, and for fixed ǫ larger λ implies larger impact times. We have also checked that the daughter singularity impact occurs well before a finite time singularity, i.e., the impact of a zero of z ζ . Thus, the effect of surface tension is significant well before the curvature in the zero surface tension solution becomes large. It is noted that the λ dependence of the daughter singularity impact is consistent with the results of steady state selection theory [15,16,17,18]. According to selection theory, for small B the possible values of λ are discretized: λ must satisfy the relation λ = λ n (B), given to leading order by where n parameterizes the branch of solutions. Note that λ n > 1/2 for all n. The steady finger shape is to leading order a Saffman-Taylor finger, with the above values of λ n substituted for the width λ. On the other hand for ǫ > 0 the asymptotic state of (7) is a Saffman-Taylor finger of width λ. From Eq. (12) it is clear there exists a steady solution with width λ n (B) close to a Saffman-Taylor finger of arbitrary width λ > 1/2. Thus the shielding of the daughter singularity, which leads to the persistence of a Saffman-Taylor solution with λ > 1/2 over long times, is consistent with steady state selection theory [41]. In contrast for λ < 1/2 there are no nearby steady solutions. Thus, a Saffman Taylor finger with λ < 1/2 cannot persist over a long time. We see that the impact of a daughter singularity provides a mechanism for the onset of finger competition, finger widening, and selection of a width λ > 1/2. For ǫ = 0 the scenario is similar, except there is an added class of exact B > 0 solutions. Magdaleno and Casademunt [36] have shown that two-finger solutions composed of steadily propagating but unequal fingers do exist for small nonzero B. The introduction of a small nonzero surface tension selects a discrete set of solutions from the continuum of fixed points of the B = 0 phase portrait. The solutions are parameterized by the total width of the fingers λ = λ 1 + λ 2 and the relative width q = λ 1 /λ, and the introduction of finite B discretizes the possible values of the parameters. In particular, they must satisfy a condition of the form λ = λ n (B) and q = q n,m (B) where n and m are integers. The expression for λ n (B) at lowest order is equivalent to Eq. (12), but with different coefficients C n . The shape of these solutions are given to leading order (in the limit t → ∞) by (7) with allowed value of λ n (B) substituted for the width λ. Again, λ n (B) > 1/2, and the consistency between daughter singularity impacts and steady state selection theory follows as above.
We conjecture that the outcome of interfacial shape evolution after the daughter singularity impinges is in general independent of the particular finger on which the impact first occurs i.e., independent of the point at which ζ d (t) impacts on |ζ| = 1. More specifically, we surmise that impact on either the shorter (trailing) or larger (leading) finger retards the velocity of that finger, and is accompanied by the widening of the leading finger, so as to maintain a constant fluid flux at infinity. The widened leading finger then shields the trailing finger, preventing it from further growth. Thus, the finger which is leading at the time of the daughter singularity impact 'wins' the competition, in the sense that it will evolve for t → ∞ to the ST steady finger. To examine this conjecture and study the dynamics of finger competition with finite (but small) surface tension we have numerically computed the evolution of an interface with initial conditions given by the conformal mapping Eq. (7) close to the planar interface (|ζ s (0)| −2 ≪ 1). The results are reported in the next section.

IV. NUMERICAL RESULTS
Numerical computations have been performed for B > 0, using an initial interface corresponding to the explicit B = 0 solutions discussed in Sec. II. The effect of positive surface tension on this class of solutions is explored for various values of ǫ and a variety of initial pole positions.
We employ the numerical method introduced by Hou et al. [32] and used in other studies of small surface tension effects in Hele-Shaw flow [19,20,33,34]. The method is described in detail in Ref. [32]. It is a boundary integral method in which the interface is parameterized at equally spaced points by means of an equal-arclength variable α. Thus, if s(α, t) measures arclength along the interface then the quantity s α (α, t) is independent of α and depends only on time. The interface is described using the tangent angle θ(α, t) and the interface length L(t), and these are the dynamical variables instead of the interface x and y positions. The evolution equations are written in terms of θ(α, t) and L(t) in such a way that the high-order terms, which are responsible of the numerical stiffness of the equations, appear linearly and with constant coefficients. This fact is exploited in the construction of an efficient numerical method, i.e., one that has no time step constraint associated with the surface tension term yet is explicit in Fourier space. We have used a linear propagator method that is second order in time, combined with a spectrally accurate spatial discretization. Results in this section are specified in terms of the scaled variables t = πt;B = π 2 B;x = πx;ỹ = πy. (13) instead of the original ones used in previous sections. The number of discretization points is chosen so that all Fourier modes of θ(α, t) with amplitude greater than round-off are well resolved, and as soon as the amplitude of the highest-wavenumber mode becomes larger than the filter level the number of modes is increased, with the amplitude of the additional modes initially set to zero. The time step ∆t is decreased until an additional decrease does not change the solution to plotting accuracy, nor lead to any significant differences in any quantities of interest. In a typical calculation 512 discretization points are initially used, and the initial time step is ∆t = 5 · 10 −4 . For small values of surface tension numerical noise is a major problem, and the spurious growth of short-wavelength modes induced by round-off error must be controlled. To help prevent this noise-induced growth at short wavelengths spectral filtering [37] is applied. Additionally, we minimize noise effects and also assess the time at which these effects become prevalent by employing extended precision calculations, as described in the next section.
Our main interest is to uncover the role of surface tension in the dynamics of finger competition. To isolate the features of finger competition from those of width selection, we will concentrate on B = 0 solutions with λ = 1/2, the value selected by surface tension in the limit B → 0. Since the B = 0 dynamics for ǫ = 0 and ǫ = 0 is quite different the numerical results for the two cases will be presented separately.

A. Solutions with ǫ = 0
We first consider parameter values λ = 1/2 and ǫ = 0. A typical set of interfacial profiles is shown in Fig. 4. The initial data is given by the mapping function Eq. (7), with λ = 1/2, ǫ = 0, d(0) = 0 and ζ 2 s (0) = 20 exp(iπ/6). With this value of ζ 2 s (0) the initial interface is well inside the linear regime. Evolutions are shown for different values ofB, and the B = 0 interface evolution is also plotted for comparison. In all these evolutions the filter level is set to 10 −13 , although later we shall make comparisons to profiles computed at higher precision.
For the largest value of surface tension the computed B > 0 and the exact B = 0 solutions first differ appreciably at the seventh curve, corresponding tot ≈ 3. At this point the velocity of the small finger (at the channel sides) begins to decrease and it is clearly left behind when compared with the small finger evolution in the B = 0 solution. Eventually, the advance of the small finger is completely suppressed and the larger finger widens to attain a width close to 1/2 of the channel. For a smaller value of surface tension, for instanceB = 0.001, the evolution displays qualitatively the same behavior.
The B > 0 interface differs appreciably from the B = 0 sightly later than before (i.e., at the eighth curve) and the region where the two solutions differ most is to some extent more localized around the small finger than for larger values of B. Additionally, for this value of surface tension the effect of numerical noise is clearly exhibited in the interfacial profiles. Here the tip-splitting and sidebranching activities are a clear effect of numerical noise, as can be easily checked redoing the computation with a different noise filter level.
In order to suppress or delay the branching induced by numerical noise that appears for small values of surface tension it is necessary to use higher precision arithmetic, e.g. quadruple precision (128-bit arithmetic). The filter level can then be reduced by a large amount and the outcome of spurious oscillations is substantially delayed. Figure 5 shows the effect of reducing the filter level to 10 −27 . The B = 0 solution is plotted, as well as the computation with double precision. ForB = 0.001 the branching is totally suppressed, at least for the times we have computed, but for smaller values ofB the use of quadruple precision is only able to delay the branching and not totally suppress it. The quadruple precision computation confirms the results observed with lower precision: the introduction of finite (but small) surface tension results in the suppression of the small finger. From Fig. 5 one can also see that for long times, when the interface is clearly affected by numerical noise (in the double precision curve), the noise-induced branching is restricted to the large finger, and the small finger is basically unaffected by noise. This observation suggests that the small finger shape, as well as its tip velocity and tip curvature, can be trusted even when the large finger has developed tip-splittings and side-branchings due to the spurious growth of round-off error. Figure 6 shows the tip velocity of both fingers versust for decreasing values of surface tension. It can be seen that the velocity of the large finger is only slightly affected by surface tension, whereas the velocity of the small finger is substantially reduced by the inclusion of finite B. AsB is decreased the tip velocity of the small finger is more faithful to the B = 0 evolution before the daughter singularity impact (shown by a cross), and clearly veers away from the B = 0 velocity later in the evolution, consistent with asymptotic theory. Note that at the smallest value ofB the tip velocity of the large finger drastically differs from the B = 0 velocity at late times. This discrepancy is a manifestation of noise effects in the neighborhood of the large finger tip. However, as previously seen, the small finger is basically unaffected by noise at the times we have plotted.
In order to further verify that the daughter singularity impact is responsible for the observed change in the small finger tip speed we follow the scheme introduced in [20]. Define t p as the time when the computed tip velocity differs by p from the B = 0 tip velocity. According to asymptotic theory this t p will be a linear function of B 1/3 in the limit B → 0 as long as p is small enough.  showst p versusB 1/3 for various values of p, and it can be seen thatt p exhibits the predicted behavior. Moreover, we have extrapolated the B = 0 value oft p using the two points of lowestB and the result is very close tõ t d , whose value is represented by a cross. We conclude that the impact of the daughter singularity is associated with the dramatic change of the B > 0 solution when compared to the zero surface tension solution, reducing the velocity of the small finger and eventually suppressing it. In contrast, for the B = 0 dynamics the small finger 'survives', propagating with the same asymptotic speed as the larger finger. Note that the average interface advances at unit velocity, and a tip velocity below 1 implies that the finger is retreating in the reference frame of the average interface.
In summary then, our numerical results show that the computed interface for B = 0 follows the B = 0 evolution for an O(1) time interval-roughly corresponding to the daughter singularity impact time-and that at further times the velocity of the small finger decreases while the large finger widens. The small finger eventually comes to a halt and the larger (leading) finger reaches an asymptotic width slightly above 1/2, the width singled out by selection theory. It is noted that for the initial condition we have studied the daughter singularity impact takes place on the tip of the small finger. Therefore, the influence of surface tension on the interface should be sig-nificant first around the impact point, that is, the small finger tip. Our numerical results show that in fact this is the case; the initial effect of the daughter singularity impact is to slow and then completely stop the growth of the small finger. Later on, as the singularity cluster centered in ζ d spreads over the unit circle, the effect of surface tension is felt by the whole interface and the large finger widens to reach the selected width.
We have also studied the finite surface tension dynamics for a more general class of initial conditions. More precisely, we have studied initial conditions of the form ζ 2 s (0) = 20 exp(i nπ/12) where n = 0, ±1, ..., ±6, and have obtained the same qualitative results as in the case previously studied, namely that the presence of small surface tension suppresses the growth of the finger which is trailing at the time of daughter singularity impact. In order to compare the B = 0 and the B = 0 dynamics in a compact and global way we have plotted the phase portrait for B = 0 using the the tip velocities v 1 , v 2 as dynamical variables. In the laboratory frame they read Now a comparison between dynamics for B = 0 and B = 0 is straightforward since the trajectories can be   plotted together and compared. In addition, the tip velocity is a useful variable because it contains geometric information; specifically the inverse of the tip velocity is equal to the width of the finger in the asymptotic (t → ∞) regime. It is important to note that (v 1 , v 2 ) are dynamical variables for the B = 0 problem, so that the plot of the zero surface tension trajectories onto the space (v 1 , v 2 ) is a true phase portrait. On the other hand (v 1 , v 2 ) are not state variables of the problem with finite surface tension, so in this case we simply obtain a projection onto the (v 1 , v 2 ) space of the original B = 0 trajectory, which is embeded in the infinite-dimensional phase space of interface configurations. Figure 8 shows the phase portrait for B = 0 together with the tip velocities obtained from the initial conditions described above forB = 0.01. From the figure it is evident that the introduction of finite surface tension has substantially changed the global phase dynamics of the problem. Only oneB = 0.01 trajectory connects the planar interface (1, 1) and the 2ST point (2, 2), corresponding to the unsteady double Saffman-Taylor finger. Any otherB = 0.01 trajectory ends in one of the two ST finger points, ST(L) at (2, 0) and ST(R) at (0, 2). In contrast, the (2, 2) point, equivalent to the continuum of fixed points present with the (α ′ , α ′′ ) or (Reζ s , Imζ s ) variables, has a finite basin of attraction for B = 0. The introduction of finite surface tension has dramatically changed the zero surface tension (v 1 , v 2 ) trajectories, to the extent that the B = 0 phase portrait and the B = 0 projection are not topologically equivalent. This result is not a complete surprise, since it was anticipated from the structural instability of the dynamical system governing the evolution of Eq. (7) for ǫ = 0 [21]. A more dramatic example of topological inequivalence of phase portraits will be given in the next subsection, when we consider the case ǫ = 0.  Although the use of the variables (v 1 , v 2 ) has allowed us to project the finite surface dynamics onto the zero surface tension phase portrait this projection has one major limitation: it only considers a local quantity, the tip velocity. We have also considered a projection that takes more global properties of the interface into account. Specifically, given a computed B = 0 solution for an initial condition of the form (7), one can use a suitable norm to define a 'distance' between the computed interface and the B = 0 interface obtained from the mapping function Eq. (7). We choose this 'distance' to be the area enclosed between the two interfaces at a given time. Additionally, we define a projection of the B = 0 interface onto the B = 0 phase space (with phase space variables (Re ζ s , Im ζ s )) by selecting the value of ζ s that minimizes the 'distance' between the two interfaces, with the restriction that the position of the two mean interfaces must be the same. The latter condition ensures that the projection satisfies mass conservation. Figure 9 shows the B = 0 phase portrait and the corresponding projected evolution for surface tensionB = 0.01. Again, the plot clearly shows that the introduction of finite surface tension modifies the phase portrait of B = 0. The projected trajectories are initially close to the B = 0 dynamics, but for well developed fingers (corresponding to |α| ∼ 1) the projection departs from the B = 0 trajectory towards the Saffman-Taylor fixed point, located at α ′ = 0, α ′′ = 1. The projected trajectory only remains close to the corresponding B = 0 trajectory when the latter evolves towards the Saffman-Taylor fixed point. More precisely, the continuum of fixed points present for B = 0 has been removed by surface tension and the Saffman-Taylor fixed point is the universal attractor of the dynamics for finite surface tension.
In Fig. 10 the projection for decreasing values ofB is plotted, using the initial condition ζ 2 s (0) = 20 exp(iπ/6). AsB is decreased the projected trajectory gets closer to the B = 0 trajectory, but as it approaches the point when the daughter singularity impinges the unit circle (this point is signaled by a cross) the projection departs from the B = 0 trajectory and approaches the Saffman-Taylor fixed point, consistent with asymptotic theory.

B. Solutions with ǫ = 0
The continuum of fixed points present for ǫ = 0 is absent for ǫ = 0, but in this case finite-time singularities in the form of zeros of z ζ impinging on the unit disk do appear for some initial conditions. Therefore, we can expect that the effect of finite surface tension will be somewhat different than for ǫ = 0. Firstly, the presence of surface tension should eliminate finite-time singularities, and secondly, finite B could modify the basin of attraction for the two attractors of the B = 0 dynamical system, namely the side Saffman-Taylor finger and the center Saffman-Taylor finger.
To explore this, we have performed computations with From the plot one can see that mostB = 0.01 velocity trajectories follow (at least qualitatively) their B = 0 counterparts, in the sense that they end up in the same fixed point. However, the second, third and fourth trajectories (counting from the upper left trajectory in clockwise direction) differ significantly from their B = 0 counterparts. The secondB = 0.01 trajectory moves apart from the B = 0 solution simply because the latter develops a finite-time singularity, which is regularized by the introduction of finite surface tension. However, the third and fourth trajectories exhibit a quite surprising behavior: the computed interface withB = 0.01 ends up in a different fixed point than the exact B = 0 solution, despite the fact that the B = 0 solution is smooth for all time and has the asymptotic width that would be selected by vanishing surface tension.
In order to get further insight into this behavior we have computed the evolution for decreasing values ofB using the specific initial pole position ζ 2 s (0) = 20 exp(−iπ/6), with λ = 1/2 and ǫ = 0.1. Quadruple precision has been used when it has been necessary. Figure 12 shows its evolution for four values of the surface tension parameter, together with the B = 0 solution. The differences between the two interfaces for long times are readily apparent. When B = 0 the finger in the central position stops growing and the side finger wins the competition, whereas for B > 0 we encounter the opposite situation-namely, the central finger surpasses the side finger and wins the competition. For the smaller values of B the finger on the sides has not quite stopped growing when the computation is stopped, although its tip speed shows a marked decrease over that for B = 0 and is less than that of the central finger. The side finger tip speed is also decreasing at the final stage of the computation. The tip speed trend in the limit B → 0 is further illustrated in Fig. (13). This figure shows the tip speed versus time of each finger for a sequence of decreasing B. The plot suggests that upon impact of the daughter singularity the side finger velocity levels off and eventually decreases, whereas the velocity of the center finger is nearly unaffected and continues to increase. The trend is indicative of the center finger "winning" the com- petition in the B > 0 dynamics, while the opposite occurs for B = 0. Finally, it is noted that the influence of surface tension is first felt by the smaller finger, which is the recipient of the daughter singularity impact. Afterwards the leading finger begins to widen, in a manner consistent with the conjecture in Sec. III. Further remarks on this point are made in Sec. V.
The projection method described in the previous section has been also applied to this case, and the results are displayed in Fig. 14 in the particular caseB = 0.01. It can be seen that for most trajectories the projection stays close to the B = 0 curves, even for long times. The daughter singularity impact still leads to O(1) differences between the B = 0 and B > 0 solutions, although the impact does not produce changes in the outcome of finger competition. However, as expected some of the trajectories (namely the third and fourth as measured clockwise from the bottom) do indicate significant qualitative differences in the long time evolution. The plot provides a simple depiction of the topological inequivalence of the B > 0 and B = 0 dynamics [42].
It has been shown that the introduction of a finite B has not changed the attractors of the problem, but it has changed their basins of attraction. Interestingly, in the B = 0 case there does not exist a single separatrix trajectory between the two Saffman-Taylor attractors, but rather a finite region, corresponding to the set of trajec-tories ending in cusps, that acts as an effective separatrix. Since for finite surface tension there are no cusps, it can be assumed that there is a single trajectory that separates the two basins of attraction. Obviously, this trajectory will depend on the value of the surface tension parameter. More precisely, the initial condition ζ 2 s (0) corresponding to the separatrix trajectory will be a function of the surface tension B. To quantitatively characterize this set of initial conditions we have studied the dependence of the separatrix trajectory in a neighborhood of the planar interface fixed point as a function ofB, using initial conditions of the form ζ 2 s (0) = 20 exp(iθ). For a given initial condition ζ s (0) introduce the parameter θ sep (B), defined as the unique value for which the evolution is attracted toward the fixed point ST(L) when θ > θ sep and to the fixed point ST(R) when θ < θ sep . Figure 15 shows the plot of θ sep versusB, and it is observed that asB decreases, θ sep saturates to a fixed value, namely θ sep (B → 0) = −0.4843 ± 0.0009. It is interesting to compare this value to the position of the separatrix region for B = 0, which is located between θ B=0 + = −0.95758 and θ B=0 − = −1.04796. The separatrix for finiteB lays outside and far away from the separatrix region for B = 0, even for vanishing surface tension. Our evidence therefore suggests that any B = 0 trajectory located between the trajectories defined by θ sep (B → 0) and θ B=0 + will not describe, even qualitatively, the reg- ularized dynamics in the limitB → 0, since the finger that will 'win' the competition under the B = 0 dynamics will 'lose' under the B → 0 dynamics. Thus, there exists a positive measure set of initial conditions of the form (7) such that the evolution with B → 0 cannot be qualitatively described by its evolution under B = 0 dynamics. This is a dramatic consequence of the singular nature of surface tension on the dynamics of finger competition which is not related to steady state selection, but confirms the ideas of the proposed dynamical solvability scenario in Ref. [22].

V. SUMMARY AND CONCLUDING REMARKS
The asymptotic theory developed in Refs. [13,20] predicts the existence of regions of the complex plane where the zero surface tension solution and the finite surface tension solution differ by O(1). These regions are the daughter singularity clusters, and their influence is felt in the physical interface when they are close to the unit circle. Daughter singularities move towards the unit circle, and when their motion is not impeded by other singularities they reach the unit circle in O(1) time. When the distance between the daughter singularity and the unit circle is O(B 1/3 ) the interface can display O(1) discrepancies with respect the interface of the B = 0 solutions. However, the asymptotic theory does not predict the nature of the discrepancies caused by daughter singularity impact.
Siegel et al. [20] showed numerically that the effect of the daughter singularity impact on the tip (in a singlefinger configuration with λ < 1/2) was a retard in the velocity of the finger accompanied by a widening of the finger. However, this provided small insight on the effect of the impact in multi-finger configurations, where finger competition could be substantially affected by the presence of finite surface tension, as suggested in Refs. [5,21,22].
Since the precise effect of the daughter singularity cannot be established by the asymptotic theory it is necessary to use numerical computation in order to establish the effects of daughter singularity on the dynamics of the interface. We have focused our efforts on uncovering the role of surface tension in the dynamics of two finger configurations, which is the simplest situation exhibiting nontrivial finger competition. Two different types of twofinger zero surface tension solutions have been studied. The first type (ǫ = 0) does not exhibit finger competition when B = 0 but rather contains asymptotic configurations consisting of two unequal steady fingers advancing with the same speed. These two-finger steady state solutions form a continuum of fixed points in the phase space of the corresponding, reduced dynamical system, which is structurally unstable. Numerical computations with small surface tension show that the introduction of a small B removes the continuum of fixed points and triggers the competition process which was absent for B = 0 by restoring the saddle-point (hyperbolic) structure of the appropriate multifinger fixed point. The second type (ǫ = 0) of two-finger solution we have studied exhibits finger competition for B = 0, but the numerical computation with small B has shown that the long time configuration of the computed interface may be qualitatively different from the B = 0 solution for a broad set of initial conditions, in the sense that the finger that 'wins' the competition is not the same with and without surface tension. Thus, the presence of surface tension seemingly can change the outcome of finger competition even in configurations that are well behaved and smooth for all time and whose asymptotic width is fully compatible with the predictions of selection theory for vanishing surface tension. This unexpected result shows that surface tension is not only necessary to select the asymptotic width and to prevent cusp formation, but plays also an essential role in multifinger dynamics through a drastic reconfiguration of the phase space flow structure.
Our calculations support the conjecture that impact on either the shorter or larger finger retards the velocity of that finger, and is accompanied by the widening of the larger finger. As a consequence, in general the outcome of finger competition is independent of the particular finger on which the impact first occurs, and the finger which is leading at the time of the daughter singularity impact 'wins' the competition. This recipe fails only for interfacial configurations with very similar fingers, when not only the position of the finger (which finger is leading) but also the tip velocities (a trailing finger can have for a certain time a larger velocity than the leading one) at the impact time may play a role.
The main conclusion of the present work is that surface tension is essential to describe multifinger dynamics and finger competition, even when the corresponding zero surface tension evolution is well behaved and compatible with selection theory. That is, we have detected singular effects of surface tension on the dynamics of finger competition that are not directly related to steady state selection. These can be properly interpreted in the context of an extended dynamical selection scenario as described in Ref. [22] where the reconfiguration of phase space flow by surface tension can be traced back to the restoring of hyperbolicity of multifinger fixed points.