Philippe Marchner
email: philippe.marchner@siemens.com

Xavier Antoine
email: xavier.antoine@univ-lorraine.fr

Christophe Geuzaine
email: cgeuzaine@uliege.be

Hadrien Bériot
email: hadrien.beriot@siemens.com

Construction and Numerical Assessment of Local Absorbing Boundary Conditions for Heterogeneous Time-Harmonic Acoustic Problems

Keywords: Dirichlet-to-Neumann operator, pseudo-differential calculus, local absorbing boundary conditions, heterogeneous time-harmonic wave propagation, finite element method

come   L'archive ouverte pluridisciplinaire

Introduction

The design of accurate spatial truncation techniques for solving wave propagation problems in unbounded domains is an active area of research, as the introduction of a fictitious boundary is mandatory for numerical methods like finite elements or finite volumes that rely on a volume discretization of the problem. The exact non-reflecting operator at the fictitious boundary is the well-known Dirichlet-to-Neumann (DtN) operator, which is unfortunately non-local and thus difficult and costly to implement in a numerical scheme [START_REF] Hasgtrom | Radiation boundary conditions for the numerical simulation of waves[END_REF]2,[START_REF] Tsynkov | Numerical solution of problems on unbounded domains. A review[END_REF]. Alternatives include Absorbing or Artificial Boundary Conditions (ABC) [START_REF] Hasgtrom | Radiation boundary conditions for the numerical simulation of waves[END_REF][START_REF] Tsynkov | Numerical solution of problems on unbounded domains. A review[END_REF][START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF][START_REF] Engquist | Radiation boundary conditions for acoustic and elastic wave calculations[END_REF][START_REF] Higdon | Radiation boundary conditions for dispersive waves[END_REF][START_REF] Bayliss | Radiation boundary conditions for wave-like equations[END_REF][START_REF] Givoli | High-order local non-reflecting boundary conditions: a review[END_REF], Perfectly Matched Layers (PML) [9,[START_REF] Johnson | Notes on Perfectly Matched Layers (PMLs)[END_REF][START_REF] Turkel | Absorbing PML boundary layers for wave-like equations[END_REF][START_REF] Bermúdez | Perfectly matched layers for time-harmonic second order elliptic problems[END_REF][START_REF] Bériot | An automatic perfectly matched layer for acoustic finite element simulations in convex domains of general shape[END_REF] and infinite elements [START_REF] Astley | Infinite elements for wave problems: a review of current formulations and an assessment of accuracy[END_REF][START_REF] Hohage | Hardy space infinite elements for scattering and resonance problems[END_REF].

Regarding local ABCs built as approximations of the DtN map, tremendous progress has been made since the pioneering work of Engquist and Madja [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF], fostered by the need for ABCs in a broader range of geometrical and physical configurations. In particular, for time-harmonic wave problems governed by the Helmholtz equation, very effective and numerically cheap ABCs are now available for homogeneous wave problems in many situations (see e.g. [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF][START_REF] Turkel | Boundary conditions and iterative schemes for the Helmholtz equation in unbounded regions[END_REF][START_REF] Kechroud | Numerical accuracy of a Padé-type non-reflecting boundary condition for the finite element solution of acoustic scattering problems at high-frequency[END_REF][START_REF] Hagstrom | A new auxiliary variable formulation of high-order local radiation boundary conditions: corner compatibility conditions and extensions to first-order systems[END_REF]20]). Fewer works have tackled the problem of building ABCs for heterogeneous media since the wave propagation problem is more complex and analytic/exact solutions are not available for general situations.

Among the works related to heterogeneous problems, the early study of Engquist and Madja [START_REF] Engquist | Radiation boundary conditions for acoustic and elastic wave calculations[END_REF] compares the "frozen" with the "variable" coefficients technique. While the former uses local values of the heterogeneous field, the latter explicitly incorporates its spatial variation to enhance the quality of the ABC. Their analytic case is revisited in this work with a high-order finite element scheme. Approximations of the DtN map are built in [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF][START_REF] Barucq | Atmospheric radiation boundary conditions for the Helmholtz equation[END_REF]23] for heterogeneous and non-linear [START_REF] Shevchenko | Absorbing boundary conditions for nonlinear acoustics: The Westervelt equation[END_REF] problems following the rules of pseudo-differential calculus. They have in common that accounting for the spatial variation of the heterogeneity in the ABC design increases its accuracy. Stability and well-posedness of such conditions are studied in [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF]25], but this topic is not addressed in the present study. In the work of Hagstrom et al. [26, section 6], an extension to the complete radiation condition [START_REF] Hagstrom | A new auxiliary variable formulation of high-order local radiation boundary conditions: corner compatibility conditions and extensions to first-order systems[END_REF][START_REF] Hagstrom | Complete radiation boundary conditions: minimizing the long time error growth of local methods[END_REF] is adapted for layered, stratified and continuous transverse heterogeneous media where the ABC tuning coefficients associated to the auxiliary equations depends on the local value of the heterogeneity. Alternative techniques combine the ABC with the PML methods [START_REF] Savadatti | Accurate absorbing boundary conditions for anisotropic elastic media. part 1: Elliptic anisotropy[END_REF][START_REF] Hagstrom | The double absorbing boundary method[END_REF]30] and have good potential to be effective in heterogeneous situations. More general methods include root-finding algorithm [START_REF] Lee | Root-finding absorbing boundary condition for poroelastic wave propagation in infinite media[END_REF] or the construction of the ABC at the discrete level [START_REF] Bélanger-Rioux | Compressed absorbing boundary conditions via matrix probing[END_REF]. A rigorous comparison between these methods in heterogeneous situations is challenging.

In this paper we focus on ABCs for heterogeneous time-harmonic acoustic problems, motivated by industrial applications linked to sound radiation in aeroengines. Such problems can be described by a scalar wave equation in which acoustic waves are convected by a potential, isentropic, compressible mean flow [START_REF] Gabard | Adaptive, high-order finite-element method for convected acoustics[END_REF]. The propagation is influenced by a variable speed of sound c 0 (x), density ρ 0 (x) and mean flow velocity v 0 (x). In addition, the work is ultimately motivated by the development of efficient and converging optimized Schwarz Domain Decomposition Methods (DDM) [START_REF] Dolean | An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation[END_REF][START_REF] Gander | Optimized Schwarz methods[END_REF][START_REF] Gander | A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods[END_REF][START_REF] Gander | Optimized Schwarz methods without overlap for the Helmholtz equation[END_REF]38,[START_REF] Boubendir | A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation[END_REF][START_REF] Vion | Double sweep preconditioner for optimized Schwarz methods applied to the Helmholtz problem[END_REF]. Indeed, in this class of hybrid solvers, it is well-known that accurate ABCs used as transmitting boundary conditions at the interfaces between the smaller sub-domains act as preconditioners for the iterative part of the DDM [START_REF] Gander | Optimized Schwarz methods[END_REF][START_REF] Gander | A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods[END_REF]38]. This is relatively well-understood for the homogeneous Helmholtz equation but needs further developments for the heterogeneous Helmholtz equation, and in particular for large scale flow acoustics simulations for industrial applications [START_REF] Lieu | A non-overlapping Schwarz domain decomposition method with high-order finite elements for flow acoustics[END_REF].

To design ABCs for the heterogeneous Helmholtz equation we adopt the Engquist and Majda methodology to build families of local ABCs that can be easily implemented in a high-order finite element solver. The strength of the approach is that an asymptotic expansion of the DtN operator for any PDE with smooth variable coefficients can be explicitly computed thanks to microlocal calculus rules for pseudodifferential operators [START_REF] Taylor | Pseudodifferential operators[END_REF]. This implies that increasing order approximations of the DtN operator can be obtained, and accurately computed locally thanks to a second level of approximation related to Taylor or Padé expansions. The method is general and applicable to arbitrary convex boundaries [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF]. However, a detailed numerical analysis is required to understand whether the method indeed leads to an accurate truncation technique, or when it can fail.

The paper is organized as follows. Section 2 introduces the generic procedure for the derivation of the DtN symbolic asymptotic expansion applied to the half-space heterogeneous Helmholtz operator. We next analyze in detail three situations of interest: the propagation in a duct with a longitudinal variation of the speed of sound in Section 3, the propagation in a non-uniform mean flow using a convected wave operator in Section 4, and the propagation in a duct with a transverse variation of the speed of sound in Section 5. For each case, we develop the construction of families of local ABCs following the process described in Section 2 and explain their implementation in a high-order finite element context. Numerical simulations allow to validate the accuracy of the ABCs, and to give recommendations for the choice of their parameters. Finally, conclusions and perspectives are given in Section 6.

DtN operator for the heterogeneous Helmholtz problem

Let us define the two-dimensional half-space Ω = {x = (x, y) ∈ R 2 ; x ≤ L} with straight vertical fictitious boundary Γ = {x = L}, L > 0. We then look for a formal derivation of the DtN operator Λ + at Γ for the Helmholtz problem in Ω with a variable density ρ 0 (x) = ρ 0 (x, y) and speed of sound c 0 (x) = c 0 (x, y). More precisely, the heterogeneous Helmholtz equation (for the e iωt convention with time frequency ω > 0) is associated to the partial differential operator

H(x, y, ∂ x , ∂ y , ω) = ρ -1 0 ∂ x (ρ 0 ∂ x ) + ρ -1 0 ∂ y (ρ 0 ∂ y ) + ω 2 c -2 0 . (2.1)
The operator Λ + is then defined as

Λ + : H 1/2 (Γ) → H -1/2 (Γ) u |Γ -→ ∂ n u |Γ = -i Λ + u |Γ , (2.2)
where n = (1, 0) is the outwardly directed unit normal vector to the straight boundary Γ, which means that ∂ n = ∂ x . All the developments that follow could readily be extended to the 3D case, but would require further algebraic computations leading to more complicated formulas to analyze. The analysis could also be extended to the case of a curved convex boundary based on the tangent plane approximation [START_REF] Antoine | An improved surface radiation condition for high-frequency acoustic scattering problems[END_REF]. We formally compute the DtN operator in the framework of microlocal analysis by using the tools of pseudo-differential operator calculus [START_REF] Taylor | Pseudodifferential operators[END_REF], as initiated by Engquist and Majda [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF] for hyperbolic systems. To this end, we consider that the density ρ 0 and the speed of sound c 0 are smooth functions of the spatial variable x. In addition, since we will work with classical pseudo-differential operators, it is well-known that the analysis is not valid for grazing waves and would require the introduction of new classes of operators (like e.g. Gevrey classes [START_REF] Lafitte | The kernel of the Neumann operator for a strictly diffractive analytic problem[END_REF]). This specific point will be clarified later when we enter into the details of the concrete cases.

Preliminary definitions

A partial differential operator may be defined as [START_REF] Taylor | Pseudodifferential operators[END_REF] P(x,

∂ x ) = |α|≤m (-i) α a α (x)∂ α x , (2.3)
with x ∈ R d (d = 2 here), α = (α 1 , . . . , α d ) ∈ N d a multi-index, |α| = d j=1
α j and m the order of the operator P. We define its symbol p to be the polynomial

p(x, ξ) = |α|≤m a α (x)ξ α , (2.4) setting ξ ∈ R d and ξ α = ξ α1 1 . . . ξ α d d . If we introduce the Fourier transform of a smooth function u of compact support in R d û(ξ) = R d u(x)e -ix•ξ dx, (2.5)
where the property (∂ α x u)(ξ) = (-iξ) α û(ξ) can be justified, a partial differential operator can be seen as the inverse Fourier representation of its symbol

P(x, ∂ x )u(x) = 1 (2π) d R d e ix•ξ p(x, ξ)û(ξ)dξ.
(2.6)

We next introduce the principal symbol σ(P) of an operator P , which is the highest order homogeneous term in ξ

σ(P)(x, ξ) = |α|=m a α (x)ξ α . (2.7)
This point of view allows to work with a wide range of operators, where the symbol is not necessarily a polynomial, but a smooth function p = p(x, ξ). This leads to the broad theory of pseudo-differential operators [START_REF] Taylor | Pseudodifferential operators[END_REF]. We will denote by P = Op(p) ∈ OPS m the pseudo-differential operator of order m associated to the symbol p. Roughly speaking, the symbol should have a polynomial growth that equals the order of the operator when |ξ| → ∞. Sets of admissible symbols were introduced by Hörmander [START_REF] Hörmander | The analysis of linear partial differential operators III: Pseudo-differential operators[END_REF] and we will use the set of classical symbols S m cl , that are defined by the estimate

∂ β x ∂ α ξ p(x, ξ) ≤ C(1 + |ξ|) m-|α| , ∀(x, ξ) ∈ K × R d , for each pair of multi-indices α, β, a compact set K of R d and a constant C = C(α, β, K) ∈ R.
We make the correspondence p ∈ S m cl ⇔ P ∈ OPS m . Finally, we define OPS -∞ = m∈R OPS m to be the set of regularizing operators. This corresponds to operators with a smooth integral kernel, for which the symbol decays faster than any polynomial. In practice, algebraic computations can be performed at the symbol level to construct approximations of increasing orders of the DtN operator. We only give the necessary formula for the computations and refer to [START_REF] Taylor | Pseudodifferential operators[END_REF] for more details on the underlying theory.

Symbols computation

The first step to derive the DtN operator for the heterogeneous Helmholtz problem consists in splitting the Helmholtz operator into two operators that characterize the forward and backward parts of the wave field u. According to Nirenberg's factorisation theorem [START_REF] Nirenberg | Pseudodifferential operators and some applications[END_REF], there exist two pseudo-differential operators Λ ± = Λ ± (x, y, ∂ y ) of order +1 such that

H = (∂ x + iΛ -)(∂ x + iΛ + ) mod OPS -∞ , (2.8)
where the equality holds modulo a regularizing pseudo-differential operator. This reformulation enlightens the role of the DtN operator, which is the trace of the outgoing (or ingoing) characteristic of the wave on the boundary Γ. Equation (2.8) can be seen as a reformulation of the Helmholtz equation as two "one-way" wave equations. In a microlocal sense, the solution to the Helmholtz problem is given by the operators (Λ + , Λ -).

For the constant coefficients case and a plane wave of the form e i(ωt-kxx) , one could directly infer and verify that

Λ ± = ±k x , k x = ω 2 c -2 0 + ∂ 2 y . (2.9)
This is a square-root, pseudo-differential operator of order +1. However, we will see that it is only an approximation in the heterogeneous case. Rather than looking for the DtN map directly in its operator form, a more suitable approach is to look for its symbol. We turn back to the Helmholtz operator and develop the factorisation (2.8) as

H = ∂ 2 x + i∂ x Λ + + iΛ -∂ x -Λ -Λ + mod OPS -∞ . (2.10)
From definition (2.6), we obtain for the wave field u

∂ x (Λ + u) = (2π) -1 ∂ x R e iyξ λ + û dξ = (2π) -1 R e iyξ ∂ x λ + û + λ + ∂ x û dξ (2.11) = Op ∂ x λ + u + Λ + ∂ x u, (2.12)
where λ + is the symbol of Λ + . The Helmholtz operator can be recast as

H = ∂ 2 x + i Λ + + Λ -∂ x + iOp ∂ x λ + -Λ -Λ + mod OPS -∞ . (2.13)
With this form, one is able to identify in (2.1) and (2.13) the first and zeroth order x-derivatives, which leads to the system

Λ + + Λ -= -iρ -1 0 ∂ x (ρ 0 ) -Λ -Λ + + iOp {∂ x λ + } = ω 2 c -2 0 + ρ -1 0 ∂ y (ρ 0 ∂ y ) , (2.14)
from which we can eliminate Λ -and obtain an equation for the outgoing operator Λ +

(Λ + ) 2 + iρ -1 0 ∂ x (ρ 0 )Λ + + iOp ∂ x λ + = ω 2 c -2 0 + ρ -1 0 ∂ y (ρ 0 ∂ y ). (2.15)
This equation is still a reformulation of the initial problem. The simplification occurs thanks to the asymptotic expansion property of classical symbols. Since Λ + is of order +1, the asymptotic expansion writes

λ + ∼ ∞ j=-1 λ + -j = λ + 1 + λ + 0 + λ + -1 + • • •
where each symbol λ + -j is homogeneous of order -j in (ω/c 0 , ξ), and the equivalence class ∼ has the meaning

∀m ≥ -1, λ + - m j=-1 λ + -j ∈ S - (m+1) cl
.

The asymptotic expansion gives a notion of "convergence" for pseudo-differential operators and allows to compute successive approximations of their symbol with respect to their homogeneity order. If we express (2.15) in its symbolic form, one obtains

(λ + ) 2 + iρ -1 0 ∂ x (ρ 0 )λ + + i∂ x λ + = ω 2 c -2 0 -ξ 2 -iξρ -1 0 ∂ y (ρ 0 ). (2.16)
There are well-defined calculus rules for classical symbols allowing them to be obtained in a recursive manner [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF]. For instance, the composition rule for pseudo-differential operators [START_REF] Taylor | Pseudodifferential operators[END_REF] can be used to compute the square of the symbol λ + , and further sort the terms by their decaying homogeneity degree

(λ + ) 2 ∼ α≥0 (-i) α α! ∂ α ξ λ + ∂ α y λ + = (λ + 1 ) 2 + 2λ + 0 λ + 1 -i∂ ξ λ + 1 ∂ y λ + 1 mod S 0 cl .
(2.17)

The identification of the second-order homogeneous terms in (2.16) yields a choice for the principal symbol

λ + 1 = ω 2 c -2 0 -ξ 2 , (2.18)
which is directly linked to the operator (2.9) in the constant coefficients case. Alternatively, the principal symbol could be defined in a inhomogeneous manner by identifying the second order terms in the left-hand side of (2.16) only

λ + 1 = ω 2 c -2 0 -ξ 2 -iξρ -1 0 ∂ y (ρ 0 ), (2.19)
from which one may associate the operator

Λ + 1 = ω 2 c -2 0 + ρ -1 0 ∂ y (ρ 0 ∂ y ). (2.20)
The golden rule is to ensure that Re(λ + 1 ) ≥ 0 in the high frequency limit ω → +∞, such that the wave is outgoing. Although both choices could be valid, we will use expression (2.18) to be the principal symbol. Once it is fixed, the lower order symbols are uniquely determined. For a transverse heterogeneous medium (no x-dependence), Nirenberg's factorization theorem suggests that expression (2.20) is the most appropriate operator to represent the DtN map. This will be emphasized and confirmed in Section 5. For the choice of λ + 1 (2.18) we identify the first-order homogeneous terms from relation (2.16) and get the zeroth order symbol

λ + 0 satisfying 2λ + 0 λ + 1 -i∂ ξ λ + 1 ∂ y λ + 1 + iρ -1 0 ∂ x (ρ 0 )λ + 1 + i∂ x λ + 1 = -iξρ -1 0 ∂ y (ρ 0 ). (2.21)
Standard calculus rules lead to

∂ ξ λ + 1 = - ξ λ + 1 , ∂ x λ + 1 = ω 2 ∂ x (c -2 0 ) 2λ + 1 , (2.22)
and similarly for the y-derivative. It results the zeroth order symbol

λ + 0 = -i ∂ x (ρ 0 ) 2ρ 0 + ξ∂ y (ρ 0 ) 2ρ 0 λ + 1 + ω 2 ∂ x (c -2 0 ) 4(λ + 1 ) 2 + ξω 2 ∂ y (c -2 0 ) 4(λ + 1 ) 3 . (2.23)
There are four terms of decreasing powers of λ + 1 . It means that the density variations have the largest impact on the wave behaviour. If Im(λ + 1 ) = 0 holds, we can interpret the zeroth order symbol as an amplitude correction to the wave. In a similar manner, the rest of the symbols can formally be obtained from the composition rule.

Microlocal regimes

The criterion Re(λ + 1 ) ≥ 0 defines the cone of propagation for outgoing waves and characterizes the hyperbolic zone in microlocal analysis [START_REF] Majda | Reflection of singularities at the boundary[END_REF]. It means that only propagative modes are modeled (low-frequency spatial modes). When the frequency is such that ωc -1 0 ≈ |ξ|, the wave is said to be in the grazing zone, and the current approach is a priori not valid. We will see that this situation introduces a singularity for higher order symbols. The last situation ωc -1 0 < |ξ| corresponds to the elliptic zone and describes evanescent modes. One may extend the present analysis to the elliptic regime by selecting the branch-cut of the principal symbol to be along the negative real axis such that

λ + 1 = -i ξ 2 -ω 2 c -2 0 .
(2.24)

Both propagative and evanescent modes can be taken into account if we rewrite the square-root as suggested by Milinazzo et al. [START_REF] Milinazzo | Rational square-root approximations for parabolic equation algorithms[END_REF]. The idea was introduced for parabolic equations in underwater acoustics, and adapted to ABCs for homogeneous Helmholtz problems in [START_REF] Antoine | An improved surface radiation condition for high-frequency acoustic scattering problems[END_REF]. It consists in rotating the branch-cut of the square-root by an angle α

λ + 1 = e iα/2 e -iα (ω 2 c -2 0 -ξ 2 ), (2.25)
such as the hyperbolic zone is fully captured when α = 0 and the elliptic zone when α = -π. The angle α can be seen as a parameter to be selected in order to find a good balance between the modeling of the propagative and evanescent modes.

Summary

We can build an approximate surface DtN operator Λ + M , called DtN M , by i) keeping the M first symbols in the sum

∂ n u = -iΛ + M u, Λ + M = M -2 j=-1 Op(λ + -j ), (2.26)
and then ii) taking the trace on the boundary Γ. In the following, to simplify the notations, we forget the tilde and the plus sign + when considering the trace on the boundary of the outgoing operator. As an example, the surface operator Λ + M is denoted by Λ M . We will now consider some simplified situations and explicitly build local ABCs for three practical situations of increasing difficulty:

• the propagation in a duct with a longitudinal variation of the speed of sound;

• the propagation in a non-uniform mean flow using the linearized potential operator;

• the propagation in a duct with a transverse variation of the speed of sound and density. These three cases allow to numerically analyze the strengths but also the inherent limitations of the developed approach.

Longitudinal heterogeneous problem

We first analyze the method in the case of a longitudinal heterogeneous duct problem. More precisely, we consider the situation of a constant density set to unity (ρ 0 = 1) and a speed of sound varying along the propagation direction, namely c 0 (x, y) = c 0 (x). We specifically revisit the case of a linear profile

c -2 0 (x) = ax + b, (3.1)
since an analytic expression for the DtN operator is available [START_REF] Engquist | Radiation boundary conditions for acoustic and elastic wave calculations[END_REF][START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF]. We always suppose a = 0, otherwise we have the homogeneous case.

Symbols and associated operators

We approximate the DtN symbol by the first two terms of the asymptotic expansion

λ + ≈ λ + 1 + λ + 0 = ω 2 c -2 0 -ξ 2 -i ω 2 ∂ x (c -2 0 ) 4(ω 2 c -2 0 -ξ 2 ) . (3.2)
The next step is to go back at the operator level by inverse Fourier transform. We propose the first-and second-order DtN M operators respectively defined by

Λ 1 = ω 2 c -2 0 + ∆ Γ mod OPS -2 , (3.3)
Λ 2 = ω 2 c -2 0 + ∆ Γ - iω 2 ∂ x (c -2 0 ) 4 (ω 2 c -2 0 + ∆ Γ ) -1 mod OPS -3 , (3.4)
where ∆ Γ = ∂ 2 y is the Laplace-Beltrami operator along the transverse direction y. By construction, the symbols of Λ 1 and Λ 2 are exactly λ + 1 and λ + 1 + λ + 0 , respectively. This is because the symbols are evaluated at x = L and thus do not depend on the transverse variable y, leading to a natural choice for the operators. These operators are still non-local but follow a hierarchic degree of regularity with respect to the exact DtN. Note that the regularity estimate for the second-order condition Λ 2 is based on the next symbol λ

+ -1 λ + -1 = 5ω 4 ∂ x (c -2 0 ) 2 32 λ + 1 5 - ω 2 ∂ 2 x (c -2 0 ) 8 λ + 1 3 , (3.5)
which corresponds to an operator in OPS -5 for a linear profile.

Construction of the ABCs

To derive an ABC that can be easily implemented in a finite element method, we need an additional approximation of the non-local operators DtN M to represent them through local partial differential operators. This can be achieved thanks to Taylor and Padé approximants, with the "angle of incidence" (ξc 0 /ω) as a small parameter.

Let us start with the Taylor expansion. From the first symbols, we can build some so-called local complete radiation boundary conditions as defined in [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF]. For example, we can compute the second-order Taylor expansion (denoted by (•) 2 ) of the first four symbols

Λ 2 4 = Op   2 j=-1 (λ + -j ) 2   , on Γ. (3.6)
With the help of the Python library SymPy [START_REF] Meurer | SymPy: symbolic computing in Python[END_REF], we formally compute the Taylor approximation with respect to (c 0 ξ/ω), leading to the local condition

Λ 2 4 = ω c 0 -i ∂ x (c -2 0 ) 4c -2 0 + 5 ∂ x (c -2 0 ) 2 32ωc -5 0 + i 15 ∂ x (c -2 0 ) 3 64ω 2 c -8 0 (3.7) + c 0 2ω + i ∂ x (c -2 0 ) 4ω 2 c -4 0 - 25 ∂ x (c -2 0 ) 2 64ω 3 c -7 0 -i 15 ∂ x (c -2 0 ) 3 16ω 4 c -10 0 ∆ Γ , (3.8)
where we have set the second-and third-order derivatives of c -2 0 to zero for a linear profile. The Python code used for the formal derivation of Λ 2 4 is given in Appendix B. Note that any lower order condition can easily be recovered by dropping some terms. The above condition is one of the most accurate ABCs that can be directly built with standard (mass and stiffness) finite element matrices, but as we will see, is less accurate than the Padé-based conditions. In the following, it will be denoted by ABC 2 4 . The Padé expansion, which is expected to be more robust than its Taylor counterpart, requires an additional computational cost. We only consider the approximation of the first two symbols, and rewrite the second-order nonlocal operator DtN 2 as

Λ 2 = k 0 √ 1 + X -i ∂ x (c -2 0 ) 4c -2 0 (1 + X) -1 , k 0 = ω/c 0 , X = ∆ Γ /k 2 0 . (3.9)
The usual Padé approximation for the square root

√ 1 + X, X → 0 leads to Λ 2 ≈ k 0 1 + N =1 (a X) (1 + b X) -1 -i ∂ x (c -2 0 ) 4c -2 0 (1 + X) -1 , (3.10)
a = 2 2N + 1 sin 2 π 2N + 1 , b = cos 2 π 2N + 1 . (3.11)
This condition is still non-local because of the inverse operators. We use auxiliary functions to give them a local, hence sparse discrete representation. A total of N auxiliary functions is required for the Padé-type ABC based on Λ 1 since it is defined by the square-root operator, and N + 1 functions for the ABC related to Λ 2 . The implementation is detailed in Section 3.3. Until now, the proposed conditions are only valid in the hyperbolic zone, whenever the square-root term of the principal symbol is strictly positive. Thanks to the branch-cut rotation introduced in (2.25), the square-root operator becomes

e iα/2 √ 1 + z, z = e -iα (1 + X) -1. (3.12)
The Taylor approximants are extended as

√ 1 + z ≈ e iα/2 N =0 1/2 (1 + z)e -iα -1 , (3.13)
and the Padé ones as

√ 1 + z ≈ C 0 + N =1 A z 1 + B z , (3.14)
C 0 = e ı α 2 R N e -ıα -1 , A = e -iα 2 a (1 + b (e -ıα -1)) 2 , B = e -ıα b 1 + b (e -ıα -1) . (3.15)
In the following, the Padé-type ABC with N terms and rotation angle α, based on M symbols, is denoted by ABC N,α M . The grazing zone leads to the situation ξ 2 ≈ ω 2 c -2 0 , which introduces an explicit singularity for higher order symbols. A workaround is to complexify the frequency ω in the denominator of the zeroth order symbol by adding a local damping term ε near the singularity, i.e. we set: ω ε = ω -iε. Thus λ + 0 is modified to be the regularized symbol

λ + 0 = -i ω 2 ∂ x (c -2 0 ) 4(ω 2 ε c -2 0 -ξ 2 ) . (3.16)
In the specific case of a linear profile, it is possible to choose ε such that it minimizes the reflection coefficient at the turning point. The optimal value ε opt is derived in Appendix A and is inspired from the approach used in [START_REF] Antoine | An improved surface radiation condition for high-frequency acoustic scattering problems[END_REF]. In practice, this zone is expected to be limited to isolated frequencies.

Description of the test case

For the numerical test case, we choose a two-dimensional straight duct problem as shown in Figure 1 with a strictly positive, linear speed of sound profile. The domain Ω represents a connected semi-infinite duct that has been truncated by a straight boundary. Since they involve the Laplace-Beltrami operator on Γ, the introduced ABCs could be applied to any smooth convex geometry by considering a tangent plane approximation. However, analytic and/or references test cases for heterogeneous Helmholtz problems are scarce and we therefore focus on this simpler situation.

wave direction For a linear profile, the heterogeneous Helmholtz equation

c -2 0 (x) = ax + b ∂ n u = 0, Γ 2 Γ 1 ∂ n u = g Γ 3 ∂ n u = -iΛu 0 Ω = [0, L] × [0, H] H L x y
∂ 2 x u + ∂ 2 y u + k 2 0 u = 0, k 0 = ω/c 0 , c -2 0 (x) = ax + b, (3.17)
can be explicitly solved by separation of variables. The exact outgoing solution u ex (with the e iωt convention) for a given mode n and a strictly positive increasing speed of sound profile is given thanks to Airy's function [START_REF] Lévy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF] as

u ex (x, y) = cos(k y y)Ai e -2iπ 3 k 2 y -k 2 0 (aω 2 ) 2/3 , k y = nπ H , n ∈ N. (3.18)
Note that the choice of a decreasing profile is also possible but is more challenging to set up numerically (see [START_REF] Lévy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF], pp. 122-125). The x-derivative of the exact solution is

∂ x u ex (x, y) = -e -2iπ 3 (aω 2 ) 1/3 cos(k y y)Ai' e -2iπ 3 k 2 y -k 2 0 (aω 2 ) 2/3 , (3.19)
such that the exact DtN operator on Γ 3 is given by

Λ = -ie -2iπ 3 (aω 2 ) 1/3 Ai'(z) Ai(z) , z = e -2iπ 3 k 2 y -ω 2 (aL + b) (aω 2 ) 2/3 , (3.20)
with L being the x-position of the fictitious boundary. Note the special situation whenever

k 2 y = ω 2 (ax + b), (3.21)
which we may solve for x or ω. When this specific situation occurs, the principal symbol cancels. The nature of the wave changes and shows a transition from cut-off to cut-on (or vice versa). The x-location where equation (3.21) holds is called a turning point [START_REF] Rienstra | An Introduction to Acoustics[END_REF] and denoted by x t . Next, we write the variational formulation of the boundary value problem for the exact DtN operator:

find u ∈ H 1 (Ω) such that ∀v ∈ H 1 (Ω), Ω {∇u • ∇v -k 2 0 u v} dΩ + i Γ3 Λu v dΓ 3 = Γ1 g v dΓ 1 , (3.22)
and g is given by the normal derivative of the exact solution at x = 0

g = e -2iπ 3 (aω 2 ) 1/3 cos(k y y)Ai' e -2iπ 3 k 2 y -bω 2 (aω 2 ) 2/3 . (3.23)
We now present how to implement the Padé-type ABCs in a finite element context. The auxiliary variables related to the principal and the zeroth order symbols are respectively denoted by (ϕ 1 , • • • , ϕ N ) and ψ. The weak form of ABC N,α 2 given by (3.9) together with (3.14)-(3.15) leads to the coupled system: for

(u, ϕ 1 , • • • , ϕ N , ψ) ∈ H 1 (Γ 3 ) × H 1 (Γ 3 ) N × H 1 (Γ 3 ) i Γ3 ABC N,α 2 u v dΓ 3 = i Γ3 k 0 C 0 u v dΓ 3 -i Γ3 A k 0 ∇ Γ3 ϕ • ∇ Γ3 v dΓ 3 + Γ3 ∂ x (c -2 0 ) 4c -2 0 ψ v dΓ 3 , (3.24) ∀v ∈ H 1 (Γ 3 ), Γ3 k 2 0 ϕ v dΓ 3 -B Γ3 ∇ Γ3 ϕ • ∇ Γ3 v dΓ 3 = Γ3 k 2 0 u v dΓ 3 , (3.25) ∀µ ∈ H 1 (Γ 3 ), Γ3 k 2 0 ψ µ dΓ 3 - Γ3 ∇ Γ3 ψ • ∇ Γ3 µ dΓ 3 = Γ3 k 2 0 u µ dΓ 3 . (3.26)
In its discrete matrix form, the complex-valued sparse linear system for the global problem is similar to (3.22) but includes the ABC N,α

2 terms (3.24)-(3.26) of size [n dof,Ω + (N + 1)n dof,Γ3 ] × [n dof,Ω + (N + 1)n dof,Γ3 ] AU = F , U =        u ϕ 1 . . . ϕ N ψ        , F =        -g 0 . . . 0 0        , (3.27)
where the FEM matrix has a quasi block-diagonal structure

A =          K -k 2 0 M + ik 0 C 0 M Γ3 -i A1 k0 K Γ3 • • • -i A N k0 K Γ3 ∂x(c -2 0 ) 4c -2 0 M Γ3 -k 2 0 M Γ3 k 2 0 M Γ3 -B 1 K Γ3 0 • • • 0 . . . 0 . . . 0 . . . -k 2 0 M Γ3 . . . 0 k 2 0 M Γ3 -B N K Γ3 0 -k 2 0 M Γ3 0 • • • 0 k 2 0 M Γ3 -K Γ3          . (3.28)
One can remove the last row and column to obtain the matrix associated to ABC N,α

1 , which is of size [n dof,Ω + N n dof,Γ3 ] × [n dof,Ω + N n dof,Γ3
]. In the above notations, the number of degrees of freedom in Ω (respectively on Γ 3 ) is denoted by n dof,Ω (respectively n dof,Γ3 ). The mass and stiffness matrices for the volume problem are M and K. For the surface Γ 3 , the mass and stiffness matrices are M Γ3 and K Γ3 , respectively.

Numerical results

Let us now study the numerical behavior of the ABCs and the effect of the various parameters. We use a high-order finite element scheme equipped with a basis of integrated Legendre polynomials [START_REF] Solin | Higher-Order Finite Element Methods[END_REF] to discretize the weak formulations. Such a scheme allows to effectively control the interpolation and dispersion errors associated to Helmholtz problems [START_REF] Bériot | Efficient implementation of high-order finite elements for Helmholtz problems[END_REF], and is well-suited for testing the accuracy limits of the ABCs.

The implementation uses the GmshFEM [START_REF] Royer | Gmsh-Fem: An Efficient Finite Element Library Based On Gmsh[END_REF] finite element library and is available online1 under open-source license. The ABCs effectiveness is measured by the relative L 2 -error (in %) in the domain Ω = [0, L] × [0, H] as

E L 2 = 100 u ex -u h L 2 (Ω) u ex L 2 (Ω) , (3.29)
where u h refers to the discretized solution. As a reference, we consider the numerical solution obtained with the exact DtN operator (3.20), such that the remaining error is only due to the finite element discretization. The speed of sound profile is assumed strictly increasing and we use c -2 0 (x) = 5x + 0.1. The duct is taken of length L = 1 and height H = 0.5. The mesh is generated by Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and postprocessing facilities[END_REF] and is composed of linear quadrangle elements Q4 of size h = 1/40. The p-FEM shape function order is fixed to p = 6, and the integration on the reference element is computed by a tensorised Gauss quadrature rule.

For the first experiment, we set up the single mode n = 3 on the input left boundary Γ 1 . We present in Figure 2 We can distinguish three frequency regimes

• if x t > L, the mode is evanescent;

• if 0 < x t < L, the mode shows an evanescent to propagative transition within the duct. The mode becomes cut-on at x = x t ;

• if x t < 0, the mode is propagative.

We report in Figure 3 the relative L 2 -error for the Taylor and Padé-based ABCs of different orders M as a function of the input frequency ω. The parameters α and ε are for now fixed to zero, and will be gradually turned on throughout the numerical experiments. As expected, adding symbols in the ABCs results in a decrease of the error. It results from a better DtN approximation which takes into account the speed of sound variation ∂ x (c -2 0 ). Regarding Taylor-based conditions, ABC 2 4 reaches a very good accuracy at no additional cost. Moreover, the contribution from the zeroth order symbol has the highest impact on the ABC accuracy. For Padé-based conditions, ABC N,0 1 reaches an error plateau as N grows. A plateau is also reached for the condition ABC N,0 2 , but is approximately two orders of magnitude lower than ABC N,0 1 . This is in contrast to the homogeneous case, where N can be increased up to the reference error precision. For such a heterogeneous case, a higher number of auxiliary fields N adds high-order terms to the ABC that do not match the full DtN operator, but only its truncated version DtN M . This is confirmed by the numerical tests, since each condition reaches an error that is consistent with the order M of the ABC. A proper localization procedure ensures that the maximal efficiency of the ABC is attained. Note that in this situation, the duct is large enough such that the low frequencies associated to evanescent modes are damped regardless of the ABC and the parameter α. ).

In the second experiment, we analyze more precisely the low frequency regime. We fix the frequency ω and vary the position of the ABC, that is the length of the duct L. The mode order is chosen to be n = 5. Thanks to the turning point relation (3.21), we select the frequency such that x t = 0.5, giving ω ≈ 19.5. The mode is evanescent at the input boundary and turns into propagative if the ABC is located at L > x t . The ABC location is varied by adding a single mesh layer to the duct for each computation. The rotation angle α is now turned on.

We focus here on the most accurate condition ABC N,α

2

. In Figure 4a, the influence of the angle α is consistent with its theoretical interpretation: a large angle improves the attenuation of evanescent modes while deteriorating the attenuation of propagative modes. The latter effect can be reduced by taking a large number of auxiliary fields N thus improving the localization of the principal symbol, as shown in Figure 4c. For a sufficiently large N (here N = 6), there is no more gain in the accuracy of ABC N,α 2 and the error plateau corresponding to the expected DtN approximation is attained.

However at L = 0.5, the ABC is located at the turning point and the condition ABC N,α 2 clearly shows a lack of accuracy. This is most likely linked to the inverse operator arising as a corrective term because such a behaviour is not observed for the condition ABC N,α 1 . As a workaround, we turn on the parameter ε in ABC N,α 2 , see Figure 4b. For conciseness, the parameter ε does not appear in the notation ABC N,α 2 but is introduced thanks to ω ε = ω -iε. The optimal value obtained theoretically in Appendix A results in a smaller error at the turning point, and improves the ABC efficiency in the propagative low frequency regime. Another value of ε could be more efficient globally, and more advanced strategies for the optimization procedure could be conducted. The choice for ε is however case dependent, and such an optimal value is in practice difficult to determine a priori.

Note that Taylor-based conditions can also be designed while rotating the branch-cut of the principal symbol. As ABCs, such complexified Taylor conditions are less robust than Padé-based conditions, which has been confirmed by our numerical experiments. Nonetheless, complexified Taylor ABCs might be very effective as a transmission condition for domain decomposition methods [START_REF] Lieu | A non-overlapping Schwarz domain decomposition method with high-order finite elements for flow acoustics[END_REF], and can be contructed numerically at no additional cost. Such conditions will be tested in a domain decomposition framework for heterogeneous problems in the future. ).

In more realistic situations, it is not always possible to predict whether the wave hitting the interface will be propagative or evanescent. Therefore, we recommend to rotate the branch-cut of the principal symbol before the localization procedure. A sufficiently large number of auxiliary fields N should be used in order to ensure a proper localization of the square-root operator, especially for a large angle α. For high frequency applications, we recommend to keep α in the range [0, -π/2] and increase N with α, e.g. α = -π/2 with N = 6.

A stability analysis remains to be conducted for the introduced ABCs. The stability of Padé-based ABC have been addressed for homogeneous problems in [START_REF] Trefethen | Well-posedness of one-way wave equations and absorbing boundary conditions[END_REF], but the precise effect of the rotation angle α remains to be examined.

The convected wave operator

We focus in this section on a second more complex situation, that is modal propagation in a non-uniform mean flow. Here, we only analyze the case of Padé-type ABCs which have been shown to be more accurate than ABCs based on Taylor expansion in the previous section. In two-dimensions, the acoustic field is described by the general convected wave operator in the frequency domain

L(x, ∂ x , ω) = D 0 Dt 1 c 2 0 D 0 Dt -ρ -1 0 ∇ • (ρ 0 ∇), D 0 Dt = iω + v 0 • ∇, v 0 = (v x , v y ), (4.1)
which is commonly used in flow acoustics [START_REF] Howe | Acoustics of fluid-structure interactions[END_REF][START_REF] Rienstra | An Introduction to Acoustics[END_REF][START_REF] Pierce | Wave equation for sound in fluids with unsteady inhomogeneous flow[END_REF][START_REF] Goldstein | Aeroacoustics[END_REF]. It governs the acoustic velocity potential u L(x, ∂ x , ω)u = 0, and the acoustic velocity is v = ∇u. The acoustic pressure variable p can be recovered through the relation

p = -ρ 0 D 0 Dt u. (4.2)
The mean flow quantities v 0 , ρ 0 , c 0 are steady and should satisfy the physical relations of a compressible, irrotational and homentropic flow:

• mass conservation: ∇ • (ρ 0 v 0 ) = 0,
• Bernoulli's equation:

v 2 0 2 + c 2 0 γ-1 = K 1 , γ = 1.
4 being the adiabatic constant in dry air,

• isentropic state equation for a perfect gas:

c 2 0 = K 2 γρ γ-1 0 .
For our purpose of testing ABCs, we would like to extend the straight duct situation of the previous section and consider a non-uniform flow along the x-direction (here the direction of propagation). Unfortunately this situation is only physically relevant for a duct of variable cross-section, and doing so would break the separable modal structure of the solution. Instead, we choose to drop the assumption of the flow to be isentropic. Although being physically incorrect, the model remains relatively simple and is still relevant to derive and test the ABCs. Future work will investigate more realistic situations with two-dimensional flows. For a variable mean flow along the x-direction, the convected wave operator takes the form

L(x, ∂ x , ω) = -ω 2 c -2 0 + 2iv x ωc -2 0 ∂ x + iωv x ∂ x (c -2 0 ) + v x ∂ x (v x c -2 0 ∂ x ) -ρ -1 0 ∇ • (ρ 0 ∇). (4.3)
We use the local wavenumber k 0 = ω/c 0 , the velocity variable M x = v x /c 0 and recast the operator as

L(x, ∂ x , ω) = (M 2 x -1)∂ 2 x -∂ 2 y + 2ik 0 M x -ρ -1 0 ∂ x (ρ 0 ) + v x ∂ x (v x c -2 0 ) ∂ x -k 2 0 + iωv x ∂ x (c -2 0 ). (4.4)
Let us remark that [START_REF] Cummings | Sound generation and transmission in flow ducts with axial temperature gradients[END_REF] presents a closely related one-dimensional model. Compared to the uniform mean flow situation, the non-uniformity adds two first-order terms to the operator and a zeroth order term that couples the frequency and the mean flow. In this model, a decrease in the flow velocity is balanced by a larger speed of sound and density. The latter can be interpreted from the mass conservation relation in one-dimension

∂ x (ρ 0 )ρ -1 0 = -∂ x (v x )v -1 x . (4.5)
We further impose the restrictions of a non-zero subsonic mean flow 0 < |M x | < 1, a positive density ρ 0 > 0 and speed of sound c 0 > 0. Then, it has been shown in e.g. [START_REF] Bonnet-Bendhia | Acoustic diffraction by a plate in a uniform flow[END_REF] that equation (4.1) with a suitable outgoing radiation condition is a problem of Fredholm type.

Symbols computation

As done throughout the document, the objective is to compute the first terms from the asymptotic expansion of the DtN symbol and further construct a hierarchical set of ABCs. Let us consider the normalized operator

L (x, ∂ x , ω) = L(x, ∂ x , ω)/(M 2 x -1), (4.6)
such that one can apply Nirenberg's factorization theorem

L (x, ∂ x , ω) = (∂ x + iΛ -)(∂ x + iΛ + ) mod OPS -∞ , (4.7)
= ∂ 2 x + i Λ + + Λ -∂ x + iOp ∂ x λ + -Λ -Λ + mod OPS -∞ . (4.8)
As in Section 2, the identification of the first and zeroth order derivatives with the normalized convected wave operator leads to a coupled system

i (Λ + + Λ -) = (A 1 + A 0 ) /(M 2 x -1) -Λ -Λ + + iOp {∂ x λ + } = k 2 0 + ∂ 2 y -iB 1 /(1 -M 2 x ) , (4.9)
where the terms

A 1 = 2ik 0 M x , A 0 = v x ∂ x (v x c -2 0 ) -ρ -1 0 ∂ x (ρ 0 ), B 1 = ωv x ∂ x (c -2 0 ), (4.10)
have been written relatively to their homogeneity degree in ω. Eliminating Λ -from (4.9) gives an operator equation for the outgoing characteristic

1 -M 2 x (Λ + ) 2 -i (A 1 + A 0 ) Λ + + i(1 -M 2 x )Op ∂ x λ + = ω 2 + ∂ 2 y -iB 1 mod OPS -∞ . (4.11)
Expanding the symbol of the DtN operator as classical symbols of decreasing orders, we can identify the second-order terms in (4.11) to find the principal symbol

λ ± 1 = -M x k 0 ± k 2 0 -(1 -M 2 x )ξ 2 1 -M 2 x , (4.12)
which has been chosen to match the DtN symbol in the uniform flow situation. We recall that the analysis is valid microlocally and does not hold a priori in the grazing region, that is when the term under the squareroot cancels. The next symbol is computed thanks to the composition rule (2.17). The identification of the first-order terms in (4.11) gives

2λ + 0 (1 -M 2 x )λ + 1 + M x k 0 = -iB 1 + iA 0 λ + 1 -i(1 -M 2 x )∂ x λ + 1 , (4.13)
or alternatively, by substitution of the principal symbol in the left-hand side

λ + 0 = i (M 2 x -1)∂ x λ + 1 + A 0 λ + 1 -B 1 2 k 2 0 -(1 -M 2 x )ξ 2 . (4.14)
After the formal computation of the principal symbol derivative, we use the mass conservation relation (4.5) and obtain the simplified expression

λ + 0 = -i ∂ x (ρ 0 ) 2ρ 0 k 2 0 -ξ 2 k 2 0 -(1 -M 2 x )ξ 2 + i ∂ x (c 0 ) 2c 0 k 2 0 + M 2 x ξ 2 k 2 0 -(1 -M 2 x )ξ 2 , (4.15)
which is of order zero with respect to (ω, ξ). Intuitively, these terms correspond to an amplitude correction factor of the wave. For M x = 0, we retrieve the zeroth order symbol from the heterogeneous Helmholtz situation (2.23).

Construction of the ABCs

To build ABCs, the next step is to associate appropriate operators to the symbols, and give them a local representation. Based on the principal symbol, we consider

Λ 1 = -M x k 0 + k 2 0 + (1 -M 2 x )∆ Γ 1 -M 2 x , (4.16)
where the inverse Fourier representation is exact in the microlocal sense since neither M x nor the boundary Γ depend on the transverse direction. It coincides with the exact DtN map in the uniform mean flow case.

For its implementation, we write the operator in the form

Λ 1 = k 0 -M x + √ 1 + X 1 -M 2 x , X = (1 -M 2 x ) ∆ Γ k 2 0 , (4.17)
such that the complex Padé approximants (3.14)-(3.15) can be used. For the next symbol λ + 0 , we consider the operator

Op(λ + 0 ) = i 2 ∂ x (c 0 ) c 0 k 2 0 -M 2 x ∆ Γ - ∂ x (ρ 0 ) ρ 0 k 2 0 + ∆ Γ k 2 0 + (1 -M 2 x )∆ Γ -1 , (4.18)
leading to Λ 2 = Λ 1 + Op(λ + 0 ). We can finally define the approximate DtN map by keeping the operators based on the first M symbols, i.e.

∂ n u = -iΛ M u on Γ, (4.19) setting Λ M = M -2 j=-1 Op(λ + -j ).
The resulting ABCs are again denoted by ABC N,α M after using the complex Padé approximants with parameters (N, α). For grazing waves, one might introduce a complexified frequency ω ε in the inverse operator arising in the definition of Λ 2 .

Numerical study

We consider the 2D duct case from Figure 1 applied to the convected wave operator (4.4). We choose a sigmoid x-velocity profile that is inspired from axial flow variations in a turbofan engine intake. The mean flow properties are shown in Figure 5. We choose the velocity profile as

v x (x) = v 0 + δ v tanh(s(x -x c )), (4.20)
where we fix the mean value to be v 0 = -0.4c ∞ , the slope strength s = {5, 10}, the deviation δ v = 0.3c ∞ and the center point x c = 1. The value c ∞ = 340 serves as a reference for large x. We also set a reference velocity v ∞ = -0.1c ∞ . The integration of (4.5) with ρ ∞ = 1.2 leads to the mean density

ρ 0 (x) = ρ ∞ v ∞ v x (x)
, and Bernoulli's relation relative to the reference point gives the speed of sound profile

c 0 (x) = c ∞ 1 + γ -1 2 v 2 ∞ -v 2 x (x) c 2 ∞ .
We recall that these profiles are used as an example and do not aim to represent any realistic physical situation. For this case, we expect the exact solution u ex to have the modal structure

u ex (x, y) = cos (k y y) A(x), k y = nπ H , n ∈ N, (4.21)
where A is an unknown oscillatory function. velocity profile has been chosen such that the initial data g can be based on the uniform flow situation. For example, a single propagative mode is enforced as

ρ 0 (1 -M 2 x )g = i -M x k 0 + k 2 0 -(1 -M 2 x )k 2 y cos(k y y). (4.22)
The variational formulation associated to the duct problem is built directly from (4.1). After multiplication by ρ 0 > 0, we look for u ∈ H 1 (Ω) such as:

∀v ∈ H 1 (Ω), Ω ρ 0 ∇u • ∇v - ρ 0 c 2 0 D 0 u D 0 vdΩ + i Γ3 ρ 0 k 0 M x u vdΓ 3 - Γ1 ρ 0 k 0 M x u vdΓ 1 (4.23) - Γ3 ρ 0 (1 -M 2 x )∂ n u vdΓ 3 = Γ1 ρ 0 (1 -M 2 x )g vdΓ 1 , (4.24)
where

D 0 = iω + v x ∂ x .
The practical finite element implementation of the ABC N,α 2 is similar to the one described in Section 3.3, in particular regarding the introduction of auxiliary functions into the variational formulation following (3.24)- (3.26). In all the simulations, we use linear Q4 elements of size h = 1/40 and the p-FEM order is fixed to p = 8.

Validation for a uniform mean flow

For a uniform mean flow, we have Λ = Λ 1 . The exact solution reads

u ex (x, y) = cos (k y y) e -ikxx , k y = nπ H , n ∈ N, (4.25)
and the exact DtN map Λ = k x is given by the wavenumber

   k x = 1 1-M 2 x -M x k 0 + k 2 0 -(1 -M 2 x )k 2 y , k 2 0 > (1 -M 2 x )k 2 y , k x = 1 1-M 2 x -M x k 0 -i (1 -M 2 x )k 2 y -k 2 0 , k 2 0 < (1 -M 2 x )k 2 y . (4.26)
We measure the L 2 -error as defined in (3.29) for M x = -0.8 and look at the performance of ABC N,α in the hyperbolic and elliptic zones, shown respectively in Figures 6a and6b. A duct of shorter length is used to evaluate the error in the elliptic zone since the exact solution is evanescent and decays exponentially in the duct. Because the mean flow is uniform, the ABC accuracy tends to the one obtained with the exact DtN operator as the number of auxiliary functions N grows. The angle α corresponds to a trade-off between the attenuation of evanescent and propagative modes. When k 0 → 0, the mode does not propagate and the square-root approximation is no longer accurate. As expected from microlocal analysis, none of the conditions are able to tackle grazing waves. They are located respectively at k 0 ≈ 23 and k 0 ≈ 15 in Figures 6a and6b. Similar results are obtained when M x > 0.

Validation for a non uniform mean flow

For a non-uniform flow, there is unfortunately no analytical solution for our problem, and we must resort to a numerical alternative. Perfectly matched layers seem a reliable option, although their efficiency for heterogeneous problems remains unclear. To bypass the effect of the medium heterogeneity in the PML, we follow an approach similar to the one proposed in [START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]. The idea is to extend the computational domain and further apply a PML in the region where the mean flow is uniform. The PML is terminated by a homogeneous Neumann boundary condition. A second issue for convected problems is the presence of inverse modes that makes the PML unstable and ineffective in practice. Fortunately stabilization techniques are available and we will use the stabilized version of the PML described in [START_REF] Marchner | Stable perfectly matched layers with Lorentz transformation for the convected Helmholtz equation[END_REF]. The strategy to compute the reference solution is illustrated in Figure 7. Three closed domains are defined: the physical domain

Ω = [0, L] × [0, H], the extended domain Ω ext = [L, L ref ] × [0, H] and the PML domain Ω pml = [L ref , L pml ] × [0, H].
The relative L 2 -error (in %) is then measured in the truncated domain Ω as ∂ In order to evaluate the efficiency of the ABCs we consider a set of frequencies and ABC positions L for a given mode n. The selected frequencies span the elliptic, grazing and hyperbolic regimes. Figure 8 shows the real part of some reference solutions in Ω for L ref = 4 and L = 2. According to the mean flow profile in Figure 5 we set k ∞ = ω/c ∞ to be the reference wavenumber. The modal behaviour has interesting features in the grazing regime: at k ∞ = 20 and k ∞ = 25, the mode is cut-on at the input and cut-off when it reaches x = 2. When this behaviour occurs within the computational domain, there is an interference pattern that does not appear in the uniform flow case. The mode has to reflect at a turning point (see e.g. [START_REF] Rienstra | An Introduction to Acoustics[END_REF] section 8.5), where the square-root term of the principal symbol (4.12) vanishes. The transition range is approximately k ∞ ∈ [START_REF] Turkel | Boundary conditions and iterative schemes for the Helmholtz equation in unbounded regions[END_REF]25] and we do not expect any ABCs to perform well in this regime. In the propagative regime (such as k ∞ = 40 and k ∞ = 60), the wavelength increases in the duct due to the mean flow variation.

E L 2 = 100 u pml|Ω -u h L 2 (Ω) u pml|Ω L 2 (Ω) , (4.27
k 0 E L 2 (%) ABC 0,0 ABC 1,0 ABC 4,0 ABC 8,0 (a)
k 0 ABC 4,0 ABC 4,-π/4 ABC 4,-π/2 ABC 4,-3π/4 ABC 4,-π (b) H = 1, L = 0.2, n = 8
n u = 0 0 Ω Ω ext Ω pml H L pml L ref L x y
For a general velocity profile, we have no a priori information on the location of the turning point. However thanks to the separable structure of the solution one may infer that it is reduced to a fixed value x t . Since the mean flow is analytically prescribed, we can find the turning point x t whenever the equality

ω 2 = c 2 0 (x t ) -v 2 x (x t ) nπ H 2 (4.28)
holds. If x t lies in the duct, its value is reported in Figure 8. We present in Figures 9 and 10 the relative L 2 -error as a function of the ABC position L and input wavenumber k ∞ for two velocity profiles of respective slopes s = 5 and s = 10. When the propagative regime is well developed (k ∞ > 25), the condition ABC N,α 2 shows in general a gain in accuracy of more than one order of magnitude compared to ABC N,α 1 , especially for the cases where the mean flow derivative is non-zero. In this model, even a small value of the velocity derivative has an impact on the ABC accuracy. It partially is less pronounced for s = 10 than for s = 5. Note that the contribution due to the density variations on the L 2 -error is more important than the one due to the variation of the speed of sound. This is directly linked to the value of the factors ∂ x (ρ 0 )/ρ 0 and ∂ x (c 0 )/c 0 in equation (4.18).

As expected, we observe challenging situations for the ABCs in the regime k ∞ ∈ [START_REF] Turkel | Boundary conditions and iterative schemes for the Helmholtz equation in unbounded regions[END_REF]25], where the wave is grazing and close to the critical situation from (4.28). More precisely when the ABC is located before the turning point L < x t it has to capture a wave that propagates along both the positive and negative xdirection, and such a behaviour is not tackled by our method. Although improvements can be found locally by tuning a dissipation parameter ε for the inverse operator of the zeroth order symbol (as it was done in Section 3.3), the results are not general enough to be included. Finally for low frequencies where the mode is evanescent, the error can be reduced with a higher value of α.

To sum up, the most efficient condition is ABC N,α 2 and the rules for tuning (N, α) are the same as in Section 3.3. It confirms the effectiveness of the operators obtained from microlocal analysis. As said before, we recommend to always rotate the branch-cut of the principal symbol, since in practice we want the attenuation of all modes and do not know the turning point location, in particular where there is a superposition of modes. We claim that the presented results stay valid in a multi-modal situation.

Transverse variation

This section describes a third situation which is more complex and highlights some limitations of the method. We consider the heterogeneous Helmholtz equation subject to a transverse variation of the speed of sound c 0 (x, y) = c 0 (y) and density ρ 0 (x, y) = ρ 0 (y) For the two-dimensional straight duct problem, we can still use the separation of variables and write the exact solution as a linear combination of modes

∂ 2 x u + ρ -1 0 ∂ y (ρ 0 ∂ y )u + k 2 0 u = 0, k 0 = ω/c 0 . ( 5
u ex (x, y) = n≥0 U n (y)e -ik n x x , (5.2)
and the exact DtN operator can be explicitly written for a single mode as

Λ ± = ±k n x on Γ 3 , (5.3)
where the k n x are the propagation constants given by the Sturm-Liouville eigenvalue problem

U (y) + ρ -1 0 ∂ y (ρ 0 )U (y) + k 2 0 (y)U (y) = (k x ) 2 U (y), U (0) = 0, U (H) = 0, (5.4)
where we consider homogeneous Neumann boundary conditions on the duct upper and lower walls. In the constant coefficient case, the eigenmodes and eigenvalues are respectively given by

U n (y) = cos (k y y) , (k n x ) 2 = k 2 0 -k 2 y , k y = nπ H , n ∈ N. (5.5)
Compared to the previous sections, the main difficulty is that all the information is contained in the eigenvalues, which moreover depend on the boundary conditions on the walls. Such problems are common in various fields of physics such as but not limited to acoustics [START_REF] Rienstra | An Introduction to Acoustics[END_REF], optics [START_REF] Okamoto | Fundamentals of Optical Waveguides[END_REF] or geophysics [START_REF] Demanet | Discrete symbol calculus[END_REF].

Computation of the dispersion relation

For some specific profiles, one may derive the general solution of the boundary value problem (5.4) and find the associated eigenvalues by looking for the zeros of a transcendental equation. We solve the problem semi-analytically thanks to a one-dimensional spectral method, which leads to highly accurate results for an arbitrary profile [START_REF] Trefethen | Software, Environments, and Tools[END_REF]. The approach allows to quickly compute the dispersion relation, which relates the propagation constants k n

x to the input frequency ω. For the speed of sound profile we choose the Gaussian profile given in [START_REF] Vion | Double sweep preconditioner for optimized Schwarz methods applied to the Helmholtz problem[END_REF] c 0 (y) = 1.25 1 -0.4e -32(y-H/2) 2 , H = 1, (

and select the density such as ρ 0 (y) = c 2 0 (y). Hence, one may identify these quantities as a refractive index for transverse electric and magnetic modes in optical waveguides [START_REF] Okamoto | Fundamentals of Optical Waveguides[END_REF]. Figure 11 compares the dispersion relations associated to propagative modes for the homogeneous case and for the Gaussian profile. A qualitative observation indicates that the speed of sound variations have the most significant impact on the modal structure. The eigenvalues associated to the Gaussian profile seem to be driven by one of the two lines of equation 1 = ω/ max(c 0 ) and 2 = ω/ min(c 0 ). The limit case would be a constant piecewise profile, for which the eigenvalues accumulate along the lines 1 and 2 in the high frequency limit [START_REF] Okamoto | Fundamentals of Optical Waveguides[END_REF]. 1. the guided modes, which are decaying close to the duct walls; 2. the radiating modes, which are fully oscillating.

Although not shown, evanescent modes are also present and have a purely imaginary propagation constant. They are highly oscillatory and decay exponentially along the x-direction.

Transverse variation of the density 5.2.1. Strategies for the choice of the principal symbol

We first consider the case of a variable density and constant speed of sound c 0 = 1 and present two strategies for selecting the principal symbol. The first one consists in choosing the classical symbol λ + 1 = ω 2 -ξ 2 , together with the zeroth order symbol given in (2.23)

λ + 0 = -iξ∂ y (ρ 0 ) 2ρ 0 ω 2 -ξ 2 .
As a consequence, the use of Λ S as an ABC should accurately represent the DtN operator. The modified square-root operator is then written as

Λ S = ω √ 1 + Z, Z = ∆ Γ + ρ -1 0 ∇ Γ ∂ y (ρ 0 )∇ Γ ω 2 , (5.13)
and localized using the complex Padé approximants (3.15). The semi-classical ABC is denoted ABC N,α S .

Weak formulation and numerical tests

We follow the methodology from the previous sections and consider the duct case problem boundary value problem depicted in Figure 1. The variational form associated to the domain equation (5.1) reads: find

u ∈ H 1 (Ω) such that ∀v ∈ H 1 (Ω), Ω {∇u • ∇v -ρ -1 0 ∂ y (ρ 0 )∂ y u v -k 2 0 u v} dΩ + i Γ3 Λu v dΓ 3 = Γ1 g v dΓ 1 . (5.14)
As input boundary condition we enforce a given mode n from the homogeneous case and set its amplitude to unity g = cos nπ H y , n ∈ N.

(5.15)

Since g is not an eigenfunction of the Sturm-Liouville problem (5.4), the solution exhibits the propagation of multiple-modes. Moreover, the transverse oscillations present a more complex pattern compared to the homogeneous situation, as shown in Figure 13a. As a reference we consider a numerical solution using a large PML. The physical domain Ω is extended by Ω pml = [L, L pml ] × [0, H] and the relative L 2 -error (in %) is measured as

E L 2 = 100 u pml|Ω -u h L 2 (Ω) u pml|Ω L 2 (Ω) . (5.16)
In order to validate the PML in the transverse heterogeneous case we have performed a mode-by-mode analysis where the semi-analytical solution for a fixed mode u ex (x, y) = U (y)e -ikxx is generated thanks to the one-dimensional spectral method used for the problem (5.4). Each eigenfunction U (y) is used as input boundary condition and the DtN map is given by the corresponding eigenvalue k x . The variational form with the PML follows from the substitution ∂ x = ∂ x -iσ(x)/ω in (5.14), where σ is zero in Ω and given by the hyperbolic function σ(x) = σ 0 /(L pml -x), σ 0 ∈ C in Ω pml . Since we later assess the ABCs accuracy in a multi-modal situation special attention must be paid to grazing modes which are characterized by a high phase velocity. After numerical experiments we conclude that the value σ 0 = 40 together with a large PML |L pml -L| = 40h yields an optimal accuracy for all modes. As a result we obtain a reference solution that has the precision of the DtN map, and by linearity we suppose it holds in the multi-modal case. Additional considerations can be found in the appendix of [START_REF] Marchner | Stable perfectly matched layers with Lorentz transformation for the convected Helmholtz equation[END_REF]. For the simulations we use again Q4 elements of size h = 1/40 and set the p-FEM order to p = 6.

For the practical implementation of the ABC N,α S and ABC N,α 2 with finite elements, we use an augmented system with auxiliary functions on the boundary as explained in Section 3. For ABC N,α S , the boundary integral on Γ 3 with Padé approximants gives, ∀v ∈ H 1 (Ω),

i Γ3 ABC N,α S u v dΓ 3 = i Γ3 ωC 0 u v dΓ 3 -i Γ3 A ω ∇ Γ3 ϕ • ∇ Γ3 v dΓ 3 + i Γ3 A ∂ y (ρ 0 ) ωρ 0 ∇ Γ3 ϕ • v dΓ 3 , (5.17)
where ϕ are given through the N auxiliary equations

∀v ∈ H 1 (Γ 3 ), Γ3 ω 2 ϕ v dΓ 3 - Γ3 B ∇ Γ3 ϕ • ∇ Γ3 v dΓ 3 + Γ3 B ∂ y (ρ 0 ) ρ 0 ∇ Γ3 ϕ • v dΓ 3 = Γ3 ω 2 u v dΓ 3 .
(5.18)

For the implementation of ABC N,α

2

, we need to introduce a total of 2N supplementary auxiliary equations. For simplicity, the same tuning parameters (N, α) are used as done in the implementation of ABC N,α 1 for the localization of the inverse square root operator. The second set of auxiliary functions is denoted (ψ

1 , • • • , ψ N ). The variational formulation reads, ∀v ∈ H 1 (Ω), i Γ3 ABC N,α 2 u v dΓ 3 = i Γ3 ωC 0 u v dΓ 3 -i Γ3 A ω ∇ Γ3 ϕ • ∇ Γ3 v dΓ 3 + i Γ3 ∂ y (ρ 0 ) 2ωρ 0 ∇ Γ3 ψ • v dΓ 3 , (5.19)
where ϕ and ψ are given through the 2N auxiliary equations We report in Figure 13b , where the error reaches a plateau when N grows.

∀v ∈ H 1 (Γ 3 ), Γ3 ω 2 ϕ v dΓ 3 - Γ3 B ∇ Γ3 ϕ • ∇ Γ3 v dΓ 3 = Γ3 ω 2 u v dΓ 3 , (5.20) ∀µ ∈ H 1 (Γ 3 ), Γ3 S ω 2 ψ µ dΓ 3 - Γ3 ∇ Γ3 ψ • ∇ Γ3 µ dΓ 3 = Γ3 R ω 2 u µ dΓ 3 . ( 5
It confirms the approximate representation of the DtN map. Although not reported, one may use a second order Taylor approximation for high frequencies. In that case, the operators Λ 2 and Λ S with the two first symbols reduce to the same expression

Λ 2 2 = ω + ∂ y (ρ 0 ) 2ωρ 0 ∇ Γ + ∆ Γ 2ω . (5.22)
Further numerical tests show that its performance is close to the one of ABC N,α 2 . In Figure 13c the length of the duct is varied at a fixed frequency ω = 30 to examine the effect of the rotation branch-cut in the elliptic regime. By contrast to the case of a longitudinal heterogeneity, evanescent modes are only present when the ABC is close enough to the input boundary. We observe as expected that the reflection of evanescent modes becomes negligible beyond a certain value of L. For the condition ABC N,α S , the rotation angle has a similar effect as in the homogeneous case. Evanescent modes are better damped when α increases but to the detriment of propagative modes. Therefore, we recommend to increase the number of auxiliary fields together with α as proposed in the previous sections.

Grazing waves deteriorate the overall ABC quality in such a multi-modal situation, especially in the low frequency regime. All the ABCs show error peaks for a discrete set of frequencies. They correspond to evanescent modes becoming cut-on and can be accurately predicted by the dispersion relation. As a workaround we may add a small dissipation parameter to the frequency within the square-root operator [START_REF] Antoine | An improved surface radiation condition for high-frequency acoustic scattering problems[END_REF]. It has the effect to remove the singular behaviour of the operator but deteriorates the attenuation of all other modes. However this strategy did not significantly improve the ABCs.

Transverse variation of the speed of sound

The effect of the heterogeneity is now examined with a transverse variation of the speed of sound and a density set to unity. Following the analysis from Section 2 we choose the following operator to represent the principal symbol

Op(λ + 1 ) = ω 2 c -2 0 + ∆ Γ mod OPS -2 , λ + 1 = ω 2 c -2 0 -ξ 2 . (5.23)
We remark that the next symbols of Op(λ + 1 ) are exactly the DtN symbolic expansion such that we have

Op(λ + ) = ω 2 c -2 0 + ∆ Γ mod OPS -∞ , (5.24)
which means that the square-root operator is an accurate representation of the DtN map from a microlocal point of view. This result can be obtained more directly from Nirenberg's factorization procedure (2.8) since the speed of sound does not depend on the direction of propagation. The situation is similar as for the operator Λ S (5.11) in the case of a variable density. However it is here more involved to find a local representation of the square-root operator, because the speed of sound variations affect the real part of λ + 1 for which the sign may change along the non-reflecting boundary.

We examine two different ways to approximate the square-root operator, denoted Λ ω and Λ k0 to allude to the factor in front of the square root

Λ k0 = k 0 1 + ∆ Γ k 2 0 , Λ ω = ω 1 + (c -2 0 -1) + ∆ Γ ω 2 , (5.25)
where we recall that k 0 = ω/c 0 . The second approximation is often used in one-way modeling for beam propagation in optical waveguides [START_REF] Lu | Some techniques for computing wave propagation in optical waveguides[END_REF]. Once again, complex Padé approximants are used for the localization procedure. The obtained ABCs are denoted ABC N,α which is nothing but the difference between the DtN and ABC principal symbols. We plot both functions in a logarithmic scale in Figure 14 for the parameters N = 8 and α = -π/2. The condition ABC N,α ω seems to be a better candidate than ABC N,α k to approximate the DtN principal symbol.

Both functions show a singularity along the characteristic line of equation ξ = ±ωc -1 0 (y) where the square root vanishes. On the contrary to the homogeneous or x-profile case, the turning region is not reduced to a single point and both oscillatory and evanescent behaviour may be present for a fixed value of ξ (a given mode). This behaviour is intimately linked to the spectral properties of the Sturm-Liouville problem. log 10 f k0 (ξ, y) log 10 f ω (ξ, y) We plot in Figure 15 an example of numerical solution and report the relative L 2 -error defined in (5.16) for different conditions. As for the case of transverse variable density, we use a large PML of width |L pml -L| = 100h as a reference solution.

The condition ABC N,α ω is more accurate than ABC N,α k0 , which can be understood from Figure 14 by a better approximation of the DtN principal symbol. The parameters (N, α) have a marginal impact on ABC N,α k0 , while we observe an improvement with N for ABC N,α ω . This suggests that the condition ABC N,α ω is a good approximation of the DtN map, although a more efficient localization of the square-root operator may be sought. As expected, grazing waves are not well tackled by the proposed ABCs. A more advanced analysis is required and is out of scope of this work. For example, a microlocal cut-off function might be used to handle the singularity [START_REF] Demanet | Discrete symbol calculus[END_REF][START_REF] Stolk | A pseudodifferential equation with damping for one-way wave propagation in inhomogeneous acoustic media[END_REF]. To conclude, Padé approximants provide a comprehensive way to design non-reflecting boundaries in the transverse heterogeneous case, but requires a careful implementation. It reaches its limitations for complex heterogeneous media, for example in the presence of corners and/or non-smooth heterogeneities. For such situations, algebraic or numerical techniques to find sparse representations of the DtN map could be more appropriate [START_REF] Lee | Root-finding absorbing boundary condition for poroelastic wave propagation in infinite media[END_REF][START_REF] Bélanger-Rioux | Compressed absorbing boundary conditions via matrix probing[END_REF][START_REF] Druskin | Compressing variable-coefficient exterior Helmholtz problems via RKFIT[END_REF].

Conclusion

In this paper, we built and numerically analyzed two families of ABCs for the heterogeneous Helmholtz equation, based on a general strategy which i) needs the computation of the first symbols of the DtN map by microlocal analysis techniques, and ii) next approximates the symbols by a Taylor or complex Padé expansion, resulting in local ABCs. These ABCs can be easily implemented in a high-order finite element solver. We analyzed in details three configurations • a longitudinal heterogeneous duct problem with linear x-dependent speed of sound c 0 (x, y) = c 0 (x);

• a convected wave duct problem with non-uniform mean flow parameters along the x-direction;

• and finally a duct problem with a transverse density ρ 0 (x, y) = ρ 0 (y) and speed of sound c 0 (x, y) = c 0 (y).

For the two first cases, our strategy led to the construction of relatively accurate local ABCs, called ABC N,α 2 , based on the two first symbols of the DtN map and Padé approximants with tuning parameters (N, α). The last situation depends on an appropriate choice for the square-root operator, which has been guided by microlocal analysis and studies on optical waveguides. It led to fairly efficient ABCs, called ABC N,α S and ABC N,α ω . Aside from specific cases, grazing waves are not correctly handled and require alternative strategies that will be studied in further works. Additional simulations, most particularly in 3D, would give more information about the capacity of the ABCs to simulate industrial situations.

A promising application of the proposed ABCs is to use them as transmitting boundary conditions in Schwarz domain decomposition methods for solving the heterogeneous Helmholtz equation. Indeed, the choice of well-designed transmitting boundary conditions provides fast converging iterative solvers for the simulation of large scale physical situations [START_REF] Boubendir | A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation[END_REF]. This constitutes a natural next step of our developments, where we will consider ABC N,α

Figure 1 :

 1 Figure 1: Sketch of the numerical case: two-dimensional acoustic propagation in a heterogeneous duct.

  the real part of the reference numerical solution for single frequencies and the location of the turning point as a function of ω. ω = 10, x t = 0.69 ω = 60, x t < 0 ω = 20, x t = 0.16 x t (ω)

Figure 2 :

 2 Figure 2: Real part of the reference numerical solution u h for fixed frequencies ω and location of the turning point xt.

Figure 3 :

 3 Figure 3: Relative L 2 -error for Taylor-based (left) and Padé-based (right) ABCs as a function of the input frequency ω for α = 0 and ε = 0. Reference solution().

2 . 2 ,

 22 Variation of the angle α for ABC 4,α 2 , ε = 0. Variation of the parameter ε for ABC 6,-π/2 Variation of the number of auxiliary fields for ABC N,-π/2 ε = 0.

Figure 4 :

 4 Figure 4: Influence of the tuning parameters for the Padé-based conditions ABC N,α 2 on the relative L 2 -error when varying the ABC location L. Reference solution().

Figure 5 :

 5 Figure 5: Mean flow profile along the x-direction for the convected wave problem.

H = 0. 5 , L = 1 , n = 6 2

 516

Figure 6 :

 6 Figure 6: Relative L 2 -error for the condition ABC N,α . Left: propagative regime for different number of auxiliary fields N and α = 0. Right: evanescent regime for different rotation angles α and N = 4. Reference solution ().

Figure 7 :

 7 Figure 7: Illustration of the strategy to compute the reference solution. The mean flow derivative is assumed to be zero at x = L ref . The ABC boundary is at x = L.

  (a) k∞ = 20, n = 4, xt ≈ 0.89 (b) k∞ = 25, n = 4, xt ≈ 1.56 (c) k∞ = 40, n = 4 (d) k∞ = 60, n = 4

Figure 8 :

 8 Figure 8: Real part of the reference solution for fixed frequencies

2 Figure 9 :

 29 Figure 9: Relative L 2 -errors for ABC 4,-π/4 1 ( ) and ABC 5,-π/4 2 ( ) in comparison with a PML reference solution. Velocity profile with slope s = 5.

Figure 10 :

 10 Figure 10: Relative L 2 -errors (%) for ABC 4,-π/4 1 ( ) and ABC 5,-π/4 2 ( ) in comparison with a PML reference solution. Velocity profile with slope s = 10.

  c 0 (y) = 1.25 (a) Homogeneous c 0 (y) = 1.25(1 -0.4e -32(y-H/2) 2 ) (b) Gaussian

Figure 11 :

 11 Figure 11: Dispersion relation associated to propagative modes for a homogeneous and Gaussian profile.

Figure 12

 12 Figure12presents the propagative eigenmodes at the fixed frequency ω = 37, where two types of modes can be identified:

. 21 )

 21 The results are first analyzed when for the fixed ABC position L = 2, where we expect propagative modes (corresponding to the hyperbolic zone) to have the largest impact on the ABC quality. The rotation angle is set to α = -π/4.

  Influence of the frequency ω for L = 2. of the ABC location at ω = 30.

Figure 13 :

 13 Figure 13: Real part of the reference solution for the mode n = 4 and relative L 2 -errors when varying the frequency and ABC location.

1 , 2 performs better than ABC N,α 1 especially

 121 the relative L 2 -error as a function of the input frequency ω for ABC N,α ABC N,α 2 and ABC N,α S . The condition ABC N,α in the high-frequency regime, and ABC N,α S outperforms the two other conditions. The microlocal theory gives a consistent interpretation to the results. While ABC N,α 2 only incorporates the contribution of the zeroth order symbol, the condition ABC N,α S encodes the full asymptotic symbolic expansion of the DtN operator thus leading to an excellent ABC. Besides, the performance of ABC N,α S can be improved by increasing the number of auxiliary functions N . It means that Padé approximants are able to accurately represent the information contained in the DtN symbols. This is not the case for ABC N,α 1 neither ABC N,α 2

  k0and ABCN,α ω . In order to analyze their potential accuracy, we compute the function of two variables at a fixed frequency ωf k0 (ξ, y) = λ + 1 -σ ABC N,α k0 , f ω (ξ, y) = λ + 1 -σ ABC N,α ω ,(5.26)

Figure 14 :

 14 Figure 14: Difference in logarithmic scale between the DtN principal symbol and its approximation by Padé approximants (N = 8, α = -π/2) for the Gaussian speed of sound profile at ω = 30.

Figure 15 :

 15 Figure 15: Real part of the reference solution at ω = 55 for the mode n = 4 and relative L 2 -error for L = 2 when varying the frequency for the Gaussian speed of sound profile.

2 ,

 2 ABC N,α S and ABC N,α ω as transmitting boundary conditions in optimized Schwarz DDMs.

https://gitlab.onelab.info/gmsh/fem

https://gitlab.onelab.info/gmsh/fem/-/tree/master/examples/helmholtz2d/waveguide, commit number ad4e0321, April 9th 2021

https://gitlab.onelab.info/gmsh/fem/-/tree/master/examples/helmholtzFlow/waveguide, commit number ad4e0321, April 9th 2021

Acknowledgements

This work was performed as part of the CIFRE contract No. 2018/1845 funded by Siemens Industry Software SAS and Association Nationale de la Recherche et de la Technologie (ANRT). This research was also funded in part through the ARC grant for Concerted Research Actions (ARC WAVES 15/19-03), financed by the Wallonia-Brussels Federation of Belgium.

Data Reproducibility

The source code to reproduce the simulation data is available online at https://gitlab.onelab.info/ gmsh/fem in the folder examples/helmholtz2d/waveguide 2 for Sections 3 and 5 and in the folder examples/helmholtzFlow/waveguide 3 for Section 4. Installation instructions are given in the corresponding README.md file.

By contrast to the case of longitudinal heterogeneity, the zeroth order symbol depends on the transverse variable and the choice of an appropriate operator is less natural. Nevertheless, we can choose here the operators such that the following correspondence holds

Note that the choice for Op(λ + 0 ) is not unique and that a different choice might be more relevant to the situation. The computation of the operator asymptotic expansion at the symbol level allows to compare it with λ + 0 and next obtain the operator regularity estimate. We build the approximate DtN maps Λ M based on the first M symbols in the high frequency limit ω → ∞ as

and use a rotated branch-cut Padé approximation for the inverse square-root [START_REF] Chaillat | Fast iterative boundary element methods for high-frequency scattering problems in 3d elastodynamics[END_REF] (

(5.9)

The resulting ABCs are again denoted ABC N,α M . The second strategy is a semi-classical approach, which enriches the information contained in the principal symbol. It consists in keeping all the terms in the right-hand side of (2.16). The modified principal symbol is

.10)

A natural choice for its operator representation would be to choose the modified square-root

(5.11)

One observes that the principal symbol of Λ S is σ (Λ S ) = ω 2 -ξ 2 , and that the next symbols coincide with the sequence {λ + -j , j ≥ -1} from the first strategy. Hence we have Op(λ + ) = Λ S mod OPS -∞ .

(5.12)

Appendix A. Optimized reflection coefficient for the linear profile

For a fixed transverse mode, the Airy differential equation

admits the outgoing and associated ingoing linearly independent solutions [START_REF] Lévy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF]

We set the outgoing wave amplitude to unity, and express the total field as the sum of the forward and reflected fields

where R is the reflection coefficient. Since the DtN operator links the solution with its normal derivative on the boundary, we obtain for the total field The next step is to analyze the accuracy of the different DtN approximations obtained by microlocal analysis. This could be done with respect to the transverse mode, the frequency or the position of the ABC. It is illustrated in Figure A.16 as a function of the frequency. As predicted theoretically, the precision of the ABC increases with its approximation order but deteriorates at the turning point.

The inverse operators arising from the asymptotic analysis introduce a singularity at the turning point x t . We propose to focus on the second order condition and add a dissipation parameter ε such that

We next look for the minimum of the reflection coefficient at the turning point, yielding a simple minimization problem

where a ε (x t , ω t ) simplifies to

The problem can be solved explicitly:

The present analysis is however only valid when an analytic solution is available. We use the value ε opt for our computations, and leave alternative strategies (e.g. purely numerical) to optimize the reflection coefficient for a future work. We plot in Figure A.17 the reflection coefficient before and after the approximation of the square-root term by Padé approximants. The condition R < 1 seems to hold if N is taken large enough. The choice of α is important to ensure R < 1 in the whole frequency range. For example, the choice α = -π/3 needs at least N = 5 to ensure R < 1 in the entire elliptic zone. Appendix B. Formal computation of the DtN approximation for the linear profile 1 from sympy import diff , series , Function , symbols , poly 2 from sympy import I , discriminant , simplify , sqrt , collect 3 from sympy import init_printing 4 init_printing () 5