GWsim: a code to simulate gravitational waves propagating in a potential well

Let $\Omega\subset\mathbb{R}^{3}$ be a bounded domain with boundary $\partial\Omega$. The boundary value problem (BVP) we aim to study in this work can be presented as

which is a second-order linear hyperbolic partial differential equation (PDE) with initial conditions

Let $\partial\Omega_{1}$ and $\partial\Omega_{2}$ are subsets of the boundary $\partial\Omega$ with $\partial\Omega_{1}\cup\partial\Omega_{2}=\partial\Omega$ and $\partial\Omega_{1}\cap\partial\Omega_{2}=\emptyset$. We assume that the boundary conditions on $\partial\Omega_{1}$ are Dirichlet and on $\partial\Omega_{2}$ is Neumann

where $\hat{n}$ is the outward pointing unit vector normal to $\partial\Omega_{2}$. $q(x,t)$ is a given function on the boundary $\partial\Omega_{2}$.

In mathematics, Eq. (1), together with the initial and boundary conditions Eqs. (2,3,4), are called the strong formulation. The strong formulation is the most commonly used format of the wave equation. It can be numerically solved using some straightforward methods such as the finite difference method (FDM). However, the strong formulation is not the only format of the wave equation. In what follows, we introduce an alternative format, namely, the weak formulation of the above equations.

Multiply Eq. (1) by a test function $\phi$ and then integrate over $\Omega$, we obtain a variational equation

where we have used the Green’s formula

The second term on the right hand represents the boundary term, which can be further split into terms according to different boundary conditions

The test function $\phi$ is chosen in such a way that it vanishes on the subset of boundaries with Dirichlet boundary conditions ${\partial\Omega_{1}}$, which is equivalent to say that we restrict $\phi$ to lie in

where $H^{1}(\Omega)=W^{1,2}(\Omega)$ is called the first order Sobolev space, which is also a Hilbert space, meaning that $\phi$ and its first order weak derivatives $\partial_{x}\phi$ are square integrable.

Given the restrictions of the test function $\phi$, the boundary term in Eq. (5) can be reduced to

Here, in the second equality, we have imposed an absorbing boundary condition

The physical meaning of the absorbing boundary condition will be explained later on.

The above we have discussed the boundary condition on $\partial\Omega_{2}$. On the other hand, the boundary condition on $\partial\Omega_{1}$,

is called the essential boundary condition, which does not affect the variational equation Eq. (5) but this condition must be imposed directly on the solution itself.

Equipped with the boundary conditions Eqs. (9,10), Eq. (5) is called the weak formulation of Eq. (1). It can be shown that if Eq. (5) holds for all $\phi\in V$, it has a unique solution $u\in V$ Lax & Milgram (1954). This solution is called the weak solution. It is obvious that the classical solutions of Eq. (1) are also weak solutions. However, conversely, weak solutions must have sufficient smoothness to be the classical solutions, which dependents both on the shape of the boundary $\partial\Omega$ and the regularity of $c^{2}f$. For example, it can be shown that if $\Omega\subset\mathbb{R}^{n}$ has $C^{k}$ boundary (smooth up to $k$-th derivative) and $c^{2}f\in H^{k}(\Omega)$ with $k>\frac{n}{2}$, then $u\in C^{2}(\bar{\Omega})$, namely $u$ is a solution in the classical sense Grisvard (1985).

In addition to the weak formulation of the wave equation, another important property of the wave function Eq. (1) is that it obeys the law of energy conservation. To see this point, we multiply Eq. (1) by $\frac{\partial u}{\partial t}$ and integrate over the domain $\Omega$

Note that the total energy carried by a wave in the domain $\Omega$ is defined by

The left hand of the equality of Eq. (11) then represents the change of the total energy of wave with time. From Eq. (11), if the boundary condition on the right hand side of the equality is a constant ($\frac{\partial{u}}{\partial t}|_{\partial\Omega}=0$) and there is no force source $f=0$ in the second term, for a constant speed of wave $\nabla c^{2}=0$, the total energy is conserved $\frac{dE(t)}{dt}=0$.

Equation (11) is an important feature of the wave function. We shall discuss this point in our following numerical analyses.

We use superscript $n$ to indicate the number of a time step and $k=t_{n}-t_{n-1}$ is the length of the present time step. We discretize the time derivative as

The discretization scheme we adopted here also involves the spatial quantities at two different time steps, which is known as the $\theta$-scheme. Equations (14) then can be represented as

When $\theta=0$, the scheme is called the forward or explicit Euler method. If $\theta=1$, it reduces to the backward or implicit Euler method. The Euler method, no matter implicit or explicit, is only of first-order accurate. Therefore we do not adopt here.

The scheme we use here is called the Crank-Nicolson Scheme, namely $\theta=1/2$, which uses the midpoint between two different time steps. This scheme is of second-order accuracy. As the midpoint rule is symplectic, the most important feature of this scheme is that it is energy-preserving. We shall discuss this point in our following numerical analyses.

We discretize the domain $\Omega\subset\mathbb{R}^{3}$ using the Finite Elements Methods (FEM). The FEM is an effective method for numerically solving partial differential equations. There are many excellent textbooks on this topic (Readers are referred to e.g. Ref. Zienkiewicz & Taylor (1989) for more details). Here, we will not go into too many mathematical details. Instead, we only summarise the main steps here.

First, the domain $\Omega$ is decomposed into subdomains $\Omega_{i}$, which consist of rectangles and triangles. This is called decomposition or triangulation. The vertices of rectangles and triangles in the domain $\Omega$ are called mesh points or nodes. Let $\Omega_{h}$ denote the set of all nodes of the decomposition. On each node, we construct a test function $\phi_{i}\in V\,,i=1,..,N$, where $V$ is the space defined in Eq. (7) and $N$ is the total number of nodes in the domain. The test function $\phi_{i}$ is required to have the property

where

Thus, $\phi_{i}$ has non-zero values only on the node with $k=i$ and its adjacent subdomains. It vanishes on other parts of the domain $\Omega$. The test function $\phi_{i}$ constructed this way is called the shape function. Clearly, $\phi_{i}\in V$ on different nodes are linearly independent. We denote the space spanned by $\phi_{i}$ as $V_{h}:=\mathrm{span}\{\phi_{i}\}_{i=1}^{N}$, which is a subspace of $V$.

Next, we expand the scalar fields $u,v,c^{2}f$ in terms of $\phi_{i}$

where $u_{h}\,,v_{h}\,,c^{2}f_{h}\in V_{h}$. These functions are approximations of the original scalar fields. The advantage of this expansion is that the coefficients are precisely the values of these functions on the node points (e.g. $u_{i}=u_{h}(p_{i})$).

Inserting the above expansions into Eqs. (16) and noting that Eqs. (16) holds for any $\phi\in V$, thus, we can choose $\phi$ as $\phi=\phi_{j}$, where $j=1,..,N$. Then we obtain $N$ different equations, which form a linear system

where $G_{1}$ and $G_{2}$ are defined by

and the elements of the matrixes are defined by

It is worth noting that the wave speed $c$ is not a constant so that the matrix $D$ is not symmetric $D_{ij}\neq D_{ji}$, which needs to be kept consistent with Eqs. (17). Equations (18,19) are called the Galerkin equations, from which the unknown coefficients $u_{i}$ and $v_{i}$ can be solved. Note that Eq. (18) explicitly contains $V^{n}$, which is unknown at the time step $n$. However, we can multiply Eq. (18) by $k\theta$ and then use Eq. (19) to eliminate $V^{n}$. Then we obtain

for $U^{n}$ and

for $V^{n}$.

Now, $U^{n}$ can be solved using the information on the time step $n-1$. Then $V^{n}$ can be solved with the knowledge of $U^{n}$. In practice, Eqs. (20,21) usually contain a very large number of equations. The rank of Eqs. (20,21) (the number of linear independent equations) is called the degrees-of-freedom (DOF) of the system, which depends both on the number of node points and the freedoms in the shape functions.

As the systems of linear Eqs. (20,21) may contain a large number of unknowns(e.g. which can be easily up to a few $10^{7}$), direct methods such as the LU decomposition is inefficient in this case. One has to use iterative methods such as the Conjugate Gradient (CG) method, which is efficient for large but sparse matrices. However, our cases are complicated due to the varying wave speed $c$. The matrix $D$ is no longer symmetric. As the CG method can be only applied to symmetric, positive-definite matrices, we can not use the CG method in this work. Instead, we adopt the GMRES (a generalized minimal residual algorithm for solving non-symmetric linear systems) method Young & Chen (1996), which does not require any specific properties of the matrices.

In this work, we use the public available code deal.ii Arndt et al. (2019); Bangerth et al. (2007); Alzetta et al. (2018) to numerically solve Eqs. (20,21). deal.ii is a C++ program library aimed at numerically solving partial differential equations using modern finite element method. Coupled to stand-alone linear algebra libraries, such as PETSc Abhyankar et al. (2018); Balay et al. (2019a, b, 1997), deal.ii supports massively parallel computing of vary large linear systems of equations. deal.ii also provides convenient tools for triangulation of various geometries of the simulation domain. Therefore, it is straightforward to implement Eqs. (20,21) in the deal.ii code. The detailed numerical analyses are presented in the next section.

The boundary-initial value problem for the propagation of waves from a point source can be presented as

where $Q(t)$ is the time-dependent waveform of the point source and $\delta$ is the Dirac delta function. The wave equation has a fundamental solution, which is the distributional solution of

The distribution $G$ is also called the Green’s function. In free-space, the explicit form of $G$ for an outgoing wave is given by

The above equation is known as the retarded Green’s function in free-space. The general solution of Eq. (22) then has the form

where $r$ is the norm of $x$, namely, $r=\left\|x\right\|$.

Although the wave function from a point source has a very simple analytic expression Eq. (25), numerically solving Eq. (22) is, indeed, non-trivial. A challenge is that the point source on the right-hand side of Eq. (22) is singular, which may cause numerical problems. In practice, in order to avoid this problem, we can make use of the Huygens-Fresnel principle to numerically solve Eq. (22) on a slightly different domain without the origin but with new boundaries

The Huygens-Fresnel principle states that the original waves from a source propagating to a distant observer can be considered as re-radiating from a wavefront of the original wave rather than directly from the source. As such, we can use the wavefront of the original wave at some radius as a new boundary for the wave equation. In this way, we can avoid the singularity in Eq. (22) at the origin. A mathematical description of the Huygens-Fresnel principle for spherical waves is provided in Appendix A.

Figure 1 shows the new boundaries and triangulation of our simulation domain. The domain consists of the bulk regions between the two spherical shells with an inner radius $R_{i}$ and an outer radius $R_{o}$, respectively. We set the boundary conditions at the inner shell $R_{i}$ as

which is exactly the wavefront of the original point source propagating there. At the outer shell, we adopt the absorbing boundary condition

The absorbing boundary condition is also called the non-reflecting boundary conditions or radiating boundary conditions. These boundary conditions can absorb and eliminate the reflections of waves on the boundaries. Therefore, with these boundary conditions, we can approximate the propagation of waves in free-space using a limited volume of the simulation domain. In addition to the boundary conditions, we set the initial conditions as

After fixing the boundary and initial conditions, we perform several numerical tests. In these tests, we choose a concrete waveform $Q(t)=-A\sin(\omega t)$. Thus, the spherical wave has an analytical expression

where $\omega=2\pi f\,$ and $f$ is the frequency of waves. $A$ is the amplitude of the wave. In our test, we set $R_{o}=10^{7}M_{\odot}=49.2535{\rm Sec}\,,$$R_{i}=9\times 10^{6}M_{\odot}=44.32815{\rm Sec}$ and $f=1$. Further, we set $A=4\pi R_{i}$, namely, the amplitude of the waves are normalized to unity at the inner boundary.

To perform our simulations, we further need to set the time step. For a time-dependent wave equation, an explicit discretization method is only stable if the time step is small enough so that the wave can have enough time to travel through the space discretization (see. Chapter 2 in Ref. Grossmann et al. (2007)). This condition is known as the Courant-Friedrichs-Levi condition

where $\sigma$ is the size of the mesh and $k$ is the size of the time step as mentioned before. However, for the implicit methods such as the one used in this work, it can be stable for arbitrary step sizes if the scheme is upwind (see. Chapter 2 in Ref. Grossmann et al. (2007)). However, a stable scheme is not sufficient to guarantee a correct solution to the wave equation. This is because of the oscillating features of waves. If the frequency of waves is too high, the wavelength will be too small and it is difficult to resolve the waveform within one period. The shape of the waves, in this case, is not smooth relative to the size of the meshes, which can lead to large errors in the numerical simulations. Therefore, in order to get a correct solution, the time step and the size of the meshes should be small enough so that the waveform within one period can be well resolved. As such, in our simulations, we set $k=\lambda/30$ and the size of the mesh $\sigma<\lambda/7$. We have tested that further reducing the size of the mesh $\sigma$ does not improve the results appreciably.

We perform a suite of simulations with different spatial resolutions. In these simulations, we adopt the linear Lagrange test function, which depends only on the values of the finite elements on their vertices but does not involve any derivatives of the unknown fields. Therefore the degree of freedom (DOF) associated with each vertex is one. In this case, the total DOF is simply the total number of the expansion coefficients of the unknown field ($u_{i}$ or $v_{i}$). This number is also equal to the total number of linear algebraic equations in the system. Table 1 lists the level of refinement and the corresponding DOF.

Figure 2 shows the waveform observed by an observer located at $r=46.63{\rm Sec}$. Dashed lines represent the numerical results and the solid lines are obtained from Eq. (29). The dashed vertical line indicates the epoch that the wavefront just arrives at the observer. Dashed lines with different colors show the results with different resolutions. With high spatial resolution $2^{7}$ (blue dashed line), the numerical results agree with the analytic one very well. This is expected as the accuracy of the FEM is strongly dependent on the resolution of elements.

After investigating the time-domain waveform for a particular local observer, we consider the integral property of the wave equation. Inserting Eq. (29) into Eq. (12), we obtain the total energy of waves in the domain

Then we measure the total energy in our simulation numerically.

Figure 3 compares the numerical results (stars) with the analytical ones (solid line) from Eq. (31). When the wavefront has not reached the outer boundary, namely, $R_{i}<t<R_{o}$, the numerical results agree with the analytic ones very well. After $t>R_{o}$ (vertical dashed line), the wavefront starts getting out of the simulation domain. The total energy inside the bulk of our simulation domain becomes steady (flat) as, in this case, the energy getting into the simulation domain equals the energy getting out of the domain.

We assume a metric $g_{\mu\nu}$ given by

and the gravitational waves as a perturbation to the background metric

where $g_{\mu\nu}^{0}$ is the background metric. Following Ref. Peters (1974), we neglect higher order non-linear terms and arrive at the propagation equation for gravitational waves $h_{ij}$

Using the eikonal approximation by Ref. (Baraldo et al. (1999)) (also see Ref. Weinberg (1972)), the gravitational wave tensor can be represented as

where $e_{ij}$ is the polarization tensor. In this work, we assume that the polarization tensor does not change during the propagation of GWs. Then we obtain the so-called scalar wave equation

where $c^{2}=1/(1-4\psi)$ is the speed of wave.

Unlike the conventional wave equation, one notable feature of this equation is that the wave speed is not a constant but is dependent on the potential well. The varying wave speed plays an important role in the propagation of GWs in a potential well. Moreover, given the fact that $\psi<0$, we have $c<1$. The presence of potential thus delays the propagation of GWs. This phenomenon is known as the Shapiro time-delay Shapiro (1964).

Figure 4 shows the domain and setups of our simulations, where $R_{i}$ is the distance from the source to the surface of the simulation domain. We assume that the source of GWs is far away from the simulation domain. The incident wave train travels along the $z$-axis and enters the surfaces of the simulation domain at $z=0$. Due to the long-range nature of gravity produced by the scatterer, the incident waves suffer the Shapiro time delay on their way from the source to the scatterer, which adds up a shift of phase to the wave function at the boundary of the simulation domain. Since the source is distant from the scatterer and the length of the simulation domain is small relative to $R_{i}$, the wave function at the simulation domain can be well approximated by plane waves with a constant phase shift $\Delta t=-2M\ln(2R_{i}/d)$, where $d$ is the distance of the scatterer to the boundary of the simulation domain. Moreover, we choose the size of the simulation domain large enough so that the scatterer does not distort the wave front initially at the boundary of the simulation domain. Such distortion effect can be accurately and robustly determined using the null geodesic equations, which are provided in appendix B. After these considerations, we set the boundary condition at $z=0$ as

where $Q(t)$ is the waveform of the source and the dot denotes the derivative of time.

In this subsection, we choose $Q$ as a sine function. For convenience, the initial phase of the source is so chosen that it cancels the change of phase due to the Shapiro time-delay when waves arrive at the simulation domain. As such, the boundary condition at $z=0$ can be represented as

where $\omega=2\pi f$ and $f$ is the frequency of the wave. The amplitude of the incident wave is normalized as unity at $z=0$. Moreover, we choose the zero point of the temporal axis at the epoch of the wavefront arriving at the surface of the simulation box. For other surfaces of the simulation domain, we adopt the absorbing boundary conditions as explained in Eq. (27). The absorbing boundary conditions aim to guarantee that GWs $h$ do not reflect into the simulation box when they arrive at the boundaries. In addition to the boundary conditions, we set the initial condition inside the simulation domain as

which means that there are no waves in the simulation domain initially.

For $\psi$, we choose it as the potential generated by a solid spherical ball with uniform density

where $R_{s}=2M$ is the Schwarzschild radius and $x_{0}$ is the center of the potential well. Note that, $\psi$ is continuous up to the first order of derivative, namely, $\psi\in C^{1}(\Omega)\subset H^{1}(\Omega)$. We choose the mass of the scatterer as

In this work, we first simulate the behavior of the leading wavefront for a finite wave train in a potential well, we choose the input wave train as a sine function but only lasting for one period

where $\lambda$ is the period. We choose the simulation domain as a cube. We choose the refinement of the simulation domain along each dimension in 3D space as $2^{9}$. The size of the simulation box is taken as $20\,[\rm sec]$ along one side. We place the scatterer at the center of the box. As shown in appendix B, in this case the scatterer is distant from the boundary surface of the incident waves. As such, the distortion effect of the scatterer at boundaries is negligible.

Figure 5 shows the waveform of the finite wave train in a potential well. We choose $\lambda=4[\rm sec]$ for illustrative purposes. From left to right, the panels show waveforms at different times. The snapshots are taken along the $z-y$ plane at $x=0$. The color bar to the right shows the amplitude of waves $h$. The black dot indicates the position of the scatterer. As the wave speed is dependent on the potential well $c^{2}=1/(1-4\psi)$, waves travel faster far away from the center than those closer to the center. The leading wavefront is bent toward the center (left panel). After the scatterer, the wavefront is warped, which leads to a complex 3D shape of the wave zones even when the waves get out of the potential well of the scatterer (right panel).

Unlike the rest part of waves, the leading wavefront represents the transport of energy from one place to another. Therefore, it is strictly subject to the constraint of causality. Moreover, the hypersurface of the wavefront (denoted as $S$) is a null hypersurface. Its normal vector $k^{a}=\nabla^{a}S$ is a null vector $k^{a}k_{a}=0$. As shown in Wald (2010) (see, Eq.(4.2.37)), the integral curves of $k^{a}$ are null geodesics, which satisfy $k^{a}\nabla_{a}k^{b}=0$ if $k^{a}$ is presented in terms of an affine parameter. As such, the behavior of the leading wavefront can be traced by null geodesics, which provide an independent way to test our numerical results. As shown in the appendix, the equations for null geodesics can be presented as (see appendix B for details)

where $\mu=1/r$ and $r$ is the distance to the center of the mass. $\varphi$ is the azimuthal angle in the spherical coordinates. The definition of $\zeta(M,\mu)$ is provided in the appendix B. As the incident wave is a plane wave at $z=0$, the spacial components of its normal vector $k^{a}$ are parallel to the $z$-axis. Therefore, we set the initial conditions for Eq. (42) at $z=0$ as

where $b$ is the impact factor. Eq. (43) simply describes a straight line in spherical coordinates along the $z$-axis, which is also the solution of Eq. (42) in the absence of the scatterer $\zeta(M,\mu)=0$.

Figure 6 shows the trajectories of the null geodesics of the normal vectors $k^{a}$ of the wavefront obtained using Eq. (42) (the bundle of thin solid black lines). The wavefront predicted by null geodesics at $t=17.40[{\rm sec}]$ (solid red lines) agrees very well with those of simulations (shaded regions). Using null geodesics, it is also clear that after the scatterer, the wavefront get warped and the signals of the leading wavefront on the $z$-axis does not come from the signals traveling on the $z$-axis itself, but from the wavefront circling around it.

In addition to the tests using null geodesics, we also test our numerical results using the analytic expressions of the Shapiro-time delay. To to this, we introduce a new notion

where $c$ is the speed of wave $c=\sqrt{1/(1-4\psi)}$. For a point mass, the Shapiro time delay $\Delta t_{\rm Shapiro}$ is given by

where $b$ is the impact factor and the integration is taken along the $z$ direction under the thin-lens assumption. ${z_{i}}$ and ${z_{f}}$ are the initial and final positions of the integration limits. In this work, we take $z_{i}=-10[{\rm sec}]$ and ${z_{f}}=7.40[{\rm sec}]$ which is the $z$ component of the position of the wavefront relative to the scatterer without delay at time $t=17.40[{\rm sec}]$. Fig. 6 shows that in regimes that are far away from the scatterer (where the thin-lens assumption is valid), our numerical results agree with the predictions of the Shapiro-time delay (blue solid line) very well.

Moreover, we demonstrate the robustness of our numerical results. To do this, we compare the results obtained from the runs with refinements of $2^{8}$ and $2^{9}$ for wave functions at different times along the $z$-axis. From Fig. 7, the waveforms for simulations with $2^{8}$ (dashed lines) agree with $2^{9}$ refinements (solid lines) very well, which means that even with $2^{8}$ refinement, it already yields good accuracy. However, a high-resolution simulation, such as $2^{9}$ refinement, can suppress the wavefront shocks, which are numerical artifacts appearing in the wave functions after the scatterer.

After testing the leading wavefront, we further test the waveform of our simulations against analytical predictions. In this new test, we take the incident wave function as a continuous sine function. Moreover, since the incident waves are of azimuthal symmetry, we choose the simulation domain as a cylinder instead of a cube to minimize the effects of boundaries. The incident waves are assumed to travel along the axis of the cylinder($z$ axis). The length of the cylinder is taken as $70[\rm sec]$ and the radius is taken as $25[\rm sec]$. To make it easier to be compared with analytical predictions, we simulate the waves with the evolution time long enough that the waves in the simulation domain can reach a steady state. In this case, the waveform can be analyzed in the Fourier space.

where $\omega=2\pi f$ is the angular frequency. Equation (36) then becomes

Compared with Eq. (36), the potential $\psi$ appears as a source term in Eq. (47). As such, the wave speed in Eq. (47) is unity and the wave number $\tilde{k}$ equals to the angular frequency $\tilde{k}=\omega$, which is different from that in Eq. (36).

Outside the Schwarzschild radius of the scatterer, the potential $\psi$ of the scatterer can be described by $-M/r$. In this case, Eq. (47) can be solved analytically(e.g.  Peters (1974)). However, unlike Eq. (36) which is hyperbolic, Eq. (47) is elliptic, which are specific to the “boundary-value" problems Grisvard (1985). The solutions of Eq. (47) are strongly dependent on the boundary conditions imposed. Since the wave functions are of azimuthal symmetry, we impose the following ansatz as boundaries

where $\xi=r-z$ and $r=\left\|x-x_{0}\right\|=\sqrt{x^{2}+y^{2}+z^{2}}$(the scatterer is placed at the origin). Since $f(\xi)$ does not depend on the direction of $\tilde{k}$, $\tilde{h}$ only has Fourier modes with wave vectors $\tilde{k}$ along the $z$-axis. As such, the above ansatz implicitly assumes that the wave vector $\tilde{k}$ of the outgoing waves at infinity is still along the $z$-axis. Under such assumption, the solution of Eq. (47) is given by

where ${{}_{1}F_{1}}$ is the confluent hypergeometric function and $C$ is a complex constant (detailed discussions are provided in Appendix C). When $M=0$

Wave function $\tilde{h}$ simply goes back to the plane wave. Moreover, it is worth noting that $\tilde{h}$ can not describe waveform on the $z$-axis because ${{}_{1}F_{1}}$ is not well defined when $\xi=0$.

If $r\gg R_{s}$ and $\xi\gg 1$, the confluent hypergeometric ${{}_{1}F_{1}}(i\omega M,1,i\omega\xi)$ can be expanded around $i\omega\xi=\infty$ as

In this case, the wave function can be written as

Note that $\tilde{h}(x,y,z)$ are complex values. The spatial waveform can then be obtained by taking the real part of $\tilde{h}(x,y,z)$. The complex constant $C$ is independent of $(x,y,z)$, which can be determined by matching the amplitude and phase of $\tilde{h}(x,y,z)$ with the incident waves. An advantage of the above equation is that it can be evaluated efficiently in the numerical processes. We therefore use it to test our simulation results.

Figure 8 shows a snapshot of our simulation along the $z-y$ plane at $x=0$. The incident wave travels along the $z$-axis with a waveform as a continuous sine function. The frequency of the incident wave is $f=0.25\,[{\rm Hz}]$. The snapshot is taken at $t=136.50[\rm Sec]$, for which the time evolution is long enough so that the waveform is in a steady state. We have checked that further increasing the evolution time does not change the waveform appreciably. The black dot indicates the position of the scatterer. The signals after the scatterer show a clear diffraction pattern with strong signals along the $z$-axis.

Figure 9 compares the wavefroms of our simulations (black lines) with the analytical predictions (red lines). The wavefroms are taken along 1D lines that are parallel to the $z$-axis. Since the wave function and simulation domain are of azimuthal symmetry, we only extract these lines from the $z-y$ plane. The distances of these lines to the $z$-axis are chosen as $y=10.0563[\rm Sec]$, $y=14.0230[\rm Sec]$ and $y=18.0946[\rm Sec]$, respectively. The waveform of our simulations agree with the analytical predictions very well.