Discrete transparent boundary conditions for the two-dimensional leap-frog scheme

We develop a general strategy in order to implement (approximate) discrete transparent boundary conditions for finite difference approximations of the two-dimensional transport equation. The computational domain is a rectangle equipped with a Cartesian grid. For the two-dimensional leapfrog scheme, we explain why our strategy provides with explicit numerical boundary conditions on the four sides of the rectangle and why it does not require prescribing any condition at the four corners of the computational domain. The stability of the numerical boundary condition on each side of the rectangle is analyzed by means of the so-called normal mode analysis. Numerical investigations for the full problem on the rectangle show that strong instabilities may occur when coupling stable strategies on each side of the rectangle. Other coupling strategies yield promising results.


Introduction
In this article, we are concerned with the construction and numerical implementation of discrete transparent boundary conditions for linear transport equations. Our goal is to derive numerical boundary conditions that minimize the parasitic wave reflections at the boundary of the computational domain. This research area has been increasingly active since the pioneering work by Engquist and Majda [19] and our goal here is to explain why following the same "small frequency approximation" strategy as in [19] can be successful at the fully discrete level when one deals with a rectangular computational domain. Namely, we shall construct fully discrete, approximate transparent boundary conditions for the two-dimensional leap-frog scheme on a rectangle. Our numerical boundary conditions will not be exactly transparent since the derivation of such conditions on a rectangle seems to be out of reach. Our strategy is to start from the exact transparent boundary conditions for a half-space and, as in [19], to localize them with respect to the tangential variable by means of a suitable "small frequency approximation". This general strategy may be applied to any finite difference scheme. In terms of wave reflection magnitude, our numerical simulations show that our approximate transparent boundary conditions Keywords and phrases. Transport equation, leap-frog schemes, transparent boundary conditions, stability.
are more efficient than more standard, purely local extrapolation strategies, hence, in our opinion, the relevance of implementing discrete time convolution operators as described below.
We focus here on the discretization of the two-dimensional transport equation by means of the leap-frog scheme. The nice feature of this approximation is that it does not dissipate high frequency signals -and therefore allows for a precise analysis of wave reflections on each side of the rectangle -and, moreover, its stencil exhibits some "dimensional splitting". The latter feature will be helpful in our construction since we shall not have to develop a specific numerical treatment for the four corners of the computational domain. Hence we shall focus on the numerical boundary conditions on each of the four sides of the rectangle and on their coupling through the eight mesh points that are closest to the four corners. More simple, purely local, outflow strategies [1,21,24,35] may be implemented for dissipative schemes (e.g., the Lax-Wendroff scheme). Even though these extrapolation boundary conditions are much cheaper from a computational point of view and are quite efficient in terms of absence of wave reflection for dissipative schemes, they only give poor intuition for more complex wave propagation phenomena modeled for instance by the Schrödinger or Airy equations. Moreover, our simulations clearly show that local extrapolation boundary conditions for the leap-frog scheme give rise to larger reflected waves than those generated by the boundary conditions we shall construct below. For inflow boundaries, high order local numerical treatments based on the inverse Lax-Wendroff method [16,20,31,33,36] have also been constructed and analyzed in the past. As for extrapolation strategies at outflow boundaries, these numerical treatments may yield good convergence properties. However, they generate wave reflections whose magnitude is larger than what we observe with our nonlocal (in time) procedure. Despite a more involved implementation and a higher computational cost, we thus believe that nonlocal strategies as those examined below may be relevant.
An alternative strategy to the one developed below would consist in starting from the exact transparent boundary condition for the underlying partial differential equation (the continuous problem) and then trying to discretize it. This strategy has been examined, at least, in [2,8,23], and the overall conclusion is that there is no general answer for the construction of such absorbing boundary conditions if one wishes to maintain stability. We therefore rather follow the approach developed in [3,18] and later used in a wide variety of contexts, and thus construct (approximate) transparent boundary conditions that are adapted to the considered finite difference approximation.
Exact Discrete Transparent Boundary Conditions (DTBC in what follows) are analytically available in two space dimensions only when the computational domain is a half-space (Z −1 × N) or the product of a half-line by a torus ((Z/ Z) −1 × N). By "analytically available", we mean that DTBC can be expressed by a welldefined operator, even though its practical implementation may be difficult. We refer for instance to [2][3][4]15] and references therein for various aspects of the theory, the main features of which are recalled below. In either case, the discrete spatial domain is a cylinder of the form × N with no tangential boundary (the variables in ) so DTBC can be derived by performing either a Fourier transform or a Fourier series decomposition with respect to the spatial tangential variables. In the end, this approach yields a one-dimensional problem that is parametrized by the tangential frequency 1 . In that framework, the two-dimensional DTBC thus take the form of a tangential pseudo-differential operator that is intrinsically nonlocal. In order to be implemented on a rectangular computational domain, one therefore has to make the DTBC local, which we achieve here as in [19] by using a suitable "small frequency approximation", which is meaningful (at least formally) for smooth solutions. The frequency cut-off amounts to retaining only finitely many terms in the Taylor expansion of the symbol of the DTBC operator and rewriting the Taylor expansion as a trigonometric polynomial (up to higher order terms). This cut-off strategy is implemented here in such a way that the stencil used for the numerical boundary conditions is not wider than that of the interior numerical scheme. Higher order strategies could be relevant but would require a specific treatment when getting closer to the corner. We do not investigate this issue here.
As in all numerical approximation problems, a crucial feature for the efficiency of the numerical scheme is stability. Following the general result of [15], the DTBC for the leap-frog scheme on a half-space do not meet the strongest possible stability properties because the leap-frog scheme exhibits glancing wave packets (even in one space dimension). On a half-space, the exact DTBC meet some kind of neutral stability which calls for special care since a slight modification in the numerical boundary conditions may yield violent instabilities. In what follows, we show that the zero, first and second order tangential frequency cut-off -as described above -maintain neutral stability when implemented on a half-space (the so-called Godunov-Ryabenkii condition is satisfied [22]). In other words, our approximate transparent boundary conditions are as spectrally stable as the exact ones on a half-space and we therefore feel rather safe to implement them on a rectangular geometry. If instabilities occur in the rectangular geometry, they will necessarily stem from the coupling between the various numerical boundary conditions on each side of the rectangle. A thorough stability analysis on the rectangle seems to be out of reach at the present time. It is not even fully understood at the continuous level! The numerical investigations that we report below show that stability in the rectangle should not be taken for granted. Some coupling strategies exhibit good stability -and convergence -properties while at least one displays violent instabilities. In a future work, we shall extend our strategy to fully discretized hyperbolic systems and to implicit difference schemes for dispersive equations.
Let us now give the plan of our work. In Section 2, we introduce the two-dimensional transport equation and the leap-frog scheme on a rectangular computational domain. The goal is then to construct and analyze approximate transparent boundary conditions. For later use and also for the ease of reading, we first go back to the corresponding one-dimensional problem in Section 3. We adapt the general strategy of [3,18], which was recently generalized in [15], to the particular case of the leap-frog scheme. Some of the tools used in [15] are complemented here by some explicit expansions, based on Legendre polynomials, which were already used in [3,12,18]. Such explicit representation formulae lower the implementation cost of our method and give a sharp asymptotic behavior for the convolution coefficients, which would not be possible by only using the tools in [15]. We then analyze the stability of the one-dimensional DTBC by means of the normal mode decomposition and show that they satisfy the so-called Godunov-Ryabenkii condition [22]. Since the numerical implementation and computational cost of DTBC is sometimes considered as heavy, we then discuss the fast implementation of approximate DTBC by the so-called sum of exponential approach. We relax some of the constraints used in [3,18] for the free parameters in this approximation, and show that the resulting approximation is efficient from a numerical cost point of view and also from the wave reflection magnitude point of view. Section 4 extends this approach to the two-dimensional problem by first considering the exact DTBC on a half-space and then using low frequency truncations in order to derive local (in the tangential space variable) approximate DTBC. We analyze, for the half-space problem, the stability of our approximate DTBC and then report on some numerical simulations on the rectangle.

The transport equation and the leap-frog approximation
In this article, we consider the linear advection equation in two space dimensions: where the space coordinates are denoted ( , ), the velocity reads c = ( , ) and ∇ = ( , ) . It is always assumed from now on that the velocity c is nonzero. To fix ideas, we shall always take and to be nonnegative, but this has no impact on the analysis below. The solution to (2.1) reads: define some integers and (that are meant to be large). Eventually, the time step > 0 will always be chosen so that the parameters are fixed (or at worse vary between two fixed positive constants). These parameters should satisfy some stability requirement, the so-called Courant-Friedrichs-Lewy condition, as explained below. The grid points are defined as ( , ), with: The interior points correspond to 1 ≤ ≤ and 1 ≤ ≤ . The four sides of the rectangle (the numerical boundary) correspond to ∈ {0, + 1} and ∈ {0, + 1}. Letting now , denote the approximation of the exact solution to (2.1) at ( , ) and time , the two-dimensional leap-frog scheme reads: which holds for 1 ≤ ≤ and 1 ≤ ≤ , that is at all interior points. The definition of the numerical scheme requires to prescribe numerical boundary conditions for +2 0, , +2 +1, , = 1, . . . , , and +2 ,0 , +2 , +1 , = 1, . . . , . A crucial observation for what follows is that the update of the numerical solution in the interior does not require to prescribe any numerical boundary condition at the corners (this is because +2 , does not depend on ±1, ±1 at any earlier time index ). In all what follows, the numerical solution ( , ) will therefore be defined for ( , ) ∈ {0, . . . , + 1} × {0, . . . , + 1} ∖ {(0, 0), (0, + 1), ( + 1, 0), ( + 1, + 1)}. We shall be careful in our definition of the numerical boundary conditions to make sure that the four corner values are not involved either.
Let us note that on the whole space Z 2 , the ℓ 2 -stability condition for (2.2) reads: where cfl is a free positive parameter to be chosen. We consider the case , ≥ 0, hence , ≥ 0, so we do not write absolute values in (2.3). The inequality (2.3) gives an upper bound for in terms of , . The bound (2.3), or its one-dimensional analogue that we shall discuss below, is always assumed to hold from now on.
We wish to construct approximate transparent boundary conditions for (2.2) on each side of the rectangle. For ease of reading and also for future use on the two-dimensional problem (2.2), we first go back to the onedimensional problem and briefly recall the construction of DTBC for the one-dimensional leap-frog scheme. We discuss the various possible implementations of the DTBC and approximations by means of sums of exponentials. This technique was used for instance in [3] in the context of the Schrödinger equation and later extended to other contexts. Once the tools have been set in one space dimension, we shall feel free to use them in two space dimensions whenever it makes sense.

Discrete transparent boundary conditions for the one-dimensional
leap-frog scheme

The numerical scheme
We now consider the one-dimensional transport equation: with an initial condition 0 that is supported in an interval [ ℓ ; ]. For concreteness, we shall assume from now on > 0 so that the initial condition is transported towards the right and eventually exits the initial support [ ℓ ; ]. We are interested here in simulating the solution to (3.1) on the fixed interval [ ℓ ; ] by imposing suitable non-reflecting boundary conditions at the boundary points ℓ , .
As in the previous section, we choose a mesh size > 0 such that is an integer. The time step > 0 is chosen in such a way that the parameter defined by: is a fixed positive constant. Keeping the notation := ℓ + , = 0, . . . , + 1, and letting now denote the approximation of the solution to (3.1) at and time , the one-dimensional leap-frog scheme reads The parameter is fixed by imposing: which is a necessary and sufficient ℓ 2 stability requirement for (3.2) on the whole real line ∈ Z, see [22,29]. Since is positive in our framework, we thus fix ∈ (0, 1) and let , vary accordingly. The numerical sheme (3.2) on the interval 1 ≤ ≤ requires a definition for the boundary values +1 0 , +1 +1 , and a definition for the initial data ( 0 ), ( 1 ). The initial condition ( 0 ) is defined by setting ∀ = 0, . . . , + 1, where 0 is the initial condition for (3.1). In the numerical simulations reported below, we have chosen [ ℓ , ] = [−3, 3] and 0 ( ) = exp(−10 2 ). The first time step value ( 1 ) is defined by imposing the second order Lax-Wendroff scheme in the interior domain, namely: The boundary values 1 0 and 1 +1 are set equal to zero for simplicity (which is almost exact for the Gaussian initial condition that we consider). We now recall how to derive the DTBC for the scheme (3.2).

Derivation of DTBC
The goal of this paragraph is to derive the DTBC for the leap-frog scheme (3.2). We follow the methodology of [18], which was originally designed for the Schrödinger equations, and therefore consider the numerical scheme: on the whole real line Z, assuming that the initial conditions ( 0 ) ∈Z , ( 1 ) ∈Z belong to ℓ 2 and vanish outside of the interval {1, . . . , }. The stability condition (3.3) is assumed to hold. The DTBC are first computed in terms of the so-called -transform of the sequences ( ) ∈N . Let us briefly recall the definition of the -transform (we do not pay too much attention here to the convergence of the Laurent series involved in the computations below and rather remain at the level of formal series).
Definition 3.1. Let ( ) ∈N ∈ C N be a sequence. The -transform of ( ) ∈N , which is denoted̂︀, is defined by:︀ , the series being well-defined and holomorphic in { ∈ C, | | > } where 1/ is the convergence radius of ∑︀ ∞

=0
(here is a placeholder for 1/ and the radius may be zero or infinite).
We now letˆ( ), ∈ Z, denote the -transform of the sequence ( ) ∈N : It is shown in [15] that the above series converges for | | > 1 thanks to the ℓ 2 -stability of the leap-frog scheme 3 . Moreover, the sequence (ˆ( )) ∈Z is square integrable. We apply the -transform to the numerical scheme (3.4) and obtain: for ≤ 0 and ≥ + 1. (On the interval {1, . . . , }, the initial conditions give rise to a nonzero source term on the right hand side of (3.5).) The equation (3.5) is a recurrence relation and we therefore look for the solutions to the characteristic equation: The two roots to (3.6) are given by with the standard definition of the square root (the branch cut is R − ). When | | → +∞, | − | → +∞ and | + | → 0. Moreover, ± do not belong to S 1 for | | > 1 so we have | + | < 1 and | − | > 1 for | | > 1. We thus relabel these two roots as a stable and unstable one and write from now on: The reason for the superscript 0 will be clear in the next section when we deal with the two-dimensional problem. The recurrence relation (3.5) and the decay at infinity of (ˆ( )) ∈Z implies: It now remains to compute the Laurent series expansion of both 0 and 1/ 0 and to perform the inversetransform in order to write the numerical boundary conditions (3.9) in the original discrete time variable .
Since we have 0 ( ) 0 ( ) = −1, see (3.6), we can equivalently rewrite (3.9) as: There are at least two ways to compute the Laurent series expansion of 0 for the leap-frog scheme, which we now detail in order to compare the possible benefits and drawbacks of each method. The first method is rather in the spirit of [18] while the second one follows [11,15].
An explicit expansion in terms of Legendre polynomials. The first method to obtain this expansion relies on the explicit formula (3.8) and on the expansion where denotes the th-Legendre polynomial [32]. This family of polynomials can be equivalently defined by the first two terms 0 ( ) := 1, 1 ( ) := and by the recurrence formula: Starting from (3.8), we have with := 1 − 2 2 ∈ (−1, 1). We therefore derive the expansion We focus below on the right boundary condition at point +1 = . Going back to (3.10) and using the above expansion (3.13) for 0 , the DTBC (3.10) reads and it is understood that −1 vanishes (in the expression of 0 ). We compute 0 = − , and thanks to the recurrence relation (3.12), we get the shorter expression and finally Recall that in (3.14), we use the notation = 1 − 2 2 . We now look for an efficient way to implement the coefficients of the above series without using the expression of the Legendre polynomials (which would be computationally rather cheap anyway). Let us define the sequence ( 0 ) ≥0 such that 0 0 := , 0 1 := (1 − 2 ), and 0 := ( −1 ( ) − +1 ( ))/((4 + 2) ) for ≥ 2. After various simplifications using (3.12), we obtain The DTBC (3.14) therefore readsˆ+ Performing the inverse -transform, we get 4 : The recurrence relation (3.15) may be efficiently implemented in a computer code, which gives access to the coefficients in the numerical boundary condition (3.17) up to any prescribed final time index . Going back to (3.10), we also derive the "left" numerical boundary condition: In view of understanding the stability properties of the numerical boundary condition (3.18), it might be important to determine the asymptotics of the sequence ( 0 ) ∈N . Let us recall that 0 equals ( −1 ( ) − +1 ( ))/((4 + 2) ) for ≥ 2 (here ∈ (0, 1) and = 1 − 2 2 ). We use the so-called Laplace formula for the Legendre polynomials, see Theorem 8.21.2 of [32]: where the remainder term is even uniform with respect to on every compact set of the form [ , − ], > 0. We fix the angle ∈ (0, ) such that cos = = 1 − 2 2 , and therefore sin = 2 √︀ 1 − 2 . We then apply the Laplace formula and obtain after a few simplifications: This asymptotic behavior can be verified on numerical experiments.
Determining inductively the expansion. As observed in [11,15], the Laurent series expansion of 0 can also be computed inductively by using the equation (3.6). The main benefit is that the method does not rely on an explicit knowledge of a power series expansion such as the one involving the Legendre polynomials (which is useful only for three point schemes). The methodology below is therefore more flexible in view of being generalized to numerical schemes with larger stencils (as the one considered in [11]). We know that 0 tends to zero at infinity so we may plug its Laurent series expansion 5 ∑︀ ≥1 0 − into (3.6) and obtain which is found to be equivalent to 0 1 = and

21)
4 It is understood in (3.17) that the summation holds over all integers between 0 and ( + 1)/2. In case is even, then the summation holds over all integers between 0 and /2. 5 It is proved in [15] that the Laurent series is convergent for | | > 1 which follows from the splitting of the two roots of (3.6) for | | > 1.

S543
where we use the convention 0 0 = 0. It is not difficult to see that if is even, then 0 = 0. Let now = 2 + 1 be odd, ≥ 0. An easy computation gives This means that the sequence ( 0 ) ∈N appearing in (3.17) and that satisfies the recurrence relation (3.15) also satisfies with 0 0 = ( 0 coincides with 0 2 +1 ). Observe that the value 0 1 = (1 − 2 ) that we have found in the previous paragraph is consistent with (3.22). It seems less clear to derive the asymptotic behavior (3.20) by starting from (3.22) rather than from the more explicit representation (3.14). Both approaches thus seem to have their own interest.

Stability analysis on a half-line
In this short paragraph, we explain why the normal mode analysis predicts "neutral stability" for any of the two half-line problems where the leap-frog scheme (3.23a) is considered on the half-line { ≤ } or { ≥ 1}, in combination with either (3.23b) or (3.23c). Let us focus on the case { ≥ 1}, since the other case is entirely similar.
The normal mode analysis consists in determining the solutions of the leap-frog scheme (3.23a) of the form = , with | | > 1, and ( ) ∈ ℓ 2 . Among such sequences, the final question is to determine whether there exist nonzero ones that satisfy the numerical boundary condition (3.23b). Following [22], we shall say that the Godunov-Ryabenkii condition is satisfied if there is no nonzero sequence ( ) ∈ ℓ 2 such that = is a solution to the leap-frog scheme that satisfies (3.23b). Otherwise, we shall say that is an unstable eigenvalue. Of course, unstable eigenvalues preclude (in the most violent way) stability estimates hence a desirable feature of any numerical boundary condition is to satisfy the Godunov-Ryabenkii condition.
Plugging the ansatz = in the leap-frog scheme, we find that the sequence ( ) should belong to ℓ 2 and satisfy ∀ ≥ 1, The computation of the ℓ 2 -solutions to this recurrence relation follows by solving the characteristic equation (3.6) and by determining among the roots to (3.6) which have modulus less than 1. This classification has already been performed in the previous paragraph, so for | | > 1, we can conclude that the sequence ( ) is given by ∀ ≥ 1, = 0 0 ( ) , 0 ∈ C. Inserting into (3.23b), we have to determine whether among such sequences (that are parametrized by their initial state 0 ), we can have Taking the limit → ∞, we thus need to determine whether there can hold where the last equality follows from the expression of 1 . In other words, we need to determine whether there can hold 0 ( ) = 0 ( ) for some complex number with | | > 1. This is clearly impossible since 0 ( ) has modulus < 1 and 0 ( ) has modulus > 1. This means that the so-called Godunov-Ryabenkii condition holds (see [22]): there does not exist any unstable eigenvalue for the half-space problem { ≥ 1} with the numerical boundary condition (3.23b). Let us observe however that the roots 0 and 0 coincide when equals one of the four values for which the characteristic equation (3.6) has a double root. These values of correspond to the glancing spatial frequencies (for which the associated group velocity vanishes). This is a neutral stability case, whose continuous counterpart has been analyzed in details in Chapter 7 of [10], see also [35] for other examples of this situation.

Fast implementation of approximate DTBC with sums of exponentials
The computation of the DTBC (3.23b) and (3.23c) at nodes ℓ and and time +2 = ( +2) requires /2 sums and /2 multiplications. Therefore, the total cost of discrete convolution computations for a full simulation up to time = is ( 2 ). Since is related to by the CFL condition (3.3), a fine space grid will make the time step small and as a consequence the number of time steps large. Since the leap-frog scheme (3.2) is explicit, each time iteration requires ( ) operations for interior nodes , = 1, . . . , , which gives ( ) operations for a full simulation. Thus, the total number of operations is ( 2 + ). If > , which corresponds to large time simulations, the main part of the computational cost is related to the computation of discrete convolutions. Moreover, it is necessary to store in memory the evolution of the solution at the two boundary nodes. It may therefore be useful to reduce this part of the computation by applying a more "local" formula.
In [3,17], a fast convolution procedure to compute approximation of discrete convolutions ( ) := ∑︀ =0 − was introduced. We assume that the convolution coefficients ( ) ∈N are given and satisfy a decay property similar to (3.20) that allows to derive their approximation by "sums of exponentials". In order to present the idea and the efficiency of the method, let us assume for a moment that we are able to compute an approximation (˜) ∈N of the convolution terms ( ) ∈N given by where the integer , the coefficients and the complex numbers , all of them satisfying | | > 1, have to be defined. Denoting˜( ) := ∑︀ =0˜− , we compute The efficiency of the method relies on the derivation of a recurrence relation to compute the terms ( ) ( ): Assuming 0 = 0, the initial value for each ( ) ( ) is given by (0) ( ) = 0. It should be understood that in a practical implementation, each value ( ) ( ) is stored in memory so the recurrence (3.25) represents only three operations at each time iteration. Accordingly, the original discrete convolution ( ) is approximated by and we have replace the computation of a nonlocal discrete convolution that involves operations by the sum of terms computed by the explicit local operations (3.25). In this case, the total number of operations is ( ( + )). It remains to derive the approximation (3.24) of the sequence ( ) ∈N by (˜) ∈N . The sum-of-exponential formula (3.24) is frequent and various techniques are available to compute the coefficients ( ) and ( ) (see e.g., [13,14]). We present here a modified version of the algorithm derived in [3]. Let us consider the power series ( ) := ∑︁ ∈N , for | | < 1, associated with the original sequence ( ) ∈N (which we assume to satisfy a polynomial bound of the form | | , ∈ R, so the power series converges on the unit disk). We then consider its [ , ] Padé approximant, see [5,7], which we denote˜( ) = ( )/ ( ), where and are respectively polynomials of degree and . We assume to be unitary in order to normalize the polynomials. Here we shall always consider the case < . The power series expansion of˜at 0 is given bỹ with˜= for 0 ≤ ≤ + (by definition of the Padé approximants). Let us now assume that has simple roots with | | > 1 for 1 ≤ ≤ . So, we have and thus Let us now define the coefficients S546 C. BESSE ET AL.
We are going to show that the polynomial −1 in the latter equality equals . Indeed, we have Since is of degree ≤ − 1, by uniqueness of the interpolating polynomial, we have −1 = . As a partial conclusion, we havẽ where the coefficients are defined in (3.27). Writing − = (1 − / ), and using | | > 1, then for | | ≤ 1 < | |, we obtain Let us recall that the power series expansion of˜at the origin reads˜( ) = ∑︀ ≥0˜s o proceeding by identification, we have and the coefficients˜satisfy the propertỹ The computation of (˜) ∈N requires to determine the [ , ] Padé approximant of , the roots of (praying for them to be simple) and the coefficients given by (3.27).
In [3], the authors suggested to take = − 1. Unfortunately, when we apply the latter approximation strategy to the function 0 (1/ ), we are not always able to verify the assumption inf | | > 1 on the roots of , even when these roots are computed with (very) high accuracy. Our choice to take "free" and not necessarily equal to − 1 allows us to find Padé approximants for which we can verify inf | | > 1. We present below some experiments to emphasize this remark. Remark 3.3. In [3], the authors justified the convergence of the approximate coefficients˜to as → ∞. Their result was based on the so-called Baker-Gammel-Wills conjecture, which is unfortunately now known to be false, see [6] and references therein. That does not mean however that the above procedure is meaningless, as the numerical results in [3] and those we present below confirm. The main difficulty in applying the above approach is that there does not seem to be any general result available to make sure that the condition inf | | > 1 holds in a systematic way (not mentioning the condition that the roots should be simple. . . ).
In [3], a Maple code is proposed to calculate the Padé approximant of . The use of Maple insures to not be limited to double precision computations. This is very important to compute accurately the coefficients of the monomials of and , and accordingly to determine the roots of with high accuracy. This process is however quite limited since it relies on formal calculus softwares and is also restricted to small values of (typically ≤ 20). We therefore present below an alternative self-contained procedure which allows us to consider much larger values of (say, ≤ 200) provided that is not close to (say ≤ 50). When is close to , the coefficients of and become extremely large and the numerical computations become rather unstable. The [ , ] Padé approximant˜of reads where we fix for instance q 0 = 1, and we restrict to the case < . (The normalization convention for is not the same as above but this does not affect the previous arguments.) We therefore have to compute + +1 unknown coefficients, which are solutions of the + + 1 equations given by where ( ) denotes the th derivative with respect to . Indeed, thanks to the property (0) ̸ = 0, (3.28) can be seen to be equivalent to the property which defines the [ , ] Padé approximant of . The equations (3.28) reduce to and min( , ) The first equation in (3.29) leads to p 0 = 0 = (0), and the remaining equations in (3.29) and (3.30) can be recast as a linear system: where ∈ M + , + (R) and , are vectors in R + . The matrix in (3.31) is given by (the first row and column below present the indices for the sake of clarity): The unknown in (3.31) is = (q 1 , . . . , q , p 1 , . . . , p ) , and the right hand side is = (− 1 , − 2 , . . . , − + ) . The linear system (3.31) thus displays block matrices The matrix is a subtriangular Toeplitz matrix. Its inverse (assuming that 0 is nonzero) is also a subtriangular Toeplitz matrix whose first column is easy to compute by solving the linear system = 1 where 1 is the first vector of the canonical basis of R . The other columns of −1 are deduced from the Toeplitz structure.
Using the Schur complement method, we can compute the unknowns by solving successively and then q = −1 ( 1 − p).
In order to get high accuracy both in the computations of p and q, but also of the roots of , we use the floating-point arithmetic with arbitrary accuracy Python library mpmath. The advantage of using the library mpmath is its ability to control the accuracy of floating-point arithmetic. For all the numerical experiments in this paper, we set the decimal accuracy to 80, meaning the library uses 266 bits to hold an approximation of the numbers that is accurate to 80 decimal places.
Going back to our numerical scheme (3.23), we wish to approximate the sequence ( 0 ) given by (3.15) by a sum of exponentials. We apply the above algorithm, which is legitimate since 0 0 ̸ = 0 ( 0 plays the above role of for all ). We present below a comparison between the thousand first elements of the sequence ( 0 ) in the complex plane, compared with the unit circle, and the evolution of˜0 compared with that of 0 for 0 ≤ ≤ 1000. We first choose = 50 with = 10 (see Fig. 1), then = 50 with = 49 (see Fig. 2), then = 100 with = 30 (see Fig. 3), and eventually = 100 with = 99 (see Fig. 4). On Figures 2 and 4, one root of the polynomial is relatively big compared with the other roots and we therefore do not represent it. They respectively take the values −4.6 10 17 and −1.2 10 29 , which means that some of the coefficients of are so large that any computation involving the roots of should be taken with much care. We see that the more is close to , the more˜0 are close to 0 for a large interval of values of . The quality of the approximations for ( , ) = (50, 49) and ( , ) = (100, 30) is equivalent ; the case ( , ) = (100, 30) is however much more stable since all roots of the polynomial remain within a "small" ball of the complex plane while has a very large root in the case ( , ) = (50, 49). When ( , ) = (100, 99), one of the roots that we computed for has a modulus less than one, which generates blow up for sequence (˜0 ) (see Fig. 4). This does not prove however that has a root of modulus less than 1, because the computation is so unstable (due to the very large coefficients in ) that the final result should be taken with care even with the decimal accuracy set to 80.

Numerical experiments
We present here some numerical tests with the following parameters: The number of grid points is +1 = 1000, and the leap-frog scheme (3.2) is implemented with a CFL parameter = 5/6 (the time step is computed accordingly). Since the initial condition is roughly concentrated in the interval [−1, 1], the exact solution to the transport equation more or less vanishes on the interval [ ℓ , ] after the time 4. We measure below the amplitude of the reflected wave at the right boundary depending on the choice of numerical boundary conditions. We present below the evolution of the logarithm log 10 | | rather than the evolution of the solution itself. This gives a much more visible representation of wave reflections and their magnitude. Let us observe that in  Figures 5 and 6 below, we always observe a left going wave emanating from the initial condition, this wave being highly oscillatory (in both space and time) and having magnitude 10 −8 . This is due to the fact that the leap-frog scheme (3.2) supports the wave (−1) + which has group velocity − and thus travels backwards [34]. In Figure 5, we first compare the leap-frog scheme with exact DTBC (3.23) (left) with the leap-frog scheme implemented with the more classical Neumann type boundary condition +2 +1 = +1 , see [35] (accordingly the numerical boundary condition at 0 is +2 0 = +1 1 ). After the theoretical exit time for the continuous solution, the DTBC strategy leaves a solution of amplitude 10 −16 in the whole domain while the easier (local) Neumann strategy exhibits multiple wave reflections.
We then show on Figure 6 the evolution of the logarithm log 10 | | when we do not implement the exact DTBC (3.23b) and (3.23c) but rather their approximation by sums of exponentials as explained in the preceding paragraph. The numerical simulation with the parameters ( , ) = (50, 6) (left of Fig. 6)    (right of Fig. 6) for the degrees of the Padé approximant. In each of these two cases, the condition that all roots of are simple and lie outside the unit circle has been checked numerically. We observe that in the numerical boundary conditions (3.23b) and (3.23c), the determination of +2 0 and +2 +1 relies on the values 1 and for either only odd or even values of . When we implement the approximate numerical boundary condition (3.26), we thus need to consider separately the cases where is either odd or even, and we therefore introduce two sequences (︁ )︁ and (︁ (2 +1) )︁ for each root of the polynomial . The numerical results presented in Figure 6 exhibit very low reflection amplitudes. In the case ( , ) = (50, 49), the overall accuracy is comparable with the exact DTBC. Note however that the condition inf | | > 1 for this approximation to make sense is not granted and the implementation of the sum of exponential approximation requires first solving for and and then identifying the roots of (which is not necessarily easier than implementing (3.15) and (3.23)).

Exact DTBC on a half-space and local approximations
In this paragraph, we go back to the leap-frog scheme (2.2) for the linear advection equation (2.1) in two space dimensions. Our first goal is to derive the exact DTBC for (2.2) when that scheme is considered on a half-space, be it { ≥ 1}, { ≤ }, { ≥ 1} or { ≤ }. This means that we shall consider the leap-frog scheme (2.2) in the whole plane ( , ) ∈ Z 2 with initial data supported, say, in the half-space { ≥ 1}, and we shall derive the numerical boundary conditions satisfied by the solution to (2.2) at = 0. This is exactly the problem considered in [15] in the general framework of multistep finite difference schemes. For ease of reading and in order to stick to the notation of Section 3, we consider from now on the case where the initial data for (2.2) are supported in { ≥ 1}, and it is understood that in (2.2) the tangential variable lies in Z. We consider the step function (·) defined by ( ) := , , ∀ ∈ [ , +1 ), ∈ Z, with = + for all ∈ Z and not only for = 0, . . . , + 1. Then (2.2) reads so applying first a partial Fourier transform with respect to and then the -transform in time, we are led (with obvious notation) to the recurrence relation which holds for all ≤ 0 (in the half-space { ≥ 1} where the initial data are supported, the analogous recurrence relation displays a nonzero source term on the right hand side). The solution to the recurrence relation (4.1) is computed by considering the characteristic equation: where is a placeholder for the rescaled frequency . For any real number and | | > 1, the roots to (4.1) do not belong to S 1 . Therefore one has modulus < 1 and the other one has modulus > 1 (see the general splitting argument in [15]). We label the two roots in such a way that: whenever | | > 1 and ∈ R. The discrete boundary condition (on the Fourier--transform side) thus readŝ which is reminiscent of (3.9). The analogous numerical boundary condition at = + 1 (for the other half space problem) isˆ+ If we now apply the inverse Fourier--transform to (4.3), this would lead to the exact, nonlocal DTBC for (2.2) on a half space. Nonlocality refers here both to time and to the tangential variable . In other words, determining +2 0, requires the knowledge of 1, ′ for all ′ ∈ Z and at all time steps up to +1. The nonlocality in time is not so harmful from a computational point of view (as we have seen in one space dimension) but the nonlocality with respect to is much more problematic since the set to which belongs is infinite. We therefore propose to follow the idea developed by Engquist and Majda [19] for continuous problems and to localize the exact DTBC (4.3) in the -direction. This is performed by means of an asymptotic expansion with respect to = in (4.3). At a formal level, this amounts more or less to assuming that the discrete solution displays only bounded tangential frequencies and to performing an asymptotic expansion with respect to (which is meant to be small in practice).
Let us therefore go back to (4.2) and observe that again the product of the two roots and equals −1.
and we now expand the stable root ( , ) to (4.2) with respect to . Up to the third order, this expansion reads where 0 ( ) is the stable root for the one-dimensional problem (which corresponds to = 0 in (4.2)) and the choice of writing sin rather than , as well as sin 2 /2 rather than 2 , has been made in order to exhibit amplification factors that are obviously linked with tangential finite difference operators (the centered first order derivative and the discrete Laplacian). At this stage, localizing the exact DTBC (4.3) with respect to amounts to using in (4.3) finitely many terms of the Taylor expansion (4.4) (e.g., either the first term, or the two/three first terms). Namely, on the Fourier--transform, our localization procedure for (4.3) yields one of the three following choices (in increasing order of approximation): After performing an inverse Fourier transform with respect to , (4.5) reads: where the hat notation refers here to the -transform only. In order to perform the inverse -transform in (4.6) and write down the approximate non-reflecting boundary conditions in the physical variables, we need to compute the Laurent series expansion of the functions 1 , 2 that appear on the right hand side of (4.6b) and (4.6c).
(The Laurent series expansion of 0 has already been derived in the analysis of the one-dimensional problem.) This is achieved below with either of the two methods that we have already used in the one-dimensional case.
Expansions based on Legendre (and Tchebychev) polynomials. Let us recall that we have already determined the Laurent series expansion of the function 0 in (4.6). We have written it under the form where the sequence ( 0 ) ≥0 is determined by either (3.15) or (3.22). The two first values are 0 0 = and 0 1 = (1 − 2 ). Let us eventually recall that the expression of 0 ( ) is given by (3.8). We are now going to determine the Laurent series expansion of 1 ( ). The root ( , ) to (4.2) is simple for ∈ R and | | > 1. It thus depends holomorphically on and is ∞ with respect to . Differentiating (4.2) with respect to at = 0, and identifying 1 = ( , 0)/(2 ), we obtain which yields by using (3.8) and (3.11). We can therefore write the Laurent series expansion of 1 under the form: with 1 0 = 0 and ∀ ≥ 1, We can, for instance, compute (recall 1 0 = 0): The expression (4.10) is only seemingly singular with respect to for ≪ 1. Indeed, we recall = 1 − 2 2 and since all Legendre polynomials satisfy (1) = 1, the expression ( ) − −1 ( ) can be written as 2 ( 2 ) for some polynomial . We therefore see from (4.10) that all coefficients 1 read ( 2 ) for some real polynomial . Using again the Laplace formula (3.19), we also get the asymptotic behavior: where the angle ∈ (0, ) is still defined by cos = 1 − 2 2 . The asymptotic behavior (4.11) can be verified by numerical experiments. Comparing with (3.20), we observe that the sequence ( 1 ) has a slower decay than ( 0 ), but it decays to zero nevertheless. This is due to the fact that 1 has (four) singularities of the form ( / 0 − 1) −1/2 , 0 ∈ S 1 , as approaches the unit circle S 1 , while 0 is continuous up to S 1 .
The Laurent series expansion of 2 can be derived by following more or less the same lines. By differentiating twice (4.2) with respect to and plugging the expression (4.7) of 1 , we first get: In other words, we have derived the relation and it only remains to use (3.11) as well as the generating function of the Tchebychev polynomials of the second kind (see [32]): Writing the Laurent series expansion of 2 under the form: we have derived the relations 2 0 = 0 and We can for instance compute It does not seem so easy to infer from (3.19) and from the relation (cos ) = sin(( + 1) ) sin , the asymptotic behavior of ( 2 ) ∈N . However, we verify below on some numerical experiments that ( 2 ) ∈N does not seem to be bounded. It actually seems to grow like √ , see Figure 7, up to an oscillating behavior similar to the one we have exhibited in (3.20) and (4.11). This growth in √ is consistent with the singularities of the form ( / 0 − 1) −3/2 displayed by 2 . The rigorous justification of this asymptotic behavior is left to a future work. . This strategy can be extended to the "correctors" 1 and 2 . Namely, the function 1 is given by the relation (4.7). Since 0 ( ) has a finite limit at infinity, we already see on (4.7) that 1 tends to zero at infinity. It actually decays like −2 . Hence its Laurent series expansion reads , and the coefficients ( 1 ) ≥1 are computed by plugging the latter expression in (4.7) (which itself has been obtained by differentiating (4.2) with respect to ). After a few simplifications, we obtain 1 1 = 0 and the recurrence relation where we use the convention 1 0 = 0. We easily deduce that all odd coefficients 1 2 +1 vanish (recall that the even coefficients 0 2 vanish), and the even coefficients 1 2 , which we also write 1 to be consistent with (4.9), satisfy 1 0 = 0 and Recalling the initial values 0 0 = and 0 1 = (1 − 2 ), we recover for instance from (4.16) the two first values 1 1 = − and 1 2 = − (2 − 3 2 ) as in the previous method. We can follow the same lines to derive the Laurent series expansion of 2 , by starting from the relation (4.12) and plugging (observing on (4.12) that 2 tends to zero at infinity): which yields 2 1 = 0 and the recurrence relation
Independently of the method we choose to compute the coefficients in the Laurent series expansion of 1 and 2 , we can implement them numerically up to any prescribed final time index . We now go back to (4.6). Once we have performed the inverse Fourier -transform, the approximate DTBC (4.6) respectively read: depending on the tangential accuracy we aim at achieving on the boundary. Of course, the asymptotic expansion that we have derived for ( , ) also provides the approximate DTBC at the right boundary = + 1: Let us observe that the boundary conditions (4.17) and (4.18) are nonlocal in time but they are local in space (unlike the original exact DTBC). The boundary conditions (4.17) and (4.18) only involve a three point stencil, as depicted on Figure 8 (left picture) on the example of the point located closest to the upper right corner. We could have pushed the Taylor expansion with respect to further by including for instance the following term 3 . The problem is that there is no linear combination of 1, cos and sin that behaves like 3 at the origin. In other words, if one writes 3 as a trigonometric polynomial ( ) up to an ( 4 ) term, then the trigonometric polynomial has degree at least 2, which means that the associated finite difference operator in will involve at least 4 points (and more likely five if one wishes to keep some symmetry). Hence this extension causes some trouble when getting close to the four corners of the rectangle because the very last boundary points close to the corner would require a specific treatment. We therefore do not pursue this higher order extension for the time being and leave it to future investigations.
At this stage, we have made the numerical boundary conditions on the "left" and "right" boundaries, resp. { = 0} and { = + 1}, explicit. The analysis can be carried out almost word for word for the "bottom" and "top" boundaries, namely { = 0} and { = + 1}. The only thing is to observe that the indices , in (2.2) can be exchanged by simply switching the roles of and . From now on, we let ( 0 ) ≥0 , ( 1 ) ≥0 , ( 2 ) ≥0 denote the sequences obtained by the same procedure as ( 0 ) ≥0 , ( 1 ) ≥0 , ( 2 ) ≥0 , but switching the roles of and . For instance, the sequence ( 0 ) ≥0 is defined by 0 0 = , 0 1 := (1 − 2 ), and Then the approximate DTBC that we consider on the bottom and top boundaries read Before reporting on several numerical simulations using the boundary conditions (4.17)-(4.20), we first study whether each subcase in these boundary conditions, for instance (4.17a)-(4.17c) satisfy some basic stability requirements. This stability analysis for half-space problems is carried out in the following paragraph.

Stability analysis
The stability analysis of initial boundary value problems for hyperbolic systems in regions with corners is a delicate matter, see for instance [9,27,28,30] and references therein for a partial account of the theory (of which much, if not all, remains to be done). The discrete counterpart is not easier and has been left (almost) undone so far. At a theoretical level, we do not claim to study the interaction of either one of the numerical boundary conditions in (4.20) with either of the conditions in (4.18) at the upper right corner of the rectangle. Such an analytical study seems to be out of reach at the present time. However we shall report below on some interesting numerical observations which call for a lot of care in the coupling of numerical strategies on each side of the rectangle. Namely, we shall show on one example that coupling two specific numerical boundary conditions on the top and right boundaries, namely (4.20c) and (4.18c) yields strong instabilities (even though the two separate half-space problems seem to be stable, as our partial proof below indicates). This phenomenon is well-known in the PDE context, see an elementary example in [26] or more elaborate examples in [9,30].
We now examine a more simple problem which consists in studying the stability of the half-space problem where the leap-frog scheme (2.2) holds in the half space {( , ) ∈ Z 2 , ≥ 1} and either of the numerical boundary conditions (4.17a)-(4.17c) is imposed on the boundary { = 0}. Following what we have done in the one-dimensional case, the stability analysis is carried out here by means of the so-called normal mode analysis. Extending slightly what we have done in the one-dimensional case, the normal mode analysis amounts here to determining the solutions to the leap-frog scheme (2.2) which are of the form: with | | > 1, ∈ R and ( ) ≥0 ∈ ℓ 2 , that also satisfy the homogeneous numerical boundary conditions, be they (4.17a), (4.17b) or (4.17c). The arguments are quite similar in each of the three cases (except that we shall not be able to carry them out completely in the last case). We deal below with (4.17a)-(4.17c) in increasing order of complexity.
The "zero order" numerical boundary condition (4.17a). First of all, we plug the expression , = exp( ) in (2.2) and find that ( ) ≥0 must satisfy the recurrence relation (Compare with (4.1).) Since ( ) ≥0 should decay to zero at infinity, we get where we recall that ( , ) is the only root of modulus < 1 to The "first order" numerical boundary condition (4.18b). We can use part of the previous analysis, namely the determination of the "stable" solutions to (2.2) of the prescribed normal mode form. However we now consider the numerical boundary condition (4.18b) or rather its Fourier--transform counterpart (4.5b), so we need to determine whether there holds or equivalently ( , ) = 0 ( ) + 2 sin 1 ( ), (4.21) for some pair ( , ) with | | > 1 and ∈ R. Let us recall that 0 is a root to (3.6) and that 1 is given by the relation (4.7). Let us first consider the case = 0 ( ), that is, sin = 0. Then ( , ) = 0 ( ) and (4.21) can not hold since the left hand side has modulus > 1 while the right hand side has modulus < 1. We therefore assume from now on sin ̸ = 0. Let us assume for a moment that the equality (4.21) holds for some pair ( , ) with | | > 1 and sin ̸ = 0. Then the right hand side of (4.21) must be a solution to the second degree polynomial equation (4.2). Plugging the right hand side of (4.21) in (4.2) and expanding, we eventually get 6 : Simplifying by sin , we get either 1 ( ) = 0 or + 1 ( ) = 0. Let us exclude the first of these two options. For | | > 1, the stable root 0 ( ) to (3.6) does not vanish. Hence we see on the relation (4.7) that 1 ( ) does not vanish either. This means that if (4.21) holds for some pair ( , ), then there necessarily holds + 1 ( ) = 0. We use again (4.7) to compute: and the latter quantity can vanish only if = 0 (recall that we have assumed ̸ = 0). At this point of the analysis, we have shown that the only possibility for (4.21) to hold is to have = 0, that is = 0. However, in that case, we have ( , ) = 0 ( ) for all , see (4.2), and 1 ( ) = 0 (see (4.7)). Hence (4.21) can never hold because for = 0, the right hand side has modulus < 1 while the left hand side has modulus > 1. In other words, the Godunov-Ryabenkiii condition holds and there does not exist any unstable eigenvalue for the half-space problem { ≥ 1} with the numerical boundary condition (4.17b). Of course, the same arguments would apply for any of the other half space problems.
The "second order" numerical boundary condition (4.18c). We can use the same approach as above. Verifying the Godunov-Ryabenkii condition amounts to determining whether there exists a pair ( , ) with | | > 1 and ∈ R such that: ( , ) = 0 ( ) + 2 sin 1 ( ) − 4 sin 2 2 2 ( ). (4.22) We verify again that (4.22) can not hold if = 0 ( ) or if = 0. Hence we assume from now on sin ̸ = 0 and ̸ = 0. Let us assume that the relation (4.22) holds, meaning that not only the right hand side of (4.22) is a root to (4.2) but also that it is the only root of modulus > 1. We first plug the right hand side of (4.22) in (4.2) and collect the terms in increasing powers of sin /2 (we use again the equation (3.6) satisfied by 0 , the relation (4.7) satisfied by 1 and we also use the relation (4.12) satisfied by 2 ). After simplifying (quite a bit. . . ), we eventually get the polynomial equation Since sin is nonzero, then sin /2 is also nonzero, and (4.23) is consequently an eight degree polynomial equation in . Changing into , (4.23) becomes a polynomial equation with real coefficients (after simplifying by ). Moreover, the roots of (4.23) are nonzero and are invariant by the transformation → −1/ . Let us also recall that we are interested in determining the roots of (4.23) that satisfy | | > 1.
We have not been able to obtain a complete proof of the facts we claim below, hence we do not have a complete proof for the verification of the Godunov-Ryabenkii condition, but repeated numerical experiments indicate the following facts: -For ∈ (0, max ), where the angle max is defined by (4.23) has a unique root verifying | | > 1. This root is a purely imaginary number (one other root of (4.23) is −1/ since the roots of (4.23) are invariant under the transformation → −1/ , and the six other roots have modulus one).
-For ∈ [ max , ), all roots to (4.23) have modulus one (the value max is obtained by considering the case where − is a (double) root to (4.23)).
For (4.23) to have a root that satisfies | | > 1, we must necessarily have ∈ (0, max ), and in that case the remaining question is to determine whether the quantity 0 ( ) + 2 sin 1 ( ) − 4 sin 2 2 2 ( ) is either the stable or the unstable root to (4.2). On the numerical experiments we have conducted, we have verified the property: for ∈ (0, max ), meaning that the right hand side of (4.22) may coincide with a root (4.2) for some specific values of ( , ) but in that case it coincides with the stable root ( , ) and not with the unstable root ( , ) (which would complete the verification of the Godunov-Ryabenkii condition). The inequality (4.24) can be shown rigorously in the regime where is small, for in that case we have → ∞, see (4.22), and all functions 0 , 1 , 2 tend to zero at infinity, so the modulus of 0 ( ) + 2 sin 1 ( ) − 4 sin 2 2 2 ( ) tends to zero as tends to zero. We can also verify (4.24) in the limit regime → max for in that case we get → − and we compute We can then verify numerically that the inequality holds for all positive values of , such that + < 1 (here we use the classical formula sin max = 2 /(1+ 2 ) and sin 2 max /2 = 2 /(1 + 2 ) with := tan max /2). Hence it does seem that (4.22) can not hold for some pair ( , ), which means that the Godunov-Ryabenkii condition is verified (though a complete analytical proof of this fact still remains open). We postpone this rigorous justification to a future work.

Numerical experiments on a rectangle
We go back to our original motivation which is the simulation of the leap-frog scheme (2.2) on the rectangle {1 ≤ ≤ , 1 ≤ ≤ } with some approximate DTBC on each side. In the numerical simulations reported below, the computations were run on the rectangle ( , ) ∈ (−3, 3) × (−2, 2) and on the time interval ∈ [0, 8].
The boundary values for 1 (all other relevant values of , ) are set equal to zero for simplicity. In the numerical simulations below, we try various coupling strategies of numerical boundary conditions on the "right" and "top" boundaries { = + 1} and { = + 1}. We first make it clear on Figure 8 that the two sides of the upper right corner are coupled through the boundary conditions if one chooses for instance to impose (4.18b) on = + 1 and (4.20b) on = + 1. The coupling is weak though since the boundary conditions (4.18b) and (4.20b) are explicit so the trace values +2 +1, and +2 , +1 are coupled only through the preceding time steps.
In order to highlight the wave reflections on each side of the rectangle, we show below on some two and three dimensional figures the logarithm of the (absolute value of the) solution. The logarithm of the solution is represented on two-dimensional figures at various time steps. The full time evolution is represented on a three dimensional figure. Logarithmic representations can be far more enlightening if one wishes to reveal the magnitude of the reflected waves. The color legend for all 2D logarithmic representations reported below is given in Figure 9 and the color legend for all 3D logarithmic representations reported below is given in Figure 10.   We first report on some cases where the initial Gaussian condition is transported towards one edge of the rectangle and does not give rise (in the considered time interval) to multiple wave reflections. The first case is depicted in Figure 11 where the chosen velocity is c = (1, 0) (and therefore all boundary conditions (4.18a)-(4.18c) coincide since = 0). We call them DTBC of order 0 since (4.18a) amounts to not retaining any tangential dependence. The results are of course as accurate as in one space dimension (because the scheme (2.2) is one-dimensional in that case and the tangential index only enters as a parameter). We now turn on the tangential dependence by considering the case c = (1, 0.1). Both and are positive then. On the right and left faces of the rectangle, we either implement: -The approximate DTBC (4.18a) and (4.17a). These are called DTBC of order 0; accordingly our choice of approximate DTBC for the top and bottom faces of the rectangle are (4.20a) and (4.19a). The results are shown in Figure 12. The reflected wave at the (right) outflow boundary has magnitude 10 −3 . -The approximate DTBC (4.18b) and (4.17b). These are called DTBC of order 1; accordingly our choice of approximate DTBC for the top and bottom faces of the rectangle are (4.20b) and (4.19b). The results are shown in Figure 13. The reflected wave at the (right) outflow boundary now has magnitude 10 −5 , which shows significant improvement. -The approximate DTBC (4.18c) and (4.17c). These are called DTBC of order 2; there is a subtlety here because our choice of approximate DTBC for the top and bottom faces of the rectangle are still (4.20b) and (4.19b), and not (4.20c) and (4.19c) as one might expect. This is further explained below. The results are shown in Figure 14. The reflected wave at the (right) outflow boundary now has magnitude 10 −8 and is therefore not visible. This shows again significant improvement.
The overall conclusion that can be drawn from these numerical results is that, up to stability issues (which can be challenging), adding more tangential terms in the Taylor expansion (4.5) seems to lead to more accurate results, which was to be expected of course since we consider smooth solutions. The numerical simulations quantify the gain of adding each new term in (4.5b) and (4.5c) with respect to (4.5a).
We now report on more intricate cases that exhibit multiple wave reflections. Namely, we now choose the velocity c = (1, 0.3). With this velocity, the initial Gaussian function should leave the rectangle through the right edge of the boundary. Because the numerical boundary conditions (4.18) are not exactly transparent, wave reflection occurs. For the leap-frog scheme, the reflection takes place according to the billiard law. A (small amplitude highly oscillating) wave then moves towards the top face of the rectangle and a second wave reflection occurs. For this second reflection, the small tangential frequency assumption that motivated the expansion (4.5) is not any longer valid. Hence it is likely that there will be little absorption through the top boundary. The numerical simulations confirm that expectation. We follow the same plan as for the case c = (1, 0.1) and report on some simulations with increasing tangential approximation. We thus implement: -The zero order DTBC (4.18a), (4.17a), (4.20a), and (4.19a) on the four sides of the rectangle. This case is shown in Figure 15; after the first reflection on the right face, the wave has magnitude 10 −3 . The magnitude remains of order 10 −3 after the second reflection on the top boundary. -The first order DTBC (4.18b), (4.17b), (4.20b), and (4.19b) on the four sides of the rectangle. This case is shown in Figure 16; after the first reflection on the right face, the wave has magnitude 10 −5 . The magnitude remains of order 10 −5 after the second reflection on the top boundary. The overall accuracy is improved. -The second order DTBC (4.18c) and (4.17c) on the right and left boundaries with the first order DTBC (4.20b) and (4.19b) on the top and bottom boundaries. This case is shown in Figure 17; after the first reflection on the right face, the wave has magnitude 10 −6 . The magnitude remains of order 10 −6 after the second reflection on the top boundary. The overall accuracy is improved again. In that case, the initial condition is transported towards the upper right corner of the computational domain. There is again an improvement when passing from zero order DTBC to first order DTBC (compare the scales between Figs. 18 and 19). Figure 20 shows that coupling approximate DTBC with different tangential approximation orders on the two sides of a corner may affect the overall accuracy. Namely, for the velocity c = (1, 2/3), using DTBC of order 1 on the four sides of the rectangle yields better results than coupling DTBC of order 2 on two opposite sides with DTBC of order 1 on the two remaining sides.
Eventually, we show in Figure 21 the evolution of the ℓ 2 norm of ( , ) as a function of the time in the case c = (1, 0.3) when the second order DTBC (4.18c), (4.17c), (4.20c), and (4.19c) are implemented on the four sides of the rectangle. The computation exhibits two stages. In the first stage, the initial condition (mostly) exits through the right boundary. As long as the solution remains small on the boundary, the effect of the numerical boundary conditions is rather invisible and the ℓ 2 norm remains approximately constant. In a second stage, the trace of the solution is no longer small and an exponential instability takes place. The logarithmic scale in Figure 21 clearly shows the exponential behavior in time. The numerical solution displays a profile that is concentrated along the top and right boundaries (right of Fig. 21). It thus seems clear that the coupling between (4.18c) and (4.20c) gives rise to an unstable eigenvalue that is associated with a "surface wave" that has exponential decay with respect to the normal directions to both the right and top boundaries. We have unfortunately not been able to obtain an analytical proof of this fact, but our conclusion is that despite good stability properties of (4.18c) and (4.20c) when considered separately on each side of the rectangle, the coupling at the corner may yield strong instabilities.    At last, we have also tested the efficiency of the sum of exponential approximation compared to the standard convolution procedure in the case of DTBC of order 1 (4.18b), (4.17b), (4.20b), and (4.19b). Namely, we have compared the computational effort in that case with the sum of exponential approximation of both sequences ( 0 ) and ( 1 ). We have chosen ( , ) = (50, 20) for the degrees of the Padé approximant. The results are presented in Table 1. Though there is little gain for the sum of exponential approximation in one space dimension, there