A non-reflecting wave equation through directional wave-field suppression and its finite difference implementation

The acoustic wave equation describes wave propagation directly from basic physical laws, even in heterogeneous acoustic media. When numerically simulating waves with the wave equation, contrasts in the medium parameters automatically generate all scattering effects. For some applications - such as propagation analysis or certain wave-equation based imaging techniques - it is desirable to suppress these reflections, as we are only interested in the transmitted wave-field. To achieve this, a modification to the constitutive relations is proposed, yielding an extra term that suppresses waves with reference to a preferred direction. The scale-factor α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} of this extra term can either be interpreted as a penetration depth or as a typical decay time. This modified theory is implemented using a staggered-grid, time-domain finite difference scheme, where the acoustic Poynting-vector is used to estimate the local propagation direction of the wave-field. The method was successfully used to suppress reflections in media with bone tissue (medical ultrasound) and geophysical subsurface structures, while introducing only minor perturbations to the transmitted wave-field and a small increase in computation time.

www.nature.com/scientificreports/ In this paper, we will build from the latter method towards a robust directional wave-field suppression theory. Using this theory and its implementation, waves are guided through the medium without reflections and without significant perturbations to the transmitted wave-field.

Theory
Inside a lossless, heterogeneous, isotropic, acoustic medium the particle velocity and the acoustic pressure obey the equations of motion and the constitutive equation (stress-strain relation). Together, they can be written as a system of first-order hyperbolic partial differential equations (PDE's): where v = v(� r, t) and p = p(� r, t) denote the particle velocity and the acoustic pressure respectively, and ρ = ρ(� r) and c = c(� r) correspond to the density and compressional wave speed at each position r in the medium. In this paper, a modification of Eq. (2) is proposed to suppress wave-fields along a preferred direction: where α = α(� r) ( m −1 ) determines the strength of suppression, and Ŝ =Ŝ(� r, t) indicates the estimated propagation direction. In order to motivate this modification, we first revert to the second order hyperbolic PDE by taking the temporal derivative of (4) and substituting Eq. (3): Thus, we have obtained the acoustic wave equation with an extra α-weighted, one-way term towards the Ŝ direction. This term corresponds to the term proposed by Fletcher et al. 5 , and has its origin in sponge-like boundary conditions 16 . This particular form was selected because it does not require any auxiliary fields.
To observe the effect of this extra term on the wave-field, we derive the solution of Eq. (5) for a plane wave traveling in the k -direction: where k and ω denote the wave-vector and angular frequency of the plane-wave, respectively. A complete derivation of this result can be found in the supplementary information. In Eq. (6), we see that the proposed modification leads to a directional suppression effect. When Ŝ coincides with the plane wave's direction, e.g: Ŝ ·k = 1 , the plane wave propagates unaltered. On the other hand, when Ŝ and k are opposite, e.g: Ŝ ·k = −1 , the plane wave is maximally suppressed by a factor exp[−αct] . In between these two extremes, suppression is proportional to the cosine of the angle between Ŝ and k . Since this method solely suppresses reflections, and does not affect the transmitted amplitude, energy is not preserved. This method must thus be viewed as a non-physical acoustic wave equation.
The observations above motivate us to define a penetration depth, δ p = α −1 , for reflecting waves propagating in a direction exactly opposite to the incident direction. Alternatively, suppression can be viewed as a temporal process by defining a time-decay constant τ , such that α(� r) = (c(� r)τ ) −1 .

Finite difference implementation
The proposed method for reflection suppression is demonstrated using a O(t 2 , x 4 ) staggered-grid FDTD implementation 17 of Eqs. (3) and (4) with Perfectly Matched Layer (PML) boundary conditions 18 . The additional term in Eq. (4) is implemented using spatial (cubic) and temporal (quadratic) interpolation. The complete FD scheme can be found in the supplementary information.
The choice of FDTD provides the added benefit of allowing one to work with the acoustic Poynting vector, which can be used to determine the local wave-field propagation direction Ŝ (� r, t) . Conveniently, the additional cost for computing and storing this quantity is low. In order to account for regions where the acoustic Poynting vector is ill-defined, we use the time-stacking technique described by Yoon et al. 7 .
Using the above method, the maximum value of α is defined based on the source frequency through a time constant τ . Additionally, the value of α is set to 0 within a small circle around the source location, since the acoustic Poynting vector is not well-defined at the time of source-injection. In theory, τ can be kept constant throughout the medium. In practice, to minimize perturbations to the transmitted wave-field, it is recommended to scale τ with respect to the local medium velocity contrast, e.g: τ (� r) = max(| � ∇c|) | � ∇c| τ min , where τ min denotes the fastest time-decay present in the medium.

Results
First, we examine a medium with a sharp velocity contrast at geophysical scale (Fig. 1a). A point-source Ricker 19 wavelet with a peak frequency of 50 Hz is injected inside the medium, after which the wave-fields arising from Eqs. (1) and (2) are compared with their reflection-suppressed counterparts (3,4), with a value of τ min = 4.68 · 10 −3 s such that α max = 0.14 m −1 (Fig. 1b). In Fig. 1c,d,e we observe the effectiveness of the reflection suppression method via a snapshot and time series display. Reflections at both small and large incident angles are fully suppressed, while the refracted wave-field propagates with little to no perturbations.
The value of the suppression constant used in Fig. 1 was determined by a sensitivity assessment for the value of τ min . Figure 2 shows the level of suppression for different levels of τ min . Complete suppression of the reflections seen in Fig. 1a is reached at a value of τ min = 1 4f 0 , or α max = 1 4cf 0 . The level of suppression can be tuned depending on the level of velocity contrasts in the medium or in the case when there exist specific regions where reflections must be suppressed.
Next, we apply the same methodology in the ultrasound regime to a human skull model 21 (Fig. 3a) with a peak frequency of 200 kHz. We use a value of τ min = 7.80 · 10 −7 s such that α max = 855 m −1 (Fig. 3b).
In Fig. 3c,d,e we once again observe that the reflected wave-fields are very strongly suppressed, while the transmitted wave-fields only exhibit small perturbations with respect to the unmodified acoustic wave equation. www.nature.com/scientificreports/ Lastly, we repeat the same procedure for the geophysical Marmousi 22 model (Fig. 4a) at a peak frequency of 10 Hz. We use a value of τ min = 7.80 · 10 −3 s such that α max = 7.50 · 10 −2 m −1 (Fig. 4b).
We once again observe a complete suppression of reflected wave-fields. The transmitted wave-field of the modified acoustic wave equation in Fig. 4e contains varying perturbations with respect to the acoustic wave equation. In part, these perturbations can be explained due to the wave-front in Fig. 4c containing both reflected and transmitted wave-fields at heterogeneous locations, where Fig. 4d only contains the transmitted wave-field, making it difficult to compare the two figures.

Discussion
The examples show that the method presented in this paper strongly suppresses internal reflections in a robust manner, even within highly heterogeneous media. In addition, the transmitted wave-fields exhibit only small perturbations. In order to fully remove these small perturbations, we recommend a combination of slowness smoothing and a contrast-dependent α to keep changes to the transmitted wave-field to a minimum. If desired, this method can also be used in conjunction with an impedance-matched wave equation, where the impedance is kept constant throughout the medium, e.g: ρ(� r) = c −1 (� r) . However, this will not allow for an independently chosen density contrast and significantly affects the transmission amplitudes. A comparison between our method and impedance matching is included in the supplementary information.
In general, the acoustic Poynting vector has shown to give an accurate estimate of the local propagation direction of wave-fields. However, problems may arise in the case of interfering waves. Firstly, we note that the use of the acoustic Poynting vector as a measure of wave-field direction breaks down for interfering waves. Secondly, and more importantly, we note that the modified Eq. (5) does not allow for reflection suppression in multiple directions simultaneously. For this reason, more sophisticated wave-field decomposition methods would not provide a solution to this issue. To remedy this, our method of contrast-dependent α allows interfering waves far away from areas exhibiting large contrasts in wave-speed to propagate unaltered. Furthermore, at high-contrast regions where interfering waves are known to appear, suppression could be turned off by setting α to zero locally.
As evident from the plane wave solution of Eqs. (3) and (4), it is possible for forward propagating components of the wave-field not exactly aligned with Ŝ to also be suppressed. The losses incurred from such misalignments, however, can be disregarded in general because of the exponential term in Eq. (6). Moreover, the results do not show any occurrence of suppression of the forward propagating wave-field.
As an alternative to the acoustic Poynting vector, a-priori ray-based methods such as Eikonal solvers can be used as a measure of the wave-field direction in the case of point-sources, by using the gradient of the shortest travel time as a time-independent propagation direction vector. It is important to note that in this way only primary arrivals are taken into account. Using this approach, the method presented here can also be applied in the frequency domain. After temporally Fourier transforming Eq. (5) we obtain a modified Helmholtz equation, which can subsequently be solved independently for each frequency component. Our experimental results using  www.nature.com/scientificreports/ be achieved by using adaptive scheme approaches 23 . Alternatively, Eq. (5) can be directly implemented using a flux-limiter 24 in a dimensional splitting approach. Results from both these methods are identical, but flux-limited schemes are significantly more expensive computationally, and thus are not preferred. The implementation of this method can also be extended naturally to the three dimensional case. Lastly, it is worth emphasising that this approach consists of an analytical modification to the acoustic wave equation. Therefore, approaches to www.nature.com/scientificreports/ numerical solutions are not limited to finite difference methods and could also be implemented using finiteelement or finite-volume methods.

Conclusion
The modification to the acoustic wave equation proposed in this paper successfully suppresses reflections within heterogeneous media, while the transmitted wave-field only exhibits small perturbations. Using a staggered grid FD scheme, the modified acoustic wave equation is implemented without significant additional computational cost. In combination with the acoustic Poynting vector, the wave-field is essentially dynamically guided through a reflection-less, heterogeneous medium. Although the solution is non-physical, this method is very suitable for use in RTM, where internal reflections often lead to imaging artifacts. Additionally, this method could serve as a tool for analysis purposes for many imaging methods, and could be used for any type of wave simulation application.