Absorbing boundary conditions for the Westervelt equation

The focus of this work is on the construction of a family of nonlinear absorbing boundary conditions for the Westervelt equation in one and two space dimensions. The principal ingredient used in the design of such conditions is pseudo-differential calculus. This approach enables to develop high order boundary conditions in a consistent way which are typically more accurate than their low order analogs. Under the hypothesis of small initial data, we establish local well-posedness for the Westervelt equation with the absorbing boundary conditions. The performed numerical experiments illustrate the efficiency of the proposed boundary conditions for different regimes of wave propagation.


1.
Introduction. Constantly growing needs of numerical simulations in science and engineering often require considering problems which are naturally formulated in unbounded domains. Typical examples can be found in many problems originating from fluid dynamics, solid mechanics, aerodynamics, electrodynamics, acoustics, etc. However, the numerical solution of such problems requires a finite region. There are basically two approaches which can be used to reformulate problems in infinite domains as problems in finite domains. The first approach is to map the originally unbounded domain to a bounded one. Simple as the problem sounds the solution in practical applications is far from known. This is mostly due to reasons which are connected with singularities of the new equation that results from the mapping. The second approach, which we follow in this work, is to impose fictitious boundaries to truncate the domain of interest. Such artificial boundaries require special boundary conditions so that the boundary value problem is well-posed and its solution is an accurate approximation to the restriction of the solution in the unbounded domain. In other words, these boundary conditions have to be transparent to or, as they are usually called, absorbing for solutions propagating outwards the artificial boundary. It is commonly recognized that absorbing boundary conditions (ABCs) play a key role in computations on unbounded domains and have a profound impact on the accuracy of numerical methods. Over the past thirty years, ABCs have developed into a vigorous research direction including a wide spectrum of methods and approaches. The description of these techniques is out of the scope of this work and therefore we restrict ourselves to referring the reader to the comprehensive review articles [10,34,14,15,11,12] and the references therein.
The focus of this work is on construction of ABCs for high-intensity focused ultrasound (HIFU) which plays an important role in many medical and industrial applications such as diagnostic ultrasound [8,31,30], thermotherapy of tumors [9,16,4], lithotripsy [1], ultrasound cleaning and sonochemistry. Linear models of wave propagation are not applicable in HIFU due to nonlinear effects requiring more sophisticated acoustic equations to be taken into account. In this work, we develop local in space and time ABCs for the Westervelt equation used as a basic acoustic model in various HIFU simulations. The Westervelt equation is one of the fundamental equations governing the propagation of acoustic waves in nonlinear regimes [35,16,4,3]: where Ω ⊆ R d , d ∈ {1, 2, 3}, u = u(·, t) is the acoustic pressure, c > 0 is the speed of sound, b > 0 is the acoustic diffusivity, > 0 is the mass density, β a = 1 + B/(2A) with B/(2A) > 0 standing for the parameter of nonlinearity of the fluid, T is the final time at which the problem is to be solved. All the parameters are assumed to be constant. We rewrite (1) in a form more convenient for further treatment c −2 u tt − ∆u − β∆u t = γ(u 2 ) tt in (0, T ) × Ω (2) with β = b/c 2 , γ = β a /( c 4 ), and complement (2) by initial conditions u(·, t = 0) = u 0 , u t (·, t = 0) = u 1 in Ω, and by inhomogeneous Neumann and absorbing boundary conditions where ∂Ω = Γ N ∪ Γ A , subscript n denotes the normal derivative on the boundary, and the operator A, on the absorbing boundary Γ A , is an annihilating operator for outgoing waves which we specify in due course.
In spite of the intensive research activity in the field of transparent boundary conditions, most results have been obtained for linear problems with constant coefficients. Wave equations with variable coefficients have received much less attention, not to mention nonlinear models. There are only few papers devoted to problems with variable coefficients [7], convective [2] and nonlinear [17,32,38,29] terms. Despite the existence of some approaches to the construction of ABCs for nonlinear wave models their application to concrete equations is rather sophisticated and still out of the scope of most research works.
In this work we design ABCs based on the theory of pseudo-differential [25,18,28] calculus. We will also address a possible approach via para-differential [21,27] calculus in the appendix. The first approach is applicable to linear wave equations with variable coefficients. Therefore it is used for the Westervelt equation linearized in a neighborhood of a reference solution. The second approach will be directly applied to the nonlinear Westervelt equation. Before going into detail with the derivation of ABCs, we remark that both theories have been already used in the construction of transparent boundary conditions. For example, the pseudo-differential calculus was exploited by Engquist and Majda in [7] to design ABCs for the linear wave equation with variable coefficients. Transparent boundary conditions for the semilinear wave equation as well as for the nonlinear Schrödinger equation with the help of para-differential operators were obtained in [33] and [32], respectively. Remark 1. The nonlinearity in the Westervelt equation comes along with a strong damping term b∆u t . In fact, this strong damping, besides being a physically imposed term, also plays a quite particular mathematical role. As already observed in [23], [24] in the context of different boundary conditions, strong damping β > 0 is essential in two and higher space dimensions in order to compensate the nonlinearity and avoid degeneracy in the equation, (whereas in 1-d we have wellposedness also in case β = 0). On the other hand, the strong damping term destroys the wave like character of the equation since it implies decay of the energy and a rather parabolic than hyperbolic behaviour of the equation, cf. [23]. This results in the observation that the (linearized) differential operator defining the Westervelt equation with strong damping is not amenable to a factorization (see, e.g., (19) below) as required for constructing absorbing boundary conditions. For this reason, we skip the strong damping term during derivation of the ABCs. Of course it has to stay in the PDE, though (for physical reasons and since otherwise wellposedness would fail, as mentioned above). It is clear that the inevitable use of integration by parts in deriving energy estimates for the PDE causes the appearance of boundary terms resulting from the presence of the term b∆u t . So these terms finally have to be taken into account in the ABC as well. As mentioned already, the factorization approach based on pseudo-or paradifferential calculus is not appropriate for doing so. Thus, incorporation of the β-term in the boundary conditions will be done as a postprocessing step after the pseudodifferential factorization, and it will be done on the basis of energy considerations. The latter will also allow to prove well-posedness of the resulting initial boundary value problems for the Westervelt equation.
1.1. Main results. The novelty of our work lies in the derivation and analysis of high-order ABCs for the Westervelt equation which have not been construct so far. We will do so for the one-and two dimensional versions of the Westervelt equation (2) first of all in a domain without corners, see Section 2. Additionally, we will provide well-posedness results for the Westervelt equation with zero and first order conditions in one and two space dimensions in Section 3.
In this section we summarize the boundary conditions derived in this paper together with the main well-posedness results.
For the case of one space dimension we will derive zero, first and second order ABC (with µ defined as in (40)) in Section 2.1. In 2-d, the zero order ABCs derived in Section 2.2 are and the first order ones are where subscript ϑ denotes the tangential derivative. Here the zero order ABC and the first line of the first order ABC are exactly what one would expect from the linear case with constant coefficients. Note that for reasons outlined in Remark 1 above we set β = 0 in these derivations. The energy considerations in Section 3.1 will allow us to appropriately take into account the third order derivative term going with β. With the according modifications in the ABCs (5), (6), (8), (9), and denoting u β := u + βu t we will obtain the following local in time well-posedness results for sufficiently small initial data u 0 , u 1 , and Theorem 1.1. For β ≥ 0, any open interval Ω = (a, b) ⊆ R and any T > 0 there exists ρ > 0 such that for all initial data u 0 , u 1 satisfying u 2 L 2 (Ω) + u 1 H 1 (Ω) < ρ, a solution u ∈ C 2 (0, T ; L 2 (Ω)) ∩ C 1 (0, T ; H 1 (Ω)) to exists and is unique.
exists and is unique.
For β > 0, any smooth and bounded domain Ω ⊆ R 2 and any T > 0 there exists ρ > 0 such that for all initial data u 0 , u 1 satisfying u 2 L 2 (Ω) + exists and is unique.
For β > 0, any smooth and bounded domain Ω ⊆ R 2 and any T > 0 there exists ρ > 0 such that for all initial data u 0 , u 1 satisfying u 2 L 2 (Ω) + exists and is unique.
The remainder of this paper is organized as follows. In Section 2 we derive absorbing boundary conditions for the Westervelt equation via (formal) pseudodifferential calculus in one and two space dimensions. Section 3 is devoted to energy estimates and the proofs of Theorems 1.1-1.4. In Section 4 we provide numerical results.
2. Derivation of absorbing boundary conditions for the Westervelt equation in one and two space dimensions. In our derivation, without loss of generality we consider the simple domains Ω = (−∞, 0] in 1-d and Ω = (−∞, 0] × R in 2-d, where x plays the role of the outward unit normal and (in 2-d) y is the tangential direction. Moreover, we will skip the term β∆u t for the reasons outlined in Remark 1.

2.1.
Absorbing boundary conditions in 1-d via linearization and pseudodifferential calculus. As it was already mentioned, the direct reformulation of (2) in terms of pseudo-differential operators is not possible because of the nonlinear term on the right hand side. Therefore we consider some linearization around a reference solution u (0) of this equation. After derivation of the ABCs from this inhomogeneous linear wave equation with variable coefficients, we will re-insert u (0) = u to arrive at ABCs for the Westervelt equation. The reason for using (14) (as was also done for the wellposedness proof in [23]) and not the standard linearization according to first order Taylor expansion, which would be (15) is that the offset terms −2γu (0) u tt would lead to problems with the commutativity of the pseudodifferential operators below.
For simplicity of exposition we first of all consider the one-dimensional version of the Westervelt equation (2) In 1-d the linearization (14) reads as where we set ν 2 = ν 2 (u (0) ) with and point out that the analysis of the Westervelt equation is based on estimates that actually make sure positivity of c −2 − 2γu, so that ν 2 > 0 is a natural assumption. In order to derive transparent boundary conditions for the linearized Westervelt equation (17) we make use of the theory of pseudo-differential calculus. For the purpose of this formal derivation, ν is assumed to be a C ∞ function both in time and space, as needed for applying pseododofferential calculus. Since we do not prove this smoothess, our derivations are only formal.
The key idea behind the derivation of ABCs is mostly based on the Nirenberg factorization of (17) written in terms of pseudo-differential operators. To construct approximate boundary conditions one can factorize the operator D 1 as where A = A(x, t, D t ) and B = B(x, t, D t ) are pseudo-differential operators with symbols a(x, t, τ ) and b(x, t, τ ) from the space The differential operator D t is defined as −i∂ t with the imaginary unit i, and R is a smoothing pseudo-differential operator with the Schwartz kernel k(x, y) ∈ C ∞ satisfying [19] ( At the symbolic level, factorization (20) reduces to with the correspondence iτ ↔ ∂ t between the frequency and the (physical) time domains, and where by a slight abuse of notation, for a function f , we denote the symbol of the zero order differential operators u → f u (multiplication operator) again by f . Now, we have to define symbols a and b in (21). For doing so, it is worth to remark that formally these symbols admit the following asymptotic expansions where a 1−j (x, t, τ ) and b 1−j (x, t, τ ) are homogeneous of degree 1 − j in τ . To proceed, one has to substitute (22) in (21) and equate symbols of the same degree of homogeneity on both sides of equality (21). However, before this substitution, we recall the reader the definition of the product of two pseudo-differential operators which are used owing to the term ab in (21).
In accordance to the theorem on the product of two pseudo-differential operators [37], has the asymptotic expansion of its symbol c(x, ξ) ∈ S m1+m2 given by for every nonnegative integer N and with the standard multi-index notation α = (α 1 , α 2 , . . . , α k ) and |α| = α Thus, the symbol c := ab of the product of the pseudo-differential operators Substitution of (22) and (24) in (21) and casting-out R lead to Evidently, the more coefficients are taken in (25) the more accurate ABCs are. However, taking more coefficients also makes the ABCs more complicated and involved to implement since they will contain higher order derivatives. Therefore, we only show how to find {a j , b j } j={1,0,−1} . In order to define the first pair of coefficients a 1 , b 1 one has to equate the symbols with the degree of homogeneity O(τ 2 ). This gives the following system of equations The solutions to (26) are given by We take to make the terms of order O(τ 2 ) vanish.
Remark 2. The choice of the sign in front of ν(iτ ) is not arbitrary. This sign defines the propagation direction of the wave.
In order to find the next pair of coefficients a 0 , b 0 we equate symbols with degree of homogeneity O(τ 1 ). In other words, we have to solve the system with the operator A 0 := ∂ x + ν∂ t .
In order to obtain more accurate boundary conditions one has to equate the symbols with degree of homogeneity O(τ 0 ) which leads to the following system The solution of (32) is given by Taking into account (28) and (31) we deduce that Note that with the Taylor linearization (15) an offset term γ(u (0) 2 ) tt would have appeared here which would have prevented the equality a −1 a 1 = a 1 a −1 . (Here, we write f for the symbol of the zero order differential operator u → f (constant mapping), which has to be strictly distinguished from the multiplication operator u → f u.) This problem is avoided by using the fixed point type linearization (14). In accordance to [26], the operator annihilates outgoing waves at {x = 0} × (0, T ). Substitution of the asymptotic expansion (22a) with the first k leading terms results in the following boundary condition  i.e., an ABC of order k is obtained by keeping the first k terms in the asymptotic expansions (22). Thus in order to construct a zero order ABC we set k = 0 and substitute the coefficient a 1 in (36) which gives Parallel to the construction of the zero order ABCs (2.1), we set k = 1 and substitute a 1 , a 0 in (36) to obtain the first order boundary conditions: For k = 2 we obtain the second order ABCs where we have multiplied with (iτ ) before converting from symbols to operators, Inserting u itself for the a priori solution u (0) , we arrive at zero first (42) and second order nonlinear ABCs. We will see in Section 5.1 that slightly different conditions result from derivation via a paradifferential approach.

2.2.
Absorbing boundary conditions in 2-d via linearization and pseudodifferential calculus. In the spatially two dimensional situation on the domain (−∞, 0) × R, where ν is defined by (18), we proceed very similarly to the 1-d case: We consider pseudo-differential operators A = A(x, y, t, D y , D t ) and B = B(x, y, t, D y , D t ) with respect to time and tangential (i.e., y) direction, but the expansion is still with respect to powers of τ , so equations (19), (20) (with A = A(x, y, t, D y , D t ) and B = B(x, y, t, D y , D t )) remain the same whereas (21), (22), (25) change to with the correspondence iη ↔ ∂ y and and respectively, where a 1−j and b 1−j are homogeneous of degree 1 − j in τ (and are additionally functions of x, y, t, and η). As in [6], in our derivations we will rely on an assumption of the type η ∼ τ or even η τ small. Considering the O(τ 2 ) terms in (47) leads us to in place of (26), which leads to At this point, a fundamental difference to the 1-d case arises, since we will have to approximate the square root in order to derive practically applicable boundary conditions. We will do so by a Taylor expansion whose order is adapted to the order of the ABCs.
The computations for a 0 , b 0 look exactly the same as in the 1-d case and yield i.e., To obtain zero order boundary conditions we use the zero order Taylor expansion (52) For our first order boundary conditions we use the first order Taylor approximations for the terms that are nonlinear with respect to τ , η in (50), (51). This yields the symbols Again we insert u itself for the a priori solution u (0) to arrive at zero order ABCs and at first order ABCs where we have multiplied the symbols with (iτ ) to obtain (54).
3. Well-posedness. In this section we will show well-posedness of the Westervelt equation with zero or first order ABC derived above in one or two space dimensions. Note that zero order ABC have already been considered in [3]. However, the conditions there do not take into account the nonlinearity in the highest order time derivative. Moreover, the proof in [3] is carried out in higher space dimensions, which necessitates the use of higher order energies. In 1-d this is not required (simply due to the fact that in 1-d already H 1 embeds into L ∞ ) and the proof is on one hand much simpler, on the other hand it enables existence also of spatially less smooth solutions and well-posedness in the absence of interior damping (i.e., with β = 0). For these reason we will also provide the well-posedness proof for the 1-d Westervelt equation with zero order ABC (41) here.
Since we derive energy estimates by only multiplying with u t for zero order ABC in 1-d, the strong damping term and the terms resulting from its integration by parts at the boundary will be easily tractable in that case. However, for the first order ABCs, carrying out energy estimates following the idea in [13], the β term yields derivatives of u on the boundary that are too high to be controllable by other boundary (or, via trace theorems, interior) terms. Therefore we will modify the first order ABC accordingly to account for the strong damping and arrive at decaying energies. Note that in the derivations of section 2 we had omitted the β terms since they would have destroyed commutativity. The terms that we insert now again in favor of energy decay are different from those omitted in section 2, though. I.e., the (formal) Nirenberg factorization from there would not have helped in obtaining energy dissipation. In fact it turns out that the ABCs derived in section 2 (plus the β-modifications made here) only allow us to show local in time well-posedness. As to be expected, the resulting ABCs coincide with the classical Engquist-Majda ones in case of constant coefficients and vanishing damping.
3.1. Energy identities for the strongly damped inhomogeneous wave equation with variable coefficients. Before proceeding to well-posedness of the nonlinear Westervelt equation with zero and first order ABCs, we will derive some energy identities (especially we will carry over the energy identities used in the well-posedness proof for first order ABCs in [13]) for inhomogeneous wave equations with variable coefficients and strong damping of the form x(, y)), β ≥ 0 and initial conditions u(t = 0) = u 0 , u t (t = 0) = u 1 . This will provide us with crucial information on how to incorporate the strong interior damping term into the ABCs and help to prove well-posedness of the Westervelt equation with ABCs in the next subsections. For simplicity of exposition we will here restrict ourselves to a geometry with x being the boundary normal and y the boundary tangential direction, respectively. The general case can be covered by applying smooth local boundary transformations.
Multiplying the PDE with u t we obtain where u β = u + βu t . This suggests to use as zero order absorbing boundary conditions √ αu t + u β n = 0 (or √ αu t + u β n = lower order terms), where "lower order terms" are expressions whose L 2 (0, T ; L 2 (∂Ω)) inner product with u t can be dominated by the energy and/or the interior dissipation and/or the boundary dissipation Similarly, if we differentiate the PDE wrt t and multiply with u tt , we arrive (after space and time integration) at the energy identity where Multiplication of the time differentiated PDE with αu tt (instead of u tt ) yields Considering the PDE that results from (55) for u β , and multiplying with u β tx , we get the energy identity For the combined higher order energy functional the identities (60), (63) yield where we have used the fact that x is the outward normal direction in our setting. This suggests to use first order boundary conditions leading to the identity αu tt + u β xx + βu β txx + 2 √ αu β tn = 0 (or lower order terms) where this time "lower order terms" are expressions whose L 2 (0, T ; L 2 (∂Ω)) inner product with u β tn can be dominated by the higher order energy E 2 [u](t) and/or the interior dissipation Using the PDE (61) for transforming second order normal (i.e., x) derivatives to tangential (i.e., y) derivatives, we can achieve the boundary identity (65) e.g. by the first order ABCs where in one space dimension the yy derivative terms are just skipped.
In the nonlinear 2d case we need to establish an L ∞ ((0, T ) × Ω) estimate of u within the coefficient c −2 − 2γu in order to guarantee nondegeneracy. We will do so via the embedding H 2 ((0, T ) × Ω) → L ∞ ((0, T ) × Ω), the Poincaré inequality applied to the domain (0, T ) × Ω with fixed Cauchy data on the boundary part {0} × Ω, as well as the energy estimate resulting from multiplication of the PDE with −∆u: 3.2. Zero and first order ABCs in 1-d; proof of Theorems 1.1, 1.2. We prove well-posedness and boundedness of the enrgy E 1 [u] as in (59) of the following initial boundary value problems with Ω = (−1, 1).
To show that T is a contraction on W, we use the fact that for v 1 , v 2 ∈ W, and u i = T v i , i = 1, 2, the functionû = u 1 − u 2 solves the following problem the estimate where we have used E 0 [û](0) = 0 as well as an estimate similar to (74) together with u 2 t (±1) L 2 (0,T ) ≤c to estimate c −1 (±1) L 2 (0,T ) . Since α 1/4 ≥ 4 c −2 − 2γm, and by estimate (72) forû t ,v t in place of u tt , we arrive at an estimate of the form with a constant C(m,m,ā,b,c) that can be made small byā,b sufficiently small. Hence, we achieve contractivity of T on W with respect to the norm induced by Thus, using Banach's Contraction Principle, we have shown Theorems 1.1, 1.2.
Remark 3. It is readily checked that replacing c −2 − 2γu by a constant c −1 , we would get rid of a couple of higher derivative terms on the boundary and end up with energy estimates enabling even global in time wellposedness even with β = 0, cf. [22]. Thus, for obtaining enhanced approximation by taking into account the full time and space dependence of the coefficient c −2 − 2γu in (68) we pay the price of losing global in time well-posedness and needing strong damping.

Remark 4.
For β > 0 it is possible to make use of maximal parabolic regularity to prove more general existence results in L p spaces even under less restrictive assumptions on the regularity of the inital data in case of pure Dirichlet boundary conditions cf. [36]. However, it seems to be at least not really straightforward to carry over these techniques to absorbing boundary conditions.

Zero and first order
in case of zero order ABC and the combined energy in case of first order ABC, with an appropriately chosen factor λ > 0. Moreover, to include pointwise bounds −m,m on u into the definition of W for avoiding degeneracy, we make use of the fact that in both cases obviously an estimate of the form with C possibly depending on T holds. Accordingly, we will use combination of the energy identities (58), (67) and (60), (63), (67), respectively for showing the self mapping property in case of zero and first order ABCs, respectively. Contractivity in case of zero order ABC will rely on the lower order energy identity (56), applied to the initial boundary value problem that holds for the difference between two solutions of the linearized problem. In case of first order ABCs we will have to use a higher order energy identity also for contractivity in order to take into account the tangential derivative terms. This is the only part of the proof that we will provide explicitely here, since the rest (self-mapping for zero and first order ABCs, contraction for zero order ABCs) goes very much along the lines of the proofs in [3], [23]: To show that the operator T mapping v ∈ W to a solution u of we use the fact that for v 1 , v 2 ∈ W, and u i = T v i , i = 1, 2, the functionû = u 1 − u 2 solves the following problem wherẽ Hence using (64) with α = α v 1 = c −2 − 2γv 1 , f = 2γv 1 t , g = 2γv t u 2 t + 2γvu 2 tt we obtain for the energy and the interior dissipation where we have used E 2 [û](0) = 0. Let us first consider the terms within the curly braces {. . .} under the integral over Ω and (0, t) on the right hand side of (79). After bounding the L ∞ norm of the α factors by α = c −2 + 2γm we see that they are all of the form where ∂ is a combination of differential operators (id + β∂ t ), ∂ x , ∂ t , ∇. Hence the time and space integrals of these terms can be estimated by products of the form where f 1 , f 2 , f 3 is an appropriate permutation of q, φ, ψ .
It is readily checked that since v 1 , u 2 ∈ W as defined in (78), all q factors can be bounded by constants depending on m,m,ā,b,c that can be made small for smallā,b,c, and all φ and ψ factors can be bounded by either the energy norms . For the terms within the brackets [. . .] under the integral over ∂Ω and (0, t) on the right hand side of (79) we have that their squared L 2 (0, t; L 2 (∂Ω)) norm is bounded by some constant depending on m,m,ā,b,c (which can be made small for small a,b,c), multiplied with Altogether we arrive at an estimate of the form which for smallā,b,c gives the desired contractivity estimate.
4. Numerical results. In this section we study the performance of the proposed boundary conditions and compare them with the first and second order Engquist-Majda ABCs [7] for different setups. In what follows, we focus on a horizontal waveguide in one and two dimensions, and study how the accuracy of ABCs is influenced by the angle of incidence in the 2-d case. Then, we consider the highintensity focused ultrasound (HIFU) problem with the physical parameters typical for simulations of thermotherapy for human liver cancer and analyze how intensively the solution is contaminated by the reflected waves. We name the ABCs as ABC d,o n , where the superscripts d and o indicate the space dimension and the order of ABC, while the subscript n takes the value PS or EM standing for the new nonlinear ABC obtained with the pseudo-differential calculus or the Engquist-Majda ABC, respectively. To approximate system (2)-(4) in time the standard Newmark scheme is applied [20]. For space discretization, the finite element method is used.
In order to compare different ABCs, a reference solution u * is computed in the domain Ω Ω which is large enough to prevent the solution in the restricted domain Ω from being polluted by reflected waves. The studied ABCs are compared in terms of the l 2 -norm relative error δ = u * − u 2 / u * 2 . In all numerical experiments the number of finite elements per wavelength is set to be 50, and the time step is chosen in such a way as to have 20 time samples per time period for each of the frequencies ω = {100 kHz, 1 MHz}. To induce a wave in the domain, we use a monofrequency transducer of the form u n = sin(2πωt). The time t as well as the acoustic pressure are normalized to their maximum values. The physical parameters used in all the computations correspond to those of human liver [16,4]: c = 1596 m · s −1 , ρ = 1050 kg · m −3 , B/A = 6.8, b = 2αc 3 /(2πω) 2 , with the acoustic absorption coefficient α = 4.5 Np · m −1 · MHz −1 .
4.1. ABC in 1-d. In this part we compare the zero and first order ABCs on a line segment x ∈ [0, 3 cm]. The transducer with a 100 kHz excitation frequency is set at x = 0 while the ABC is prescribed at x = 3 cm. The results of the comparison are presented in figure 1.
As it can be seen from figure 1, the behaviour of the boundary conditions brings no surprise: the higher the order of ABC the more accurate solution we have. In contrast to ABC 1,0 PS and ABC 1,1 PS , ABC 1,1 EM is of much less accuracy. This result is expectable and reconfirms the attention one has to pay to the ABCs for the Westervelt equation. 4.2. ABC in 2-d. The zero order versions of the proposed ABCs, as well as the Engquist-Majda ABCs, provide the best absorption of the wave hitting the boundary at normal incidence. The higher order the ABC are, the better a deviation from this specific angle should be taken into account. Therefore, it is worth studying how the ABCs react to different angles θ of incidence. In this respect, we consider 100 kHz waves traveling from left to right in a rectangular waveguide Ω = [0, 3 cm] × [0, 0.05 cm] on the right wall of which one of the studied ABCs is set. We begin from showing the well-known effect: the more the deviation of the incidence angle from zero the more the solution is contaminated by reflected waves. The first example presented in this series is a wave impinging the boundary at θ = 0 • . In this simple case all the ABCs should work the best. The results are shown in figure 2. As in the 1-d case, ABC 1,1 PS outperforms ABC 1,0 PS . However, ABC 2,2 EM , which coincides with ABC 1,1 EM for θ = 0 • , shows low accuracy. In the next setup, we increase the angle θ by 15 • , which is a assumed to be a harder trial for the ABCs (see figure 3). Indeed, all the ABCs revealed to be quite sensitive to the angle of incidence and exhibit higher errors introduced by the reflected waves into the solution. The second order Engquist-Majda ABC gives more accurate results compared to the first order condition, but the error is quite large. The new boundary condition of the first order demonstrates the lowest error while the zero order ABC is less efficient. In the last example, the ABCs are studied in a much more realistic situationthe HIFU problem which is routinely used in computational setups to simulate the thermotherapy for human liver cancer. We consider a concave transducer, with a much higher excitation frequency ω  At the very beginning of the simulation (t < 0.3), the first and second order Engquist-Majda ABCs work equally well. However, the situation gets worse as time advances: the discrepancy between the boundary conditions grows and the solution becomes substantially contaminated by the reflected waves. The second order Engquist-Majda ABC does not dramatically affect the situation, and the numerical solution is still quite poor. The proposed ABCs demonstrate much lower errors. However, the difference between them is less pronounced compared to the waveguide example. Another remarkable feature of the new ABC is that the error exhibit a much less fluctuating behavior, which suggests that the new boundary conditions are robust with respect to the wave propagation regime.

5.
Conclusions. In this work we proposed zero and first order ABCs, based on pseudo-differential calculus, for the Westervelt equation in one and two space dimensions. Well-posedness of the boundary value problem with the new ABCs is stated and proven. All our numerical results reconfirm the fact that using the ABCs which are not especially tailored for the Westervelt equation lead to poor numerical solutions. The zero order ABCs are computationally easier than the first order conditions, however, more prone to the regimes of the wave propagation and less accurate. It is important to remark that the application of the self-adapting technique [?] to the developed ABCs will result in further improvements. defined as where F is the Fourier transformation, χ ∈ C ∞ (R d × {R d \ {0}}) is a function of homogeneity degree zero satisfying χ(ζ, η) = 1 if |ζ| ≤ ε 1 |η|, where 0 < ε 1 < ε 2 . Let us consider a nonlinear differential equation of order N defined by the superposition operator induced by Φ with Φ ∈ C ∞ and x ∈ R d . In accordance to [21], the para-linearization of (82) with Φ(x, ·) vanishing at 0 is given by T ∂Φ ∂λα (·,u,...,∂ α u,...) 0≤|α|≤N ∂ α u + R(u), where T a is a para-differential operator with symbol a, and R(u) is a smooth error. More precisely, for all u ∈ H s (R d ) with s > d/2 equation (83) implies R(u) ∈ H 2s−d/2 (see [27]). Equation (83) is often referred to as the para-linearization formula of Bony. Before the derivation of ABCs for the Westervelt equation (16), we develop the nonlinear term on its right hand side γ(u 2 ) tt = 2γ((u t ) 2 + uu tt ) and recast (16) in the form ν 2 (u)u tt − u xx − βu txx = 2γ(u t ) 2 with ν 2 (u) = c −2 − 2γu. Based on (83) and taking into account that (see [21]) where T f and T g are para-differential operators with symbols f and g, we get a para-differential equation instead of the nonlinear Westervelt equation (16). Acting similar to the previous derivation, we can apply Nirenberg's factorization, analogous to (19), and rewrite (86) in the form where A and B are para-differential operators with symbols a and b, respectively. A similar argument as for the linearized Westervelt equation yields Note that this equation differs from (21) and it will also lead to different ABCs. Again we will skip the β terms for the same reason as in section 2.1 The equation for the O(τ 2 ) terms remains the same, namely (26) Thus we get The equation for the O(τ 1 ) terms is different from (29), namely a 0 + b 0 = 0, which upon setting β = 0 yields (differently from (31) with the operator A 0 := ∂ x + ν∂ t . Finally the O(τ 0 ) equation is (differently from (32) so that we get Parallel to the ABCs for the linearized Westervelt equation from section 2.1, we obtain the zero, first and second order boundary conditions where B 2 :=μ(u), which contains multiplication with u tt , as opposed to (39).
As in the pseudo-differential case, we do not consider higher order boundary conditions, although their derivation follows the same lines.
Remark 5. It is probably due to the approximation by Taylor linearization (83) for an actually quadratic nonlinearity that the paradifferential approach yields a different and in our view most likely worse ABC than the pseudodifferential one from section 2, since the linearization (14) seems to be better capable to capture the quadratic nonlinearity than just Taylor linearization (cf. (15)). E-mail address: i.shevchenko@imperial.ac.uk E-mail address: Barbara.Kaltenbacher@uni-klu.ac.at