Bruno Issenmann

Régis Wunenburger

B I S S E N M A N N

H. C H R A I B I , M. G A N D I L , A N D R W U N E N B U R G E R

J. -P D E L V I L L E

Unsteady deformations of a free liquid surface caused by radiation pressure

come

Unsteady deformations of a free liquid surface caused by radiation pressure

Introduction

The measurement of interfacial and bulk properties of liquids, soft solids and films has motivated the development of several non-contacting techniques in order to overcome the drawbacks of the classical tensiometers and related dynamometric techniques that are their intrusive nature, the associated possible chemical contamination of the samples and the difficulty to study liquid-liquid and soft interfaces. In this context, the analysis of the spatio-temporal dynamics of either thermally or artificially excited small amplitude deformations (so called ripples) of free liquid surfaces, of liquid-liquid interfaces or of surfaces of soft solids has been shown for a long time to be an efficient way for determining in a non invasive manner several interfacial and bulk properties such as surface [START_REF] Sakai | Measurement of liquid surface properties by laserinduced surface deformation spectroscopy[END_REF]) and interfacial [START_REF] Mitani | Measurement of ultralow interfacial tension with a laser interface manipulation technique[END_REF]) tension, bulk viscosity [START_REF] Yoshitake | Measurement of high viscosity with laser induced surface deformation technique[END_REF]), surface [START_REF] Bonfillon-Colin | Contribution l'étude des mousses et des émulsions : dynamique interfaciale[END_REF]) and bulk elasticity [START_REF] Monroy | Direct experimental observation of the crossover from capillary to elastic surface waves on soft gels[END_REF]), adsorption kinetics [START_REF] Sakai | Rapid ripplon spectroscopy with ms time resolution[END_REF]),... While the analysis of the light scattered by thermally excited ripples is the earliest non-contacting technique [START_REF] Langevin | Light scattering by liquid surfaces and complementary techniques[END_REF]; [START_REF] Madsen | Capillary waves at the transition from propagating to overdamped behavior[END_REF]), other powerful methods have been more recently developed. They are based on the controlled generation of ripples by the Maxwell stress due to an electric field produced by a charged blade or needle (Sohl et al. (1978); [START_REF] Langevin | Light scattering by liquid surfaces and complementary techniques[END_REF]; [START_REF] Stenvot | Study of viscoelasticity of soluble monolayers using analysis of propagation of excited capillary waves[END_REF]; [START_REF] Sakai | Electric field tweezers for characterization of liquid surface[END_REF]), or by the radiation pressure (RP) of a laser [START_REF] Komissarova | Light pressure-induced deformations of a free liquid surface[END_REF]; [START_REF] Grigorova | Coherent correlation spectroscopy of capillary waves[END_REF]; [START_REF] Sakai | Measurement of liquid surface properties by laserinduced surface deformation spectroscopy[END_REF]) or of an acoustic beam [START_REF] Khuri-Yakub | Ultrasonic excitation and detection of capillary waves for the measurement of surface film properties[END_REF]; [START_REF] Cinbis | A noncontacting technique for measuring surface tension of liquids[END_REF]). The appeal of all these techniques relies on the key fact that the spatio-temporal dynamics of small amplitude ripples is independent of their amplitude. This indeed makes the accuracy of these techniques independent of dynamometric calibration of the excitation setups.

These non-contacting techniques have already proved their usefulness not only in soft matter science (Mitani & Sakai (2005); [START_REF] Monroy | Direct experimental observation of the crossover from capillary to elastic surface waves on soft gels[END_REF]) but also for the in situ monitoring of the quality of inks used in printers [START_REF] Stockhausen | Method and apparatus for measuring and/or monitoring the surface tension of a fluid[END_REF]), for the size calibration of droplets produced by liquid sample dispensers used in biotechnological applications [START_REF] Williams | Non-contact techniques for measuring viscosity and surface tension information of a liquid[END_REF]), and should constitute promising technologies for applications at high pressure, high temperature or in natural environments [START_REF] Cinbis | Noncontacting techniques for measuring surface tension of liquids[END_REF]).

Given the growing use of the above listed non-contact techniques for bulk and surface metrology, a precise knowledge of the generation, propagation and attenuation of the ripples is required in order to accurately determine the interfacial properties of interest from the monitored surface dynamics. Concerning liquid free surfaces and liquid-liquid interfaces, inertia and surface tension (and gravity to a less relevant extent in this context) are the causes of the propagation of the ripples, while viscosity is responsible for their attenuation and to their eventual over-damping at large dissipation. Various models, either empirical or with a limited range of validity, have been published together with the description of the associated non-contact measurement techniques: simplified analytical [START_REF] Khuri-Yakub | Ultrasonic excitation and detection of capillary waves for the measurement of surface film properties[END_REF]; [START_REF] Cinbis | A noncontacting technique for measuring surface tension of liquids[END_REF]) and numerical [START_REF] Cinbis | Effect of surface tension on the acoustic radiation pressure-induced motion of the water-air interface[END_REF]) models of inviscid dynamics of free liquid surfaces, analytical model of weakly damped dynamics of free liquid surfaces (Ostrovskaya (1988a); [START_REF] Sakai | Measurement of liquid surface properties by laserinduced surface deformation spectroscopy[END_REF]), model of over-damped dynamics of free liquid surfaces with zero Bond number [START_REF] Yoshitake | Measurement of high viscosity with laser induced surface deformation technique[END_REF]), model of free liquid surface over-damped dynamics applied to liquid-liquid interfaces over-damped dynamics [START_REF] Mitani | Measurement of ultralow interfacial tension with a laser interface manipulation technique[END_REF][START_REF] Williams | Non-contact techniques for measuring viscosity and surface tension information of a liquid[END_REF]).

On the other side, since Rayleigh's and Kelvin's seminal works, the physics of ripples propagation has been the subject of several theoretical efforts: Summarizing and extending previous works, Lamb (1932) and [START_REF] Levich | Physicochemical hydrodynamics[END_REF] treated the case of one-dimensional (1D), harmonic waves at the free surface of a liquid, both without and with viscous effects. [START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF] determined the exact solution of the initial-value problem of the propagation of 1D, freely decaying waves at a free surface of a viscous liquid.

The aim of the present work is threefold:

• to derive from this theoretical corpus an exact description of the dynamics of ripples with circular symmetry caused by a spatially localized pulse of acoustic or electromagnetic RP applied on a free liquid surface. We focus our attention on these two kinds of surface forcing because they have already been implemented in efficient and versatile non contact surface measurement techniques [START_REF] Sakai | Measurement of liquid surface properties by laserinduced surface deformation spectroscopy[END_REF]; [START_REF] Cinbis | A noncontacting technique for measuring surface tension of liquids[END_REF]).

• to validate this analytical model through a quantitative comparison with experimental measurements of the dynamics of a free liquid surface forced by pulsed acoustic radiation pressure.

• to validate a numerical code describing the spatio-temporal dynamics of a fluid-fluid interface by comparison with the analytical model. In § 2 we present the characteristic features of the RP excitations under consideration. In § 3 we present and solve the problem of the spatio-temporal dynamics of axially symmetric deformations of a free liquid surface. In § 4 we give a simple physical interpretation of the analytical solution of the problem. In § 5 we present simplified solutions of the problem in the asymptotic regimes of zero, weak and strong dissipation. In § 6 we illustrate the model in the particular case of an acoustic excitation. In § 7 we compare the predictions of the model with measurements of the motion of the hump produced at a free liquid surface by a pulse of acoustic radiation pressure. Finally, in § 8 present a numerical code describing the spatio-temporal dynamics of a fluid-fluid interface and based on the Boundary Element applied to Brinkman's equation and compare its predictions to the analytical model.

Radiation pressure fields

Background assumptions

In this work we focus on the deformation of liquid free surfaces induced by the radiation pressure of either a focused electromagnetic beam with Gaussian profile or an acoustic beam produced by a spherical ultrasonic transducer. Since both kinds of beams are axially symmetric, the resulting interface deformations are also axially symmetric. The focal plane of the exciting beam is assumed to coincide with the initially flat liquid free surface or interface. Cylindrical coordinates (r, θ, z) centered on the vertical beam axis with origin O coinciding with the beam focus and the associated local base (e r , e θ , e z ) will be used throughout this work, see figure 1. Free surface or interface deformations are described by their position and time dependent height h(r, t). In the following, small amplitude deformations are assumed, i.e. (i) their slope with respect to the horizontal is small: ∂h ∂r t (r, t) 1 ∀r, ∀t and (ii) they are much smaller than the characteristic length scale of divergence of the electromagnetic or acoustic field along z. Consequently, (i) the beams are always assumed to impinge on the interface at normal incidence and (ii) the radiation pressure field does not depend on h (r, t).

We now detail the characteristic features of the radiation pressure applied by each kind of beam.

Acoustic excitation

Basically, the radiation stress exerted by an acoustic field on the irradiated interface separating two fluids results from the imbalance between the time averaged Lagrangian stresses exerted on both sides of the interface. According to [START_REF] Borgnis | Acoustic radiation pressure of planar compressional waves[END_REF]; [START_REF] Chu | Acoustic radiation pressure produced by a beam of sound[END_REF]; [START_REF] Lee | Acoustic radiation pressure[END_REF]; Landau & Lifshitz (1987), the acoustic radiation stress exerted by a continuous plane wave of finite width impinging on the interface between a lower fluid 1 and an upper fluid 2 (of density ρ i and sound speed c i , i = 1, 2) is normal to the interface (this is the reason why we call it radiation pressure hereafter) and can be written as π rad = AI inc e z in the case of a horizontal interface irradiated at normal incidence, where I inc is the intensity of the incident acoustic wave:

I inc = p 2 inc t ρ i c i (2.1)
where p inc is the variation of the pressure with respect to the hydrostatic pressure associated with the incident acoustic wave (acoustic pressure), t denotes time averaging, i = 1 (resp. i = 2) if the incident beam propagates upward (resp. downward) and:

A = ε 2 c i ρ 2 i c 2 i + ρ 2 j c 2 j -2ρ i ρ j c 2 i (ρ i c i + ρ j c j ) 2 , (2.2) where ε = 1, i = 1, j = 2 (resp. ε = -1, i = 2, j = 1) if the incident beam propagates upward (resp. downward).
Since intense acoustic fields are required for observing noticeable mechanical effects of the acoustic radiation pressure, we consider axially symmetric acoustic beams emitted by spherically focused ultrasonic transducers. Considering a spherical ultrasonic transducer with focal length z 0 and radius a harmonically excited at frequency f , in its focal plane the emitted acoustic field has the structure of an inhomogeneous plane wave and its harmonic acoustic pressure field p inc (r, 0, t) in the transducer focal plane can be written as [START_REF] Kino | Acoustic waves: devices, imaging and analog signal processing[END_REF]):

p inc (r, 0, t) = P f (r) cos(2πf t + ϕ(r)) (2.3)
where P is the amplitude of oscillation of the acoustic pressure at the focus, ϕ(r) is a r dependent phase, f bs (r) is the beam shape factor at z = 0:

f bs (r) = 2J 1 (k c r) k c r (2.4)
with k c = 2πa λz 0 , λ the acoustic wavelength of the incident beam and J n the n th order Bessel function of the first kind. In the following, we consider a harmonic wave train of central frequency f , duration ∆t, and of square envelope. Thus, in the focal plane of the transducer the radiation pressure field π rad can be written as:

π rad (r, t) = AI inc,0 φ(r)D(t) (2.5)
where I inc,0 = P 2 2ρ i c i is the acoustic intensity of the incident beam at the focus (i is the index of the fluid in which the incident beam propagates), φ(r) = f 2 bs (r) and D(t) = H(t)H(∆t -t) is a door function of duration ∆t, H being the Heaviside function. The Hankel transform φ(k) of φ(r), defined as φ(k) = ∞ 0 rJ 0 (kr)φ(r) dr, has the following expression:

φ(k) = 2 k 2 c 1 -k πk c 1 -k 2 4k 2 c -2 π arcsin k 2k c for k < 2k c φ(k) = 0 for k 2k c (2.6)
The variations of the normalized acoustic intensity distribution in the transducer focal plane k c r → φ(r) and of its normalized Hankel transform in figure 2 (a,b).

Electromagnetic excitation

A photon of frequency ν propagating in a dielectric medium with refractive index n possesses a momentum nhν/c, where c is the speed of light in vacuum and h the Planck constant. The radiation stress exerted by an electromagnetic wave on the irradiated interface separating two dielectric, isotropic media of different refractive indices results from the change of the momentum of the photons crossing the interface [START_REF] Wunenburger | Light-induced deformation and instability of a liquid interface. i. statics[END_REF]). The electromagnetic radiation stress exerted by a wave impinging on the interface between a lower dielectric fluid 1 and an upper dielectric fluid 2 (of respective refractive index n i , i = 1, 2) is normal to the interface [START_REF] Wunenburger | Light-induced deformation and instability of a liquid interface. i. statics[END_REF]) (hence it is called radiation pressure) and can be written as π rad = AI inc e z in the case of a horizontal interface irradiated at normal incidence, where I inc is the intensity of the incident electromagnetic wave, and :

0 2 4 6 0 0.5 1 k c r φ(r) (a) 0 1 2 3 0 0.5 1 k c r φ(r) (c) 0 2 4 0 1 2 k/k c k c 2 φ(k) (b) 0 2 4 0 0.5 1 k/k c k c 2 φ(k) (d)
A = ε 2n i c n j -n i n j + n i (2.7) ε = 1, i = 1, j = 2 (resp. ε = -1, i = 2, j = 1
) if the incident beam propagates upward (resp. downward).

Since continuous wave Gaussian laser beams are usually used (Ostrovskaya (1988b); [START_REF] Sakai | Measurement of liquid surface properties by laserinduced surface deformation spectroscopy[END_REF]; [START_REF] Mitani | Measurement of ultralow interfacial tension with a laser interface manipulation technique[END_REF]; [START_REF] Sakai | Laser excitation of high frequency capillary waves[END_REF][START_REF] Sakai | Rapid ripplon spectroscopy with ms time resolution[END_REF]; [START_REF] Yoshitake | Measurement of high viscosity with laser induced surface deformation technique[END_REF]; Mitani & Sakai (2005); [START_REF] Wunenburger | Light-induced deformation and instability of a liquid interface. i. statics[END_REF]), we assume for I inc the intensity profile of a continuous laser beam in the TEM 00 Gaussian mode in its focal plane (i.e. at its waist):

I(r) = I inc,0 φ(r) (2.8)
where I inc,0 is the electromagnetic intensity of the incident beam at the focus,

φ(r) = exp - 2r 2 ω 2 0 , (2.9)
with ω 0 the beam waist radius. I inc,0 is related to the beam power by the equality

I inc,0 = 2P πω 2 0 .
As for the acoustic case, we consider a monochromatic electromagnetic wave train of central frequency ν, duration ∆t and of square envelope. Thus, the electromagnetic radiation pressure field exerted on the horizontal interface located in the beam focal plane at normal incidence can be written as:

π rad (r, t) = AI inc,0 φ(r)D(t)e z
(2.10)

Defining k c as k c = 2 ω 0 , the Hankel transform φ(k) of φ(r) is:

φ(k) = 1 k 2 c exp - k 2 2k 2 c .
(2.11)

The variations of the normalized electromagnetic intensity distribution at the beam waist and of its normalized Hankel transform in figure 2 (c,d ).

Model of free surface dynamics

Flow equations

We consider a liquid free surface, i.e. a liquid in contact with a gas of viscosity and density both sufficiently small compared to the ones of the liquid so that during the surface deformation, both its inertia and the viscous stresses it exerts on the liquid can be considered as negligible compared to the liquid inertia and the liquid viscous stresses.

The liquid under consideration in this study is Newtonian and has homogeneous density ρ, dynamic viscosity η, and surface tension σ. It is submitted to the acceleration of Earth's gravity g = -ge z , and occupies the lower half-space defined by z < 0 at rest. The upper half-space, defined by z > 0, is filled with a gas whose pressure is assumed to be uniform and constant equal to P 0 . The liquid flow associated with the deformation of its free surface under the effect of radiation pressure is assumed to be characterized by a small Re number and to be incompressible. Notice that this assumption implies that, in the case of an acoustic excitation of the surface, the compressible flow associated with the propagation of the acoustic beam is assumed to be decoupled from the incompressible liquid motion associated with the deformation of the liquid surface, an assumption usually made in studies of streaming flows induced by bulk acoustic absorption (see [START_REF] Kamakura | Acoustic streaming induced in focused gaussain beams[END_REF] and references therein). Therefore, the liquid velocity field v(r, z, t) associated with the surface deformation is solenoidal:

∇ • v = 0 (3.1)
and satisfies the linearized Navier-Stokes equation:

ρ ∂v ∂t = -∇P s + η∇ 2 v + ρg. (3.2)
where P s is the liquid pressure. The fluid is put out into axially symmetric motion by the axially symmetric pulse of radiation stress π rad(r,t) exerted on the free surface of altitude h(r, t). The kinematic boundary condition at the free surface is :

∂h ∂t (r, t) = v(r, h(r, t)) • n(r, t) (3.3)
where n(r, t) is the vector normal to the interface oriented from the liquid to the gas.

The stress balance at the free surface is written as:

-P 0 n(r, t) -T(r, h(r, t), t) • n(r, t) + π rad (r, t) + σκn(r, t) = 0 (3.4)
where T = -P s I + 2ηD(v) is the hydrodynamic stress tensor expressed for the incompressible flow of a Newtonian fluid, I the identity tensor, D(v) = 1 2 (∇v + t ∇v), and κ(r, t) is the surface curvature. Introducing the dynamic pressure p = P s -P 0 + ρgz in (3.2) and (3.4), both become:

ρ ∂v ∂t = -∇p + η∇ 2 v (3.5)
and (p(r, h(r, t), t) -ρgh(r, t))n(r, t) -2ηD(v)(r, h(r, t), t) • n(r, t) + π rad (r, t) + σκn(r, t) = 0, (3.6) respectively. Finally, reminding that the radiation pressure is exerted from t = 0, the liquid is initially at rest and its free surface is initially flat.

Implications of the assumption of small amplitude deformations

The assumption of small amplitude deformations, as defined in § 2.1, allows to also considerably simplify the expressions of the boundary conditions at the free surface, since then: (i) the surface curvature κ can be made linear in h:

κ(r, t) 1 r ∂ ∂r r
∂h ∂r (r, t).

(3.7)

(ii) n can be approximated by e z and (iii) evaluating (3.6) at the free surface is equivalent to evaluating it at z = 0. Thus, since v is assumed to be independent of θ and to lie in the meridian plane (v θ = 0), the projections of (3.6) along e r and e z , respectively, are:

η ∂v r ∂z + ∂v z ∂r (r, 0, t) = 0, (3.8) p(r, 0, t) -ρgh(r, t) -2η ∂v z ∂z (r, 0, t) + π rad (r, t) + σ 1 r ∂ ∂r r
∂h ∂r (r, t) = 0. (3.9) (iii) the kinematic boundary condition (3.3) can also be simplified:

∂h ∂t r (r, t) = v z (r, 0, t) (3.10)
3.3. Flow decomposition Following Lamb (1932) and [START_REF] Levich | Physicochemical hydrodynamics[END_REF], we solve this problem by expressing the flow as the superposition of a solenoidal, irrotational flow whose velocity v (p) satisfies the Euler equation and of a solenoidal, rotational flow whose velocity field v (r) satisfies an equation of pure diffusion of momentum involving no pressure gradient, i.e.

v = v (p) + v (r) (3.11) with ρ ∂v (p) ∂t = -∇p (3.12) and ρ ∂v (r) ∂t = η∇ 2 v (r) . (3.13)
Expressing p as the sum of irrotational and rotational contributions p (p) and p (r) , respectively, and assuming ρ ∂v (p) ∂t = -∇p (p) instead of (3.12) and ρ ∂v (r) ∂t = -∇p (r) + η∇ 2 v (r) instead of (3.13), [START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF] evidenced the pressure correction induced by the rotational flow to the pressure of the potential flow. But since only the sum p of both pressure terms actually matters for predicting the free surface shape, following Lamb (1932) and [START_REF] Levich | Physicochemical hydrodynamics[END_REF] we do not split p and we use (3.12) and (3.13), despite the fact that the method used in [START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF] is physically more enlighting.

In order to solve this set of equations, we introduce the velocity potential φ associated with v (p) and the stream function ψ associated with v (r) . Since φ is defined by

v (p) = ∇φ (3.14)
it satisfies:

ρ ∂φ ∂t = -p. (3.15)
Furthermore, the solenoidal property of v (p) leads to:

∇ 2 φ = 0. (3.16)
On the other hand, the vorticity r) also satisfies a viscous diffusion equation:

w = ∇ ∧ v = ∇ ∧ v (
∂w ∂t = ν∇ 2 w. (3.17)
where ν = η ρ is the liquid kinematic viscosity. Because of the axial symmetry of v, w is orthoradial: w = we θ . The stream function ψ associated with v (r) satisfies:

v (r) r = ∂ψ ∂z (3.18) v (r) z = - 1 r ∂(rψ) ∂r .
(3.19) because of the solenoidal property of v (r) . From (3.18-3.19) and the definition of vorticity we deduce the following relation between w and ψ:

w = ∂ ∂r 1 r ∂(rψ) ∂r + ∂ 2 ψ ∂z 2 .
(3.20)

Hankel-Laplace transform

Since the response of the surface to the pulse of radiation pressure is expected to be harmonic neither in time nor in space, it involves a continuum of wavelengths and of frequencies. Therefore, following Ostrovskaya (1988a), we apply Hankel and Laplace transforms to the hydrodynamic fields. Assuming k ∈ R + , we define:

V (α) r (k, z, s) = ∞ 0 e -st ∞ 0 rJ 1 (kr)v (α) r (r, z, t) dr dt, α = p, r (3.21a) V (α) z (k, z, s) = ∞ 0 e -st ∞ 0 rJ 0 (kr)v (α) z (r, z, t) dr dt, α = p, r (3.21b) Φ(k, z, s) = ∞ 0 e -st ∞ 0 rJ 0 (kr)ϕ(r, z, t) dr dt (3.21c) W (k, z, s) = ∞ 0 e -st ∞ 0 rJ 1 (kr)w(r, z, t) dr dt (3.21d ) Ψ(k, z, s) = ∞ 0 e -st ∞ 0 rJ 1 (kr)ψ(r, z, t) dr dt (3.21e) P (k, z, s) = ∞ 0 e -st ∞ 0 rJ 0 (kr)p(r, z, t) dr dt (3.21f ) H(k, s) = ∞ 0 e -st ∞ 0 rJ 0 (kr)h(r, t) dr dt (3.21g) Π(k, s) = ∞ 0 e -st ∞ 0 rJ 0 (kr)π rad (r, t) dr dt (3.21h)
These transformations lead to rewrite (3.14), (3.16), (3.15), (3.17), (3.18), (3.19), (3.20), (3.10), (3.8), (3.9) in the form:

V (p) r = -kΦ (3.22a) V (p) z = ∂Φ ∂z (3.22b) -k 2 Φ + ∂ 2 Φ ∂z 2 = 0 (3.22c) -sΦ = P ρ (3.22d ) sW = ν -k 2 + ∂ 2 ∂z 2 W (3.22e) V (r) r = ∂Ψ ∂z (3.22f ) V (r) z = -kΨ (3.22g) W = ∂ 2 ∂z 2 -k 2 Ψ (3.22h) sH = V z (k, z = 0, s) (3.22i ) ∂V r ∂z (k, z = 0, s) -kV z (k, z = 0, s) = 0 (3.22j ) P (k, z = 0, s) -ρgH(k, s) -2η ∂V z ∂z (k, z = 0, s) + Π(k, s) = σk 2 H(k, s) (3.22k )
(3.22e) and (3.22h) imply that, to a gauge function Ψ g satisfying ∂ 2 ∂z 2 -k 2 Ψ g = 0, i.e. describing an irrotational, hence irrelevant flow, Ψ also satisfies a viscous diffusion equation:

sΨ = ν ∂ 2 ∂z 2 -k 2 Ψ (3.23)
The solutions of (3.22c) and (3.23) are of the form Φ = A(k)e kz and Ψ = C(k)e lz , respectively, with and Ψ in (3.22d, 3.22i, 3.22j, 3.22k) allows to express H as a function of Π:

l 2 = k 2 (1 + s νk 2 ) and A(k), C(k) ∈ R. Injecting these expressions of Φ
H(k, s) = k (νk 2 ) 2 1 F (X) Π ρ , (3.24)
where

X 2 (k, s) = l 2 k 2 = 1 + s νk 2 , (3.25) F (X) = (X 2 + 1) 2 -4X + ∆, (3.26) with ∆(k) = Ω(k) 2 (νk 2 ) 2 , (3.27)
where

Ω(k) = gk + σ ρ k 3 . (3.28)
From a physical point of view, Ω(k) is the pulsation of plane surface waves of real wave number k propagating at the surface of an inviscid fluid, and (νk 2 ) -1 is the characteristic timescale of viscous diffusion of momentum over the length scale k -1 (Lamb (1932); [START_REF] Levich | Physicochemical hydrodynamics[END_REF]). F (X(k, ω)) = 0 where ω = is is the dispersion relation of plane surface waves of real wave number k and complex pulsation ω at the surface of a viscous fluid.

Analytical solution

The surface deformation height h(r, t) can be determined by inverting (3.21g):

h(r, t) = ∞ 0 kJ 0 (kr)L -1 (H(k, s)) dk (3.29)
where L -1 (H(k, s)) is the inverse Laplace transform of H(k, s). Given (3.24), (2.5) and

(2.7), h(r, t) can be rewritten as:

h(r, t) = AI inc,0 1 ρν 2 ∞ 0 J 0 (kr) k 2 φ(k)L -1 L(D)(s) F (X(k, s)) dk (3.30)
where

L(D)(s) is the Laplace transform of D(t). Given the definition of D(t), L -1 L(D)(s) F (X(k,s))
is calculated using the convolution theorem:

L -1 L(D)(s) F (X(k, s)) = t 0 D(τ )L -1 1 F (X(k, s)) (t -τ ) dτ, (3.31) = t 0 L -1 1 F (X(k, s)) (τ ) dτ if t < ∆t, (3.32) = t t-∆t L -1 1 F (X(k, s)) (τ ) dτ if t ∆t. (3.33)
In order to determine L -1 1 F (X(k,s)) , following [START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF], we exploit the polynomial form of F (X(k, s)) of the variable X by performing a partial fraction expansion of 1/F (X):

1 F (X) = 4 i=1 a i (k) X -X i (k) (3.34)
where X i (k), i = 1 -4 are the roots of the polynomial F (X),

X i (k) = α i 2 - 4 3 + 8 + 6∆(k) 3G(k) + 2 3 G(k)+ β i 2 - 8 3 - 8 + 6∆(k) 3G(k) - 2 3 G(k) + γ i 8 -4 3 + 8+6∆(k) 3G(k) + 2 3 G(k) , (3.35) with G(k) = 19 -9∆(k) + 3 33 -54∆(k) -3∆(k) 2 -3∆(k) 3 1/3 , (3.36) α 1,2,3,4 = (-1, -1, 1, 1), β 1,2,3,4 = (-1, 1, -1, 1), γ 1,2,3,4 = (-1, -1, 1, 1
), and

a i (k) = 1 Π j =i (X i (k) -X j (k)) . (3.37) As L -1 1 √ s-a = 1 √ πt + a exp(a 2 t) 1 + erf(a √ t) , L -1 (L(Q)(s + 1)) = exp(-t)Q(t) and L -1 (L(Q)( s ω )) = ωQ(ωt), L -1 1 F (X(k, s)) = ω(k) exp(-ω(k)t) πω(k)t 4 i=1 a i (k) + ω(k) exp(-ω(k)t) 4 i=1 a i (k)X i (k) exp X i (k) 2 ω(k)t 1 + erf(X i (k) ω(k)t) .
Since 4 i=1 a i (k) = 0, as shown in [START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF], the general integral

G(k, t) of L -1 1 F (X(k,s)) is: G(k, t) = 4 i=1 a i (k)X i (k) X 2 i (k) -1 exp(ω(k)t(X 2 i (k) -1)) × 1 + erf(X i (k) ω(k)t) -1 -X i erf( ω(k)t) . (3.38) Finally, h(r, t) = AI inc,0 1 ρν 2 ∞ 0 J 0 (kr) k 2 φ(k)(G(k, t) -G(k, 0)) dk if t < ∆t = AI inc,0 1 ρν 2 ∞ 0 J 0 (kr) k 2 φ(k)(G(k, t) -G(k, t -∆t)) dk if t ∆t (3.39)
The computation of h(r, t) is performed using the dimensionless wave number k = k/k c :

h(r, t) = AI inc,0 σk 2 c ∞ 0 J 0 (k c r k)(k 2 c φ(k c k)) 1 k + Bo c k ∆( k)(G(k, t) -G(k, 0)) d k if t < ∆t = AI inc,0 σk 2 c ∞ 0 J 0 (k c r k)(k 2 c φ(k c k)) 1 k + Bo c k ∆( k)(G(k, t) -G(k, t -∆t)) d k if t ∆t (3.40) where ∆( k) = σ ρ k 3 c (νk 2 c ) 2 1 k + Boc k3 , Bo c = ( 2 k 2 c ) -1
being the Bond number associated with the characteristic wave number k c and = σ ρg the capillary length associated with the liquid surface. We now give a physical interpretation of this analytical solution.

Physical interpretation

Flow type

The variations of X i , i = 1 -4 as a function of ∆ are shown in figure 3 (a,b). The association of these roots in one or two pairs of conjugate complex roots is due to the fact that the coefficients of F (X) are real. While X 1 and X 2 always form a pair of conjugate complex numbers, X 3 and X 4 are real for ∆ < ∆ 0 with ∆ 0 0.582 and complex for ∆ > ∆ 0 .

The corresponding dependence of the growth rates 3 (c,d ). For ∆(k) < ∆ 0 , s 1 and s 2 are conjugate complex numbers with (negative) real part and imaginary part both scaling as νk 2 . The corresponding flow is thus oscillatory and overdamped since it is characterized by an oscillation period comparable to its damping time scale. For ∆(k) > ∆ 0 , the imaginary parts of s 1 and s 2 both scale as νk 2 √ ∆ = Ω(k), while their negative real part still scales as νk 2 . The corresponding flow is thus oscillatory and weakly damped: circular surface waves of wave number k propagate with small damping. For ∆(k) < ∆ 0 , both s 3 and s 4 are negative reals. The corresponding flow and surface motion are simply damped. Since s 3 scales as νk 2 and s 4 scales as Ω 2 (k)/(νk 2 ), in the ∆(k) → 0 limit |s 4 | is asymptotically much smaller than |s 3 |, |Re(s 1 )| and |Re(s 2 )|. Thus, the flow associated of characteristic damping time scale s -1 4 is damped asymptotically much more slowly than the three other flows of characteristic damping time scales s -1 i , i = 1 -3. For ∆(k) > ∆ 0 , s 3 and s 4 behave asymptotically like s 1 and s 2 . The corresponding surface motion is propagative and weakly damped.

s i = νk 2 (1 -X 2 i ) on ∆ is shown in figure
We conclude from this analysis that for ∆(k) < ∆ 0 the surface motion is globally damped with Ω 2 (k)/(νk 2 ) as the leading damping time scale while for ∆(k) > ∆ 0 circular surface waves with wave number k propagate with pulsation Ω(k) and are weakly damped with damping time scale (νk 2 ) -1 .

Flow spatial extension

The flow is the superposition of an irrotational flow and a viscous flow, whose velocity fields are (V

(p) r = -kA(k) exp(kz), V (p) z = kA(k) exp(kz)) and (V (r) r = lC(k) exp(lz),

V (r) z

= -kC(k) exp(lz)), respectively. The resolution of the hydrodynamic problem performed in 3.4 implies that the viscous to irrotational velocity ratios corresponding to each root X i of F (X) have the following expressions:

V (r) r /V (p) r = 2X i /(1 + X 2 i ) and V (r) z /V (p) z = 2/(1 + X 2 i ).
The variations versus ∆ of these velocity ratios are shown in figure 4 (a). When ∆ → ∞, all the velocity components of the viscous flow become small compared to the velocity components of the irrotational flow. Simultaneously, k/l = 1/X i → 0 for all X i , as shown in figure 4 (b), i.e. the irrotational flow spreads more deeply within the liquid than the viscous flow: for ∆ → ∞ the flow is almost everywhere irrotational except in a thin layer close to the surface, where a flow contribution of small amplitude due to the viscous diffusion of momentum arises. On the contrary, when ∆ → 0, all the velocity fields are of comparable amplitude and spatial extend: during the deformation of the surface, the viscous diffusion of momentum has enough time to occur down to the penetration depth of the inertial flow within the liquid, therefore affecting it noticeably.

Asymptotic behaviours

Identifying the dynamic regime to which a given experiment corresponds requires to evaluate the typical value of ∆ characterizing the surface deformation. Since k c is the typical value of k within the surface wave packet excited both acoustically and electromagnetically, as shown by (2.6) and (2.11), the typical value of ∆ is ∆ c = ∆(k c ). ∆ c , which is ratio of the frequency of inviscid surface waves of wavenumber k c and of the attenuation rate of surface waves due to viscosity, is usually called the frequency parameter. If ∆ c 1 (resp. 1), the dynamics of the surface deformation is expected to be overdamped (resp. propagative and weakly damped). One may advance that, since viscous dissipation increases with k, evaluating ∆ at k = k c underestimates viscous effects. But since φ(k) vanishes for large k in both the acoustic and electromagnetic cases (in the acoustic case, it is even strictly equal to zero for k > 2k c ), the influence of contributions corresponding to values of k much larger than k c on the flow is expected to be negligible.

The inviscid case

For low viscosity liquids and broad beams such that ∆ c 1, we expect the liquid surface response to be well described by a model of inviscid flow. Considering ν = 0 in § 3, the flow reduces to a potential flow. The corresponding surface deformation height can be straightforwardly deduced:

h(r, t) = AI inc,0 1 ρ ∞ 0 k 2 J 0 (kr) φ(k)(G(k, t) -G(k, 0)) dk if t < ∆t, = AI inc,0 1 ρ ∞ 0 k 2 J 0 (kr) φ(k)(G(k, t) -G(k, t -∆t)) dk if t ∆t. (5.1)
where

G(k, t) = 1 Ω 2 (k) (1 -cos(Ω(k)t)).
(5.

2)

The spatio-temporal evolution of the liquid surface height appears as the result of the propagation of a wave packet of undamped surface eigenmodes of wave number k and pulsation Ω(k) with k ranging from 0 to typically k c . Using the dimensionless wave number k:

h(r, t) = AI inc,0 σk 2 c ∞ 0 J 0 (k c r k)(k 2 c φ(k c k)) 1 k + Bo c k (1 -cos(Ω( k)t)) d k if t < ∆t = AI inc,0 σk 2 c ∞ 0 J 0 (k c r k)(k 2 c φ(k c k)) 1 k + Boc k (cos(Ω( k)(t -∆t)) -cos(Ω( k)t)) d k if t ∆t (5.3)
where Ω( k) = Considering the case of a pulse wise excitation of the surface (i.e. such that ∆t Ω -1 c ), for t > ∆t (5.3) can be rewritten as:

h(r, t) = AI inc,0 σk 2 c ∞ 0 J 0 (k c r k)(k 2 c φ(k c k)) 1 k + Bo c
k Ω(k)∆t sin(Ω( k)t) d k.

(5.4)

Consequently, since the integrand in (5.4) is of the order of Ω(k)∆t and k varies between 0 and O(1), in the range of small values of the Bond number, the characteristic height of the unsteady surface deformation h is predicted to scale as

AIinc,0 σk 2 c Ω c ∆t = AI inc,0 ∆t 1 ρσk c .

The case of weak viscous dissipation

Still trying to grasp the effect of viscous dissipation on the surface deformation dynamics in the ∆ c 1 limit, i.e. in situations for which the typical damping time constant of the surface motion ω -1 c , where ω c = ω(k c ) = νk 2 c , is much larger than the typical surface oscillation period Ω -1 c , we can simplify the resolution of the flow problem by exploiting the asymptotic expressions of the physical roots of F (X). As shown in figure 3 (a,b), in the ∆ c 1 limit, both the real and imaginary parts of X 3 and X 4 scale as ∆ 1/4 . Thus, following Ostrovskaya (1988a), since in the ∆ 1 limit

|X i | n |X i | 1, n = 2
, 4 and i = 3, 4, we simplify (3.26) to a bi-quadratic form:

F (X) (X 2 + 1) 2 + ∆.
(5.5)

In the ∆ 1 limit, the roots of (5.5) are:

X 3 = -1 -i √ ∆ = ∆ 1/4 √ 2 (1-i)+O(∆ -1/4
) and

X 4 = -1 + i √ ∆ = ∆ 1/4 √ 2 (1 + i) + O(∆ -1/4
), in agreement with their direct calculation shown in the insets of figure 3 (a,b). In the ∆ 1 limit, the corresponding values of the growth rates have the following expressions s 3 = -2νk 2 -iΩ(k) and s 4 = -2νk 2 + iΩ(k) . The corresponding surface deformation dynamics is therefore expected to be mainly inertial, i.e. oscillatory and weakly damped, with a characteristic inertial time scale of the order Ω(k c ) -1 and a damping time scale of the order of (νk 2 c ) -1 Ω(k c ) -1 . The major interest of this simplification is to transform F (X) into a polynomial function of s:

F (X(k, s)) (s + 2) 2 + ∆ = (s -s 3 )(s -s 4 ) (5.6)
i.e. an analytic function of s, so that the resulting surface deformation dynamics can be approximated with a simple superposition of exponential functions of time. This simplification means that wa can neglect the memory term in the harmonic oscillator-like equation describing the evolution of the amplitude of each surface eigenmode established in [START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF] (equation (18) of [START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF]). By proceeding in the same way as in § 3.5 and using L -1 1 s-a = exp(-at), we deduce that the surface deformation height satisfies (3.39) with:

G(k, t) = ω 2 (k) Ω(k) exp(-2ω(k)t)(-2ω(k) sin(Ω(k)t) -Ω(k) cos(Ω(k)t)) + Ω(k) 4ω 2 (k) + Ω 2 (k) (5.7)
The spatio-temporal evolution of the liquid surface height appears as the result of the propagation of a wave packet of weakly damped surface eigenmodes of wave number k, pulsation Ω(k) and damping time constant ω(k) with k ranging between 0 and typically k c . Using the dimensionless wave number k:

h(r, t) = AI inc,0 σk 2 c ∞ 0 J 0 (k c r k)(k 2 c φ(k c k)) 1 k + Boc k ( G(k, t) -G(k, 0)) d k if t < ∆t = AI inc,0 σk 2 c ∞ 0 J 0 (k c r k)(k 2 c φ(k c k)) 1 k + Boc k ( G(k, t) -G(k, t -∆t)) d k if t ∆t (5.8)
where

G(k, t) = 1 -exp(-2ω(k)t)(cos(Ω(k)t) + 2 √ ∆(k) sin(Ω(k)t)) 1 + 4 ∆(k)
(5.9) Thus, still considering Bo c 1, the characteristic time scale of variation of the surface deformation is Ω -1 c .

The case of strong viscous dissipation

In the ∆ c 1 limit, the surface motion is expected to be essentially damped. As shown in figure 3 (c,d ), since |s 1 |, |s 2 | and |s 3 | scale as νk 2 (more precisely, to zeroth order in ∆ in the ∆ → 0 limit, s j νk 2 (-3.5437-(-1) j 2.2303 i), j = 1, 2 and s 3 -0.9126 νk 2 ), the corresponding damped flows are expected to vanish with a characteristic time constant τ that scales as ω -1 c . On the contrary, since in the ∆ c 1 limit -s 4 = νk 2 ∆ 2 νk 2 , the corresponding flow is expected to be damped with a characteristic time constant τ that scales as

Ω 2 c ω(k c ) -1
∝ η/(σk c ) in the Bo c 1 limit. Since τ τ , as soon as t τ , the surface motion is expected to be characterized by a single time constant, namely τ . More precisely, since in the |z| → ∞ limit (with

|arg(z)| < 3π 4 ) exp(z 2 )erf(z) ∼ exp(z 2 )-1 √ π 1 z (1+O(|z| -4 )),
G(t) = 4 i=1 a i X i X 2 i -1 exp(s i t) 2 - 1 √ π 1 X i √ ωt (1 + O(|X i √ ωt| -4 )) -1 -X i erf( √ ωt) .
(5.10) Thus, since in the ∆ → 0 limit

|s 4 | ∝ Ω 2 ω |s i | ∝ ω, i = 1 -3, in the range t ω -1
, G(t) can be rewritten as:

G(t) = 3 i=1 -a i X i X i -1 + ω a 4 X 4 s 4 (2 exp(s 4 t) -1 -X 4 ).
(5.11) Thus,

h(r, t) = AI inc,0 σk 2 c ∞ 0 J 0 (k c r k)(k 2 c φ(k c k)) 1 k + Bo c k (1 -exp(s 4 (k)t)) d k if t < ∆t = AI inc,0 σk 2 c ∞ 0 J 0 (k c r k)(k 2 c φ(k c k)) 1 k + Boc k ∆(k) 2 exp - ∆(k) 2 ω(k)t ω(k)∆t d k
if t ∆t and ∆t τ (5.12)

Consequently, since the integrand in (5.12) is of the order of ∆(k)ω(k)∆t, in the range of small values of the Bond number, the characteristic height of the unsteady surface deformation h is predicted to scale as

AI inc,0 σk 2 c Ω 2 c νk 2 c ∆t = AI inc,0 1 ηkc ∆t.

Illustration of the surface dynamics

We now illustrate the above given theoretical predictions and demonstrate their consistency in the particular case of an acoustic surface excitation. In the measurement techniques cited above, the surface deformation dynamics is detected by measuring (i) the surface height by confocal microscopy [START_REF] Cinbis | A noncontacting technique for measuring surface tension of liquids[END_REF]), acoustic timeof-flight measurement [START_REF] Khuri-Yakub | Ultrasonic excitation and detection of capillary waves for the measurement of surface film properties[END_REF]; [START_REF] Williams | Non-contact techniques for measuring viscosity and surface tension information of a liquid[END_REF]) or holography (Ostrovskaya (1988b)), (ii) the deflection of a light beam reflected by the liquid surface acting as a tilted mirror (Sohl et al. (1978); [START_REF] Grigorova | Coherent correlation spectroscopy of capillary waves[END_REF]; [START_REF] Bonfillon-Colin | Contribution l'étude des mousses et des émulsions : dynamique interfaciale[END_REF]; [START_REF] Mitani | Measurement of ultralow interfacial tension with a laser interface manipulation technique[END_REF]; [START_REF] Sakai | Laser excitation of high frequency capillary waves[END_REF]; [START_REF] Sakai | Electric field tweezers for characterization of liquid surface[END_REF]) that is proportional to the surface slope, (iii) the defocus of a light beam reflected by the deformation tip acting as a spherical mirror [START_REF] Sakai | Measurement of liquid surface properties by laserinduced surface deformation spectroscopy[END_REF] (2005)) that is proportional to the tip curvature. Thus, we will study the behavior of both the surface height h(r, t) and its curvature κ(r, t). Given (3.7), expressions for κ(r, t) are straightforwardly derived from those for h(r, t) by multiplying the integrand by -k 2 in (3.39, 5.1) and by -k 2 c k2 in (3.40, 5.3, 5.8, 5.12). We consider a beam emitted by a focused transducer of numerical aperture NA = 2a z0 = 1 working at frequency f = 2.25 MHz and propagating in a fluid of sound speed c = 1500 m.s -1 , for which k c = π/λ = 4712 m -1 . The chosen wave train duration is ∆t = 8.88 µs, i.e. 20 acoustic periods. The liquid density is ρ = 1000 kg.m -3 and has a variable viscosity ranging from 3.92 × 10 -3 to 1.24 × 10 -6 m 2 • s -1 . The Bond number is constant and equal to 6.1 × 10 -3 .

We define a characteristic surface deformation height h c as h c = AIinc,0 σk 2 c and a characteristic surface deformation curvature κ c as κ c = k 2 c h c . Since the surface curvature is directly related to the radiation pressure forcing through (3.9), we begin our illustration of the surface dynamics by analyzing the spatio-temporal variations of the dimensionless curvature κ(r, t)/κ c in the inviscid case. The r variation of κ(r, t)/κ c for several values of t are shown in figure 5. At the bottom of the graph the acoustic beam shape factor f bs (r) is compared with the vertically expanded surface curvature profile at early stage of the surface deformation r → κ(r, t = 30µs). Notice the coincidence between the distance between two consecutive side lobes of f (r) and the pseudo-wavelength of the wavy curvature profile (this initial wavy curvature profile of the liquid surface is not observed in the case of a Gaussian, electromagnetic forcing). Also notice the propagation of the humps of the curvature profile away from the beam axis as divergent circular waves and the oscillation of the sign of the curvature at r = 0.

In figure 6 tip t → κ(0, t)/κ c in the inviscid case is shown. Given our last remark on figure 5, we interpret its oscillations as the consequence of the successive arrivals at r = 0 of convergent circular waves induced at t = 0 by the side lobes of f (r). Figure 6 displays the time dependence of the dimensionless curvature of the deformation tip for several values of ∆ c ranging from 10 -3 to 10 4 predicted using the general model (3.40). These predictions are confronted to their overdamped approximation in the ∆ c 1 limit and to their weakly viscous (5.8) and inviscid (5.3) approximations in the ∆ c 1 limit. Notice the good agreement between the general model and (i) the overdamped approximation for ∆ c 10 -3 , as shown in figure 6 (a,b), (ii) the weakly viscous and inviscid approximations for ∆ c 10 3 , as shown in figure 6 (d ), that defines a broad range 10 -3 < ∆ c < 10 3 for the relevance of the general model. Also notice the damping of the curvature oscillations as ∆ c decreases leading eventually to their disappearance for ∆ c < 10 2 .

We repeat this analysis for the surface deformation height. Figure 7 displays the time dependence of the dimensionless height of the surface deformation tip, h(0, t)/h c , for several values of ∆ c ranging from 10 -3 to 10 4 predicted using the general model (3.40). These predictions are confronted to their overdamped approximation in the ∆ c 1 limit and to their weakly viscous (5.8) and inviscid (5.3) approximations in the ∆ c 1 limit. Notice again the good agreement between the general model and (i) the overdamped approximation for ∆ c 10 -3 , as shown in figure 7 (a), (ii) the weakly viscous and inviscid approximations for ∆ c 10 3 , as shown in figure 6 (d ), that confirms the broad range 10 -3 < ∆ c < 10 3 of relevance of the general model for h(0, t) too.

Figure 8a displays the ∆ c dependence of the maximal dimensionless tip height h max /h c reached by h(0, t)/h c and of the maximal dimensionless tip curvature κ max /κ c reached by κ(0, t)/κ c . As predicted by the inviscid model ( § 5.1), in the ∆ c 1 range both h max and κ max scales as Ω c ∆t and are thus independent of ∆ c in the ∆ c 1 range. Figure 8 (b) displays the ∆ c dependence of the time t h max (resp. t κ max ) at which h(0, t) (resp. κ(0, t)) attains its maximum h max (resp. κ max ). As predicted by the inviscid model ( § 5.1), t h max and t κ max both scale as Ω -1 c and are thus independent of ∆ c in the ∆ c 1 range. Furthermore, as predicted by the overdamped model ( § 5.3), in the ∆ c 1 range both h max and κ max scale as 1 ηkc and thus scale as √ ∆ c , as shown in figure 8 (a).

The overdamped model ( § 5.3) also predicts a late stage exponential relaxation of the surface deformation with a single characteristic time scale τ = η σk c . Figure 9 displays the variations of h(0, t)/h c and of κ(0, t)/κ c (in logarithmic scale) as a function of t/η for several values of ∆ c falling in the overdamped regime. As expected, log(h(0, t)/h c ) and log(κ(0, t)/κ c ) decrease linearly in t/η with a slope that is independent of ∆ c . Finally, notice that in the ∆ c 1 range the time t h max at which h(0, t) attains its maximum, that corresponds to the transition from the regime of viscous diffusion of vorticity, of time constant τ = ω -1 c , to the regime of relaxation of the surface deformation, of time constant τ = ∆ -1 c ω -1 c , slowly increases as ∆ c decreases, whereas the time t κ max at which κ(0, t) attains its maximum slowly decreases as ∆ c decreases.

Comparison with experiments

The goal of the experiments presented in this section is to test the validity of the analytical model presented above. Our setup is based on an acoustic excitation of the free liquid surface and on the optical detection of the surface curvature at the tip of the deformation. The bumped liquid surface acts as a partially reflecting spherical mirror that focuses/defocuses a laser beam impinging on it. The focusing/defocusing of the reflected beam is detected by simple spatial filtering. First, we present the experimental setup, then we compare our measurements of the temporal evolution of the curvature of the deformation tip generated at the surface of several liquids of variable viscosity with our analytical model.

Experimental setup

The liquid sample is contained in a glass tank closed by a PTFE bottom cap, see figure 10. A vertically immersed Imasonic spherical ultrasonic transducer (UT) of focal length F = 38 mm, numerical aperture 1, central frequency f = 2.25 MHz, bandwidth 600 kHz, focal zone length 3.2 mm, is mounted through the cap. It periodically emits sine wave trains toward the free liquid surface located in its focal plane. The acoustic wave trains are composed of 20 cycles at carrier frequency f and are emitted at a typical repetition rate of 50 Hz. After a propagation delay τ = F/c, where c is the sound speed in the liquid, each wave train induces an axially symmetric surface deformation (hump) whose characteristic diameter is of the order of 1 mm. The hump raises then relaxes over a characteristic time scale of the order of 1 ms. The waist formed by a set of lenses (not shown in figure 10) of an initially circularly polarized 17 mW cw TEM 00 He-Ne laser beam (wavelength in vacuum λ 0 = 0.6328 µm) is located at the tip of the interface deformation. A Glan prism (GP), a beam splitter cube (BSC) and mirrors (M) allow to linearly polarize the beam and to direct it toward the liquid surface. Since the incident laser beam impinges on the surface at normal incidence, a quarter-wave plate (λ/4) is used to extract the part of the laser beam that has been reflected by the tip of the acoustic RP induced surface deformation. This reflected laser beam is then filtered through a tunable iris diaphragm (ID) of typical diameter 1 mm distant of D = 0.79 m from the liquid surface. Finally, it is focused on an amplified photodiode (PD) of bandwidth 170 kHz. The photodiode voltage signal is acquired in the DC mode and averaged over 500 sweeps using a digital oscilloscope of bandwidth 300 MHz. We now explicit the relation between the hump tip curvature κ(0, t) and the optical power P(t) of the reflected and filtered laser beam which is measured using the photodiode. The condition at which the illuminated area of the hump can be considered as a surface of homogeneous (r-independent) curvature is that the diameter of the laser beam at the altitude of the liquid surface is much smaller than the characteristic diameter of the surface deformation, namely λ for the transducer used in this study [START_REF] Issenmann | Bistability of a compliant cavity induced by acoustic radiation pressure[END_REF]). If this condition is fulfilled, the hump tip can be actually considered as a spher-ical mirror of curvature κ = κ(0, t). Since the beam waist radius ω 0 , defined by (2.9), is equal to 44 µm, this condition is satisfied not only when the waist of the laser beam is located at the liquid surface (z = 0), but also at close distance from it, say a distance comparable to its natural diffraction length (Rayleigh range) z 0 = πω 2 0 λ0 = 10 mm. So, considering that the beam waist of the laser beam impinging on the liquid surface is distant of z from it (z is positive if the laser beam waist is located in air, i.e. in front of the liquid surface), the power P of the part of the reflected beam that passes through the diaphragm of radius R I and located at D + z from the mirror can be written as [START_REF] Saleh | Fundamentals of photonics[END_REF]):

P = P 0 1 -exp - 2R 2 I ω 2 (D + z -z ) (7.1)
where P 0 is the total optical power of the reflected laser beam, z

= -κ -1 + M 2 (z + κ -1 )
is the distance from the waist of the reflected beam to the mirror,

M (z) = -κ -1 z+κ -1 (1+ z 0 z+κ -1 ) 1/2
, and ω 0 (z) = M (z)ω 0 is the waist radius of the reflected beam. The dependence of ∆P * = (P(κ)-P(0))/P(0) versus κ for several values of z ranging from -2z 0 to 2z 0 is represented in figure 11. The sensitivity of ∆P to κ is observed to strongly depend on z (amplitude and sign of variation) and displays a non monotonic behavior. Nevertheless, is it found to always depend linearly on κ for |κ| 10 m -1 . This means that in the range of I inc,0 in which κ(0, t) is linear in I inc,0 and ∆P * is linear in κ, (i) the photodiode voltage is linear in I inc,0 , or equivalently, the shape of the photodiode voltage signal is independent of I inc,0 , and (ii) ∆P * (t) is simply proportional to κ(0, t). This range of values of I inc,0 is called hereafter the fully linear regime. Consequently, a simple experimental criterium for tuning I inc,0 to sufficiently small values for being in the fully linear regime is that the shape of the photodiode voltage signal shall not depend on I inc,0 (i.e. on the UT supply power). Moreover, in this regime the measured photodiode voltage signal can be directly compared to the predicted hump tip curvature κ(0, t). Notice also that the sensitivity of ∆P * to κ increases with |z|. Consequently, we have chosen z z 0 in order to satisfy the condition of sphericity of the mirror formed by the hump and to make ∆P * sensitive to the κ variations.

Experimental results

In order to test the validity of the analytical model described for all dynamic regimes, i.e. over a large range of values of ∆ c , with a special attention to the intermediate regime, for which no analytical prediction was formerly available, we performed surface deformation experiments using several liquids (water and Silicone oils) having different viscosities ranging from 1 cSt to 100 cSt. The density and the viscosity of the Silicone oils could be retrieved from the manufacturer (ABCR, Germany), while those of water could be found in [START_REF] Weast | Handbook of Chemistry and Physics[END_REF]. The sound speed has been measured in-house using the time of flight method. Since the surface tension is very sensitive to impurities, we measured the surface tension of each liquid sample immediately after the surface deformation experiment using the drop weight method and the correlation of [START_REF] Zhang | Formulation of the harkins-brown correction factor for dropvolume description[END_REF]. All the measurements have been performed in the temperature range T = 23 -26 • C. The physical properties of these liquids are presented in Table 1.

In figure 12 the time-dependent relative variation of the power of the reflected and filtered beam ∆P * measured in the fully linear regime is compared with the hump tip curvature evolution κ(0, t) predicted using the above presented analytical model (3.39). When implementing the radiation pressure excitation into the model, we took into account the additional contribution of the acoustic echoes impinging on the surface after their reflection on the UT surface with a delay of nτ , n being the number of reflections on the UT surface encountered by the echo. The amplitude of the echoes could be independently and precisely measured by visualizing their electric signature on the UT voltage signal as they reflect on the UT emitting surface. The associated pressure reflection coefficient was found to be equal to 0.45 for water and to decrease monotonously from 0.43 for 5 cSt to 0.32 for 100 cSt Silicone oil. Given these values, only the two first echoes noticeably contribute to the overall radiation pressure excitation. As shown in figure 12, a quantitative agreement between the hump tip curvature evolution and the time-dependent measured optical power, both normalized by their maximum, is obtained without any adjustable parameter for all the experiments, demonstrating the accuracy of our analytical model for each dynamic regime.

8. Validation of a numerical procedure for computing the spatio-temporal dynamics of a fluid-fluid interface

The aim of this section is to present in details a numerical procedure describing the spatio-temporal dynamics of a fluid-fluid interface and to validate it by comparing its predictions in the particular case of a liquid-gas interface to the analytical model presented above. Such a procedure is indeed required for predicting the behavior of forced interfaces in more complex situations such as free [START_REF] Zhang | Shape distorsion of a signle water droplet by laser-induced electrostriction[END_REF]) or sessile drops [START_REF] Chraïbi | Stretching and squeezing of sessile dielectric drops by the optical radiation pressure[END_REF]), meniscii, wetting films, for which the assumption of infinite liquid phases and of initial flat interface are not valid. This method solves a two-phase unsteady Stokes problem where the resolution of stress and velocity field is performed in both fluids. The Boundary Element Method (BEM) is used for this purpose. It is indeed very accurate for solving interfacial problems, as demonstrated in many references [START_REF] Sherwood | Breakup of fluid droplets in electric and magnetic fields[END_REF]; [START_REF] Manga | Low reynolds number motion of bubbles, drops and rigid spheres through fluid-fluid interfaces[END_REF]; B. [START_REF] Berkenbusch | Liquid interfaces in viscous straining flows: numerical studies of the selective withdrawal transition[END_REF]). Since the hydrodynamic problem is unsteady, The BEM is applied to Brinckman's equation, as presented hereafter.

Numerical implementation

We consider two fluid phases labeled i, i = 1, 2. The mass conservation and unsteady Stokes equations are for each phase, respectively:

∇ • v i = 0 ; i = 1, 2 (8.1) ρ i ∂v i ∂t = -∇P i + η i ∆v i ; i = 1, 2 (8.2)
For sake of generalization, we make the problem dimensionless considering v 0 , k -1 c , P i0 , τ 0 respectively as a reference velocity, length, pressure (in each fluid) and time such as:

v 0 = σ η 1 ; P i0 = η i v 0 k c ; τ 0 = 1 k c v 0 (8.3)
The temporal discretization of equation (8.2) follows a standard implicit approach,

∆ i v (t+∆t) i -v (t) i ∆t = -∇P (t+∆t) i + ∆v (t+∆t) i ; i = 1, 2 (8.4)
where

∆ 1 = σ ρ 1 ν 2 1 k c ; ∆ 2 = σ ρ 2 ν 2 2 k c ζ (8.5)
and ζ = η 2 η1 is the dynamic viscosity ratio. ∆ 1 , which is the relative magnitude of inertial effects to viscous effects in fluid 1, corresponds to its frequency parameter.

Defining the parameter α i and the body force b i in each phase as:

α i = ∆ i ∆t ; b i = α 2 i v (t) i (8.6)
we find the following Brinkman's equation, which can successfully be solved with the BEM [START_REF] Occhialini | Boundary intergral/spectral element approaches to the navier-stokes equations[END_REF]; [START_REF] Pozrikidis | A practical guide to boundary element methods with the software library BEMLIB[END_REF]).

α 2 i v (t+∆t) i = -∇P (t+∆t) i + ∆v (t+∆t) i + b i ; i = 1, 2 (8.7)
In addition to the bulk flow equations, we consider the stress balance condition at the interface:

T 1 • n 12 -ζT 2 • n 12 = (κ(r) -Bo c z + π rad (r, t))n 12 (8.8)
where T i = -p i I + 2D(v i ) and π rad (r, t) is a Gaussian excitation:

π rad (r, t) = π 0 e -r 2 2 H(t)H(δt -t). (8.9)
π 0 is a dimensionless number giving the relative magnitude of radiation to Laplace pressure, δt the dimensionless pulse duration.

The continuity of velocity at the interface together with a Lagrangian description of its motion give:

v 1 = v 2 (8.10) dx dt = v(x) (8.11)
Since the solutions to Brinkman's problem can be formulated in terms of Green's functions [START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF][START_REF] Pozrikidis | A practical guide to boundary element methods with the software library BEMLIB[END_REF])), we can rewrite the governing equations as a system of integral equations over the boundaries of the computational domain. Once boundary conditions on S I , S W 1 and S W 2 are used (see figure 13), and n e being a vector normal to the surface and pointing outward, the velocity field in the fluid bulk can be written as:

v i (x) = S I +S W i U • T i • n e dS y - S I +S W i n e • K • v i dS y + V i U • b i dV y i = 1, 2 (8.12)
Adding the latter equation written for fluid 1 to the same written for fluid 2 multiplied by ζ and using (8.8) and n = n 12 leads to the following equality:

1 + ζ 2 v(x) = S I U • n(κ(r y ) -Bo c z(r y ) + π rad (r y , t))dS y + (ζ -1) S I n • K • vdS y + S W 1 U • (T 1 • n)dS y -ζ S W 2 U • (T 2 • n)dS y + V1 U • b 1 dV y + ζ V2 U • b 2 dV y (8.13)
U and K are Brinkman's Green kernels for velocity and stress respectively and are given by ((Pozrikidis 1992)):

U ij = 1 8πd Aδ ij + B d i d j d 2 (8.14) K ijk = 1 8π - 2 d 3 (δ ij d k + δ kj d i )[e -D (D + 1) -B] - 2 d 3 δ ik d j [1 -B] -2 d i d j d k d 5 [5B -2e -D (D + 1)] (8.15)
where d = y -x, y(r y , z y ) is the integration point and D = αd. A and B are defined as:

A(D) = 2e -D 1 + 1 D + 1 D 2 - 2 D 2 (8.16) B(D) = -2e -D 1 + 3 D + 3 D 2 + 6 D 2 (8.17)
and A(0) = B(0) = 1. In (8.13), the first term in the right hand side describes the flow contribution due to interfacial tension, radiation pressure and gravity, whereas the second term accounts for the shear rates contrast on the interface. The latter term vanishes when ζ = 1. The third and fourth terms account for shear occurring along S W 1 and S W 2 as a result of the no-slip boundary condition. The two last terms account for the contribution of the body force. Unlike the other terms, the latter are evaluated over the bulk of fluids 1 and 2. The velocity along the interface as well as the stress exerted on all the boundaries S I , S W 1 and S W 2 are determined by solving the discrete form of this equation using a numerical procedure. This procedure requires first the reduction of surface boundaries S I , S W 1 and S W 2 and volumes V 1 and V 2 into line boundaries and surfaces, respectively. Due to this integral formulation and to the axial symmetry, the problem is reduced to a 1D problem. Line boundaries, shown in figure 13, need to be discretized with line mesh while surfaces S 1 and S 2 are discretized in quadrilaterals. The numerical integration on these quadrilaterals is performed using 2D Gauss quadratures [START_REF] Davis | Methods of numerical integration[END_REF]) and isoparametric transformations. The line mesh is made of constant boundary elements i.e. line segments with centered nodes. Both azimuthal and line integrations of (8.13) are performed using Gauss quadratures. The fluid-fluid interface is parameterized in terms of arc length and is approximated by local cubic splines, so that the curvature can be accurately computed.

The distribution and the number of points are adapted to the shape of the interface, so that the concentration of elements is larger in regions where the interface curvature gradient is large. The motion of the interface is followed using the kinematic condition (8.11) which is discretized using an Euler time scheme, such as:

r (t+∆t) = r (t) + v (t) r + v (t+∆t) r 2 ∆t (8.18) z (t+∆t) = z (t) + v (t) z + v (t+∆t) z 2 ∆t (8.19)
At the beginning of the first iteration, t = 0 and v i = 0 and therefore there is no bulk term b i . At the end of this iteration, the velocity field is solved at the interface but η 1 (Pa.s) ν 1 (cSt) ∆ 1 ∆ 2 0.054 54 5 17.9 0.077 80 2.5 12.5 0.17 170 0.5 5.7 0.55 550 0.05 1.76

Table 2. Values of the fluid properties used in the computations.

is also calculated at the interior of the domains on a N r × N z uniform grid with N r and N z adapted such that the velocity spatial variations is accurately described. For next iterations, the bulk term b i , which is always calculated from the previous velocity field, is interpolated using a four-points inverse distance weighting method with power parameter equal to one.

A typical computation begins with a flat interface at rest. The excitation is switched on at t = 0. The direction of deformation of the interface depends on the sign of π 0 . The time step is chosen to provide an accurate description of the temporal variation of velocity. Variations of the time step did not show significative change in the results.

Comparison with analytical predictions

7 We now compare the results given by the numerical resolution with the analytical model. Since this model is suited to a liquid-gas interface, we consider that phase 2 is air (ρ 2 = 1.3 kg/m 3 , η 2 = 2 10 -5 Pa.s, ν 2 = 15 cSt) and that phase 1 is a liquid (ρ 1 = 1000 kg/m 3 ) whose dynamic viscosity η 1 can be freely varied. Surface tension is chosen to be equal to σ = 70 mN/m and the characteristic wave number k c is fixed equal to 4712 m -1 in order to correspond to a water-air interface in the experimental configuration presented on §.

Table 2 summarizes the numerical values used in the computations. Figure 14 shows a comparison between the numerical prediction using the BEM and the analytical model. A quantitative agreement is observed concerning both the axial velocity at the hump tip and the temporal variation of its curvature in the frequency parameter range ∆ c 5 (since Bo 1, ∆ 1 ∆ c ). We have noticed that the mismatch between the numerical predictions and the analytical model increases as ∆ c increases. This can be explained by noting that our numerical method is better suited to flows where the viscosity has important to moderate effects. For purely inertial or weakly viscous regimes (∆ c 1), a resolution based on a potential flow formulation (i.e. Laplacian equation) and BEM would be more adequate. The observed quantitative agreement between the numerical simulation and the analytical model in the viscous and intermediate regimes demonstrates the accuracy of this numerical scheme.

Conclusion

To summarize, we developed an analytical model of the time dependent, small amplitude deformation of a free liquid surface caused by a spatially localized, axisymmetric, pulsed of continuous, acoustic or electromagnetic radiation pressure exerted on the surface. By exactly solving the unsteady Stokes equation, we have been able to predict the surface dynamics in all dynamic regimes, namely inertial, intermediate and strongly ) 2 , i.e. in all dynamic regimes, by comparing its prediction to experiments consisting in optically measuring the time dependent curvature of the tip of a hump created at a free liquid surface by the radiation pressure of an acoustic pulse. Finally, we have developed a numerical scheme simulating the behavior of a fluid-fluid interface submitted to a time-dependent radiation pressure, and we have showed its accuracy by comparing the numerical predictions with the analytical model in the intermediate and strongly damped regimes. Since radiation pressure based non-contacting techniques of measurement of liquid properties usually involve all the dynamic regimes, this new general solution of this initial-value problem of spatiotemporal dynamics of a liquid surface should constitute a key for the improvement of the accuracy and versatility of these techniques. Moreover, the analytical model presented in this paper can also be straightforwardly applied to the case of unsteady deformations of a free liquid surface caused by an axisymmetric electric field. Finally, note that the seemingly similar problem of the deformation of a liquid-liquid interface caused by radiation pressure, which is also of considerable practical interest [START_REF] Mitani | Measurement of ultralow interfacial tension with a laser interface manipulation technique[END_REF][START_REF] Williams | Non-contact techniques for measuring viscosity and surface tension information of a liquid[END_REF]), presents some additional difficulties [START_REF] Chandrasekhar | The character of the equilibrium of an incompressible heavy viscous fluid of variable density[END_REF]) which we plan to address using the same tools in the near future.

Figure 1 .

 1 Figure 1. Geometry and notations used in the model.

Figure 2 .

 2 Figure 2. (a) Variations of the normalized acoustic intensity distribution in the transducer focal plane k c r → φ(r). (b) Corresponding variation of its normalized Hankel transform k/k c → k 2 c φ(k). (c) Normalized electromagnetic intensity distribution at the beam waist. (b) Corresponding normalized Hankel transform.

∆Figure 3 .

 3 Figure 3. Variations versus ∆ of the real part (a) and the imaginary part (b) of the four roots X i , i = 1 -4 of (3.35). Variations versus ∆ of the real part (c) and the imaginary part (d) of the associated growth rates s i , i = 1 -4, as defined by (3.25). Upper inset of (a): Variation of ∆ versus k for a water free surface (black curve) and a 100 cSt Silicone oil (red curve). Other insets: same data as the main figure in logarithmic scales.

Figure 4 .

 4 Figure 4. (a) Variations versus ∆ of the ratio of the radial components 2X i /(1 + X 2 i ) and of the axial components 2/(1 + X 2 i ) of the irrotational part and of the viscous part of the flow associated with each root X i , i = 1 -4 of (3.35). (b) Variations versus ∆ of the ratio k/l the characteristic lengthscales of penetration of the irrotational part and of the viscous part of the flow associated with each root X i of (3.35), i = 1 -4.

  Bo c k2 . Thus, assuming Bo c 1 (which is the case in all experiments reported in the introduction), the characteristic time scale of variation of the surface deformation is Ω -1 c , where Ω c = σk 3 c ρ .

Figure 5 .

 5 Figure5. Upper curves: radius dependent curvature of the free surface of an inviscid liquid at various instants after a pulse of acoustic radiation pressure. Lower, dashed curve: radius dependence f (r) of the axisymmetric radiation pressure pulse exerted on the surface. Lower, solid curve: radius dependent curvature of the free surface at t = 30 µs, magnified four times for clarity purpose. Both latter curves have been shifted vertically.

Figure 6 .

 6 Figure 6. Dimensionless curvature of the hump tip κ(0, t)/κ c as function of time t for various values of the frequency parameter ∆ c .

Figure 7 .Figure 8 .Figure 9 .

 789 Figure 7. Dimensionless height of the hump tip h(0, t)/h c as function of time t for various values of the frequency parameter ∆ c .

Figure 10 .

 10 Figure 10. Experimental setup. UT: ultrasonic transducer. M: mirror. λ/4 : quarter-wave plate. BSC: beam splitter cube. GP: Glan Prism. ID: iris diaphragm. PD: photo-detector. cw laser: continuous wave laser.

Figure 11 .

 11 Figure11. Predicted variations versus the curvature of the hump tip of the relative variation of the optical power measured by the photodetector for several values of the distance between the laser beam waist and the liquid surface (z is positive if the laser beam waist is located in air, i.e. in front of the liquid surface). z 0 is the Rayleigh range of the laser beam.

Figure 12 .

 12 Figure 12. Comparison between the measured time-dependent power of the reflected and filtered beam ∆P * (t) that is proportional to the hump tip curvature (solid curves) and the hump tip curvature evolution κ(0, t) predicted using our analytical model (symbols): (a) water, (b) 5 cSt, (c) 10 cSt, (d) 20 cSt, (e) 50 cSt, (f) 100 cSt Silicone oil. All the signals have been normalized by their maximum. Note the progressive change of the horizontal scale.

Figure 13 .

 13 Figure 13. Typical configuration of the computational domain showing the two fluids and the liquid interface initially at rest.

Figure 14 .

 14 Figure 14. Temporal variation of interface tip curvature of a liquid-air interface for various values of the frequency parameter ∆ 1 ∆ c (symbols: numerical results, solid curve: theory). Pulse duration is 40 µs, π0 = 1 and Boc = (ρ 1 -ρ 2 )g σk 2 c = 0.006. Inset: Temporal variation of tip vertical velocity (symbols: numerical results, solid curve: theory).

Table 1 .

 1 Physical properties of the liquid samples used for the experimental test of the analytical model.

		viscosity density sound speed ν ρ c (cSt) (kg.m -3 ) (m.s -1 )	surface tension σ (mN.m -1 ) (drop weight method)	∆ c
	water	1	1000	1490	72	15, 000
	Silicone oil	5	918	968	19.6	100
	Silicone oil	10	935	980	19.8	30
	Silicone oil	20	950	995	19.9	6
	Silicone oil	50	970	997	20.2	1
	Silicone oil	100	966	1000	20.2	0.3