Systematic weakly nonlinear analysis of interfacial instabilities in Hele-Shaw flows

We develop a systematic method to derive all orders of mode couplings in a weakly nonlinear approach to the dynamics of the interface between two immiscible viscous fluids in a Hele-Shaw cell. The method is completely general. It includes both the channel geometry driven by gravity and pressure, and the radial geometry with arbitrary injection and centrifugal driving. We find the finite radius of convergence of the mode-coupling expansion. In the channel geometry we carry out the calculation up to third-order couplings, which is necessary to account for the time-dependent Saffman-Taylor finger solution and the case of zero viscosity contrast. Both in the channel and the radial geometries, the explicit results provide relevant analytical information about the role that the viscosity contrast and the surface tension play in the dynamics. We finally check the quantitative validity of different orders of approximation against a physically relevant, exact time-dependent solution. The agreement between the low order approximations and the exact solution is excellent within the radius of convergence, and reasonably good even beyond that.


I. INTRODUCTION
The morphological instability of fluid interfaces in Hele-Shaw flows [1,2] has become a paradigm of interfacial pattern formation in nonequilibrium systems [3][4][5]. As opposed to the most commonly studied "bulk" pattern forming systems [6], the inherent difficulties of the free-boundary problems associated to interfacial growth makes the latter even more elusive to analytical treatment. As a prototype of interfacial instabilities in diffusion-limited growth problems (including for instance dendritic growth, solidification of mixtures, chemical electrodeposition, flame propagation, etc.) the Saffman-Taylor problem [7] is a relatively simple case, both theoretically and experimentally, well suited to gain insight into generic dynamical features in the broad context of nonlinear interface phenomena.
The intrinsic difficulty of free-boundary problems is manifest in the fact that the interface dynamics is highly non-local. Furthermore, the nature of the instability (except in some cases, such as in directional solidification of binary alloys) usually produces non-saturated growth, which inevitably results in a highly nonlinear dynamics. In bulk instabilities, when the control parameter is near threshold, the traditional weakly nonlinear techniques lead to a universal description of patterns in terms of amplitude equations, based on center manifold reduction [6,8]. These techniques, however, are not so useful for interfacial problems in which nonlinearities do not saturate the growth. This is the case for instance of viscous fingering. In the case of the channel geometry (Saffman-Taylor problem) [1,7] the interface restabilizes in a nontrivial morphology (the Saffman-Taylor finger) which keeps growing at a finite rate. For circular geometries [9,10] the patterns do not reach an equivalent steady state and the interplay of tip-splitting events and screening effects may result in a variety of complicated morphologies. In these problems all the weakly nonlinear techniques apply only to a transient in the early nonlinear regime. Nevertheless, compared to the more traditional ones [6,8], the weakly nonlinear analysis developed in this paper is not restricted to situations near the instability threshold, where a separation of scales is exploited. Instead, we expand on the amplitudes of the whole spectrum of modes.
In this paper we will deal both with channel and radial geometry Hele-Shaw flows. In the traditional Saffman-Taylor problem (channel geometry) the pressure-and gravity-driven instabilities can be formally mapped into each other in the appropriate reference frames, so there is really no different interface dynamics for the two physical situations. The problem, then, contains two independent dimensionless parameters, namely, a dimensionless surface tension B and the viscosity contrast or Atwood ratio A [11]. In the radial geometry, though, there is no such formal mapping. We will see that the injection and the centrifugal forcing are not equivalent and three independent parameters must be considered.
The situation most commonly studied in the literature is the high viscosity contrast limit, A = 1, where one of the two fluids is non-viscous [1] (typically air displacing a viscous fluid). The singular perturbation character of the surface tension B has received most of the attention as responsible for the subtle mechanism of steady-state selection, namely the fact that surface tension "selects" a single finger solution out of a continuum of solutions for B = 0 [12][13][14]. More recently, the crucial role of surface tension in the dynamics of fingering patterns has been pointed out. Siegel and Tanveer [15][16][17] have shown that the exact, nonsingular time-dependent solutions known for the case with B = 0 may differ significantly from the corresponding ones with B → 0 + after a time which is of order one (B 0 ). In practice, this implies that exact solutions of the problem with B = 0 (including those with no finite time singularities) may lead to completely incorrect asymptotic behaviour as compared to the regularized ones, with B ≪ 1. A careful analysis of these questions may be found in Ref. [18]. Notice however that such analysis restricts to small B, while in many cases (for instance, for fingers emerging naturally from the linear instability, with the characteristic length scale of the linearly most unstable mode) the effective dimensionless surface tension is necessarily B ∼ 1. Understanding the dynamics of finger competition in typical experimental conditions thus requires considering relatively large values of B, for which the perturbative techniques of [15][16][17] fail.
On the other hand, an important role of viscosity contrast A in the dynamics of finger competition has been observed both numerically [11,19] and experimentally [20][21][22]. A careful characterization of the interface evolution has shown that for A = 0 the finger competition process is ineffective, and that the system does not approach the usual single finger Saffman-Taylor attractor [23,24]. Although the nature or existence of other attractors is still an open question, it seems that the basin of attraction of the Saffman-Taylor solution does depend on A, and is particularly sensitive to A in the neighborhood of A ≃ 1. In any case, it is clear that the viscosity contrast plays also a crucial role in the highly nonlinear regime, and that tuning A in its full range is necessary to elucidate some of the important open questions.
Finally, an interesting interplay between B and A in connection with the selection problem is apparent in that, despite the fact that single finger stationary solutions of any width do exist for B = 0 regardless of viscosity contrast A, the only single finger time dependent solution of the A = 1 (B = 0) problem which is also a solution for any viscosity contrast A is the one that fills one half of the channel [25,26], which is precisely the solution selected by surface tension in the limit B → 0. Whether deeper consequences concerning the selection problem can be drawn from this fact is also an interesting open question.
In order to gain analytical insight into these dynamical questions we propose here a systematic weakly nonlinear expansion of the problem of viscous fingering in Hele-Shaw flows, including all traditional setups and the most recent one of rotating flows [27]. The basic motivation is to be able to extract information which is nonperturbative in any of the two basic parameters, which are taken as completely arbitrary. The expansion parameter will be basically the mode amplitudes.
In this paper we will focus on unstably stratified flows, for which the approach is necessarily restricted to the early evolution of the interface. Although some of the nontrivial dynamic effects mentioned above are associated to the highly nonlinear regime, it may be useful to know within a controlled approximation to what extent those or other effects show up already at the early stages of nonlinear mode coupling.
For the stably stratified case the weakly nonlinear analysis is obviously valid for long times since all mode amplitudes decay with time. Although this configuration may seem trivial, this is not the case in some situations, for instance when some external source systematically drives the interface out of its equilibrium state (planar or circular). An example of this is the presence of noise sources, such as the quenched noise associated to a porous medium [28]. In a study of the long-wavelength, low-frequency scaling properties of the interface fluctuations, the knowledge of the lowest order nonlinear terms and their dependence on parameters such as viscosity contrast is crucial. In this context the weakly nonlinear expansion is the starting point of any Renormalization Group analysis of the relevant terms and of the fixed points of the problem. This line of research is clearly beyond the scope of this paper and will not be pursued here.
The basic ideas of weakly nonlinear analysis developed here were first applied to viscous fingering by Miranda and Widom [29,30] to both channel and circular geometry (with fluid injection). The present work is in part an extension of those previous contributions in several directions, and in part a detailed study of selected particular situations to assess the validity and limitations of the approach.
First of all, we provide a fully systematic methodology which may be carried out to arbitrary order. We explicitly calculate the results up to third order to include important situations for which the second order couplings, first discussed in Refs. [29,30], vanish identically, such as for A = 0 or for configurations with up-down symmetry (which includes the physically relevant time dependent single finger solution with width one half). We give the explicit predictions both in real and Fourier space.
In our formulation we also address the case of centrifugal forcing of Hele-Shaw flows [31]. The experimental study of rotating Hele-Shaw flows has revealed a rich variety of new phenomena [27,32,33]. From a theoretical point of view, new classes of exact solutions with B = 0 have been found [34,35]. The role of rotation in the possible suppression of finite time cusp singularities in the absence of surface tension has been discussed in [36]. From an experimental point of view, important differences in pattern morphology and new dynamical effects have been found for low viscosity contrasts [37]. It seems thus important to have this case included in the weakly nonlinear formalism.
In addition, our study extends the earlier results of Miranda and Widom [29,30] with the discussion of the convergence of the weakly nonlinear analysis. We find the explicit exact criterion to assure uniform convergence of the series. Beyond that condition the series is asymptotic and different resummation schemes are also explored. The different orders of approximation, including possible resummations, are carefully compared to exact solutions for the case of a single finger configuration. We find that in some cases the agreement even at relatively low orders is quite remarkable.
The layout of the rest of the paper is as follows: in Section II we introduce the formalism. Section III deals with the derivation of the weakly nonlinear equations and their application to Hele-Shaw flows in channel geometry. The analysis of the radial geometry is carried out in Section IV, and Section V presents a numerical analysis of exact and approximate solutions. The main results and the conclusions are summarized in Section VI.

A. Channel Geometry
Let us first consider the Hele-Shaw problem in the channel geometry. We consider fluid 1 (viscosity µ 1 , density ρ 1 ) below fluid 2 (µ 2 , ρ 2 ) (Fig. 1). Theẑ axis is perpendicular to the cell. A velocity V ∞ is imposed at infinity in thê y direction. Gravity points from fluid 2 to fluid 1. The width of the cell is L, the gap between plates is b, and the surface tension between the fluids is σ.
The equations of motion for the interface and the boundary conditions are well known [7]. Here we will use the formulation of Tryggvason and Aref [11]. We introduce the velocity U = ( u 1 + u 2 )/2 as the mean of the two limiting values of the velocities from both sides of the interface ( u 1 , u 2 ) at a given point. This velocity U can be expressed in terms of the vortex sheet distribution γ at the interface as: where: with s the arclength, κ the curvature, and For the purposes of this work it is convenient to rewrite these equations in terms of the interface height h(x), following [38]. The equations read: where:γ The dependence of h andγ on time is not written explicitly.
To complete the definition of the moving boundary problem the continuity of the normal velocity at the interface is required. This means that the velocity in the y direction, dh/dt, projected along the normal direction, is equal to the normal component of the average velocity of the interface, U ·n: In Section III we will consider the interface in the comoving frame (h = 0) and will look for the equation of evolution of h(x, t).

B. Radial geometry
We proceed to write the corresponding equations for a circular cell rotating with angular velocity Ω. The initial condition is R = R 0 (constant) and we label the outer (inner) fluid as the fluid 1 (2) which have known viscosities µ 1 , µ 2 and densities ρ 1 , ρ 2 (Fig. 2).
In the presence of sinks or sources the velocities u 1 and u 2 which define U include only the solenoidal part of the total velocity field. For this reason, in the presence of injection (Q > 0) or withdrawal (Q < 0) of the inner fluid, Eqs. (8) and (9) 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 the local equations known for this problem [27,31]: and the boundary conditions: After some algebra, we can write an expression for the vorticity as a function of U : with The last step is to derive the equation of the continuity of the normal velocity. Now the projection of the radial velocity along the normal directionn has two contributions: the solenoidal part of the average velocity, U , projected alongn, and the irrotational part of the mean interface velocity, U irrot : This completes the system of equations for the two geometries considered. For an interface with no overhangs, Eqs. (4), (5), and (7) are the starting point for the generalized weakly nonlinear analysis in the rectangular cell. Equations (8), (9), (12), and (14) play the same role for the analysis in the rotating circular cell.

III. SYSTEMATIC WEAKLY NONLINEAR ANALYSIS. CHANNEL GEOMETRY
Our goal in this section is to introduce a systematic method to derive an evolution equation of the interface in real space to a given order in nonlinear couplings, in the channel geometry. 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 form: where F [h], G[h], . . . are nonlocal operators on the function h(x, t), including nonlinearities of order n + 1 in the term of order ε n . The small parameter ε is defined as the ratio of two lengths, ε = w/L. We take w as a measure of the characteristic scale of variation of the interface h(x), while L is either the width of the cell or, alternatively, the characteristic scale of variation in the x direction. The weakly nonlinear regime is defined by the condition w ≪ L.
The order ε 0 in Eq. (15) corresponds to the linearized equation. The order ε 1 , when written in Fourier space, corresponds to the result of Miranda and Widom [29]. Here we will perform the explicit calculation to one order higher in ε (up to I[h]) which is the leading nonlinear contribution in several important cases, such as the one discussed in Section III D. The systematics is presented in the channel geometry for simplicity.

A. Dimensionless equations, expansion and convergence
We scale the interface height with w, the coordinate x with L, and the time with L/C, where the velocity C has been defined in (3).
The equations (4), (5), (6), and (7) become: where and The starting point of our approach is an expansion of U ,γ and κ in powers of ε. Notice that the term 2 between curly brackets in (16) is bounded provided that h x does not diverge. Equation (19) thus takes the form: Before pursuing the calculation in detail, let us point out that the ε-expansion in Eq.(20) has a finite radius of convergence. This is guaranteed by the properties of uniform convergence of both the expansion of the denominator in Eq.(16) and the curvature. These properties allow us to commute the expansion with the integral in Eq.(16) and yield a convergent series of the form (20). The radius of convergence of the expansion of the denominator of Eq.(16) is given by the condition while the convergence of the curvature expansion is assured by the condition ε 2 h 2 x < 1. In the nonscaled original variables, the two conditions for convergence coincide, and read M AX(|h x |) < 1. If |h x | < 1 in the whole domain of integration, then the ε-expansion converges. If the condition is not fulfilled in some intervals, then Eq. (20) is an asymptotic expansion. Even in this case, the expansion contains useful information about the original problem.
An interesting case is A = 0 which makes the vorticity independent of U in Eq. (16), so that Eq.(19) becomes a closed equation for h(x, t). Then, in Eq. (20) it is easy to show that U (n) y = 0 when n is even and U (n) x = 0 when n is odd, a property that makes the even power terms of the expansion in Eq. (20) to vanish.
B. First and second order expansion (ε 0 and ε 1 ) Following the scheme introduced in the previous section we obtain from Eq.(16) that: With this result we can study the first order term in the vorticity equation:γ Taking into account the definition of the Hilbert Transform: Eq. (20) up to order ε 0 reads: The linear operator F [h] in Eq.(15) thus reads explicitly: Writing h(x, t) as a superposition of Fourier modes in Eq.(25) we recover the linear dispersion relation: We will take k as an integer but we should keep in mind that, upon restoring dimensions, k should become (2π/L)n, with n integer. Let us pursue the systematics of the method by computing the next order in ε in Eq. (19): The computation of U x , which must be computed from the first order term of the vorticity. We get: Taking the Fourier transform, we obtain: which coincides with the result of Miranda and Widom in Ref. [29].
C. Third order expansion (ε 2 ) The expansion to order ε 2 is necessary to account for the lowest order nonlinearities in the case of zero viscosity contrast, and for other relevant situations such as the time dependent Saffman-Taylor finger solutions (section III D). We now have: We have already computed U in Eq. (29). On the other hand: Integrating twice by parts, the last integral can also be written as a Hilbert Transform. After some algebra we obtain the explicit form of the operator I[h] containing the cubic nonlinearities in Eq. (15), which reads: with: and: The same result in Fourier space takes the form: with: It is clear that the viscosity contrast in Eq.(37) appears squared because of the reflection symmetry (the simultaneous change A → −A and h → −h is a dynamical symmetry of the problem). Symmetry reasons alone, however, do not allow to discard a three mode coupling contribution when A = 0. We see from our calculation that a three mode coupling is indeed present independently of A.
Following this scheme, the fourth order will carry a contribution proportional to A, and another proportional to A 3 , for symmetry reasons. The fifth order will carry a contribution independent of A and two others, proportional to A 2 and A 4 , respectively. This scheme will continue for subsequent orders. At this point it is worth recalling that the case A = 1 allows for a continuum of Saffman-Taylor finger solutions corresponding to different finger widths. A continuum of time dependent exact solutions, leading to those stationary states, is also known for A = 1. However, for A = 1 only the one of width λ = 1 2 remains a solution, λ being the ratio of the finger width to the width of the channel. This result, which has been recently addressed in Ref. [26] although it was first discovered in Ref. [25], points out an intriguing connection between the width selection problem and the dynamical role of viscosity contrast. Here we will analyze the interplay between A and λ in the early nonlinear regime, and elucidate at what stage of the nonlinear dynamics does the viscosity contrast A = 1 prevents the possibility of having λ = 1/2.
From now on we consider L = 2π, B = 0, C = 1, and we take as scaling velocity (1 − λ)C = (1 − λ) ≡ η 2 . Conformal mapping techniques enable us to write the single finger solution of this problem in the form [25,26]: where f (w, t) = y + ix, is an analytic function inside the unit disk in the w-complex plane, which maps that disk into the physical region occupied by the more viscous fluid. The interface is obtained in a parametric form by setting w = e iθ . The functions α(t), d(t) verify:ḋ To obtain an expression for h(x) in the weakly nonlinear regime of the evolution, α(t) and d(t) are expanded in powers of a small parameter ν: Introducing these expansions in Eqs.(42) and (43) we obtain: which will be useful later. From (41) and (42) we obtain: This last equation can be inverted in a systematic way to get: The relation between θ and x can be inverted in Eq.(47) and the cosine functions expanded. We obtain in this way the following expression for h(x, t): In order to follow the scheme developed in the previous section, we must measure h(x, t) in units of its characteristic amplitude ν. In this way, the previous equation (50) takes the form: where h(x, t) represents now h(x, t)/ν, and the small parameter ν is directly comparable to the small parameter ε of the previous section. The expression for h(x, t) can be regarded as a series of modes of decreasing amplitude in our expansion (15) (written in the units of this problem), and matching the corresponding powers in either ε or ν we obtain: to order ε 0 . This identity can be verified (independently of A) by substituting the expression of h (0) and using Eqs.(46) and (25) for the left and right hand side respectively. At order ε we have: The second term of the right hand side gives no contribution, and we obtain an interesting result: the equation at order ε is verified independently of η and A. Hence, we cannot establish the difference between A = 1 (compatible with any η) and A = 1 (compatible only with η = 1) at this order. The difference between these two situations arises at order ε 2 : Since the left hand side of the equation involves only the modes cos x, cos 2x, and cos 3x, the right hand side of the equation must include these same modes with the same coefficients. The coefficients for cos 2x and cos 3x are easily matched. For cos x the left hand side reads: where we have used Eq. (46), and the right hand side reads: Hence, the requirement for matching the coefficients on both sides is: which is always true if η = 1 (λ = 1/2), but for η = 1 (λ = 1/2) it requires A = 1. This result shows that the nontrivial relationship between A and λ which is known from exact solutions is already manifest at the early nonlinear stages of the dynamics. This clearly illustrates the potential usefulness of the weakly nonlinear expansion at a purely analytical level, in that a dynamical property of the problem which must be satisfied at all orders (in this case the incompatibility of A = 1 and λ = 1/2) may be detected at a low order in the expansion, without having to know the exact solution to all orders. If we pursue the expansion to higher orders, the general expression for h (n) takes the form: where each coefficient is a function of α(t). The even modes of the solution with η = 1 have zero coefficients because they must remain invariant under a change of sign and translation by π. So far in this section we have restricted ourselves to B = 0 in order to compare with exact results. However, we can carry out the analysis also for B = 0. It can be shown that the structure of the expression (58) is preserved in this case. The coefficients β (n) k (t) are obtained as solutions of a set of differential equations which contain the surface tension B. For instance, to third order we have: which, introduced in Eq. (15), leads to the following equations: once the dimensions are reintroduced. In this way we can also see how surface tension perturbs the dynamics at the weakly nonlinear regime. Clearly, at these early stages of the nonlinear evolution there is no sign of the singular perturbation character of surface tension, which, at the late stages of the evolution will be responsible for the selection of the steady state [1].

IV. APPLICATION OF THE METHOD TO THE ROTATING HELE-SHAW CELL
A. Dimensionless equations in the radial geometry The ε expansion in the channel geometry originated in the different scaling of lengths in the x and y directions. To perform an equivalent expansion in the radial geometry we first split the function r(φ, t) in two parts, namely the radius of the unperturbed circle and the departure from this circle: The largest length scale in the problem is given by R(t), which is only a constant in the absence of injection. For finite injection rate, R(t) defines an evolving (unstable) solution. The proper scaling of r(φ, t), which will define the weakly nonlinear dynamics now, is not the naive extension of the scaling defined in the channel geometry (namely ε = w/R, w being the typical departure from circularity). In fact, this scaling would mix different orders of the weakly nonlinear analysis. The reason for this was already pointed out by Miranda and Widom [30]: 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: Following this observation, we split the deviation from R in two terms, so that the radial function reads r = R+r+r 0 , where r 0 represents the zero mode (an expansion or contraction of the circle), andr the other modes. We scale them in the formr → wr, r 0 → (w 2 /R)r 0 , and write the position of the interface r(φ) as: where: To simplify the notation, we have dropped out the time dependence. The velocities will be scaled by: implying that time will be scaled by R/V 0 . Once equations (8), (9), (12), and (14) are made dimensionless, we have: for the vorticity, where: For the integrals: with: and for the evolution of the deviation from R(t): Eq.(71) above is the only one which is sensitive to the different possible scalings of the deviation from circularity, due to the time derivatives. The choice we propose is the only one which is truly systematic. Other possibilities are not incorrect but will mix different orders of the weakly nonlinear expansion at a given order in ε.

B. The linear dispersion relation
Our goal now is to obtain the linear dispersion relation and the leading weakly nonlinear corrections, in both real and Fourier space, and to discuss the interplay between rotation and injection at the early stages of the instability.
We will use the following definitions for the average of a function and for the Hilbert transform in the unit circle: The zero order of the vorticity and the two components of the velocity are: and thus:γ This is different from its counterpart in the channel geometry. 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 linearized equation for the evolution of the interface deviation reads: with:γ We can now perform a Fourier transform of the equations and, after reintroducing the adequate dimensions, we find: This is the linear dispersion relation found in [27]. We now return to the question of the different scaling alternatives introduced at the end of Section IV A. Had we used the scaling ε = w/R 2 (t), the term Q/(2πR 2 ) in (77) 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 (63) and (64) becomes δ n /R in the first of the scaling alternatives above.
The term Q/(2πR 2 ) at the end of Eq. (77) 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 [36] the term Q/(2πR 2 ), which can be scaled out, must be distinguished from the rest, which do describe the intrinsic instability of the problem. Disregarding that last term, it is then clear that injection and rotation (CV 0 and DV 0 respectively) play an equivalent role in the linear dispersion relation, since both appear multiplying n. 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 all the nonlinear evolution (in the appropriate reference frame). This is not the case in radial geometry, as shown in the next section.

C. Weakly nonlinear theory with rotation
The purpose of this section is to derive the leading order nonlinear contributions for the general problem with both rotation and injection. We first recall that mass conservation related the zero mode to 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 U can be obtained from the expansion of Ur: where:γ κ (2) and U (2) φ can be computed from (67) and (69) respectively. The latter involves onlyγ (1) . We obtain finally: with: By introducing a superposition of Fourier modes inr we directly recover Eq.(78). For the other modes (k = 0) we get the result:δ where: For the particular case Ω = 0, the expression above reproduces the result of Miranda and Widom [30]. We find that the presence of rotation does not change the formal structure of the equations, but introduces new terms.
To emphasize the roles of injection and rotation, we define the quantity H(k, n) as the part of the coupling matrix in the r.h.s. of Eq.(84) which contains Q and/or Ω. This quantity reads: 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Ω andQ is different in the linear part and in the leading nonlinear part of the dynamics. This clearly shows that the effects of injection and rotation cannot be interchanged by simple changes in parameters, and that the intrinsic new dimensionless parameter related to rotation already shows up at the early nonlinear regime. It is also remarkable that, for Q = 0, R(t) is time dependent, and the relative weight ofΩ andQ changes with time. This may have remarkable consequences in the highly nonlinear regime. For instance, rotation can prevent the formation of finite time singularities in the case with zero surface tension [36].
Finally, notice that at nonlinear order there is also a geometrical term, Q/(2πR 3 ), independent of n. This term cannot be scaled out by the same procedure that eliminated its counterpart at linear order.

V. NUMERICAL STUDY OF AN EXACT SOLUTION AND ITS WEAKLY NONLINEAR APPROXIMATION
The purpose of this section is to put the weakly nonlinear analysis to the test, as a quantitative approximation. We will check it against an exact time dependent solution of the case B = 0 which evolves towards the Saffman-Taylor finger of width λ = 1/2 in the channel geometry. This will allow us to check how fast is the convergence of the expansion and how accurate it may be even beyond its radius of convergence, when it becomes an asymptotic series.
We define F as the ratio between the maximum height of the interface (at the finger tip) and half the width of the cell (in our case L = 2π). F is a dimensionless amplitude which is proportional to ε and thus measures how deep the system is into the nonlinear regime. Furhermore, since the system is described by the evolution of a curve, it may also be convenient to have a more global characterization of its configuration in order to compare the exact result and the different approximations. We propose the use of the flux Φ, defined as the total amount of fluid 1 per unit length and unit time flowing across the horizontal line located at the mean interface position [18]). An explicit solution of the problem without surface tension (B = 0) and valid for any viscosity contrast A, describing the growth of a single finger which asymptotically fills half of the channel (λ = 1/2), can be obtained from Eqs.(41), (42), and (43). This reads: with a(0) and b(0) as initial constants, and b(0) ≫ 1. Since this solution is valid for any A, particularly A = 0, all terms of the weakly nonlinear equation (15) depending on A will necessarily have no contribution. We will take the following values for the parameters: (91) Fig. 3 shows the evolution of the interface for this solution. When F ≃ 1 the finger shape is indistinguishable from the stationary Saffman-Taylor solution (Fig. 4), except near the points x = π/2 and x = 3π/2 when the asymptotes of the finger must develop at infinite time. When F = 1/π, the slope of the interface becomes 1 at these two points and the series (15) loses its convergence, according to the discussion of Section III A.
In Fig. 5 we present the amplitude of the Fourier modes of the exact solution. It is clear that relatively few modes suffice to accurately reproduce the stationary finger shape in the tip region.

B. Expansion of the exact solution
In Section III D we expanded the exact solution (41) in a small parameter ν. This expansion revealed a hierarchy of amplitudes of the different modes for an initial condition close to the planar interface. The solution (90) can be expanded similarly and yields: In agreement with (50), and using the same parameters as in the previous section, we obtain up to third order: h(x, t) = νe t/2 cos x + ν 3 1 2 e t/2 − 1 8 e 3t/2 cos x − 1 24 e 3t/2 cos 3x .
The exact flux Φ can be easily computed from the stream function, which in its turn can be related to the form of the mapping Eq.(90) (see for instance Ref. [18]). In Fig. 6 we compare the flux Φ of the approximate result above with the exact result coming from (90). We clearly observe that the first correction to the linear regime improves the flux (to a 5 % accuracy) from F ≃ 0.12 to F ≃ 0.17. The amplitude of the mode k = 1 (Fig. 8) is fairly well reproduced by the linear approximation up to F ≃ 0.20, and the correction of order ν 3 gives a value for the mode k = 3 reasonably good up to F ≃ 0.20 (Fig. 9) and improves the mode k = 1 almost until F ≃ 0.25. Above these values of F the deviation from the exact interface is exponential. Notice that the expansion up to ν 3 considered here does not account for the complete three mode coupling of Eq.(37) which contains part of the higher orders in ν. This explains why going from order ν (linear) to order ν 3 improves only slightly the shape of the finger, the flux Φ, and the amplitude of the mode k = 1.

C. Weakly nonlinear expansion
In this section we study the mode coupling equation (37). We solve this equation in cases where only modes k = 1 and k = 1, 3 are present, and compare the result with both the exact solution and the approximation to order ν 3 introduced in the two previous sections.
We consider the initial condition β 1 (0) = ε = 0.001 for the first mode and zero for the rest of modes. We use Eq. (37) recalling that, since β i represent the amplitudes of the cosine functions, δ i = β i /2. The result is that the modes k = 1 remain zero and: where we see that the mode k = 1 saturates to a finite value as t → ∞. Returning to Fig. 6 and 7, we see that the flux Φ approximates the exact value to a 5 % uncertainty) until F ≃ 0.23. Fig. 8 shows that the amplitude of the first mode starts deviating significantly from the exact solution around F ≃ 0.40. Next, we add the third mode to the equation (37) with the initial condition β 1 = ε + 3ε 3 /8 and β 3 = ε 3 /24. The set of equations to be solved numerically is: As shown in Figs. 6 and 9 the flux does not improve significantly, while the amplitude of the mode k = 3 remains acceptable almost until F ≃ 0.25, and improves the previous result by more than a 50 % almost until F ≃ 0.33. Beyond this value of F the third mode falls to zero due to the second term in the r.h.s. of Eq.(96), proportional to β 2 1 β 3 but negative. This is a signature that higher order terms, for instance of order β 4 1 β 3 , become important when β 1 is of order one.
From this study we can conclude that using only two modes and the lowest order nonlinear correction for this case (involving three mode couplings), both the shape of the interface and the total flux are reasonably well described within the whole range of convergence of the series (up to F ≃ 1/π). This indicates that the weakly nonlinear description of the problem works remarkably well at the quantitative level, at least in some physically relevant cases.
To what extent this conclusion holds for more complicated situations, for instance those involving competition of different fingers, remains an open question.

D. Partial resummation
According to Eq. (37), we can write a closed equation for the mode k = 1 which is correct up to third order couplings, of the form:β where the mode k = 3 has been set initially to zero. The presence ofβ 1 in the right hand side gives rise to terms of order higher than β 3 1 . Thus, solving this equation amounts to doing a partial resummation of these higher order terms. This kind of closed differential equation was already obtained in Ref. [29]. In our approach this differential equation is obtained in a systematic way by identifying the terms λ(n)β n and replacing them forβ n . For example, in Eq. (37) the term T (k, s, l) could have been partially resummated in the two mode coupling differential equation.
The solution of Eq.(97) yields the transcendental form: A numerical computation of the flux based on this equation shows a dramatic improvement up to F ≃ 1/π. It is surprising and remarkable that the equation for a single mode, including the resummation suggested by the first nonlinear correction, reproduces the exact solution to a great accuracy in the full range of convergence and defines an excellent approximation further inside the nonlinear regime. This suggests that this kind of resummation is not arbitrary, and might have a deeper physical justification. In addition, we obtain an amplitude for the mode k = 1 within 10 % accuracy up to F ≃ 0.8, and it is still quite good until F = 1. Although the exact spectrum shows the increasing importance of the mode k = 3 at these values of F , the agreement between the exact result and this approximation for a single mode k = 1 is quite remarkable.
We recall that the series converges up to F ≃ 1/π. Beyond this value of F , increasing the order of mode coupling in Eq. (97) does not guarantee that the approximation improves. Accordingly, there is an optimal finite order of approximation for a given value of ε, as usual in asymptotic series. In this sense, our results above seem to indicate that Eq.(97) is the optimal approximation for the mode k = 1 in the asymptotic (nonconvergent) region (i.e., for moderately long times).
In the same spirit of Eq. (97), we can also derive a closed system of differential equations which includes third order couplings for the modes k = 1, 3. Its solution shows no qualitative differences with the case of the single mode k = 1 above, as far as the flux and the amplitude of the first mode are concerned. The improvement in the amplitude of the third mode and consequently in the shape of the finger does not reach the high values of F that the first mode reaches (F ≃ 0.8), as we can see in Fig. 9. In a short range of F the mode k = 3 is better than the case shown in the previous section, but deviates from the exact solution much earlier than the mode k = 1 for the same reasons than in the previous section. Specifically the mode k = 3 is extremely good until F ≃ 0.25, starts deviating by more than a 10 % at F ≃ 0.33 and is about 50 % of the value at F ≃ 0.40.
We conclude that the weakly nonlinear formalism, apart from its nominal range of validity for very small amplitudes where it converges very fast to the exact dynamics, can describe regimes beyond those expected in principle with reasonable accuracy. In the case of the growth of a single finger, for instance, we have explicitly seen that simply two modes together with an appropriate partial resummation, yields a remarkably good approximation in the full range of convergence of the expansion (F < 1/π).

VI. CONCLUSION AND PERSPECTIVES
We have developed a systematic scheme to derive the successive orders of mode couplings in a weakly nonlinear regime, adapted to the study of interfacial dynamics in Hele-Shaw flows. This has been done with full generality, including both the channel geometry driven by gravity and pressure, and the radial geometry with arbitrary injection and centrifugal driving, which includes the case of a rotating Hele-Shaw cell. The method could also be applied to even more general (position dependent) drivings. The formulation in real space has enabled us to address the issue of convergence of the mode coupling expansion. We have found that the exact condition for convergence in the channel geometry is |h x | < 1 in every point of the interface.
We have carried out the explicit derivation of nonlinear couplings up to third order in the case of channel geometry, in order to obtain the leading nonlinear contributions in cases where the second order contribution vanishes. These include the case of zero viscosity contrast and the time dependent single finger solution of width λ = 1/2.
On the analytical side, we have also illustrated the usefulness of the weakly nonlinear analysis in elucidating the role of the different parameters on the dynamics, in some examples. In the case of single finger solutions growing in a channel, we have shown how an exact nontrivial property of the problem (the relationship between λ and viscosity contrast A for zero surface tension) can be extracted from an analysis order by order, without really knowing the exact solution. In the case of rotating Hele-Shaw flows, we have discussed the asymmetry between the roles of injection and rotation at the lowest nonlinear order, as opposed to what happens in the channel geometry, where the equivalence between gravity and injection driving is exact.
On the numerical side we have checked the predictions of the weakly nonlinear analysis, at different orders, against exact solutions. We have found that the convergence is quite fast for small amplitudes. Furthermore, in the case of a single finger growing in a channel we have explicitly seen that the analysis to low orders (up to three mode couplings) yields a good description of the interface dynamics even close to the radius of convergence. With an appropriate resummation scheme, the prediction is relatively good even beyond that point.
Some of the most interesting applications of the systematic weakly nonlinear approach have not been developed here since they deserve a separate, in-depth analysis. One of them is the study of the scaling properties of fluctuations in stably stratified Hele-Shaw flows with external noise [28]. A more direct application of our formalism which also deserves a separate analysis is the derivation of amplitude equations for the region near the instability threshold through a center manifold reduction. This would made contact with what is most commonly referred to as "weakly nonlinear" approach in the literature of pattern formation. A detailed study of this point with a careful discussion of the nature of the bifurcation and its implications will be presented elsewhere [39].

VII. ACKNOWLEDGEMENTS
We acknowledge financial support from the Dirección General de Enseñanza Superior (Spain) under Projects PB96-1001-C02-02 and PB97-0906, E. Alvarez-Lacalle 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 Project ERB FMRX-CT96-0085 (Training and Mobility of Researchers).