Jean-Luc Achard
email: jean-luc.achard@hmg.inpg.fr

Sanda-Carmen Georgescu

Quasi-steady formation of bubbles and drops viewed as processes that break bifurcation

Keywords: bifurcation curve, bubbles and drops formation, cusp, meniscus stability, Rayleigh-Taylor

Instabilities occurring during the quasi-steady formation of bubbles and drops at a circular orifice through a thin plate, which separates a cylindrical upper vessel of quiescent liquid from a lower air chamber maintained at constant pressure, are considered in connection with the Rayleigh-Taylor instability of a flat meniscus. The control parameters of the problem are the E ötv ös number E ö, and the excess-pressure number , which characterizes the pressure difference between gas and liquid across the orifice. The Rayleigh-Taylor case i.e., = 0, can be viewed as a perfect bifurcating problem. A subcritical bifurcation emerges from the critical point, E öc = 5•783186, beyond which the flat meniscus is unstable for axisymmetric perturbations. Bubbles ( > 0), and drops ( < 0) appear as solutions that break bifurcation. When an appropriate measure of their magnitude ε is introduced, it can be shown analytically that the equilibrium surface at the critical point is a cusp; its intermediate sheet is stable, while its two upper and lower sheets are unstable. The analytical bifurcation set onto the control parameters plane is valid only around the critical point.

Introduction

In this paper, we study the well-known problem of the Rayleigh-Taylor instability in connection with the more general problems regarding the formation of emerging bubbles and pendant drops. The Rayleigh-Taylor instability refers to the instability of a plane interface between two immiscible fluids of different density in an external gravitational field. While examining the stability of drops and bubbles, at least in their quasi-steady formation regime, we are led to consider interfaces that are not necessarily flat, but upward (bubbles), or downward (drops) curved. The purpose of this paper is to present a unified theoretical treatment of those two types of instability in a well-defined physical system.

Rayleigh-Taylor instability has generated a large theoretical literature [START_REF] Sharp | An overview of Rayleigh-Taylor instability[END_REF], tracing back to Maxwell [2, pp. 56-71], long before Rayleigh [START_REF] Lord Rayleigh | Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density[END_REF] and Taylor [START_REF] Taylor | The instability of liquid surfaces when accelerated in a direction perpendicular to their planes[END_REF]. As noted by Jacobs and Catton [START_REF] Jacobs | Three-dimensional Rayleigh-Taylor instability. Part 2. Experiment[END_REF], there is a comparable amount of work devoted to the experimental validation of the theory, beginning with the pioneering studies of Duprez [6, Plate 1], [7, Plate 1], and Lewis [START_REF] Lewis | The instability of liquid surfaces when accelerated in a direction perpendicular to their planes[END_REF]. The formation of emerging bubbles and pendant drops is not less studied in many cases, both theoretically and experimentally [START_REF] Michael | Meniscus stability[END_REF]; summaries of the pertinent literature are presented by Tsuge [START_REF] Tsuge | Hydrodynamic of bubble formation from submerged orifices[END_REF]Chapter 9,, and Sadhal et al. [START_REF] Sadhal | Transport Phenomena with Drops and Bubbles[END_REF]Chapter 7]. Facing the huge amount of the above mentioned studies, it is useful to classify them according to some simple criteria, which will characterize either peculiarities of the physical systems, or assumptions of the theoretical models. This classification will help us to specify the framework of our own study.

Physical systems can be distinguished by the following three main conditions: Control condition -one of the fluids is maintained either at a constant pressure, or at a constant flow; Lateral constraints -the two superposed fluids being contained into a finite apparatus, the interface is attached at some edge (orifice through a thin plate, opening mouth of a nozzle), or it joins the vertical and more or less smooth bounding walls; Geometrical configuration of the bounding walls -it may be two-dimensional (parallel and vertical plates separated by a gap), axisymmetric (vertical circular cylinder), or three-dimensional (rectangular vessel).

We are interested in a meniscus attached to the edge of a circular orifice through a thin plate. This plate separates an upper cylindrical vessel of quiescent liquid, from a lower air chamber maintained at constant pressure.

Regarding the Rayleigh-Taylor instabilities, the theoretical models can be distinguished by at least six basic categories of assumptions that form their basis: (I) The interface is of infinite, or finite lateral extent. Many theoretical works, but not all of them [START_REF] Pimbley | Stationary solutions of the problem of Rayleigh-Taylor instability[END_REF], have adopted the former condition, leaving aside boundary conditions at the interface. On a smooth, rigid, solid surface, the contact line is free to move, and one is faced with the difficult problem of removing the singularity resulting from the usual hydrodynamic assumptions, i.e., Newtonian fluid behaviour, and no slip on the wall [START_REF] Dussan | Stability in systems with moving contact lines[END_REF]. Fortunately, there are simple situations where the contact line can be considered as prescribed at a fixed location (e.g. sharp edge of an orifice, wire rim); (II) Upper and lower fluids are infinitely thick, or have a finite thickness [START_REF] Elgowainy | The Rayleigh-Taylor instability of viscous fluid layers[END_REF]; (III) Fluids may be inviscid and irrotational, or viscous effects are introduced [START_REF] Bellman | Effects of surface tension and viscosity on Taylor instability[END_REF]; (IV) The density of the lighter fluid is negligibly small, or a finite density ratio is allowed [START_REF] Chandrasekhar | The character of the equilibrium of an incompressible heavy viscous fluid of variable density[END_REF]; (V) Axisymmetric, two-dimensional, or fully three-dimensional initial disturbances are considered [START_REF] Jacobs | Three-dimensional Rayleigh-Taylor instability. Part 1. Weakly nonlinear theory[END_REF]. Oblong apparatuses promote the second case; (VI) The growth rates of instabilities are examined according to linear or nonlinear methods. The former are analytic, being restricted to describe the initial stage of the process, while the latter are generally numerical [START_REF] Daly | Numerical study of two-fluid Rayleigh-Taylor instability[END_REF], describing the subsequent stages (spikes, bubbles, drops,. . . ). Nevertheless, some analytical studies [START_REF] Pimbley | Stationary solutions of the problem of Rayleigh-Taylor instability[END_REF][START_REF] Tan | On the steady solutions of the problem of Rayleigh-Taylor instability[END_REF] are devoted to these stages.

In any of these six categories of assumptions, the first option corresponds to the simplest, and to the more manageable assumption, while the second option corresponds to a more realistic, or more general model. Each time, one among the first authors who tried to perform this extension has been quoted.

In our study, we will adopt systematically the simplest option, except for the first category of assumptions where, as stated above, the meniscus is attached to the edge of an orifice. The fixed end-point condition applies on it. As we shall see, there are some loose restrictions on the contact angle at this point.

To some extent our model is modest. Our objective is not to refine the analysis of the Rayleigh-Taylor instability, which was already thoroughly studied, but to show that, if we enlarge the set of control parameters for this phenomenon, we address simultaneously other phenomena as the slow (quasi-steady) formation of bubbles and drops. This enlargement means basically to take into account initial and unavoidable pressure steps, which exist at the interface. In order to provide a new global picture of previously separated instabilities, a simple model is sufficient.

Formulation

Consider a physical system composed of a cylindrical vessel of liquid over a gas chamber. The vessel is assumed open at the upper side. A thin plate perforated with a single small orifice of radius R separates both fluids. At the orifice, a meniscus forms (Figure 1).

The meniscus may be flat (Rayleigh-Taylor case), upward oriented (emerging bubble), or downward oriented (pendant drop). The gas density is very small with respect to the liquid density ρ L , so that the gas pressure p G is uniform, while the vertically oriented gravity g imposes a hydrostatic profile to the liquid pressure p L . The cylindrical vessel containing the liquid is assumed large enough, so that the effect of the sidewalls is negligible. Also, the depth of the liquid is so big that the free surface is ignored. Mathematically, the liquid above the plate will be assumed to extend to infinity. The steady initial interface shape, as well as its time of evolution are determined by R, ρ L , p G , g, by the hydrostatic pressure on the plate p L0 , and also by the surface tension σ . It is possible to identify two independent dimensionless parameters, which are supposed to control this evolution. The E ötv ös number, which is also known as the Bond number [START_REF] Sadhal | Transport Phenomena with Drops and Bubbles[END_REF], measures the importance of the gravitational force compared to the surface tension force:

E ö = gρ L R 2 /σ.
(

The second control parameter is the dimensionless excess pressure:

= (p G -p L0 )R/σ, (2)
which characterizes the pressure difference between gas and liquid across the orifice. The Rayleigh-Taylor case, where the meniscus is initially flat, is one situation given by = 0; by analogy to the bending of an initially straight column, this case will be called the perfect problem.

Bubble and drop formation is defined by > 0, and < 0, respectively. They correspond to initially bent column problems, and will be termed imperfect problems.

In some literature devoted to equilibrium shapes of drops and bubbles [START_REF] Chesters | An analytical solution for the profile and volume of a small drop or bubble symmetrical about a vertical axis[END_REF], it is quite common to introduce the radius of curvature at the apex of the bubble or drop, instead of the orifice radius in Equation [START_REF] Sharp | An overview of Rayleigh-Taylor instability[END_REF]. Then the differential equations giving the sought profiles are further modified, by shifting the origin to this apex. These two operations give equations containing a single resulting parameter, the shape factor. As mentioned by Rienstra [START_REF] Rienstra | The shape of a sessile drop for small and large surface tension[END_REF], the simplification is only apparent since the distance of the apex to the plate is a part of the solution. Furthermore, the shape factor is not convenient to use when the profile is flat, and the radius of curvature is infinite.

In addition to the set of two parameters {E ö, } controlling the evolution of the interface, we might have introduced some Reynolds number expressing the influence of liquid viscosity. Neglecting such a refinement is justified, provided that the Ohnesorge number is Oh = µ L / √ 2ρ L σ R 1, where µ L is the liquid viscosity. We choose a cylindrical polar coordinate system (r, θ, z), with r = 0 on the orifice axis, and z = 0, the horizontal reference plane, at the plate level. Gravity is acting in the negative zdirection. Throughout this paper we restrict our attention to axisymmetric configurations with respect to the z-axis. For bubble and drop configurations close to the flat meniscus, interfaces have a simple (one-to-one) projection on the reference plane, so that a representation z = ξ (r, t) is convenient.

The appropriate boundary conditions that govern the continuously changing contact line, which at time t = 0 is attached to the orifice edge, will be simplified due to the following common hypotheses [START_REF] Marmur | Equilibrium shapes and quasi-static formation of bubbles at submerged orifice[END_REF]. First, as the orifice edge is rounded, the contact line can find a location on this edge, respecting the local physics of contact angles. Basically, this physics stipulates that, for one contact location, a region of equilibrium contact angles θ eq exists [START_REF] Dussan | On the spreading of liquids on solid surfaces: static and dynamic contact lines[END_REF] with: θ r < θ eq < θ a , where θ r and θ a are the limits of the receding and advancing contact angles. Otherwise a complex dynamic process occurs through which the contact line looks for another place of contact, where the above inequalities can be satisfied.

The interface configurations that will be considered in this study are those having a contact line that does not leave the round part of the edge. The corresponding permissible limits can be expressed in terms of the azimuthal angle ψ c at the contact point. When a bubble emerges over the plate, it can occur that no equilibrium contact angle can be found on the edge as soon as ψ c decreases to θ r : the contact line is then on the point of retreating. On the contrary, when a drop hangs under the plate, it can occur that no equilibrium contact angle can be found on the edge as soon as (ψ cπ ) increases to θ a : the contact line is then on the point of advancing.

Secondly, we assume that the radius of curvature of the orifice edge is much smaller than the radius of the orifice. Thus, adopting a macroscopic point of view, the contact line can be considered as attached to the seemingly sharp-edged orifice defined by: ξ (R, t) = 0. All the detailed local physics of contact angles can be by-passed, provided that the studied interface configurations are such that θ r < ψ c < (π +θ a ).

A small disturbance is applied to the system, and we wish to study the subsequent motion of the interface, and of the liquid. The imposed axial symmetry of the boundary conditions favours mainly axial disturbances to which our analysis will be restricted. Now that we have a better definition of our model, we can list the questions that our study will address: (I) What is the stability of the null solution (perfect problem)? This is not more than the classical Rayleigh-Taylor problem, which will be given for sake of completeness. (II) What are the steady solutions that bifurcate from the null solution? More specifically, we shall give the values of E ö (or branch-points) from which at least a bifurcation emanates; a theoretical description of bifurcations (or branches) in the neighbourhood of the branch-points, i.e., number, dependence of the bifurcations on E ö. (III) What is the stability of the bifurcating solution? (IV) What are isolated solutions that break bifurcation? They correspond, of course, to emerging bubbles and pendant drops, and we shall give a theoretical description of these solutions in terms of in the neighbourhood of the branch-points. (V) What is the stability of isolated solutions?

Governing equations

The motion in the incompressible inviscid liquid above the interface is assumed to be irrotational initially and remain so, with the velocity potential φ, and velocity u = ∇φ. The velocity potential φ satisfies the Laplace equation:

∇ 2 φ = 0, for    0 ≥ r < ∞ ξ (r, t) ≥ z < ∞, if r ≤ R 0 ≥ z < ∞, if r > R , ( 3
)
and its evolution is governed by the Bernoulli equation:

ρ L ∂φ ∂t + 1 2 |∇φ| 2 + p L + ρ L gz = C (t) , for    0 ≤ r < ∞ ξ (r, t) ≥ z < ∞, if r ≤ R 0 ≤ z < ∞, if r > R . ( 4
)
The velocity potential must be finite along the z-axis, and it vanishes at infinity for the considered axisymmetric problem:

∂φ/∂r = 0, for (r = 0, ξ (0, t) ≤ z < ∞) , ( 5
)
φ = 0, for (r → ∞, 0 ≥ z < ∞) . ( 6
)
The pressure there becomes the steady hydrostatic pressure. Thus, on the plate it is:

p L (r, 0, t) = C (t) = p L0 . ( 7
)
On the interface, where 1 + (∂ξ/∂r)2 1/2 n = (ξ r , -1), the so-called kinematic condition applies:

∇φ • n = - ∂ξ ∂t 1 + ∂ξ ∂r 2 -1/2 , or ∂φ ∂z = ∂ξ ∂t + ∂φ ∂r ∂ξ ∂r , for 0 ≤ r < R z = ξ (r, t) , ( 8
)
and the boundary condition on the impermeable plate is given by

∂φ/∂z = 0, for (R ≤ r < ∞, z = 0) . ( 9
)
Here the momentum balance on the interface reduces to its normal component (the additional pressure jump due to surface tension is the Laplace-Young boundary condition):

p G -p L = 2σ H, for (0 ≤ r ≤ R, z = ξ (r, t)) , ( 10
)
where H is the mean curvature given by:

H = - 1 2 ∇ • n = - r∂ 2 ξ/∂r 2 + 1 + (∂ξ/∂r) 2 ∂ξ/∂r 2r 1 + (∂ξ/∂r) 2 3/2 . ( 11
)
Two additional boundary conditions for the interface equation must be specified:

∂ξ/∂r = 0 at (r = 0) , and ξ = 0 at (r = R) . [START_REF] Pimbley | Stationary solutions of the problem of Rayleigh-Taylor instability[END_REF] Initial conditions are also prescribed:

ξ (r, 0) = ξ i (r) , and ∂ξ ∂t (r, 0) = ξi (r) , for (0 ≤ r ≤ R, t = 0) . ( 13
)
To get a better idea of the mathematical structure of the above problem, it would be interesting to express it solely in terms of ξ (r, t). One way to do this, is first to solve the potential problem defined by Equations ( 3), ( 5), ( 6), [START_REF] Lewis | The instability of liquid surfaces when accelerated in a direction perpendicular to their planes[END_REF], and ( 9), showing explicitly the dependence of φ on the various parts of the data. Assuming that the Green function G ξ exists for this partial problem [24, Volume II, pp. 130-145], and initial conditions for φ being put equal to 0, we may show that φ is a highly nonlinear functional of ξ r , t :

φ (r, z, t) = -2π R 0 r G ξ r, z r , ξ r , t ∂ξ ∂t r , t 1 + ∂ξ/∂r 2 -1/2 dr . ( 14
)
Combining momentum balance and Bernoulli equation evaluated at the interface gives the following equation, written in dimensionless form:

-E ö ∂φ * ∂t * + |∇ * φ * | 2 2 = ∂ 2 ξ * ∂r * 2 1 + ∂ξ * ∂r * 2 -3 2 + 1 r * ∂ξ * ∂r * 1 + ∂ξ * ∂r *
where the length scale R, and velocity scale U = √ gR have been adopted. The dimensionless variables have been taken as:

ξ * = ξ/R, r * = r/R, φ * = φ/ (UR) , and t * = Ut/R. ( 16
)
Substituting the dimensionless form of Equation ( 14) for φ * evaluated at z * = ξ * (r * , t * ) in Equation ( 15), yields a complex nonlinear integro-differential evolution equation governing ξ * (r * , t * ). The left-hand side contains the transient part. It is well-defined only if the Green function G ξ is determined explicitly. In fact, such a task cannot be performed analytically, excepted if we limit ourselves to the study of small disturbances of the interface from a steady flat profile. This will be done in the next section, where the Rayleigh-Taylor stability problem will be addressed. The right-hand side contains the steady part; we shall denote this part by F (E ö, ξ * , ). We may observe that the flat profile ξ * = 0 (steady solutions are noted with a tilde) is the solution of:

F E ö, ξ * , 0 = 0. ( 17
)

Stability of the null solution (Rayleigh-Taylor stability problem)

The particular static interface, for which the stability is tested in this section, is the flat interface ξ * = 0; it is given by = 0. We linearize the general dynamic problem stated above by considering small perturbations ξ * , φ * about a given static equilibrium configuration:

ξ * r * , t * = ξ * r * + ξ * r * , t * , ( 18
)
φ * r * , z * , t * = φ * r * , z * , t * . ( 19
)
Following the procedure defined at the end of Section 3, a small disturbance ξ * about ξ * is shown to satisfy the linear integro-differential equation:

2π E ö 1 0 r * G 0 r * , 0 r * , 0 ∂ 2 ξ * ∂t * 2 r * , t * dr * = F ξ E ö, 0, 0| ξ * , ( 20
)
where we have defined the linearized steady operator:

F ξ E ö, 0, 0| ξ * ≡ ∂ 2 ξ * ∂r * 2 + 1 r * ∂ ξ * ∂r * + E ö ξ * , (21)
which is the Bessel operator of order 0. The boundary and initial conditions ( 12) and ( 13) become:

∂ ξ * /∂r * = 0, and ξ * = ∞, for r * = 0, t * > 0 , ( 22
)
ξ * = 0, for r * = 1, t * > 0 , ( 23
)
ξ * = ξ * i , and ∂ ξ * /∂t * = ξ * i , for 0 ≤ r * ≤ 1, t * = 0 , ( 24
)
For convenience asterisks will be suppressed in the remainder of this paper. The Green function G 0 is defined by the following expression (see Appendix A):

G 0 r, z| r , z = 1 4π ∞ 0 J 0 (γ r) J 0 γ r e -|z-z |γ + e -|z+z |γ dγ, (25)
where J 0 (γ r) is the Bessel function of first kind and order zero, with parameter γ 2 .

Using an appropriate eigenfunction expansion solves the above stability problem. The corresponding eigenvalue problem [24, Volume I, p. 282] is based on the zero-order Bessel equation. Assuming a uniformly converging series, the solution of Equation ( 20) can be expanded in a complete set, {Z k (r)}, with coefficients depending parametrically on t. Thus,

ξ (r, t) = ∞ k=1 ξ (k) (t) Z k (r), where ξ (k) (t) = 1 0 r ξ (r, t) Z k (r) dr = ξ, Z k . ( 26
)
The set of functions:

Z k (r) = √ 2J 0 (λ k r)/J 1 (λ k ) ( 27
)
is orthonormal over the range 0 < r < 1. Here λ k is the kth zero of J 0 . In our case, the eigenvalue problem is self-adjoint. As a consequence, {Z k } is an orthonormal system of eigenfunctions. Assuming uniform convergence, initial values of ξ (k) (t) result from [START_REF] Stakgold | Boundary Value Problems of Mathematical Physics. I[END_REF]:

ξ (k) (0) = ξi , Z k , and ∂ξ (k) (0) ∂t = ξ i , Z k . ( 28
)
We multiply both sides of the integro-differential Equation ( 20) by (rZ k (r)) , and integrate over the range [0, 1] to obtain:

2πE ö ∞ j =1 d 2 ξ (j ) (t) dt 2 1 0 r G 0 (r, 0|r , 0), Z k (r) Z j (r )dr = 1 0 ∂ ∂r r ∂ξ ∂r Z k dr + E ö ξ, Z k . ( 29
) Since 1 0 ∂ ∂r r ∂ξ ∂r Z k dr = 1 0 ξ ∂ ∂r r ∂Z k ∂r dr = -λ 2 k 1 0 rξ Z k dr = -λ 2 k ξ (k) (t) , ( 30
)
and using the definition of Bessel functions, we finally have the differential equation:

E ö 2 [J 1 (λ k )] 2 d 2 ξ (k) (t) dt 2 + λ 2 k -E ö ξ (k) (t) = 0, (31)
subject to the initial conditions [START_REF] Padday | Transient surface properties in liquid bridge and pendant drop menisci in gravity and low gravity[END_REF]. The coefficient of d 2 ξ (k) /dt 2 is viewed as an added mass. For small motions, the upper fluid imposes on the kth interface component some inertia, which can be characterized by a single mass coefficient

m k = (E ö[J 1 (λ k )] 2 )/2.
Since the upper fluid has been assumed inviscid, there is no friction term. The coefficient of ξ (k) has the role of a restoring force.

The most unstable component is obtained for k = 1. Thus, the general condition for loss of stability must be derived from the equation controlling this component:

m 1 d 2 ξ (1) dt 2 -µξ (1) = 0, (32)
where µ = E öλ 2 1 = E ö -5•783186, which corresponds to the first eigenvalue of the axisymmetric problem; λ 1 = 2•40482556. Let us define ξ (1) = e σ (1) t x (1) as a typical solution of the Equation (32). The corresponding eigenvalues σ (1) have to satisfy:

m 1 σ (1) 2
x (1) µx (1) = 0. (33) Thus, there are two eigenvalues:

Im( ) Re(s) O s + (2) (m) s (1) (m) s + (1) (m) s (1) (m) - - - (1) (m)
σ (1) ± (µ) = ± µ/m 1 . ( 34
)
They are imaginary when µ < 0, and real when µ ≥ 0. In the latter case, there is always a positive eigenvalue, and we conclude that ξ = 0 is unstable. If µ ≥ 0 implies instability, stability may be obtained when µ < 0, i.e., when √ µ is purely imaginary. In fact, in this case, an extra negative real term should be present to dampen oscillations. This term would obviously come from viscous effects, which have been neglected. At the critical value µ = 0, σ

± = 0 is a double real eigenvalue. We observe that φ and ξ (1) = ξ i , Z 1 t + ξi , Z 1 are given via the Green function as in Equation ( 14). The behaviour of σ [START_REF] Sharp | An overview of Rayleigh-Taylor instability[END_REF] ± (µ), when µ increases around 0, is presented in Figure 2. When σ [START_REF] Sharp | An overview of Rayleigh-Taylor instability[END_REF] ± coalesce at 0, other eigenvalues σ (k) ± are still on the imaginary axis. They follow the same path as σ (1) ± when µ keeps increasing.

Bifurcating solutions

There is a close connection between loss of stability of the linearized problem (i.e., between some critical values of the control parameters), and the condition for existence of a branch-point [START_REF] Stakgold | Branching of solutions of nonlinear equations[END_REF].

So we look for steady bifurcating solutions in the neighbourhood of E ö = λ 2 1 = 5•783186 (or µ = 0), ξ = 0, and = 0. Let us expand the nonlinear operator F(λ 2 1 + µ, ξ, ) around the steady solution ξ = 0, which satisfies Equation (17):

F λ 2 1 + µ, ξ, = F ξ λ 2 1 + µ, 0, | ξ + 1 2 F ξξ λ 2 1 + µ, 0, |ξ |ξ + + 1 6 F ξξξ λ 2 1 + µ, 0, |ξ |ξ |ξ + O |ξ | 4 = 0. ( 35
)
To simplify the notation throughout the paper, we shall use f for F, and for all partial derivatives of F with respect to the indicated subscript, when these quantities are evaluated at the point (λ 2 1 , 0, 0). We obtain the following relations: = 0. Steady solutions that bifurcate will be constructed as a power series in the solution amplitude ε, that is:

f = F (λ 2 1 , 0, 0) = 1, f ξ (u) = F ξ (λ 2 1 , 0, 0|u) = d 2 u dr 2 + 1 r du dr + λ 2 1 u, f ξµ (u) = F ξµ (λ 2 1 , 0, 0|u) = u, f ξ (u) = F ξ (λ 2 1 , 0, 0|u) = 0, f ξξ (u|v) = F ξξ (λ 2 1 , 0, 0|u|v) = 0, f ξξµ (u|v) = F ξξµ (λ 2 1 , 0, 0|u|v) = 0, f ξξ (u|v) = F ξξ (λ 2 1 , 0, 0|u|v) = 0, f ξξξ (u|v|w) = F ξξξ (λ 2 1 , 0, 0|u|v|w) = -3 d 2 u
ξ(ε) = ε ξ1 + ε 2 2 ξ2 + ε 3 6 ξ3 + O |ε| 4 , ( 37
)
µ(ε) = εµ 1 + ε 2 2 µ 2 + ε 3 6 µ 3 + O |ε| 4 . ( 38
)
The amplitude

ε = ξ(r), Z 1 = 1 0 r ξ(r)Z 1 (r)dr = √ 2 1 0 r ξ J 0 (λ 1 r) J 1 (λ 1 ) dr ( 39
)
is defined by a projection on the eigensubspace associated with the adjoint eigenvector Z 1 (r), belonging to the eigenvalue µ = 0 of the problem:

f ξ ( ξ) + µf ξµ ( ξ) = 0. ( 40
)
The amplitude definition may be regarded as a normalization, which implies some constraints on ξk appearing in Equation (37):

ξ1 , Z 1 = 1 and ξk , Z 1 = 0 for k > 1. ( 41
)
Physically, ε can be regarded as a scalar, which characterizes the bubble (ε > 0), or the drop (ε < 0) magnitude. We may try to obtain ξk , and µ k (k ≥ 1), by expanding the left-hand side of F(λ 2 1 + µ(ε), ξ(ε), 0) = 0 in powers of ε, and identifying independent powers of ε, that is:

f ξ ( ξ1 ) = 0, ( 42
)
f ξ ( ξ2 ) + 2µ 1 f ξµ ( ξ1 ) = 0, ( 43
)
f ξ ( ξ3 ) + 3µ 1 f ξµ ( ξ2 ) + 3µ 2 f ξµ ( ξ1 ) + f ξξξ ξ1 | ξ1 | ξ1 = 0, ( 44
)
where relations (36) have been taken into account. Since ξ1 , Z 1 = 1, the eigenvalue problem (42) has a single solution, namely

ξ1 = Z 1 = √ 2J 0 (λ 1 r)/J 1 (λ 1 ). ( 45
)
The Fredholm alternative applies for the solvability of the other problems ( 43) and (44). Consider first Equation (43). Since µ 1 f ξµ ( ξ1 ), Z 1 = µ 1 Z 1 , Z 1 = 0, and Z 1 , Z 1 = 1, we have µ 1 = 0. The resulting equation f ξ ( ξ2 ) = 0 may then be solved. It has also a single solution: ξ2 = 0, since ξ2 , Z 1 = 0. For the resulting problem (44), the compatibility condition gives:

µ 2 = - 1 3
f ξξξ (Z 1 |Z 1 |Z 1 ) , Z 1 f ξµ (Z 1 ), Z 1 = -73• 4685 (46)
To obtain this result we have used:

f ξξξ (Z 1 |Z 1 |Z 1 ), Z 1 = 12λ 3 1 J 4 1 (λ 1 ) -2 1 0 J 0 (λ 1 r)J 3 1 (λ 1 r)dr + 3λ 1 1 0 rJ 2 0 (λ 1 r)J 2 1 (λ 1 r)dr = 220• 4055, ( 47
)
while the denominator is

f ξµ (Z 1 ), Z 1 = Z 1 , Z 1 = 1.
The associated equation becomes with ξ3 , Z 1 = 0:

f ξ ( ξ3 ) = -3µ 2 Z 1 -9Z 1 dZ 1 dr 2 - 6 r dZ 1 dr 3 . ( 48
)
So, two distinct branches emanate symmetrically from E ö = λ 2 1 (or µ = 0):

µ = 1 2 µ 2 ε 2 + O |ε| 4 = -36.73425ε 2 + O |ε| 4 . ( 49
)
It is easy to show that only even powers of ε can occur in Equation ( 49). Conversely, concerning the ξ -expansion, only odd powers of ε occur, showing that bubble profiles and drop profiles are symmetrical with respect to the orifice. As µ 2 is negative, the branching is to the left; as this branching is to the stable side, the system is subcritical. Obviously it is quite difficult to compute further terms in Equation (48) analytically.

A bifurcation diagram restricted to the neighbourhood of ε = 0, including the analytical results (Equation 49), is given in Figure 3. The upper branch corresponds to bubbles, the lower to drops. The flat meniscus (ε = 0) is stable before the critical point P 1 , where µ < 0. The solution ξ = 0 loses its stability when µ = (E öλ 2 1 ) is increased beyond zero. New solutions µ(ε 2 ), corresponding to curved menisci (bubbles, or drops), undergo double-point bifurcation at the singular point (µ, ε) = (0, 0) denoted by P 1 .

Stability of the bifurcating solution

We study the stability of the steady bifurcating solution we have just found. In Section 4 we showed that, if µ is increased along the branch ε = 0, criticality is reached at µ = 0; a conjugate pair of purely imaginary eigenvalues gives rise to two real eigenvalues, one having a strictly positive part (instability). Instead of increasing µ from µ = 0, let us increase |ε| along the other two branches, and examine the behaviour of this pair. Such bifurcating solutions (µ(ε), ξ(ε)) satisfy:

F(λ 2 1 + µ(ε), ξ(ε), 0) = 0. ( 50
)
Instead of formulating the fully linearized problem for the small perturbations ξ, φ about the above solutions as in Section 4, let us only consider the associated spectral problem for X and F defined by: ξ(ε; r, t) = ξ(ε; r) + ξ(ε; r, t) = ξ(ε; r) + e σ (ε)t X(ε; r), (51)

φ(ε; r, z, t) = φ(ε; r, z, t) = e σ (ε)t F (ε; r). ( 52
)
This spectral (eigenvalue) problem can be decomposed into two parts, one for F, the other for X. The first part for F is a classical Laplace problem with Neumann boundary conditions; it is complicated by the dependence of F on the shape of the mathematical domain, which is modified by the moving boundary z = ξ(ε; r, t). The second part, involving F and X, results from the interfacial momentum balance incorporating the Bernoulli equation evaluated at the moving boundary. Following a procedure similar to the one defined at the end of Section 3, which collects all the preceding equations, we obtain an extended eigenvalue problem for the single function X(r):

σ 2 (ε)K(ε; X(r )) = F ξ λ 2 1 + µ(ε), ξ(ε; r), 0|X(r) , (53)
where

K(ε; X(r )) = 2π(λ 2 1 + µ(ε))G ξ r, ξ(ε; r)|r , ξ(ε; r ) × ×   1 + ∂ ξ(ε; r ) ∂r 2   -1 2 , X(r ) . ( 54
)
The sign of the real part of σ (1) (ε), which is the eigenvalue with the largest real part, determines the stability of ξ(ε; r). It can be shown that there is always a positive eigenvalue in the neighbourhood of ε = 0, i.e., σ (1)

(ε) = ±9•71ε + O |ε| 2 , ( 55
)
so we conclude that ξ(ε; r) is unstable as the solution bifurcates at the critical point.

Isolated solutions that break bifurcation and their stability

We consider here the imperfect problem, i.e., = 0. Using the simplified notations defined in Section 5, we can write the steady equation in the neighbourhood of (λ 2 1 , 0, 0) as

F(λ 2 1 + µ, ξ, ) = f + f ξ ( ξ) + µf ξµ ( ξ) + 1 6 f ξξξ ( ξ | ξ | ξ) + µ 6 f ξξξµ ( ξ | ξ | ξ) + +R(µ, ξ) = 0. (56)
In this expansion we have used the Equations (35), and (36). The remainder R is:

R(µ, ξ) = (λ 2 1 + µ) ∞ n=2 a n ∂ ξ ∂r 2n ξ, (57)
where a 2 = 3/8, and a n = 3(-1) n (2n -5)!/2 n n!, with n ≥ 3. We have seen that f ξ (u) = 0 at (µ, ε) = (0, 0), where F(λ 2 1 + µ, ξ, 0) = 0. Due to the Equation (45), the following key condition for the remainder of this section is guaranteed:

f , Z 1 = √ 2 J 1 (λ 1 ) 1 0 rJ 0 (λ 1 r)dr = √ 2 λ 1 = 0. ( 58
)
The amplitude ε that we have already defined for = 0 in Equation ( 39) is extended here for = 0. In that context of imperfection theory, we follow Iooss and Joseph [START_REF] Iooss | Elementary Stability and Bifurcation Theory[END_REF]Chapter 6,, and invoke condition (58), in conjunction with the implicit-function theorem, claim the existence of a smooth solution ξ(µ, ε), and (µ, ε) of:

F λ 2 1 + µ, ξ (µ, ε), (µ, ε) = 0, ( 59
) ε = ξ(µ, ε), Z 1 , ( 60
)
as well as the values of these functions at ε = 0, i.e., ξ(µ, 0) = 0, and (µ, 0) = 0. We expand ξ = ξ(µ, ε), and = (µ, ε) as series around (0, 0) up to the third order:

ξ = ξµ µ + ξε ε + 1 2 ξµµ µ 2 + 2 ξµε µε + ξεε ε 2 + + 1 6 ξµµµ µ 3 + 3 ξµµε µ 2 ε + 3 ξµεε µε 2 + ξεεε ε 3 + O (|µ| + |ε|) 4 , (61)
= µ µ + ε ε + 1 2 µµ µ 2 + 2 µε µε + εε ε 2 + + 1 6 µµµ µ 3 + 3 µµε µ 2 ε + 3 µεε µε 2 + εεε ε 3 + O (|µ| + |ε|) 4 . ( 62
)
Owing to the above-mentioned values of these functions at ε = 0, we can write:

ξµ (0, 0) = ξµµ (0, 0) = ξµµµ (0, 0) = 0, ( 63
) µ (0, 0) = µµ (0, 0) = µµµ (0, 0) = 0. ( 64
)
Inserting Equations (61), and (62) into (56), and identifying powers of ε i , µ i , ε i µ j , (i ≥ 1, j ≥ 1), we find to first order:

f ξ ξε = -ε , (65)
to second order:

f ξ ξµε = -µε -f ξµ ξε , ( 66
)
f ξ ξεε = -εε , ( 67
)
and to third order:

f ξ ξµµε = -µµε -2f ξµ ξµε , ( 68
)
f ξ ξµεε = -µεε -f ξµ ξεε , ( 69
)
f ξ ξεεε = -εεε -f ξξξ ξε | ξε | ξε . ( 70
)
Again the Fredholm alternative is used to insure the solvability of the Equation (65), and successively of all the other differential Equations (66)-(70). Moreover, the normalization Equation (39) imposes ξε , Z 1 =1, ξµε , Z 1 = ξεε , Z 1 =0, and

ξµµε , Z 1 = ξµεε , Z 1 = ξεεε , Z 1 = 0. ( 71
)
Thus, considering the first-order Equation (65), we obtain:

ξε = Z 1 , and ε = 0. ( 72
)
In the second-order Equation (66), we have ξµε = 0, and the insertion of (72) into the solvability condition gives:

µε = - f ξµ (Z 1 ), Z 1 1, Z 1 = - Z 1 , Z 1 1, Z 1 = - λ 1 √ 2 = -1•70046. ( 73
)
The solution of this equation may be obtained by the eigenfunction method

ξµε = - ∞ k=2 λ 1 J 1 (λ k )Z k λ k (λ 2 k -λ 2 1 )J 1 (λ 1 ) . ( 74
)
In the second-order Equation (67), we observe that ξεε = 0, and εε = 0. Now the solvability conditions of the three third-order equations are simply, on account of the preceding results:

µµε = - 2 f ξµ ( ξµε ), Z 1 1, Z 1 = - 2 ξµε , Z 1 1, Z 1 = 0, ( 75
) µεε = - f ξµ ξεε , Z 1 1, Z 1 = - ξεε , Z 1 1, Z 1 = 0, ( 76
) εεε = - f ξξξ (Z 1 |Z 1 |Z 1 ), Z 1 1, Z 1 = - λ 1 √ 2 f ξξξ (Z 1 |Z 1 |Z 1 ), Z 1 = 12λ 4 1 √ 2J 4 1 (λ 1 ) 2 1 0 J 0 (λ 1 r)J 3 1 (λ 1 r)dr -3λ 1 1 0 rJ 2 0 (λ 1 r)J 2 1 (λ 1 r)dr = -374•79173. (77)
Computing ξµµε , ξµεε , and ξεεε does not present any difficulty, and will not be presented here. What is worth mentioning is the final form of , which reduces to:

(µ, ε) = µε µε + εεε 6 ε 3 + O |ε|(|µ| + |ε|) 3 . ( 78
)
With ( 73) and (77), the Equation (78) becomes:

(µ, ε) = -1•70046µε -62•46528ε 3 . ( 79
)
This analytical steady solution defines the equilibrium 3D surface, which is plotted in Figure 4. In the terminology of Catastrophe Theory we have identified a cusp. The upper cusp sheet corresponds to bubbles, the lower to drops. For some typical µ-values, cross-sections of the equilibrium surface (µ, , ε) are plotted in Figure 5. For µ < 0, there is a stable branch that crosses the µ-axis between the turning point P - 1 for drops, and P + 1 for bubbles (solid line in Figure 5). For ε-values larger than the one corresponding to P + 1 , and for ε-values smaller than the one corresponding to P - 1 , the branches are unstable (dashed lines in Figure 5). For the particular value µ = 0, the whole bubble branch (ε > 0), and the whole drop branch (ε < 0) are unstable, the cross-section passing beyond the critical point P 1 .

The turning points P + 1 and P - 1 extend into > 0, and < 0 regions as singularity lines (fold curves) of the cusp. So, the intermediate cusp sheet is stable, being delimited by those two singularity lines. The upper and lower cusp sheets are unstable. The bubble and drop configurations, corresponding to the intermediate sheet, are the ones that are usually observed. The intersection of the 3D equilibrium surface with the plane = 0 returns the afore-described Rayleigh-Taylor problem.

The stability study of isolated solutions ( = 0), which perturb bifurcating solutions ( = 0), will be based on the stability study of the latter, presented in Section 6. The relations µ = µ(ε), and ξ = ξ(ε) for = 0, will be replaced by the relations = (µ, ε), and ξ = ξ(µ, ε), in which µ is a fixed parameter such that µ < 0.

The bifurcation curve near the critical point P 1 is given by:

= (µ, ε), and 0 = ∂ ∂ε (µ, ε), ( 80
)
where ε plays the role of a curve parameter. This curve is composed of two branches through which the number of steady solutions varies when a turning point, P + 1 for bubbles, or P - 1 for drops, is reached. These two branches merge at the bifurcation point P 1 at µ = 0, ε = 0, and = 0. Each branch is given by the system issued from Equations (79) and (80): The analytical bifurcation curve near the critical point P 1 is expanded in Figure 6 on the control parameter plane (µ, ).

= 124•93056ε 3 µ= -110•20302ε 2 . ( 81)

Conclusions

A systematic theoretical analysis has afforded a unified picture of the Rayleigh-Taylor instability, and of the quasi-steady formation of bubbles and drops, under constant-pressure conditions. Via , the excess pressure number ( > 0 for bubbles, and < 0 for drops), these latter processes are viewed in this paper as perturbing the Rayleigh-Taylor physical situation, which appears as a perfect bifurcating problem with = 0. In this idealized problem, the flat interface is stretched according to a value measured by the E ötv ös number E ö (or µ = E ö -5• 783186). When E ö increases, the stretching decreases to some critical value corresponding to E öc = 5•783186 (or µ = 0). Beyond E öc , the flat interface is unstable, and in its neighbourhood no steady stable solution is available (a transient process necessarily takes place). For a fixed gas/liquid couple, the critical value E öc has a practical meaning, allowing computation of the maximum orifice radius, which sustains a meniscus in equilibrium. On the contrary, before reaching the critical value E öc , at = 0, two sets of steady symmetric deflections ξ(r), one upward (bubble), and one downward (drop), are theoretically available to the interface. These solutions, whose magnitudes can be described by some appropriate scalar ε = ξ, Z 1 , emerge from the critical point P 1 to the stable side. They correspond to a subcritical bifurcation.

An explicit picture of the 3D equilibrium surface near the critical point has been given when the excess pressure number is added to E ö as an extra control parameter. A cusp has been identified. The turning points P + 1 and P - 1 extend into > 0, and < 0 regions as singularity lines (fold curves) of that cusp. The intermediate sheet of the cusp is stable, while its two upper and lower sheets are unstable. The intersection of the 3D equilibrium surface with the plane = 0 returns the afore-described Rayleigh-Taylor problem.

The bubble and drop configurations corresponding to the intermediate cusp sheet are the ones that are usually observed in experiments [START_REF] Hugues | The formation of bubbles at simple orifices[END_REF][START_REF] Padday | Transient surface properties in liquid bridge and pendant drop menisci in gravity and low gravity[END_REF]. They are steady, stable curved menisci and represent the first stage of the interface growing evolution [START_REF] Marmur | Equilibrium shapes and quasi-static formation of bubbles at submerged orifice[END_REF][START_REF] Davidson | Bubble formation at an orifice in a viscous liquid[END_REF]. Within this first stage, at a given E ötv ös number (or µ-value), the gas pressure increases for bubbles from the hydrostatic pressure p L0 (or = 0), to its maximum value where = max . For drops, the gas pressure decreases from p L0 to its minimum value where = min . The gas excess pressure has equal absolute extreme values for bubbles and drops, max = | min |, at the turning points P + 1 and P - 1 (Figure 5). The practical value of the proposed description is limited since all the above results are known from the literature. However, it has two important aspects:

-it has an educational value by allowing a unified treatment of separate classical descriptions of interface instabilities. The role of the nondimensional parameters is clearly identified; -it allows a sound basis for future studies. In fact, it can be shown numerically [START_REF] Georgescu | Évolution d'une Bulle: Formation à Partir d'un Orifice et Éclatement à la Traversée d'une Surface Libre[END_REF]Chapter 3] that below the critical E ötv ös number the upper and lower sheets keep on zigzagging away, and that between the turning points, stable equilibrium states exist, and correspond to bubbles or drops having an undulatory structure. In [START_REF] Georgescu | Multiple equilibrium states during the quasi-steady formation of bubbles and drops at a circular orifice[END_REF] (in preparation), such computations are discussed, and experimental evidence of the physical reality of some mathematical solutions we have exhibited presented.

Figure 1 .

 1 Figure 1. Emerging bubble (left), flat meniscus (centre), and pendant drop (right).

Figure 2 .Figure 3 .

 23 Figure 2. Behaviour of the eigenvalues around the critical value µ = 0.

Figure 4 .

 4 Figure 4. Cusp near the critical point P 1 . Analytical results for µ ∈ [-3•9, 1•2], ε ∈ [-0•32, 0•32].

Figure 5 .Figure 6 .

 56 Figure 5. Cross sections of the equilibrium surface for µ = 0, µ = -1•5, and µ = -3.

-1 2 + + E ö ξ * ,[START_REF] Bellman | Effects of surface tension and viscosity on Taylor instability[END_REF]

Appendix A. Green's function for the potential of the liquid velocity field

In cylindrical co-ordinates, the Green function G 0 (r, z|r 0 , z 0 that gives the dependence of φ on the flat interface velocity evolution, satisfies the equation:

with boundary conditions:

Multiplying the Equation (A1) by J 0 (γ r), and integrating from r = 0 to ∞, we obtain:

where

is the Hankel transform of order zero of G 0 , with r, r 0 , and z 0 as parameters. Integrating by parts twice, we obtain:

where the integrated terms are each time equal to zero.

Since J 0 (γ r) satisfies the Bessel equation of order zero with parameter γ 2 :

then (A6) becomes:

The boundary conditions for this simple differential equation are:

according to Equations (A2), (A4), and (A7).

The method of images can be used to construct the transformed Green's function for a semi-infinite positive axis, when the left end of this axis is insulated. Thus, G 0H is found by placing a positive image source at z = -z 0 :