Systematic weakly nonlinear analysis of radial viscous fingering

We present a weakly nonlinear analysis of the interface dynamics in a radial Hele-Shaw cell driven by both injection and rotation. We extend the systematic expansion introduced in @E. Alvarez-Lacalleet al., Phys. Rev. E 64, 016302~2001!# to the radial geometry, and compute explicitly the first nonlinear contributions. We also find the necessary and sufficient condition for the uniform convergence of the nonlinear expansion. Within this region of convergence, the analytical predictions at low orders are compared satisfactorily to exact solutions and numerical integration of the problem. This is particularly remarkable in configurations ~with no counterpart in the channel geometry ! for which the interplay between injection and rotation allows that condition to be satisfied at all times. In the case of the purely centrifugal forcing we demonstrate that nonlinear couplings make the interface more unstable for lower viscosity contrast between the fluids.


I. INTRODUCTION
The evolution of an unstable fluid interface in a Hele-Shaw cell ͓1,2͔ has been widely studied, both theoretically and experimentally.While it exhibits the inherent difficulties of a free-boundary problem with nonsaturated, nonlinear, and nonlocal growth, yet it is simple enough to allow for analytical insight ͓3-6͔.The channel and radial geometries are the two basic configurations.
In the channel geometry the interface between two fluids is driven by an external pressure ͑with different viscosities of the fluids͒ or by gravity ͑with different densities͒ ͓2,3,7͔.The gravity-driven and the pressure-driven problems can be mapped into each other in the appropriate reference frame, leaving only two dimensionless parameters, viscosity contrast A and dimensionless surface tension B. A thorough study of the problem in channel geometry has been pursued since Saffman and Taylor ͓3͔ explained the instability mechanisms and obtained the analytical shape of steady finger solutions for Bϭ0.In particular, the singular perturbation character of surface tension ͓8͔ has been identified as source of different subtle effects.The most celebrated one is the mechanism of steady-state selection ͓9-11͔, but more recently, surprising singular effects on the dynamics of finger competition have also been unveiled ͓12-16͔.On the other hand, viscosity contrast also plays an important role in the fingering dynamics.For instance, it has been demonstrated that the finger competition is inhibited when the viscosity contrast is low ͓7,17-20͔ and that consequently, the basin of attraction of the single-finger solution does depend on viscosity contrast ͓21,22͔.
The radial configuration has also been object of intensive experimental and theoretical analysis.In the typical radial configuration ͓23,24͔, the circular Hele-Shaw cell is filled with two fluids.The less viscous fluid is injected at the center of the cell and displaces the more viscous one.The interface is unstable due to viscosity contrast, as in the pressure-driven channel instability.To provide an instability independent of the viscosity contrast in radial geometry, a counterpart of the gravity-driven experiment is realized by introducing a centrifugal force.This situation has been extensively studied experimentally ͓25-27͔ and also theoretically ͓28 -31͔.
A basic and important feature distinguishes the radial geometry from the channel one.In the former, the injection and the rotation forcings are not equivalent due to the different dependence of viscous and centrifugal forces with radial distance.Injection and rotation drivings thus give rise to different interface dynamics and the problem involves three dimensionless parameters.Viscosity contrast A and dimensionless surface tension B are now supplemented with the ratio of the two driving forces.It has been shown experimentally that complicated finger morphologies appear due to tip-splitting and screening effects ͓23,24͔ when injection is the unique destabilizing force.On the contrary, a purely centrifugally driven experiment leads to a different morphology of the fingers and the enhancement of pinch-off singularities, especially when the viscosity contrast is low ͓25,32͔.
The rich experimental phenomena and the specific features of the radial geometry indicated above suggest the relevance of extending to this geometry the weakly nonlinear analysis of the viscous fingering problem developed in Ref.

͓33͔.
The main reason for developing the weakly nonlinear analysis is the possibility to extract analytical information which is not perturbative in the relevant parameters of the problem.For instance, the region of finite B, which is relevant in the typical configuration emerging from the linear instability, can be explored naturally although, in the unstable configuration, the analysis is limited to a transient.The interplay of rotation and injection in radial geometry can also be best understood by analyzing how the corresponding terms interact in the early nonlinear regime.Similarly, the weakly nonlinear analysis provides explicit nonperturbative information about the connection between centrifugal forces and viscosity contrast at the nonlinear level in the purely centrifugally driven case.
Remarkably enough, in the radial geometry one can find situations which combine rotation and injection in such a way that the early unstable growth is stabilized at long times.It is then possible to apply the weakly nonlinear analysis to the entire interface evolution, giving rise to a reduced description of the whole Saffman-Taylor dynamics in terms of a set of ordinary differential equations.This situation has no counterpart in channel geometry.
Finally, the weakly nonlinear analysis is also relevant in the analysis of a stably stratified interface which is driven by an external force such as modulations of the gap or of the wetting condition ͑see Ref. ͓34͔ for quenched disorder͒.Specifically, the information provided by the weakly nonlinear analysis in channel geometry has already proven useful in morphologically stable interfaces subject to a local forcing in the form of random cell gap variations ͓35͔.In such cases, the lowest order of the nonlinear terms is required to obtain the long-wavelength, low-frequency scaling properties of the interface fluctuations.While nonlinear couplings that are local in space are relatively easy to obtain, this is not the case for nonlocal terms.The weakly nonlinear analysis is then a useful tool to obtain the complete set of nonlinearities to the desired order, to discuss, for instance, the universal properties of interface growth ͓35͔.The present analysis will be useful to extend these studies to the radial geometry.
The aim of this paper is to work out the details of the expansion in radial geometry.We will explicitly carry out the nonlinear expansion up to second-order couplings, including both rotation and injection, which extends the result of Ref.
͓36͔ to the case with rotation.We will extract analytical information concerning the interplay between injection, rotation, and viscosity contrast at different orders of nonlinear couplings.We will also be interested in obtaining an exact criterion of convergence of the nonlinear expansion in the radial geometry, extending the analysis of the channel geometry.
The approximate evolution obtained with the nonlinear analysis will be compared with numerical simulations with finite surface tension.The nonlinear expansion is also applied to the zero surface tension case which is special in that explicit time-dependent solutions are known, even for the combined case of rotation and injection ͓29͔, and therefore provide an additional ground for testing the method.Some of the physical anomalies of such solutions will be clearly manifested within the weakly nonlinear scheme.
The layout of the rest of the paper is as follows: in Sec.II we introduce the formalism and obtain the weakly nonlinear equations for Hele-Shaw flows in radial geometry.We also address the convergence criterion in this section.Its mathematical proof is given in Appendix.Section III presents a numerical analysis of exact, simulated, and approximate solutions.The main results and the conclusions are summarized in Sec.IV.

A. Vortex sheet formalism in radial geometry
Let us consider the Hele-Shaw problem in the radial geometry.The initial interface is a constant circle RϭR 0 which separates an outer fluid ͑labeled 1͒ and an inner fluid ͑labeled 2͒, which have known dynamic viscosities 1 , 2 and densities 1 , 2 .The gap between the plates is b, Q is the areal injection rate, and ⍀ is the angular velocity of the cell ͑Fig.1͒.
The equations of motion for the interface and the boundary conditions are well known ͓3͔.Following the same approach as that in Ref. ͓33͔ and using the formulation of Tryggvason and Aref ͓7͔, we introduce velocity U ជ ϭ(u ជ 1 ϩu ជ 2 )/2 where u ជ 1 and u ជ 2 are the two limiting values ͑from both sides of the interface͒ of the solenoidal part of the velocity at a given point.This velocity U ជ can be expressed in terms of vortex sheet distribution ␥ at the interface as where ␥ ˜ϭͱ1ϩ(r /r) 2 (u ជ 1 Ϫu ជ 2 )•s ˆ, and we have used nota- tion r( 1 ,t)ϵr 1 , r( 2 ,t)ϵr 2 .
In the presence of sinks or sources, velocities u ជ 1 and u ជ 2 which define U ជ include only the solenoidal part of the total velocity field.For this reason, when injection QϾ0 or withdrawal QϽ0 of the inner fluid is present, Eqs.͑2.1͒ and ͑2.2͒ must be supplemented with the corresponding irrotational part of the velocity field.
In order to obtain the expression for the vorticity as a function of u ជ i , we use Darcy's law for the gap-averaged velocity ͓25͔: ٌ ជ p i ϭϪ and boundary conditions After some algebra, we can write an expression for the vorticity as a function of U ជ : with the curvature and the viscosity contrast as We now impose the continuity of the normal velocity.The projection of the radial velocity along normal direction n ˆhas two contributions: the solenoidal part of the average velocity, U ជ , and the irrotational part of the mean interface velocity, U ជ irrot :

B. Systematic weakly nonlinear analysis: Radial geometry
Our goal in this section is to introduce a systematic method to derive an evolution equation of the interface in real space, up to a given order in nonlinear couplings, in the radial geometry.As in Ref. ͓33͔, the different orders of mode couplings will be ordered as powers of a ''book-keeping'' parameter , to be defined below.The evolution of the interface will thus take the following form: where F͓r͔,G͓r͔, . . .are nonlocal operators on function r(,t), including nonlinearities of order nϩ1 in the term of order n .
To define small parameter , we first split function r(,t) in two parts, namely, the radius of the unperturbed circle and the deviation from this circle: r͑,t ͒ϭͱR 0 2 ϩ Qt ϩr͑,t ͒ϭR͑ t ͒ϩr͑ ,t͒.͑2.9͒ The largest length scale in the problem is given by R(t), which is a constant when Qϭ0.For finite injection rate, R(t) defines an evolving ͑unstable͒ solution.We could naively define as ϭw/R, w being the typical scale of r(,t).This straightforward generalization of the scaling procedure used in the channel geometry is not appropriate here.The reason was already pointed out by Miranda and Widom ͓36͔.In the radial geometry, mass conservation implies that the zero mode ͑which, in the channel geometry, is decoupled from the rest and drops out of the formulation͒ has a higher-order nonzero amplitude, since it must satisfy a mass conservation relation which, to lowest order, reads: Since the zero mode is explicitly one order higher than the others, we must consider a scaling that does not mix different orders of the weakly nonlinear analysis.We use the following scaling: with r͑ ͒ϭr ˜͑ ͒ϩr 0 , ͑2.12͒ where r 0 stands for the zero mode and r ˜() is the sum of all the other modes.To simplify the notation we have dropped out the time dependence.
We define the characteristic scaling velocity as: implying that time will be scaled with R/V.Once Eqs.͑2.1͒, ͑2.2͒, ͑2.5͒, and ͑2.7͒ are made dimensionless, we have for the vorticity, where

͑2.15͒
B is the dimensionless surface tension, C is the scaled injection, and D the scaled rotation.For the integrals we have and the dimensionless equation for the evolution of r()

͑2.19͒
This last Eq.͑2.19͒ is the only one which is sensitive to the different scalings in r().This equation has a time derivative which produces different results depending on how R(t) is involved in the scaling.The scaling that we propose is the only one that is truly systematic.Other possibilities are not incorrect but will mix different orders of the weakly nonlinear expansion at a given order in .
Finally, to complete the theoretical analysis we address the convergence of the weakly nonlinear expansion ͑2.8͒.For the channel geometry, we showed in Ref. ͓33͔ that the expansion in converges uniformly if and only if ͉h x ͉Ͻ1 in the whole domain of integration.Here, we find a similar condition in terms of nonscaled variable r.Whenever an interface fulfills the nonlinear expansion converges uniformly.This result is proved in Appendix and it reduces to ͉h x ͉Ͻ1 in the limit of channel geometry (R→ϱ).This result provides a practical criterion to assess the validity of the expansion in different situations.

C. The linear dispersion relation
We want to obtain the linear dispersion relation, in both real and Fourier space, to discuss the interplay between rotation and injection at the linear stage of the instability.
We will use the following definitions for the average and the Hilbert transform of a function f in the unit circle: where subindex in operator H indicates the argument of the Hilbert transform, not a derivative.The zero orders of vorticity and U ជ are: and thus This is different from its counterpart in the channel geometry ͓33͔.Here, the equations for the vorticity at consecutive orders require the knowledge of the average vorticity.The equations are solved by averaging on both sides.For example, at zero order, A 1 will lead to ␥ ˜(0) ϭ0 but for A ϭ1, ␥ ˜(0) can be any constant.The presence of an arbitrary constant is not a problem, since it is eliminated at each order in the equation for the evolution of the interface.
Expanding the equations up to first order and using ␥ ˜(0) ϭ0, the linear equation for the evolution of the interface deviation reads We can now perform a Fourier transform of the equations and, after reintroducing the adequate dimensions, we reproduce the linear dispersion relation found in Ref. ͓25͔: We now return to the question of the different scaling alternatives introduced at the end of Sec.II B. Had we used scaling ϭw/R 2 (t), term Q/(2R 2 ) in Eq. ͑2.26͒ would have disappeared before reintroducing the dimensions.On the contrary, for ϭw this term becomes multiplied by a factor 2. These changes in the dimensionless linear dispersion relation come from the fact that a mode ␦ n in the scaling adopted in Eqs.͑2.11͒ and ͑2.12͒ becomes ␦ n /R in the first of the scaling alternatives above.
Term Q/(2R 2 ) at the end of Eq. ͑2.26͒ is stabilizing for positive injection rate (QϾ0) and destabilizing for QϽ0.In a sense, this is a purely geometric effect, which contributes to the growth or decay of modes due to the expansion or contraction of the base state.As pointed out in Ref. ͓29͔, term Q/(2R 2 ), that can be scaled out, must be distinguished from the other terms that describe the intrinsic instability of the problem.Disregarding that last term, it is clear that injection and rotation (CV 0 and DV 0 , respectively͒ play an equivalent role in the linear dispersion relation, since both appear multiplied by k.By choosing the proper viscosity and density contrasts, they can produce exactly the same stabilizing or destabilizing effect.This is the counterpart, in the radial geometry, of the equivalence between the roles of injection and gravity in the channel geometry.In the channel geometry, however, the equivalence is exact and can thus be extended to the whole nonlinear evolution ͑in the appropriate reference frame͒.This is not the case in radial geometry, as shown in the following section.

D. First weakly nonlinear order
The purpose of this section is to derive the leading nonlinear contribution for the radial Hele-Shaw cell with both rotation and injection.We first recall that mass conservation relates the zero mode with the other modes.In the original nonscaled variables this reads To reproduce this result and the nonlinear couplings for the other modes we start with the equation for the interface at order 2 : where we have used that (2) can be obtained from the expansion of U r where (2) Ϫ4Cr ˜r ˜ .͑2.30͒ (2) and U (2) can be computed from Eqs. ͑2.15͒ and ͑2.17͒, respectively.The latter involves only ␥ ˜(1) .We obtain finally ϩ2r ˜r ˜ , s ϭ͓B͑r ˜ϩr ˜ ͒ϩ͑ DϪC ͒r ˜͔ .

͑2.32͒
By introducing a superposition of Fourier modes in r ˜we directly recover Eq. ͑2.27͒ for mode kϭ0.For the other modes (k 0) we get the following result: where where sgn͑ ͒ is the sign function.For the particular case ⍀ ϭ0, the expression above reproduces the result of Miranda and Widom ͓36͔.We find that the presence of rotation does not change the formal structure of the equations.
To emphasize the roles of injection and rotation, we define quantity H(k,p) as the part of the coupling matrix in the right hand side ͑rhs͒ of Eq. ͑2.33͒ which contains Q and/or ⍀.This quantity reads:

͑2.37͒
We observe that experimental parameters occur only in two groups, namely, Both of them turn out to multiply the same functions of the wave numbers, except for changes in sign.However, notice that the relative sign of ⍀ ˜and Q ˜is different in the linear part and in the leading nonlinear part of the dynamics.
Considering the approach taken in the preceding section, where only the morphological instabilities are taken into account, there are only two important parameters in the linear regime: viscosity contrast A and ratio B/(CϩD) between the capillary forces and the sum of the injection and rotation forces.On the contrary, the weakly nonlinear solution shows that another important parameter must be taken into account, namely, coefficient Q ˜/⍀ ˜.The overall dynamics presents three different parameters A, B/(CϩD), and C/D.This clearly shows that the effects of injection and rotation are not equivalent, and that the intrinsic dimensionless parameter related to rotation already shows up at the early nonlinear regime.
Finally, notice that at nonlinear order there is also a geometrical term Q/(2R 3 ) independent of k.This term cannot be scaled by the same procedure used at the linear order with analogue term Q/(2R 2 ).

III. ANALYSIS OF EXACT SOLUTIONS AND THEIR WEAKLY NONLINEAR APPROXIMATIONS
We proceed to compare the approximate results that the weakly nonlinear equations ͓Eq.͑2.33͔͒ provide, with exact solutions of the full problem obtained analytically or numerically.Our purpose is to gain further information about the implementation of the weakly nonlinear analysis and to discuss its accuracy.
We will focus on three comparisons.In Sec.III A we will check the weakly nonlinear analysis against exact analytical solutions with zero surface tension Bϭ0.In Secs.III B 1 and III B 2 we will do the same with numerical simulations of the full problem with finite surface tension B 0 and finite angular velocity ⍀.In Sec.III B 1 we will consider the presence of injection as a destabilizing force while the centrifugal force will be stabilizing.In Sec.III B 2 injection will not be present and the centrifugal force will be destabilizing.

A. Zero surface tension
In Ref. ͓33͔ we showed that the weakly nonlinear analysis at low orders, for the channel geometry, did approximate accurately the zero surface tension evolution of single-finger configurations.However, a word of caution is necessary in general when testing approximations on the zero surface tension case which, being an ill-posed problem ͓8͔, may exhibit different types of pathologies.The most apparent ones are the generation of finite-time singularities, which are regularized by surface tension, but other more subtle singular effects of surface tension have been recently unveiled ͓15,16͔, which show that the integrable dynamics of the zero surface tension problem may be dramatically different from the regularized one even for smooth nonsingular solutions.
These series of results point out precisely to the necessity of having analytical tools that are nonperturbative in surface tension, such as the present weakly nonlinear analysis.A detailed study of the small surface tension limit is obviously beyond the scope of the present paper but we will illustrate how such pathologies of the integrable dynamics show up in the systematic nonlinear analysis, in an example.
Consider the exact solutions of high viscosity contrast (AϭϪ1) discussed in Ref. ͓29͔ and first obtained in Ref.
͓30͔, based on the conformal mapping formalism.The solutions are written in the form where f (,t)ϭzϭxϩiy is an analytical function inside the unit disk in the complex plane, which maps the disk onto the viscous fluid domain in the physical plane.The interface is obtained in a parametric form by setting ϭe i ͑with 0 рϽ2).Real functions a 0 (t) and a n (t) satisfy

͑3.2͒
Constants of motion, k n and k 0 , are fixed by the initial conditions.The injection (QϾ0) destabilizes the interface and the rotation (⍀ ˜) is stabilizing.
It is important to recall the basic properties of this family of solutions.The interface described by Eq. ͑3.1͒ presents nfold symmetry and its Fourier expansion contains only harmonics of the basic n mode.If a n /a 0 is small, the amplitude of mode kϭn satisfies ␦ kϭn Ӎa n and the harmonics are hierarchically ordered as ␦ kϭ2n ϳ␦ n 2 , ␦ kϭ3n ϳ␦ n 3 , . . . .Furthermore, these solutions can present finite-time singularities depending on the initial parameters ͑see Ref. ͓29͔ for de-tails͒.
We proceed with the case of basic periodicity nϭ3.Integrating the linearized dynamics, we obtain a nonexponential growth of mode kϭ3:

͑3.3͒
Using and Eq.͑3.2͒, Eq. ͑3.3͒ reads Next, we compare this linear analysis with the exact evolution of the basic mode kϭ3, which is given by the Fourier transform where we omit the time dependence for simplicity.Using that ϭg(), where g() can be computed from Eq. ͑3.1͒, and integrating by parts, we can write which yields The surprising result here is that the full exact nonlinear evolution of amplitude ␦ 3 of the basic mode kϭ3 coincides exactly with the linearized evolution at all times.This means that the rest of the series cancels out exactly at all times.This remarkable dynamical symmetry does reflect the awkward character of the integrable zero surface tension dynamics.This property is obviously missed by the weakly nonlinear analysis.Notice, however, that the exact knowledge of a specific mode amplitude does not yield the best possible description of the full interface, so even in this case, through the description of the harmonics of the dominant mode, the weakly nonlinear analysis may be useful to approximate the interface evolution on the early stages.In neighborhood of a cusp singularity it will fail again because all orders are necessary then.
We have also checked our scheme with a less pathological exact solution with zero surface tension, namely, the petallike configuration of Ref. ͓37͔, which has injection but no rotation, and is thus the direct analogue of the single-finger configuration of Ref. ͓33͔.In this case, the weakly nonlinear analysis typically performs as well as in the channel geometry.

B. Nonzero surface tension
For the numerical integration of the interface we use the numerical method introduced and described in detail by Hou et al. ͓38͔, as also in Refs.͓14,39,40͔.The interface is parametrized at equally spaced points by means of an equalarclength variable ␣.As a consequence, since s measures arclength along the interface, quantity s ␣ (␣,t) is independent of ␣ and a function of time only.The interface is described using tangent angle (␣,t) and interface length L(t).These dynamical variables replace cartesian coordinates x, y of the interface.The evolution equations are written in terms of (␣,t) and L(t) in such a way that the high-order terms dominant at small scales appear separated from the other terms and in a linear manner.These terms are responsible for the numerical stiffness of the equations, which is significantly reduced through the use of this method.Once the high-order terms appear linearly ͑and with constant coeffi-cients͒ it is straightforward to apply an implicit time integration method to these terms.We have used a linear propagator method of second order in time.We use spectrally accurate spatial discretization.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-wave-number mode becomes larger than the filter level, the number of modes is increased and the amplitude of the additional modes is set to zero.In a typical calculation, 768 discretization points are initially used, and this number is a multiple of three to account for the threefold symmetry of the computed interfaces.Time step ⌬t is decreased until an additional decrease of ⌬t does not change the solution to plotting accuracy, and none of the other physical quantities are significantly different.In all the following numerical solutions we take the initial radius as the length scale, and set R 0 ϭ1 in dimensionless units.We also set the dimensionless time with the characteristic velocity at tϭ0, V 0 ϭV(tϭ0)ϭ1.We define the dimensionless surface tension at tϭ0 as B 0 ϭB(tϭ0).Dimensionless injection and rotation are defined by Q*ϭQ ˜R0 /V 0 and ⍀*ϭ⍀ ˜R0 /V 0 .We will consider values of these parameters which make mode kϭ3 the most unstable at tϭ0 in the linear dispersion relation.All the evolutions considered will have the symmetry ␦ m ϭ␦ Ϫm ϭ␤ ͉m͉ /2, where ␤ ͉m͉ is the amplitude of the cosine function cos(mx) and mϭ3p, where p is an integer.

Configuration with rotation and injection
We now address a configuration with rotation and injection, with no counterpart in channel geometry, which is well described at all times with a weakly nonlinear approximation.We would like to discuss to what extent a few orders can account for the whole interface evolution and hence, a certain truncation of expansion ͑2.8͒ can be regarded as a good model of the system.We will deal with interfaces initially unstable due to injection (QϾ0, AϽ0) and with rotation present as stabilizing force (⌬Ͻ0).As the interface evolves and mean radius R(t) increases, the stabilizing effect of the angular velocity also increases, producing a circular interface at long times.In this configuration the convergence condition is fulfilled at both the initial and final stages.Depending on the values of Q ˜, ⍀ ˜, and ˜and on the initial condition, the intermediate stages of the interface evolution may or may not fulfill the convergence condition.For the case when the convergence criterion is not met during a certain time window, we ask ourselves how well the weakly nonlinear analysis approaches the interface in this intermediate nonconvergent regime.
We use the following parameters: ⍀*ϭ0.001,Q* ϭ0.999, AϭϪ1, and B 0 ϭ1/30.The time window, where the convergence condition is not fulfilled, depends on these parameters and also on the initial condition.First, we consider an initial condition given by mapping f ͑ ,0͒ϭ ͩ 1.0ϩ 0.01 3 ͪ .
Several snapshots of the interfaces defined by conditions ͑3.9͒ and ͑3.10͒ are presented in Figs. 2 and 3. We use dimensionless time ϭV o /͓R o (kϭ3,tϭ0)͔(R o /⌬), where ⌬ is the maximum perturbation of the initial condition.The interface configurations that do not satisfy the convergence criterion of Sec.III are indicated in the figures with dotted lines.
The evolution in Fig. 2 does not fulfill the convergence condition in a very short period of time, while the evolution in Fig. 3 has a longer window of nonconvergence.The basic configuration of fingers is also clearly different, depending on the initial condition.In the first case, the basic pattern consists of three fingers that grow and later on, vanish.Mode kϭ3 is not the unique relevant mode in this configuration.Although mode kϭ6 decays at tϭ0, it starts growing when a larger R(t) is reached because the effective surface tension is reduced.Its amplitude becomes even larger that mode k ϭ3 after a certain period of time.For the evolution shown in Fig. 3, the contribution of mode kϭ6 becomes much larger than that of kϭ3 much earlier, making the basic pattern to be a six-finger configuration with a significant kϭ9 contribution.
Having defined these two configurations, we compare the different linear and weakly nonlinear approximations with the exact evolution.We begin with configuration Eq. ͑3.9͒ and plot ␤ 3 /R(t), ␤ 6 /R(t) in Figs. 4 and 5, respectively.The time window where the convergence condition is not met is indicated with vertical lines.
Mode kϭ3 is approximated fairly well by the linear evolution ͓␤ ˙3ϭ(3)␤ 3 and ␤ ˙6ϭ(6)␤ 6 ], even in the nonconvergence regime.On the contrary, mode kϭ6 is not well described, not only in the nonconvergence regime but in the initial stages as well.A hierarchy of modes in the initial condition, which results in a different order of magnitude for each mode amplitude, such as the hierarchy ␤ 3m ϳ␤ 3 m considered here, makes the order of magnitude of the terms coupling n modes not to correspond with order ␤ n ͑equivalent to order nϪ1 ) of the weakly nonlinear expansion.For example, when ␤ 6 ϳ␤ 3 2 and ␤ 9 ϳ␤ 3 3 , a three-mode coupling ␤ 3 ␤ 3 ␤ 3 has a magnitude similar to that of a two-mode coupling ␤ 3 ␤ 6 or a single mode ␤ 9 .Therefore, linear equation ␤ ˙6ϭ(6)␤ 6 is not systematic because it does not take into account all orders up to ␤ 3 2 in Eq. ͑2.33͒.The proper equations for the initial stages of modes kϭ3 and kϭ6, which are referred to as ''two-mode hierarchy'' in the figures, read with C(6,3)ϭ1/2͓F(6,3)ϪS( 6,3)ϩ(3)J(6,3)͔.Figure 5 shows that these equations correctly approximate mode k ϭ6 in the linear stages and improve the result also in the nonlinear regime.
On the other hand, the weakly nonlinear approximation, Eq. ͑2.33͒, with all the couplings involving ␤ 3 and ␤ 6 , reads 6)J(3,6)͔.This system of equations is not systematic when ␤ 3m ϳ␤ 3 m because term ␤ 3 ␤ 6 of order ␤ 3 3 is considered, while three-mode coupling ␤ 3 ␤ 3 ␤ 3 is not.Nevertheless, this is a partial resummation ͑referred to as ''weakly nonlinear'' in Figs. 4 and 5͒ which improves significantly the evolution of mode kϭ3 while leaving mode kϭ6 almost unchanged with respect to the evolution obtained with Eq. ͑3.11͒.It should be mentioned that a partial resummation does not always lead to a similar improvement, as we will see in the following section.
We now consider initial condition ͑3.10͒.The main plots in Figs. 6 and 7 present the evolution of modes kϭ3 and k ϭ6, respectively, up to the nonconvergence regime.The whole evolution of the modes is plotted in the insets.In this case, linear equations ␤ ˙3ϭ(3)␤ 3 and ␤ ˙6ϭ(6)␤ 6 provide the consistent description of the initial stages, both for modes kϭ3 and kϭ6, as shown in the figures.Since ␤ 3 ϳ␤ 6 , the first proper correction to the linear evolution is the two-mode coupling ͑3.12͒ and not the ''two-mode hierarchy'' ͑3.11͒.Equation ͑3.12͒ is not a partial resummation now since it includes all the terms of order ␤ 3 2 ͓41͔.As shown in Figs. 6  and 7, it improves the linear evolution of modes kϭ3 and kϭ6.
Besides the analysis indicated above, we have also investigated the resummation scheme introduced in Ref. ͓33͔, where (k)␤ k is substituted by ␤ ˙k in the left hand side of Eq. ͑3.12͒.In the present case this scheme does not produce any significant improvement in the approximation during the convergence regime, and the differences are extremely small beyond that point.
From the previous results and from the comparison between the two different initial conditions, we can extract two main practical conclusions for the use of the weakly nonlinear analysis.First, for a given initial condition we can provide specific rules on the appropriate truncation scheme.If the order of magnitude of the initial modes is equal (␤ m ϳ␤ for any m integer͒, each order of the weakly nonlinear expansion is self-consistent regarding the order of magnitude of the couplings involved, i.e. an l-mode coupling is of order ␤ l ͑or lϪ1 ) and therefore, each additional order of the expansion improves the approximation.In practice, one has to deal with a reduced number of modes.Typically, a good criterion is to include only the linearly unstable modes.The stable modes will only play a role when activated through nonlinear couplings at later stages.It is worth remembering that when there is a finite injection rate Q 0, a stable mode at tϭ0 does not necessarily remain stable during all the linear evolution.When the amplitudes of modes are not uniform in the initial condition, the weakly nonlinear equations must be changed into a set of equations consistent with the corresponding hierarchy of amplitudes.To obtain these equations all the modes and couplings which give a contribution of order ␤ p r must be taken into account, regardless of the number of modes involved in the coupling, i.e., regardless of the order of the expansion.As an example, the weakly nonlinear evolution of mode kϭ9, considering initial condition Eq. ͑3.9͒, makes the approximation worse unless all the terms of order ␤ 3 3 are included in the approximation.The second main practical conclusion is that remarkably, a few modes to low orders do approximate the exact solution accurately, up to the nonconvergence regime.Particularly, when the convergence condition is fulfilled along the whole evolution, a low-order weakly nonlinear approximation provides a satisfactory reduced description of the full nonlocal equations.In any case, when there is a small time window of nonconvergence, the approximation may still be remarkably good.

Configuration with rotation as the only destabilizing driving
We compare, finally, the weakly nonlinear approximation with the exact numerical solution of the problem where rotation is the only destabilizing force and injection is not present.Our main purpose is to obtain analytical information, using the weakly nonlinear analysis, about the interplay of rotation and viscosity contrast.This is an intrinsic nonlinear effect, which already shows up in the first weakly nonlinear correction of the linear dispersion relation ͓see Eq. ͑2.33͔͒.The viscosities of the fluids have no influence in the initial mechanism of the instability ͑while densities must fulfill 1 Ͻ 2 ).The lack of injection makes the linear evolution a simple exponential growth ͑or decay͒ independent of viscosity contrast A. We will study the two limiting cases: highest viscosity contrast, Aϭ1, and lowest viscosity contrast, Aϭ0.We also have ⍀*ϭ1 by our definition of parameters and take B 0 ϭ1/45, which makes kϭ3 the most unstable mode during the whole linear evolution.
The numerical simulations of the exact evolution for viscosity contrasts Aϭ0 and Aϭ1 are presented in Fig. 8.The initial condition is set by Eq. ͑3.10͒.The interfaces that do not verify the convergence condition are indicated again with dotted lines.The evolution of the interfaces for the two values of A are very similar in the convergence regime.The beginning of the nonconvergence regime happens almost at the same time for both values of A. The viscosity contrast, however, has a stronger effect in the later stages of the evolution when the outward growing fingers have developed.The shape of the fingers and the width of their necks depend strongly on the viscosity contrast.In particular, the width of the finger necks are significantly smaller for Aϭ0.
We present in Figs. 9 and 10 the evolution of mode k ϭ3 and kϭ6, respectively, for both viscosity contrasts and initial condition ͑3.10͒.Following the discussion of the preceding section on the influence of the hierarchy of modes in the nonlinear approximation, the linear dispersion relation describes the initial stages correctly and two-mode coupling equation ͑2.33͒ ͑''weakly nonlinear'' in the figures͒ improves the linear evolution, including only modes kϭ3,6.
Observing Figs. 9 and 10, we see that mode kϭ3 is always the largest mode and its linear approximation is better for Aϭ0 than for Aϭ1.Actually, the linear evolution and the exact Aϭ0 solution are visually indistinguishable in the convergence regime ͑inset of Fig. 9͒.This comes from the fact that the linear evolution does not depend on A and is therefore well suited to approximate an evolution where viscosity contrast Aϭ0 eliminates all the terms depending on A.
The first weakly nonlinear approximation does depend on A: In the two cases Aϭ0 and Aϭ1, this approximation reproduces the exact evolution in the convergence regime for the two modes kϭ3 and kϭ6 ͑see insets͒.As it is expected, the weakly nonlinear approximation of kϭ3 for Aϭ0 further improves the almost exact linear evolution of kϭ3.It is remarkable that the two leading modes, which basically define the whole interface, are obtained accurately in the convergence regime using only the first nonlinear approximation From Eq. ͑2.33͒ we see how A modifies the effect of rotation at the nonlinear level.For low viscosity contrast the stabilizing role of rotation, for instance, becomes less effective.We thus conclude that the coupling of viscosity contrast and rotation in the first order of approximation has the important effect of making low-viscosity-contrast interfaces more unstable.This prediction of the weakly nonlinear analysis is confirmed in Figs. 9 and 10.The amplitude of modes kϭ3 and kϭ6 is always larger in the case Aϭ0.The effect of the first nonlinear coupling is persistent along the whole evolution, even in the nonconvergent regime.
This result has also been confirmed studying the exact evolution and the weakly nonlinear approximation of modes kϭ3,6, for both Aϭ0 and Aϭ1, with the initial condition f ͑ ,0͒ϭ ͩ 1.0ϩ 0.1 3 ͪ ,

͑3.13͒
where the modes follow a hierarchy of form ␤ 3m ϳ␤ 3 m .Furthermore, the analysis also shows that the partial resummation Eq. ͑3.12͒ worsens the approximation obtained with Eq. ͑3.11͒ in the present configuration, in contrast with the improvement obtained in the preceding subsection.

IV. CONCLUSION AND PERSPECTIVES
We have extended to radial geometry the systematic scheme discussed in Ref. ͓33͔, to derive successive orders of mode couplings in the weakly nonlinear regime of the Saffman-Taylor problem.We have found that the nonlinear expansion converges uniformly in the radial geometry whenever r 2 ϩr 2 Ͻ2R 2 is fulfilled at every point of the interface.We have tested the weakly nonlinear approach against exact solutions with zero surface tension and numerical integration of the full problem in several representative situations.The comparison is satisfactory in general, as in the channel geometry.Difficulties in the approach appear only in classes of solutions which exhibit the ill posedness of the zero surface tension problem in the form of finite-time cusp singularities.The small surface tension region is known to be a very delicate limit which can be studied perturbatively until times of order one, well into the highly nonlinear regime ͓16͔.One of the advantages of the present scheme, though, is that it can be used for arbitrarily large surface tension, which is generically relevant to physical situations.In this case, the method provides an accurate and controlled analytical approximation of the dynamics.
The explicit knowledge of the successive orders of the approximation and of the convergence criterion will be particularly useful in morphologically stable arrangements subject to external perturbations.For instance, our scheme is relevant to work out the systematic nonlinear, nonlocal terms in the problem of fluid invasion of random media ͓34,35͔ if it is ever extended to radial geometry.
We have shown that one of the specific features of the radial geometry, with no counterpart in the channel geometry, is the fact that an appropriate combination of the stabilizing effect of rotation and the destabilizing effect of injection can yield situations where the interface is initially unstable and yet, the weakly nonlinear analysis is useful during the whole evolution.The weakly nonlinear analysis gives rise to a reduced description of the whole Saffman-Taylor dynamics in terms of ordinary differential equations, valid when the convergence criterion is fulfilled at all times.
We have also shown that the weakly nonlinear analysis can provide useful and nonperturbative information in the case when rotation is the unique destabilizing force.Using the first weakly nonlinear correction, we have demonstrated that the main nonlinear coupling between A and ⍀ is stabilizing.This shows that the low-viscosity-contrast case ͑e.g., two similar liquids in the cell͒ is more unstable than the high-viscosity-contrast case ͑e.g., one liquid inside and air outside͒.This result has been confirmed numerically using viscosity contrasts Aϭ0 and Aϭ1.

ACKNOWLEDGMENTS
We acknowledge financial support from the Direccio ´n General de Ensen ˜anza Superior ͑Spain͒ under Project No. DGES BXX2000-0638-C02-02. J.O. acknowledges the Generalitat de Catalunya for additional financial support.E.A.L. also acknowledges a grant from the Comissionat per a Universitats i Recerca of the Generalitat de Catalunya.The work was also supported by the European Commission through Project No. HPRN-CT-2002-00312.

APPENDIX: PROOF OF THE CONVERGENCE CONDITION OF THE NONLINEAR EXPANSION
In order to obtain the different orders of the modecoupling equation, we have carried out two power series expansions.First, we have expanded the inverse of the denominator in Eqs.͑2.1͒ and ͑2.2͒, which corresponds to the expansion of the average velocity and second, we have expanded the inverse of the denominator of curvature ͑2.15͒.If both expansions are uniformly convergent then Eq. ͑2.8͒ is also a uniformly convergent series.
Written in terms of the nonscaled variables, the convergence condition for the curvature is r 2 ϩr 2 Ͻ2R 2 , ͑A1͒ at every point of the interface.The convergence condition for the average velocity is for any two points 1 and 2 considered.Our purpose is to demonstrate that Eqs.͑A1͒ and Eqs.͑A2͒ are fully equivalent and thus, Eq. ͑A1͒ provides the necessary and sufficient convergence condition.First, it is clear that Eq. ͑A2͒ reproduces Eq. ͑A1͒ when 1 ϭ 2 .Therefore, it remains to be proved that condition ͑A1͒ implies Eq. ͑A2͒.

͑A3͒
Eq. ͑A1͒ will imply Eq. ͑A2͒ if and only if any extreme of Z with 1 2 is smaller or equal to the maximum of r 2 ϩr 2 for any value of , i.e., Max"Z( 1 , 2 1 )… рMax(r 2 ϩr 2 ).To simplify the notation we use the following definitions: We also write r( 1 )ϵr 1 , r( 2 )ϵr 2 , and take r 1 уr 2 without loss of generality and hence ⌬у0 ͑taking this prescription, the sign of r ¯ and ␦ cannot be fixed a priori͒.Using this notation, Z is written as We differentiate Z with respect to 1 and 2 and set them to zero to obtain information about the extremes: By adding both equations we obtain the first condition that any extreme must fulfill:

FIG. 2 .
FIG. 2. Numerical simulation of the problem with rotation and injection, using initial condition ͑3.9͒.The initial interface is the tiny circle in the center of the image.Snapshots in interval t/ϭ0 to t/ϭ15.6, separated by ⌬t/ϭ1.3.See text for more details.

FIG. 5 .
FIG. 5. Comparison of the exact evolution of the scaled amplitude of mode kϭ6, against different approximations, using initial condition ͑3.9͒ with both injection and rotation.See text for details.FIG. 6.Comparison of the exact evolution of the scaled amplitude of mode kϭ3, against different approximations, using initial condition ͑3.10͒ with both rotation and injection.The main graph displays the initial convergence regime and the inset the whole evolution.See text for details.

FIG. 8 .
FIG. 8. Numerical simulation of the problem with rotation as the only destabilizing force.The upper evolution corresponds to Aϭ0 and the lower to Aϭ1.The initial condition in both cases is Eq.͑3.10͒.Snapshots in interval t/ϭ0 to t/ϭ0.93, separated by ⌬t/ϭ0.133.See text for more details.

FIG. 9 .
FIG. 9. Comparison of the exact evolution of the scaled amplitude of mode kϭ3, against different approximations, using initial condition ͑3.10͒ with rotation as the only destabilizing force.The end of the convergence regime is indicated by the vertical line.The inset is an enlargement of the main graph in the weakly nonlinear regime.See text for more details.
Introducing this condition in Z we obtain the value of any maxima of Z: To complete the demonstration we need to show that, if r ¯ and ␦ have opposite signs, Max(Z)рMax(r 2 ϩr 2 ) too.We begin considering case ͉␦͉у͉r ¯͉.The following inequalities hold:͑A12͒and use this to write Max(Z) as Max͑Z ͒ϭr ¯2ϩr ¯ 2 cos 4 ͩ 1 Ϫ 2 2 ͪͫ 1Ϫ ͩ ␦ r ¯ ͪ 2 ͬ 2 .