Changgui Zhang

ON THETA-TYPE FUNCTIONS IN THE FORM (x; q) ∞

, a function is said to be of thetatype when its asymptotic behavior near any root of unity is similar to what happened for Jacobi theta functions. It is shown that only four Euler infinite products have this property. That this is the case is obtained by investigating the analyticity obstacle of a Laplace-type integral of the exponential generating function of Bernoulli numbers.

Introduction

In his last letter to Hardy, Ramanujan wrote that he had discovered very interesting functions that he called mock ϑ-functions. As was said in Watson's L.M.S. presidential address [START_REF] Watson | The final problem: An account of the mock theta functions[END_REF], the first three pages, where Ramanujan explained what he meant by "mock ϑ-functions", are very obscure. Therefore, Watson quoted the following comment of Hardy:

A mock ϑ-function is a function defined by a q-series convergent when |q| < 1, for which we can calculate asymptotic formulae, when q tends to a "rational point" e 2rπi/s of the unit circle, of the same degree of precision as those furnished for the ordinary ϑ-functions by the theory of linear transformation.

In our previous work [START_REF] Zhang | On the mock theta behavior of Appell-Lerch Series[END_REF], we proposed definitions of what we call theta-type, false theta-type and mock theta-type functions, following directly from the abovementioned comment of Hardy. The main goal of this paper is to determine the possible values of x for which the Euler q-exponential function (x; q) ∞ is of thetatype, where (1.1) (x; q) ∞ = n≥0

(1 -x q n ) .

1.1. Statement of Main Theorem. Let ζ = e 2πir = e(r) be any root of unity, with r ∈ Q. As in [START_REF] Zhang | On the mock theta behavior of Appell-Lerch Series[END_REF], a function f (q) is said to be of theta-type as q a.r.

-→ ζ, and one writes f ∈ T ζ , if there exists a quadruplet (υ, λ, I, γ), composed of a couple (υ, λ) ∈ Q × R, a strictly increasing sequence I ⊂ R and a C * -valued map γ defined on I, such that the following relation holds for any N ∈ Z ≥0 as τ a.v.

-→ r:

(1.2) f (q) = i τ υ e(λτ )   k∈I∩(-∞,N ] γ(k) q k 1 + o(q N 1 )
  .

Date: September 24, 2021. The author was supported by Labex CEMPI (Centre Européen pour les Mathématiques, la Physique et leurs Interaction).

Here and in what follows, q = e 2πiτ = e(τ ), τ > 0, τ = τ -r and q 1 = e(-1 τ ). The symbol "q a.r.

-→ ζ" could be read as "q almost radially converges to the root of unity ζ". At the same time, "τ a.v.

-→ r" means that "τ almost vertically converges to the rational point r"; see §1.2 (ii), below.

The above-named variable q 1 may be considered as the modular variable with respect to the root ζ. By considering the respective modular formulae, one can easily see that the ordinary ϑ-functions, as well as the famous Dedekind η-function, satisfy the condition in (1.2) for any root of unity ζ.

Let U denote the set of the roots of unity. One remembers that η(τ ) = q 1/24 (q; q) ∞ , where q = e(τ ). Thus, one can observe that (q; q) ∞ ∈ T ζ for any ζ ∈ U. Furthermore, an elementary calculation shows that the following identities hold: (-q; q) ∞ = (q 2 ; q 2 ) ∞ (q; q) ∞ , ( √ q; q) ∞ = ( √ q; √ q) ∞ (q; q) ∞ and (-√ q; q) ∞ = (q; q) 2 ∞ ( √ q; √ q) ∞ (q 2 ; q 2 ) ∞ .

As one may check directly from (1.2), the set ∩ ζ∈U (T ζ \ {0}) constitutes a multiplicative group which is stable by the ramification operator q → q ν for all ν ∈ Q >0 . Therefore, one obtains from the above the following property:

Remark 1.1. Given x ∈ {q, -q, √ q, -√ q}, one has (x; q) ∞ ∈ T ζ for all ζ ∈ U.

The following result can be viewed as being converse to the above statement:

Theorem 1.1 (Main Theorem). Let (x 0 , β) ∈ C × R be such that |x 0 | = 1 and β = 0, and consider x = x 0 q β . Then, the following conditions are equivalent:

(1) one has (x; q) ∞ ∈ T 1 , with ζ = 1 = e(0);

(2) there exists a root of unity ζ = e 2πir = e(r) such that (x; q) ∈ T ζ ;

(3) one has x ∈ {q, -q, √ q, -√ q}.

It should be noted that the four Euler q-exponential functions obtained in Theorem 1.1 are intimately linked to the quadruplet (φ(q), ψ(q), f (-q), χ(q)), situated at the heart of the Ramanujan's theory about theta-functions and modular equations [START_REF] Berndt | Ramanujan's notebooks. Part III[END_REF]Ch. 16]. Indeed, one can observe that f (-q) = (q; q) ∞ , χ(q) = (-q; q 2 ) ∞ and, furthermore, that φ(q) = (-q; q 2 ) ∞ (q 2 ; q 2 ) ∞ (q; q 2 ) ∞ (-q 2 ; q 2 ) ∞ , ψ(q) = (q 2 ; q 2 ) ∞ (q; q 2 ) ∞ .

In particular, these relations imply that (q; q 2 ) ∞ = f (-q 2 ) ψ(q) , (-q 2 ; q 2 ) ∞ = χ(q) ψ(q) φ(q) .

1.2. Some ideas for the proof of Main Theorem. The following notational conventions will be used through the whole paper:

(i) For simplicity, we will write e(z) = e 2πiz for all z ∈ C, and this gives, in particular, a map τ → q = e(τ ) from the Poincaré half-plane H onto the unit disc D. Furthermore, one has e(Q ∩ [0, 1)) = U, where U denotes the set of the roots of unity.

(ii) Let ζ = e(r) ∈ U and r ∈ Q ∩ [0, 1). One writes q a.r.

-→ ζ if there exists ∈ (0, π 2 ) such that q → ζ in the sector | arg(q -ζ) + r| < inside the unit disc D. Similarly, one writes τ a.v.

-→ r if there exists ∈ (0, π 2 ) such that τ → r in the sector | arg(τ -r) -π 2 | < in the upper half-plane H. (iii) Given ζ = e(r) ∈ U, one says that f is exponentially small as q a.r.

-→ ζ or τ a.v.

-→ r, and one writes f

∈ A ≤-1 ζ if there exists (C, κ) ∈ R 2 >0 such that |f (q)| ≤ C e -κ/|τ -r| for all q = e(τ ) in some sector {|τ -r| < ρ, | arg(τ - r) -π 2 | < }, where ∈ (0, π 2 ) and ρ > 0. (iv) For any given (z, τ ) ∈ C×H, we set (z | τ ) 0 = (x; q) 0 = 1, (z | τ ) ∞ = (x; q) ∞ and, for N ∈ Z >0 , (1.3) (z | τ ) N = (x; q) N = N -1 n=0 (1 -x q n ) ,
where x = e(z) and q = e(τ ). This is in line with (x; q) ∞ given in (1.1). Letting k = min I and c = γ(k) in (1.2) implies that

f (q) = c i τ υ e λ τ - k τ (1 + f 1 (q)) , f 1 ∈ A ≤-1 ζ .
If one takes the principal branch of the logarithm for both members of the above equation, one can observe the following fact:

Remark 1.2 (Asymptotic form of the theta-type functions). Given ζ ∈ U and f (q) ∈ T ζ , there exists a quadruplet (υ, c ∞ , c 0 , c 1 ) ∈ Q × (iR) × C × (iR) such that (1.4) log f (q) = υ log i τ + c ∞ τ + c 0 + c 1 τ mod A ≤-1 ζ .
The formula stated in (1.4) can be viewed as a necessary condition for any function to be of theta-type. Furthermore, let C{z} be the set of analytic functions at z = 0. One remembers that A ≤-1 ζ ∩ C{τ } = {0}. By replacing c 0 + c 1 τ with any convergent power series of τ in the relation in (1.4), we will introduce a larger class of functions as follows: [START_REF] Andrews | Special functions, Encyclopedia of Mathematics and its Applications[END_REF]. One says that f (q) admits an exponential-convergent expansion as q a.r.

Definition 1.1. Let ζ = e(r) ∈ U and r ∈ Q ∩ [0,
-→ ζ or τ a.v.

-→ r and one writes f

∈ C ζ if there exists (υ, c ∞ ) ∈ Q × (iR) such that (1.5) log f (q) = υ log i τ + c ∞ τ mod C{τ } ⊕ A ≤-1 ζ . It is obvious that T ζ ⊂ C ζ .
In this way, one will see that Theorem 1.1 can be easily deduced from Theorem 1.2. Let (x 0 , β) ∈ C × R be such that |x 0 | = 1 and β = 0, and consider x = x 0 q β . Then, the following conditions are equivalent:

(1) one has (x; q) ∈ C 1 , with ζ = 1 = e 2πi0 ;

(2) there exists a root of unity ζ = e 2πir such that (x; q) ∈ C ζ ;

(3) one has x 0 ∈ {1, -1} and β ∈ 1 2 Z \ {0}. A modular-like formula has been found for (x; q) ∞ in [START_REF] Zhang | On the modular behaviour of the infinite product (1-x)(1-xq)(1-xq 2 )(1-xq 3 )[END_REF]Th. 3.2] and [START_REF] Zhang | A modular-type formula for (x; q)∞[END_REF]Th. 2.9], by means of one certain perturbed factor named P (z, τ ), where x = e(z) and q = e(τ ). Thus, it suffices to understand the analyticity obstacle of P (α + β τ, τ ) around each given rational point τ = r ∈ Q ∩ [0, 1). We shall obtain the condition for this function to be analytically continued at τ = 0 by a Stokes analysis, with the help of the Ramis-Sibuya Theorem [START_REF] Malgrange | Sommation des séries divergentes[END_REF], [START_REF] Ramis | Gevrey Asymptotics and Applications to Holomorphic Ordinary Differential Equations[END_REF]. This analysis will be generalized for every r ∈ Q ∩ (0, 1) by means of a series of transformations associated to the continued fraction of r; transformations often used in the classical theory of the modular functions.

1.3.

Plan for the paper. The rest of the paper is divided into three sections. In Sec. 2, we define a family of integrals involving the exponential generating function associated with the Bernoulli numbers. These integrals can be seen as being of Laplace type, and they will be used for stating an equivalent version of the abovementioned result on (x; q) ∞ ; see Theorems 2.1 and 2.2.

Sec. 3 is essentially devoted to the part ζ = 1 of Theorem 1.2; see Theorem 3.1. By means of Theorems 2.1 and 2.2, we will see that the fact that a Euler q-exponential function, modulo some exponentially small term, can be analytically continued at τ = 0 and may be interpreted as one problem of the analytic continuation inside the theory of the Gevrey asymptotic expansions; see Theorem 3.4 and the proof of Theorem 1.2 given in §3.3.

Sec. 4 aims to obtain Theorem 1.2 for an arbitrary root ζ of unity; see Theorem 4.1. Lemma 4.1 will play a key role, especially in terms of permitting us to make use of both continued fractions and modular transforms. Finally, a complete scheme for proving our main result, Theorem 1.1, will be outlined at the end of the paper.

A Laplace-type integral involving Bernoulli exponential generating function

The goal of this section is to develop appropriate means for properly understanding the following result obtained in [START_REF] Zhang | On the modular behaviour of the infinite product (1-x)(1-xq)(1-xq 2 )(1-xq 3 )[END_REF] and [START_REF] Zhang | A modular-type formula for (x; q)∞[END_REF]: Theorem 2.1 ( [13, Th. 3.2], [15, Th. 2.9]). Let (z, τ ) ∈ U and let s = z/τ . If s / ∈ (-∞, 0], then

(2.1) (z | τ ) ∞ = 2πs(1 -e(z)) Γ(s + 1) e(- τ 24
) e s(log s-1)+ Li 2 (e(z))

2πiτ

+P (z,τ ) ( z -1 τ | - 1 τ ) ∞ ,
where P (z, τ ) denotes the analytic function in U defined by the following integral:

(2.2)

P d (z, τ ) = ∞e id 0 sin( zt τ ) e it/τ -1 cot t 2 - 2 t dt t (-π < d < 0) .
In the above, U denotes the domain defined in C × H by the relation

(2.3) U = ∪ δ∈(0,π) C δ × H δ ,
where

H δ = {τ ∈ H : arg τ ∈ (0, δ)} and C δ = C \ (1 + H δ ) ∪ (-1 -H δ ) .
The functions Γ, log and Li 2 are the Euler Gamma function, the principal branch of the complex logarithm function and the dilogarithm function, respectively. In §2.1, we will introduce a family of Laplace-like integrals denoted as b d , involving the exponential generating function of Bernoulli numbers. It will be shown that the term ((z -1)/τ | -1/τ ) ∞ in the right-hand side of (2.1) can be obtained from comparing these integrals in different directions; see Theorem 2.2. One will also see that the same integrals are closely linked to the function P (z, τ ) used in Theorem

I(x) := ∞ 0 B(t) e -xt dt t = log Γ(x) -(x - 1 2 ) log x + x - 1 2 log 2π .
Here, first of all, one supposes that x > 0, so the integration path is the half-axis (0, +∞). By using an open interval (0, ∞e id ) in the half-plane t > 0, this integral representation can then be valid for all x ∈ C \ (-∞, 0]. The integral I(x) stated in (2.6) is the Laplace transform of the function t → B(t)/t. In what follows, we shall consider a modified Laplace-type integral b

d (z, τ ) associated to each d ∈ (-π 2 , π 2 ):
(2.7)

b d (z, τ ) = ∞e id 0 e -zu -1 e u -1 B(τ u) du u .
To be brief, b d will be called a Bernoulli integral.

Let us determine the values (z, τ ) ∈ C 2 where the integral b d (z, τ ) is well-defined. From (2.5), it follows that B(t) = O(t) at t = 0 in C. Thus, the above integral in (2.7) converges at u = 0 for all (z, τ ) ∈ C × C. With regard to the convergence at infinity, we define (2.8)

V d,+ = {τ ∈ C : (τ e id ) > 0} , V d,-= {τ ∈ C : (τ e id ) < 0} and
(2.9)

U d = {z ∈ C : z + 1 ∈ V d,+ } = V d,+ -1 .
By using (2.4), one finds that

B(t) → ±1/2 when t → ±∞. Therefore, b d (z, τ ) is defined in two separated domains U d × V d,+ and U d × V d,-in C 2 .
Geometrically, U d represents the half-plane containing the point at origin and delimited by the straight-line -1 + e i(-d+ π 2 ) R, while V d,± are half-planes separated by the straight line e i(-d+ π 2 ) R. One can find that the interval (-1, ∞) belongs to

U d for every argument d ∈ (-π 2 , π 2 ); see Figure 1. Let (2.10) W + = ∪ d∈(-π 2 , π 2 ) U d × V d,+ , W -= ∪ d∈(-π 2 , π 2 ) U d × V d,-. Since ∪ d∈(-π 2 , π 2 ) V d,+ = C \ (-∞, 0], ∪ d∈(-π 2 , π 2 ) V d,-= C \ [0, ∞), it follows that (2.11) (-1, ∞) × (C \ (-∞, 0]) ⊂ W + , (-1, ∞) × (C \ [0, ∞)) ⊂ W -. E -1 0 1 ∞ ¨¨¨B r r r r j j j d -d -d + π 2 V d,+ V d,- U d ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ Figure 1. Half-planes U d , V d,-and V d,+ Definition 2.1. We define b + (z, τ ) and b -(z, τ ) in W + and W -, respectively, by applying the analytic continuation procedure to b d (z, τ ) from U d × V d,+ and U d × V d,-as d runs through (-π 2 , π 2
). We shall make use of the following result to express the difference between b + (z, τ ) and b -(z, τ ) in their common domain

W + ∩ W -: Lemma 2.1. If τ ∈ H and z ∈ H, then (2.12) n≥1 1 n e(nz) 1 -e(nτ ) = -log ((z | τ ) ∞ ) .
Proof. This follows from [2, p. 36, (

.

By (2.11), one finds that ((-1, ∞) × (C \ R)) ⊂ W + ∩ W -.
In what follows, we will write

C \ R = H ∪ H -, where H -= -H = {τ ∈ C : τ < 0}. Theorem 2.2. Let (z, τ ) ∈ W + ∩ W -.
The following assertions hold:

(1) if τ ∈ H, then

(2.13) b + (z, τ ) -b -(z, τ ) = -log (-(z + 1)/τ | -1/τ ) ∞ (-1/τ | -1/τ ) ∞ ; (2) if τ ∈ H -, then (2.14) b + (z, τ ) -b -(z, τ ) = -log ((z + 1)/τ | 1/τ ) ∞ (1/τ | 1/τ ) ∞ .
Proof. (1) By the standard argument of analytical continuation, it suffices to prove (2.13) for (z, τ ) ∈ (-1, ∞) × H. Thus, one chooses

d 1 ∈ (-π 2 , 0) and d 2 ∈ (0, π 2 ) such that τ ∈ V d1,+ ∩ V d2,-. The contour integral in (2.7) allows one to write that b + (z, τ ) -b -(z, τ ) = b d1,+ (z, τ ) -b d2,-(z, τ ) = ∞e id 1 0 - ∞e id 2 0 e -zu -1 e u -1 B(τ u) du u . (2.15) Since both d 1 and d 2 belong to (-π 2 , π 2
), the two half straight-lines used in the contour-integral (2.15) are separated in the u-plane by the half straight-line τ defined by the relation τ = {u ∈ C * : (τ u) = 0, u > 0}; see Figure 2.

By observing that the function B(τ u) admits simple poles u = 2nπi/τ (n ∈ Z >0 ) on the line τ , applying the Residues Theorem to (2.15)

yields that b + (z, τ ) -b -(z, τ ) = n≥1 1 n e(-nz/τ ) -1 e(n/τ ) -1 E 0 £ £ £ # • τ $ $ $ $ $ $ τ : poles of B(τ u) w d 1 G G V d1,+ Q d 2 s s V d2,-
+ (z, τ ) -b -(z, τ ) = -log (- z + 1 τ | - 1 τ ) ∞ + log (- 1 τ | - 1 τ ) ∞ ,
so that one obtains (2.13).

(2) When τ ∈ H -, the above proof can be adopted as follows: choose d 1 ∈ (0, π 2 ) and d 2 ∈ (-π 2 , 0), and observe that the simple poles of B(τ u) to which the Residues Theorem is applied become u = -2nπi/τ (n ∈ Z >0 ). A direct calculation implies (2.14), which ends the proof of Theorem 2.2. Now, consider τ ∈ H, with arg τ = δ ∈ (0, π). By (2.10), it follows that (z, τ ) ∈ W + if and only if z ∈ V d,+ for some suitable d ∈ (-π 2 , π 2 -δ). Thus, one obtains the equivalence

(2.17) (z, τ ) ∈ W + ⇐⇒ z ∈ H ∪ Z τ ,
where Z τ is the half-plane associated with τ in the following manner:

Z τ = {z ∈ C : z + 1 τ < 0} .
One may see that if z ∈ Z τ , then z + z 0 ∈ Z τ for all z 0 ∈ R >0 ; see Figure 3.

E -1 0 1 ∞ ¨¨¨¨B $ $ X•τ £ £ £ £ £ £ y y d H π 2 -d π 2 -δ V d,+ Z τ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ δ = arg τ -π 2 < d < π 2 -δ Figure
P d (z, τ ) = ∞e id 0 e zt/τ -e -zt/τ e t/τ -1 B(t) dt t ,
where we write

d = d + π 2 ∈ (-π 2 , π 2
). Let W + be as in (2.10), and let U be as in (2.3). A simple computation shows that

U = {(z, τ ) ∈ W + : τ ∈ H, (-z, τ ) ∈ W + } . Furthermore,
H d = V d,+ - 1 2 ∩ -V d,+ + 1 2 .
It is easy to see that the integral

(2.22) B d (z, τ ) = ∞e id 0 e zu -e -zu e u/2 -e -u/2 B(τ u) du u is well-defined for any (z, τ ) ∈ H d × (V d,+ ∪ V d,-
). Furthermore, by noticing that

e zu -e -zu e u/2 -e -u/2 = e (z+ 1 2 )u -1 e u -1 - e -(z-1 2 ) -1 e u -1 , comparing (2.22) with (2.7) yields that (2.23) B d (z, τ ) = b d (-z - 1 2 , τ ) -b d (z - 1 2 , τ ) . Definition 2.2. Let (2.24) Ω + = ∪ d∈(-π 2 , π 2 ) H d × V d,+ , Ω -= ∪ d∈(-π 2 , π 2 ) H d × V d,-.
We will denote by B + (z, τ ) and B -(z, τ ) the respective functions defined in Ω + and Ω -by the integral (2.22).

If there is no possible confusion, we will simply write Ω and B(z, τ ) instead of Ω + and B + (z, τ ), respectively.

From (2.23), it follows that (2.25) B ± (z, τ ) = b ± (-z - 1 2 , τ ) -b ± (z - 1 2 , τ ) .
Furthermore, combining this last equality with (2.20) and (2.18) yields that

(2.26) B(z, τ ) = P (z + 1 2 , τ ) -I( z + 1/2 τ ) .
Theorem 2.3. The following relation holds for all (z, τ ) ∈ Ω with τ ∈ H:

(2.27) (z + 1 2 | τ ) ∞ = e(- τ 24 ) 1 + e(z) e Li 2 (-e(z)) 2πiτ +B(z,τ ) ( z -1/2 τ | - 1 τ ) ∞ .
Proof. By considering (2.6), the formula in (2.1) can be put into the following form:

(z | τ ) ∞ = 1 -e(z) e(- τ 24
) e

Li 2 (e(z)) 2πiτ

-I( z τ )+P (z,τ ) ( z -1 τ | - 1 τ ) ∞ .
Thus, one obtains (2.27), with the help of (2.26).

3. Conditions for a Euler q-exponential function to be of theta-type at one Let x 0 = e(α) = e 2πiα with α ∈ R, and let β ∈ R. By using (1.3) with N = ∞, we write (x 0 q β ; q) ∞ = (α + β τ | τ ) ∞ , where τ ∈ H and q = e(τ ). The goal of this section is to establish the next result, which will be useful for the proof of Theorem 1.2.

Theorem 3.1. Let (α, β) ∈ [0, 1)×(0, 1], and consider f (q) = (α+β τ | τ ) ∞ . Then, f ∈ C 1 if and only if α ∈ {0, 1 2 } and β ∈ { 1 2 , 1}.
The main idea here will consist of using Theorems 2.1 and 2.3 to rewrite f (q) in such a way that log f (q) = B(...) or = P (...) mod C{τ } ⊕ A ≤-1 1 , the exponential small term being furnished by an infinite product of e(-1/τ ). Thus, we will be led to consider the analytic continuation of B or P around τ = 0 in the complex plane; see Theorem 3.4 in §3.2. In this way, we will obtain the condition for (α, β) required by Theorem 3.1, whose proof will be completed in §3.3.

We shall make use of the Gevrey asymptotic expansions for understanding the analytic obstacle at τ = 0 of the above-mentioned functions B and P . This is linked to the so-called Stokes' phenomenon. One tool to treat this problem may be Ramis-Sibuya Theorem, which will be briefly explained in §3.1 in what follows.

3.1. Ramis-Sibuya's Theorem on Gevrey asymptotic expansions. Let x 0 ∈ C and let Cx0 be the Riemann surface of the function x → log(x -x 0 ); let I = (α 1 , α 2 ) ⊂ R and let R > 0. We let V x0 (I; R) denote the sector of a vertex at x 0 in Cx0 , with an opening in I and a radius R; that is to say,

(3.1) V x0 (I; R) = {x 0 + re iα : α ∈ I, r ∈ (0, R)} .
By definition, a proper sub-sector of V x0 (I; R) will be any domain of the form V x0 (J; ρ) such that J ⊂ I and ρ < R.

If the length of the open interval I is smaller than or equal to 2π, any sector V x0 (I; R) is not overlapped in Cx0 ; in this case, one will consider V x0 (I; R) as a sector in C. When x 0 = 0, we will remove the sub-index 0 and simply write V (I; R) instead of V 0 (I; R).

Let V = V (I; R) be a sector in C at 0. By definition ( [START_REF] Malgrange | Sommation des séries divergentes[END_REF], [START_REF] Ramis | Gevrey Asymptotics and Applications to Holomorphic Ordinary Differential Equations[END_REF], . . . ), a given function f defined and analytic in V is said to have a power series n≥0 a n x n , a n ∈ C, as a Gevrey asymptotic expansion at 0 in V , if, for any proper sub-sector U = V (J; ρ), one can find C > 0 and A > 0 such that the following estimates hold for all n ∈ Z ≥0 :

(3.2) sup x∈U (f (x) - n-1 m=0 a m x m )x -n ≤ C A n n! .
As a typical example, the Borel-sum function of a given divergent series, if it exists, admits this series as a Gevrey asymptotic expansion. A Gevrey type asymptotic expansion is also called an exponential asymptotic expansion, due to the following fact: Remark 3.1. [6, p. 175, Th. 1.2.4.1 1)] A function f admits the identically null series as a Gevrey asymptotic expansion at 0 in V if and only if f is exponentially small there, which means that, for all proper sub-sectors U in V , there exists C > 0 and κ > 0 such that, for all

x ∈ U , |f (x)| ≤ C e -κ/|x| .
In what follows, we will denote by A ≤-1 (V ) the space of all functions that are exponentially small in V as indicated in Remark 3.1. More generally, when V = V x0 (I; R), we will say that f ∈ A ≤-1 (V ) when f is exponentially small as

x → x 0 in V . Theorem 3.2. [6, p. 176, Th. 1.3.2.1] Let V 1 , ..., V m , V m+1 be a family of open sectors at 0 in C such that V m+1 = V 1 , V j ∩ V j+1 = ∅ for 1 ≤ j ≤ m
and that the whole union ∪ m j=1 V j contains a neighborhood of 0 in C. For every j, let f j be a given analytic and bounded function in V j . If

f j+1 = f j mod A ≤-1 (V j ∩ V j+1 ) ,
then all f j 's admit the same Gevrey asymptotic expansion at 0. The above result is currently called Ramis-Sibuya's Theorem. We shall make use of the following statement deduced from Theorem 3.2: Corollary 3.1. Let R > 0, and let I 1 and I 2 be open intervals such that

[-, π -] ⊂ I 1 ⊂ (-π, π), [π -, 2π -] ⊂ I 2 ⊂ (0, 2π)
for some ∈ (0, π). Let V 1 = V (I 1 ; R), V 2 = V (I 2 ; R), and consider two analytic and bounded functions f 1 and f 2 defined, respectively, in

V 1 and V 2 . If f 1 -f 2 ∈ A ≤-1 (V 1 ∩ V 2 )
, then f 1 and f 2 have the same Gevrey asymptotic expansion and, moreover, the following conditions are equivalent:

(1) one of the functions f 1 and f 2 can be continued into an analytic function at 0 in C; (2) both f 1 and f 2 can be continued into an analytic function at 0 in C;

(3)

f 1 ≡ f 2 in V 1 ∩ V 2 .
Proof. The existence of a Gevrey asymptotic expansion for f 1 and f 2 follows immediately from Theorem 3.2.

Let f = n≥0 a n x n be the common asymptotic expansion of f 1 and f 2 . Since the length of I 1 and that of I 2 are larger that π, one finds that f 1 and f 2 are the respective Borel-sum functions of f in V 1 and V 2 . Thus, the above statement in [START_REF] Andrews | Special functions, Encyclopedia of Mathematics and its Applications[END_REF] implies that f is really a convergent series, so that their two Borel-sums are equal to each other. In this way, one obtains that (1) implies all other statements.

On the other hand, if the statement in is true, then both f 1 and f 2 equal to a same analytic and bounded function in the punctuated disc {0 < |x| < R}. By the Riemann removable singularities Theorem, one finds the statements (1) and (2).

Asymptotic expansion of Bernoulli integrals.

From now on, we will identify the upper half-plane H as the sector V (I; R) with I = (0, π) and R = ∞. Thus, A ≤-1 (H) will be the space of all analytic functions in H that are exponentially small as τ → 0. It is easy to see that A ≤-1 (H) ⊂ A ≤-1

1

, where ζ = 1 = e 2πi0 with r = 0; see §1.2 (iii). Proposition 3.1. Let (α, β) ∈ R 2 , and consider

f (τ ) = log (- α + β τ τ | - 1 τ ) ∞ for τ → 0 in H. If α > 0, then f ∈ A ≤-1 (H).
Proof. Thanks to Euler [1, p. 490, Corollary 10.2.2 (b)], one can write that, for all x ∈ C,

(3.3) (x; q) ∞ = n≥0
q n(n-1)/2 (q; q) n (-x) n .

Letting x = e(-(α + β τ )/τ ) and q = e( -1/τ ) into (3.3), one gets that

(- α + β τ τ | - 1 τ ) ∞ = 1 + n≥1 (-1) n e(-n(α + β τ )/τ ) ( -1/τ | -1/τ ) n e n(n -1) 2 ( - 1 τ ) ,
where (. | .) n is defined as in (1.3). Since e( -ν τ ) ∈ A ≤-1 (H) for any ν > 0, it follows that, when α > 0,

(- α + β τ τ | - 1 τ ) ∞ = 1 - e(-(α + β τ )/τ ) 1 -e( -1/τ ) mod A ≤-1 (H) = 1 mod A ≤-1 (H) .
This finishes the proof.

Proposition 3.2. Let (α, β) ∈ R 2 and let I = (-π, π). If α > -1, then b + (α + β τ, τ ) is well-defined and analytic in V (I; R), and is bounded in every proper subsector of V (I; R) with R > 0.

Proof. For all τ ∈ C * , let D τ be the sector containing 0 that is bounded by (-∞, -1] ∪ [-1, -1 -∞τ ), where [-1, -1 -∞τ ) denotes the half straight-line starting from -1 to ∞ with the direction -τ . By combining (2.8) together with (2.9), one can find that, for all fixed τ ∈ C\(-∞, 0], the function b + (z, τ ) is defined and analytic for z ∈ D τ . If α > -1, one can easily see that α + β τ belongs to this half-plane D τ when τ / ∈ R -. This implies that b + (α + β τ, τ ) is well-defined and analytic in any sector V (I; R).

The boundedness of this function over any proper sub-sector comes from direct estimates done for (2.7).

In a similar way, one can find that the statement of Proposition 3.2 remains true if b + (z, τ ) and I are replaced with b -(z, τ ) and (0, 2π), respectively. Thus, one obtains the following:

Theorem 3.3. Let (α, β) ∈ R 2 and let b + (z, τ ) as in Definition 2.1. If α > -1, then b + (α + β τ, τ
) admits a Gevrey asymptotic expansion in any sector V (I; R) with I = (-π, π) and R > 0.

Moreover, b + (α + β τ, τ ) can be continued into an analytic function at τ = 0 if and only if α = 0 and β ∈ Z.

Proof. Fix R > 0, and let

V 1 = V ((-π, π); R), V 2 = V ((0, 2π); R). Define f 1 (τ ) = b + (α + β τ, τ ) = b + (α + β τ, τ ), f 2 (τ ) = b -(α + β τ, τ )
for τ ∈ V 1 and V 2 , respectively. By putting z = α + β τ into both relations (2.13) and (2.14) of Theorem 2.2, it follows that

(3.4) f 1 (τ ) -f 2 (τ ) = -log (-(α + β τ + 1)/τ | -1/τ ) ∞ (-1/τ | -1/τ ) ∞ , if τ ∈ H ∩ V 1 ∩ V 2 , and that (3.5) f 1 (τ ) -f 2 (τ ) = -log ((α + β τ + 1)/τ | 1/τ ) ∞ (1/τ | 1/τ ) ∞ , if τ ∈ H -∩ V 1 ∩ V 2 . One observes that V 1 ∩ V 2 ∩ R = ∅.
Therefore, by considering Proposition 3.1, relation (3.4) together with (3.5) imply that f 1 (τ ) -f 2 (τ ) is exponentially small in the common domain V 1 ∩ V 2 . This allows us to apply Corollary 3.1 to get a particularly the common Gevrey asymptotic expansion of both f 1 and f 2 .

Furthermore, Corollary 3.1 implies that f 1 can be extended into an analytic function at τ = 0 in C if and only if

(- α + β τ + 1 τ | - 1 τ ) ∞ = (- 1 τ | - 1 τ ) ∞
for all τ ∈ H, or, equivalently,

- α + β τ + 1 τ = - 1 τ mod Z .
In this way, one finds the necessary and sufficient condition α + β τ ∈ τ Z in order to have an analytic function b + (α + β τ, τ ) at τ = 0 in C. This ends the proof of Theorem 3.3. Now, consider the functions P (z, τ ) and B(z, τ ) appearing in Theorems 2.1 and 2.3. In keeping with the spirit of Theorem 3.3, one finds the following result: Theorem 3.4. Let (α, β) ∈ R 2 , and let V = V (I; R) with I = (-π, π) and R > 0.

Then

(1) if α ∈ (-1, 1), the function P (α + β τ, τ ) admits a Gevrey asymptotic ex-

pansion as τ → 0 in V ; (2) if α ∈ (-1 2 , 1 2
), the function B(α + β τ, τ ) admits a Gevrey asymptotic expansion as τ → 0 in V . Furthermore, P (α + β τ, τ ) or B(α + β τ, τ ) can be continued into an analytic function at τ = 0 in C if and only if α = 0 and β ∈ 1 2 Z. Proof. We shall consider only the function P (z, τ ), and the case of B(z, τ ) is very similar.

By using (2.20), P (α + β τ, τ ) can be expressed in terms of b + (α + β τ, τ ) as follows:

P (α + β τ, τ ) = -b + (α + β τ, τ ) + b + (-α -β τ, τ
) . Thus, Theorem 3.3 implies that P (α + β τ, τ ) remains analytic and has a Gevrey asymptotic expansion as τ → 0 in V when α ∈ (-1, 1).

Furthermore, combining the expression of P (z, τ ) in (2.20) with Stokes's relations (2.13) and (2.14) allows one to obtain the following equation: for all τ ∈ H,

P (z, τ ) -P -(z, τ ) = log (-(z + 1)/τ | -1/τ ) ∞ (-(-z + 1)/τ | -1/τ ) ∞ ,
where P -(z, τ ) denotes the function defined by (2.19) with d ∈ (0, π). Thus, P (α + β τ, τ ) can be continued into an analytic function at τ = 0 if and only if

-(α + β τ + 1)/τ = (α + β τ -1)/τ mod Z .
This achieves the proof of Theorem 3.4.

Proof of Theorem 3.1.

In what follows, we will denote by C{τ } the space of the germs of analytic functions at τ = 0 in H. One knows that the dilogarithm is well-defined and analytic in the universal covering of C \ {1}. Thus, u → Li 2 (e(u)) represents an analytic function on the Riemann surface of logarithm, i.e, the universal covering C0 of C \ {0}.

Lemma 3.1. The following relation holds for all u ∈ C0 :

(3.6) Li 2 e(ue 2πi ) -Li 2 (e(u)) = 4π 2 u .
Proof. Let x = e(u) for u ∈ C0 . When u makes a complete rotation along a circle around u = 0, the corresponding x forms a circle around x = 1. By using a relation between Li 2 (x) and Li

2 (1 -x) [12, §2], Li 2 (1 -x) = -Li 2 (x) + π 2 6 -log x log(1 -x),
one finds that the monodromy of Li 2 around x = 1 can be expressed as follows:

Li 2 (1 + x e 2πi ) = Li 2 (1 + x) -2πi log(1 + x) .

Therefore, one gets that

Li 2 e ue 2πi -Li 2 (e (u)) = -2πi log (e (u)) , which implies the desired relation (3.6).

Proof of Theorem 3.1. First of all, suppose that (α,

β) ∈ {0, 1 2 ) × { 1 2 , 1}. It follows from Remark 1.1 that f ∈ T 1 , so f ∈ C 1 , also.
Now, consider the "only if" part, and suppose that f ∈ C 1 . The rest of the proof will be divided into two parts, according whether α may be null or not.

• Case 1: α ∈ (0, 1). Let α = α -1 2 , and observe that α ∈ (-1 2 , 1 2 ). By putting z = α + β τ into (2.27) of Theorem 2.3, it follows that

(3.7) log f (q) = A(τ ) + L(τ ) + B(α + β τ, τ ) + R(τ ),
where one introduces the following notation:

A(τ ) = log e - τ 24 1 + e(α + β τ ) = π 12 τ i + 1 2 log (1 -e(α + β τ )) , L(τ ) = Li 2 (e(α + β τ )) 2πiτ , R(τ ) = log( α + β τ -1 τ | - 1 τ ) ∞ .
On the one hand, as α -1/2 < 0, Proposition 3.1 implies that R(τ

) ∈ A ≤-1 (H) ⊂ A ≤-1 1
. On the other hand, it is easy to see that

A(τ ) ∈ C{τ }, L(τ ) = Li 2 (e(α)) 2πiτ mod C{τ } . Thus, comparing (1.5) (ζ = 1, r = 0, τ = τ ) with (3.7) yields that (3.8) B(α + β τ, τ ) = ν log τ i + λ τ mod C{τ } ⊕ A ≤-1 1 ,
where ν ∈ Q and λ ∈ C. By Theorem 3.4, it follows that B(α +β τ, τ ) has a Gevrey asymptotic expansion as τ → 0 in any sector V = V (I; R), where I = (-π, π) and R > 0. This implies that ν = 0 and λ = 0 in (3.8). Furthermore, the exponentially small term used in (3.8) will be bounded in any proper sub-sector of V . As the openness of V is larger than π, a classical argument such as the Phragemen-Lindeloff Theorem implies that this term is identically null; see [START_REF] Fruchard | Remarques sur les développements asymptotiques[END_REF] for more on this matter. Thus, one gets that B(α + β τ, τ ) can be really continued into an analytic function at τ = 0. Applying Theorem 3.4 (2) implies that α = 0 and β ∈ 1 2 Z, so it follows that α = 1 2 and

β ∈ { 1 2 , 1}. • Case 2: α = 0. Putting z = β τ and s = z/τ = β > 0 into (2.1) of Theorem 2.1 gives that f (q) = 2πβ(1 -e(β τ )) Γ(β + 1) e(- τ 24 ) e β(log β-1)+L(τ )+P (β τ,τ ) R 1 (τ ) , so that (3.9) log f (q) = 1 2 log τ i -I(β) + A(τ ) + L(τ ) + P (β τ, τ ) + R(τ ) .
In the above, I denotes the function given by (2.6),

A(τ ) = π 12 τ i + 1 2 log e (β τ ) -1 i β τ , L(τ ) = Li 2 (e(β τ )) 2πiτ , and
R(τ ) = log R 1 (τ ), R 1 (τ ) = ( β τ -1 τ | - 1 τ ) ∞ .
One can easily see that A(τ ) ∈ C{τ }. Letting u = β τ in (3.6) gives that L(τ e 2πi ) -L(τ ) = -2βπi. In view of the equality Li 2 (1) = π 2 6 , the function L can be put into the following form:

L(τ ) = c ∞ τ + β log τ i mod C{τ } , c ∞ = - π 12 i .
In addition, by Proposition 3.1, one gets that R ∈ A ≤-1

1

. Thus, it follows from (3.9) that

(3.10) log f (q) = c ∞ τ + ( 1 2 + β) log τ i + P (β τ, τ ) mod C{τ } ⊕ A ≤-1 1 .
One knows that P (β τ, τ ) admits a Gevrey asymptotic expansion as τ → 0 in C \ (-∞, 0]. As in Case 1 for B(..., τ ), comparing (3.10) with (1.5) gives that P (β τ, τ ) can be continued into an analytic function at τ = 0. Thus, applying Theorem 3.4 (1) implies, finally, that β ∈ 1 2 Z, whis gives that β ∈ { 1 2 , 1}. • In summary, one finds that f ∈ C 1 implies that α ∈ {0, 1 2 } and β ∈ { 1 2 , 1}. This ends the proof of Theorem 3.1.

Asymptotic behavior at an arbitrary root via continued fractions

With regard to an arbitrary root ζ of unity, we shall establish the following result, which, together with Theorem 3.1, will imply Theorem 1.2: 1 2 } and β ∈ { 1 2 , 1}. First, one will observe, in §4.1, that the corresponding functions B and P used in Theorems 2.3 and 2.1 are analytic at each non-zero rational point τ = r. This allows us to establish one key lemma, Lemma 4.1, in §4.2, that permits us to pass an arbitrary rational value r to an other r 1 . By iterating this procedure, one arrives at the case of r = 0, to which case Theorem 3.1 can be applied. This is realized in terms of the continued fractions relative to r and related modular transforms; see Theorem 4.4 in §4.3. We complete the proofs of Theorems 4.1, 1.2 and 1.1.

Theorem 4.1. Let r ∈ Q ∩ (0, 1), ζ = e(r) and (α, β) ∈ [0, 1) × (0, 1], and consider f (q) = (α + β τ | τ ) ∞ . Then f ∈ C ζ if and only if α ∈ {0,

4.1.

Bernoulli integral and associated functions on a real axis. We will discuss the degenerate case τ ∈ R >0 for the functions b + (z, τ ), B(z, τ ) and P (z, τ ). In what follows, we will make use of the notational convention

(4.1) ∈ (0, π 2 ), W = V ((-. ); ∞), W c = C \ W ,
and the letter r always denotes a given positive number.

First of all, we consider the function b + (z, τ ). It should be noted that the relation stated in (2.18) is valid for any (z, τ ) ∈ H × H. If z / ∈ H, we have to avoid the poles of the Gamma function, and the right-hand side of (2.18) continues to be well-defined over the Riemann surface of log while s / ∈ Z ≤0 . Thus, Lemma 2.2 allows one to make the analytic continuation of the function b + at (z, τ ) provided that (z + n)/τ / ∈ Z ≤0 for all n ∈ Z ≥1 . This yields the following observation:

Remark 4.1. For any fixed τ ∈ H, z → b + (z, τ ) can be continued into an analytic function on the universal covering of C \ ∆ τ , where

(4.2) ∆ τ = Z ≤-1 ⊕ τ Z ≤0 .
By Definition 2.1 and Remark 4.1, b + (z, τ ) is well-defined and analytic in the domain (-W c -1) × W , where -W c -1 = C \ -W -1 ; see Figure 4 below. In particular, we note the following fact:

H W -W -1 u -W c -1 Figure 4. b + (z, τ ) is analytic for z ∈ (-W c -1) and τ ∈ W When → 0 + , W becomes (0, ∞) and -W c -1 is reduced into C \ (-∞, -1].
By replacing τ with r in the partial lattice ∆ τ given by (4.2) for all τ ∈ H, we will continue to write ∆ r = {n + mr : n ∈ Z ≤-1 , m ∈ Z ≤0 }. It is easy to see that ∆ r is discrete on the real axis if and only if r is a rational number. In this way, we shall make use of the following remark:

(4.3) ∆ r = -1 + 1 d Z ≤0 ⊂ (-∞, -1] .
For any α ∈ R, let [α] denote the integral part of α and {α} denote the corresponding fractional part, that is α -

[α]. If z 0 ∈ (-∞, -1] \ ∆ r , we define (4.4) b(z 0 , r) =   [-z 0 ] -1 -2 [-z0]-1 k=1 z 0 + k r   πi .
Theorem 4.2. Let r ∈ (0, 1) ∩ Q and let ∆ r , b be as in the above. Then b + (z, r) is analytic in C \ (-∞, -1] and can be continued to be an analytic function over the universal covering of C \ ∆ r in such a way that the following relation holds for all

z 0 ∈ (-∞, -1] \ ∆ r : (4.5) lim →0 + b + (z 0 + i , r) -b + (z 0 -i , r) = b(z 0 , r) .
Proof. As τ → r in H, the limit set ∆ r of the singularities of b + (z, τ ) is discrete, as stated in (4.3). By considering Remark 4.1, one obtains that b + (z, r) is analytic over the universal covering of C \ ∆ r . Now, let D(z 0 ) denote the expression in the left-hand side of (4.5). By putting τ = r into (2.18), one finds that, if z / ∈ (-∞, 0], then b + (z -1, r) = b + (z, r) + I(s), where s = z/r and I(s) is given in (2.6). Thus, one can write

D(z 0 -1) = D(z 0 ) + lim →0 + I( z 0 r + i) -I( z 0 r -i)
for all z 0 ∈ (-∞, 0] \ ∆ r . By using (2.6), one gets that

(4.6) D(z 0 -1) = D(z 0 ) + (1 - 2z 0 r )πi + lim →0 + log Γ( z 0 r + i) -log Γ( z 0 r -i) .
Let n be any negative integer, say, n = -m, m > 0, and let s ∈ (n, n + 1). From the relation Γ(x) = Γ(x + m) (x) m and the fact that log Γ(x + m) is well-defined and analytic for (x) > n, it follows that

lim →0 + (log Γ(s + i) -log Γ(s -i)) = -lim →0 + (log(s + i) m -log(s -i) m ) = 2πin .
Therefore, (4.6) yields that (4.7)

D(z 0 -1) = D(z 0 ) + 2 n - z 0 r πi + πi .
Letting z 0 = (n + δ 0 )r with n 0 ∈ Z <0 and δ 0 ∈ (0, 1), the above expression in (4.7) becomes D(z 0 -1) -D(z 0 ) = -2πiδ 0 + πi. By replacing z 0 with z 0 + 1 and iterating this process, one obtains the finite sequences (n k ) and (δ k ) associated with the pair (z 0 , r) in the following manner:

z 0 + k = (n k + δ k )r, where n k ∈ Z <0 and δ k = δ k (z 0 , r) = { z0+k r }. Since D(z 0 + k) = 0 for k > -1 -z 0 , one finds that D(z 0 -1) = -2 δ 0 + ... + δ [-z0]-1 πi + [-z 0 ] πi .
Replacing z 0 with z 0 + 1 in this last relation gives D(z 0 ) = b(z 0 , r), where b is as given in (4.4), so one obtains the expected relation (4.5) and Theorem 4.2.

By using (2.25), one finds that B(z, τ ) is analytic in the domain

((1/2 + W c ) ∩ (-1/2 -W c )) × W .
In addition, from (2.20) one obtains that P (z, τ ) can be continued to be analytic in the domain ((-1 -W c ) ∩ (1 + W c )) × W ; see Figure 5 Moreover, things are the same for τ → P (α + β τ, τ ) when α + βr ∈ (-1, 1).

By letting → 0 + , one sees that, for any r > 0, B(z, r) is analytic for all z ∈ C\(-∞, -1/2)∪(1/2, ∞), while P (z, r) is analytic for z ∈ C\(-∞, -1]∪[1, ∞). Theorem 4.3. Let r ∈ (0, 1) ∩ Q, and let ∆ r and b be as in (4.3) and (4.4). Then

(1) the function B(z, r) can be continued to be analytic in the universal covering of C \ ((1/2 + ∆ r ) ∪ (-1/2 -∆ r )), and the following relations hold for all

z 0 ∈ R \ ((1/2 + ∆ r ) ∪ (-1/2 -∆ r )): (4.8) lim →0 + (B(z 0 + i , r) -B(z 0 -i , r)) = -b(-|z 0 | - 1 2 , r) ;
(2) the function P (z, r) can be continued to be analytic in the universal covering of C \ (∆ r ∪ (-∆ r )) in such a way that, for all z 0 ∈ R \ (∆ r ∪ (-∆ r )), (4.9) lim →0 + (P (z 0 + i , r) -P (z 0 -i , r)) = -b(-|z 0 |, r) . Proof. This follows directly from Theorem 4.2 together with relations (2.25) and (2.20).

4.2.

One key lemma. As in the definition of b(z, r) in (4.4), we will let [a] and {a} denote the integral and fractional part, respectively, of any given real a. Given each non-zero real r, consider the associated one-to-one map T r defined on [0, 1) × (0, 1] as follows:

(4.10) T r : (α, β) → ({(1 -α) - 1 r + β}, 1 -α) .
One finds easily that

(4.11) T r ({0, 1 2 } × { 1 2 , 1}) = {0, 1 2 } × { 1 2 , 1} .
f (q) = (z(τ ) | τ ) ∞ , g(q 1 ) = (z 1 (τ 1 ) | τ 1 ) ∞ ,
where q = e(τ ), τ 1 = -1 τ --1 r , q 1 = e(τ 1 ) and (4.12) -→ r ⇐⇒ τ 1 a.v.

z 1 (τ 1 ) = α 1 + β 1 τ 1 , (
-→ r 1 .

By observing that τ = -1/(τ 1 +N ), one gets that (z(τ )-1)/τ = (1-α)(τ 1 +N )+β, so it follows from (4.12) that z 1 (τ 1 ) = (z(τ ) -1)/τ mod Z. One remembers that

(z | τ ) ∞ = (z | τ ) if (z, τ ) = (z , τ ) mod Z 2 ; see (1. 3
). Thus, one finds that g(q 1 ) = g(q) if one defines

(4.14) g(q) = z(τ ) -1 τ | - 1 τ ∞ .
By noticing the relation τ 1 -r 1 = (τ -r)/(rτ ), it follows from (4.13) that g(q 1 ) ∈ C ζ1 if and only if g(q) ∈ C ζ . Thus, we shall use Theorems 2.3 and 2.1 to link f (q) with g(q) in the following fashion: (4. [START_REF] Zhang | A modular-type formula for (x; q)∞[END_REF] f (q) = H(q) g(q) .

By hypothesis, (α, β) / ∈ Z×{0}, so f (q) is not identically null. As C ζ \{0} constitutes a multiplicative group, Lemma 4.1 says exactly that H ∈ C ζ , which is what we need to establish.

As before, write τ = τ -r. As in the proof of Theorem 3.1, we shall distinguish two cases: α + βr / ∈ Z and α + βr ∈ Z. • Case 1: α + βr -n ∈ (0, 1) with some n ∈ Z. Let α = α -n -1 2 , and observe that α + βr ∈ (-1 2 , 1 2 ). By applying Theorem 2.3 to z = z(τ ) -n -1 2 , the factor H(q) defined by (4.15) can be written as follows: H(q) = H 1 (q) H 2 (q) H 3 (q), where

H 1 (q) = 1 -e(z(τ )) e(- τ 24
) e

Li 2 (e(z(τ ))) 2πiτ

, H 2 (q) = e B(α +β τ,τ ) and (4.16)

H 3 (q) = ((z(τ ) -n -1)/τ | -1/τ ) ∞ ((z(τ ) -1)/τ | -1/τ ) ∞ .
When τ → r, it follows that e (z(τ )) → e(α + βr) = 1 for α + βr / ∈ Z. Thus Li 2 (e(z(τ ))) is really analytic at τ = r in C. As r = 0, one finds, finally, that H 1 ∈ C{τ }. In addition, Remark 4.4 implies that B(α + β τ, τ ) is analytic at τ = r. Furthermore, one can express H 3 as

H 3 (q) = z(τ ) -n -1 τ | - 1 τ -n or H 3 (q) = 1 ((z(τ ) -1)/τ | -1/τ ) n
for n ∈ Z ≤0 or Z ≥0 , respectively; see (1.3). This also shows that H 3 ∈ C{τ }. Thus, one gets that H ∈ (C{τ } \ {0} ⊂ C ζ .

• Case 2: α + βr = n ∈ Z. Let s = (z(τ ) -n)/τ , and notice that z(τ ) -n = β τ , so s = β τ /τ . As β = 0, one gets that s ∈ C \ R for all τ ∈ H. Moreover, one has C{s} = C{τ } for r = 0.

Putting z = z(τ ) -n = β τ into Theorem 2.1 gives that the factor H(q) defined by (4.15) can be read as follows: H(q) = H 0 (q) H 1 (q) H 2 (q) H 3 (q), where H 0 (q) = 2πs(1 -e(τ s)) Γ(s + 1) e(-τ 24 ) , H 1 (q) = e s(log s-1)+ Li 2 (e(τ s)) 2πiτ , H 2 (q) = e P (β (τ -r),τ ) , and where H 3 (q) is given by (4.16). It is immediately apparent that H 0 ∈ C{s}, so H 0 ∈ C{τ } also. Also, from Remark 4.4, one knows that H 2 ∈ C{τ }; it is the same for H 3 , as explained in the above. -→ r ∈ (0, 1) ∩ Q or, equivalently, when q a.r.

-→ ζ = e(r). Our strategy is to use continued fractions in order to reduce the general case τ a.v.

-→ r to the known case τ a.v.

-→ 0. Indeed, the above operation (r, τ, z) → (r 1 , τ 1 , z 1 ), considered in Lemma 4.1, allows one to link two rational numbers: r and r 1 . By iterating this process, one arrives at the case where τ a.v.

-→ 0. This iteration procedure requires one to write r into continued fraction. Thus, to any given r = p m ∈ Q ∩ (0, 1) will be associated the sequences r j ∈ Q ∩ [0, 1) and d j ∈ Z >0 in the following manner: (4.17) r 0 = p m , d 0 = 0; r j = -1 r j-1

, d j = -1 r j-1

( 1 ≤ j ≤ ν) ,

where ν denotes the smallest index such that r ν = 0, i.e., 1/r ν-1 ∈ Z >0 . With the standard notation for the continued fractions, one can notice that (4.18) r = [0, d 1 , -d 2 ..., (-1) ν-

1 d ν ] = 1 d 1 -| 1 d 2 -| • • • 1 |d ν .
Now, given (α 0 , β 0 ) ∈ [0, 1) × (0, 1] and r ∈ (0, 1) ∩ Q as in (4.18), define the r-depending sequence (z j , τ j ) 0≤j≤ν as follows: τ 0 = τ , z 0 = z 0 (τ 0 ) = α 0 + β 0 τ , (4.19) τ j = -1 τ j-1 -d j , z j = z j (τ j ) = α j + β j τ j , where (α j , β j ) = T rj-1 (α j-1 , β j-1 ), T rj-1 being the transform obtained by substituting r j-1 to r in (4.10). If τ = τ 0 a.v.

-→ r in H, then τ j a.v.

-→ r j in H, particularly with τ ν a.v.

-→ 0. Furthermore, it is easy to see that j ∈ H, with

τ j = M j τ j-1 , M j = -d j -1 1 0 ∈ SL(2; Z) ,
where τ → M τ denotes the classic modular transform associated with a modular matrix M . Thus, one can find that τ ν = M τ with M = M ν ... M 1 ∈ SL(2; Z).

Theorem 4.4. Let (ν, r), τ j and z j (τ j ) be given as in (4.18) and (4.19), with (α, β) ∈ [0, 1) × (0, 1]. Let ζ j = e(r j ) and q j = e(τ j ) for j from 0 to ν. Consider f (q) = (z 0 (τ 0 )|τ ) ∞ , f j (q j ) = (z j (τ j )|τ j ) ∞ .

Then, the following conditions are equivalent:

(1) f ∈ C ζ ;

(2) f ν ∈ C 1 ;

(3) f j ∈ C ζj for all j from 1 until ν;

(4) α ∈ {0, 1 2 } and β ∈ {1, 1 2 }. Proof. For simplicity, write ∆ = {0, 1 2 }×{1, 1 2 }. By (4.11), it follows that (α, β) ∈ ∆ if and only if (α j , β j ) ∈ ∆ for (one of) all indices j from 0 to ν. In addition, applying Theorem 3.1 to f ν implies that (α ν , β ν ) ∈ ∆ if and only if f ν ∈ C 1 . Thus, by considering Lemma 4.1, one finds that all conditions (1) -( 4) stated in Theorem 4.4 are equivalent. Now, we are ready to finish, successively, the proofs for Theorems 4.1 and 1.2 and, therefore, the proof for the main theorem.

Proof of Theorem 4.1. This follows directly from Theorem 4.4.

Proof of Theorem 1.2. In view of Theorems 3.1 and 4.1, it suffices to notice that, given ζ ∈ U, one has (x 0 q β ; q) ∞ ∈ C ζ if and only if the same holds by replacing β with β + 1. This last equivalence can be deduced from the relation (x 0 q β ; q) ∞ = (1 -x 0 q β ) (x 0 q β+1 ; q) ∞ and the fact that (1 -x 0 q β ) ∈ C ζ , for C ζ \ {0} constitutes a multiplicative group.

Figure 2 .

 2 Figure 2. τ belongs to the common domain V d1,+ ∩ V d2,-while the directions d 1 and d 2 are separated by the half-line τ

Remark 4 . 2 .

 42 Given r > 0, α ∈ (-1, 0) and β ∈ R, there exists ρ ∈ (0, r) such that the function τ → b + (α + β(τ -r), τ ) is well-defined and analytic in the open disc {τ ∈ C : |τ -r| < ρ}.

Remark 4 . 3 .

 43 If n and d are co-prime positive integers such that r = n d ∈ (0, 1), then

  below, where W and W c are given in (4.1).

Figure 5 .Remark 4 . 4 .

 544 Figure 5. P (z, τ ) is analytic for z ∈ (-W c -1) ∩ (W c + 1) and τ ∈ W

Lemma 4 . 1 .

 41 Let r ∈ (0, 1) ∩ Q, ζ = e(r), ζ 1 = e({-1 r }), and let z(τ ) = α + β τ with (α, β) ∈ [0, 1) × (0, 1]. Consider

α 1 , β 1 )

 11 being the transform of (α, β) defined by (4.10). Then f ∈ C ζ if and only if g ∈ C ζ1 . Proof. For simplicity, write r 1 = {-1 r } and N = [-1 r ]. As -1 r = r 1 + N , the following equivalence holds on the upper half-planes τ > 0 and τ 1 > 0

4 . 3 .

 43 Furthermore, replacing u with τ s in (3.6) of Lemma 3.1 gives that Li 2 (e(τ se 2πi )) -Li 2 (e(τ s)) = 4π 2 τ s, so one gets the identitys log(se 2πi ) + Li 2 (e(τ se 2πi )) 2πiτ = s log s + Li 2 (e(τ s)) 2πiτ , which implies that H 2 ∈ C{s} = C{τ }. Finally, one finds that H ∈ (C{τ } \ {0}) ⊂ C ζ ,which ends the proof. Continued fractions and modular transforms. Let us consider the asymptotic behavior of the Euler function (z | τ ) ∞ when τ a.v.

  Functions related with Bernoulli integrals. First, let P (z, τ ) be as in (2.1) and consider how to express it by using b + (z, τ ). In view of the fact that

	Proof. Since (z -1, τ ) ∈ W + , relation (2.17) implies that either z ∈ H ∪ (0, ∞)
	or z ∈ H -but (z/τ ) < 0. This implies that s / ∈ (-∞, 0]. Also, relation (2.18)
	follows immediately, by comparing (2.7) with (2.6).
	2.2. cot	t 2	-	2 t	=	e it/2 + e -it/2 e it/2 -e -it/2 i -	2 t	= 2i	1 e it -1	+	1 2	-	1 it	,
	it follows that cot	t 2	-	2 t	= 2i B

3. The half-plane Z τ contains both the point τ and the segment (-1, ∞) Lemma 2.2. Let b + (z, τ ) be as in Definition 2.1, and let s = z/τ . One supposes that τ ∈ H. If (z -1, τ ) ∈ W + , then s / ∈ (-∞, 0] and (2.18) b + (z -1, τ ) -b + (z, τ ) = I(s) , where I(s) is the Laplace integral stated by (2.6). (it); see (2.4) for B(t). Thus, replacing the integration path (0, ∞e id ) with (0, i∞e id ) in (2.2) yields that (2.19)

Proof of Theorem 1.1. By taking into account Remark 1.1 and Theorem 1.2, one needs only to observe that, for any positive integer n ∈ Z >0 and any root ζ ∈ U, any finite product of the form (x 0 q β ; q) n does not belong to T ζ , although the same function belongs to the larger class C ζ .

Addendum. After having finished a first version of our paper, we learned that the interesting work [START_REF] Katsurada | Asymptotic expansions of certain q-series and a formula of Ramanujan for specific values of the Riemann zeta function[END_REF] is closely related to the present paper. Indeed, let α > 0 and µ ∈ [0, 1) be as in [4, Theorem 1]. By combining [4, (3.2) & (3.3)] together with (1.2) and (1.4), one can observe the following result: Remark 4.5. One has (e(µ)q α ; q) ∞ ∈ T 1 only if the following conditions are satisfied for all integers k ≥ 2:

In the above, B k (α) denotes the usual Bernoulli polynomials, and B k (α, y) are the rational functions defined by the Taylor series expansion

Proposition 4.1. The following condtions are equivalent for any (µ, α) ∈ [0, 1) × (0, +∞):

(1) relation (4.20) holds simultaneously for all integers k ≥ 2;

(2) relation (4.20) holds simultaneously for k = 2 and k = 3;

(3) (µ, α) ∈ {(0, 1 2 ), (0, 1), (

2 )(α -1). Thus, (4.20) holds simultaneously for k = 2 and k = 3 if and only if α ∈ { 1 2 , 1}. Next, let µ ∈ (0, 1) and y = e(µ) = 1; putting α = 0 into (4.21), one can get that B 0 (0, y) = 0, B 1 (0, y) = 1 y-1 , B 2 (0, y) = -2 y (y-1) 2 and B 3 (0, y) = 3y(y+1) (y-1) 3 . Again, one can notice that B 2 (0, y) = 0, so relation (4.20) holds for k = 2 only if α ∈ { 1 2 , 1}. Since B 4 = -1 30 , one obtains, applying (4.22), that B 4 (α) = α 4 -2α 3 + α 2 -1 30 . One can see that B 4 ( 1 2 ) = 0 and B 4 (1) = 0. Therefore, (4.20) holds simultaneously for both k = 2 and k = 3 if and only if α ∈ { 1 2 , 1} and e(µ) = -1. This implies the equivalence between the conditions (2) and (3) stated in Proposition 4.1. Finally, suppose that condition (3) holds. It follows from Theorem 1.1 that (e(µ) q α ; q) ∞ ∈ T 1 , which, together with Remark 4.5, implies [START_REF] Andrews | Special functions, Encyclopedia of Mathematics and its Applications[END_REF]. In this way, we complete the proof of Proposition 4.1.


Laurent Debreu
email: laurent.debreu@inria.fr

Eric Blayo

Two-way embedding algorithms: a review submitted to Ocean Dynamics: Special Issue on Multi-Scale Modelling: Nested Grid and Unstructured Mesh Approaches

Keywords: Two way embedding, Mesh refinement, Structured grids

Local mesh refinement features have now been added to a number of numerical ocean models. In its crudest form, a high resolution grid is embedded (or nested) in a coarse resolution grid, which covers the entire domain, and the two grids interact. The aim of this paper is to review existing two way grid embedding algorithms. The basic algorithms and specificities related to ocean modelling are first described. Then we address several important issues: conservation properties, design of interpolation/restriction operators and noise control techniques.

Introduction

An increase of the horizontal resolution of a numerical ocean model still remains a key point in the improvement of the realism of its solutions, mainly through a better representation of small scales and domain geometry. This increase of resolution is generally not performed everywhere in the domain, both because it is not necessary from a physical point of view and because of limited computational resources. That is why a number of ocean models include local mesh refinement features. The idea is to refine the mesh where (and potentially when) necessary according to the objective of the simulation. For models based on a structured grid, a possible way to locally increase the resolution is to use a grid with variable resolution. Starting from a uniform grid, a mapping is introduced to produce a grid with increased resolution in areas of interest [START_REF] Zhuo | The Application of Time-Dependent Adaptive Mesh Model in Typhoon Track Prediction[END_REF]. This approach has the advantage that the model is still written on a single grid and thus there is almost no additional coding complexity and properties like conservation are easily handled. The disadvantages lie mainly in the difficulty of generating grids with good numerical properties, and in the fact that the time step constraint relies on the smallest grid cell (in the case of explicit time stepping schemes). The second approach, which is the subject of the present paper, is to embed a high resolution (HR) version of the same model at a specific location. The information is exchanged in two ways between the coarse resolution (CR) solution and the HR solution: the coarse grid provides the boundary conditions of the child grid, while the high resolution solution is used to update the coarse grid solution in the common domain. The major interest of this approach is that the coarse grid provides accurate boundary conditions for incoming information along the HR grid boundaries while the update (or feedback) step enables outgoing information to leave the fine grid much better than in one way interaction (i.e. without feedback from the HR to the CR solution). In this approach, even if the same model is applied, numerical schemes and/or physical parameterizations can be adapted to the grid resolution. We will discuss here only refinement of the mesh along the horizontal directions. A few applications of vertical mesh refinement have also been reported. In [START_REF] Fox | Two-way interactive nesting of primitive equation ocean models with topography[END_REF], the vertical refinement is applied to the whole water column. For similar reasons to the non decomposition of the vertical direction in parallel applications (use of implicit schemes and particular treatment of the barotropic mode), local mesh refinement on the vertical introduces additional difficulties, some of them very close to the ones reported in this paper for horizontal refinement. Concerning horizontal refinement, we will focus here on several important aspects including conservation, intergrid transfer operators and noise control techniques. We begin in section 2 with a description of the basic algorithms underlying the two way mesh refinement method and then give an overview of difficulties that can arise in practical numerical simulations. Conservation issues are addressed in section 3. In section 4, intergrid transfer operators (interpolations and updates) are examined and important details of implementations are given. Noise control techniques, which aim at making the algorithm more robust, are presented in section 5.

Two-Way nesting: Algorithms

In this paper, we will consider only one high resolution grid embedded in the coarse resolution grid. Extension to more than one grid and/or more than one level of refinement does not present additional complexities. So let's consider Ω the domain covered by the coarse grid and a subdomain ω, covered by the fine grid. γ is the boundary of ω. At the discrete level, the numerical grids are denoted by Ω H and ω h (cf. Fig. 1). ω H is the part of the grid Ω H corresponding to the domain ω.

Basic Algorithm

The ratio between the coarse and fine horizontal mesh sizes is an integer ρ, the mesh refinement factor (ρ = 2 on Fig. 1). Typical values of ρ in actual applications range from 2 to 5. Using larger values introduces a too strong change in the resolved scales on the different grids. In that case, it may not be appropriate to try to enforce a strong coupling between the two grids and methods with weaker interactions than the ones we will discuss in this paper could be of interest [START_REF] Sheng | A new two-way nesting technique for ocean modeling based on the smoothed semi-prognostic method[END_REF].

Ω H ω h γ h
Fig. 1 Local mesh refinement of a structured grid with a mesh refinement factor of 2. The high resolution grid ω h is embedded in the coarse resolution grid Ω H .

Choosing an odd mesh refinement factor simplifies grid interactions since in that case a coarse grid point has always one underlying high resolution point. Figure 2 shows, for B and C grids (using the Arakawa classification), that if a variable is staggered (cell or face centered) and if the mesh refinement factor is even, then a coarse grid point does not have a corresponding point on the fine grid. Associated with this spatial mesh refinement, time refinement may also be applied. Most of the ocean models now use explicit time integration algorithms both for computational reasons (efficiency on distributed memory parallel computers) and for ac-curacy reasons (smaller dispersion errors than implicit schemes). The model is then subject to a given CFL (Courant Friedrichs Lewy) stability condition and the ratio Δt/Δx must be kept smaller than a given value on the whole grid hierarchy. As we shall see, refining in time induces several difficulties but is required for some applications, especially when only a small portion of the coarse grid domain is refined or when there are several levels of refinement. In these cases, keeping the same time step over the whole grid hierarchy greatly increases the computational cost. The integration algorithm for a time refinement factor ρ t of 3 is depicted on Fig. 3.

1 2 ω h Ω H 3 4 ∆t c ∆t f ∆t f ∆t f Interpolation Update

Model Integration

Fig. 3 Integration algorithm for a time refinement factor of 3

The model is first integrated on the coarse grid with a time step equal to Δt c and then the fine grid is integrated ρ t times with Δt f = Δt c /ρ t . Interpolations of coarse grid boundary data occur at the end of each fine grid time step and updates occur at the end of coarse time step, when the two solutions have been advanced to the same physical time.

Let L represent the model integration from time t n to time t n+1 and L c , L f its discretizations on the CR and HR grids. Let P denote the interpolation (or prolongation) operator from Ω H to γ h , the boundary of ω h , and R the update (or restriction) operator from ω h to ω H . Then, assuming that the model is fully explicit, the algorithm can be written in the following simplified form :

1. u n+1 c = L c (u n c ) 2. For m =1...ρ t do u n+ m ρ t f = L f u n+ (m-1) ρ t f u n+ m ρ t f | γ h = P (u n c ,u n+1 c ) 3. u n+1 c | ω H = R(u n+1 f )
P , the interpolation operator, includes a time interpolation. This time interpolation is done linearly in most of the applications and thus requires only the knowledge of u c at time t n and t n+1 . To our knowledge, the effect of increasing the interpolation order in time has not been studied in the literature. This is because, as for numerical models on a uniform grid, the error relative to time is expected to be lower than the one relative to spatial interpolation. However, as for numerical models, it may not be the case anymore due to the increase of horizontal resolution and to the use of higher order spatial schemes. New experiments are probably required to test the validity of the linear interpolation. Note finally that interpolation and restrictions can also be applied to time tendencies instead of instantaneous values [START_REF] Jones | A nested grid for a three-dimensional model of a tropical cyclone[END_REF][START_REF] Spall | A nested primitive equation model for oceanic applications[END_REF][START_REF] Oey | A nested-grid ocean model: with application to the simulation of meanders and eddies in the Norwegian Coastal Current[END_REF].

Complexities associated with time stepping

The numerical schemes are seldom as simple as previously described and this can lead to several complications. In a numerical ocean model, some of them are linked to the particular treatment of the gravity waves associated with the evolution of the barotropic mode. This treatment leads either to the use an implicit time stepping, thus leading to the solution of an elliptic system, or to the use of a time splitting method.

Elliptic systems

Elliptic systems can come from a 2D implicit solve (e.g. implicit treatment of free surface [START_REF] Dukowicz | Implicit free-surface method for the Bryan-Cox-Semtner ocean model[END_REF] or rigid lid approximation). A 3D system has also to be solved if the hydrostatic assumption is removed. If we denote the associated linear system by :

Av = B
then after discretization, on the nested grid, the simplest choice is to force the high resolution solution with an interpolation of the coarse grid solution :

a) A c v c = B c on Ω H ,b ) A f v f = B f on ω h v f | γ h = Pv c (1)
This implies continuity of the solution at the interface γ h but the resulting gradient across the boundary is discontinuous. Additionally a loss of accuracy is introduced because errors produced by the coarse grid resolution a) propagate inside the HR domain through the boundary forcing in step b). A first step toward better accuracy is to update the right hand side of the coarse grid equation so that Eq. 1 is replaced by Eq. 2:

a) A c v c = RB f in ω H B c in Ω H \ω H ,b ) A f v f = B f on ω h v f | γ h = Pv c (2)
In [START_REF] Spall | A nested primitive equation model for oceanic applications[END_REF], this technique is applied to a rigid lid ocean model. The linear system corresponds to the computation of the streamfunction tendency and the right hand side is the barotropic vorticity tendency. The barotropic vorticity tendency on the coarse grid is updated by the one computed on the fine grid. However, despite this improvement, the gradient of the solution across the interface is still discontinuous. To be properly solved the two grid system should be handled as a truly local multilevel system. This can be done using local defect correction methods (see [START_REF] Laugier | Nested grid methods for an ocean model: a comparative study[END_REF] for an application to a rigid lid model and a comparison with results obtained in [START_REF] Spall | A nested primitive equation model for oceanic applications[END_REF]. The system can also be seen as discretized on a single composite grid. In this approach, the operator A is modified at the boundary and involves both coarse and fine grids variables. The resulting system is:

⎛ ⎝ A c 0 0 A f A cγ A fγ ⎞ ⎠ v c | Ω H \ω H v f = ⎧ ⎨ ⎩ B c in Ω H \ω H B f on ω h B γ in γ h (3) A cγ v c | Ω H \ω H +A fγ v f = B
γ is an approximation of the original system on the interface which leads at convergence to continuity of both the solution and its gradient across the boundary. See [START_REF] Martin | Solving Poisson's equation using adaptive mesh refinement[END_REF] for applications of the local multilevel method to the solution of a Poisson equation.

Remarks:

-Algorithms 2 and 3 are more difficult to apply when time refinement is used. Indeed the right hand side B c and B f should be computed at the same time and the way to write the system at intermediate time steps on the fine grids is not obvious. -Since we are using two way techniques, the right hand side B c , which depends on coarse grid variables that have been updated, is never significantly different from RB f . In some cases, the update schemes used for the primary variables are precisely chosen in such a way that B c = RB f in ω H . However, this leads most of the time to the use of very simple restriction operators (e.g. average, see [START_REF] Clark | Severe downslope windstorm calculations in two and three spatial dimensions using anelastic interactive grid nesting: A possible mechanism for gustiness[END_REF] that do not have good filtering properties and thus can affect the quality of the solution (see paragraph 4.1).

Free surface and Time splitting

Similar problems are also inherent when a time-splitting is used for the treatment of the external mode. A number of models now use the time-splitting approach of [START_REF] Blumberg | A description of a three-dimensional coastal ocean circulation model[END_REF]; [START_REF] Killworth | The development of a free-surface Bryan-Cox-Semtner ocean model[END_REF]. The barotropic time step is set to a fraction of the baroclinic time step. The barotropic mode is then integrated as a solution of a shallow water model forced by the vertical integral of the right hand side of the 3D equations. At the end of this integration, the newly computed barotropic quantities are used to replace the barotropic part of the 3D fields. In order to prevent aliasing errors, two approaches can be taken. The first one, as in [START_REF] Killworth | The development of a free-surface Bryan-Cox-Semtner ocean model[END_REF], is to use in the barotropic time steps a time integration scheme that sufficiently damps high temporal frequencies (e.g. Euler Backward). In this case, the fine and coarse grid can interact during the small barotropic time steps throughout the barotropic quantities (see [START_REF] Oey | A nested-grid ocean model: with application to the simulation of meanders and eddies in the Norwegian Coastal Current[END_REF]. This is illustrated on Fig. 4 with a time refinement factor of 2.

The second approach consists in applying a filter to the barotropic quantities to filter out the high temporal frequencies. The integration extends to a time larger than t+Δt c and a filtering formula, the result of which is centered at time t + Δt c , is applied. For example, as shown on Fig. 5, one may extend the integration period to t +2Δt c and use flat weights over [t, : t +2Δt c ] (constant weights equals to 1/N where N is the number of barotropic time steps) to compute the average at time t + Δt c . When filtering is applied, it is relatively easy to interact between grids when no time refinement is done as on Fig. 6 [START_REF] Barth | Two-way nested model of mesoscale circulation features in the Ligurian Sea[END_REF]. However, when both filtering and time refinement are used (cf. Fig. 7) a new problem arises, due to the fact that some coarse grid points do not have fine grid points located at the same time. For this reason, the coupling between coarse and fine grids cannot be done anymore through the instantaneous barotropic quantities. In that case, a simple choice is to perform the coupling at the baroclinic level [START_REF] Penven | Evaluation and application of the ROMS 1-way embedding procedure to the central california upwelling system[END_REF], i.e. the exchange occurs only between filtered quantities. However this is like having the system in a one-way mode for the external mode and, as for implicit solvers, errors produced on the coarse grid can then propagate inside the HR grid. One way to perform the coupling at the barotropic level is to exchange information between another set of variables corresponding to intermediate quantities (Debreu et al. 2008). Using intermediate weights, quantities centered at times between t and t + Δt c are computed so that coupling as without filtering (Fig. 4) is permitted.

t t+∆t c t t+∆t c /2 t+∆t c G c G f
t G c G f t+∆t c t+2∆t c t t+∆t c t+∆t c /2 t+3∆t c /2 ?
Partial conclusion As we have seen, besides the simple algorithm, several additional complexities arise when implicit schemes or time splitting schemes are used. Moreover, it is clear that, in both cases, use of time refinement leads to additional difficulties. Looking at the problem as a truly multilevel system in the case of elliptic systems or doing the coupling at the barotropic level in the case of a time splitting scheme are probably the best choices. Additionally, if more simple techniques are used, then it brings strong constraints on the update operators in order to maintain the coarse grid right hand side in agreement with the fine grid right hand side. Most of the time, these constraints leads to the use of very simple update operators that do not have good numerical properties (see 4.1).

Conservation

On a uniform grid, conservation is guaranteed when internal numerical schemes are written in flux form. Concerning two-way embedded models, it is not so common that the resulting system is conservative. Ensuring conservation leads to several computational issues and imposes strong requirements on intergrid transfer operators that can lead to a loss of accuracy. However it is recommended for long term integration. After reviewing the basic requirements for conservation on an embedded grid, two familiar approaches for the preservation of conservation properties are described.

Definition and discretization

Let us consider a two dimensional domain and q, a solution of the following equation written in conservative form ∂q ∂t

+ ∂f ∂x + ∂g ∂y =0
where f and g may contain both advective and diffusive fluxes. Then, assuming that fluxes f and g cancel on the boundaries of Ω, Q Ω , integral of q over the domain Ω, is constant in time (at the continuous level).

Q Ω (t)= Ω q(x, y, t)dx dy ⇒ dQ Ω (t) dt = ∂Ω f ds + ∂Ω g ds =0
In the nested grid system, the quantity Q Ω is defined by the summation over the high resolution domain ω and its complement in Ω:

Q Ω = Q ω + Q Ω\ω
Now let us make the following assumptions for the sake of simplicity:

-As shown on Fig. 8, we consider a two dimensional domain infinite in both x and y directions. The left (resp. right) part of the domain is at coarse (resp. high) resolution. -The variable q n i,j is cell centered

Ω H ω h i c i f j j 1 j 1 +1 j 1 +2
Fig. 8 The CR domain Ω H on the left and the HR domain ω h on the right for a C-grid with a mesh refinement factor of 3 -The time stepping scheme is an explicit Euler scheme

q n+1 i,j = q n i,j - Δt ΔxΔy F i,j -F i-1,j - Δt ΔxΔy G i,j -G i,j-1 (4)
where F i,j ,G i,j are volumetric fluxes: F i,j = f i,j Δy, G i,j = g i,j Δx. At discrete level, Q Ω at time t = t n is given by

Q n Ω = i≤i c ,j Δx c Δy c q c,n i,j + i≥i f ,j Δx f Δy f q f,n i,j
where, as shown on Fig. 8, i c and i f denote the first coarse and fine indices close to the interface. Then after one time step, according to Eq. 4 we obtain

Q n+1 Ω = Q n Ω -Δt c j F c,n i c ,j + Δt f ρ t -1 p=0 j 1 +ρ-1 j f =j 1 F n+p/ρ t i f -1,j f = Q n Ω -Δt c j ⎛ ⎝ F n i c ,j - 1 ρ t ρ t -1 p=0 j 1 +ρ-1 j f =j 1 F n+p/ρ t i f -1,j f ⎞ ⎠ (5)
In general, there will be a misfit between the coarse and fine grid fluxes so that the flux differences of the right hand side do not cancel and conservation is artificially lost:

Q n+1 Ω = Q n Ω .
Conservation can be achieved if the fine grid fluxes derive from a conservative interpolation (in time and space) of the coarse grid flux so that the following relation holds:

1 ρ t ρ t -1 p=0 j 1 +ρ/2 j f =j 1 -ρ/2 F n+p/ρ t i f -1,j f = F n i c ,j (6)
This approach is applied in several ocean models for volume conservation (q being the free surface and F, G the vertically integrated transport in x and y directions). It is indeed sufficient to produce a conservative interpolation of the transport for the global volume to be conserved. A number of classical conservative interpolation schemes can then be used (see 4.2). Let's consider more generally the computation of a boundary flux on a C-grid as indicated on Fig. 9 and let's suppose than the we are trying to compute the flux of a tracer on the boundary so that: The computation of fine grid flux on the boundary involves the result of the interpolation of the velocity u on the boundary (grey ellipses on Fig. 9). This interpolation has to be conservative in order not to create an artificial loss of conservation. The fine grid flux also involves the computation of the tracer's value on the interface that makes use of both interpolated values (grey circles on Fig. 9) and internal values in the fine grid domain. Thus, we may suppose that the resulting fine grid fluxes are actually more accurate than the coarse grid ones, so that trying to enforce Eq. 6 is not the best choice.

f = uT Ω H ω h i c i f j j 1 j 1 +1 j 1 +2
Two other approaches to conservation, that share several similarities, are now introduced: the flux correction algorithm and the Kurihara method.

Flux correction algorithm

The flux correction algorithm comes from the adaptive mesh refinement community. It follows from the algorithm of [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF]; [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF].

The idea is to apply a modification of the coarse grid variables that take into account the misfit between coarse and fine grid fluxes. Starting from Eq. 5, a correction is applied to the coarse grid variable at time n +1 near the boundary as follows :

q n+1, * i c ,j = q n+1 i c ,j + Δt c Δx c Δy c ⎛ ⎝ F n i c ,j - 1 ρ t ρ t -1 p=0 j+ρ/2 j f =j-ρ/2 F n+p/ρ t i f -1,j f ⎞ ⎠ (7)
Basically it means that the coarse grid variable has been integrated using, on the right interface, fluxes computed by the fine grid :

q n+1,⋆ i c ,j = q n i c ,j - Δt c Δx c Δy c ⎛ ⎝ 1 ρ t ρ t -1 p=0 j+ρ/2 j f =j-ρ/2 F n+p/ρ t i f -1,j f -F i c -1,j ⎞ ⎠ (8)
The algorithm is easy to implement assuming that it is possible to write the time evolution of q in term of flux divergences as it was indeed the case for the Euler time scheme. The program stores the fine grid fluxes at the boundary and makes a summation in time and in space over the fine grid cells. With other time stepping schemes, e.g. leap-frog associated to an Asselin filter, this is not achievable (even without the Asselin filter, it is more complex since two arrays of fluxes should be stored, see [START_REF] Herrnstein | Structured adaptive mesh refinement using leapfrog time integration on a staggered grid for ocean models[END_REF].

When using this approach, care must be taken concerning potential stability issues. Indeed, if the original fluxes are computed in a centered way, the resulting scheme is always biased and instabilities can occur [START_REF] Olsson | Stability of interpolation on overlapping grids[END_REF]Debreu et al. 2008). One remedy is to compute the interfacial fluxes using information along characteristics of the flow (Part-Enander and Sjogreen 1994, and Blayo and Debreu 2005 for an interpretation of usual open boundary conditions in term of characteristic variables). This is actually the case in the adaptive mesh refinement community where, due to the hyperbolic nature of the problems, the internal numerical schemes already make use of decompositions into characteristic variables. Moreover, there is an associated loss of accuracy at the boundary. Indeed it can be shown that if the fluxes are computed with a second order accurate approximation, use of Eq. 8 leads to a first order only approximation of the original equation (Debreu et al. 2008). These two remarks also apply to the Kurihara method that we now describe.

Kurihara method

The Kurihara method [START_REF] Kurihara | Design of a movable nested-mesh primitive equation model[END_REF]) has been used by several authors to enforce conservation [START_REF] Kurihara | Use of a movable nested-mesh model for tracking a small vortex[END_REF][START_REF] Ginis | Design of a multiply nested primitive equation ocean model[END_REF]. It is important to note that it applies to a model with all variables cell centered. To understand the method, it is first necessary to explain that the scheme makes use of an intermediate area, composed of two coarse grid cells (grey circles on Fig. 10), where the space resolution is the one of the coarse grid and the time resolution is the one of the fine grid. In this method, conservation is achieved thanks to two principles :

-The coarse grid computes its fluxes which are then interpolated in time in a conservative way at ✷ points to force the separation area.

Dynamic Interface Mesh Interface

t n t n +∆t f t n +2∆t f t n +∆t c
Fig. 10 The three zones of the Kurihara method: on the left the coarse grid, on the right the fine grid. In between, an area composed of two coarse grid cells with the space resolution of the coarse grid and the time resolution of the fine grid.

-A spatial linear interpolation of variables is performed at △ points and the computed flux is used both for integrating the coarse grid cell on the left and the fine grid cell on the right. In two dimensions, a summation is applied along the y axis just as for the flux correction approach.

Thus the resulting scheme is conservative. The problems of loss of conservation due to time refinement and to space refinement has been treated separately.

It is clear that if there is no time refinement, then the flux correction approach is equivalent to the Kurihara method if the boundary variables used in the flux computation have also been linearly interpolated.

As already said, the original method is designed for a non staggered grid with all variables located at the cells center. Like for the flux correction approach, the Kurihara method is more difficult to implement for non cell-centered variables (for example in order to preserve conservation of momentum on a C-grid, cf. [START_REF] Sobel | Nested grids in numerical weather prediction and an application to a mesoscale streak[END_REF].

Partial conclusion

We have presented two different approaches for enforcing conservation. Due to the presence of the transition area, the Kurihara method is more difficult to implement than the flux correction approach. Additionally, although possible, it is more difficult to implement such methods for quantities not defined at the centers of the cells. The flux correction approach is easier to implement; however care about potential stability issues must be taken, especially when the original model computes its fluxes using second order centered differences. When conservation properties are not preserved in a two way simulation, it is important, especially for long term integration, to have an idea of the artificial loss of conservation. That is why the misfit between coarse and fine grid fluxes on the interface should always be diagnosed.

Intergrid Transfer operators

For very simple equations (e.g. 1D advection-diffusion), the properties of the two grids methods can be analyzed with simple techniques related to study of numerical schemes on non uniform grids (Harrisson and Elsberry 1972;[START_REF] Berger | Stability of Interfaces with Mesh Refinement[END_REF][START_REF] Olsson | Stability of interpolation on overlapping grids[END_REF]. In the linear case, it can even be extended to more complex systems of equations using matrix stability analysis [START_REF] Heggelund | A method for analysing nesting techniques for the linearized shallow water equations[END_REF]. However in the more general case, the design of intergrid transfer operators must rely on more crude considerations.

Update schemes

Obviously, the restriction operator has a crucial role in two way nesting algorithms. If we just think in terms of resolved scales, the two following properties should hold :

1. The transfer of information should be maximum for scales well resolved on the coarse grid.

Otherwise the order of the approximation will be lowered. 2. Small scales should be strongly filtered.

Otherwise, by aliasing, noise will be produced on the coarse grid.

In several applications, authors have made the choice to use the average operator for the feedback step : the coarse grid value in a cell is replaced by the area weighted sum of the fine grid values in the same cell :

u c = 1 Δx c Δy c i,j Δx f Δy f u f
Note that this is often seen as a requirement for conservation but as it has been previously seen, it is not. In an adaptive mesh refinement context, the restriction is still made with this average formula because the grid can move from one time step to another, so that an average restriction is necessary to maintain globally the conservation before and after the regridding step (periodic modification of the grid hierarchy following a refinement criterion). Note that this was also the main reason for using the average restriction operator in the Kurihara paper [START_REF] Kurihara | Design of a movable nested-mesh primitive equation model[END_REF], where the methodology is applied to moving mesh methods. Other familiar restriction schemes are written on Fig. 11 along with their corresponding transfer functions. Recently in Debreu et al. (2008), it is shown that use of higher order restriction operators, using a larger stencil (full weighted schemes), can lead to strong improvements. Even if the full weighted operator damps the well resolved scales a little more, it also produces a strong damping of the small scales which is required to inhibit aliasing and noise on the coarse grid.

Interpolation

Numerous interpolation operators have been used in two way nesting methods. See [START_REF] Koch | A survey of Nested Grid Techniques and Their Potential for Use Within the MASS Weather Prediction Model[END_REF]; [START_REF] Zhang | A Two-Way Interactive Nesting Procedure with Variable Terrain Resolution[END_REF] for the details in several models. The Copy

u c = u f Average u c = 1 3 u f j-1 + u f j + u f j+1 Shapiro u c = 1 4 u f j-1 +2u f j + u f j+1 Full Weighting u c = 1 9 u f j-2 +2u f j-1 +3u f j +2u f j+1 + u f j+2
Fig. 11 Transfer functions for usual one dimensional restriction operators with a mesh refinement factor ρ = 3 for typical filters (k is the wavenumber). The corresponding formula are also given.

basic idea is that the order of the interpolation order should be high enough not to lower the global order of approximation and, as usual when using high order interpolation schemes, care must be taken about the potential oscillations produced by the schemes (this problem is well known, for example, when using Lagrange interpolation operators). As previously seen, it is also advantageous to have conservative interpolations (in particular for the interpolation of normal velocities on the boundaries). To fulfill these two conditions several methods can be used: parabolic conservative interpolation or approximation [START_REF] Clark | Severe downslope windstorm calculations in two and three spatial dimensions using anelastic interactive grid nesting: A possible mechanism for gustiness[END_REF], parabolic interpolation with minimization of second order derivative [START_REF] Barth | Two-way nested model of mesoscale circulation features in the Ligurian Sea[END_REF], PPM (Piecewise Parabolic Method), WENO (Weighted Essentially Non Oscillatory), advection equivalent interpolation schemes [START_REF] Alapaty | Intercomparaison of Spatial Interpolation Schemes for Use in Nested Grid Models[END_REF] . . . Several authors also mention the fact that interpolation and restriction operators must be devised in agreement one with the other. This is mainly justified when implicit solvers are used. The main reason is that, after feedback, the right hand side of the elliptic equation, computed on the coarse grid, must be in agreement with the one of the fine grid [START_REF] Clark | Severe downslope windstorm calculations in two and three spatial dimensions using anelastic interactive grid nesting: A possible mechanism for gustiness[END_REF]. In that case, assuming a conservative interpolation has been used, the average restriction operator is the most natural. Let's take the example of the free surface equation written under the following form (Eq. 9):

∂η ∂t + ∇•(U, V ) = 0 (9)
where η is the free surface and (U, V ) the barotropic transports in x and y directions. If the transport U has been conservatively interpolated on the fine grid boundary, then a simple update (average in one direction and simple copy in the other direction as shown on Fig. 12) of the transports leads to a barotropic divergence in the coarse grid cell equal to the average of the barotropic divergence on the fine grid cell:

∇•(U c ,V c )= 1 9 i,j ∇•(U f i ,V f j ) (10) Ω H ω h i c i f j j 1 j 1 +1 j 1 +2 V c = 1 3 V f i f +V f i f +1 +V f i f +2 U c = 1 3 U f j 1 +U f j 1 +1 +U f j 1 +2 1 3 U f j 1 +U f j 1 +1 +U f j 1 +2 = U c
Fig. 12 Update and interpolation schemes for the preservation of the barotropic divergence on a C-grid with a mesh refinement factor of 3. The transport (U ) are conservatively interpolated on the boundary (grey ellipses) and then simple update operators are used at coarse grid points (black ellipses).

However, as we already saw (paragraph 2.2.1), the constraint given by Eq. 10, which actually corresponds to an equality between tendencies of the free surface equation instead of instantaneous values, can be removed if either the right hand side of the coarse grid is directly updated by the right hand side of the fine grid or if a truly multilevel solver is used. Removing this constraint allows the use of more scale selective update operators than those used here.

Inconsistency at boundary

We now look at a particular problem that also induces differences between implementation of intergrid transfer operators. Let's rewrite the explicit algorithm given in paragraph 2.1

1. u n+1 c = L c (u n c ) 2. For m =1...ρ t do u n+ m ρ t f = L f u n+ (m-1) ρ t f u n+ m ρ t f | γ h = P (u n c ,u n+1 c ) 3. u n+1 c | ω H = R(u n+1 f (u n+1 c
))

where we insist on the fact that u n+1 f is indeed dependent on u n+1 c through the boundary interpolation so that an inconsistency (sometimes called overspecification) is inherent to the third step. Note that this does not apply to variables entirely lying on the interface, like normal velocities on a C-grid, because the coarse grid variables used in the interpolation (along the interface) are not updated in the feedback step. On a C-grid, this inconsistency occurs for tracers and tangential velocities. One possible way to remove this inconsistency is to modify the interpolation step 2 by doing intermediate updates so that step 2 is modified as follows

u n+ m ρ t f | γ h = P (u n c ,u n+1 c ) → u n+ m ρ t f | γ h = P (u n c , u n+1 c in Ω H \ω H R(u n+ m ρ t f | γ h )i n ω H )
eventually using a fixed point iteration if the operators P and/or R are non linear.

To our knowledge, this method has not been explored. In the literature, several other methods have been developed to try to remove this contradiction:

-Interpolation operators It is possible to write the interpolations so that no coarse grid point, that will be modified during the feedback step (i.e. inside the HR grid domain) are used in the interpolation process. In order to do this, a tangential interpolation is first applied, followed by a normal interpolation. The idea is also to use the HR solution as soon as available thus lowering interpolation errors (both in space and time).

Note that here too, care must be taken about stability issues. As an example, using a high order interpolation, which involves several fine grid points inside the fine grid domain, biases the interpolation schemes and can produce unstable solutions if the flow is actually entering the fine grid. So that in practice the normal interpolation should be computed using the direction of the flow. Again this is especially true when the model computes its fluxes using second order centered differences.

Interpolation

Interpolation Tangential Normal Fig. 13 Tangent interpolation followed by a normal interpolation. A set of fictitious points (squared boxes) are first obtained by interpolation along the tangential direction, this is followed by a normal interpolation that, in addition to these points, uses fine grid points inside the fine grid domain. Note that the coarse grid points inside the fine grid domain are not used in the interpolation process.

-Update operators -Separation of dynamic and feedback interface Another remedy to remove this inconsistency is to modify the feedback step. Several authors have proposed to separate the feedback interface from the dynamic interface (sometimes also called input interface) where the boundary values are interpolated [START_REF] Phillips | On the strategy of combining coarse and fine grids meshes in numerical weather prediction[END_REF][START_REF] Zhang | A Two-Way Interactive Nesting Procedure with Variable Terrain Resolution[END_REF][START_REF] Oey | A nested-grid ocean model: with application to the simulation of meanders and eddies in the Norwegian Coastal Current[END_REF].

Similarly, the Kurihara method can also be seen as introducing a separation between forcing and feedback interfaces. Figure 14 shows a separation of the dynamic and feedback interface by one coarse grid cell. This separation can also be composed of two coarse grid cells. In these cases, the contradiction has been removed since the coarse grid points used in the interpolations are not updated. Another reason for using a mesh separation between dynamic and feedback interfaces is that if noise is produced, it will be larger near the dynamic interface so that it is safer not to use the fine grid values near the dynamic interfaces. However, as mentioned by [START_REF] Spall | A nested primitive equation model for oceanic applications[END_REF], setting the feedback interface closer to the dynamic interface maximizes the information transferred to the coarse grid. In addition since there is no feedback in this area, the consistency between coarse and fine solutions will be less than with feedback and this inconsistency has to be lowered by use of noise control techniques that will indeed recombine and smooth the HR and CR solutions in the area between dynamic and feedback interfaces.

Partial conclusion

From the authors' experience, update operators have a greater impact on the quality of the solution than interpolation operators (assuming that these last ones are properly built, e.g. conservative when they have to be). The usual update operators (average, Shapiro filter . . . ) do not have sufficiently good numerical properties to lead to stable simulations. In particular, use of update operators that produce a stronger damping of small scales is required. The technique of separation of dynamic and feedback interfaces can be used associated with noise control techniques to make robust two way embedding simulations. However errors introduced by these methods may be found to be very large when quantitatively evaluated.

Noise Control

Two principles can motivate the use of noise control techniques :

-Maintain a strong consistency between HR and CR solutions in the area where solutions interact (i.e. near the common interface). -Prevent waves reflection: waves unresolved on the coarse grid reflect inside the fine grid domain and have to be filtered in some way.

It is generally accepted that one criterion for choosing the noise control technique is that the associated modifications should cancel if the refinement ratio is equal to 1, so that if the simulation is made with a mesh refinement factor of 1, the result is identical to the result computed with a uniform grid.

Of course, use of noise control techniques modify the original differential operator near the interface and thus lower the order of approximation. Two-way embedding models can make use of noise control techniques that are traditionally used in open boundary and/or one way embedded simulations. We briefly recall here some of them:

-Use of a time stepping algorithm or other internal numerical schemes that damps small wavelengths. This is the crudest form of noise control technique: the original model solution is sufficiently damped so that small scales are strongly removed. -Relaxation methods -Blending (Flow Relaxation Scheme)

q ⋆ f =(1-μ x,∂ω )q f + μ x,∂ω Pq c μ x,
∂ω is a coefficient varying from 1 on the interface ∂ω t o0a w a yt ot h e interface [START_REF] Davies | A lateral boundary formulation for multi-level prediction models[END_REF]). This is used in [START_REF] Oey | A nested-grid ocean model: with application to the simulation of meanders and eddies in the Norwegian Coastal Current[END_REF] in the separation area between feedback and dynamic interfaces.

-Nudging ∂q f ∂t = ...- (q f -Pq c ) T d
The discrete form of nudging is equivalent to the flow relaxation scheme (Martinsen and Engedahl 1987). -Increase of dissipation coefficient -Sponge Layer When using a sponge layer it is advantageous to apply it to the difference between external (coarse resolution) and internal (high resolution) fields ∂q f ∂t = ...+(-1) n+1 (Δ) n μ x,∂ω (q f -Pq c ) (11)

Here the objective is clearly to damp small scales as defined by the restriction operator. Since the coarse grid variable q c has been updated (q c = Rq f ), the diffusion term mostly acts on the scales lying in the kernel of the restriction operator (Rq f = 0). Note that the preceding relaxation methods can be put under this category just by choosing n = 0 in Eq. 11 so that they produce a less selective damping than for n ≥ 1. The advantage of this formulation is also that the diffusion term naturally cancels if there is no refinement. Using differences between coarse and fine grid variables can also be at the basis of the derivation of open boundary schemes suitable for embedding methods (e.g. Perkins et al. 1997, Oddo and[START_REF] Oddo | Lateral open boundary conditions for nested limited area models: A scale selective approach[END_REF].

Additionally, in a flux correction approach of conservation, it is easy to maintain conservation just by adding the diffusive flux of Eq. 11.

Moreover, due to the inability to build ideal restriction operators, it may also be useful to smooth the coarse grid solution inside the fine grid domain. The filtering is only applied to the coarse grid points inside the fine grid domain in order not to violate conservation principles. As an example, the smoothing-desmoothing algorithm of [START_REF] Bender | Improvements in tropical cyclone track and intensity forecasts using the GFDL initialization system[END_REF] can be applied on the first coarse grid points near the interfaces.

Partial conclusion Two way grid embedding algorithms should be first evaluated without using noise control techniques in order not to hide fundamental problems of the algorithm. Because the primary requirement is that the associated modifications of the original model should cancel if the refinement factor is equal to 1, most of the methods are formulated in term of differences between coarse and fine fields as in Eq. 11.

Choosing n = 1 in Eq. 11 (laplacian diffusion) is a good compromise between a low computational cost and a selective damping of the difference. For n ≥ 1, the correction can be easily written in term of flux divergences thus allowing the conservation properties to be maintained (e.g. using the flux correction approach).

Conclusion

This paper reviews existing two way embedding techniques. The emphasis has been put on what are the sources of potential approximations and errors. In particular, we have shown how the use of implicit solver and/or time splitting algorithms and the use of refinement in time can introduce several problems in comparison with an idealized model. Conservation issues can be tackled in a simple way, using the flux correction approach (at least for quantities based on cell centered variables: mass, tracers . . . ).

There is also a large choice of intergrid transfer operators, the choice between one of them is a trade-off between conservation, accuracy, robustness and coding complexities. Some guidelines have been given as partial conclusions of each section. Several issues related to grid embedding have not been covered. Computational implementation of embedding can be done by hand or by using dedicated software (e.g. AGRIF Debreu et al. 2008, RSL Michalakes 1998). In the context of operational oceanography and production of high resolution forecasts at particular areas of interest, high resolution grid initialization schemes which are in agreement with the two way interactions have to be derived. As most of the operational systems are now using data assimilation methods, one possible way to look at the problem is to introduce the high resolution grid initialisation in the assimilation process and this is an active research subject.

Fig. 2

 2 Fig. 2 Refinement of a cell with even (ρ = 2) and odd (ρ = 3) refinement factors for a B and a C grid. Only one coarse cell is shown. When ρ = 3, a coarse grid point has always a corresponding fine grid point.

Fig. 4 Fig. 5 Fig. 6

 456 Fig. 4 Time splitting algorithm with a time refinement ratio equal to 2. Intermediate times (due to time refinement) on the fine grid are not shown. G c is the coarse grid and G f the fine grid.

Fig. 7

 7 Fig. 7 Time splitting algorithm with time refinement equal to 2 and with filtering

Fig. 9

 9 Fig. 9 Computation of a boundary flux on a C-grid

Fig. 14

 14 Fig.14Separation of dynamic and feedback interface on a C-grid for a mesh refinement factor of 3. Without a separation of these interfaces, the tracers and tangential velocities would also be updated in the light grey area.

Acknowledgements This work is a contribution to the MERSEA project. Partial support from the European Commission under contract SIP3-CT-2003-502885 is gratefully acknowledged. We would like to thank two anonymous reviewers for their comments and careful reading of the manuscript.