Babich's expansion and the fast Huygens sweeping method for the Helmholtz wave equation at high frequencies

Abstract

In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation so that the Helmholtz equation can be viewed as an evolution equation in one of the spatial directions. With such applications in mind, starting from Babich's expansion, we develop a new high-order asymptotic method, which we dub the fast Huygens sweeping method, for solving point-source Helmholtz equations in inhomogeneous media in the high-frequency regime and in the presence of caustics. The first novelty of this method is that we develop a new Eulerian approach to compute the asymptotics, i.e. the traveltime function and amplitude coefficients that arise in Babich's expansion, yielding a locally valid solution, which is accurate close enough to the source. The second novelty is that we utilize the Huygens–Kirchhoff integral to integrate many locally valid wavefields to construct globally valid wavefields. This automatically treats caustics and yields uniformly accurate solutions both near the source and remote from it. The third novelty is that the butterfly algorithm is adapted to accelerate the Huygens–Kirchhoff summation, achieving nearly optimal complexity O(Nlog⁡N), where N is the number of mesh points; the complexity prefactor depends on the desired accuracy and is independent of the frequency. To reduce the storage of the resulting tables of asymptotics in Babich's expansion, we use the multivariable Chebyshev series expansion to compress each table by encoding the information into a small number of coefficients. The new method enjoys the following desired features. First, it precomputes the asymptotics in Babich's expansion, such as traveltime and amplitudes. Second, it takes care of caustics automatically. Third, it can compute the point-source Helmholtz solution for many different sources at many frequencies simultaneously. Fourth, for a specified number of points per wavelength, it can construct the wavefield in nearly optimal complexity in terms of the total number of mesh points, where the prefactor of the complexity only depends on the specified accuracy and is independent of frequency. Both two-dimensional and three-dimensional numerical experiments have been carried out to illustrate the performance, efficiency, and accuracy of the method.