Roberto Feola

Felice Iandoli
email: felice.iandoli@sorbonne-universite.fr

Federico Murgante
email: fmurgant@sissa.it

Laboratoire Jacques-

LONG-TIME STABILITY OF THE QUANTUM HYDRODYNAMIC SYSTEM ON IRRATIONAL TORI

We consider the quantum hydrodynamic system on a d-dimensional irrational torus with d = 2, 3. We discuss the behaviour, over a "non-trivial" time interval, of the H s -Sobolev norms of solutions. More precisely we prove that, for generic irrational tori, the solutions, evolving form ε-small initial conditions, remain bounded in H s for a time scale of order O(ε -1-1/(d-1)+ ), which is strictly larger with respect to the time-scale provided by local theory. We exploit a Madelung transformation to rewrite the system as a nonlinear Schrödinger equation. We therefore implement a Birkhoff normal form procedure involving small divisors arising form three waves interactions. The main difficulty is to control the loss of derivatives coming from the exchange of energy between high Fourier modes. This is due to the irrationality of the torus which prevent to have "good separation" properties of the eigenvalues of the linearized operator at zero. The main steps of the proof are: (i) to prove precise lower bounds on small divisors; (ii) to construct a modified energy by means of a suitable high/low frequencies analysis, which gives an a priori estimate on the solutions.

INTRODUCTION

We consider the quantum hydrodynamic system on an irrational torus of dimension 2 or 3    ∂ t ρ = -m∆φ -div(ρ∇φ)

∂ t φ = -1 2 |∇φ| 2 -g(m + ρ) + κ m + ρ ∆ρ - κ 2(m + ρ) 2 |∇ρ| 2 , (QHD)
where m > 0, κ > 0, the function g belongs to C ∞ (R + ; R) and g(m) = 0. The function ρ(t, x) is such that ρ(t, x) + m > 0 and it has zero average in x. The space variable x belongs to the irrational torus

T d ν := (R/2πν 1 Z) × • • • × (R/2πν d Z) , d = 2, 3 , (1.1
) with ν = (ν 1 , . . . , ν d ) ∈ [START_REF] Antonelli | Analysis of acoustic oscillations for a class of hydrodynamic systems describing quantum fluids[END_REF][START_REF] Audiard | Global Well-Posedness of the Euler-Korteweg System for Small Irrotational Data[END_REF] d . We assume the strong ellipticity condition g (m) > 0.

(1.2)

We shall consider an initial condition (ρ 0 , φ 0 ) having small size ε 1 in the standard Sobolev space H s (T d ν ) with s 1. Since the equation has a quadratic nonlinear term, the local existence Felice Iandoli has been supported by ERC grant ANADEL 757996.

theory (which may be obtained in the spirit of [START_REF] Berti | Local well posedness of the Euler-Korteweg equations on T d[END_REF][START_REF] Feola | Local well-posedness for the Hamiltonian quasi-linear Schrödinger equation on tori[END_REF]) implies that the solution of (QHD) remains of size ε for times of magnitude O(ε -1 ). The aim of this paper is to prove that, for generic irrational tori, the solution remains of size ε for longer times.

For φ ∈ H s (T d ν ) we define

Π 0 φ := 1 (2π) d ν 1 • • • ν d T d ν φ(x)dx , Π ⊥ 0 := id -Π 0 . (1.3)
Our main result is the following.

Theorem 1.1. Let d = 2 or d = 3. There exists s 0 ≡ s 0 (d) ∈ R such that for almost all ν ∈ [1, 2] d , for any s ≥ s 0 , m > 0, κ > 0 there exist C > 0, ε 0 > 0 such that for any 0 < ε ≤ ε 0 we have the following. For any initial data

(ρ 0 , φ 0 ) ∈ H s 0 (T d ν ) × H s (T d ν ) such that ρ 0 H s (T d ν ) + Π ⊥ 0 φ 0 H s (T d ν ) ≤ ε , (1.4)
there exists a unique solution of (QHD) with (ρ(0), φ(0)) = (ρ 0 , φ 0 ) such that

(ρ(t), φ(t)) ∈ C 0 [0, T ε ); H s (T d ν ) × H s (T d ν ) C 1 [0, T ε ); H s-2 (T d ν ) × H s-2 (T d ν ) , sup t∈[0,Tε) ρ(t, •) H s (T d ν ) + Π ⊥ 0 φ(t, •) H s (T d ν ) ≤ Cε , T ε ≥ ε -1- 1 d-1 log -d+1 2 1 + ε 1 1-d . (1.5)
Derivation from Euler-Korteweg system. The (QHD) is derived from the compressible Euler-Korteweg system

∂ t ρ + div(ρ u) = 0 ∂ t u + u • ∇ u + ∇g(ρ) = ∇ K(ρ)∆ρ + 1 2 K (ρ)|∇ρ| 2 , (EK)
where the function ρ(t, x) > 0 is the density of the fluid and u(t, x) ∈ R d is the time dependent velocity field; we assume that K(ρ), g(ρ) ∈ C ∞ (R + ; R) and that K(ρ) > 0. In particular, in (QHD), we assumed

K(ρ) = κ ρ , κ ∈ R + . (1.6)
We look for solutions u which stay irrotational for all times, i.e.

u = c(t) + ∇φ , c(t) ∈ R d , c(t) = 1 (2π) d ν 1 • • • ν d T d ν udx , (1.7)
where φ : T d ν → R is a scalar potential. By the second equation in (EK) and using that rot u = 0 we deduce

∂ t c(t) = - 1 (2π) d ν 1 • • • ν d T d ν u • ∇ udx = 0 ⇒ c(t) = c(0) .
The system (EK) is Galilean invariant, i.e. if (ρ(t, x), u(t, x)) solves (EK) then also ρ c (t, x) := ρ(t, x + ct) , u c (t, x) := u(t, x + ct) -c , solves (EK). Then we can always assume that u = ∇φ for some scalar potential φ : T d ν → R. The system (EK) reads

∂ t ρ + div(ρ∇φ) = 0 ∂ t φ + 1 2 |∇φ| 2 + g(ρ) = K(ρ)∆ρ + 1 2 K (ρ)|∇ρ| 2 .
(1.8)

Notice that the average

1 (2π) d ν 1 • • • ν d T d ν ρ(x)dx = m ∈ R , (1.9)
is a constant of motion of (1.8). Notice also that the vector field of (1.8) depends only on Π ⊥ 0 φ (see (1.3)). In view of (1.9) we rewrite ρ m + ρ where ρ is a function with zero average. Then, the system (1.8) (recall also (1.6)) becomes (QHD).

Phase space and notation. In the paper we work with functions belonging to the Sobolev space

H s (T d ν ) := u(x) = 1 (2π) d/2 j∈Z d ν u j e ij•x : u(•) 2 H s (T d ν ) := j∈Z d ν j 2s |u j | 2 < +∞ , (1.10) where j := 1 + |j| 2 for j ∈ Z d ν with Z d ν := (Z/ν 1 ) × • • • × (Z/ν d ). The natural phase space for (QHD) is H s 0 (T d ν ) × Ḣs (T d ν )
where Ḣs (T d ν ) := H s (T d ν )/ ∼ is the homogeneous Sobolev space obtained by the equivalence relation

ψ 1 (x) ∼ ψ 2 (x) if and only if ψ 1 (x) -ψ 2 (x) = c is a constant; H s 0 (T d ν )
is the subspace of H s (T d ν ) of functions with zero average. Despite this fact we prefer to work with a couple of variable

(ρ, φ) ∈ H s 0 (T d ν ) × H s (T d ν )
but at the end we control only the norm

Π ⊥ 0 φ H s (T d ν )
which in fact is the relevant quantity for (QHD). To lighten the notation we shall write

• H s ν to denote • H s (T d ν )
. In the following we will use the notation A B to denote A ≤ CB where C is a positive constant depending on parameters fixed once for all, for instance d and s. We will emphasize by writing q when the constant C depends on some other parameter q.

Ideas of the proof. The general (EK) is a system of quasi-linear equations. The case (QHD), i.e. the system (EK) with the particular choice (1.6), reduces, for small solutions, to a semi-linear equation, more precisely to a nonlinear Schrödinger equation. This is a consequence of the fact that the Madelung transform (introduced for the first time in the seminal work by Madelung [START_REF] Madelung | Quanten theorie in Hydrodynamischer Form Zeit[END_REF]) is well defined for small solutions. In other words one can introduce the new variable ψ := √ m + ρe iφ/ (see Section 2 for details), where = 2 √ k, obtaining the equation

∂ t ψ = i 2 ∆ψ - 1 g(|ψ| 2 )ψ .
Since g(m) = 0, such an equation has an equilibrium point at ψ = √ m. The study of the stability of small solutions for (QHD) is equivalent to the study of the stability of the variable z = ψ -√ m. The equation for the variable z reads

∂ t z = -i |D| 2 ν 2 + mg (m) z -i mg (m) z + f (z),
where f is a smooth function having a zero of order 2 at z = 0, i.e. |f (z)| |z| 2 , and |D| 2 ν is the Fourier multiplier with symbol

|ξ| 2 ν := d i=1 a i |ξ j | 2 , a i := ν 2 i , ∀ ξ ∈ Z d . (1.11)
The aim is to use a Birkhoff normal form/modified energy technique in order to reduce the size of the nonlinearity f (z). To do that, it is convenient to perform some preliminary reductions. First of all we want to eliminate the addendum -i mg (m) z. In other words we want to diagonalize the matrix

L = 2 |D| 2 ν + 1 mg (m) 1 mg (m) 1 mg (m) 2 |D| 2 ν + 1 mg (m)
.

(1.12)

To achieve the diagonalization of this matrix it is necessary to rewrite the equation in a system of coordinates which does not involve the zero mode. We perform this reduction in Section 2.2: we use the gauge invariance of the equation as well as the L 2 norm preservation to eliminate the dynamics of the zero mode. This idea has been introduced for the first time in [START_REF] Faou | Sobolev stability of plane wave solutions to the cubic nonlinear Schrödinger equation on a torus[END_REF]. After the diagonalization of the matrix in (1.12) we end up with a diagonal, quadratic, semi-linear equation with dispersion law

ω(j) := 2 4 |j| 4 ν + mg (m)|j| 2
ν , where j is a vector in Z d \ {0}. At this point we are ready to define a suitable modified energy. Our primary aim is to control the derivative of the H s -norm of the solution

d dt z(t) 2 H s , (1.13)
where z is the variable of the diagonalized system, for the longest time possible. Using the equation, such a quantity may be rewritten as the sum of trilinear expressions in z. We perturb the Sobolev energy by expressions homogeneous of degree at least 3 such that their time derivatives cancel out the main contribution (i.e. the one coming from cubic terms) in (1.13), up to remainders of higher order. In trying to do this small dividers appear, i.e. denominators of the form ±ω(j 1 ) ± ω(j 2 ) ± ω(j 3 ) .

It is fundamental that the perturbations we define is bounded by some power of z H s , with the same s in (1.13), otherwise we obtain an estimate with loss of derivatives. Therefore we need to impose some lower bounds on the small dividers. Here it enters in the game the irrationality of the torus ν. We prove indeed that for almost any ν ∈ [1, 2] d , there exists γ > 0 such that [START_REF] Bambusi | Birkhoff normal form for some nonlinear PDEs[END_REF][START_REF] Bambusi | Birkhoff normal form for partial differential equations with tame modulus[END_REF], that the power of µ 3 is not dangerous if we work in H s with s big enough. Unfortunately we have also a power of the highest frequency µ 1 which represents, in principle, a loss of derivatives. However, this loss of derivatives may be transformed in a loss of length of the lifespan through partition of frequencies, as done for instance in [START_REF] Delort | On long time existence for small solutions of semi-linear Klein-Gordon equations on the torus[END_REF][START_REF] Ionescu | Long-time existence for multi-dimensional periodic water waves[END_REF][START_REF] Feola | Long time solutions for quasi-linear Hamiltonian perturbations of Schrödinger and Klein-Gordon equations on tori[END_REF][START_REF] Bernier | Long-time existence for semi-linear beam equations on irrational tori[END_REF].

| ± ω(j 1 ) ± ω(j 2 ) ± ω(j 3 )| ≥ γ µ d-1 1 log d+1 1 + µ 2 1 µ M (d) 3 , if ±j 1 ± j 2 ± j 3 = 0,
Some comments. As already mentioned, an estimate on small divisors involving only powers of µ 3 is not dangerous. We may obtain such an estimate when the equation is considered on the squared torus T d , using as a parameter the mass m. In this case, indeed, one can obtain better estimates by following the proof in [START_REF] Faou | Sobolev stability of plane wave solutions to the cubic nonlinear Schrödinger equation on a torus[END_REF]. This is a consequence of the fact that the set of differences of eigenvalues is discrete. This is not the case of irrational tori with fixed mass, where the set of eigenvalues is not discrete. Having estimates involving only µ 3 one could actually prove an almost-global stability. More precisely one can prove, for instance, that there exists a zero Lebesgue measure set N ⊂ [1, +∞), such that if m is in [1, +∞) \ N , then for any N ≥ 1 if the initial condition is sufficiently regular (w.r.t. N ) and of size ε sufficiently small (w.r.t. N ) then the solution stays of size ε for a time of order ε -N . The proof follows the lines of classical papers such as [START_REF] Bambusi | Birkhoff normal form for some nonlinear PDEs[END_REF][START_REF] Bambusi | Birkhoff normal form for partial differential equations with tame modulus[END_REF][START_REF] Bambusi | Almost global existence for Hamiltonian semi-linear Klein-Gordon equations with small Cauchy data on Zoll manifolds[END_REF] by using the Hamiltonian structure of the equation. More precisely, the system (QHD) can be written in the form

∂ t ρ φ = X H (ρ, φ) = ∂ φ H(ρ, φ) -∂ ρ H(ρ, φ) , (1.14)
where ∂ denotes the L 2 -gradient and H(ρ, φ) is the Hamiltonian function

H(ρ, φ) = 1 2 T d ν (m + ρ)|∇φ| 2 dx + T d ν 1 2 κ m + ρ |∇ρ| 2 + G(m + ρ) dx (1.15)
where G (ρ) = g(ρ).

We do not know if the solution of (QHD) are globally defined. There are positive answers in the case that the equation is posed on the Euclidean space R d with d ≥ 3, see for instance [START_REF] Audiard | Global Well-Posedness of the Euler-Korteweg System for Small Irrotational Data[END_REF]. Here the dispersive character of the equation is taken into account. For recent developments in this direction see [START_REF] Antonelli | Analysis of acoustic oscillations for a class of hydrodynamic systems describing quantum fluids[END_REF] and reference therein. It is worth mentioning also the scattering result for the Gross-Pitaevsii equation [START_REF] Gustafson | Scattering for the Gross-Pitaevskiiequation[END_REF]. Since we are considering the equation on a compact manifold, the dispersive estimates are not available.

It would be interesting to obtain a long time stability result also for solutions of the general system (EK). In this case the equation may be not recasted as a semi-linear Schrödinger equation. Being a quasi-linear system, we expect that a para-differential approach, in the spirit of [START_REF] Ionescu | Long-time existence for multi-dimensional periodic water waves[END_REF][START_REF] Feola | Long time solutions for quasi-linear Hamiltonian perturbations of Schrödinger and Klein-Gordon equations on tori[END_REF] should be applied. However, in this case, the quasi-linear term is quadratic, hence big. In [START_REF] Ionescu | Long-time existence for multi-dimensional periodic water waves[END_REF][START_REF] Feola | Long time solutions for quasi-linear Hamiltonian perturbations of Schrödinger and Klein-Gordon equations on tori[END_REF] the quasi-linear term is smaller. Therefore new ideas have to be introduced in order to improve the local existence theorem. By using para-compositions (in the spirit of [START_REF] Berti | Almost global solutions of capillary-gravity water waves equations on the circle[END_REF][START_REF] Feola | Long time existence for fully nonlinear NLS with small Cauchy data on the circle[END_REF][START_REF] Feola | A non-linear Egorov theorem and Poincaré-Birkhoff normal forms for quasi-linear pdes on the circle[END_REF]), in the case d = 1, i.e. on the torus T 1 , it is possible to obtain stronger results. This is the argument of a future work of one of us with other collaborators [START_REF] Berti | Long time existence results for the Euler-Korteweg equations[END_REF].

2. FROM (QHD) TO NONLINEAR SCHR ÖDINGER 2.1. Madelung transform. For λ ∈ R + , we define the change of variable (Madelung transform)

ψ := M ψ (ρ, φ) := √ m + ρe iλφ , ψ := M ψ (ρ, φ) := √ m + ρe -iλφ . (M)
Notice that the inverse map has the form

m + ρ = M -1 ρ (ψ, ψ) := |ψ| 2 , φ = M -1 φ (ψ, ψ) := 1 λ arctan -i(ψ -ψ) ψ + ψ . (2.1)
In the following lemma we provide a well-posedness result for the Madelung transform.

Lemma 2.1. Define 1 4λ 2 = κ , := 1 λ = 2 √ κ . (2.
2)

The following holds.

(i) Let s > 1 and

δ := 1 m ρ H s ν + 1 √ κ Π ⊥ 0 φ H s ν σ := Π 0 φ . There is C = C(s) > 1 such that, if C(s)δ ≤ 1, then the function ψ in (M) satisfies ψ - √ me iσ H s ν ≤ 2 √ mδ . (2.3) (ii) Define δ := inf σ∈T ψ - √ me iσ H s ν . There is C = C (s) > 1 such that, if C (s)δ ( √ m) -1 ≤ 1, then the functions ρ, 1 m ρ H s ν + 1 √ κ Π ⊥ 0 φ H s ν ≤ 8 1 √ m δ . (2.4)
Proof. The bound (2.3) follows by (M) and classical estimates on composition operators on Sobolev spaces (see for instance [START_REF] Moser | A rapidly convergent iteration method and non-linear partial differential equations -I[END_REF]). Let us check the (2.4). By the first of (2.1), for any σ ∈ T, we have

ρ H s ν ≤ √ m(ψe -iσ - √ m) H s ν + √ m(ψe iσ - √ m) H s ν + (ψe -iσ - √ m)(ψe iσ - √ m) H s ν (2.5) ≤m 2 √ m ψ - √ me iσ H s ν + ( 1 √ m ψ - √ me iσ H s ν ) 2 . (2.6)
Therefore, by the arbitrariness of σ and using that (

√ m) -1 δ 1, one deduces 1 m ρ H s ν ≤ 3 1 √ m δ .
Moreover we note that

(ψ -ψ)(ψ + ψ) -1 H s ν ≤ 2 1 √ m ψ - √ m H s ν .
Then by the second in (2.1), (2.2), composition estimates on Sobolev spaces and the smallness condition

( √ m) -1 δ 1, one deduces, for any σ ∈ T such that ( √ m) -1 ψ - √ me iσ H s ν 1, that 1 √ κ Π ⊥ 0 φ H s ν + 2 Π ⊥ 0 arctan -i(ψ -ψ) ψ + ψ H s ν = 2 Π ⊥ 0 arctan -i(ψe -iσ -ψe iσ ) ψe -iσ + ψe -iσ H s ν ≤ 5 √ m ψ - √ me iσ H s ν . Therefore the (2.4) follows.
We now rewrite equation (QHD) in the variable (ψ, ψ).

Lemma 2.2. Let (ρ, φ) ∈ H s 0 (T d ν ) × H s (T d ν ) be a solution of (QHD) defined over a time interval [0, T ], T > 0, such that sup t∈[0,T ) 1 m ρ(t, •) H s ν + 1 √ κ Π ⊥ 0 φ(t, •) H s ν ≤ ε (2.7)
for some ε > 0 small enough. Then the function ψ defined in (M) solves

∂ t ψ = -i -2 ∆ψ + 1 g(|ψ| 2 )ψ ψ(0) = m + ρ(0)e iφ(0) . (2.8)
Remark 2.3. We remark that the assumption of Lemma 2.2 can be verified in the same spirit of the local well-posedness results in [START_REF] Feola | Local well-posedness for the Hamiltonian quasi-linear Schrödinger equation on tori[END_REF] and [START_REF] Berti | Local well posedness of the Euler-Korteweg equations on T d[END_REF].

Proof of Lemma 2.2. The smallness condition (2.7) implies that the function ψ is well-defined and satisfies a bound like (2.3). We first note

∇ψ = ( 1 2 √ m + ρ ∇ρ + iλ √ m + ρ∇φ)e iλφ , (2.9
)

1 2λ 2 |∇ψ| 2 = 1 2 1 4λ 2 1 m + ρ |∇ρ| 2 + 1 2 (m + ρ)|∇φ| 2 .
(2.10)

Moreover, using (QHD), (2.2), (M) and that

div(ρ∇φ) = ∇ρ • ∇φ + ρ∆φ ,
we get

∂ t ψ = ie iλφ ∆ρ 4λ √ m + ρ - |∇ρ| 2 8λ(m + ρ) 3 2 + i √ m + ρ∆φ 2 - √ m + ρλ|∇φ| 2 2 + i∇ρ • ∇φ 2 √ m + ρ -iλg(|ψ| 2 )ψ .
(2.11)

On the other hand, by (2.9), we have

i∆ψ = ie λφ ∆ρ 2 √ m + ρ - |∇ρ| 2
4(m + ρ)

3 2 + iλ √ m + ρ∆φ -λ 2 √ m + ρ|∇φ| 2 + iλ∇ρ • ∇φ √ m + ρ . (2.12)
Therefore, by (2.11), (2.12), we deduce

∂ t ψ = i 2λ ∆ψ -iλg(|ψ| 2 )ψ , where 1 λ = , (2.13)
which is the (2.8).

Notice that the (2.8) is an Hamiltonian equation of the form

∂ t ψ = -i∂ ψ H(ψ, ψ) , H(ψ, ψ) = T d ν 2 |∇ψ| 2 + 1 G(|ψ| 2 ) dx , (2.14)
where

∂ ψ = (∂ Re ψ + i∂ Im ψ )/2. The Poisson bracket is {H, G} := -i T d ν ∂ ψ H∂ ψ G -∂ ψ H∂ ψ Gdx . ( 2
ψ 0 = αe -iθ α ∈ [0, +∞) , θ ∈ T ψ j = z j e -iθ j = 0 , (2.16)
which are the polar coordinates for j = 0 and a phase translation for j = 0. Rewriting (2.14) in Fourier coordinates one has

i∂ t ψ j = ∂ ψ j H(ψ, ψ) , j ∈ Z d , (2.17)
where H is defined in (2.14). We define also the zero mean variable

z := (2π) -d 2 j∈Z d \{0} z j e ij•
m := j∈Z d |ψ j | 2 = α 2 + j =0 |z j | 2
is a constant of motion for (2.8). Using (2.16), one can completely recover the real variable α in terms of {z j } j∈Z d \{0} as

α = m - j =0 |z j | 2 .
(2.20)

Note also that the (ρ, φ) variables in (2.1) do not depend on the angular variable θ defined above. This implies that system (QHD) is completely described by the complex variable z. On the other hand, using

∂ ψ j H(ψe iθ , ψe iθ ) = ∂ ψ j H(ψ, ψ)e iθ , one obtains i∂ t α + ∂ t θα = Π 0 (g(|α + z| 2 )(α + z)) i∂ t z j + ∂ t θz j = ∂H ∂ψ j (α + z, α + z) . (2.21)
Taking the real part of the first equation in (2.21) we obtain

∂ t θ = 1 α Π 0 1 g(|α + z| 2 )Re (α + z) = 1 2α ∂ α H(α, z, z) , (2.22)
where (recall (1.11))

H(α, z, z) := 2 T d |D| 2 ν z • zdx + 1 T d G(|α + z| 2 ) dx . ( 2
∂ ψ j H(α + z, α + z) = ∂ z j H(α, z, z) ,
one obtains

i∂ t z j =∂ z j H(α, z, z) - z j 2α ∂ α H(α, z, z) = ∂ z j K m (z, z) , j = 0 , (2.24)
where

K m (z, z) := H(α, z, z) |α= √ m-j =0 |z j | 2 .
We resume the above discussion in the following lemma.

Lemma 2.4. The following holds.

(i) Let s > d 2 and

δ := 1 m ρ H s + 1 √ κ Π ⊥ 0 φ H s , θ := Π 0 φ .
There is

C = C(s) > 1 such that, if C(s)δ ≤ 1, then the function z in (2.18) satisfies z H s ≤ 2 √ mδ . (2.25) (ii) Define δ := z H s .
There is

C = C (s) > 1 such that, if C (s)δ ( √ m) -1 ≤ 1, then the functions ρ, 1 m ρ H s + 1 √ κ Π ⊥ 0 φ H s ≤ 16 1 √ m δ .
(2.26)

(iii) Let (ρ, φ) ∈ H s 0 (T d ν ) × H s (T d ν )
be a solution of (QHD) defined over a time interval [0, T ], T > 0, such that

sup t∈[0,T ) 1 m ρ(t, •) H s + 1 √ κ Π ⊥ 0 φ(t, •) H s ≤ ε
for some ε > 0 small enough. Then the function z ∈ H s 0 (T d ν ) defined in (2.18) solves (2.24). Proof. We note that

z H s = Π ⊥ 0 ψ H s ≤ ψ - √ me iθ Hs (2.3) ≤ 2 √ mδ , (2.27)
which proves (2.25). In order to prove (2.26) we note that ) one can study the system (QHD) near the equilibrium point (ρ, φ) = (0, 0) by studying the complex hamiltonian system

inf σ∈T ψ - √ me iσ H s ≤ ψ - √ me iθ H s = α - √ m + z H s ≤ m -z 2 L 2 - √ m + z H s ≤ 2δ ,
i∂ t z = ∂ z K m (z, z) (2.28)
near the equilibrium z = 0. Note also that the natural phase-space for (2.28) is the complex Sobolev space H s 0 (T d ; C), s ∈ R, of complex Sobolev functions with zero mean. 2.3. Taylor expansion of the Hamiltonian. In order to study the stability of z = 0 for (2.28) it is useful to expand K m at z = 0. We have

K m (z, z) = 2 T d |D| 2 ν z • z dx + 1 T d G m - j =0 |z j | 2 + z 2 dx = (2π) d G(m) + K (2) m (z, z) + N -1 r=3 K (r) m (z, z) + R (N ) (z, z) , (2.29)
where

K (2) m (z, z) = 1 2 T d 2 |D| 2 ν z • z dx + g (m)m T d 1 2 (z + z) 2 dx , (2.30)
for any r = 3,

• • • , N -1, K (r)
m (z, z) is an homogeneous multilinear Hamiltonian function of degree r of the form

K (r) m (z, z) = σ∈{-1,1} r , j∈(Z d \{0}) r r i=1 σ i j i =0 (K (r) m ) σ,j z σ 1 j 1 • • • z σr jr , |(K (r) m ) σ,j | r 1 , and X R (N ) (z) H s s z r-1 H s , ∀ z ∈ B 1 (H s 0 (T d ; C)) .
(2.31)

The vector field of the Hamiltonian in (2.29) has the form (recall (1.14))

∂ t z z = -i∂ z K m i∂ z K m = -i |D| 2 ν 2 + mg (m) mg (m) -mg (m) -|D| 2 ν 2 -mg (m) z z + N -1 r=3 -i∂ z K (r) m i∂ z K (r) m + -i∂ z R (N ) i∂ z R (N ) .
(2.32)

Let us now introduce the 2 × 2 matrix of operators

C := 1 2ω(D)A(D, m) A(D, m) -1 2 mg (m) -1 2 mg (m) A(D, m) , with A(D, m) := ω(D) + 2 |D| 2 ν + 1 2 mg (m)
, and where ω(D) is the Fourier multiplier with symbol

ω(j) := 2 4 |j| 4 ν + mg (m)|j| 2 ν . (2.33)
Notice that, by using (1.2), the matrix C is bounded, invertible and symplectic, with estimates

C ±1 L(H s 0 ×H s 0 , H s 0 ×H s 0 ) ≤ 1 + √ kβ, β := mg (m) k . ( 2
K m (w, w) := K (2) m (w, w) + K (3) m (w, w) + K (≥4) m (w, w) , K (2) m (w, w) := K (2) m C w w := 1 2 T d ω(D)z • zdx , K (3) m (w, w) := K (3) m C w w , K (≥4) m (w, w) := N -1 r=4 K (r) m C w w + R (N ) C w w .
(2.36)

Therefore system (2.32) becomes

∂ t w = -iω(D)w -i∂ w K (3) m (w, w) -i∂ w K (≥4) m (w, w) .
(2.37)

SMALL DIVISORS

As explained in the introduction we shall study the long time behaviour of solutions of (2.37) by means of Birkhoff normal form approach. Therefore we have to provide suitable non resonance conditions among linear frequencies of oscillations ω(j) in (2.33). This is actually the aim of this section.

Let a = (a 1 , . . . ,

a d ) = (ν 2 1 , . . . , ν 2 d ) ∈ (1, 4) d , d = 2, 3. If j ∈ Z d \ {0} we define |j| 2 a = d k=1 a k j 2 k . (3.1)
We consider the dispersion relation

ω(j) := k|j| 4 a + mg (m)|j| 2 a , (3.2)
we note that ω(j)

= √ k |j| 2 a + β 2 -β 2 8 1 |j| 2 a +O( β 3 |j| 4 a
) for any j big enough with respect to β := mg (m) k .

Throughout this section we assume, without loss of generality,

|j 1 | a ≥ |j 2 | a ≥ |j 3 | a > 0,
for any j i in Z d , moreover, in order to lighten the notation, we adopt the convention ω i := ω(j i ) for any i = 1, 2, 3. The main result is the following. Proposition 3.1. (Measure estimates). There exists a full Lebesgue measure set A ⊂ (1, 4) d such that for any a ∈ A there exists γ > 0 such that the following holds true. If σ 1 j 1 + σ 2 j 2 + σ 3 j 3 = 0, σ i ∈ {±1} we have the estimate The proof of this proposition is divided in several steps and it is postponed to the end of the section. The main ingredient is the following standard proposition which follows the lines of [START_REF] Bambusi | Birkhoff normal form for some nonlinear PDEs[END_REF][START_REF] Bernier | Long-time existence for semi-linear beam equations on irrational tori[END_REF]. Here we give weak lower bounds of the small divisors, these estimates will be improved later.

k -1 2 σ 3 ω 3 + σ 2 ω 2 + σ 1 ω 1 d γ |j 1 | d-1 log d+1 (1+|j 1 | 2 )|j 3 | M (d) , if σ 1 σ 2 = -1 1, if σ 1 σ 2 = 1 . (3.3) for any |j 1 | a ≥ |j 2 | a ≥ |j 3 | a , j i ∈ Z
Proposition 3.2. Consider I and J two bounded intervals of R + \ {0}; r ≥ 2 and j 1 , . . . , j r ∈ Z d such that j i = ±j k if i = k, n 1 , . . . n r ∈ Z \ {0} and h : J d-1 → R measurable. Then for any γ > 0 we have

µ (p, b) ∈ I × J d-1 : h(b) + r k=1 n k |j k | 4 (1,b) + p|j k | 2 (1,b) ≤ γ I,J,d,r,n γ 1 2r ( j 1 • • • j r ) 1 r , with (1, b) = (1, b 1 , . . . , b d-1 ) ∈ R d .
Remark 3.3. We shall apply this general proposition only in the case r = 3, however we preferred to write it in general for possible future applications.

Proof of Prop. 3.2. For simplicity in the proof we assume

|j 1 | (1,b) > . . . > |j r | (1,b) .
Since by assumption we have j i = j k for any i = k then one could easily prove that for any η > 0 (later it will be chosen in function of γ) we have

µ(P i,k η ) < ηµ(J d-2 ) , P i,k η := {b ∈ J d-1 : |j i | 2 (1,b) -|j k | 2 (1,b) < η} . We define P η = ∪ i =k P i,k
η , and

B γ := (p, b) ∈ I × J d-1 : h(b) + r k=1 n k |j k | 4 (1,b) + p|j k | 2 (1,b) ≤ γ , then we have µ(B γ ) ≤ µ(B γ ∩ P η ) + µ(B γ ∩ (P η ) c ) ≤ µ(I)µ(P η ) + µ(J d-1 ) sup i =k b / ∈Pη µ {p ∈ I : h(b) + r k=1 n k |j k | 4 (1,b) + p|j k | 2 (1,b) ≤ γ} r µ(I)µ(J d-2 )η + µ(J d-1 ) sup i =k b / ∈Pη µ {p ∈ I : h(b) + r k=1 n k |j k | 4 (1,b) + p|j k | 2 (1,b) ≤ γ} .
We have to estimate from above the measure of the last set. We define the function

g(p) := h(b) + r k=1 n k |j k | 4 (1,b) + p|j k | 2 (1,b) .
For any ≥ 1 we have

d dp g(p) = c r k=1 n k |j k | (1,b) (p + |j k | 2 (1,b) ) 1 2 -, c := i=1 ( 1 2 -i).
Therefore we can write the system of equations

   c -1 1 ∂ 1 p g(p) . . . c -1 r ∂ r p g(p)    =    (p + |j 1 | 2 (1,b) ) 0 . . . (p + |j r | 2 (1,b) ) 0 . . . . . . . . . (p + |j 1 | 2 (1,b) ) 1-r . . . (p + |j r | 2 (1,b) ) 1-r       n 1 |j 1 | (1,b) (p + |j 1 | 2 (1,b) ) -1/2 . . . n r |j r | (1,b) (p + |j r | 2 (1,b) ) -1/2    .
We denote by V the Vandermonde matrix above. We have that V is invertible since

| det(V )| = 1≤i<k≤r 1 p + |j i | 2 (1,b) - 1 p + |j k | 2 (1,b) ≥ 1≤i<k≤r |j i | 2 (1,b) -|j k | 2 (1,b) (p + |j i | 2 (1,b) )(p + |j k | 2 (1,b) ) 1≤k≤r η (p + |j k | 2 (1,b) ) 2 η r 1 j 1 2 • • • j r 2 ,
where in the penultimate passage we have used that b / ∈ P η and

|j i | (1,b) ≤ |j k | (1,b) if i > k. Therefore we have r max =1 |c ∂ p g(p)| r | det(V )| r max =1 n |j | (1,b) (p + |j |) -1 2 (1,b) r,n η r |j 1 | 1/2 (1,b) j 1 2 • • • j r 2 r,n η r j 1 2 • • • j r 2 .
At this point we are ready to use Lemma 7 in appendix A of the paper [START_REF] Procesi | Reducible quasi-periodic solutions of the non linear Schrödinger equation Boll[END_REF], we obtain

µ {p ∈ I : h(b) + r k=1 |j k | 4 (1,b) + p|j k | 2 (1,b) ≤ γ} ≤ ( γ j 1 2 . . . j r 2 η r 1 r .
Summarizing we obtained

µ(B γ ) I,J,d,r,n η + η -1 γ 1 r ( j 1 2 . . . j r 2 ) 1 r ,
we may optimize by choosing η = γ

1 2r ( j 1 • • • j r ) 1 r
and we obtain the thesis.

As a consequence of the preceding proposition we have the following.

Corollary 3.4. Let r ≥ 1, consider j 1 , . . . , j r ∈ Z d such that j k = j i if i = k and n 1 , . . . , n k ∈ Z \ {0}. For any γ > 0 we have µ a ∈ (1, 4) :

r i=1 n i k|j i | 4 a + mg (m)|j i | 2 a ≤ γ d,r,n γ √ k 1 2r ( j 1 . . . j r ) 1 r .
Proof. We write

r i=1 n i k|j i | 4 a + mg (m)|j i | 2 a = √ ka 1 r i=1 n i |j i | 4 (1,b) + β a 1 |j i | 2 (1,b) ,
where we have set

β := mg (m) k , b := ( a 2 a 1 , . . . , a d a 1
).

(3.4)

The map (a 1 , . . . ,

a d ) → ( 1 a 1 , b) is invertible onto its image, which is contained in ( 1 4 , 1)×( 1 4 , 4) d-1
. The determinant of its inverse is bounded by a constant depending only on d. Therefore the result follows by applying Prop. 3.2 and the change of coordinates (a 1 , . . . , a d ) → ( 1 a 1 , b). Owing to the corollary above we may reduce in the following to the study of the small dividers when we have 2 frequencies much larger then the other.

Lemma 3.5. Consider Λ := √ k|j 1 | 2 a - √ k|j 2 | 2 a -ω 3 and β defined in (3.4). If there exists i ∈ {1, . . . , d} such that |j 3,i | 1 + β |j 3 | 2 a ≤ 1 2 |j 1,i + j 2,i | , (3.5)
then for any γ > 0 we have

µ a ∈ (1, 4) d : | Λ| ≤ γ ≤ 2γ √ k|j 1,i + j 2,i | .
Proof. We give a lower bound for the derivative of the function Λ with respect to the parameter a i .

|∂ a i Λ| ≥ √ k |j 3,i (j 1,i + j 2,i )| -j 2 3,i 1 + β |j 3 | 2 a ≥ √ k 1 2 |j 3,i ||j 1,i + j 2,i | ≥ √ k 1 2 |j 1,i + j 2,i |.
Therefore a i → Λ is a diffeomorphism and applying this change of variable we get the thesis.

Proposition 3.6. There exists a set of full Lebesgue measure A 3 ⊂ (1, 4) d such that for any a in A 3 there exists γ > 0 such that

σω 3 + ω 2 -ω 1 ≥ √ kγ |j 1 | d-1 log d+1 (1 + |j 1 | 2 )|j 3 | d+1 , for any σ ∈ ±1, for any j 1 , j 2 , j 3 in Z d satisfying |j 1 | a > |j 2 | a ≥ |j 3 | a , the momentum condition σj 3 + j 2 -j 1 = 0 and J(j 1 , β) = min ||j 1 | 2 -4d 2 β| 2d , min γ 4β 2 1 d+2 , γ 2β 3 1 d+1 |j 1 | 4-d log(1 + |j 1 |) d+1 1 d+2 > |j 3 | , (3.6)
where β is defined in (3.4).

Proof. We suppose σ = 1, we set Λ := ω 1 -ω 2 -ω 3 and

L(γ) := √ kγ |j 3 | d+1 |j 1 | d-1 log d+1 (1 + |j 1 |) .
From the first condition in (3.6) we deduce that β/|j 1 | 2 < 1, therefore, by Taylor expanding the (3.2), we obtain

Λ = √ k |j 1 | 2 a -|j 2 | 2 a + β 2 8 |j 1 | 2 a -|j 2 | 2 a |j 2 | 2 a |j 1 | 2 a + R -ω 3 , (3.7)
where |R| ≤ 1 8

β 3 |j 2 | 4 a . We define Λ := √ k|j 1 | 2 a - √ k|j 2 | 2
a -ω 3 and the following good sets

G γ := {a ∈ (1, 4) d : |Λ| > L(γ), ∀j 1 , j 3 ∈ Z d } , Gγ := {a ∈ (1, 4) d : Λ > 3L(γ) , ∀j 1 , j 3 ∈ Z d } .
We claim that, thanks to (3.6), we have the inclusion Gγ ⊂ G γ . First of all we have

Λ ≥ ω 3 + √ k|j 2 | 2 a - √ k|j 1 | 2 a - √ k β 2 8 |j 1 | 2 a -|j 2 | 2 a |j 1 | 2 a |j 2 | 2 a - √ k|R| . (3.8)
From the momentum condition j 1 -j 2 = j 3 and the ordering

|j 1 | a > |j 2 | a ≥ |j 3 | a we have that |j 1 | a ≤ 2|j 2 | a , which implies |j 1 | 2 a -|j 2 | 2 a |j 1 | 2 a |j 2 | 2 a = d k=1 a k j 3,k (j 1,k + j 2,k ) |j 1 | 2 a |j 2 | 2 a ≤ 2 |j 3 | a |j 1 | a |j 1 | 2 a |j 2 | 2 a ≤ 2 |j 3 | a |j 1 | a |j 2 | 2 a ≤ 32 |j 3 | |j 1 | 3 , (3.9)
where we used

| • | < | • | a < 4| • |.
Therefore from (3.6), more precisely from

( γ 4β 2 ) 1 d+2 |j 1 | 4-d log(1 + |j 1 |) d+1 1 d+2 > |j 3 |, we deduce that √ k β 2 8 |j 1 | 2 a -|j 2 | 2 a |j 1 | 2 a |j 2 | 2 a < L.
Analogously one proves that √ k|R| < L. We have eventually proved that Gγ ⊂ G γ using (3.8). We define the bad sets Bγ := (1, 4) d \ Gγ ⊃ B γ := (1, 4) d \ G γ and we prove that the Lebesgue measure of ∩ γ Bγ equals to zero, this implies the thesis.

We want to apply Lemma 3.5 with γ L. We know that there exists i ∈ {1, . . . , d} such that d|j 1,i | ≥ |j 1 |. We claim that, thanks to (3.6), we satisfy condition (3.5) for the same index i. Let us suppose by contradiction that

|j 3,i | 1 + β |j 3 | 2 a > 1 2 |j 1,i + j 2,i | = 1 2 |2j 1,i -j 3,1 | ≥ |j 1,i | -1 2 |j 3,i | > |j 1 | d -1 2 |j 3,i | , from which we obtain |j 1 | ≤ 2d|j 3 | 1 + β/|j 3 | 2 a .
Taking the squares we get

|j 1 | 2 ≤ 4d 2 |j 3 | 2 + 4d 2 β |j 3 | 2 |j 3 | 2 a , which, recalling that | • | < | • | a < 4| • |, contradicts (3.6).
Therefore, by using Lemma 3.5, we have

µ Bγ = µ a ∈ (1, 4) d |∃j 1 , j 3 ∈ Z d : | Λ| ≤ √ kγ|j 3 | -d-1 |j 1 | 1-d log(|j 1 |) -d-1 ≤ j 3 ∈Z d 1 |j 3 | d+1 j 1 ∈Z d γ |j 1 | d-1 |j 1,i | log(|j 1 |) d+1 d γ .
This implies that meas(∩ γ B γ ) = 0, hence we can set A 3 = ∪ γ G γ .

We are now in position to prove Prop. 3.1.

Proof of Prop. 3.1. The case σ 1 σ 2 = 1 is trivial, we give the proof if σ 1 σ 2 = -1. From Prop. 3.6 we know that there exists a full Lebesgue measure set A 3 and γ > 0 such that the statement is proven if |j 3 | ≤ J(j 1 , γ). Let us now assume |j 3 | > J(j 1 , γ). Let us define

B γ := j 1 ,j 3 ∈Z d a ∈ (1, 4) d : |σ 3 ω 3 + ω 2 -ω 1 | ≤ √ k γ |j 3 | M (d) ,
where γ will be chosen in function of γ and M (d) big enough w.r.t. d. (3.6)) and Corollary 3.4 with r = 3, we have

Let us set

p := ( M (d) 6 -d -1) 1 d+2 suppose for the moment γ/4β 2 1 d+2 ≤ γ/2β 3 1 d+1 . From |j 3 | > J(j 1 , β) (see
µ(B γ ) d ( √ k) -1 6 j 1 ,j 3 ∈Z d γ 1 6 |j 3 | M (d)/6 j 1 d ( √ k) -1 6 γ1/6 γ -p (4β 2 ) p j 1 ∈Z d log p(d+1) (1 + |j 1 |) |j 1 | (4-d)p-1 j 3 |j 3 | -d-1 .
If the exponent M (d) (and hence p) is chosen large enough we get the summability in the r.h.s. of the inequality above. We now choose γ1/6 γ -p = γ m , we eventually obtain µ(B γ )

γ m . If γ/4β 2 1 d+2 > γ/2β
:= A 3 ∩ (∪ γ>0 B c γ ).

ENERGY ESTIMATES

In this section we construct a modified energy for the Hamiltonian K m in (2.36). We first need some convenient notation.

G 3 (w) = σ i ∈{-1,1}, j i ∈Z d \{0} σ 1 j 1 +σ 2 j 2 +σ 3 j 3 =0 (G 3 ) σ,j w σ 1 j 1 w σ 2 j 2 w σ 3 j 3 , (G 3 ) σ,j ∈ C , σ := (σ 1 , σ 2 , σ 3 ) j := (j 1 , j 2 , j 3 ) (4.1)
with symmetric coefficients (G 3 ) σ,j (i.e. for any ρ ∈ S 3 one has (G 3 ) σ,j = (G 3 ) σ•ρ,j•ρ ) and where we denoted w σ j := w j , if σ = + , w σ j := w j , if σ = -.

The Hamiltonian in (2.36) has the form (see (2.33))

K m := K (2) m + K (3) m + K (≥4) m , K (2) m = j∈Z d \{0} ω(j)w j w j , (4.2)
where K

m is a trilinear Hamiltonian in L 3 with coefficients satisfying

|( K (3) m ) σ,j | 1 , ∀ σ ∈ {-1, +1} 3 , j ∈ (Z d ) 3 \ {0} , (4.3)
and where

K (≥4) m satisfies for any s > d/2 X K (≥4) m (w) H s s w 3 H s , if w H s < 1 . (4.4)
The main result of this section is the following.

Proposition 4.3. Let A and M given by Proposition 3.1. Consider a ∈ A. For any N > 1 and any s ≥ s0 , for some s0 = s0 (M ) > 0, there exist ε 0 s,δ log -d-1 (1 + N ) and a trilinear function E 3 in the class L 3 such that the following holds:

• the coefficients (E 3 ) σ,j satisfies

|(E 3 ) σ,j | s N d-2 log d+1 (1 + N )µ 3 (j) M +1 µ 1 (j) 2s , (4.5
)

for σ ∈ {-1, 1} 3 , j ∈ (Z d ) 3 \ {0}; • for any w in the ball of radius ε 0 of H s 0 (T d ; C) one has |{N s + E 3 , K m }| s N d-2 log d+1 (1 + N ) w 4 H s + N -1 w 3 H s . (4.6)
where N s is defined as

N s (w) := w 2 H s = j∈Z d j 2s |w j | 2 , (4.7)
and K m in (4.2).

In subsection 4.1 we study some properties of the Hamiltonians in L 3 of Def. 4.2. Then in subsection 4.2 we give the proof of Proposition 4.3. Finally, in subsection 4.3, we conclude the proof of the main theorem. 4.1. Trilinear Hamiltonians. We now recall some properties of trilinear Hamiltonians introduced in Definition 4.2. We first give some further definitions.

Definition 4.4. Let N ∈ R with N ≥ 1. (i) If G 3 ∈ L 3 then G >N 3 denotes the element of L 3 defined by (G >N 3 ) σ,j := (G 3 ) σ,j , if µ 2 (j) > N , 0 , else . (4.8) We set G ≤N 3 := G 3 -G >N 3 . (ii) We define G (+1) 3 ∈ L 3 by (G (+1) 3
) σ,j := (G 3 ) σ,j , when ∃i, p = 1, 2, 3, s.t.

µ 1 (j) = |j i | , µ 2 (j) = |j p | and σ i σ p = +1 . We define G (-1) 3 := G 3 -G (+1) 3 .
Combining Lemmata 4.5 and 4.6 we deduce

|{N s , ( K (3,-1) m ) (>N ) }(w)| s N -1 w 3 H s , (4.12)
for s large enough with respect to M . We now prove the estimate (4.6). We have

{N s + E 3 , K m } (4.2) = {N s + E 3 , K (2) m + K (3) m + K (≥4) m } (4.13) = {N s , K (2) m } (4.14)
+ {N s , K (3) m } + {E 3 , K (2) m } (4.15) + {E 3 , K (3) m + K (≥4) m } + {N s , K (≥4) m }. ( 4
t w = ∂ w K m (w, w) w(0) = w 0 ∈ H s 0 (T d ; C) , (4.17)
for some s > 0 large. We shall prove the following.

Lemma 4.8. (Main bootstrap) Let s 0 = s 0 (d) given by Proposition 4.3. For any s ≥ s 0 , there exists ε 0 = ε 0 (s) such that the following holds. Let w(t, x) be a solution of (4.17) with t ∈ [0, T ), T > 0 and initial condition w(0 Proof of Lemma 4.8. We study how the equivalent energy norm E s (w) defined in (4.20) evolves along the flow of (4.17 Since ε can be chosen small with respect to s, with this choices we get for times T as in (4.25). This implies the thesis.

, x) = w 0 (x) ∈ H s 0 (T d ; C). For any ε ∈ (0, ε 0 ) if w 0 H s ≤ ε , sup t∈[0,T ) w(t) H s ≤ 2ε , T ≤ ε -1-1 d-1 log -d+1 2 1 + ε 1 1-d , ( 4
Proof of Theorem 1.1. In the same spirit of [START_REF] Feola | Local well-posedness for the Hamiltonian quasi-linear Schrödinger equation on tori[END_REF], [START_REF] Berti | Local well posedness of the Euler-Korteweg equations on T d[END_REF] we have that for any initial condition (ρ 0 , φ 0 ) as in (1.4)

.15) 2 . 2 .

 22 Elimination of the zero mode. In the following it would be convenient to rescale the space variables x ∈ T d ν ν • x with x ∈ T d and work with functions belonging to the Sobolev space H s (T d ) := H s (T d 1 ), i.e. the Sobolev space in (1.10) with ν = (1, . . . , 1). By using the notation ψ = (2π) -d 2 j∈Z d ψ j e ij•x , we introduce the set of variables

  d and where M (d) is a constant depending only on d.

Definition 4 . 1 .

 41 If j ∈ (Z d ) r for some r ≥ k then µ k (j) denotes the k st largest number among |j 1 |, . . . , |j r | (multiplicities being taken into account).

Definition 4 . 2 .

 42 (Formal Hamiltonians). We denote by L 3 the set of Hamiltonian having homogeneity 3 and such that they may be written in the form

  ). Notice that∂ t E s (w) = -{E s , H}(w) .Therefore, for any t ∈ [0, T ], we have thatT 0 ∂ t E s (w) dt

( 4 . 0 ∂

 40 6),(4.18)s T N d-2 log (d+1) (1 + N )ε 4 + N -1 ε 3 . Let 0 < α and set N := ε -α . Hence we have T 0 ∂ t E s (w) dt s ε 2 T ε 2-α(d-2) log (d+1) (1 + ε -α ) + ε 1+α . (4.23) We choose α > 0 such that 2 -α(d -2) = 1 + α , i.e. t E s (w) dt s ε 2 T ε d d-1 log d+1 (1 + ε -α ) .

T 0 ∂ 0 ∂

 00 t E s (w) dt ≤ ε 2 /4as long as T ≤ ε -d/(d-1) log -the equivalence of norms (4.21) and choosing c 0 > 0 small enough, we havew(t) 2 H s ≤ (1 + c 0 )E 0 (w(t)) ≤ (1 + c 0 ) E s (w(0)) + T t E s (w) dt ≤ (1 + c 0 ) 2 ε 2 + (1 + c 0 )ε 2 /4 ≤ ε 2 3/2 ,

  we denoted by M (d) a positive constant depending on the dimension d and µ i the i-st largest integer among |j 1 |, |j 2 | and |j 3 |. It is nowadays well known, see for instance

  3 1 d+1 one can reason similarly. The wanted set of full Lebesgue measure is therefore obtained by choosing A

  Proof of the main result. Consider the Hamiltonian K m (w, w) in (4.2) and the associated Cauchy problem i∂

	.16)
	We study each summand separately. First of all note that (recall (4.7), (4.2)) the term (4.14)
	vanishes. By (4.4), (4.5) and Lemma 4.5-(ii) we obtain
	|(4.16)| s N d-2 log d+1 (1 + N ) w 4 H s .
	Moreover, by (4.11), (4.12), we deduce
	|(4.15)| s N -1 w 3 H s .
	The discussion above implies the bound (4.6).
	4.3.

  First of all we show that the energy N s + E 3 constructed by Proposition 4.3 provides an equivalent Sobolev norm. Lemma 4.9. (Equivalence of the energy norm). Let N ≥ 1. Let w(t, x) as in (4.18) with s 1 large enough. Then, for any 0 < c 0 < 1, there exists C = C(s, d, c 0 ) > 0 such that, if we have the smallness conditionεCN d-2 log (d+1) (1 + N ) ≤ 1 , (4.19) the following holds true. Define E s (w) := (N s + E 3 )(w) (4.20) with N s is in (4.7), E 3 is given by Proposition 4.3. We have1/(1 + c 0 )E s (w) ≤ w 2 H s ≤ (1 + c 0 )E s (w) , ∀t ∈ [0, T ] . (4.21)Proof. Fix c 0 > 0. By (4.5) and Lemma 4.5, we deduce|E 3 (w)| ≤ C w 3 H s N d-2 log (d+1) (1 + N ) ,(4.22)for some C > 0 depending on s. Then, recalling (4.20), we get|E s (w)| ≤ w 2 H s (1 + C w H s N d-2 log (d+1) (1 + N )) H s (1 + c 0 ), where we have chosen C in (4.19) large enough. This implies the first inequality in (4.21). On the other hand, using (4.22) and (4.18), we havew 2 H s ≤ E s (w) + CN d-2 log (d+1) (1 + N )ε w 2 H s .Then, taking C in (4.[START_REF] Moser | A rapidly convergent iteration method and non-linear partial differential equations -I[END_REF]) large enough, we obtain the second inequality in (4.21).

	(4.19) ≤ w 2

.

[START_REF] Madelung | Quanten theorie in Hydrodynamischer Form Zeit[END_REF]

then we have the improved bound sup t∈[0,T ) w(t) H s ≤ 3/2ε .

  there exists a solution of (QHD) satisfying •) H s ≤ 2ε for some T > 0 possibly small. The result follows by Lemma 4.8. By Lemma 2.4 and estimates (2.34) we deduce that the function w solving the equation (2.37) is defined over the time interval [0, T ] and satisfies sup As long as ν ∈ [1, 2] d (defined as at the beginning of section 3) belongs to the full Lebesgue measure set given by Proposition 3.1, we can apply Proposition 4.3 if ε is small enough. Then by Lemma 4.8 and by a standard bootstrap argument we deduce that the solution w(t) is defined for t ∈ [0, T ε ], T ε as in (1.5), and Using again Lemma 2.4 and (2.34) one can deduce the bound (1.5). Hence the thesis follows.

	sup t∈[0,T )	1 m 0 φ(t, t∈[0,T ] ρ(t, •) H s + 1 √ k Π ⊥ w(t) H s ≤ 4 √ √ m(1 +	kβ)ε .
		sup	w(t) H s ≤ 8 √ m(1 +	√ kβ)ε .
		t∈[0,Tε]

Consider the quadratic Hamiltonian K

(2) m in (4.2). Given a trilinear Hamiltonian G 3 in L 3 we define the adjoint action ad K (2) m G 3 := { K (2) m , G 3 } (see (2.15)) as the Hamiltonian in L 3 with coefficients

The following lemma is the counterpart of Lemma 3.5 in [START_REF] Bernier | Long-time existence for semi-linear beam equations on irrational tori[END_REF]. We omit its proof.

Assume that the coefficients (G i 3 ) σ,j satisfy (recall Def. 4.1)

for some

In particular one has the following estimates:

for any w ∈ H s 0 (T d ; C). (ii) (Poisson bracket). The Poisson bracket between G 1 3 and G 2 3 satisfies the estimate

We have the following result. (recall Definition 4.2). Assume also that the coefficients of G 3 satisfy

for some β > 0, C > 0 and q ≥ 0. We have that the Hamiltonian Q

3 }, η ∈ {-1, 1}, belongs to the class L 3 and has coefficients satisfying

Proof. One can reason as in the proof of Lemma 4. ), provides a decay property of the coefficients w.r.t. the highest index. This implies (by Lemma 4.5-(ii)) a smoothing property of the Hamiltonian {N s , G


Exploring multi-programming for quantum algorithms

Siyuan Niu and Aida Todri-Sanial LIRMM, University of Montpellier, CNRS, 34090, Montpellier, France Today's quantum computing is performed on Noisy Intermediate-Scale Quantum (NISQ) hardware. The NISQ era allows for a limited number of qubits with several physical limitations. As NISQ devices have unavoidable noisy quantum operations, we can execute only small circuits on a quantum machine to get reliable results. But, this leads to the quantum hardware under-utilization issue. With the increase of hardware qubit number and the improvement of error rates, it becomes possible to execute multiple circuits on a quantum chip simultaneously, enhancing the throughput and the utilization of NISQ hardware. Thus, we address a timely problem on how to execute reliably multiple quantum programs on given quantum hardware. This is the first attempt to propose a complete multi-programming process flow for executing an optimal number of workloads in parallel, ensuring output fidelity by analyzing the hardware limitations.

The multi-programming mechanism can be particularly interesting for some specific quantum algorithms like VQE and VQLS, which can prepare several ansatz states simultaneously by executing them in parallel in one quantum processor. The proposed flow is also general enough to be applied to other quantum circuits requiring parallel sub-problem executions regardless of applications or algorithms.

Our main contributions can be listed as follows.

1). We introduce a parallelism manager that can optimally select the number of circuits being executed on the quantum hardware simultaneously.

2). We propose two different qubit partition algorithms, which are part of a hardware-aware multiprogramming compiler. One is greedy and aims to provide the optimal partitions for multiple circuits. The other one is a heuristic that can give nearly optimal results and significantly reduce the time complexity. 3). We use the Simultaneous Randomized Benchmarking protocol to characterize the crosstalk properties and consider them in the qubit partition process to avoid the crosstalk effects during simultaneous executions. 4). We enhance the mapping transition algorithm included in the scheduler by reducing the number of inserted gates to make quantum circuits executable on real quantum devices. 5). We execute circuits of different sizes simultaneously on IBM quantum hardware and compare it with other methods found in the literature. To the best of our knowledge, this is the first attempt to report less than 10% fidelity loss when executing simultaneously four independent circuits on a quantum chip. We also apply multi-programming to the VQE algorithm to show its potential interest.

We execute four different quantum circuits on IBM Q 65 Manhattan and obtain circuit fidelity improvement by 51.8% and the reduction of the additional gate by 53.9 % compared to state of the art. We also investigate our multi-programming approach on the VQE algorithm and perform experiments on IBM Q 65 Manhattan to estimate deuteron's ground state energy. We execute four optimizations (eight circuits) simultaneously with an error rate of 7% compared to the theory result. The number of required circuits is reduced by eight times compared to the Pauli grouping measurement method and 16 times compared to the naive measurement approach. It demonstrates the possibility of enabling a multi-programming mechanism for quantum algorithms.


Palas Kumar Farsoiya

Stéphane Popinet

Luc Deike

Bubble-mediated transfer of dilute gas in turbulence

Bubble-mediated gas exchange in turbulent flow is critical in bubble column chemical reactors as well as for ocean-atmosphere gas exchange related to air entrained by breaking waves. Understanding the transfer rate from a single bubble in turbulence at large Péclet numbers (defined as the ratio between the rate of advection and diffusion of gas) regimes is important as it can be used for improving the models on a larger scale. We characterize the mass transfer of dilute gases from a single bubble in a homogeneous isotropic turbulent flow in the limit of negligible bubble volume variations. We show that the mass transfer occurs within a thin diffusive boundary layer at the bubble-liquid interface, whose thickness decreases with an increase in turbulent Péclet number, Pe. We propose a suitable time scale θ for Higbie (1935) penetration theory, θ = d 0 /ũ, based on d 0 the bubble diameter and ũ a characteristic turbulent velocity, here ũ = √ 3 u rms , where u rms is the large-scale turbulence fluctuations. This leads to a non-dimensional transfer rate Sh = 2(3) 1/4 Pe/π, from the bubble in the isotropic turbulent flow. The theoretical prediction is verified by direct numerical simulations of mass transfer of dilute gas from a bubble in homogeneous and isotropic turbulence, and very good agreement is observed as long as the thin boundary layer is properly resolved.

Introduction

The evolution of bubbles in turbulence has multiple applications in environmental and industrial contexts, from exchange at the ocean-atmosphere interface [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF][START_REF] Deike | Gas transfer by breaking waves[END_REF], as bubble-mediated gas exchange accounts for a significant part of carbon dioxide (CO 2 ) uptake by the ocean [START_REF] Reichl | Contribution of sea-state dependent bubbles to air-sea carbon dioxide fluxes[END_REF], to bubble column chemical reactors [START_REF] Risso | Agitation, mixing, and transfers induced by bubbles[END_REF], while providing an effective pathway to bring oxygen in lakes [START_REF] Karn | Mass transfer studies across ventilated hydrofoils: A step towards hydroturbine aeration[END_REF]. In an engineering context, bubble-mediated gas exchange is controlled by the level of agitation in the flow, which is itself influenced by the presence of bubbles, and in turn affects the bubble velocity statistics controlling the gas transfer from individual bubbles [START_REF] Risso | Agitation, mixing, and transfers induced by bubbles[END_REF][START_REF] Mathai | Bubble and buoyant particle laden turbulent flows[END_REF].

While empirical formulae (e.g. [START_REF] Karn | Mass transfer studies across ventilated hydrofoils: A step towards hydroturbine aeration[END_REF]; Colombet et al. (2015)) have been proposed to describe the gas transfer in bubble swarms, their general applicability and theoretical foundation remain an active research topic. Bubble-mediated gas transfer models in the context of ocean-atmosphere interaction often use formulae based on gas transfer theory [START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF] based on bubble rising in quiescent flow [START_REF] Woolf | Bubbles and the air-sea exchange of gases in nearsaturation conditions[END_REF][START_REF] Keeling | On the role of large bubbles in air-sea gas exchange and supersaturation in the ocean[END_REF][START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF][START_REF] Deike | Gas transfer by breaking waves[END_REF]. The diffusive gas transfer by a single bubble rising in a quiescent flow has been described theoretically by [START_REF] Boussinesq | Calcul du poivoir refroidissant des courants fluide[END_REF] and [START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF], within the assumption of negligible variation in bubble volume. The non-dimensional transfer rate from a single bubble is the Sherwood number Sh = k L d 0 /D l , where k L is the transfer rate, d 0 the bubble size and D l the gas diffusivity in the liquid, and reads, Sh = 2 √ π √ Pe where Pe = U d 0 /D l is the bubble Péclet number and U the bubble rise velocity. Several studies dealing with mass transfer from bubbles rising in a quiescent liquid have developed numerical techniques to resolve interphase mass transfer and validated their results against diffusive mass transfer theory [START_REF] Haroun | Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film[END_REF][START_REF] Marschall | Numerical simulation of species transfer across fluid interfaces in free-surface flows using openfoam[END_REF][START_REF] Deising | A unified single-field model framework for volume-of-fluid simulations of interfacial species transfer applied to bubbly flows[END_REF][START_REF] Fleckenstein | A volume-of-fluid-based numerical method for multi-component mass transfer with local volume changes[END_REF][START_REF] Deising | Direct numerical simulation of mass transfer in bubbly flows[END_REF]. Direct numerical simulations (DNS) of mass transfer of dilute gas from rising bubble swarms by [START_REF] Roghair | Direct numerical simulations of hydrodynamics and mass transfer in dense bubbly flows[END_REF] concluded that the transfer rate increases as the gas hold-up (ratio of gas volume to the total volume) increases.

The gas transfer at the interface between two fluids involving a turbulent flow has been described as presenting two regimes depending on the turbulent Reynolds number [START_REF] Theofanous | Turbulent mass transfer at free, gasliquid interfaces, with applications to open-channel, bubble and jet flows[END_REF]. At low-Reynolds number, the transfer rate can be written as k L /u 0 ∝ Sc -1/2 Re -1/2 , where Sc = ν/D l is the Schmidt number (ratio of kinematic viscosity over mass diffusivity), u 0 a characteristic turbulent velocity and Re a turbulence Reynolds number. The high-Reynolds number regime can be described through the action of the smallest eddies at the Kolmogorov scale enhancing the transfer, leading to k L ∼ Sc -1/2 (ν ) 1/4 , where is the turbulence dissipation rate. This can be expressed in terms of the turbulence Reynolds number and reads k L /u 0 ∝ Sc -1/2 Re -1/4 (Theofanous et al. 1976), further discussed by [START_REF] Magnaudet | Turbulent mass transfer through a flat shear-free surface[END_REF]; [START_REF] Katul | Multiple mechanisms generate a universal scaling with dissipation for the air-water gas transfer velocity[END_REF]. These two regimes and their crossover have been observed experimentally and numerically by [START_REF] Herlina | Isotropic-turbulence-induced mass transfer across a severely contaminated water surface[END_REF][START_REF] Elghobashi | Direct numerical simulation of turbulent flows laden with droplets or bubbles[END_REF], following earlier experimental work by [START_REF] Fortescue | On gas absorption into a turbulent liquid[END_REF] (who used the root mean square of the fluctuating velocity as characteristic velocity).

The same reasoning can be applied in the context of bubbles in a turbulent flow. [START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF] provides a brief discussion on gas dissolution from a bubble of diameter d 0 in a turbulent stream with a characteristic velocity u 0 . The gas transfer rate can then be estimated as Sh ∝ Re 3/4 Sc 1/2 , Sc = ν/D l is the Schmidt number (ratio of kinematic viscosity over mass diffusivity) and Re = d 0 u 0 /ν a turbulence Reynolds number, based on the bubble size d 0 and a turbulence velocity u 0 (the velocity scale considered by [START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF] is the maximum velocity of the eddies in the liquid that flows past the bubble). This regime is equivalent in terms of scalings to the high Reynolds number regime described by [START_REF] Theofanous | Turbulent mass transfer at free, gasliquid interfaces, with applications to open-channel, bubble and jet flows[END_REF]; [START_REF] Magnaudet | Turbulent mass transfer through a flat shear-free surface[END_REF] and its applicability in the context of mass exchange by a bubble swarm is discussed by Colombet et al. (2015).

Numerical methods for interfacial mass transfer started with [START_REF] Sato | Direct simulation of droplet flow with mass transfer at interface[END_REF][START_REF] Davidson | Volume-of-fluid calculation of heat or mass transfer across deforming interfaces in two-fluid flow[END_REF] where the dilute gas concentration is continuous across the interface. [START_REF] Bothe | Direct numerical Farsoiya, Popinet and Deike simulation of mass transfer between rising gas bubbles and water[END_REF] introduced a method to simulate the discontinuous concentration due to solubility. A three-dimensional front tracking model with mass transfer was presented by [START_REF] Darmana | Detailed 3d modeling of mass transfer processes in two-phase flows with dynamic interfaces[END_REF]. A one-fluid formulation for the algebraic volume of fluid method was presented independently by [START_REF] Haroun | Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film[END_REF] and [START_REF] Marschall | Numerical simulation of species transfer across fluid interfaces in free-surface flows using openfoam[END_REF]. [START_REF] Bothe | A volume-of-fluid-based method for mass transfer processes at fluid particles[END_REF] introduced a two-field approach using a geometrical volume of fluid method for multicomponent conjugate mass transfer. Sub-grid scale models to simulate high-Schmidt-number, bubble-mediated mass exchange have been developed by [START_REF] Bothe | A volume-of-fluid-based method for mass transfer processes at fluid particles[END_REF]; [START_REF] Weiner | Advanced subgrid-scale modeling for convectiondominated species transport at fluid interfaces with application to mass transfer from rising bubbles[END_REF]; [START_REF] Claassen | An improved subgrid scale model for front-tracking based simulations of mass transfer from bubbles[END_REF] In the present work, we adapt the classic penetration theory [START_REF] Higbie | The rate of absorption of a pure gas into a still liquid during short periods of exposure[END_REF] to describe the mass transfer from a bubble in a turbulent flow, by considering a turbulent time scale and predict the transfer rate of low-solubility gases, presented in §2. We next present a framework in §3 which combines recent advances in numerical algorithms for interfacial mass transfer, together with progress in turbulent multiphase flow modeling to characterize the diffusive mass transfer of dilute gas from a bubble in turbulence. Dilute gas diffusion from a bubble is similar to conjugate mass transfer from a spherical droplet [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate reynolds number[END_REF]. However, in the present study, the fluid inside the spherical cavity is considerably less dense than the surrounding fluid and hence is called a bubble. We use the Basilisk flow solver [START_REF] Popinet | 2020 Basilisk[END_REF](Popinet & collaborators -2020) ) which uses adaptive mesh refinement, a momentum-conserving scheme for velocity and a geometric volume-of-fluid method to capture the interface. We implement a module for the concentration advection and diffusion using one-fluid formulation without phase change [START_REF] Haroun | Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film[END_REF][START_REF] Bothe | A volume-of-fluid-based method for mass transfer processes at fluid particles[END_REF][START_REF] Taqieddin | Modeling of bubbles hydrodynamics and mass transfer in electrochemical gas-evolving systems[END_REF][START_REF] Yang | Direct numerical simulation of mass transfer and mixing in complex two-phase systems using a coupled volume of fluid and immersed boundary method[END_REF]. The mass transfer module is validated by comparing the numerical results of diffusion from a static bubble with a solution using inverse Laplace transform, as well as for the classic diffusion from a rising bubble in a quiescent liquid. Finally, in §4, we perform DNS of bubble-mediated mass transfer in a homogeneous isotropic turbulent (HIT) flow, solving the three-dimensional, incompressible, two-phase Navier-Stokes equations coupled with an advection-diffusion equation for the gas concentration. We consider bubbles that can deform but at Weber number below the break-up threshold. The theoretical model is confirmed by the DNS results for a wide range of turbulent Péclet numbers.

The present configuration of mass exchange of a dilute component from a bubble to the surrounding turbulent water is especially relevant for bubble mediated CO 2 gas transfer at the ocean-atmosphere interface, as CO 2 is present in a small concentration in the atmosphere so that its exchange with the surrounding water does not change the overall bubble volume, while gases such as N 2 and O 2 which contribute for most of the volume of the bubbles have a much lower solubility and exchange over longer times (see detailed discussions in [START_REF] Woolf | Bubbles and the air-sea exchange of gases in nearsaturation conditions[END_REF]; [START_REF] Keeling | On the role of large bubbles in air-sea gas exchange and supersaturation in the ocean[END_REF]

Theory of bubble mass transfer in turbulence

We consider the mass transfer of dilute gas from a bubble of diameter d 0 , the diffusivity of gas inside and outside the bubble is given by D g and D l respectively. The ratio of momentum to mass diffusivity defines the Schmidt number Sc = ν l /D l . The mass transfer predicted by the two-film theory [START_REF] Whitman | The two-film theory of gas absorption[END_REF]) assumes a steady-state diffusion through the film, while in turbulent flows the fluctuations in the velocity field keep the diffusion transient [START_REF] Treybal | Mass transfer operations[END_REF]. [START_REF] Higbie | The rate of absorption of a pure gas into a still liquid during short periods of exposure[END_REF] penetration theory states that when eddies of liquid in a turbulent flow are exposed to bubbles for a timescale θ, then the mass transfer rate is given by,

k L = 2 √ π D l θ .
(2.1)

In a turbulent flow, the eddies of various sizes will interact with the bubble and advect the gas present in the surrounding thin boundary layers. These interactions will be characterized by a turbulent velocity ũ. Several choices could be considered, in particular the large scale fluctuation velocity, the velocity fluctuations at the bubble scale or the small scale fluctuations. We consider here that the large scale fluctuations of the flow will drive the transfer process, and can be characterized by the magnitude of velocity fluctuations |ũ|= √ 3 u rms . We propose that the time scale of exposure to eddies is given by θ = d 0 /( √ 3 u rms ). The transfer rate, Eq. 2.1, is then,

k L = 2(3) 1/4 √ π D l u rms d 0 (2.2)
For a bubble-mediated mass transfer in turbulent flow, we define the non-dimensional transfer rate, or Sherwood number, Sh = k L d 0 /D l , and the turbulent Péclet number, Pe = u rms d 0 /D l . Finally, Eq. 2.2, can be written,

Sh = 2(3) 1/4 √ π Pe. (2.3)
This equation is analogous to that for the transfer rate for a bubble rising in a quiescent fluid, where the Péclet number would be defined based on the bubble rise velocity [START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF] and mass transfer for bubble swarms where the mean rise velocity of the bubbles in the swarm is considered (Colombet et al. 2015).

Numerical Framework

The Basilisk solver

We solve the three-dimensional, incompressible, two-phase Navier-Stokes equations using the open-source solver Basilisk [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Popinet | A quadtree-adaptive multigrid solver for the serre-green-naghdi equations[END_REF][START_REF] Van Hooft | Towards adaptive grids for atmospheric boundary-layer simulations[END_REF],

∂ t u + ∇ • (uu) = 1 ρ -∇p + ∇ • (µ(∇u + ∇u T )) + γ ρ κδ s n, (3.1) ∇ • u = 0, (3.2) ∂T ∂t + u • ∇T = 0. (3.3)
where u, p, γ, µ, ρ, κ, n and T are the velocity, pressure, surface tension coefficient, viscosity, density, curvature, interface normal and volume fraction fields respectively. The solver has been extensively validated for complex interfacial flows [START_REF] Popinet | A quadtree-adaptive multigrid solver for the serre-green-naghdi equations[END_REF][START_REF] Farsoiya | Axisymmetric viscous interfacial oscillations-theory and simulations[END_REF][START_REF] Gumulya | Dynamics of bubbles rising in pseudo-2d bubble column: Effect of confinement and inertia[END_REF][START_REF] Ruth | Bubble pinch-off in turbulence. PNAS . Sander, Rolf 2015 Compilation of henry's law constants (version 4.0) for water as solvent[END_REF][START_REF] Mostert | Inertial energy dissipation in shallow-water breaking waves[END_REF][START_REF] Berny | Role of all jet drops in mass transfer from bursting bubbles[END_REF]. It uses the projection method to compute the velocity and pressure and the geometric volume of fluid (VoF) method for the evolution of the interface between two immiscible fluids [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]. The Piecewise Linear Interface Calculation (PLIC) geometric interface and flux reconstruction ensures a sharp representation of the interface [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF][START_REF] Marić | Unstructured un-split geometrical volume-of-fluid methods-a review[END_REF] and is combined with an accurate Height-Function curvature calculation and a well-balanced, continuum surface tension model [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF][START_REF] Popinet | Numerical models of surface tension[END_REF].

One fluid formulation for mass transfer of dilute gas

The continuous formulation for mass transfer which we use in this study has been independently developed by [START_REF] Haroun | Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film[END_REF] and [START_REF] Marschall | Numerical simulation of species transfer across fluid interfaces in free-surface flows using openfoam[END_REF]. We have implemented the concentration diffusion of dilute gas [START_REF] Haroun | Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film[END_REF]) using the harmonic mean diffusion coefficient as verified by [START_REF] Deising | A unified single-field model framework for volume-of-fluid simulations of interfacial species transfer applied to bubbly flows[END_REF]. Note that the present study investigates the mass transfer of dilute gas present in the bubble where the effect of loss of volume and phase change are ignored. Numerical methods for phase change require changes in the volume of fluid advection equation and continuity equation [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF][START_REF] Fleckenstein | A volume-of-fluid-based numerical method for multi-component mass transfer with local volume changes[END_REF][START_REF] Maes | A unified single-field volume-of-fluid-based formulation for multi-component interfacial transfer with local volume changes[END_REF][START_REF] Scapin | A volume-of-fluid method for interfaceresolved simulations of phase-changing two-fluid flows[END_REF]. The time evolution of the jth gas concentration, c l/g,j , for the liquid phase l or the gas phase g is given by [START_REF] Standart | The mass, momentum and energy equations for heterogeneous flow systems[END_REF][START_REF] Haroun | Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film[END_REF])

∂c l/g,j ∂t + ∇ • (uc l/g,j ) = -∇ • (J l/g,j ). (3.4)
The continuity of normal fluxes across the interface Σ [START_REF] Standart | The mass, momentum and energy equations for heterogeneous flow systems[END_REF][START_REF] Bothe | A volume-of-fluid-based method for mass transfer processes at fluid particles[END_REF],

(c j (u -u Σ ) + J j ) • n Σ = 0, (3.5)
where u Σ is the interface velocity. As we assume that the transfer of a dilute component does not cause volume change, (3.5) reduces to,

J l,j • n Σ = J g,j • n Σ . (3.6)
The standard assumption of continuous chemical potentials which is good for most applications at interface Σ results in Henry's law, (we refer interested readers to [START_REF] Bothe | A volume-of-fluid-based method for mass transfer processes at fluid particles[END_REF] for the discussion on local chemical equilibrium and the generalized Henry's law),

c l,j = c g,j α j , (3.7)
where the dimensionless ratio of the liquid phase concentration to the gas phase concentration of the component transferred α j , is called Henry's law solubility constant (Sander 2015) (solubility hereafter). The problem investigated is isothermal and there is no bubble breakup or large deformation of the interface which can change the pressure inside the bubble significantly. Hence, solubility which at least depends on temperature and pressure [START_REF] Bothe | A volume-of-fluid-based method for mass transfer processes at fluid particles[END_REF], is assumed to be constant. Introducing variables for the one-fluid formulation which is valid for both phases l and g and for gas j, we get,

c j = T c l + (1 -T )c g , J j = T J l + (1 -T )J g , (3.8)
and the flux is given by

J j = -(T D jl ∇c l + (1 -T )D jg ∇c g ).
(3.9)

The diffusivity for the interfacial cells is calculated using the harmonic mean of the two diffusivities of the gas inside and outside the bubble [START_REF] Haroun | Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film[END_REF],

D j = D jg D jl D jl (1 -T ) + D jg T , (3.10)
to get a single equation for both phases,

∂c j ∂t + ∇ • (uc j ) = ∇ • D j ∇c j -D j c j (α j -1) α j T + (1 -T ) ∇T . (3.11)
The coefficient of c j in the second term containing the right hand side of Eq. (3.11) can be written as β j ,

∂c j ∂t = ∇ • (D j ∇c j + β j c j ) .
(3.12)

Farsoiya, Popinet and Deike

Using a time-implicit Euler discretisation,

c n+1 j -c n j ∆t = ∇ • D j ∇c n+1 j + β j c n+1 j (3.13)
and rearranging the implicit terms, gives

∇ • D j ∇c n+1 j + β j c n+1 j - c n+1 j ∆t = - c n j ∆t , (3.14)
Eq. (3.14) is a set of linear equations which is solved efficiently using the multigrid method [START_REF] Popinet | A quadtree-adaptive multigrid solver for the serre-green-naghdi equations[END_REF].

The solubility boundary condition (3.7) presents a discontinuity for the concentration field across the interface similar to the volume fraction field T . As discussed by [START_REF] Bothe | A volume-of-fluid-based method for mass transfer processes at fluid particles[END_REF], a non-consistent advection leads to artificial diffusion for concentration. For a consistent advection two tracer fields φ g and φ l associated with the volume of fluid T are defined ,

φ l = c l T , φ g = c g (1 -T ) (3.15)
Using Eq. (3.7) and Eq. (3.8), we get,

φ l,j = c n+1 j α j T α j T + (1 -T ) , φ g,j = c n+1 i (1 -T ) α j T + (1 -T ) (3.16)
The advection equation for φ g/l reads,

∂φ g/l ∂t + ∇ • (uφ g/l ) = 0 (3.17)
and is solved using the volume-of-fluid associated fields (López-Herrera et al. 2015) which guarantees strictly non-diffusive transport close to the interface. The concentration tracer is updated after advection using c j = φ g,j + φ l,j .

(3.18)

Validation

We validate the numerical methods and their implementation for static and rising bubble cases.

Diffusion from a static bubble

We propose a new test case of diffusion from a constant-size, static, spherical bubble. Test cases are available for concentration profiles in the case of a planar interface in [START_REF] Bird | Transport phenomena[END_REF]; [START_REF] Haroun | Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film[END_REF]. We provide the solution for the transient concentration both inside and outside the spherical bubble in the form of integrals which are then evaluated numerically.

Consider a static axisymmetric bubble of radius R 0 = d 0 /2, the diffusivities of the gas are D g and D l for the gas and liquid phase respectively. The one-dimensional transient concentration diffusion for a spherical geometry for both inside and outside the bubble is given by

∂c g ∂t = D g 1 r 2 ∂ ∂r r 2 ∂c g ∂r , ∂c l ∂t = D l 1 r 2 ∂ ∂r r 2 ∂c l ∂r .
(3.19)

Applying the Laplace transform to the above equations, we get, c(s) = L[c(t)], given boundary conditions, rc g → 0 as r → 0, cl → 0 as r → ∞, D∂c/∂r = 0 and cl /c g = α at r = R 0 , the concentrations (Laplace transformed) inside the bubble cg and outside cl are given by cg

d 2 cg dr 2 + 2 r dc g dr - s D g cg + 1 D c g0 = 0, d 2 cl dr 2 + 2 r dc l dr - s D l cl + 1 D c l0 = 0, ( 3
(s, r) = c g0 s 1 - 2 ζ(s)r sinh [λ g (s)r] , cl (s, r) = ξc g0 sζ(s)r exp [-λ l (s)r] (3.21)
where,

λ g/l = s D g/l , ζ(s) = ξ exp [-λ l R 0 ] α R 0 + 2 R 0 sinh(λ g R 0 ) ξ(s) = 2D g (λ g R 0 cosh(λ g R 0 ) -sinh(λ g R 0 )) D l exp[-λ l R 0 ](1 + λ l R 0 ) , c g0 = c g (0, r), c l0 = 0. (3.22)
Eq. 3.21, can be inverted using Cauchy residue theorem (using keyhole contour as discussed in Farsoiya et al. (2020)),

c g (t, r) = - 2c g0 πr ∞ 0 1 x Im sinh(λ g (-x)r) ζ(-x) exp[-xt] dx, (3.23) c l (t, r) = c g0 πr ∞ 0 1 x Im ξ(-x) ζ(-x) exp [-λ l (-x)r] exp[-xt] dx, (3.24)
where Im(•) is the imaginary part of a complex number. We validate the results of the

Farsoiya, Popinet and Deike

Basilisk solver against (3.23) and (3.24) (using the numerical integration functions based on double exponential quadrature in Wolfram Research, Inc. ( 2020)) in Fig. 1.

A static bubble of diameter d 0 /L = 0.2, diffusivity ratio D g /D l = 10 and solubility α = 10 -3 to 10 -1 is considered for this test. Figure 1(a) and 1(b) show the transient concentrations inside the bubble at r/d 0 = 0.25 and outside the bubble at r/d 0 = 0.75 respectively. The accuracy of the numerical solution is very good above 25 cells per diameter. The stringent solubility condition α = 10 -3 is achieved at the interface r/d 0 = 0.5 even at a low resolution of six cells per diameter (figure 1(c)). To quantify the convergence, we compute the norm, max|c 11 -c n |/c 11 at different resolutions where c 11 and c n are numerical solutions at uniform resolutions 2 11 × 2 11 and lower respectively. Figure 1(d) shows first-order of convergence with respect to the grid size ∆x. The effect of solubility on the error has not been discussed in earlier studies. We show in Figure 1(d) that decreasing the solubility at the interface increases the error in the solution for a given resolution, but that the error remain small at high resolution. [START_REF] Clift | Pascal 2015 Dynamics and mass transfer of rising bubbles in a homogenous swarm at large gas volume fraction[END_REF][START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF]) and numerical results [START_REF] Roghair | Direct numerical simulations of hydrodynamics and mass transfer in dense bubbly flows[END_REF][START_REF] Deising | A unified single-field model framework for volume-of-fluid simulations of interfacial species transfer applied to bubbly flows[END_REF].

Diffusion from a rising bubble

For the validation of the advection-diffusion scheme, we consider a bubble in a quiescent fluid rising due to buoyancy. We set up the test cases with the parameters considered in earlier studies [START_REF] Darmana | Detailed 3d modeling of mass transfer processes in two-phase flows with dynamic interfaces[END_REF][START_REF] Roghair | Direct numerical simulations of hydrodynamics and mass transfer in dense bubbly flows[END_REF][START_REF] Deising | A unified single-field model framework for volume-of-fluid simulations of interfacial species transfer applied to bubbly flows[END_REF][START_REF] Jia | Investigation of a free rising bubble with mass transfer by an arbitrary lagrangian-eulerian method[END_REF], which encompass significant variations in bubble conditions, through variations of the Bond and Archimedes (or Morton) numbers. The rise velocity of a bubble in a quiescent liquid is indeed determined by the bubble-liquid physical parameters, summarized by the Archimedes number Ar = gd 3 0 ρ l (ρ l -ρ g )/µ 2 or Morton number Mo = gµ 4 l /(ρ l γ 3 ) and Bond number Bo = ρ l gd 2 0 γ (Moore 1965; [START_REF] Maxworthy | Experiments on the rise of air bubbles in clean viscous liquids[END_REF][START_REF] Clift | Pascal 2015 Dynamics and mass transfer of rising bubbles in a homogenous swarm at large gas volume fraction[END_REF][START_REF] Cano-Lozano | Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability[END_REF]). The terminal rise velocity U can be computed theoretically [START_REF] Moore | The velocity of rise of distorted gas bubbles in a liquid of small viscosity[END_REF][START_REF] Clift | Pascal 2015 Dynamics and mass transfer of rising bubbles in a homogenous swarm at large gas volume fraction[END_REF], and expressed as a non-dimensional bubble Reynolds number Re = ρ l U d 0 /µ l . We consider four bubble configurations leading to Re from 5 to 100, for Bond numbers ranging from 1 to 40, and Archimedes numbers ranging from 100 to 8000. The mass transfer for these bubble conditions is then computed, considering the ratio of momentum to mass diffusivity given by the Schmidt number Sc = ν l /D l = 1 and the gas solubility α = 1/30, as in [START_REF] Deising | A unified single-field model framework for volume-of-fluid simulations of interfacial species transfer applied to bubbly flows[END_REF]. The fourth configuration is further tested for high Schmidt numbers of 10 and 100 in D 2 and D 3 respectively. Our results are shown in figure 2 and validated against previous work as shown in Table 1. Note that to achieve the correct mass transfer, the correct rise velocity must be obtained, since [START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF]. Very good agreement between the theoretical and numerical mass transfer is observed except for case C where the terminal shape is far from spherical. The scripts sufficient to reproduce these results are provided in Farsoiya et al. (2020b).

Farsoiya, Popinet and Deike the mass transfer will directly depend on the bubble velocity. Figure 2 (a to d) shows examples of the shape of the bubble and the surrounding adaptive grid once the terminal velocity is reached. The average gas concentrations inside and outside the bubble are computed as

cg,j = 1 V g Vg c j dV, cl,j = 1 V l V l c j dV (3.25)
where V g and V l are the volumes of bubble and liquid respectively. The mass transfer rate k L is calculated as,

k L = cn+1 g,j -cn g,j A g ∆t(α cn+1/2 g,j - cn+1/2 l,j ) (3.26)
where A g is the instantaneous surface area of the bubble. The non-dimensionalized transfer rates are computed for both the axisymmetric and three-dimensional simulations with adaptive mesh refinement and a maximum refinement corresponding to a resolution on the bubble of d 0 /∆x = 100. [START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF] derived the mass transfer rates from a rising bubble with constant size and surface concentration, which reads, Sh = 2 √ π

√

Pe where Pe = Re Sc, is the bubble Péclet number. Upon reaching terminal velocity the transfer rates are steady after t U/d 0 ≥ 2 as shown in figure 2(e) and are accurately predicted by [START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF] as shown in figure 2(e) and 2(f). The values of the non-dimensional transfer rate, scaled by the diffusion velocity scale, or Sherwood number Sh = kd 0 /D l are also provided in Table 1. Overall, the agreement between theoretical and numerical transfer rates is very good. The deviation from the spherical shape in case C, visible in Figure 2(c), causes a difference with the predicted transfer rate which assumes a spherical bubble shape.

Mass transfer in homogeneous and isotropic turbulence

We now present DNS of diffusion of dilute gas from a bubble inside a surrounding homogeneous and isotropic turbulent flow (HIT). We follow recent work on bubble deformation in turbulence [START_REF] Perrard | Bubble deformation by a turbulent flow[END_REF][START_REF] Rivière | Sub-hinze scale bubble production in turbulent bubble break-up[END_REF] first prepare the turbulent flow by solving the momentum equation with a forcing term, and then insert a bubble at the center of the HIT flow, once a statistically stationary state has been reached. A similar approach has been used to study the interaction of droplets with isotropic turbulence [START_REF] Dodd | On the interaction of taylor length scale size droplets and isotropic turbulence[END_REF]Elghobashi 2019a). As the bubble is inserted, the volume of fluid advection Eq. ( 3.3) and mass advection-diffusion Eq. (3.11) is solved coupled with the momentum equations Eq. (3.1). We then calculate the mass transfer rate from the bubble and compare it with the theory presented in §2.

Case

Resolution

Precursor simulation for isotropic turbulence

The turbulent flow is generated by adding a linear volumetric forcing term f = Au(x, t) in Eq. 3.1. This approach has been introduced by [START_REF] Rosales | Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties[END_REF] and yields similar turbulence properties as using a forcing in spectral space and leads to a wellcharacterized homogeneous and isotropic turbulent flow. Such an approach has been applied to study bubble rising in turbulence by [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF], and more recently we used this approach within the Basilisk solver to study bubble deformation in turbulence [START_REF] Perrard | Bubble deformation by a turbulent flow[END_REF][START_REF] Rivière | Sub-hinze scale bubble production in turbulent bubble break-up[END_REF]. We consider a 3D periodic box of size L, for a precursor simulation to achieve isotropic turbulence. We use adaptive mesh refinement on the velocity field, and the maximum level of refinement can be used to compare the resolution with that of a fixed grid. The turbulent flow is generated for increasing resolutions with the maximum level of refinement going from 6 to 8, corresponding to an equivalent of 64 3 to 256 3 grid points. The resolution will be increased once we insert the bubble. The turbulence state is characterized by the kinetic energy density K, turbulence dissipation rate and Taylor microscale Reynolds number Re λ , which are given by [START_REF] Pope | Turbulent flows[END_REF])

K = 1 V V 1 2 ρ l |u (x, t)| 2 dV, = 1 V V ν l ∂u i ∂x j ∂u i ∂x j dV, Re λ = 2K 3ν l 15ν l (4.1)
and are computed over time to characterize the turbulent flow. The root mean square of the velocity is u rms = 2K/3ρ l , and the eddy turn over time at the scale of the bubble of diameter d 0 is given by [START_REF] Pope | Turbulent flows[END_REF][START_REF] Perrard | Bubble deformation by a turbulent flow[END_REF]. Figure 3(a) shows the evolution of the turbulent kinetic energy with time for increasing Reynolds numbers. It shows that at t/t c ≈ 25 the flow has reached a statistically stationary state. We show that the state is grid-independent when using an adaptive mesh refinement of L7≡ 2 7 and L8≡ 2 8 for the three cases. The turbulence Reynolds number and the turbulence dissipation rate have similar time evolutions. The range of turbulent Taylor Reynolds numbers is Re λ ≈ 38 to 77 which is a typical value for current two-phase simulations of turbulent flow [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF]Elghobashi 2019b). We characterize the turbulent stationary state using the second-order structure functions Farsoiya, Popinet and Deike in the longitudinal D LL (r) and the transverse direction D N N (r), given by D LL (r) =

t c = d 2/3 0 -1/3 (
1 3 i (u i (r, t) -u i (r + dr i , t)) 2 , and D N N (r) = 1 6 i =j (u i (r, t) -u i (r + dr j , t)) 2 ,
where ri is the unit vector along the i th direction. Figure 3(b) shows that the scaled structure functions plateau at C = 2 [START_REF] Pope | Turbulent flows[END_REF] in the inertial range. The relation D LL = 3/4D N N is verified and the inertial range is relatively limited due to the relatively coarse resolution and limited turbulence Reynolds number. The bubble is inserted once the turbulent stationary state is reached, and is of a size within the inertial range, the turbulence at this scale being reasonable, as described in [START_REF] Perrard | Bubble deformation by a turbulent flow[END_REF][START_REF] Rivière | Sub-hinze scale bubble production in turbulent bubble break-up[END_REF].

Bubble insertion

The bubble is inserted at the center of the box after reaching isotropic turbulence, i.e. for t 0 > 25t c (See Fig. 3(a)). The bubble is of diameter d 0 , viscosity µ b , and density ρ b surrounded by a liquid of viscosity µ l and density ρ l . The solubility of dilute gas is α j = 0.3 which is transferred across the interface. The Weber number, We ≡ ρ l u 2 rms d 0 /γ = 1.3 is below the critical number for bubble breakup [START_REF] Perrard | Bubble deformation by a turbulent flow[END_REF][START_REF] Rivière | Sub-hinze scale bubble production in turbulent bubble break-up[END_REF], so that all results are for bubbles that can deform but do not break. Theoretical discussion and experimental data in [START_REF] Theofanous | Turbulent mass transfer at free, gasliquid interfaces, with applications to open-channel, bubble and jet flows[END_REF] provide the approximation that the transfer rate is a weak function of surface tension within 3 orders of magnitude (0.001 W e 1) and depends mainly on the bulk turbulence properties. Diffusion rates for six different gases are calculated corresponding to Schmidt numbers ranging from 1 to 100. The bubble size with respect to the Taylor micro-scale length (λ = 15νu 2 rms / ) is in the range 1.72-2.82, and the box size is L = 7.5d 0 . The turbulence properties and simulation parameters are given in Table 2. Th higher resolutions are used to properly resolve the diffusion of mass and the bubble deformation dynamics, with effective resolution using an adaptive mesh refinement of L10≡ 2 10 , L11≡ 2 11 and L12 ≡ 2 12 . As will be discussed in detail later, higher Schmidt numbers lead to thinner diffusive boundary layers which require smaller grid sizes. The resolution of smallest momentum and mass length scales are discussed in appendix A.

Mass transfer from the bubble to the surrounding turbulent flow

As the bubble moves with the flow, vortices of outer liquid come in contact with the interface, as seen in Figures 4(a)-4(b). Molecular diffusion of gas around the bubble interacts with the eddies of the flow and unsteady boundary layers are formed. The corresponding concentration field for Sc = 1 and Sc = 10 are shown in figures 4(c)-4(d) and 4(e)-4(f). It can be observed that the concentration fields have followed the flow field due to advection. The gas with high diffusivity (Sc = 1) has thicker boundary layers (Figures 4(c) and(d)) compared to the low diffusivity gas (Sc = 10, Figures 4(e) and(f)) subjected to the exact same flow. Figure 5 shows a 3D rendering of the same times of the simulations, displaying the 3D concentration field around the bubble, as the bubble moves in the flow, with the lower Sc = 1 in figures 5(a)-5(b), and higher Sc = 10 in figures 5(c)-5(d). Again, higher Schmidt numbers lead to thinner boundary layers around the bubble for the same flow field.

We next compute the total transfer rates k L of gases corresponding to increasing Schmidt numbers, and show their evolution with time in Fig. 6 Finally, Figure 7 shows the Sherwood numbers obtained from the bubble mass transfer in homogeneous and isotropic turbulence simulations for a wide range of turbulent Péclet numbers (50 < Pe ≤ 2 × 10 4 ) for three different resolutions (simulations conditions are summarized in table 2). Very good agreement is observed between Eq. (2.3), up to Pe ∼ 10 4 . As such, Figure 7 validates Eq. (2.3), which predicts the transfer rates from the bubble in homogeneous and isotropic turbulent flow. We note that for Pe ≈ 2 × 10 4 , an effective resolution L12 brings the simulated transfer rate closer to our theoretical prediction, thanks to increased resolution in the turbulent-diffusive boundary layer next to the bubble interface. We recall that Figure 7 shows the scaling for the transfer rate k L , as the non-dimensional Sh number, which reads Sh ∝ Re 0.5 Sc 0.5 = Pe, and considering the definition of the Sh number, this corresponds to k L ∝ Sc -0.5 urmsν l d0 . We note that [START_REF] Dodd | Analysis of droplet evaporation in isotropic turbulence through droplet-resolved dns[END_REF] presents simulations of a droplet evaporating in turbulence and their evaporation rates at early times before significant volume change has occurred are in agreement with our model, as shown in figure 7.

The computed Sherwood numbers show convergence to the predicted value as the grid resolution increases. For the thinner boundary layers, corresponding to Re λ = 77 and high Schmidt number at d 0 /∆x = 273 resolution, the total number of cells in the domain reached 27.5 million and used 2.16 × 10 4 CPU hours for 0.5 eddy turnover time when the transfer rates approach a steady state. A corresponding constant-resolution direct numerical simulation would have required (2 11 ) 3 ≈ 8.6 billion cells, which illustrates the gain in efficiency brought by adaptive mesh refinement. It highlights the potential of adaptive mesh refinement methods for complex two-phase flows at high Peclet number, complementary to sub-grid scale modeling approaches.

We note that in §2, we have proposed the large-scale velocity fluctuations as the controlling velocity scale. As already mentioned, another choice could have been the velocity fluctuations at the scale of the bubble. Such a choice would not have changed the results obtained in figure 7, as these velocities are similar at the Reynolds number at the Taylor-scale we consider.

Finally let us comment that at high Reynolds numbers, in the case of a flat interface, Eq. ( 2.3) Sh = 2(3) 1/4 Pe/π. For higher Péclet numbers (1.02 × 10 4 , 2.05 × 10 4 ), higher resolutions may converge to the prediction as visible when using L12 simulations. The transfer rates from [START_REF] Dodd | Analysis of droplet evaporation in isotropic turbulence through droplet-resolved dns[END_REF] corresponds to the transfer rates from simulations of an evaporating drop at early times and agree with our data and model.

the scaling of the gas transfer is shifted from Re 1/2 to Re 3/4 due to the role of small eddies in enhancing the gas exchange [START_REF] Theofanous | Turbulent mass transfer at free, gasliquid interfaces, with applications to open-channel, bubble and jet flows[END_REF][START_REF] Herlina | Isotropic-turbulence-induced mass transfer across a severely contaminated water surface[END_REF]. Whether such a regime could be observed in the case of bubble mass exchange at high Peclet or Reynolds number remains to be investigated.

Concluding remarks

We propose a simple theoretical formula for mass transfer rate of dilute gas from bubbles for a dilute component within a turbulent flow, extending the classic formulation from [START_REF] Higbie | The rate of absorption of a pure gas into a still liquid during short periods of exposure[END_REF] to a homogeneous and isotropic turbulent flow, by considering that the renewal time-scale is controlled by the turbulence fluctuations at large scale and the bubble size:

Sh = 2(3) 1/4 √ π Pe.
We developed and presented a versatile numerical framework to perform direct numerical simulations of multiphase flows, and tested the theoretical prediction against numerical simulations. The agreement between theory and numerical results is very good, as long as the diffusive boundary layer thickness is correctly resolved. The numerical framework can now be used for more complex physical configurations such as bubbles under a breaking wave and bubble swarms. The proposed theoretical scaling for dilute mass transfer is remarkable as it could be leveraged for multi-scale systems, such as bubbles entrained by breaking waves in the upper ocean, or bubbles evolving in turbulent swarms, and avoid solving for the very high Schmidt number necessary in practical problems (for example the Schmidt number of CO 2 at 20 degrees is ∼660). Finally, our formulation could also be used to improve recent bubble-mediated mass transfer models which have, up to now, used the bubble transfer rate in a quiescent flow [START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF][START_REF] Deike | Gas transfer by breaking waves[END_REF].

As part of this work, a video has been published in the gallery of fluid motion (Farsoiya et al. 2020c).

  . Recent advances in numerical methods by Tanguy et al. (2014); Fleckenstein & Bothe (2015); Maes & Bubble-mediated transfer 3 Soulaine (2020); Scapin et al. (2020) allow to simulate problems of mass transfer with local volume changes.

Figure 1 :

 1 Figure 1: Diffusion from a static bubble, comparing the numerical results with (3.23) and (3.24). (a) Concentration inside the bubble at r/d 0 = 0.25. (b) Concentration outside the bubble at r/d 0 = 0.75. (c) Radial profile at time tD l /d 2 0 = 1.2. (d) Maximum relative error at different resolutions max|c 11 -c n |/c 11 , where c 11 and c n are numerical solutions at resolution 2 11 and lower respectively, and for different solubilities α, displaying first order convergence. The scripts sufficient to reproduce these results are provided in Farsoiya et al. (2020a).

Figure 2 :

 2 Figure 2: Mass transfer of dilute gas from a bubble rising in a quiescent flow. (a)-(d)Concentration c and bubble interface (three dimensional) at t U/d 0 ≈ 4 for the four cases considered (increasing bubble Re number). (e) Evolution with time of transfer rates for axisymmetric and three-dimensional simulations and from[START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF] using the computed terminal velocity, (f) Steady state transfer rate (for t U/d 0 ≥ 2) as a function of Péclet number compared against[START_REF] Levich | Modeling bubbles and dissolved gases in the ocean[END_REF]. Very good agreement between the theoretical and numerical mass transfer is observed except for case C where the terminal shape is far from spherical. The scripts sufficient to reproduce these results are provided inFarsoiya et al. (2020b).

Figure 3 :

 3 Figure 3: Properties of the homogeneous and isotropic turbulent flow. (a) Turbulent kinetic energy as a function of time. After a short transient, a statistically stationary state is reached. The bubble is inserted once the statistically stationary state is reached. (b) Second-order structure function, D LL and D N N in the longitudinal and transverse directions respectively, compensated by the homogeneous and isotropic turbulence scaling (r ) -2/3 and D LL = 3/4D N N . Turbulence theory 4/3D LL (r)(r ) 2/3 is superimposed in red dashed line. η is the Kolmogorov length scale. The bubble has a size comparable to the Taylor turbulence scale, within the inertial range.

Figure 4 :FarsoiyaFigure 5 :

 45 Figure 4: Mass diffusion from the bubble in turbulence at two different times, (t-t 0 )/t c = 1.5 (top row) and 2 (bottom row), showing a 2D planar cut of the magnitude of the vorticity field (a)-(b) and concentration field for Re λ = 77, Sc = 1 (P e (t) = 205) (c)-(d) and Re λ = 77, Sc = 10 (P e (t) = 2050) (e)-(f). The wake in the vorticity field presents similarities with the structure of the concentration field. Higher Schmidt numbers lead to a thinner boundary layer around the bubble and a thinner wake structure.

Figure 6 :

 6 Figure 6: Non-dimensional mass transfer rates Sh as a function of time, as the bubble is exposed to the turbulent flow. (a,b) Lower resolution d 0 /∆x = 136 (level 10), for Re λ = 38 (a) and Re λ = 77 (b). (c,d) Higher resolution (solid line for level 11 d 0 /∆x = 273, dotted line for level 12 d 0 /∆x = 546), for Re λ = 38 (c) and Re λ = 77 (d). Dashed lines represent the theoretical prediction, Eq. 2.3 for the different Schmidt numbers. At steady state, very good agreement between simulations and theory is achieved for Pe = u rms d 0 /D l = Sc(u rms d 0 /D l ) ≤ 10 4 , which corresponds to a diffusive boundary layer δ k resolved with more than ∼ 4 grid points.

  .20)

	(a)	1.000					Eq. (3.23)	(b)
		0.995					d0/Δx ≈ 6 d0/Δx ≈ 25		0.0006
		0.990					d0/Δx ≈ 102
	cb/cb0	0.985							cl/cb0	0.0004
		0.980								0.0002					Eq. (3.24) d0/Δx ≈ 6
		0.975								0.0000					d0/Δx ≈ 2Δ d0/Δx ≈ 102
		0.0	0.2	0.4	0.6	0.8	1.0	1.2		0.0	0.2	0.4	0.6	0.8	1.0	1.2
					t l/d 2 0								t l/d 2 0
		1.0				Eq. (3.23), (3.24)
		0.8				d0/Δx ≈ 6 d0/Δx ≈ 2Δ
		0.6				d0/Δx ≈ 102
	c/c0	0.4
		0.2
		0.0
		0.0	0.2	0.4	0.6		0.8	1.0
					r/d0

Table 1 :

 1 Reynolds and Sherwood numbers (for t U/d 0 ≥ 2) for axisymmetric and 3D rising bubbles obtained from our numerical study, indicated as present work (PW Axi, and PW 3D), compared with existing theoretical

	Mo	Bo			Re						Sh
	Ar		Clift	Roghair	Deising	PW:	PW:	Levich	Roghair	Deising	PW: Axi		PW: 3D
	≈ (100 -		et al.	(2012)	et al.	Axi	3D	(1962)	(2012)	et al.
	8000)		(2005)		(2016)					(2016)
	A 10 -4	1.0	5	5.53	6.4	5.6	5.6	2.5	3.2	3.6	3.07-3.19		3.17-3.25
	B 5 • 10 -4	3.125 10	11.31	10.5	10.9 10.9	3.7	4.4	4.2	3.99-4.09		4.07-4.09
	C 9.2•10 -3 40	33	32.35	32.5	32.8 32.9	6.4	6.3	6.4	4.63-5.42		5.25-6
	D1 5 • 10 -7	3.125 103 102.78 97.7	102.4 102.9 11.5 12.03	11.5	10.52-10.98 11.15-11.33
	D2 5 • 10 -7	3.125 103 102.78 97.7	102.4 102.9 36.5 -	-	33.12-37.10 35.07-36.38
	D3 5 • 10 -7	3.125 103 102.78 97.7	102.4 102.9 114.7 -	-	108.10	-	116.91-
											122.13		121.18

Table 2 :

 2 Simulation parameters (with adaptive mesh refinement) of the turbulence simulation of mass transfer. Three Reynolds numbers are used, with two effective resolutions, and a range of Schmidt numbers. The Weber number, density and viscosity ratio are kept constant.

		(2 L )	d 0 /λ We ρ r µ r	Sc	Re λ
	1	L10 ≡ 2 10 , L11 ≡ 2 11	1.72 1.3 850 25 (1, 2, 10, 20, 50, 100) 38
	2	L10 ≡ 2 10 , L11 ≡ 2 11	2.30 1.3 850 25 (1, 2, 10, 20, 50, 100) 55
	3 L10 ≡ 2 10 , L11 ≡ 2 11 , L12 ≡ 2 12 2.82 1.3 850 25 (1, 2, 10, 20, 50, 100) 77

  Re λ ≈ 38 d 0 /Δx = 136 Re λ ≈ 55 d 0 /Δx = 136 Re λ ≈ 77 d 0 /Δx = 136 Re λ ≈ 38 d 0 /Δx = 273 Re λ ≈ 55 d 0 /Δx = 273 Re λ ≈ 77 d 0 /Δx = 273 Re λ ≈ 77 d 0 /Δx = 546 Dodd et. al. Δ2021) Figure 7: The transfer rates (for (t-t 0 )/t c ≥ 0.2) with respect to Péclet number. Symbols are for different Re λ and resolution. For Pe ≤ 10 4 , the numerical results are grid converged (first order) between L10 and L11, and agree very well with the theoretical prediction,

						Eq. 2.3
	Sh	10 2
		10 1
		10 2	10 3	Pe	10 4	10 5

† Email address for correspondence: ldeike@princeton.edu

Acknowledgments

This work was supported by the NSF CAREER award 1844932, the Catalysis Initiative at Princeton and the Cooperative Institute for Modeling the Earth System (CIMES) between Princeton and NOAA-GFDL. We would like to acknowledge highperformance computing support from Tiger, Princeton research computing and Cheyenne (doi:10.5065/D6RX99HX) provided by NCARs Computational and Information Systems Laboratory, sponsored by the National Science Foundation. Declaration of Interests. The authors report no conflict of interest.

Appendix A

The accepted resolution criteria for DNS in the literature [START_REF] Overholt | Direct numerical simulation of a passive scalar with imposed mean gradient in isotropic turbulence[END_REF][START_REF] Pope | Turbulent flows[END_REF][START_REF] Schumacher | Very fine structures in scalar mixing[END_REF][START_REF] Dodd | Analysis of droplet evaporation in isotropic turbulence through droplet-resolved dns[END_REF]) is typically k max η > 1.5 and k max η B > 1.5, where k max , η and η B are the maximum resolved wavenumber k max = πN/L, the Kolmogorov and Batchelor scale respectively. The Kolmogorov length scale η = (ν 3 l / ) 1/4 defines the length scale at which viscous dissipation becomes dominant while the Batchelor scale is defined as η B = η/ √ Sc. Figure 8 shows these characteristic length scales as functions of the turbulent Peclet number for all our simulations. In all of the cases the Kolmogorov length scale is well resolved, with k max η > 8, in agreement with the fact that convergence is already achieved in lower resolutions as shown in figure 3a. For the highest Péclet number, the Batchelor length scale is resolved up to k max η B ≈ 3.2 (refinement d 0 /∆x = 546). The boundary layer thickness δ ν > 2.5η and δ k > 2.5η B .