Franc ¸ois-Régis Jasnot

Louis-Anne De Vaulchier

Yves Guldner

Gérald Bastard

Angela Vasanelli

Christophe Manquest

Carlo Sirtori

Mattias Beck

Jérôme Faist

Direct surface cyclotron resonance terahertz emission from a quantum cascade structure

A strong magnetic field applied along the growth direction of a semiconductor quantum well gives rise to a spectrum of discrete energy states, the Landau levels. By combining quantum engineering of a quantum cascade structure with a static magnetic field, we can selectively inject electrons into the excited Landau level of a quantum well and realize a tunable surface emitting device based on cyclotron emission. By applying the appropriate magnetic field between 0 and 12 T, we demonstrate emission from a single device over a wide range of frequencies (1-2 THz and 3-5 THz).

In 1960 well before the advent of diode lasers, it was proposed that semiconductors immerse in a static magnetic field were likely candidates for use as quantum amplifier and oscillators 1 . A follow up on this idea was implemented in the 80's on a highly purified bulk germanium lightly doped with holes 2,3 . In this system laser action in the terahertz (THz) region was obtained, at cryogenic temperatures, by cross combining intense magnetic and electric fields. More recently, semiconductor THz lasers have been demonstrated using quantum cascade laser concepts 4,[START_REF] Liu | Intersubband Transitions in Quantum Wells : Physics and Device Applications I[END_REF] . Owing to band structure engineering and the use of guided optics a much higher degree of control on the spectral and spatial properties have been obtained in these last devices. Using quantum cascade technology distributed feedback lasers 6 , photonic crystal devices 7 and external cavities have been demonstrated 8,9 . However, an inherent property of these lasers is that the spectral position of the optical gain is fixed, at a given current, and therefore to realise wavelength tunable devices it is necessary to employ mechanical systems.

In this paper we propose to combine the quantum engineering of a quantum cascade structure (QCS) with a static magnetic field to realise frequency tunable surface emitting electroluminescent devices. The specific strategy is the use of a resonant injector which allows to populate selectively a given Landau level (LL) of a quantum well. Besides, the fundamental concept of the QCS allows to magnify even further the inter-LL emission signal by increasing the number of periods, the extraction miniband of a period being the injection miniband of the next one. In our device electrons are injected into an excited state of a transition between two LLs arising from the extra confinement imposed by the magnetic field. It is in fact known that a magnetic field parallel to the growth axis of a quantum well structure breaks the two-dimensional parabolic energy dispersion of each subband ε i (k) into a ladder of discrete LLs with an energy ε i,j = E i + (j + 1/2) ω c and separated by the cyclotron energy ω c = eB/m * , where i is the subband index, E i the energy of the subband edge at zero magnetic field, j the LL integer index, B the magnetic field and m * the energy dependent electron effective mass. Experiments with magnetic field have been widely used to evaluate the different contributions of scattering mechanisms in com-plex QCSs [10][11][12][13][14][15][16][17][18] .

To achieve electrical injection into an excited LL we have conceived a QCS in which electrons are injected into the first excited subband |2 of a large quantum well. Under the influence of a magnetic field the excited LLs arising from the first subband, |1, j , get in energy proximity with the first LL of the subband |2, 0 . It has been theoretically and experimentally demonstrated that in the vicinity of crossing LLs which belong to different subbands, a disorder activated transition becomes significant and allows a redistribution of the electron population between them 19 . By exploiting this mechanism electrons are therefore transferred, via an elastic scattering, from |2, 0 to |1, j and subsequently generate photons by spontaneous emission between |1, j and |1, j -1 . The optical dipole between LLs is perpendicular to the growth axis Oz and thus photons can be emitted through the surface of the sample with a polarization ε x contrary to intersubband transition where light has a polarization ε z and is emitted through the side of the device. [START_REF] Liu | Intersubband Transitions in Quantum Wells : Physics and Device Applications I[END_REF] The sample used is a GaAs/Al 0.15 Ga 0.75 As THz quantum cascade structure, emitting at 4.3 THz. The active region is made of 80 identical periods, based on a main quantum well, where the radiative transition takes place between levels |2 and |1 (E 2→1 ≈ 18 meV → 4.3 THz), and an injection/extraction miniband (Figure 1(a)). The detailed succession of wells and barriers (in Angstrom) for a single period is the following: 44 / 266 / 40 / 164 / 22 / 154 / 22 / 148 / 24 / 146 / 26 / 144 / 28 / 140 / 34 / 134 (Al 0.15 Ga 0.75 As barriers in bold). The 148 Å and 146 Å GaAs layers are Si-doped 2 × 10 16 cm -3 . Two devices, specially conceived to observe surface emission, have been processed : a mesa and a ridge. A schematic of the mesa is shown in the inset of Figure 1(b). The top of the mesa is a 300 µm × 300 µm GaAs free surface surrounded by an edge of gold of 50 µm width to realise electrical connection. The gold layer is separated from the mesa by a 0.5 µm layer of silicon nitride preventing electrical shortcircuit. As all the sides of the device are covered by gold, only surface emission is possible from the device. A 300 µm wide and 3.6 mm long ridge was also processed. The top surface of this device consists of two 50 µm wide bands of gold separated by a 200 µm wide GaAs free surface. This geometry allows the observation of a luminescence signal from both the surface and the facets.

Electroluminescence spectra of the two devices at zero magnetic field are shown in Figure 1(b). The spectra are measured for the same current density with a Bruker Vertex 80V FTIR in the step scan mode on a liquid-helium cooled Sibolometer, the signal being detected with a lock-in amplifier. In our setup, the detector is turned toward the surface of the sample. The spectrum measured from the ridge consists of a narrow intersubband electroluminescence peak centered at an energy of ∼ 18 meV in agreement with the calculated energy E 2→1 . This indicates that, although the detector is turned towards the surface of the sample, the ε z polarized light emitted through the facet is nonetheless detected. This is due to the diffraction of the light on the side of the ridge and to its multi-reflection on the tube around the sample. In addition, a weak broad peak is highlighted between 4 and 14 meV and is attributed to the black body of the sample. By using the Wien law, we find a temperature of ∼ 17 K. On the contrary, no peak is visible on the mesa spectrum. No black body radiation is observed either because the current flowing through the device is lower. The size of the sample is indeed reduced by a factor of 7. These two spectra demonstrate the efficiency of our mesa processing to block the ε z polarized emission.

The mesa device was mounted inside an insert at the center of a superconducting coil capable of producing fields up to 12 T such that the magnetic field lines are perpendicular to the plane of the quantum wells. Figure 2(a) shows the evolution of the voltage as a function of the magnetic field at fixed current density. The voltage exhibits oscillations as a function of B superimposed on a continuous increasing background which is attributed to the magnetoresistance of the contact 15,17 . Three main minima are located at B = 3.5 T, 5.2 T and 10.5 T. An extra minimum, less pronounced, is present at B = 7.2 T. These features are a consequence of the electron lifetime modulation on the level |2, 0 of the thick quantum well. Figure 2 Two elastic mechanisms can occur in such situation: interface roughness or electron-electron (Auger like) scattering. The two mechanisms can occur at B = 3.5, 5.2 and 10.5 T whereas at B = 7.2 T only electron-electron scattering can occur. Interface roughness has been widely reported as the main mechanism in GaAs/AlGaAs THz QCS 14,16 . Signature of electron-electron scattering has been also reported in transport 21,22 and laser power measurements 23 .

Experimental spectra of the mesa at different magnetic fields are presented in Figure 3. For the sake of clarity the residual background and the black body appearing at high magnetic field have been removed. Two main series of peaks can be observed, one around 8 meV and another one around 18 meV. These two characteristic energies correspond to the cyclotron energies in the vicinity of 5.2 and 10.5 T respectively. It has to be noted that around 6 and 12 meV, which corresponds to 3.5 and 7.2 T respectively, no such series could be measured. This is related to the weak features of the transport measurement minima at these magnetic fields as reported in Figure 2(a). It can also be seen that several spectra show two or three peaks rather than one only. In order to analyse the data we have represented in Figure 4 the peak energy of each emission spectra in Figure 3 as a function of magnetic field. The data reported in red squares correspond to the cyclotron emission and follow the ω c law (m * = 0.067m 0 for GaAs) indicated in dashed line. As expected red points only show up in the grey areas where LLs are populated with electrons (see caption).

In addition to the cyclotron emission peaks, several others are measured. They are reported in black squares and two series can be identified; one at a roughly constant energy of 18 meV between 5 and 9 T, and another one, in the similar B range, 2 meV below. These two characteristic energies are close to the intersubband transitions E 2→1 and E inj→1 respectively. As shown in the first part of this paper, light emission with a polarization ε z is blocked in our device. Therefore, even if these peaks correspond to the intersubband energy, they cannot originate from |2, 0 → |1, 0 emission, that is polarized perpendicular to the surface. It could originate from carrier localization due to disorder 23 . Also, it is worth noting that theoretical considerations have demonstrated that normally forbidden transitions in ideal structures could be activated by disorder. Specifically, Regnault et al. have shown that the disorder induced coupling between LLs could be responsible for transitions in the wrong polarization. This mechanism could explain the two series at 18 and 16 meV as electron-electron interaction can couple LLs. This observation implies that the widely used picture of energy-degenerate LLs, crossing without interaction has to be revisited.

In conclusion we report a magnetic field tunable surface emitting device in the 1 -2 and 3 -5 THz range based on cyclotron emission. The study of the CR emission allows the identification of the dominant scattering mechanisms occurring in the QCS. Our experiments show that QCS are good candidates for realizing Landau-level laser 24 or tunable THz VCSEL. In order to obtain a population inversion between LLs |j and |j + 1 , non equidistant LLs are necessary in order to dissymetrise the absorption and emission transition energies between LLs. This situation could occur in our sample by taking advantage of the disorder inducing an inhomogeneous broadening of the LLs at large magnetic field. Additionally, this situation could occur in QCS based on low gap materials, such as InAs, in which the non-parabolicity could be exploited.

  FIG. 1: (a) Schematic conduction band diagram and squared wave functions of one QCS period under an applied electric field. The radiative transition occurs between the first excited state |2 and the fundamental state |1 of the thick quantum well. Electrons are injected into |2 via the injector level |inj . The grey shaded regions correspond to the injection and extraction minibands. (b) Zero magnetic field electroluminescence spectra at a current density of 40 A.cm -2 for the mesa (red) and the ridge (black) devices. The temperature is 4.5 K. Inset : schematic of the QCS with the specific processing to block intersubband emission with a polarization εz.

FIG

  FIG. 2: (a) Voltage as a function of the magnetic field at constant current density for the mesa device. (b) Plot of the LLs |2, 0 and |1, j . Two elastic scattering mechanisms are displayed : electron-electron interaction (red arrows) and interface roughness (green circles).

  (b) shows the LLs |2, 0 and |1, j as a function of the magnetic field. Dashed vertical lines are used to indicate the values of the magnetic field for which |2, 0 is in resonance with one of the |1, j and/or halfway between two LLs of |1 .

FIG. 3 :

 3 FIG.3: Emission spectra for different magnetic field from B = 4.7 to 12.0 T at fixed current density J = 40 A.cm -2 . The resolution is 1.2 meV and the temperature is 4.5 K.

FIG. 4 :

 4 FIG. 4: Emission peak energy at different magnetic field and fixed current density. Error bars are due to the experimental resolution of our spectra (1.2 meV). Red squares correspond to the cyclotron emission peaks. The black squares correspond to disorder activated forbidden transitions. The grey areas are proportional to the number of injected electrons in the LLs as deduced from the experimental transport data of Figure 2(a).

  B. Lax, Cyclotron resonance and impurity levels in semiconductors : Proceedings of the International Symposium on Quantum Electronics, edited by C. H. Townes (Columbia Univ. Press, New York, 1960), p. 428 2 Y. L. Ivanov and Y. V. Vasil'ev, Sov. Tech. Phys. Lett. 9, 264 (1983) 3 K. Unterrainer, C. Kremser, E. Gornik, C. R. Pidgeon, Yu. L. Ivanov, and E. E. Haller, Phys. Rev. Lett. 64, 2278 (1990)

1


E Maitre

C Misbah

P Peyla

A Raoult

Comparison between advected-field and level-set methods in the study of vesicle dynamics

Keywords: numbers: 87.16.D-83.50.Ha 87.17.Jj 83.80.Lz 87.19.rh Keywords: level-set, phase-field, vesicle dynamics

Phospholipidic membranes and vesicles constitute a basic element in real biological functions. Vesicles are viewed as a model system to mimic basic viscoelastic behaviors of some cells, like red blood cells. Phase field and level-set models are powerful tools to tackle dynamics of membranes and their coupling to the flow. These two methods are somewhat similar, but to date no bridge between them has been made. This is a first focus of this paper, where we show how the phasefield methods developed in [1-3] for immersed vesicles could be considered as a level-set method for a particular strain-stress relationship. The main conclusion is that the two methods share several common features and we shall provide the correspondence between the two methods. Furthermore, a constitutive viscoelastic law is derived for the composite fluid: the ambient fluid and the membranes. We present two different approaches to deal with the membrane local incompressibility, and point out differences. Some numerical results following from the level-set approach are presented.

I. INTRODUCTION

Vesicles are closed membranes which are suspended in an aqueous solution. The membrane is made of a bilayer of phospholipid molecules. These molecules have a polar hydrophilic head which points towards the solvent, and hydrophobic tails which point towards the interior of the membrane. Both at room and physiological temperatures the bilayer is a two dimensional incompressible fluid.

Human red blood cells (RBCs) are among the simplest animal cells. They are made, like vesicles, of phospholipid bilayer, plus a protein network (called spectrin), known also under the generic name of cytoskeleton. The internal solution of RBC is made of a hemoglobin solution (a newtonian fluid). The RBC is devoid of a nucleus and organella, and this is why it may be viewed as a simple cell. Its main function consists in oxygen supply to tissues. It is hoped that vesicles may represent a simplistic starting point to understand viscoelastic properties, dynamics and rheology of bio-fluids, such as blood.

The problem of vesicles at equilibrium (equilibrium shapes) is now fairly understood [4]. The study of vesicles under non-equilibrium conditions has constituted a major focus in recent years: [1,2,. Vesicles under flow have revealed several fascinating non-equilibrium behaviors (tank-treading, tumbling, vacillating-breathing-aka swinging, trembling, and so on). RBCs exhibit also similar kinds of dynamics (see recent review [31]). This research activity knows nowadays an increasing interest in the field of physics, mechanics, applied mathematics, engineering science, and so on. This field has known during the past decade analytical, numerical and experimental progresses (see review [31]). Even at zero Reynolds number, where the hydrodynamics equations are linear (the Navier-Stokes equations reduce to the Stokes equations), the problem remains highly nontrivial due to the fluid/membrane coupling. This coupling triggers nonlinearity (inherent to any moving boundary problem, even if the bulk equations are linear), and nonlocality (a disturbance of the membrane shape at some point is felt by other regions of the membranes due to the fact that hydrodynamics have long range effects; very much like Coulomb interactions in electrostatics). In mathematical terms, the velocity field in the fluid can be integrated out in favor of a closed nonlinear integro-differential equations. This is the so-called boundary integral formulation based on the use of the Greens function technique [START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF]. This method has been used for vesicles by several groups: [5,8,[START_REF] Veerapaneni | [END_REF]. This method, despite its efficiency and precision, has some limitations. For example, in a situation where the internal (or external) fluid is non newtonian, and/or if the non zero Reynolds number limit (a situation encountered for blood flows in veins and arteries, for example) has to be treated, then the resulting nonlinearities of the fluid equations clearly rule out the use of the Greens function techniques. We have thus to resort to other methods.

Phase-field and level-set approaches can, in principle, handle the above problems quite naturally. Phase field methods have known some popularity in the scientific community, and have been used in various topics where interfacial phenomena are present (crystal growth, fluid/fluid interfaces, and so on). More recently, a phase-field method has been adapted to vesicles [1,3]. Later, other phase-field formulations have been presented for the same system [34][35][36].

In classical interfacial problems (crystal growth, vapor/liquid phase transitions, fracture...), the total length of the interface can grow (or shrink) without bound. In contrast, for membranes a new important ingredient comes into play: the membrane is a two dimensional incompressible fluid, and thus its area does not change in the course of time. More precisely, not only should the total area remain constant in the course of time, but also the area must be conserved locally. This question is not easy to deal with in practice (a reason why only global area conservation have been imposed in some phase-field models as in [34,35]). Like the phase-field, the level-set approach treats the problem in the same spirit, but the precise formulations differ in each case, as we shall see here. How does the level-set formulation compare with the phase-field one is a question which has not yet been treated in the literature, and especially regarding the problem of vesicles and membranes. This paper addresses this question as a first focus.

For the sake of completeness, it is worthwhile to cite other alternative methods which have been adapted for the study of biological membranes. These are (i) techniques based on dynamically triangulated models [37] or particlebased mesoscale solvent, multi-particle collision dynamics [38], or their combination [12], (iv) the so-called immersed boundary method [39], used to model some features of red blood cells [40], and (v) Lattice-Boltzmann methods [41].

In this paper we shall draw a parallel between the phase-field and level-set methods. It will be shown that, despite their apparent differences, the two methods have a quite number of links. The level-set that we shall present here will be compared to the phase-field models previously introduced in two and three dimensions [1,3]. We shall make a bridge between the two methods. We shall also give here a derivation of the postulated equation in [1] and [3] regarding the treatment of the local membrane incompressibility. Finally, we shall show that the membrane force derived therein can be written as a divergence of some rank two tensor, and this will allow us to write a closed form for the constitutive law of the fluid/membrane system. This law exhibits a viscoelastic behavior. We shall then briefly discuss the numerical method used to solve the level-set equations. Some numerical results will be presented.

The scheme of the paper is as follows. In section II we present some essential preliminaries. In section III we recall briefly the phase-field equations for vesicles (or membranes). In section IV we introduce the level-set formulation. Both formulations are Eulerian methods relying on an auxiliary field to capture the interface. We adopt the notation φ for the phase-field and ϕ for the level-set function. We then discuss the numerical scheme used for the solution of the level-set model, and present some numerical results. In section V we shall make the bridge between the two methods (level-set and phase-field). We shall show in that section how one can write down a constitutive law for the composite fluid (fluid+membrane). Section VI is devoted to conclusion and discussion of some future research directions. Some technical details are relegated into annexes.

II. SOME PRELIMINARIES

Both phase-field and level-set methods have the common feature that they define the interface in an implicit manner; that is the phase-field also adopts the level-set notion. A significant difference exists, however. In a phasefield approach the phase-field values in the two fluids are forced to remain parallel (or have constant values away from the interface) by imposing a double well potential, while in a level-set approach the (color) function is just advected by the flow. Both methods have therefore potential advantages and drawbacks. Enforcing the phase-field function to given values and controlling the interface width by appropriate potential brings a very stable interface description. However this might introduce volume loss and additional terms in the potential are usually added to control this leak. The philosophy of phase-field in this context is therefore to give an Eulerian description of the sharp interface by a diffuse interface of constant and controllable width. The phase-field is a function φ which is only used to give geometrical information on the interface. By contrast, in the level-set method the auxiliary field ϕ is merely advected by the velocity field of the continuous medium, starting initially from a signed distance to the interface. There is no potential to keep the field stuck to a distance function (further details will be given below), and as the diffuse interface is typically considered to lie between two given values of ϕ, e.g. -ε and +ε, the interface width varies. The gain to move the level-set function by advection is that some mechanical information is recorded in ϕ. For example |∇ϕ| records the area change of interface. This makes it therefore possible to write elastic energy merely in terms of the auxiliary field and to obtain an attractive complex fluid model of the fluid-structure coupling. The non-constant width of interface has however to be taken care of. But this is made as post-processing each time the level-set function is needed to compute the distance to interface.

III. PHASE-FIELD FORMULATION

In this section let us recall the main results of the phase-field model introduced earlier [1][2][3]. In that model the vesicle is described by an advected field φ which goes smoothly from -1 to +1 while crossing the membrane. Let -1 represent the interior of the vesicle and +1 its exterior, the membrane being localized by the zero level-set of φ. In order to ensure a constant width for the interface, one considers the following functional

E intrinsic [φ] = 1 4 (1 -φ 2 ) 2 + ε 2 2 |∇φ| 2 dx (1)
when minimized in the special case of a flat interface leads to the interfacial profile given by:

φ(r) = tanh r ε √ 2
The above energy is called intrinsic in the sense that its role is only to define the boundary, and it should not affect the physics (like the forces acting on the membrane). The physical forces are accounted for by the introduction of the configurational energy which is given by

E conf ig = κ 2 c(φ) 2 |∇φ| 2 dx + ξ |∇φ| 2 dx (2)
In this expression the first term stands for the curvature energy with modulus κ (which is a phase-field expression of the Helfrich energy, c(φ) = div ∇φ |∇φ| being the mean curvature), while the second one is a penalization of the membrane local length variation, ξ being therefore a local (i.e. it depends on the given position and on time) Lagrange multiplier (we do not use ζ as in [1] and [3], since this symbol will denote a cut-off function thereafter). It was shown in [42] that the asymptotic limit of the above phase field equations recovers the sharp interface formulation. Other formulations of the Helfrich energy have been introduced (see [43] and references therein) that give the same sharp interface limit. Note that |∇φ| is nothing but a Dirac-like function whose role is just to localize the action of the curvature term c at the interface. Moreover a thermodynamically consistent picture of these equations has been studied in [44].

It is classical that the normal to the interface and projector on its tangent plane are given by

n = ∇φ |∇φ| , P = I -n ⊗ n = I - ∇φ ⊗ ∇φ |∇φ| 2 ,
where I is the identity tensor, while the curvature field is expressed as c = -div n. Differentiating the above energy with respect to the φ gives an explicit expression of the external force

F conf ig = -κ c 2 (c 2 -4G) + (P ∇) 2 c n + ξcn + P ∇ξ |∇φ| 2 (3)
where G = det((t 1 • ∇)n, (t 2 • ∇)n, n), is the Gaussian curvature and (t 1 , t 2 , n) an orthonormal trihedron. The nonstationnary Stokes equations are used to find the velocity field:

ε u u t -div η(∇u + ∇u t ) -∇p = F conf ig , div u = 0 (4)
where ε u is a relaxation parameter (taken small enough to mimic the zero Reynolds number limit) and η = η out (1 + φ)/2 + η in (1 -φ)/2 the smeared viscosity (that accounts for a viscosity contrast between the interior and exterior of the vesicle; η in and η out are the viscosities of the internal and external fluids). As the membrane is simply advected by flow, a simple transport equation with velocity u (i.e. dφ/dt = 0 where d/dt is the material derivative) should be used to find φ. For stability reasons, and to guarantee that the advected field function minimizes the intrinsic energy, and thus it represents the interface in the course of time, the phase-field is taken to obey the equation dφ/dt = -ε φ δE intrinsic /δφ, where δ/δφ is the functional derivative and ε φ a kinetic constant fixing the time scale.

Using (1) one obtains the following equation:

φ t + u • ∇φ = ε φ (φ(1 -φ 2 ) + ε 2 (∆φ + c|∇φ|)) (5)
where the term cε 2 |∇φ| has been added by "hand" in order to cancel the wall free energy of the membrane associated with ε 2 |∇φ| 2 (which is known to lead to a surface tension-like term in the force). Since the membrane does not exchange matter with its surrounding environment (unlike drops where molecules from bulk can migrate to surface and vice versa, leading to surface variation), its surface energy is zero. The added terms guarantees the absence of a surface tension [1,3] (recently another method to deal with this problem has been suggested [44]). Since the membrane is locally incompressible, one has to impose that the surface projected divergence of the velocity field must be zero. In the sharp interface picture [8,45] the local area incompressibility can be handled by introducing a space and time dependent Lagrange multiplier. This amounts to writing the contribution of the membrane energy related to the local membrane incompressibility condition as

E inc = ζ(r m , t)dA (6)
where the integration is performed over the vesicle membrane, and r m is the vector position of the membrane. The Lagrange multiplier ζ is then determined by imposing

(I -n ⊗ n) : ∇u = 0. (7)
The above expression is nothing but the divergence along the membrane (note that I -n ⊗ n is the projector). The Lagrange multiplier ξ does not appear in the above equation which is the associated constraint, just like the pressure does not appear in the solenoidal constraint in incompressible hydrodynamics. However, like the pressure, we show in the following that ξ appears in the other equations of the model (it couples to the velocity field). In the phase-field spirit [1,3] the idea is to define ξ everywhere (and not only along the membrane) but confining its action to the membrane region. We then write

E inc = ξ(r, t)|∇φ|dx ( 8
)
where dx is the volume element of the total domain. Because φ(r) ∼ tanh(r/ √ 2 ), it is clear that |∇φ| is a Dirac-like function of width . This implies that the energy acts in the membrane region only, as it should be.

In [1,3] the tension field was postulated to obey the following equation (apart from the Laplacian term which was included in [3] for some numerical regularization)

dξ dt ≡ ξ t + u • ∇ξ = T (I -n ⊗ n) : ∇u (9)
T is a tension-like parameter which is chosen large enough, so that we expect this to enforce (In ⊗ n) : ∇u 0, that is a local membrane quasi-incompressibility. This field was built on the basis of intuition and knowledge of the sharp interface problem [8,45]. It was shown in [3] by asymptotic techniques that in the limit → 0 the above equation (9) recovers the sharp interface limit (free surface divergence of the velocity field). In [1,3] ξ is said to be proportional to the local extension of the membrane, and indeed the right hand side of ( 9) is a measure of the surface variation, as we shall see. It will be shown in the next section how equation ( 9) can be derived.

IV. LEVEL-SET FORMULATION

A. Introduction

Let us now introduce the level-set approach. In the level-set formulation [START_REF] Osher | Level set methods and Dynamic Implicit Surfaces[END_REF], one still introduces a function ϕ which is negative inside the membrane and positive outside. This function is not constrained to take values between -1 and +1. Rather, it is initially the signed distance to the interface, and is then advected by a transport equation

ϕ t + u • ∇ϕ = 0 (10)
It was shown in [START_REF] Cottet | [END_REF]48] that when u is divergence free, |∇ϕ| records the stretching of the interface {ϕ = 0}; this opened the way to a complete formulation in terms of ϕ of the membrane forces. In order to localize the interface, one introduces a cut-off function ζ, i.e. a non-negative function with compact support included into [-1, 1], of unit mass such that 1 ε ζ( r ε ) converges to the Dirac mass when ε goes to zero. We used in our simulations the following expression

ζ(r) = 1 2 (1 + cos(πr)) on [-1, 1], ζ(r) = 0 elsewhere. This function is of unit mass. Considering ζ( r ε )
, this gives a function with support in [-ε, ε] which has a mass of ε as readily seen by integration. This is why the scaling by 1 ε is necessary to maintain a unit mass and ensure the convergence of 1 ε ζ( r ε ) to the Dirac mass. By composition with the level-set function, we obtain

1 ε ζ( ϕ ε ) whose support is localized in the strip -ε < ϕ < ε. If ϕ is a distance function 1 ε ζ( ϕ ε )
converges to the Dirac mass on the curve ϕ = 0 when ε → 0. In general 1 ε ζ( ϕ ε )|∇ϕ| converges to this Dirac mass as ε → 0.

Let us first consider the case where the membrane energy would depend on |∇ϕ| only (we shall see that |∇ϕ| is a measure of local surface variation). We shall call the corresponding energy elastic energy (in contrast to curvature or bending energy). A vesicle membrane is inextensible. In our spirit we shall allow for a finite, albeit quite small, area variation. If we have in mind a capsule, for example (see review [49]), then the membrane can undergo a certain local area variation. Therefore our energy may be used to various cases. We define the energy in the entire space provided that the density energy is multiplied by the appropriate localized function introduced above. We then have

E m (ϕ) = Ω E(|∇ϕ|) 1 ε ζ( ϕ ε )dx (11)
and differentiating (i.e. taking functional derivative with respect to ϕ) it along (10) leads to the following elastic force:

F m = ∇[E (|∇ϕ|)] -div E (|∇ϕ|) ∇ϕ |∇ϕ| ∇ϕ |∇ϕ| |∇ϕ| 1 ε ζ( ϕ ε ) (12)
Note that E stands for a constitutive law for the membrane. A trivial example is a linear law E (r) = Λ(r -1). In [50] it was shown that E (r) = Λ max(r -1, 0), is a more appropriate choice, since it can be obtained from a variational asymptotic derivation. Λ is a parameter. For vesicles Λ is taken quite large in order to enforce quasi-inextensibility.

A prior study [START_REF] Cottet | [END_REF]48] has already considered these membrane elastic forces, but not the curvature energy, which we consider here. We set

E c (ϕ) = Ω G(c(ϕ))|∇ϕ| 1 ε ζ( ϕ ε ) ( 13
)
where G is again a constitutive law for the curvature energy. A standard choice, which is compatible with [1] and [3], is G(r) = κ 2 r 2 , which is nothing but the Helfrich curvature energy density. Recall that c is the divergence of the normal vector, and hence c(ϕ) = div ∇ϕ |∇ϕ| (note that we have the +, while in [1,3] the opposite was chosen), which is thus positive for convex vesicles. Our strategy to derive the force from the energy is the following: we compute the time derivative of the energy and equal it to minus the power of the curvature force (this is sometimes called the virtual power -note that a direct functional derivative can be used as well, as in [3]):

d dt E c (ϕ) = - Ω F c • udx
From differential calculus and using the transport equation (10) solved by ϕ, we have

d dt E c (ϕ) = dE c (ϕ)(ϕ t ) = dE c (ϕ)(-u • ∇ϕ)
where dE c (ϕ)(δ) means the differential of E c at point ϕ applied to the increment δ. Therefore we only need to compute the differential of curvature energy and apply it to the increment -u • ∇ϕ. We show in annex I that

dE c (ϕ)(δ) = Ω div -G(c(ϕ)) ∇ϕ |∇ϕ| + 1 |∇ϕ| P ∇ϕ ⊥ (∇[|∇ϕ|G (c(ϕ))]) 1 ε ζ( ϕ ε )δdx
which by identification leads to

F c = div -G(c(ϕ)) ∇ϕ |∇ϕ| + 1 |∇ϕ| P ∇ϕ ⊥ (∇[|∇ϕ|G (c(ϕ))]) 1 ε ζ( ϕ ε )∇ϕ (14)
The equations for u in the level-set formulation is similar to the phase-field formulation (see ( 4)), with F m + F c as source term. In fact we rather use the full Navier-Stokes equations (i.e. by including u.∇u) but this induces non significant differences as long as the Reynolds number is small enough. Note that in this formulation there is no need to introduce a Lagrange multiplier and to postulate some corresponding evolution equation (as presented in the last section), since the stretching is encoded the ϕ variable. Note also that this implies (unlike in the case of the phase-field approach), that we have to solve a transport equation for ϕ, since it is under this evolution that |∇ϕ| records the membrane stretching (see next section and section V B). Before discussing the bridges between the phase-field and level-set methods, we shall first present the numerical method used for the level-set method along with some numerical results.

B. Numerical procedure

The level-set model amounts to solve the following set of equations for (u, ϕ):

ρ(u t + u • ∇u) -div(2ηD(u)) + ∇p = F m + F c (15) div u = 0 (16) ϕ t + u • ∇ϕ = 0 (17)
where D(u) = 1 2 (∇u + ∇u T ), with appropriate initial and boundary conditions. As there is no additional term in the transport equation, ϕ does not remain a distance function for t > 0. To see this, we use the fact that distance function are functions with unit gradient modulus. Differentiating (10), one easily find that ∇ϕ verifies

∇ϕ t + u • ∇(∇ϕ) + ∇u T ∇ϕ = 0.
Taking the scalar product with ∇ϕ leads to

1 2 (|∇ϕ| 2 ) t + 1 2 u • ∇(|∇ϕ| 2 ) = -∇ϕ T ∇u T ∇ϕ
Next we remark that ∇ϕ T ∇u T ∇ϕ = ∇ϕ T D(u)∇ϕ = (∇ϕ ⊗ ∇ϕ) : D(u) so that expanding the derivatives of squared norms leads to for all time. This means that the variation of u along the direction normal to the interface is zero, which is not true in general. As a consequence the support of the smeared delta function 1 ε ζ( ϕ ε ) will vary, making the smeared membrane width locally proportional to the inverse of extension. To circumvent this effect, one uses the following trick: we use ϕ |∇ϕ| as an approximation of the distance function to the interface, thus replacing

|∇ϕ| t + u • ∇|∇ϕ| = -|∇ϕ| ∇ϕ ⊗ ∇ϕ
|∇ϕ| 1 ε ζ( ϕ ε ) by 1 η ζ( ϕ |∇ϕ|η ). Indeed ϕ |∇ϕ| is a distance at first order in a neighborhood of interface (ϕ = 0) since ∇( ϕ |∇ϕ| ) = ∇ϕ |∇ϕ| + O(ϕ), thus ∇( ϕ |∇ϕ| ) = 1 + O(ϕ)
. This order one approximation could seem brutal, but several numerical evidences [48] proved its efficiency and precision in comparison with the usual reinitialization step used in level-set methods by e.g. [START_REF] Osher | Level set methods and Dynamic Implicit Surfaces[END_REF]. Especially in terms of volume conservation, it avoids redistancing usually involved in level-set methods and acknowledged to be responsible for volume leaks.

The numerical scheme used is a Chorin projection method on a MAC mesh, which enforces the exact divergence free condition, and therefore the volume conservation, at the discrete level. This is of high importance for our problem since change in the volume affect the shape of the minimizer of the curvature energy.

C. Dimensionless parameters

The level-set main equation is given by

ρ(u t + u • ∇u) -div(ηD(u)) + ∇p = ∇[E (|∇ϕ|)] -div E (|∇ϕ|) ∇ϕ |∇ϕ| ∇ϕ |∇ϕ| + div -G(c(ϕ)) ∇ϕ |∇ϕ| + 1 |∇ϕ| P ∇ϕ ⊥ (∇[|∇ϕ|G (c(ϕ))]) ∇ϕ |∇ϕ| |∇ϕ| 1 ε ζ ϕ ε Let L, U ,
u t = U 2 L u t , ∇u = U L ∇ u , D(u) = U L D (u ), div(ηD(u)) = U η ref L 2 div (rD (u )), ∇p = η ref U L 2 ∇ p , ∇ϕ = ∇ϕ , c(ϕ) = 1 L c (ϕ )
In dimensionless variables (dropping the ), and for the particular case E (r) = Λ(r -1) and G(r) = κ 2 r 2 we get:

Reρ(u t + u • ∇u) -div(rD(u)) + ∇p = 1 W e ∇[E (|∇ϕ|)] -div E (|∇ϕ|) ∇ϕ |∇ϕ| ∇ϕ |∇ϕ| + 1 C k div -G(c(ϕ)) ∇ϕ |∇ϕ| + 1 |∇ϕ| P ∇ϕ ⊥ (∇[|∇ϕ|G (c(ϕ))]) ∇ϕ |∇ϕ| |∇ϕ| 1 ε ζ ϕ ε
where

R e = LU ρ ref η ref , W e = η ref U Λ , C k = η ref U L 2 κ
are the Reynolds, Weissenberg and capillary numbers. In a shear flow one important quantity is the shear rate γ, from which we can express the characteristic velocity U = γL. Substituting we finally get

W e = η ref γL Λ , C k = η ref γL 3 κ
which are the dimensionless parameters from [1,3]. C k is a measure of bending distortion of the vesicle. The higher C k is the easiest the bending mode is (note that bending rigidity appears in the denominator. W e measures the stretching (or dilatation) modes. Thus Λ has to be taken large enough in order to prevent significant membrane extensibility.

There are two additional dimensionless parameters, namely the viscosity contrast λ and the reduced volume (or reduced area in two dimensions) ν, defined by

λ = η in η out , ν 3D = 6 √ πV A 3/2 , ν 2D = A π[P/2π] 2 (18)
where V (respectively A) is the vesicle volume (respectively area) and P is the perimeter in the two dimensional problem. In three dimensions (respectively two dimensions), ν corresponds to the ratio between the actual volume (respectively area) over the volume (respectively area) of the a sphere (respectively circle) having the same area (respectively perimeter). For a sphere (or circle) ν = 1, and it is less otherwise. As an example a human RBC has a reduced volume ν 0.65. The numerical results presented below correspond to two dimensions, so ν will refer to ν 2D introduced in (18).

D. Numerical results

Below we present simulation results corresponding to generic situations where the vesicle makes the classical tanktreading or tumbling motion, and less classical motion called vacillating-breathing (an intermediate regime which has been described theoretically by one the authors [17], and experimentally by Podgorsky and Steinberg). The vesicle is put in a shear flow by prescribing opposite horizontal velocities on the top and bottom boundary, with free (Neumann) condition on the pressure, periodic boundary condition for the horizontal velocity and pressure on the vertical boundaries, and free condition on vertical velocity. We set Re = 0.0001 (quite close to the Stokes limit), W e = 0.000025, and C k = 0.25. The numerical parameter ε is chosen typically as 1.5 times the grid step, so that the spreaded interface roughly contains 3 points. Note that C k larger than W e means that bending modes (which relax at the time scale given typically by η ref L 3 /κ) are slower than the elastic modes (which time scale like η ref LΛ). In other words the vesicle elastic response is quasi-instantaneous (on the scale on the physical bending response) and tries to keep the local area (perimeter in two dimensions) as close as possible to the initial one. Typical variations of the observed variations is of the order of percent. The volume is conserved with a much better accuracy (typically 1% for a grid size of N = 64 and 0.1% for N = 128 at the end of the presented simulations). We first fix the viscosity contrast λ = 1. We observe that the vesicle reaches a stationary angle (1), as is expected. The membrane (which is fluid) undergoes a tank-treading motion.

Figure 2 corresponds to the same parameters as before, but with a viscosity ratio λ = 8 (the inner viscosity 8 times larger than in the former test), we observe a tumbling motion. It is interesting to see what happens if the external viscosity is lowered by a factor 8. The viscosity contrast is still the same, but in this case the external viscosity is lowered (rather than increasing the internal one, as done above). Tumbling motion still prevails, but here the vesicle is more deformed than in the previous case (the peanut shape of the vesicle is more evident (Fig. 3)). This result is understood by noting that decreasing η out is equivalent to the previous situation provided that one increase C k (as well as W e , while R e is reduced by a factor of 8; actually we have multiplied R e by 8 in this simulation, since we do not want to vary the relative effect of inertia, as we have in mind the pure Stokes limit) by the same amount (a factor of 4), due to the similarity properties of the two situations (actually the two situations are equivalent). Since C k is a measure of bending modes, its increase allows for more flexibility of the vesicle.

Figure 4 shows the variation of the angle of the vesicle main axis with respect to the flow direction (when the angle is zero this means that the vesicle is elongated and aligned along the flow direction). At low enough viscosity contrast the vesicle shows a tang-treading (Fig. 4a) At large enough λ tumbling prevails (Fig. 4c). We have found that at intermediate regime there is a vacillating-breathing regime as shown on Fig. 4b. Actually it can be shown from general arguments in two dimensions that for a fully inextensible vesicle, that a vacillating-breathing mode can not take place in the small deformation regime (i.e. close to a circular shape where only the second order mode is taken into account; the mode of a deformation about a circular shape is proportional to e imψ where ψ is the polar angle and m is an integer. The second mode corresponds to m = 2). On the one hand, at larger deformation the argument does not hold. In addition, recent numerical studies based on the boundary integral formulation [51] did not observe the vacillating-breathing, but only as a transient. However, in that work the inextensibility condition was preserved up to a relative variation of about 10 -3 for a resolution of 128 2 . If the variation is larger, then the vesicle may behave as a capsule, and in that case it has been recently shown that the vacillating-breathing mode precedes indeed the tumbling one on increasing λ. A systematic future study is needed in order to ascertain this more quantitatively.

Finally we draw the phase diagram representing the three types of motion in the plane of the reduced volume and the viscosity contrast, as has been done in [10]. The results are reported on Fig. 5. The line of transition towards tumbling is quantitatively close to that reported in [10], albeit the present line is a bit lower. This can be attributed to confinement (in [10] the most quantitative results are obtained from the boundary integral formulation in an unbounded domain). A phase diagram in two dimensions was also reported recently in [52]. In that work numerical simulations are presented (their Fig. 2) for a fixed reduced area and a fixed viscosity contrast, and the precedes indeed the tumbling one on increasing λ. A systematic future study is needed in order to ascertain this more quantitatively.

Finally we draw the phase diagram representing the three types of motion in the plane of the reduced volume and the viscosity contrast, as has been done in [9]. The results are shear rate is varied instead. Our diagram (Fig. 5) is represented by varying both the reduced area and the viscosity contrast at a fixed shear rate. Considering a given shear rate, the transition from tank-treading, Vacillating-breathing (called swinging in [52]) and tumbling agree qualitatively with ours. A quantitative comparison can not be made at present. Indeed their Reynolds number differs significantly from ours, and in addition thermal noise (included in Ref.

(Fig. 5), as stated by the authors (in their section C): Obviously, thermal noise has a large impact on the vesicle dynamics. Nevertheless, we compared our transition curve with the phase-field simulation developed in [2]. Their transition curve from tank-treading to tumbling, is in good agreement with our, taking into account that we are not, in our case, in a strictly inextensible regime. This also explains the appearance of the vacillating breathing zone in reported on Fig. 5. The line of transition towards tumbling is quantitatively close to that reported in [9], albeit the present line is a bit lower. This can be attributed to confinement (in [9] the most quantitative results are obtained from the boundary integral formulation in an unbounded domain). A phase diagram in two dimensions was also reported recently in [50]. In that work numerical simulations are presented (their Fig. 2) for a fixed reduced area and a fixed viscosity contrast, and the shear rate is varied instead. Our diagram (Fig. 5) is represented by varying both the reduced area and the viscosity contrast at a fixed between tumbling and tank-trading zones. 5: Transition between different modes depending on the reduced volume for a confinement of 0.3 (confinement is defined as the ratio of the vesicle effective diameter (the diameter of a circle having the same perimeter as the actual vesicle) over the channel width. Comparison with numerical results by the phase field method of [2], for ε = 0.035R.

V. SOME BRIDGES BETWEEN THE METHODS

A. Divergence + gradient form of the membrane forces

Let us first show that the total membrane force (the bending and elastic contribution) can be written as divergence of tensor plus a gradient. After some algebra one gets (see Annex II)

F m = ∇ E (|∇ϕ|)|∇ϕ| 1 ε ζ( ϕ ε ) -div E (|∇ϕ|) ∇ϕ ⊗ ∇ϕ |∇ϕ| 1 ε ζ( ϕ ε ) (19)
This form of F m is interesting since the gradient term is useless when plugged as a source term to incompressible (Navier-)Stokes equation (it is absorbed by the pressure gradient). Another way consists in including the gradient term as the divergence of a diagonal tensor:

F m = div E (|∇ϕ|)|∇ϕ|(I - ∇ϕ ⊗ ∇ϕ |∇ϕ| 2 ) 1 ε ζ( ϕ ε ) (20)
The remaining divergence term may be grouped with the divergence term of the fluid equation, turning the fluidstructure problem into the study of a fluid containing a viscoelastic band having the following constitutive law:

σ = -pI + η(∇u + ∇u t ) + E (|∇ϕ|) ∇ϕ ⊗ ∇ϕ |∇ϕ| 1 ε ζ( ϕ ε )
This form is interesting since in the particular case E (r) = λr (which corresponds to an elastic interface with zero area surface at rest) we recover the stress tensor of a Korteweg fluid (setting Φ = Z(ϕ/ε) with Z = ζ). This remark has been used in [53] to prove existence results for fluid-structure coupling problems formulated in Eulerian coordinates.

Of course this remark holds also for the phase-field formulation [3]. In this case the membrane force is given by

F m = ∇(ξ|∇φ|) -div ξ ∇φ ⊗ ∇φ |∇φ| = div ξ|∇φ|(I - ∇φ ⊗ ∇φ |∇φ| 2 ) (21)
Note that I -∇φ⊗∇φ |∇φ| 2 is the projector on the tangent subspace to the vesicle membrane.

B. Link between ξ and |∇ϕ|

The first apparent difference between the level-set and the phase-field methods is that in one case (the level-set) we have a pure advection, while in the second one we enforce the values of φ to be ±1 inside and outside the vesicle by introducing a double well potential. Remember, however (see section IV B), that ϕ if it is initially a distance function, this properties is not preserved under dynamics. One has thus usually to resort to a renormalization such as to keep it a distance function [START_REF] Osher | Level set methods and Dynamic Implicit Surfaces[END_REF] (note that in section IV B we propose an alternative as compared to the usual method). Thus this procedure is somewhat similar to forcing φ in the phase-field to remain constant in the bulk phases. More precisely, first observe that introducing the function Z such that Z(r) = 2 r -∞ ζ(s)ds -1, the phase function Z( ϕ ε ) has a behaviour similar to φ: it is equal to -1 when ϕ < -ε, inside the vesicle, and to 1 when ϕ > ε, outside the vesicle.

Disregarding the gradient terms in the two expressions ( 19) and ( 21) of F m , which play no role as stated above, we remark that thanks to the chain rule we have ∇ϕ ⊗ ∇ϕ |∇ϕ|

1 ε ζ( ϕ ε ) = ∇Z( ϕ ε ) ⊗ ∇Z( ϕ ε ) |∇Z( ϕ ε )|
Thus in order to compare the formulations, regarding the membrane forces, we must compare ξ (see equation 3) from the phase-field method with E (|∇ϕ|) (see equation 12) of the level-set method. It is easy to prove that if ϕ verifies a transport equation with divergence free velocity field, then |∇ϕ| verifies

|∇ϕ| t + u • ∇|∇ϕ| = |∇ϕ|t • (t • ∇u)
and dividing by |∇ϕ| and multiplying by T to obtain

(T log |∇ϕ|) t + u • ∇(T log |∇ϕ|) = T t • (t • ∇u)
i.e. T log |∇ϕ| verifies the same equation as the postulated equation for ξ of [1,3] (equation 9). Of course we can not conclude that ξ = T log |∇ϕ|, since the velocity fields in the two methods should not coincide. Indeed the φ equation in the phase-field approach is not a simple transport equation as in level-set, and this fact affects the elastic forces and thus the dynamics. But up to these technical misfits, one may consider that ξ is more or less proportional to the log of the local extension of the membrane instead of being simply proportional to this stretching. Another consideration is that, roughly speaking, and up to the Ginzburg-Landau intrinsic energy, the phase-field method can be considered as the level-set method where we choose as constitutive law

E(r) = T r(log r -1)
in the elastic membrane energy.

C. Divergence + gradient form of the curvature force

The level-set curvature force ( 14) (a similar form could be extracted in the phase-field formulation) can be put into divergence form. After a lengthy calculation (see Annex III) we can write

F LS c under divergence form F LS c = div σ ε c
with the tensor σ ε c given by

σ ε c = |∇ϕ| G(c(ϕ))I + ∇ϕ |∇ϕ| ⊗ ∇G (c(ϕ)) I - ∇ϕ ⊗ ∇ϕ |∇ϕ| 2 -G (c(ϕ)) I - ∇ϕ ⊗ ∇ϕ |∇ϕ| 2 D 2 φ I - ∇ϕ ⊗ ∇ϕ |∇ϕ| 2 1 ε ζ( ϕ ε ) (22)
In this expression the higher order term under the divergence contains an extension of the second fundamental form of the surface, to the whole space.

D. Comparison of curvature forces with that obtained in [3] The above derivation was general without specifying the form of the function G. Let us take G(r) = κ 2 r 2 (this is the Helfrich form) in ( 14) which now reads

F LS c = κ div - c(ϕ) 2 2 ∇ϕ |∇ϕ| + 1 |∇ϕ| P ∇ϕ ⊥ (∇[|∇φ|c(ϕ)]) 1 ε ζ( ϕ ε )∇ϕ
while the curvature force of (3) is

F AF c = - κ 2 c 3 2 + t • ∇(t • ∇c) ∇φ
which is equivalent to that derived in [3] (equation A13, for c 0 = 0, where c 0 is the spontaneous curvature considered in [3]; recall that the definition of the curvature there has the opposite sign here).

In two dimension, P ∇ϕ ⊥ (v) = (t • v)t thus the formula for F LS c is identical to formula F AF c of [1,3], up to a factor 2 and the convention for the sign of curvature. Note that in [3] it was not realized that the force can be written as a divergence, as shown here. Instead the force was in the numerical study substituted by an approximate expression, leading to F AF c . To arrive to that expression in [1] and [3] the property div t = 0 was used, which is not valid in general. Indeed, t = ∇×φ |∇φ| thus div t = -∇×φ • ∇|∇φ| |∇φ| 2 = -t • ∇|∇φ| |∇φ| = 0. In [1,3] the curvature force (equation ( 8) of [3]), could be obtained from (6) only to leading order of the interface width. The present study shows in fact that expression (8) is in fact exact. While this remark also holds in three dimension, the algebra to prove the coincidence of level-set and phase field expressions is somehow cumbersome and is therefore omitted here.

E. Comparison of curvature forces with those of [44] In [44] the authors found the following elastic stress for the curvature stress (we use their notation):

τ C = -(Φ -div T ) ⊗ ∇ϕ -T D 2 ϕ
where (considering the case C 0 = 0, i.e. no spontaneous curvature) they set

T = αc(ϕ) I - ∇ϕ |∇ϕ| ⊗ ∇ϕ |∇ϕ| ,
and

Φ = - 5 2 c(ϕ) 2 ∇ϕ |∇ϕ| - 2αC(ϕ) |∇ϕ| 2 D 2 ϕ∇ϕ -∆ϕ∇ϕ
We show in Annex IV that up to the cut-off ζ term which is included in the phase-field function, this is the same as the transpose of σ ε c , in the case G(r) = α r 2 2 , which proves that the stresses are the same (recall that all geometric quantity express the same in terms of the phase-field or level-set equation).

VI. CONCLUSION AND FORTHCOMING WORKS

In this paper we investigated links between the phase-field and level-set modeling of immersed elastic membranes subject to curvature energy. We proved that these two formulations are equivalent from a theoretical point of view and validated the level-set formulation by recovering known behavior of phospholipid vesicles under shear flow. Undergoing works concern the generalization of our model to full membrane energy for the elastic membrane. Indeed while phospholipid vesicles only react to local change of area and curvature, red blood cells also react to shear in the tangent plane to their membrane. This is due to the spectrin network underneath. Therefore there is a need for the modeling of the full membrane energy of an immersed interface. This question is currently under investigation.

VII. ANNEXES A. Annex I: differential of curvature energy

Differentiating this energy with respect to ϕ gives

dE c (ϕ)(δ) = Ω G (c(ϕ)) div ∇δ |∇ϕ| - ∇ϕ • ∇δ |∇ϕ| 3 |∇ϕ| 1 ε ζ( ϕ ε )dx + Ω G(c(ϕ)) ∇ϕ • ∇δ |∇ϕ| 1 ε ζ( ϕ ε ) + G(c(ϕ))|∇ϕ| 1 ε 2 ζ ( ϕ ε )δdx (23)
Integration by part of the second term yields:

-

Ω G(c(ϕ))c(ϕ) 1 ε ζ( ϕ ε )δ + G(c(ϕ)) ∇ϕ |∇ϕ| 1 ε 2 ζ ( ϕ ε )∇ϕδ + ∇G(c(ϕ)) • ∇ϕ |∇ϕ| 1 ε ζ( ϕ ε )δ.
The second term cancels with the last one of equation (23). Denoting by P ∇ϕ ⊥ the linear projection operator on ∇ϕ ⊥ ,

P ∇ϕ ⊥ (v) = v -(v • ∇ϕ |∇ϕ| ) ∇ϕ |∇ϕ| .
we thus have

dE c (ϕ)(δ) = Ω G (c(ϕ)) div P ∇ϕ ⊥ (∇δ) |∇ϕ| |∇ϕ| 1 ε ζ( ϕ ε ) -G(c(ϕ))c(ϕ) 1 ε ζ( ϕ ε )δ -∇G(c(ϕ)) • ∇ϕ |∇ϕ| 1 ε ζ( ϕ ε )δdx
which thanks to the expression for c(ϕ) also reads

dE c (ϕ)(δ) = Ω G (c(ϕ)) div P ∇ϕ ⊥ (∇δ) |∇ϕ| |∇ϕ| 1 ε ζ( ϕ ε ) -div(G(c(ϕ)) ∇ϕ |∇ϕ| ) 1 ε ζ( ϕ ε )δdx Since P ∇ϕ ⊥ (∇δ) • ∇ϕ = 0 the first term integrates by part into - Ω ∇ [|∇ϕ|G (c(ϕ))] • P ∇ϕ ⊥ (∇δ) 1 |∇ϕ| 1 ε ζ( ϕ ε ) = - Ω P ∇ϕ ⊥ [|∇ϕ|∇G (c(ϕ))] • ∇δ |∇ϕ| 1 ε ζ( ϕ ε )dx
where we used the symmetry property of the projection on ∇ϕ ⊥ . Integrating by parts once again, there holds

dE c (ϕ)(δ) = Ω div -G(c(ϕ)) ∇ϕ |∇ϕ| + 1 |∇ϕ| P ∇ϕ ⊥ (∇[|∇ϕ|G (c(ϕ))]) 1 ε ζ( ϕ ε )δdx B.
Annex II: divergence form of membrane energy

Let us start with the elastic force (12) and rewrite it as

F m = ∇[E (|∇ϕ|)]|∇ϕ| 1 ε ζ( ϕ ε ) -div E (|∇ϕ|) ∇ϕ |∇ϕ| ∇ϕ 1 ε ζ( ϕ ε ) =: A -B
we use the following notations and elementary results from tensorial analysis:

(i) For two vectors a and b, a ⊗ b is the matrix with coefficient a i b j on line i, column j.

(ii) This matrix is such that (a ⊗ b)c = (b • c)a; thus n ⊗ n is the projector on n.

(iii) For two vector fields a and b, div(a ⊗ b) = (∇a)b + (div b)a.

(iv) For a scalar function f and a vector field a, ∇(f a) = a ⊗ ∇f + f ∇a.

The B term in F m may be transformed by taking a

= 1 ε ζ( ϕ ε )∇ϕ and b = E (|∇ϕ|) ∇ϕ |∇ϕ| in (iii) to get B = div E (|∇ϕ|) ∇ϕ ⊗ ∇ϕ |∇ϕ| 1 ε ζ( ϕ ε ) -∇( 1 ε ζ( ϕ ε )∇ϕ)E (|∇ϕ|) ∇ϕ |∇ϕ| Observe that ∇( 1 ε ζ( ϕ ε )∇ϕ) = ∇ϕ ⊗ ∇( 1 ε ζ( ϕ ε )) + 1 ε ζ( ϕ ε )∇ 2 ϕ = ∇( 1 ε ζ( ϕ ε )) ⊗ ∇ϕ + 1 ε ζ( ϕ ε )∇ 2 ϕ
(in general a ⊗ b = b ⊗ a but equality occurs when a and b are collinear). Moreover from (ii) one has

(∇( 1 ε ζ( ϕ ε )) ⊗ ∇ϕ) ∇ϕ |∇ϕ| = |∇ϕ|∇( 1 ε ζ( ϕ ε ))
and

∇ 2 ϕ ∇ϕ |∇ϕ| = ∇|∇ϕ|.
Thus B may be written as

B = div E (|∇ϕ|) ∇ϕ ⊗ ∇ϕ |∇ϕ| 1 ε ζ( ϕ ε ) -E (|∇ϕ|)∇(|∇ϕ| 1 ε ζ( ϕ ε ))
and combined with A so that F m reduces to (19).

C. Annex III: divergence form of curvature energy

To get that divergence form we start from ( 14) that we recall for convenience

F c = div -G(c(ϕ)) ∇ϕ |∇ϕ| + 1 |∇ϕ| P ∇ϕ ⊥ (∇[|∇ϕ|G (c(ϕ))]) 1 ε ζ( ϕ ε )∇ϕ
and first use the tensor identity (iii) in V A which gives

F c = div -G(c(ϕ)) ∇ϕ ⊗ ∇ϕ |∇ϕ| + 1 |∇ϕ| ∇ϕ ⊗ P ∇ϕ ⊥ (∇[|∇ϕ|G (c(ϕ))]) 1 ε ζ( ϕ ε ) -∇ 1 ε ζ( ϕ ε )∇ϕ × -G(c(ϕ)) ∇ϕ |∇ϕ| + 1 |∇ϕ| P ∇ϕ ⊥ (∇|∇ϕ|G (c(ϕ)))
that we denote by div A -B. Note that × denotes a matrix-vector product in the above formula. Working on the second term and using (iv

) of V A, ∇ 1 ε ζ( ϕ ε )∇ϕ = ∇ϕ ⊗ ∇ 1 ε ζ( ϕ ε ) + 1 ε ζ( ϕ ε )D 2 ϕ = 1 ε 2 ζ ( ϕ ε )∇ϕ ⊗ ∇ϕ + 1 ε ζ( ϕ ε )D 2 ϕ = ∇ 1 ε ζ( ϕ ε ) ⊗ ∇ϕ + 1 ε ζ( ϕ ε )D 2 ϕ. Moreover, recall that D 2 ϕ∇ϕ |∇ϕ| = ∇|∇ϕ| and ∇ 1 ε ζ( ϕ ε ) ⊗ ∇ϕ P ∇ϕ ⊥ (u) = 0
for any u using (ii) from V A. Thus

that we can arrange like in (22).

D. Annex IV: comparison of curvatures forces

From section V E, we first write Φ in another way. Indeed we note that after an expansion of the curvature expression (which is c(ϕ) = div(∇ϕ/|∇ϕ|)) we obtain

c(ϕ) = 1 |∇ϕ| ∆ϕ - D 2 ϕ∇ϕ • ∇ϕ |∇ϕ| 2 .
One also has which is a spherical tensor, thus not modifying the dynamics in an incompressible flow. Therefore the two curvature stress tensors are identical.

Φ • ∇ϕ |∇ϕ| = - 5α 2 c(ϕ) 2 + 2αc(ϕ) 2 = - α 2 c(ϕ)

|∇ϕ| 2 :

 2 D(u) therefore, starting from |∇ϕ 0 | = 1 initially, we have |∇ϕ| = 1 for t > 0 if and only if ∇ϕ ⊗ ∇ϕ |∇ϕ| 2 : D(u) = 0 or equivalently ∇ϕ ⊗ ∇ϕ |∇ϕ| 2 : ∇u = 0

  ρ ref and η ref represent characteristic length, velocity, density and viscosity scales. Accordingly we set x = Lx , u = U u , t = (L/U )t , ρ = ρ ref ρ , η = η ref r, p = η ref (U/L)p , φ = Lφ , and ε = Lε . Differentiating we find

FIGURE 1 :

 1 FIGURE 1: Stationary stream lines for the tank-treading motion of a vesicle in a shear flow. Viscosity contrast: λ = 1, reduced area ν = 0.7.

FIGURE 2 :

 2 FIGURE 2: Stream lines for the tumbling motion of vesicle in a shear flow. Viscosity contrast: 8 at times 4, 12, 20, 28 (left to right, and top to bottom). Reynolds: 10 -4 .

FIGURE 2 :

 2 FIGURE 2: Stream lines for the tumbling motion of vesicle in a shear flow. Viscosity contrast: 8 at times 4, 12, 20, 28 (left to right, and top to bottom). Reynolds: 10 -4 .

FIGURE 3 :

 3 FIGURE 3: Stream lines for the Tumbling motion of vesicle in a shear flow, at time 4, 12, 20, and 28 (left to right, and top to bottom). Viscosity contrast: 8. R e = 8.10 -4 .

FIGURE 3 :

 3 FIGURE 3: Stream lines for the Tumbling motion of vesicle in a shear flow, at time 4, 12, 20, and 28 (left to right, and top to bottom). Viscosity contrast: 8. Re = 8.10 -4 .

  (a) Tank treading for λ = 2 (b) Vascillating breathing for λ = 5.6 (c) Tumbling for λ = 7.

FIGURE 4 :

 4 FIGURE 4: Angle variation depending on the viscosity ratio, for ν = 0.8.

FIGURE

  FIGURE 5: Transition between different modes depending on the reduced volume for a confinement of 0.3 (confinement is defined as the ratio of the vesicle effective diameter (the diameter of a circle having the same perimeter as the actual vesicle) over the channel width. Comparison with numerical results by the phase field method of[2], for ε = 0.035R.

2 2 which 2 + 2 Finally 2 -

 22222 From the definition of projector,P ∇ϕ ⊥ (u) = I -∇ϕ ⊗ ∇ϕ |∇ϕ| 2 uand by using u⊗ Av = (u ⊗ v)A T ,Let us try to write B as a divergence term minus a remainder term. To start with,B = -∇ G(c(ϕ)computing ∇c(ϕ) leads to ∇c(ϕ) = ∇ div ∇ϕ |∇ϕ| = div ∇ ∇ϕ |∇ϕ| T Therefore B = -∇ (...) + G (c(ϕ))|∇ϕ| div ∇ ∇ϕ |∇ϕ| term is zero thanks to the expression of ∇ ∇ϕ |∇ϕ| Twith the projector on ∇ϕ ⊥ and the fact that ∇ 1 ε ζ( ϕ ε ) is colinear to ∇ϕ. Writing the gradient term as the divergence of a diagonal tensor we get a first expression in divergence form:F c = div -G(c(ϕ))|∇ϕ| I -∇ϕ ⊗ ∇ϕ |∇ϕ| 2 + ∇ϕ |∇ϕ| ⊗ P ∇ϕ ⊥ ∇ (|∇ϕ|G (c(ϕ)))We may arrange terms further by using u ⊗ Av = (u ⊗ v)A T , with v = ∇|∇ϕ|G (c(ϕ)) gives ∇ϕ |∇ϕ| ⊗ P ∇ϕ ⊥ (∇|∇ϕ|G (c(ϕ))) = ∇ϕ |∇ϕ| ⊗ ∇ (|∇ϕ|G (c(ϕ))) ∇ϕ ⊗ ∇G (c(ϕ)) I -∇ϕ ⊗ ∇ϕ |∇ϕ| 2 and -G (c(ϕ))|∇ϕ|∇ ∇ϕ |∇ϕ| T = -G (c(ϕ))D ϕ I -∇ϕ ⊗ ∇ϕ |∇ϕ| we get the expression F c = div σ ε c with σ ε c = G(c(ϕ))|∇ϕ| I -∇ϕ ⊗ ∇ϕ |∇ϕ| 2 + (∇ϕ ⊗ ∇G (c(ϕ))) I -∇ϕ ⊗ ∇ϕ |∇ϕ|

2

 2

2 ϕ∇ϕ-D 2

 22 The transpose of the level-set stress is given by that each term in the expression of τ C is hidden in this expression. Consider the term T D 2 ϕ which is included in the last term of (25), if we consider only the identity I in the last projection term. The remainder is αc(ϕ) on the term div T ⊗ ∇ϕ. We first compute div T = div αc(ϕ) ϕ)2 ∇ϕ |∇ϕ|Thus taking the tensorial product givesdiv T ⊗ ∇ϕ = α I -The first term is exactly the second term of(25). The second term has the wrong sign to match(26). But combined with the new expression for -Φ ⊗ ∇ϕ, namely ϕ∇ϕ ⊗ ∇ϕ this gives the remaining terms, up to the extra term in the LS stressα c(ϕ) 2 2 |∇ϕ|I

Acknowledgements

This work was supported by the French Ministry of Education through ANR MOSICOB and by University Joseph Fourier by PPF DYSCO. E.M. was also supported by ANR COMMA. C.M. Acknowledges financial support from CNES.


Phil Hoole

Lasse Bombien

Barbara Kühnert
email: bkuhnert@univ-paris3.fr

Christine Mooshammer

Intrinsic and Prosodic Effects on Articulatory Coordination in Initial Consonant Clusters

EMA was used to study the coordination of the articulatory gestures for C1 and C2 in onset clusters, firstly as a function of the segmental make-up of the clusters, and secondly as a function of stress and prosodic boundary conditions. The segmental results, which compared German and French, indicated a much lower degree of overlap of C1 and C2 for C2=/n/ compared to C2=/l/ (with C1=/p, b, k, g/). Overlap was also less for voiceless compared to voiced C1 for German. Prosodic boundary strength affected the duration of the C1 articulatory constriction more strongly than that of C2, while, conversely, differences in lexical stress affected C2 more strongly than C1. The coordination relations between C1 and C2 were, however, not systematically affected by prosodic conditions. Nor did the different clusters studied (/kl, kn, sk/) appear to differ in their internal cohesiveness. The results are discussed with respect to how articulatory coordination is constrained to allow acoustic recovery of segmental information by the listener; the possibility of cross-language differences in these constraints is raised.

Introduction

This study examines how articulatory coordination in heterorganic initial consonant clusters is modulated firstly by the segmental make-up of the cluster itself (see Section 2), and secondly by the prosodic conditions in which the cluster is spoken (see Section 3).

Segmental structure of clusters 2.1 Background

The main impetus for this first area comes from the contention advanced in [START_REF] Chitoran | Gestural overlap and recoverability: Articulatory evidence from Georgian[END_REF] that generalizations about preferred syllable structures (such as the sonority hierarchy) reflect at a more fundamental level the outcome of successful compromises between parallel transmission of segmental information, i.e a high degree of overlap among gestures (assumed to be efficient for the speaker) and clear modulation of the acoustic signal (assumed to be efficient for the hearer). One illustration of this is the place-order effect in stop-stop clusters reported for Georgian by these authors: C1(back) + C2(front) clusters show less overlap than the reverse order, presumably to ensure recoverability of C1 by the listener.

Nonetheless, relative to the very large number of complex syllabic onsets found across languages there has still been little detailed exploration of these principles. We focus here on a potentially revealing comparison among onset clusters found, for example, in German and French, but not in English, namely plosive + /l/ vs. plosive + /n/. The latter appear to be less stable diachronically than the former (cf. [START_REF] Vennemann | Triple-cluster reduction in Germanic: Etymology without sound laws? Historische Sprachforschung[END_REF], and occur less frequently in extant languages. Why might they then be less suitable for the formation of complex syllable onsets? One possibility is that C1 and C2 need to overlap less (than in the case of plosive + /l/) to ensure that acoustic properties of the C1 burst are not compromised by velar lowering.

As a further manipulation of the phonetic properties of the syllable onset, the study also systematically varies the voicing category of the C1 plosive. This can be assumed to have a substantial effect on the nature of the acoustic information available to the hearer particularly in the C1 to C2 transition. But whether this may in turn lead to different preferred coordination patterns for C1 with C2 is, to our knowledge, essentially unknown. German and French were considered to be particularly appropriate for study here, since they are known to differ substantially in the timing of voice onset following voiceless plosives.

Material and Recordings

The main structure aimed for in the material for both German and French was to use real monosyllabic words (CCVC) with all combinations of C1 = /p, b, k, g/, C2 = /l, n/, V=/i, a/. Not all of these combinations are lexically attested; for example, /bn/ combinations are missing in both languages and /pn/ is lexically somewhat marginal in German.

Since preliminary analysis indicated no consistent effects of vowel category on overlap between C1 and C2, and since gaps in the lexicon occasionally required the use of vowels other than /i/ and /a/, for the purposes of the present paper we will present the results averaged over vowels. In short, results can be presented for all four combinations of C1 with /l/, and three with /n/ (/bn/ missing).

Each target word was embedded in a carrier phrase with additional material not of concern here and recorded 10 times in randomized order.

Movements of the articulators were measured by means of a 3D EMA system (AG500, Carstens Medizinelektronik). Of the three sensors located on the tongue the frontmost one was used to analyze the tongue-tip gesture for /l, n/, the rearmost one for the tongue-dorsum gesture for /k, g/. The sensors located on lower and upper lip were combined to a lip aperture measure for analysis of the labial gesture for /p, b/. The key dependent variable is the measure of overlap between C1 and C2. This is illustrated in Fig. 1 ). The one used here is based on the constriction phases of C1 and C2. Specifically, the timepoint of the onset of the constriction phase of C2 (phase 4 in the figure) is subtracted from the timepoint of the offset of the constriction phase of C1 (phase 2 in the figure). This time difference is normalized by the time interval from the onset of C1 constriction to offset of C2 constriction and expressed as a percentage. The important point to bear in mind for the interpretation of the results is that, with this definition, increasingly positive values correspond to increasing overlap of C1 and C2 (negative values in effect mean that the constriction phases do not overlap).

Results

Regarding the manner of articulation of C2 (lateral vs. nasal), for German there is a clear effect in the form of the anticipated lower degree of overlap for the plosive-nasal sequences (see Fig. 2). This appears to be a very robust effect because the same effect is observed (for the cluster pair /kl/, /kn/) in the different corpus used for the prosodic investigation reported in section 3 below. Moreover, with the prosody corpus used in section 3 the effect has been observed not only for the speakers analyzed there with EMA but also for an additional 7 speakers recorded by electropalatography.

For French the pattern is less consistent, the speakers showing strong, weak and non-existent trends in the expected direction for AM, NN and CG, respectively (see Fig. 3). For CG the overlap between plosive and /n/ oftened seemed large enough to lead to nasal release of the plosive.

Regarding voicing of C1 there is a very consistent effect for German showing less overlap between C1 and C2 when C1 is voiceless. For French no consistent effect of voicing is seen. On the basis of only three speakers per language it is, of course, too early to be confident about possible cross-linguistic differences. Nonetheless, it is intriguing that voicing-related differences are observed in that language, namely German, where acoustic properties of C2 are more strongly affected by the voicing status of C1. Perhaps this fits in with the first finding too: The constraint to maintain low overlap between plosive and /n/ may be stronger in languages where glottal opening is wide at the release of C1.

Prosodic influences on clusters

Background and material

The second main area of the study extends to clusters the wellestablished experimental paradigm (e.g Fougeron & Keating, 1997; [START_REF] Keating | Domain-initial articulatory strengthening in four languages[END_REF] of investigating articulatory strengthening at prosodic boundaries, as part of an overall concern to better understand what information on prosodic structure the speaker makes available to the listener. A plausible hypothesis -consistent, for example, with the pi-gesture approach of [START_REF] Byrd | The elastic phrase: Modeling the dynamics of boundary-adjacent lengthening[END_REF] -is that prosodic boundary effects on articulation decrease with increasing distance from the boundary. However, with a few exceptions (e.g Byrd & Choi, in press), most work within this paradigm has used rather simple syllable structures. Thus it is currently difficult to estimate whether boundary effects are indeed best modelled as a simple function of distance from the boundary, or whether structural properties of the syllable (e.g consonantal onset versus vocalic nucleus) have an independent influence.

The material investigated included the clusters /kl, kn, sk/; they formed the onset of the target words, i.e the words immediately following the prosodic boundaries, varied in four steps from utterance initial to no boundary, as shown in the following examples for the /kl/ cluster at the onset of "Claudia": U Thomas studiert in Fulda. Claudia geht noch zur Schule.

(Thomas goes to college in Fulda. Claudia is still in school.) P Olga sagt immer, Claudia sei zu jung.

(Olga always says that Claudia is too young.) L Thomas, Peter, Claudia und Elke fahren in den Süden.

(Thomas, Peter, Claudia and Elke are travelling south.) W Gestern war Claudia noch gesund.

(Yesterday, Claudia was still well.) (U=Utterance initial, P=Phrase initial, L=List element, W=Word boundary (i.e no prosodic boundary expected))

The basic expectation was that C1 would exhibit more boundary-induced variation than C2.

In addition to the boundary condition the experiment also contained orthogonal variation of a second prosodic variable, namely word stress: The target clusters formed the onsets of words with lexical stress on either the first or second syllable (e.g 'Claudia, vs. Klau'sur). For this prosodic manipulation the literature gives even less to go on than for the boundary condition. A tentative hypothesis was based on the assumption that stress effects are centered on the nucleus of the stressed syllable and decline towards the syllable margins. Accordingly, stress effects were expected to be more visible in C2 than C1.

Analyses focussed on whether C1 or C2 showed durational or spatial enhancement at higher boundary levels and in the stress condition. In addition it was hypothesized that overlap between C1 and C2 would be reduced at higher boundary levels and in the stress condition.

In addition to the interest in prosodic effects in their own right, the manipulation of prosodic conditions in the present experiment can also be regarded as a probe to elucidate the internal cohesion of clusters (cf. Byrd & Choi, in press). We have seen above that /kn/ had less overlap of the component oral gestures than /kl/. Could it also be the case that the components are also less cohesive in /kn/, such that the coordination relations (as captured by the overlap measure) are more strongly affected by 'external' conditions? (An interesting distributional fact about /kn/ in French is that unlike /kl/ and /kr/ it cannot be followed by a glide. [START_REF] Rialland | The phonology and phonetics of extrasyllabicity in French[END_REF] has argued on this basis for the extrasyllabicity of /k/ in /kn/.)

Recordings and analysis procedures

EMA recordings have been analyzed for the three German speakers who participated in section 2 above. For this part of the recording session, in addition to the EMA data, respiratory activity of the speakers was monitored by means of Respitrace in order to give additional information on prosodic divisions into breath groups.

As a preliminary stage in the analysis the actual prosodic realization at the target word of the four boundary categories given above was assigned to one of three prosodic categories, which will be referred to as follows (cf. [START_REF] Cho | Prosodic influences on consonant production in Dutch: Effects of prosodic boundaries, phrasal accent and lexical stress[END_REF]: Big Boundary (BG):

Pause and boundary tone Small Boundary (SM):

No pause but boundary tone Prosodic Word (WD):

No pause, no boundary tone

In the present dataset only very few big boundaries (BG) were found, so for the following results only the SM and WD categories will be considered. This is also the reason why the results of a previous experiment in which 7 speakers were recorded with basically the same corpus by means of EPG will not be examined here: In the EPG experiment a large number of BG realizations were found (essentially for the utterance-initial condition), which would unduly complicate the comparison with the present EMA results. The results of the EPG experiment will be considered in detail in a separate publication (Mooshammer et al., in preparation).

Results

We will first present the results for the durations of the target phases of C1 and C2, i.e the durations of the phases labelled '2' and '4' in the example utterance given in Fig. 1 above.

The results for C1 are shown in Fig. 4. There is a very consistent trend in the hypothesized direction of longer durations at the more prominent prosodic boundary (i.e the 'SM' boundary compared to the 'WD' boundary): Of 18 SM-WD pairs matched for stress level (and cluster and speaker), 17 are in the predicted direction. In ANOVAs carried out separately for each combination of speaker and cluster type the main effect of boundary is significant at p<0.01 in 6 out of 9 possible cases. On the other hand, the effect of stress on C1 appears to be very inconsistent. Of the eighteen stress-unstress pairs matched for boundary type immediately adjacent bars in the figure), 10 are in the direction of longer duration for the stress condition, and in the ANOVAs only one of nine cases showed a main effect of stress in this direction that was significant at p<0.01 (/sk/ for speaker MO).

The results for C2 are effectively a mirror-image of those for C1: The effect of stress is highly consistent, with 16 of 18 matched pairs of bars in Fig. 5 showing the expected trend of longer durations for the stress condition (main effect was significant at p<0.01 in 3 of 9 cases). On the other hand there are very few clear effects of boundary condition: Even if the balance of the overall trend is in the direction of longer durations at higher boundary levels (14 of 18 pairs), only one of the 9 ANOVAs gave a main effect that reached p<0.01 (/sk/ for speaker IK). Thus taken together, the results for C1 and C2 durations correspond quite closely to the original hypotheses that boundary effects would be more salient in C1 and stress effects more salient in C2. Regarding the maximum temporal extent of the effects, it seems at the moment possible that boundary effects may extend into C2 more often than stress effects extend to C1, but judgement on this must be reserved: Several more speakers have been recorded, but not yet analyzed

The other main parameter relating to the temporal organization of the clusters is the overlap between C1 and C2 (defined exactly as in section 2.2 above). This is shown in Fig. 6. The main point to note is that there is a strong contrast with the durational results for C1 and C2 just shown above: No trends with respect to either boundary or stress condition can be observed. In fact, counting matched pairs for the boundary comparison and the stress comparison gives 9/18 in both cases.

There is effectively only one clear case where the hypothesized greater degree of overlap in the unstressed condition is found (/kl/ for speaker IK), and even this is counterbalanced by one clear case in the opposite direction (/sk/ for speaker MO).

In short, while prosodic conditions clearly -and not surprisingly -affect the durational properties of gestures they appear to have relatively little influence on the coordination relations between gestures in onset clusters. This result is consistent with that of Byrd & Choi (in press), at least in the sense that in their data it appeared much less easy to demonstrate a consistent effect of boundary condition on onset clusters compared to clusters in other syllabic positions.

This absence of prosodic effects on overlap contrasts markedly with the effect of the segmental make-up of the clusters: Fig. 6 reinforces the results of the first part of the investigation with the radically different coordination patterns between the two superficially very similar clusters /kl/ and /kn/.

This leads to the final issue to consider with respect to overlap, namely whether the varying prosodic conditions can reveal differences in the internal cohesion of different clusters. Byrd & Choi (in press) considered this issue for /sp, sk/ vs. /kl/, hypothesizing that the fricative-onset clusters could represent particularly cohesive gestural 'molecules'. However, no clear differences were found. Similarly, here, there is little evidence from Fig. 6 that the three clusters /kl, kn, sk/ differ in the degree to which the prosodic manipulation affects the degree of overlap (we had speculated, in particular, that /kn/ might show less cohesion than /kl/).

The final aspect of the results to consider is whether prosodic influences manifest themselves not only temporally but also spatially. While this has not yet been exhaustively explored for the EMA data presented here, preliminary analysis indicates negligible effects of the boundary and stress conditions on positions and velocities of the articulators. This fits in with the results of our EPG investigation mentioned above, which also found little consistent evidence for articulatory strengthening in the form of increased tongue-palate contact at prosodically stronger locations. While this may appear to contrast with some earlier findings in the literature it should be noted that dorsal consonants and clusters have scarcely been investigated (unlike single coronal consonants). Perhaps activation in the dorsal region tends to simultaneously reduce the susceptibility to variation of an adjacent coronal articulation. Our corpus actually includes a wider range of clusters than those presented here. This issue may thus become clearer when this material has been analyzed, too.

Conclusion and Outlook

Concluding with an attempt to tie together the two different areas of this study, it was interesting to note that the prosodic effects on clusters with respect to patterns of overlap, for example, were much more subtle than the robust effects found when contrasting the different segmental structures plosive + /l/ vs. plosive + /n/. It is tempting to speculate that the robust effects emerge because intergestural timing forms an integral part of the phonological representation of words, whereas the weaker prosodic effects may reflect the greater lee-way at the disposal of speakers for implementing prosodic structure, and may also reflect greater variability (particularly, perhaps, in the context of a lab-speech experiment) in how speakers interpret the communicative demands of specific utterances. This pattern of findings may be a further facet of the often observed stability of onset clusters: this is probably the syllable-position where lexically-specified timing patterns dominate overall temporal organization most strongly.

Although the difference in timing patterns between /kl/ and /kn/ appears robust we have not provided a direct answer as to why this pattern occurs. The obvious next step would be to use articulatory synthesis to systematically vary gestural organization: Assuming oral and velar gestures for /n/ are organized more or less in synchrony then earlier timing of /n/ relative to /k/ would eventually lead to a nasal release of /k/. Even before this stage the acoustic salience of the release phase of /k/ may be reduced, leading to a situation that speakers try to avoid. This line of thought may also be relevant for the lower degree of overlap for voiceless vs. voiced C1. In languages such as German that need to distinguish a very large number of syllable onsets and in which the timing of peak glottal opening for voiceless plosives is close to the time-point of the oral release then the clarity of the acoustic information on the release of C1 may be improved by not overlapping too strongly with C2.

could simultaneously also enhance the salience of the C2 to vowel transition by increasing the chances that it is voiced rather than voiceless. The latter constraint, at least, would be of less concern in a language such as French where voice onset is timed to occur soon after the release of voiceless plosives. Currently, although more subjects have been recorded, not enough have been analyzed to be confident yet that crosslanguage differences in laryngeal timing may be accompanied by different preferred patterns in overlap of supraglottal gestures. But the results are suggestive enough to make this avenue appear well worth exploring further.

A final area to follow up with the present data is whether the clusters show a stable pattern of timing with the following vowel, such as the c-center effect, despite the differences in within-cluster timing (see Goldstein et al., 2007, for a revealing example of how cross-linguistic differences in syllable structure are reflected in different patterns of timing of consonant clusters with the following vowel). For two of the French speakers presented here [START_REF] Kühnert | Gestural overlap and C-center in selected French consonant clusters[END_REF] found quite stable c-center locations; this is currently being followed up for the German speakers.

C1=K C1=G C1=P C1=B

Figure 2 :Figure 3 :

 23 Figure 2: C1-C2 overlap for three German subjects (see text for definition; expressed here as percentage). From left to right bars are /kl, kn, gl, gn, pl, pn, bl/ (/bn/ not in corpus). Error bars give standard error of mean.

Figure 4 :

 4 Figure 4: C1 constriction durations for three German subjects for the onset clusters /kl, kn, sk/. Each pair of bars gives the results for one prosodic boundary condition; left (darker) bar the stressed (S), right (lighter) bar the unstressed (U)condition. each cluster, small (SM) and prosodic word (WD) prosodic boundary conditions are given in adjacent pairs of bars.Error bars give standard error of mean.

Figure 5 :

 5 Figure 5: C2 constriction durations. Other details as for Fig. 4.

Figure 6 :

 6 Figure 6 : C1-C2 relative overlap. A value of +1.0 would indicate 100% overlap. Other details as for Fig. 4.

  . Different measures of overlap have been

	used in previous investigations (Chitoran et al., 2002, Gafos et
	al., in press
	[k	Nnnnnn	l	a	k]
		n
				Tongue
	3	4	Tip
			Ynnnn
		Nn	Nnnyy	Tongue
				Back
	2
	1
		Time (s)
	Figure 1 Illustration of determination of C1-C2 overlap for a
	/kl/ sequence. The constriction phases for the two consonants
	(labelled '2' and '4') are determined using a 20% velocity
	criterion in the tongue-back and and tongue-tip signals
	respectively. For illustrative purposes only the vertical
	movement component is shown.

Acknowledgement

Work supported by German Research Council (DFG) grants HO 3271/3 and MO 1687/1 to Phil Hoole & Tine Mooshammer for the project "Articulatory cross-language study of initial consonant clusters in varying prosodic conditions" within DFG Schwerpunktprogramm 1234 "Phonological and phonetic competence: between grammar, signal processing, and neural activity".