Juliette Chabassier

Marc Duruflé

Victor Péron

Equivalent boundary conditions for acoustic media with exponential densities. Application to the atmosphere in helioseismology

Keywords: Multiscale Expansions, Equivalent Boundary Conditions, Numerical Methods for Helioseismology Analyse Asymptotique, Conditions aux Limites Équivalentes, Méthodes Numériques pour l'Héliosismologie

We present equivalent boundary conditions and asymptotic models for the solution of a transmission problem set in a domain which represents the Sun and its atmosphere. This problem models the propagation of an acoustic wave in time-harmonic regime. The specific non-standard feature of this problem lies in the presence of a small parameter δ which represents the inverse rate of the exponential decay of the density in the atmosphere. The atmosphere can be approximated by local equivalent boundary conditions for small values of this parameter. We derive rigorously equivalent conditions up to the fourth order of approximation with respect to δ for the exact solution u. The construction of equivalent conditions is based on a multiscale expansion in power series of δ for u. Numerical results illustrate the good behavior of these equivalent boundary conditions for realistic values close to those observed for the Sun. Finally we measure the boundary layer phenomenon by introducing a "skin depth" function that turns out to depend on the mean curvature of the interface between the subdomains.

Introduction

The concept of Equivalent Boundary Conditions (also called approximate, effective, or impedance conditions) is rather classical in the modeling of wave propagation phenomena. Equivalent Conditions (ECs) are introduced to reduce the computational domain, for instance for scattering problems from highly absorbing obstacles [START_REF] Rytov | Calcul du skin effect par la méthode des perturbations[END_REF][START_REF] Leontovich | Approximate boundary conditions for the electromagnetic field on the surface of a good conductor[END_REF][START_REF] Lafitte | Équations de Maxwell et opérateur d'impédance sur le bord d'un obstacle convexe absorbant[END_REF][START_REF] Senior | Institution of Electrical Engineers[END_REF][START_REF] Ida | Impedance boundary conditions for transient scattering problems[END_REF][START_REF] Antoine | High-frequency asymptotic analysis of a dissipative transmission problem resulting in generalized impedance boundary conditions[END_REF][START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF][START_REF] Dauge | Comportement asymptotique à haute conductivité de l'épaisseur de peau en électromagnétisme[END_REF][START_REF] Caloz | On the influence of the geometry on skin effect in electromagnetism[END_REF][START_REF] Buret | Eddy currents and corner singularities[END_REF][START_REF] Dauge | Corner asymptotics of the magnetic potential in the eddy-current model[END_REF]. The main idea consists in replacing an "exact" model inside a part of the domain of interest by an approximate boundary condition. This idea is relevant when the effective condition can be readily handled for numerical computations, for instance when this condition is local [START_REF] Senior | Institution of Electrical Engineers[END_REF][START_REF] Bendali | The effect of a thin coating on the scattering of a time-harmonic wave for the Helmholtz equation[END_REF][START_REF] Haddar | Generalized impedance boundary conditions for scattering problems from strongly absorbing obstacles: the case of Maxwell's equations[END_REF].

The main application of this work concerns helioseismology. It requires a precise knowledge of the acoustic wave field both in the Sun and in the atmosphere of the Sun. We intend to work in the context of this application for which we consider that the region of interest consists of two subdomains : a first spherical region (the Sun) which is surrounded by a second region (the atmosphere) where the density is exponentially decreasing and behaves like exp (-α r) (here α is a large parameter and r represents the distance to the surface of the Sun). Then the rapid decay of the density inside the atmosphere limits the penetration of the acoustic wave inside a boundary layer which occurs inside the atmosphere in a vicinity of the Sun. The characteristic size of this boundary layer is proportional to the small parameter δ = 1/α. This boundary layer phenomenon is an obstacle that prevents from reaching the required accuracy for the acoustic field. It raises the difficulty of applying a discretization method (e.g. a finite element method FEM) on a mesh that combines fine cells in the second region and much larger cells in the first region. To overcome this difficulty and to solve this problem, we adopt an asymptotic method which consists in "approximating" the second region by an EC. It is then possible to solve a simpler problem resulting of the Helmholtz equation (set in the first region) with this new boundary condition. The boundary layer phenomenon which is described in this paper looks like the so-called skin effect phenomenon which occurs inside a highly absorbing obstacle in electromagnetism [START_REF] Rytov | Calcul du skin effect par la méthode des perturbations[END_REF][START_REF] Leontovich | Approximate boundary conditions for the electromagnetic field on the surface of a good conductor[END_REF][START_REF] Maccamy | A skin effect approximation for eddy current problems[END_REF][START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF][START_REF] Haddar | Generalized impedance boundary conditions for scattering problems from strongly absorbing obstacles: the case of Maxwell's equations[END_REF][START_REF] Péron | Modélisation mathématique de phénomènes électromagnétiques dans des matériaux à fort contraste[END_REF][START_REF] Dauge | Comportement asymptotique à haute conductivité de l'épaisseur de peau en électromagnétisme[END_REF][START_REF] Caloz | On the influence of the geometry on skin effect in electromagnetism[END_REF]. However these boundary layer phenomena are different.

In this paper we derive ECs up to the fourth order of approximation with respect to the parameter δ. These ECs are of "Neumann-to-Dirichlet" (NtD) nature for the wave field u since a local impedance operator links the Neumann traces of the solution u and the Dirichlet traces. From NtD ECs, we deduce "Dirichlet-to-Neumann" (DtN) ECs up to the fourth order of approximation. The construction of ECs relies on a derivation of a multiscale expansion for u in power series of δ. This expansion exhibits profile terms (defined in the second region) rapidly decaying and which describe the boundary layer phenomenon. One difficulty to validate the equivalent conditions lies in the proof of uniform energy estimates for the exact solution. We overcome this difficulty using a compactness argument and removing a discret set of resonant frequencies which may appear in the first region. As a consequence of these estimates we develop an argument for the convergence of the multiscale expansion by proving error estimates. We prove stability results for DtN ECs up to the fourth order of approximation. The proof is based on a compactness argument. Finally we can deduce convergence results for ECs from both error estimates and stability results.

The outline of the paper proceeds as follows. In Section 2 we introduce the mathematical model, we present briefly a formal derivation of equivalent conditions and we state uniform estimates for the solution of the transmission problem. In section 3 we prove uniform estimates for the solution of the transmission problem. In section 4, we exhibit the first terms of a multiscale expansion for the solution of the transmission problem and we prove error estimates. Then, we state a hierarchy of equivalent "Neumann-to-Dirichlet" boundary conditions up to the fourth order of approximation, we derive "Dirichlet-to-Neumann" ECs. In section 5 we present and prove stability results for ECs. In section 6, we present numerical results. In Appendix A we present elements of derivations for the multiscale expansion. In Appendix B, we present additional results and applications of this expansion. We measure the boundary layer phenomenon by introducing a "skin depth" function that turns out to depend on the Inria mean curvature of the interface and we make the comparison with the skin effect in electromagnetism.

The mathematical model, equivalent conditions and main results

In this section, we introduce the model problem and the framework, §2.1. We present a formal derivation of the approximate boundary conditions. Then we define the order of accuracy for equivalent conditions and we state uniform estimates for the solution of the exact problem.

The model problem. Framework

In this paper, we denote by (r, θ, φ) the classical spherical coordinates. Let Ω be a smooth bounded domain in R 3 , with boundary ∂Ω, and Ω -:= {r < r t } ⊂ Ω be a ball of radius r t , with boundary Σ. We denote by Ω + := Ω \ {r r t } the complementary of Ω -in Ω, see Figure 1. In this paper we investigate the Helmholtz equation

Ω - Ω + Σ n
-div 1 ρ ∇u - ω 2 ρc 2 (1 + iν)u = f in Ω .
where ρ is the density, c the celerity of acoustic waves, ω the pulsation, ν a damping parameter. Here the density ρ satisfies a specific assumption in the domain Ω + , that is :

ρ + := ρ| Ω+ is a radial function defined as ρ + (r) = γ exp (-α(r -r t )) ,
where the parameter α > 0 is large with respect to the wave number at the surface k = ω c and γ is a real constant. The coefficient ν takes two different values in Ω ± : ν = (0, ν) in Ω -× Ω + where ν ∈ R \ {0} (i.e. the domain Ω + corresponds to a dissipative media). For the sake of simplicity, we assume that the right-hand side f (which represents an acoustic source) is a smooth function and its support does not touch the domain Ω + . We will work under usual assumptions (regularity and positiveness) on the density and the celerity Hypothesis 2.2 (i) The density ρ and the celerity c are real valued smooth functions in Ω.

(ii) The density ρ and the celerity c are strictly positive in Ω.

A consequence of the hypothesis 2.2 is that ρ and c are bounded with strictly positive constants:

ρ 0,-≤ ρ -≤ ρ 1,-, c 0,-≤ c -≤ c 1,-, c 0,+ ≤ c + ≤ c 1,+
We introduce the small parameter δ

:= α -1 > 0. Then u = (u - δ , u + δ ) satisfies the following problem                        -div 1 ρ ∇u - δ - ω 2 ρc 2 u - δ = f in Ω -, -∆u + δ - 1 δ ∂ r u + δ - ω 2 c 2 (1 + iν)u + δ = 0 in Ω + , u + δ = u - δ on Σ , ∂ n u + δ = ∂ n u - δ on Σ , u + δ = 0 on ∂Ω . (2.1)
Our goal is to derive approximate models for the restriction u - δ of u δ to the interior domain Ω -when δ → 0. The limit problem (when δ tends to zero) has the following solution : u + 0 = 0 in the domain Ω + and u - 0 (in the domain Ω -) which satisfies the problem

   -div 1 ρ ∇u - 0 - ω 2 ρc 2 u - 0 = f in Ω -, u - 0 = 0 on Σ .
(2.2) Remark 2.1 The limit problem (2.2) with f = 0 may exhibit eigenfrequencies ω.

In this paper we prove stability and convergence results for approximate models (which are denoted by u δ k in Sec. 2.3) when δ → 0 by removing resonant frequencies which may appear in the domain Ω -, Rem 2.1. Therefore we work under the following spectral assumption.

Hypothesis 2.3

The angular frequency ω is non-zero and is not an eigenfrequency of the problem

   div 1 ρ ∇u + ω 2 ρc 2 u = 0 in Ω - u = 0 on Σ . (2.3)
In the framework above, we prove that it is possible to replace the domain Ω + by approximate boundary conditions set on Σ provided that the parameter δ is small enough.

Uniform estimate

We introduce a suitable variational framework for the solution of the problem (2.1) with more general right-hand sides. This framework is useful to prove error estimate (2.18).

Inria

Weak solutions

For given data (f -, f + , g) we consider the boundary value problem

                         -div 1 ρ ∇u - δ - ω 2 ρc 2 u - δ = f - in Ω - -div 1 ρ ∇u + δ - ω 2 ρc 2 (1 + iν)u + δ = f + in Ω + u + δ = u - δ on Σ ∂ n u + δ = ∂ n u - δ + g on Σ u + δ = 0 on ∂Ω .
(2.4)

Hereafter, we explicit a weak formulation of the problem (2.4). The variational problem writes : Find

u δ ∈ V = H 1 0 (Ω) such that ∀v ∈ V, a δ (u δ , v) = F, v V ,V , (2.5)
where the sesquilinear form a δ is defined as

a δ (u, v) = Ω± 1 ρ ± ∇u ± • ∇v ± dx -ω 2 Ω- 1 ρ -c 2 - u -v-dx -ω 2 (1 + iν) Ω+ 1 ρ + c 2 + u + v+ dx ,
(here the data ρ + depends on the parameter δ: ρ + (r) = γ exp -δ -1 (r -r t ) ) and the right-hand side F is defined as

F, v V ,V = - Ω fv dx - Σ gv dσ .
We assume that the data (f, g) satisfies the regularity assumption f ∈ L 2 (Ω) and g ∈ L 2 (Σ) .

(2.6)

Statement of uniform estimate

In the framework of Sec. 2.1 we prove δ-uniform a priori estimates for the solution of problem (2.5). The following theorem is the main result in this section.

Theorem 2.4 Under Hypothesis 2.2-2.3, there is a constant C > 0 such that for all δ > 0 the problem (2.5) with data (f, g) satisfying (2.6) has a unique solution u δ ∈ V which satisfies the uniform estimates

u - δ 1,Ω-+ 1 √ ρ+ u + δ 0,Ω+ + 1 √ ρ+ ∇u + δ 0,Ω+ C( f 0,Ω + g 0,Σ ) . (2.7)
This result is proved in Sec. 3. As an application of uniform estimate (2.7), we develop in Sec. 4.2 an argument for the convergence of the asymptotic expansion introduced in Sec. 2.3.

Formal derivation of equivalent conditions

In this sub-section, the main results of the paper are presented : equivalent boundary conditions (ECs) up to the fourth order of approximation. In the framework above, it is possible to derive high-order boundary conditions set on the interface Σ, when the parameter δ is small enough. We shall derive for all k ∈ {0, 1, 2, 3} an EC on Σ which is associated with the problem (2.1) and satisfied by u k δ , i.e. u k δ solves the problem

   -div 1 ρ ∇u k δ - ω 2 ρc 2 u k δ = f in Ω -, u k δ + D k,δ (∂ n u k δ ) = 0 on Σ .
(2.8)

Here D k,δ is a surface differential operator acting on functions defined on Σ and which depends on δ. In this section, we first present a multiscale expansion ansatz for the solution. Then Equivalent Conditions (ECs) up to the fourth order are derived for the solution of the exact problem. These conditions are of "Neumann-to-Dirichlet" (NtD) nature. Then we introduce "Dirichlet-to-Neumann" (DtN) conditions.

First step : a multiscale expansion

The first step consists in deriving a multiscale expansion for the solution (u - δ , u + δ ) of the model problem (2.1): it possesses an asymptotic expansion in power series of the small parameter δ

u - δ (x) = u - 0 (x) + δu - 1 (x) + δ 2 u - 2 (x) + δ 3 u - 3 (x) + • • • in Ω -, u + δ (x) = u + 0 (x; δ) + δu + 1 (x; δ) + δ 2 u + 2 (x; δ) + δ 3 u + 3 (x; δ) + • • • in Ω + , with u + j (x; δ) = χ(r)U j (θ, φ, r -r t δ ) .
The term U j (θ, φ, S) is a "profile" defined on Σ×R + , and (r, θ, φ) are the classical spherical coordinates.

The function x → χ(r) is a smooth cut-off with support in Ω + and equal to 1 in a tubular neighborhood of Σ. This cut-off function is needed to satisfy Dirichlet condition u = 0 on the external boundary of Ω.

Each profile U j tends to 0 when S → +∞. We make explicit the first asymptotics (u - j , U j ) for j = 0, 1, 2, 3 in Section 4.1.

Second step : construction of Neumann-to-Dirichlet (NtD) equivalent conditions of order k + 1

The second step consists in identifying for all k ∈ {0, 1, 2, 3} a simpler problem satisfied, up to a residual term in O(δ k+1 ), by the truncated expansion

u - k,δ := u - 0 + δu - 1 + • • • + δ k u - k
We define u δ k as the solution of the simpler problem:

   -div 1 ρ ∇u δ k - ω 2 ρc 2 u δ k = f in Ω -, u δ k + D k,δ (∂ n u δ k ) = 0 on Σ .
This simpler problem is the problem set in Ω -with the equivalent boundary conditions (EC) set on Σ.

Here f is the right-hand side of problem (2.1) and D k,δ is a differential operator acting on functions defined on the surface Σ and which depends on δ:

D 0,δ = 0 , (2.9)
D 1,δ = δ I , (2.10)
D 2,δ = δ 1 - 2δ r t I , (2.11)
D 3,δ = δ 1 - 2δ r t + δ 2 6 r 2 t + ω 2 c 2 0 (1 + iν) + ∆ Σ .
(2.12) c 0 is the Dirichlet trace of c on Σ. ∆ Σ is the Laplace-Beltrami operator along the surface Σ. It writes in spherical coordinates (θ, φ) as

∆ Σ = 1 r 2 t sin θ ∂ θ (sin θ∂ θ ) + 1 r 2 t sin 2 θ ∂ 2 φ . (2.13) Inria Remark 2.
2 For all k ∈ {0, 1, 2}, D k,δ is a scalar operator whereas D 3,δ is a second order differential operator.

This second step and the construction of equivalent conditions are detailed in Section 4.3.

Third step : Dirichlet-to-Neumann (DtN) conditions

A Dirichlet-to-Neumann equivalent condition of order k + 1 writes

N k,δ (v δ k ) + ∂ n v δ k = 0 on Σ ,
where N k,δ denotes a local approximation of the operator N δ , and the simpler problem writes

     -div 1 ρ ∇v δ k - ω 2 ρc 2 v δ k = f in Ω -, N k,δ (v δ k ) + ∂ n v δ k = 0 on Σ . (2.14)
We compute in section 4.4 the expressions of N 1,δ , N 2,δ and N 3,δ

N 1,δ = δ -1 I , (2.15)
N 2,δ = δ -1 1 + 2δ r t I , (2.16
)

N 3,δ = δ -1 (1 + 2δ r t ) I -δ 2 2 r 2 t + ω 2 c 2 0 (1 + iν) + ∆ Σ . (2.17)
For shortness of notations, we define a DtN condition of order 1 as the Dirichlet BC on Σ.

Stability and convergence of NtD Equivalent conditions

We present well-posedness and convergence results for NtD conditions, Sec. 2.4. Elements of derivation and mathematical validations for ECs are presented in Sec. 4 and Sec. 5.

Our goal in the next sections is to validate ECs set on Σ (Sec. 4.4) proving δ-estimates for u δ -u k δ (k ∈ {0, 1, 2, 3}) where u k δ is the solution of the approximate model (2.14) , and u δ satisfies the problem (2.1).The functional setting for u k δ is described by the Hilbert space

V k : Notation 2.5 V k denotes the space H 1 0 (Ω -) when k = 0; V k = H 1 (Ω -) when k = 1, 2, and V 3 = {u ∈ H 1 (Ω -) | u| Σ ∈ H 1 (Σ)}
The main result of this section is the following statement, that is for all k ∈ {0, 1, 2, 3} the problem (2.14) is well-posed in the space V k , and its solution satisfies uniform H 1 error estimates.

Theorem 2.6 Under Hypothesis 2.2-2.3, for all k ∈ {0, 1, 2, 3} there are constants δ k , C k > 0 such that for all δ ∈ (0, δ k ), the problem (2.14) with a data f ∈ L 2 (Ω -) has a unique solution u k δ ∈ V k which satisfies the uniform estimates: u δ -u k δ 1,Ω- C k δ k+1 . (2.18)
The stability result for the problem (2.8) is stated in Thm. 5.1 and it is proved in Sec. 5.2. It appears nontrivial to work straightforwardly with the difference u δ -u k δ . A usual method consists in using the truncated series u k,δ introduced in Sec. 4.3 as intermediate quantities [START_REF] Haddar | Generalized impedance boundary conditions for scattering problems from strongly absorbing obstacles: the case of Maxwell's equations[END_REF][START_REF] Péron | Equivalent boundary conditions for an elasto-acoustic problem set in a domain with a thin layer[END_REF]. Then, the error analysis is split into two steps :

1. We prove a uniform estimate for the remainder u δ -u k,δ in Thm. 4.1, Sec. 4.2 , 2. The second step consists of proving a uniform estimate for the difference

u k,δ -u k δ u k,δ -u k δ 1,Ω- C k δ k+1 .
The first step is detailed in the next sections. The second step can be deduced straightforwardly from stability estimates (5.1a)-(5.1b) stated in Th. 5.1 (Sec. 5.1) since by construction the quantity u k,δ -u k δ satisfies Problem (2.14) with a right-hand side of order δ k+1 and which has a support on the surface Σ. Estimates (5.1a)-(5.1b) are proved in Section 5. We refer the reader to [17, Sec. 5.2]- [START_REF] Haddar | Generalized impedance boundary conditions for scattering problems from strongly absorbing obstacles: the case of Maxwell's equations[END_REF] where the authors adopt the same method and use similar arguments.

Uniform estimate for the reference solution

In this section, we prove the Thm. 2.4, that is a uniform estimate for the solution of the reference problem. We consider the problem (2.4) at a fixed pulsation ω satisfying Hypothesis 2.3. We are going to prove the following Lemma Lemma 3.1 Under Hypothesis 2.2-2.3, there is a constant C > 0 such that for all δ > 0 any solution

u δ ∈ H 1 0 (Ω) of problem (2.5
) with a data f ∈ L 2 (Ω) satisfies the uniform estimate with g = 0 :

u δ 0,Ω C f 0,Ω (3.1)
We conjecture that the following estimate holds when g = 0:

u δ 0,Ω + u δ 0,Σ C( f 0,Ω + g 0,Σ )
This Lemma is going to be proved in this section. The proof of this result involves both a compactness argument and the spectral hypothesis 2.3. Then as a consequence of estimate (3.1), we deduce estimate (2.7). Finally Thm. 2.4 is obtained as a consequence of the Fredholm alternative since the problem (2.5) is of Fredholm type.

Proof of Lemma 3.1 : Uniform estimate of u δ

We prove this lemma by contradiction : We assume that there exists a sequence {u m } ∈ V , m ∈ N, of solutions of problem (2.5), each solution being associated with a parameter δ m > 0 and a right-hand side

f m ∈ L 2 : ∀v ∈ V, Ω± 1 ρ ± ∇u ± m • ∇v ± dx -ω 2 Ω- 1 ρ -c 2 - u - m v-dx -ω 2 (1 + iν) Ω+ 1 ρ + c 2 + u + m v+ dx = Ω f m v dx , (3.2)
satisfying the following conditions

u m 0,Ω = 1 for all m ∈ N, (3.3a) f m 0,Ω → 0 as m → ∞. (3.3b) Choosing tests functions v = u m in (3.
2) and taking the imaginary part, since ν = 0, we obtain with the help of conditions (3.3a)-(3.3b) the following convergence result :

Inria (i) the sequence { 1 √ ρ+ u + m } converges to 0 in L 2 (Ω + ).
Then taking the real part, we obtain as a consequence of this convergence result the following uniform bounds on the gradients :

∇u - m 0,Ω-+ 1 √ ρ+ ∇u + m 0,Ω+ C (3.4)
Since ρ + ρ 0,-and ρ -is bounded, we deduce ∇u m 0,Ω C and the sequence {u m } is bounded in

H 1 .

Limit of the sequence and conclusion

The domain Ω being bounded, the embedding of H 1 in L 2 is compact. Hence, since the sequence {u m } is bounded in H 1 we can extract a subsequence of {u m } (still denoted by {u m }) which is strongly converging in L 2 (Ω) and weakly converging in

H 1 ∇u m ∇u in L 2 (Ω) = (L 2 (Ω)) 3 u m → u in L 2 (Ω). (3.5)
A consequence of the strong convergence in L 2 (Ω) and (3.3a) is that u 0,Ω = 1. Using Hypothesis 2.3, we are going to prove that u = 0, which will contradict u 0,Ω = 1, and finally prove estimate (3.1).

Since

u + m 0,Ω+ √ ρ 0,- 1
√ ρ+ u + m 0,Ω+ , letting m → +∞ and using (i) we get u + 0,Ω+ = 0. Hence u = 0 in Ω + . (3.6)
In particular, (3.6) implies that u -:= u| Ω-belongs to the space

H 1 0 (Ω -). Let w ∈ H 1 0 (Ω -).
The extension w 0 of w by 0 on Ω + defines an element of H 1 0 (Ω). We can use w 0 as test function in (3.2) and we obtain

Ω- 1 ρ - ∇u - m • ∇w -dx -ω 2 Ω- 1 ρ -c 2 - u - m w -dx = Ω f - m w dx .
According to (3.5) and (3.3b), taking limits as m → +∞, we deduce from the previous equalities

Ω- 1 ρ - ∇u -• ∇w dx -ω 2 Ω- 1 ρ -c 2 - u -w dx = 0 , (3.7)
i.e., u -∈ H 1 0 (Ω -) satisfies (3.7) for all w ∈ H 1 0 (Ω -). Then integrating by parts we find

   div( 1 ρ - ∇u -) + ω 2 ρ-c 2 - u -= 0 in Ω - u -= 0 on Σ . (3.8)
By Hypothesis 2.3, we deduce

u -= 0 in Ω -.
Hence, with (3.6), we have u = 0 in Ω, which contradicts u 0,Ω = 1 and ends the proof of Lemma 3.1.

Derivation and stability of equivalent conditions

It is possible to exhibit series expansions in powers of δ for u δ in (2.1) :

u - δ (x) ≈ j 0 δ j u - j (x) , (4.1)
u + δ (x) ≈ j 0 δ j u + j (x; δ) with u + j (x; δ) = χ(r)U j (θ, φ, r -r t δ ) . (4.2)
The term U j (θ, φ, S) is a "profile" which is defined on Σ × R + , it tends to 0 when S → +∞. The function χ is a smooth cut-off, see Sec. 2.3. Elements of derivations for this expansion are presented in detail in Appendix A. In this section, we explicit the first terms in asymptotics, §4.1. Then we construct formally equivalent conditions, §4.3. In §4.2 we validate the asymptotic expansion with estimates for the remainders.

First terms of the asymptotic expansion

Straightforward calculations lead to the first-order terms U j and u - j . In this section we explicit these terms when j ≤ 3 . We give elements of proof in Appendix A, Section A.3. There holds

u + 0 = 0 , (U 0 = 0).
Then u - 0 solves the problem

   -div( 1 ρ ∇u - 0 ) -ω 2 ρc 2 u - 0 = f in Ω -, u - 0 = 0 on Σ . (4.3)
The next term which is determined in the asymptotics is the profile

U 1 U 1 (θ, φ, S) = -e -S ∂ n u - 0 X(r t θ, r t φ) . (4.4)
Here, X(r t θ, r t φ) = (r t sin θ cos φ, r t sin θ sin φ, r t cos θ) ∈ Σ .

Then u -

1 solves    -div( 1 ρ ∇u - 1 ) -ω 2 ρc 2 u - 1 = 0 in Ω -, u - 1 = -∂ n u - 0 on Σ . (4.6)
The next term which is determined is the profile

U 2 U 2 (θ, φ, S) = (a 2 (θ, φ) + Sb 2 (θ, φ))e -S , (4.7)
where

a 2 (θ, φ) = -∂ n u - 1 + 2 r t ∂ n u - 0 X(r t θ, r t φ) and b 2 (θ, φ) = 2 r t ∂ n u - 0 X(r t θ, r t φ) .
Then u - 2 solves the problem

     -div( 1 ρ ∇u - 2 ) -ω 2 ρc 2 u - 2 = 0 in Ω -, u - 2 = -∂ n u - 1 + 2 r t ∂ n u - 0 on Σ . (4.8)
The next term which is determined is the profile

U 3 U 3 (θ, φ, S) = (a 3 (θ, φ) + Sb 3 (θ, φ) + S 2 c 3 (θ, φ))e -S , (4.9)
Inria where the functions a 3 , b 3 , c 3 are given by

                   a 3 = -∂ n u - 2 + 2 r t ∂ n u - 1 - 6 r 2 t + ω 2 c 2 0 (1 + iν) + ∆ Σ ∂ n u - 0 , b 3 = 2 r t ∂ n u - 1 - 6 r 2 t + ω 2 c 2 0 (1 + iν) + ∆ Σ ∂ n u - 0 , c 3 = - 3 r 2 t ∂ n u - 0 .
Then, u - 3 solves the problem

       -div( 1 ρ ∇u - 3 ) -ω 2 ρc 2 u - 3 = 0 in Ω -, u - 3 = -∂ n u - 2 + 2 r t ∂ n u - 1 - 6 r 2 t + ω 2 c 2 0 (1 + iν) + ∆ Σ ∂ n u - 0 on Σ . (4.10)
Existence and regularity results in Sobolev spaces for the asymptotics U j and u - j are stated in Prop. A.1, Sec. A.4.

Estimates for the remainders

The validation of the asymptotic expansion (4.1)-(4.2) consists in proving estimates for remainders r N δ defined in Ω as

r N δ = u δ - N n=0 δ n u n (4.11)
The convergence result is the following statement.

- n ∈ H 1 (Ω -) and U n ∈ H 1 Σ × R + . For each N ∈ N, the remainders r N δ satisfy r N,- δ 1,Ω-+ δ -1/2 1 √ ρ+ r N,+ δ 0,Ω+ + δ 1/2 1 √ ρ+ ∇r N,+ δ 0,Ω+ C N δ N +1 (4.12)
with a constant C N independent of δ.

PROOF. The error estimate (4.12) is obtained through an evaluation of the right-hand sides when the Helmholtz operator is applied to the remainder r N δ . By construction, the remainder is solving the problem:

                         -div 1 ρ ∇r N,- δ - ω 2 ρc 2 r N,- δ = 0 in Ω - -div 1 ρ ∇r N,+ δ - ω 2 ρc 2 (1 + iν)r N,+ δ = f N,δ in Ω + r N,+ δ = r N,- δ on Σ ∂ n r N,+ δ = ∂ n r N,- δ + g N,δ on Σ r N δ = 0 on ∂Ω , (4.13)
RR n°8954

which is nothing but the boundary value problem (2.4) with right-hand sides (f -, f + , g) = (0, f N,δ , g N,δ ).

The right hand sides g N,δ and f N,δ are residuals which are roughly of the order δ N . By construction of the expansion, we obtain :

g N,δ = δ N ∂ n u - N on Σ (4.14)
and

f 0,δ = 0 and for all N 1, f N,δ = -χρ -1 + δ N -1 A 0 U N +1 (θ, φ, r -r t δ ) + O(δ) = O(δ N -1 )
(4.15) where the operator A 0 is defined in appendix (A.2a). The terms with derivatives of χ should appear with higher powers of δ (see [START_REF] Caloz | Uniform estimates for transmission problems with high contrast in heat conduction and electromagnetism[END_REF]Eq. (5.6)] and [16, Prop. 7.4, Eq. (7.17) 2 ]). We can apply Theorem 2.4 to the remainder r N δ :

r N,- δ 1,Ω-+ 1 √ ρ+ r N,+ δ 0,Ω+ + 1 √ ρ+ ∇r N,+ δ 0,Ω+ C N f N,δ 0,Ω+ + g N,δ 0,Σ (4.16)
We deduce from (4.14) and (4.15) the estimates

g N,δ 0,Σ Cδ N and f N,δ 0,Ω+ Cδ N -1/2 since A 0 U N +1 ∈ H 1 (Σ × R + ).
Here C may depend on N . Combining this estimate with (4.16), we infer

r N,- δ 1,Ω-+ 1 √ ρ+ r N,+ δ 0,Ω+ + 1 √ ρ+ ∇r N,+ δ 0,Ω+ C N δ N -1/2
we therefore have:

r N,- δ 1,Ω-+ δ -1/2 1 √ ρ+ r N,+ δ 0,Ω+ + δ 1/2 1 √ ρ+ ∇r N,+ δ 0,Ω+ Cδ N -1
In order to prove optimal estimates for r N δ , we apply the previous estimates to r N +2 δ . Using the formula

r N +2 δ = r N δ + δ N +1 u N +1 + δ N +2 u N +2
we finally obtain the wanted estimates since we have for any n ∈ N

u - n 1,Ω-+ δ -1/2 1 √ ρ+ u + n 0,Ω+ + δ 1/2 1 √ ρ+ ∇u + n 0,Ω+ C .

Construction of Neumann to Dirichlet (NtD) equivalent conditions

In this section, we derive formally equivalent boundary conditions (Sec. 2.3).

Order 1

Since u - 0 solves the problem (4.3), the condition of order 1 is the homogeneous Dirichlet boundary condition, see (2.9)

u 0 = 0 on Σ Inria Order 2
According to (4.3)-(4.6), the truncated expansion u - 1,δ := u - 0 + δu - 1 solves the Helmholtz equation

-div( 1 ρ ∇u - 1,δ ) -ω 2 ρc 2 u - 1,δ = f in the domain Ω -together with the Robin boundary condition set on Σ u - 1,δ + δ∂ n u - 1,δ = δ 2 ∂ n u - 1 on Σ
Neglecting the term of order 2 in the previous right-hand side, we infer the equivalent boundary condition, see (2.10)

u δ 1 + δ∂ n u δ 1 = 0 on Σ Order 3
According to (4.3)-(4.6)-(4.8), the truncated expansion u - 2,δ := u - 0 + δu - 1 + δ 2 u - 2 solves the Helmholtz equation in Ω -together with the boundary condition set on Σ

u - 2,δ + δ 1 - 2 r t δ ∂ n u - 2,δ = - 2 r t δ 3 ∂ n u - 1 + δ 3 1 - 2 r t δ ∂ n u - 2 on Σ
We neglect the terms of order 3 in the previous right-hand side and we infer the equivalent boundary condition of order 3, see (2.11)

u δ 2 + δ 1 - 2δ r t ∂ n u δ 2 = 0 on Σ Order 4
According to (4.3)-(4.6)-(4.8)-(4.10), the truncated expansion u - 3,δ := u - 0 + δu - 1 + δ 2 u - 2 + δ 3 u - 3 solves the Helmholtz equation in Ω -together with a boundary condition set on Σ

u - 3,δ + δ 1 - 2δ r t + δ 2 6 r 2 t + ω 2 c 2 0 (1 + iν) + ∆ Σ ∂ n u - 3,δ = O(δ 4 )
We neglect the terms of order 4 in the previous right-hand side to infer the equivalent boundary condition of order 4, see (2.12)

u δ 3 + δ 1 - 2δ r t + δ 2 6 r 2 t + ω 2 c 2 0 (1 + iν) + ∆ Σ ∂ n u δ 3 = 0

Dirichlet-to-Neumann (DtN) conditions

We introduce the operator N δ = (D δ ) -1 where the operator D δ is the exact Neumann-to-Dirichlet operator that characterizes u + δ . Then the exact boundary condition for u - δ can be rewritten as

N δ u - δ + ∂ n u - δ = 0 on Σ . A Dirichlet-to-Neumann equivalent condition of order k + 1 writes N k,δ (v δ k ) + ∂ n v δ k = 0 on Σ ,
where N k,δ denotes a local approximation of the operator N δ , and the simpler problem writes

     -div 1 ρ ∇v δ k - ω 2 ρc 2 v δ k = f in Ω -, N k,δ (v δ k ) + ∂ n v δ k = 0 on Σ .
(4.17)

RR n°8954

Remark 4.1 There holds N 1,δ = δ -1 I, i.e the DtN condition of order 2 coincides with the NtD condition of order 2 (2.10).

The expression of N k,δ , for any k 2, can be obtained by a formal Taylor expansion of

(D k,δ ) -1 such that D k,δ = (N k,δ ) -1 + O(δ k+1 ) .
Straightforward calculi lead to the following expansions

(D 2,δ ) -1 = δ -1 1 + 2δ r t I + O(δ) , (D 3,δ ) -1 = δ -1 (1 + 2δ r t ) I -δ 2 2 r 2 t + ω 2 c 2 0 (1 + iν) + ∆ Σ + O(δ 2 ) .
We infer the expressions of N 2,δ and N 3,δ

N 2,δ = δ -1 1 + 2δ r t I , N 3,δ = δ -1 (1 + 2δ r t ) I -δ 2 2 r 2 t + ω 2 c 2 0 (1 + iν) + ∆ Σ .
For shortness of notations, we define a DtN condition of order 1 as the Dirichlet BC on Σ (2.9).

Analysis of DtN Equivalent Conditions

The main result of this section is Theorem 5.1 in §5.1 which proves the stability of equivalent conditions.

Main result for DtN equivalent conditions (ECs)

In this section we address the issue of stability results for ECs. We focus on DtN equivalent condition which are defined for k = 1, 2, 3 (2.14). For shortness of notations we define for k = 0 the DtN EC of order 1 as the Dirichlet BC on Σ. We consider the problem (2.14) at a fixed non-zero frequency ω satisfying Hypothesis 2.3. The functional setting for u k δ is described by the Hilbert space V k (see Notation 2.5). The main result of this section is the following statement, that is for all k ∈ {0, 1, 2, 3} the problem (2.14) is well-posed in the space V k , and its solution satisfies uniform H 1 estimates. Theorem 5.1 Under Hypothesis 2.2-2.3, for all k ∈ {0, 1, 2, 3} there are constants δ k , C k > 0 such that for all δ ∈ (0, δ k ), the problem (2.14) with a data f ∈ L 2 (Ω -) has a unique solution u k δ ∈ V k which satisfies the uniform estimates:

u k δ 1,Ω- C k f 0,Ω-for all k ∈ {0, 1, 2} ,
(5.1a)

u 3 δ 1,Ω-+ √ δ ∇ Σ u 3 δ 0,Σ + √ δ -1 u 3 δ 0,Σ C 3 f 0,Ω-.
(5.1b)

The key for the proof of Thm. 5.1 is the following Lemma.

Lemma 5.2 Under Hypothesis 2.2-2.3, for all k ∈ {0, 1, 2, 3} there exists constants δ k , C k > 0 such that for all δ ∈ (0, δ k ), any solution u k δ ∈ V k of problem (2.14) with a data f ∈ L 2 (Ω -) satisfies the uniform estimate:

u k δ 0,Ω- C k f 0,Ω- (5.2)
Inria Remark 5.1 For k = 0, the Theorem 5.1 and the Lemma 5.2 hold for all δ > 0 since the EC is nothing but the Dirichlet BC on Σ. For k = 1, 2, 3, one proves uniform estimate when δ is small enough. The proof is based on a compactness argument.

The Lemma 5.2 is proved in Section 5.2. As a consequence of this Lemma, any solution of the problem (2.14) satisfies uniform H 1 -estimates (5.1a), (5.1b), respectively when k = 0, 1, 2, 3. Then, the proof of the Thm. 5.1 is obtained as a consequence of the Fredholm alternative since the problem (2.14) is of Fredholm type.

In this section, one first proves the Lemma 5.2, i.e. uniform L 2 -estimate (5.2) for any solution of problem (2.14), Sec. 5.2. For k = 0, the proof of Lemma 5.2 is a consequence of the Hypothesis 2.3 together with the Fredholm alternative. We focus on the proof of Lemma 5.2 for k = 3 since the proofs when k = 0, 1, 2 are simpler. Hence we consider the problem (here

u = u 3 δ )    -div( 1 ρ ∇u) -ω 2 ρc 2 u = f in Ω -, ∂ n u + N δ u = 0 on Σ , (5.3)
where the operator N δ is defined as

N δ = δ -1 J δ I -δ∆ Σ
Here, J δ is a scalar defined as

J δ = (1 + 2δ r t ) -δ 2 2 r 2 t + ω 2 c 2 0 (1 + iν) .
To prepare for the proof, we introduce the variational formulation for u. If u ∈ V 3 is a solution of (5.3), then it satisfies for all v ∈ V 3 :

Ω- 1 ρ ∇u • ∇v -ω 2 ρc 2 uv dx + δ ρ(rt) Σ ∇ Σ u • ∇ Σ v dσ + δ -1 ρ(rt) J δ Σ u v dσ = Ω- fv dx . (5.4)

Proof of Lemma 5.2

Proof by contradiction: We assume that there is a sequence (u m ) ∈ V 3 , m ∈ N, of solutions of the problem (5.3) associated with a parameter δ m and a right-hand side

f m ∈ L 2 (Ω -): -div( 1 ρ ∇u m ) -ω 2 ρc 2 u m = f m in Ω -, (5.5a)
∂ n u m + δ -1 m J m u m -δ m ∆ Σ u m = 0 on Σ , (5.5b)
(with J m := J δm ) satisfying the following conditions

δ m → 0 as m → ∞ , (5.6a) u m 0,Ω-= 1 for all m ∈ N , (5.6b) f m 0,Ω-→ 0 as m → ∞ .
(5.6c)

Estimates of the sequence {u m }

We first prove that the sequence {u m } is bounded in H 1 (Ω -) (not in V 3 ). We particularize the weak formulation (5.4) for the sequence {u m }: For all v ∈ V 3 :

Ω- 1 ρ ∇u m • ∇v -ω 2 ρc 2 u m v dx + δm ρ(rt) Σ ∇ Σ u m • ∇ Σ v dσ + δ -1 m ρ(rt) J m Σ u m v dσ = Ω- f m v dx . (5.7)
Then choosing v = u m in (5.7) and taking the real part, we obtain with the help of condition (5.6b) the uniform bound (since

1 ρc ∈ L ∞ (Ω -)) Ω- 1 ρ - ∇u m 2 dx + δm ρ(rt) ∇ Σ u m 2 0,Σ + δ -1 m ρ(rt) (J m ) u m 2 0,Σ ω 2 1 ρc ∞ + f m 0,Ω-.
We infer that the sequence

{u m }, resp. { δ m ∇ Σ u m }, is bounded in H 1 (Ω -) (since 1 ρ - D 0 > 0), resp. in L 2 (Σ); and since J m tends to 1, the sequence {δ -1/2 m u m } is bounded in L 2 (Σ) : u m 1,Ω-C , (5.8a)
(δ m ) 1 2 ∇ Σ u m 0,Σ C , (5.8b)
(δ m ) -1 2 u m 0,Σ C . (5.8c)
Another consequence of (5.8a) is that the sequence {u m } is bounded in

H 1 2 (Σ). u m 1 2 ,Σ
C .

(5.9)

Limit of the sequence and conclusion

The domain Ω -being bounded, the embedding of H 1 (Ω -) in L 2 (Ω -) is compact. Hence as a consequence of (5.8a) we can extract a subsequence of {u m } (still denoted by {u m }) which is converging in L 2 (Ω -) and we can assume that the sequence {∇u m } is weakly converging in L 2 (Ω -). As a consequence of (5.9), up to the extraction of a subsequence, we can also assume that the sequence {u m } is strongly converging in L 2 (Σ) : We deduce that there is u ∈ L 2 (Ω -) such that

u m → u in L 2 (Ω -) (5.10a) ∇u m ∇u in L 2 (Ω -) (5.10b)
u m → u in L 2 (Σ) . (5.10c)
A consequence of the strong convergence result in L 2 (Ω -) (5.10a) and (5.6b) is that u 0,Ω-= 1. Furthermore, as a consequence of (5.10c) and (5.8c) one can deduce that u = 0 on Σ. Using Hypothesis 2.3, we are going to prove that u = 0 in Ω -, which will contradict u 0,Ω-= 1 and finally prove estimate (5.2). Let v ∈ C ∞ c (Ω -), it defines an element of V 3 . Hence we can use v as test function in (5.7). Since v = 0 on Σ and ∇ Σ v = 0 on Σ we infer

Ω- 1 ρ ∇u m • ∇v -ω 2 ρc 2 u m v dx = Ω- f m v dx . (5.11)

Inria

As a consequence of (5.10c), taking limits as m → ∞, we can deduce from the previous equality and convergence results (5.10a)-(5.10b)- (5.11)

that u ∈ H 1 0 (Ω -) satisfies for all v ∈ C ∞ c (Ω -) Ω- 1 ρ ∇u • ∇v -ω 2 ρc 2 uv dx = 0 .
Then integrating by parts we find that u satisfies the problem

   -div( 1 ρ ∇u) -ω 2 ρc 2 u = 0 in Ω - u = 0 on Σ .
According to Hypothesis 2.3, we infer u = 0 in Ω - which contradicts u 0,Ω-= 1 and ends the proof of Lemma 5.2.

6 Numerical illustration

Numerical method

We solve Helmholtz equation in a radial configuration (r, θ, φ), where the data only depends on the radius r. Details and validations of the numerical method described in this section are given in [START_REF] Chabassier | High Order Finite Element Method for solving Convected Helmholtz equation in radial and axisymmetric domains. Application to Helioseismology[END_REF]. The interval [0, r t ] is subdivided into sub-intervals:

[0, r t ] = ∪[x i , x i+1 ]
One-dimensional finite elements will be used in r-coordinate, while spherical harmonics will be used in θ, φ. Since spherical harmonics are orthonormal and diagonalize the laplacian, we will obtain a decoupled sequence of 1-D problems to solve. The 1-D finite element space is equal to

V h = u ∈ H 1 (Ω) such that u| [xi,xi+1] ∈ P r
where P r is the space of polynomials of degree lower or equal to r, r being the order of the approximation. The solution u is then searched under the form

u(r, θ, φ) = N h j=0 L =0 m=- u ,m j ϕ j (r) Y m (θ, φ)
where ϕ i are basis functions generating the finite element space V h of dimension N h , and Y m spherical harmonics given as

Y m (θ, φ) = (-1) m (2 + 1) 4π ( -m)! ( + m)! P m (cos θ)e imφ
where P m are the associated Legendre polynomials. L is the maximal degree of spherical harmonics used in the expression of u. The variational formulation solved by u ,m is given as

-ω 2 rt 0 1 ρ 0 c 2 0 r 2 u ,m ϕ i dr + rt 0 r 2 ρ 0 ∂u ,m
∂r ∂ϕ i ∂r dr

+ ( + 1)

rt 0 1 ρ 0 u ,m ϕ i dr - 1 ρ 0 r 2 ∂u ,m ∂r ϕ i rt 0 = f ,m i RR n°8954
where

f ,m i = rt 0 π 0 2π 0 r 2 f (r, θ, φ) Ȳ m (θ, φ) sin θ dφ dθ ϕ i dr
The boundary term in square brackets of the variational formulation is replaced by the correct term depending on the boundary condition imposed at r = r t . For the discretization of the finite element space V h , tenth order polynomials (r = 10) are used. Gauss-Lobatto points are used both as interpolation and quadrature points. The source f ,m i is computed with Gauss-Legendre integration formulae.

Radial source and solution

We consider an artificial data of a star of interior radius equal to 1 with constant celerity and piecewise exponential density. More precisely, ρ 0 (r) = D 1 e (r-1)/δin if r < 1 D 1 e -(r-1)/δ if r ≥ 1 (6.1)

D 1 = 2 × 10 -4 kg/m 3 , c 0 = 8 × 10 3 m/s, δ in = +∞ (6.2) ω = 1 + 2iγ ω 0 ω 0 , γ = ω 0 100 (6.3)
The parameters values have been chosen of the same order of magnitude than the Sun's values after adimensionalization with respect to the radius, in order to face the same numerical difficulties. δ in = +∞ means that the interior density is constant. Two typical frequencies are chosen : ω 0 = 2πf 0 with f 0 = 7 × 10 4 and 7 × 10 5 Hz (which corresponds for the Sun, before the adimensionalisation procedure to f sun 0 = 0.1 mHz and 1 mHz). They are related to the interior wavelengths λ = c 0 /f 0 : 1.14 × 10 -1 m and 1.14 × 10 -2 m. A radial gaussian source is located at the radius r 0 = 0.8 m :

f (r) = e -α(r-r0) 2
As a result the computations are performed for a single mode = 0. The parameter α is chosen such that f (r 0 + 0.02) = 10 -6

We choose to evaluate the performances of our equivalent atmosphere conditions for values of δ related to the wavelength of the considered problem. For f sun The value of ∂w ∂r (r t ) will be the exact impedance Z exact , and we impose

0 = 0.1 mHz, δ ∈ [5 • 10 -5 , 10 -2 ] m. For f sun 0 = 1 mHz, δ ∈ [5 • 10 -6 , 10 -3 ] m.

Reference solutions

∂u ∂r (r t ) = Z exact u(r t )
as a boundary condition to compute the reference solution. A plane wave analysis provides two complex wavelengths:

2k ± = iα ± 4k 2 ∞ -α 2 (6.5)
The right-going wave will oscillate (for r large enough) as e ik + r whereas a left-going wave will oscillate as e ik -r . The problem (6.4) is solved in the domain [r t , R max ] where R max is chosen such that

|e i(k + -k -)(R max -rt) | = 10 -16
since the left term of this relation corresponds to the amplitude of the reflected wave. This exponential is always decreasing as soon as γ is positive. The problem (6.4) is solved by searching w as

w = v exp(ik + r)
The advantage of this approach is that less degrees of freedom are needed to discretize v (since the phase is removed), and no underflow occurs.

Equivalent Atmosphere Boundary Conditions

For each value of δ, each of the equivalent atmosphere boundary conditions associated with the problem (2.14) is imposed at r = r t , leading to four approximate numerical solutions (called Dirichlet, Atmo1, Atmo2 and Atmo3). Dirichlet corresponds to condition (2.9), Atmo1 to condition (2.15), Atmo2 to the condition (2.16) and Atmo3 to condition (2.17

Axisymmetric source and solutions

In a more realistic setting, we consider smooth approximations of ρ sun , c sun obtained by fitting eighthorder B-splines with the data given by the Standard Solar Model [START_REF] Christensen-Dalsgaard | The current state of solar modeling[END_REF] (c and log ρ are fitted). The approximation of ρ sun is plotted in figure 3, R is the radius of the Sun. The density is therefore given

ρ 0 (r) = ρ sun if r/R ≤ 1.0007126 D 1 e -(r/R -1)/δ otherwise (6.6) D 1 ≈ 3.292 × 10 -6 kg/m 3 , (6.7
)

ω = 1 + 2iγ ω 0 ω 0 , γ = ω 0 100 (6.8)
The velocity is given as

c 0 (r) = c sun if r/R ≤ 1.0007126 C 1 otherwise (6.9) C 1 ≈ 6.865 × 10 3 m/s (6.10)
The source is a 3-D gaussian:

f (x) = exp(-α||x -x 0 || 2 )
where x 0 ≈ (0, 0, 0.99) and α chosen such that exp(-0.007 α) = 10 -6 .

The computation is done with spherical harmonics with ≤ 100 for a frequency f sun 0 = 3 mHz. The reference solution is displayed in figure 4 on a plane Oxz for a given value of δ. In figure 5(a), the convergence of the different equivalent boundary conditions is represented versus δ. For this frequency and these parameters, we observe a convergence in O(δ 3 ) for equivalent boundary conditions Atmo1 and Atmo2, and in O(δ 5 ) for Atmo3. We think that we are in a preasymptotic regime. In figure 5(b), the different equivalent boundary conditions are compared for different frequencies from 0.3 mHz until 9 mHz. The condition Atmo3 is quite efficient for low frequencies. We observe that all the conditions (Dirichlet, Atmo1, Atmo2 and Atmo3) are no longer accurate for high frequencies (f Sun 0 ≥ 5 mHz). It is expected since for this range of frequencies, propagative modes exist inside the atmosphere and they are no longer restricted to a small tubular neighborhood of the atmosphere. Future works will address the derivation of accurate equivalent boundary conditions in high-frequency regime.

Inria

A Elements of derivations for the multiscale expansion and equivalent conditions

The aim of this appendix is to present elements of derivations for the multiscale expansion and equivalent conditions. We make the following ansatz

u - δ ∼ j 0 δ j u - j (x)
u + δ ∼ j 0 δ j u + j (x; δ) with u + j (x; δ) = χ(r)U j (θ, φ, r -r t δ )

for the solution of the problem (2.1). We present all the details for the calculi of the first terms u ± j (j = 0, • • • , 3). As a consequence, we exhibit formally equivalent boundary conditions up to the fourth order.

A.1 Expansion of the operators in power series of δ

To derive the formal expansion of u + δ we perform the scaling

S = r -r t δ (A.1)
in the equations set in Ω + and Σ in order to make appear the small parameter δ in these equations. Then it is possible to expand the Helmholtz operator in power series of δ with coefficient intrinsic operators : Performing an asymptotic expansion when δ tends to zero, we infer the asymptotic behavior of the skin depth function (B.2).

∆ + 1 δ ∂ r + ω 2 c 2 (1 + iν)I = δ -2 N -1 n=0 δ n A n + δ N R N,
Comparison with the skin effect in electromagnetism

Figure 1 -Notation 2 . 1

 121 Figure 1 -A cross-section of the domain Ω and the subdomains Ω -and Ω +

Theorem 4 . 1

 41 Under Hypothesis 2.2-2.3, the solution u δ of problem (2.1) has a two-scale expansion which can be written in the form (4.1)-(4.2), with u

For each value of δ = 1 α 1 ( 6 Figure 1 :

 1161 Figure 1: Modulus of the reference solution and solutions using atmosphere boundary conditions.

= 1

 1 ). In Fig 1aredisplayed the obtained solutions. The relative L 2 errors with the reference solution on the inside interval [0, 1] as δ tends to zero are displayed in Fig.2. For the low frequency case (f sun 0 = 0.1 mHz), a preasymptotic regime is observed with a convergence in O(δ 3 ) for Atmo1 and Atmo2 and O(δ 5 ) for Atmo3, before observing the theoretical orders of convergence. For the high frequency case (f sun 0 mHz), the preasymptotic regime extends to the explored range of δ.

Figure 2 :

 2 Figure 2: Relative L 2 errors between the reference solution and solutions using atmosphere boundary conditions.

Figure 3 :

 3 Figure 3: Exact value and approximated value of ρ sun versus the relative radius r/R .

Figure 4 := 7000 and f sun 0 = 3

 403 Figure 4: Real part of u/ √ ρ in the plane Oxz for 1 δ = 7000 and f sun 0 = 3 mHz and the parameters of the Sun.

  Relative L 2 error between exact solution and approximate solution obtained with different boundary conditions versus δ for parameters of the Sun and f sun 0 = 3 mHz.

  Relative L 2 error between exact solution and approximate solution obtained with different boundary conditions versus f sun 0 for parameters of the Sun and δ = 1/7000.

Figure 5 :

 5 Figure 5: Convergence in the realistic context of the Sun.

2 φ 2 .

 22 δ (∂ S , ∂ θ , ∂ φ ; S, θ, φ) for all N ∈ N.The remainder R N,δ has smooth coefficients in (S, θ, φ) which are bounded in δ. Since the Laplace operator ∆ writes in spherical coordinates (r, θ, φ) as∆ = 1 r 2 ∂ r (r 2 ∂ r ) + 1 r 2 sin θ ∂ θ (sin θ∂ θ ) + 1 r 2 sin 2 θ ∂ Inria with m(θ, φ; s) = 1 -2 r t s + O (δ + s)This formula (B.3) is a consequence of equations (4.4) and (4.7). According to Definition B.1, the skin depth function L(δ; .) satisfies |u δ (θ, φ; L(δ; θ, φ))| = |u δ (θ, φ; 0)| e -1 . (B.4)Hence, m(.; L(δ; .)) e -L(δ; .)/δ = e -1 .

According to [ 4 , 2 ωµ0σ

 42 Th. 3.2] the skin depth function L 0 (σ; y α ) associated with the skin effect phenomenon has the following behavior for large conductivity σ,L 0 (σ; y α ) ≈ (σ) 1 -H(y α ) (σ) + O(σ -1 ) , σ → ∞ . (B.5)Here y α , α = 1, 2 are tangential coordinates on Σ, H is the mean curvature of the surface Σ and(σ) = is the classical skin depth parameter which corresponds to δ in this paper. Hence, the δ-asymptotic behavior of the skin depth function (B.2) coincides only up to the first order of approximation with the behavior of the skin depth function L 0 (σ; y α ) characterizing the skin effect phenomenon since the mean curvature of the sphere Σ = {r = r t } is H = 1 rt . Similarly the equivalent boundary conditions associated with these two boundary layer phenomena coincide only up to the second order of approximation, (2.9)-(2.10).RR n°8954RESEARCH CENTRE BORDEAUX -SUD-OUEST200 avenue de la Vieille Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr

RR n°8954

Chabassier & Duruflé & Péron

Chabassier & Duruflé & PéronInria

we observe that

where ∆ Σ is the Laplace-Beltrami operator along the surface Σ, see (2.13). Similarly the normal derivative writes

on the interface Σ. Define v δ by v δ (θ, φ, S) = u + δ (x) where x = X(r t θ, r t φ) + sn (X(r t θ, r t φ)) , and X(r t θ, r t φ) is defined by (4.5).

Remark A.1 The coordinate system (r t θ, r t φ, s) with s = r -r t is a "normal coordinate system" to the surface Σ on the manifold Ω + [START_REF] Caloz | On the influence of the geometry on skin effect in electromagnetism[END_REF].

Then, after the scaling s → S = s δ (where s = r -r t ), the problem (2.1) writes

Moreover v δ tends to 0 when S tends to ∞ since u + δ tends to 0 when r tends to ∞.

A.2 Equations for the terms U n and u - n

We insert the Ansatz

3) and (A.4), and we perform the identification of terms with the same power in δ. Then the terms U n and u - n satisfy the following family of problems coupled by their conditions on the interface Σ (corresponding to S = 0):

In (A.5), we use the convention u - -1 = 0 and in (A.6) δ 0 n denotes the Kronecker symbol. Hereafter, we make explicit the first asymptotics U n and u - n for n = 0, 1, 2, 3 by induction.

A.3 First terms of the asymptotics

The unique solution of (A.7) such that U 0 → 0 when S → ∞ is

Hence, according to (A.5) for n = 0, u - 0 solves the problem

The next term which is determined in the asymptotics is the profile U 1 . According to (A.5) for n = 1 and since

(A.9)

The unique solution of (A.9) such that

Remind that X(r t θ, r t φ) ∈ Σ is defined by (4.5). Hence, according to (A.6

The next term which is determined is the profile U 2 . According to (A.5) for n = 2 and since U 0 = 0, U 2 solves the ODE

(A.12)

According to (A.10), the unique solution of (A.12) such that

where

Hence, according to (A.6) for n = 2, u + 2 solves the problem (see (4.8))

The next term which is determined is the profile U 3 . According to (A.5) for n = 3 and since U 0 = 0, U 3 solves the ODE for S ∈ (0, +∞)

where the operator A 2 is given by (A.2c). According to (A.10)-(A.13), the right hand side in (A.15) writes

since (remind that the operator ∆ Σ is given by (2.13))

Then the unique solution of (A.15) such that U 3 → 0 when S → ∞ is (see (4.9))

where the functions a 3 , b 3 , c 3 are given by

Then, according to (A.6) for n = 3, u - 3 solves the problem (see (4.10))

(A.17)

A.4 Existence and regularity of the asymptotics

We claim without proving the next proposition which ensures existence and regularity results for the asymptotics.

Then for all n ∈ {0, 1, • • • , k + 1} there exists asymptotics U n and u - n satisfying

• U 0 = 0 and (A.8), when n = 0,

• (A.10) and (A.11), when n = 1,

• (A.13) and (A.14), when n = 2,

• (A. [START_REF] Péron | Modélisation mathématique de phénomènes électromagnétiques dans des matériaux à fort contraste[END_REF]) and (A.17), when n = 3,

• (A.5) and (A.6), when n 4.

There holds

The regularity result for the asymptotics is obtained as a consequence of classical elliptic regularity results in Sobolev spaces.

B A measure of the boundary layer phenomenon. Comparison with the skin effect in electromagnetism

As another application of the multiscale expansion, we make the link in this section between the parameter δ and the "pointwise energy" of the solution of the problem (2.1) across the boundary layer. We measure this boundary layer phenomenon by introducing a "skin depth" function that turns out to depend on the mean curvature of the boundary of the domain Ω -. An asymptotic expansion when δ → 0 for this function shows the influence of the geometry of the interface : the skin depth is larger for small values of δ when the mean curvature of the interface Σ is smaller. We refer also to the works [START_REF] Dauge | Comportement asymptotique à haute conductivité de l'épaisseur de peau en électromagnétisme[END_REF][START_REF] Caloz | On the influence of the geometry on skin effect in electromagnetism[END_REF] where the authors measure similarly the skin effect phenomenon by introducing a suitable skin depth function.

Introduction of a skin depth function

For a data f of the problem (2.1), let us define u δ (θ, φ, s) = u + δ (x) for all s 0 and (θ, φ) ∈ Σ (i.e. x ∈ Ω + ). Definition B.1 Let f be a smooth data of problem (2.1) such that for all (θ, φ) ∈ Σ, u δ (θ, φ, 0) = 0. The skin depth is the length L(δ; θ, φ) defined on Σ and taking the smallest value such that

Thus this length L(δ; θ, φ) is the distance from the interface where the field has decreased of a fixed rate e. It depends on δ, of each point (θ, φ) ∈ Σ and a priori on the data f . However it is possible to exhibit the asymptotic behavior of the skin depth function L for small values of δ independently of both (θ, φ) and the data f .

Asymptotic behavior of the skin depth function when δ → 0 Theorem B.1 Let f be a smooth data of problem (2.1). We assume that ∂ n u - 0 X(r t θ, r t φ) = 0 on Σ. Then the skin depth function has the following behavior for small values of δ L(δ; θ, φ) ≈ δ 1 -


Dora Karvouniari

Lionel Gil

Olivier Marre

Serge Picaud

Bruno Cessac

Classifying the spatiotemporal patterns within stage II retinal waves through dynamical systems analysis

Classifying the spatiotemporal patterns within stage II retinal waves through dynamical systems analysis.

Dora Karvouniari, 1 Lionel Gil, 2 Olivier Marre, 3 Serge Picaud, 3 Bruno Cessac 1

(1) Biovision team, Universite Cote d'Azur, INRIA, France, (2) INLN, Sophia-Antipolis, France, (3) Vision Institute, Paris, France theodora.karvouniari@inria.fr; bruno.cessac@inria.fr

Abstract

Retinal waves are bursts of activity occurring spontaneously in the developing retina of vertebrate species, contributing to the shaping of the visual system organization. They are characterized by localized groups of neurons becoming simultaneously active, initiated at random points. Based on our previous modelling work [START_REF] Karvouniari | Modeling the emergence of stage II retinal waves in immature retina[END_REF], we now propose a classification of stage II retinal waves patterns as a function of acetylcholine coupling strength and a possible mechanism for waves generation. Our model predicts that spatiotemporal patterns evolve upon maturation or pharmacological manipulation and that there is a regime of cholinergic coupling, only for which, waves are characterized by power-law distributions.

Context & Motivation

Necessary components for the emergence of retinal waves [START_REF] Zheng | A transient network of intrinsically bursting starbust cells underlies the generation of retinal waves[END_REF] • Spontaneous rhythmic bursting activity of isolated SACs • Refractory mechanism modulating the silent period of the bursting activity (slow After Hy-perPolarisation current, sAHP) • Coupling via cholinergic synapses to ensure the necessary level of synchrony

Isolated sacs

Proposed biophysical process

SACs are in a regime where they can oscillate spontaneously. As they oscillate, the calcium load increases, so the effect of sAHP increases up to a point where oscillations stop, reaching a steady state where the level of the voltage is quite lower. Then, intracellular calcium concentration unloads, I sAHP decreases, until the effect of sAHP is small and oscillations start again.

Bifurcation Analysis

3D Bifurcation diagram of the mechanism generating rhythmic intrinsic bursts in SACs

Network interactions

• SACs are placed on a square lattice connected to their four nearest neighbours. • Synapses are excitatory through cholinergic transmission and their characteristic refractory time follows a Gaussian distribution [START_REF] Ford | Cellular mechanisms underlying spatiotemporal features of cholinergic retinal waves[END_REF].

Synchrony

Two-neurons example

• Introducing a variability in the bursting period, we end up with a network of coupled bursters. • By increasing the coupling, Arnold tongues appear for two neurons. • We expect therefore to see roughly 3 regimes: weak coupling (disordered phase), moderate coupling (waves), strong coupling (full synchrony).

Network Simulations

Simulated local Voltage of 625 neurons on a square lattice. Blue and red colours correspond to low and high activity respectively.

• Weak coupling: localised bumps of activity.

Strong coupling: complete synchrony standing waves. Moderate coupling: propagating patterns.

• Pairwise correlations with respect to distance

show an intermediate regime where we observe anticorrelations. Anticorrelations correspond to a region of hyperpolarization linked to the boundary of the wave.

Calcium Waves

Weak Moderate Strong

Simulated Calcium waves of 4096 neurons on a square lattice.

Black and white colours correspond to low and high activity respectively.

Characterizing Waves

Average Population Bursting Rate • For region B, in all cases, we observe powerlaw distributions of wave characteristics (duration, size).

• This type of distribution, has been primarily observed in [START_REF] Hennig | Early-Stage Waves in the Retinal Network Emerge Close to a Critical State Transition between Local and Global Functional Connectivity[END_REF] in experimental data of stage II retinal waves for different species.

Conclusion

• Cholinergic coupling evolves upon maturation. According to our model, spatiotemporal dynamics depends on Ach coupling. This suggests that, within stage II, there may exist different periods of spatio temporal activity, with possible consequences on visual system development. • In our model, there is a regime of cholinergic coupling where waves are power-law distributed. • Biophysical parameters (e.g. conductances) could vary upon maturation or pharmacological manipulations, affecting the characterstics of emerging waves.

  Distribution of waves duration (top) and size (bottom) for weak, moderate and strong coupling. • Low coupling: exponential distribution. Moderate coupling: power-law distributions. Strong coupling: large waves. Distribution of waves duration and size moderate coupling (region B).

Acknowledgements

This work was supported by the French ministry of Research and University of Nice (EDSTIC). We warmly acknowledge Matthias Hennig and Evelyne Sernagor for their invaluable help.