Computation of High Frequency Wave Diffraction by a Half Plane via the Liouville Equation and Geometric Theory of Diffraction

We construct a numerical scheme based on the Liouville equation of geo- metric optics coupled with the Geometric Theory of Diffraction (GTD) to simulate the high frequency linear waves diffracted by a half plane. We first introduce a condition, based on the GTD theory, at the vertex of the half plane to account for the diffractions, and then build in this condition as well as the reflection boundary condition into the numerical flux of the geometrical optics Liouville equation. Numerical experiments are used to verify the validity and accuracy of this new Eulerian numerical method which is able to capture the moments of high frequency and diffracted waves without fully resolving the high frequency numerically.


Introduction
In this paper, we construct a numerical scheme for the high frequency wave equation in two-dimension: here c(x) is the local wave speed and ǫ ≪ 1. When the essential frequencies in the wave field are relatively high, and thus the wavelength is short compared to the size of the computational domain, direct simulation of the standard wave equation will be very costly, and approximate models for wave propagation based on geometric optics (GO) are usually used [9,12]. We are concerned with the case when there are some wedges in the computational domain, i.e. the tips and discontinuity in the boundary. When waves hit the wedges, there will be reflections and diffractions.
One of the approximate models for high frequency wave equation is the Liouville equation, which arises in phase space description of geometric optics (GO) [9,32]: (1.4) where the Hamiltonian H possesses the form The derivation of GO does not take into account the effects of geometry of the domain and boundary conditions, which give rise to GO solutions that are discontinuous. Diffractions are lost in the infinite frequency approximation such as the Liouville equation. In this case, correction terms can be derived, as done in Geometric Theory of Diffraction (GTD) by Keller in [25]. GTD provides a systematic technique for adding diffraction effects to the GO approximations. The methods for computing the GO solution can be divided into Lagrangian and Eulerian methods.
Lagrangian methods are based on the ODEs (1.6). The simplest Lagrangian method is the ray tracing method where the ODEs in (1.6) together with ODEs for the amplitude are solved directly with numerical methods for ODEs. This approach is very popular in standard free space GO, [6], and the diffractions, [2,8]. The ray tracing method gives the phase and amplitude of a wave along a ray tube, and interpolation must be applied to obtain those quantities everywhere when rays diverge. Such interpolations can be very complicated for diverging rays.
In the last decade, Eulerian methods based on PDEs have been proposed to avoid some of the drawbacks of the ray tracing method [1]. Eulerian methods discretize the PDEs on fixed computational grids to control errors everywhere and there is no need for interpolation. The simplest Eulerian methods solves the eikonal and transport equations in GO. This technique has been used in standard GO [12]. However, the eikonal and transport equations pick up the so-called viscosity solution [15]. When there are multivalued solutions, more elaborate schemes must be devised. Recently several phase space based level set methods for high frequency waves, in particular the multivalued solutions in GO are based on the Liouville equations, see [5,10,13,18,19,30].
There are very few results on Eulerian methods for diffractions. In this direction, we mention recent numerical methods for creeping waves [28,29,37]. For curved interfaces, the authors [24] constructed Eulerian method for diffraction at interfaces that takes into consideration of partial transmissions, reflections and diffractions. The idea was to revise the transmission/reflection interface condition used by Jin and Wen [20,22] for the Liouville equation in the case of critical and tangent angles to account for diffractions. The diffraction coefficients and decay rates derived in GTD are used in the interface condition. These interface conditions are then built into the numerical fluxes for the Liouville solver. Such an Eulerian computational method is able to capture the moments of high frequency waves without-at least away from the interfaces-numerically resolving the high frequencies, yet still captures the correct interface scatterings and diffractions.
This paper is to revise our previous work [24] to a different geometry, namely, waves through a half plane. When a wave hits a half plane, it usually reflects. However, at the vertex of the half plane, it generates diffracted waves into all directions. In particular, the diffracted waves can reach the shadow zone-the zone that the GO theory cannot cover. We provide a diffraction condition, based on the GTD theory, at the vertex to reflect this diffraction nature. We then build this condition, as well as the reflection boundary condition, into the numerical flux of the Liouville solver, in order to capture the diffractions. This paper is organized as follows. The GO approximations by the Wigner transform for wave equation are sketched in Section 2. In Section 3, we present the behavior of waves at a half plane based on the GTD theory, and provide the conditions for (1.4) that account for reflections at the half plane and diffractions at the vertex of the half plane. In Section 4, the diffraction conditions derived in the previous section is built into the numerical flux in the two space dimension. Numerical examples are given in Section 5 to validate the model and to verify the accuracy of the scheme against the full simulation based on the original wave equation (1.1)-(1.3). Finally, we make some concluding remarks in Section 6.

Geometric optics approximation of the wave equation in the phase space
Consider the two dimensional wave equation We introduce the new dependent variables The energy density is given by , while each of the matrices D i is constant and symmetric with entries either 0 or −1.
To study the GO limit of solution of (2.5), we assume that the coefficients of the matrix A(x) vary on a scale much longer than the scale on which the initial data vary. Let ǫ be the ratio of these two scales. Rescaling space and time coordinates (x,t) by x → ǫx,t → ǫt, one obtains Note that the parameter ǫ does not appear explicitly in (2.6). It enters through the initial data (2.7). We are interested in the initial data of the standard GO form Following [31], one can study the GO limit of (2.6) by using the Wigner distribution matrix W ǫ : n e ik·y w ǫ (t,x−ǫy/2)w ǫ (x+ǫy/2) t dy, (2.9) where n is the space dimension and w t is the conjugate transpose of w. Although W ǫ is not positive definite, it becomes so as ǫ = 0. The energy density for (2.6) is given by As ǫ → 0, the high frequency limit of E ǫ (t,x) is x,k))dk = a + (t,x,k)dk, (2.11) where the amplitude a ± (t,x,k) is given by (2.14) and therefore one needs only to keep track of a + (t,x,k) which satisfies the Liouville equation [31] ∂a + ∂t +c(x)k·∇ x a + −|k|∇ x c(x)·∇ k a + = 0.
(2.15) Therefore, a + can be interpreted as phase space energy density. It solves the Liouville equations (1.4)-(1.5), with the zeroth moment giving the spatial energy density E (0) (t,x) as in (2.11).
The GO approximation is good when ǫ is very small. For moderately small ǫ, diffraction can not be ignored. The Liouville equation (2.15), valid at ǫ=0, does not contain any information about reflection, which occurs even for ǫ=0, nor diffraction which occurs for ǫ > 0. It is not valid near the half plane. In the next section, we will discuss the behavior at an wedge.

Wave behavior on the half plane
In GO, a wave moves with its energy distribution governed by the Liouville equation (1.4). When a wave hits the half plane, it will be completely reflected. However, according to GTD, when the wave hits the vertex of the half plane, it will produce diffracted waves into all directions (see Fig. 1).
In the vicinity of the vertex Q of the half plane, there exist boundary layers [4]. A boundary layer is a narrow zone through which the waves undergo rapid variations. Its thickness is O(ǫ 1/2 ). Diffraction coefficients D o away from the vertex of the half plane (outside the boundary layer) is given [26] by with α is the incident angle, θ is the diffracted angle, and r is the distance from Q. The "+" sign applies when the boundary condition on the half plane is u = 0 (the soft boundary condition), while the "-" sign applies if it is ∂u/∂n = 0 (the hard boundary condition). In GTD, the considered wave propagation phenomena are the incident direct illumination, reflection, and diffraction by vertex of the half plane. The surrounding area of the half plane can be divided into three regions ( Fig. 2): • Region I, where the field is composed of incident, reflected and diffracted waves.
• Region II, where the diffracted waves are added to the incident waves.
• Region III, where only the diffracted waves, generated at the vertex of the half plane exist. No incident and reflected wave reach this region .
Clearly, the above diffraction coefficient is not valid when θ = π +α or θ = 2π −α, i.e., near the shadow boundary (boundary between Region II and Region III) and the reflection boundary (boundary between Region I and Region II). There are boundary layers near the shadow boundary and the reflection boundary with thickness of order √ ǫ.
The Uniform Geometric Theory of Diffraction (UTD) [27] can overcome the singularities at the shadow and reflection boundaries by introducing the transition functions. The uniform diffraction coefficient for UTD is given by with a(β) = 2cos 2 β 2 and the transition function in which one takes the principle (positive) branch of the integral, and β = θ ±α.
The diffraction coefficient in the boundary layer of the vertex of the half plane is [4] We match D U o and D l by finding r 0 , which is the r satisfying D U o =D l . And the matched new diffraction coefficient for the half plane diffraction is with D − for the soft boundary condition and D + for the hard boundary condition. From the above formula, we have D(θ,α) = D(α,θ). The relation between the diffraction coefficient and the relative wavelength ǫ is shown in Fig. 3, for two representative angles. Near the shadow boundaries . We now discuss the wave behavior in more details. Assume the half plane is (x,y) x = x 0 ,y ≤ y 0 , and the vertex of the half plane is at (x 0 ,y 0 ).
Let x=(x,y) and v=(ξ,η). Assume the incident wave hits the half plane with velocity (ξ,η). There are two possibilities: 1. The wave hits the half plane at point (x 0 ,y) with y < y 0 . In this case, the wave will be completely reflected with a new velocity (−ξ,η).
2. The wave hits the vertex (x 0 ,y 0 ) of the half plane. In this case, according to the GTD, the wave can partly diffract and partly travel in the original direction. Introduce the polar coordinates by To describe the correct reflection and diffraction at the half plane, we provide suitable reflection boundary conditions at the half plane and a diffraction condition at the vertex. This has been the strategy for transmission and reflection done in the work of Jin and Wen [20,22] and for diffractions through curved interfaces by the authors [24]. The diffraction condition here is different from the previous works, and is new.
For the wave hitting the half plane below the vertex, it will be completely reflected, At the vertex of the half plane, we use the following condition to account for the diffraction: f + (t,x 0 ,y 0 ,r,α) . These conditions corresponds to the following physical picture. When a particle hits the half plane besides the vertex, it will be completely reflected with a negative momentum. But when a particle hits the vertex of the half plane with an incident angle α, it will be diffracted at angle θ with diffraction coefficient D(θ,α); and the energy of the particle c(x)|v| = c(x) ξ 2 +η 2 = r will not change.
In (3.5), the density function of waves f + (t,x 0 ,y 0 ,ξ,η) is a superposition of the incident wave that passes through the vertex, and all diffracted waves, generated by other incident waves, that move in the direction of v = (ξ,η). These conditions will be used in the next section to construct the numerical flux on the half plane.

The numerical flux
Consider the 2D Liouville equation . Assume the half plane is on the line x = x i 0 +1/2 , (see Fig. 4) and the vertex of the half plane is at the point (x i 0 +1/2 ,y j 0 ).
Let ∆t be the time step, t n = n∆t. The cell average of f is defined as while f n ijkl = f ijkl (t n ). We approximate c(x,y) by a piecewise bilinear function, and for simplicity, we always provide two values of c(x,y) at half plane. Let the half plane value of c(x,y) be c ± i±1/2,j = 1 ∆y y j+1/2 y j−1/2 c(x ± i±1/2 ,y)dy, where the superscripts ± indicate the right and left limits of the quantity at the half plane. Note c + The 2D Liouville equation (4.1) can be semi-discretized as 2 ,jkl and all the numerical fluxes are defined using the upwind discretization, except for f ± i 0 + 1 2 ,jkl (for j ≤ j 0 ). Here we treat the half plane as an interface, and then, similar to [20,24], define the numerical flux here using conditions (3.4) and (3.5).
By (3.5), we can define the numerical flux at (x i 0 + 1 2 ,y j 0 ) in an upwind fashion. Firstly, we divide the interval [− 1 2 π, 3 2 π] into 2I subinterval [θ m ,θ m+1 ],θ m = m∆θ,∆θ = π/I,m = 0,1,··· ,2I, then approximate the diffraction condition (3.5) by with ξ k =rcosα,η l =rsinα. The first two terms are the diffraction from other waves hitting the vertex with an incident angle θ m ,m = 0,1,··· ,2I. When m ≤ I, − 1 2 π ≤ θ m ≤ 1 2 π, the diffraction term comes from the waves hitting the vertex from the left to right, and when m ≥ I, 1 2 π ≤ θ m ≤ 3 2 π, the diffraction comes from the waves hitting the vertex from the right to left. Since ξ ′ m = rcosθ m ,η ′ m = rsinθ m may not be grid points, we have to define them approximately. One can first locate the cell centers that bound these velocities, and then use a bilinear interpolation to evaluate the value at (ξ ′ m ,η ′ m ). Next we consider the numerical flux at half plane below the vertex point ( where the superscripts ± indicate the right and left side of the quantity at the half plane.
Similarly, if ξ k < 0, (4.7) After the spatial discretization is specified, one can use any time discretization for the time derivative.

Numerical examples
In this section we present numerical examples to demonstrate the validity of our scheme and to show its numerical accuracy. In the numerical computations the second order Runge-Kutta time discretization is used.
Since it is difficult to get the exact solution for this problem, as in [24], we use the numerical solution with the mesh size small enough to represent the exact solution. The two-dimensional Lax-Wendroff method with space meshsize ∆x =∆y= ǫ 20 and ∆t= ∆x/2 are used to solve the system (1.1) in the form with s = ∂u/∂t,r = ∇u to get the energy density distribution The numerical energy density is defined as and r ij can be defined similarly.
The discrete wave equation is quite dispersive [7], so one needs many grid points per wavelength to compute it. The mesh size h = ǫ/20 is the biggest mesh size we can get satisfactory numerical results for the discrete wave equation.
The limit energy density is the zeroth moment of the density distribution of Liouville equation The computational tool we used is the super computer in Tsinghua National Laboratory for Information Science and Technology, 512 Itanium 2 64 bit processor, the peak computational speed is of 2.662×10 13 , the total EMS memory is 1024G, the storage space is 26T.
In the computation, we first approximate the delta function initial data of the Liouville equation by the product of a discrete delta function in 1-D [11]: with ω = ∆ξ = ∆η (For more recent numerical studies on the approximations of the delta function, see [33][34][35][36]). Then the energy density distribution are recovered by We use the L 1 -error in the cumulative distribution function (cdf), i.e., the antiderivative of energy density [14] +∞ (5.5) which can be expected to flatten as ǫ is decreased, to measure the weak convergence in the semiclassical limit. Lemma 2.1 in [3] ensures that (5.5) going to zero is equivalent to the weak convergence of E (0) (x,y,t) For a more through discussion about the model error and numerical discretization error of this approach we refer to our previous work [24].
with ǫ = 1/3000 and some suitable boundary conditions (to be specified later) on Γ = (x,y)| x = 0.2,y ≤ 0.1 . The approximate Liouville equation is with initial data The choice of initial data is based on (2.12)-(2.15). So 1.28(x+y) 2 comes from ∇u, and 2 comes from ∂u/∂t. For convenience, we denote our scheme by GTD and the pure geometrical optics method by GO. Fig. 5 shows the numerical energy densities E (0) of GTD and GO compared with E at t = 0.2,0.3. One can see that there are some diffracted waves behind the half plane-the shadow zone. The numerical results of GTD is very close to the solution of the wave equation, which shows that our method can capture the main feature of the diffracted waves. But the GO can't capture diffraction in the shadow region.  mesh type 50 2 ×50 2 100 2 ×100 2 200 2 ×200 2 400 2 ×400 2 t=0.1 3.0032e-2 1.1623e-2 5.7521e-3 2.8282e-3 t=0. 2 3.3943e-2 1.3945e-2 7.0934e-3 3.5082e-3 t=0. 4 4.2894e-2 1.6402e-2 8.1006e-3 4.0097e-3 Table 1 presents the relative errors of the numerical energy density E (0) computed with different meshes in the phase space at t = 0.1,0.2 and 0.4. The l 1 error is very small. The convergence rate is about first order. Table 2 shows the relative errors of the numerical energy density E (0) in the shadow  zone (x ≥ 0.2,y ≤ 0.1). The GTD solution is a good approximation to the solution of the wave equation in the shadow zones. Notice that the convergence rate in the shadow region is smaller than first order. This is partly because that there is a boundary layer near the shadow boundary, which is harder to resolve numerically then elsewhere. Next, we consider the problem with the hard boundary condition on Γ, i.e. ∂u ∂n Γ = 0. We use the extrapolation boundary condition for the Lax-Wendroff method in the fully resolved simulation of the high frequency wave equation. The physically relevant values for the diffraction coefficient D + (θ,α) are given by (3.3). Fig. 6 shows the numerical energy densities E (0) and E at t = 0.2,0.3. One can see that the energy of the diffracted waves behind the half plane-the shadow zone is stronger than the case of the soft boundary conditions. The numerical results of GTD is very close to the solution of the wave equation and the numerical results of GO are not right in the shadow region. Table 3 presents relative errors of the numerical energy density E (0) computed with different meshes in the phase space at t = 0.1,0.2 and 0.4. The convergence rate is of first order. mesh type 50 2 ×50 2 100 2 ×100 2 200 2 ×200 2 400 2 ×400 2 t=0.1 2.4953e-2 1.1299e-2 5.5003e-3 2.6706e-3 t=0. 2 3.1014e-2 1.2596e-2 6.1806e-3 3.0848e-3 t=0. 4 3.6141e-2 1.5024e-2 7.1363e-3 3.5453e-3 Table 4 shows the relative errors of the numerical energy density E (0) in the shadow zone. The GTD solution is a good approximation to the solution of the wave equation in the shadow zones.
The solution of GTD and GO depend on wavelength ǫ, Fig. 7 gives when t = 0.2, the relation between the relative errors of GTD and GO and the wavelength. One can see that the error of solution of GO and GTD is of same order-near O(ǫ). The reason that both GTD and GO are first order in ǫ is that the diffraction effect is only important near the vertex, and along RB and SB which are lower dimensional things. So even if GTD is more accurate in these places, the L 1 error is based on the error of the entire domain which is basically O(ǫ) for both GTD and GO.
The corresponding Liouville equation is Firstly, we simulate the problem with the soft boundary condition. The physically relevant values for the diffraction coefficient D − (θ,α) is given by (3.3). Fig. 8 shows the numerical energy densities E (0) and E at t = 0.3,0.4. The numerical results of GTD is very close to the solution of the wave equation, which shows that our method can capture the main feature of the diffracted waves. On the contrary, diffractions are lost in the GO solutions. Table 5 presents the relative errors of the numerical energy density E (0) computed with different meshes in the phase space at t = 0.15,0.3 and 0.4. The error is very small. The convergence rate is about first order.  mesh type 50 2 ×50 2 100 2 ×100 2 200 2 ×200 2 400 2 ×400 2 t=0.15 2.0044e-2 9.2131e-3 4.6428e-3 2.3042e-3 t=0. 3 2.5404e-2 1.1043e-2 5.4954e-3 2.7751e-3 t=0. 4 3.1988e-2 1.2408e-2 6.1198e-3 3.0087e-3 Table 6 shows the relative errors of the numerical energy density E (0) in the shadow zone. The GTD solution is a good approximation to the solution wave equation in the shadow zones (x > 0.11,|y| < 0.2).
Finally, we consider the problem with the hard boundary condition on Ω. The physically relevant values for the diffraction coefficient D + (θ,α) are given by (3.3).   Fig. 9 shows the numerical energy densities E (0) and E at t =0.3,0.4. GTD can capture the diffractions, but GO can not capture the diffractions. Table 7 presents the relative errors of the numerical energy density E (0) computed with different meshes in the phase space at t = 0.15,0.3 and 0.4. The convergence rate is about first order. Table 8 shows the relative errors of the numerical energy density E (0) in the shadow   It is important also to point out that only near the vertex we need to impose ∆x,∆y ∼ O( √ ǫ). Away from it we can use ∆x,∆y,∆ξ,∆η = o(1) if we program the method in the adaptive mesh framework. Even though the time step is still required to be of O( √ ǫ), this is still a tremendous saving compared with the full wave simulation.

Conclusion
In this paper, we revise our previous work [24] to a different geometry, namely, high frequency waves through a half plane. When a wave hits a half plane, it usually reflects.
However, at the vertex of the half plane, it generates diffracted waves into all directions. In particular, the diffracted waves can reach the shadow zone-the zone that the GO theory cannot cover. We provide a diffraction condition, based on the GTD theory, at the vertex to reflect this diffraction nature. We then build this condition, as well as the reflection boundary condition, into the numerical flux of the Liouville solver, in order to capture the diffractions. This gives an Eulerian computational method for high frequency waves through a half plane, which is able to capture wave reflection and diffractions at a half plane without fully resolving the high frequency waves in the entire computational domain.
An immediate application of this method is to the problem of high frequency waves diffracted by a wedge. Similar ideas, including those in our previous work [24], can also be applied to other geometries, and to elastic and electromagnetic waves, which will be the subjects of future research.