Far-field analytical solutions of the non-homogeneous Helmholtz and wave equations for spatially non-localized sources

Analytical solutions of partial differential equations (PDEs) are employed in many scientific disciplines, as they provide a basic understanding of physical phenomena. Furthermore, in the last few decades, with the rapid development of computational physics,^1,2 they have become an essential tool for evaluating the accuracy of numerical methods. In acoustic propagation models, explicit analytical solutions are commonly derived for PDEs with point sources³ that are formally defined as Dirac distributions. However, real sources in physical applications have a natural spatial extent and, consequently, exhibit specific directivity patterns. In addition, unlike computational approaches based on the integral formulation of a PDE, such as finite elements,² numerical methods based on a direct discretization of a PDE, for instance, finite-difference (FD) algorithms,¹ are intrinsically unable to handle point sources arbitrarily located within a grid.^4–6 When using an FD method, a point source is replaced by a regularized approximant of the Dirac distribution (a kernel function), and the numerical solution is then sought for the modified PDE. However, approximating a point source with a continuous function unavoidably introduces a regularization error that limits the accuracy of the numerical scheme irrespectively of the formal convergence order of the computational algorithm.⁵ To cope with this issue, spatially extended sources that do not require any regularization may be employed.¹ 

This study focuses on far-field analytical solutions to the non-homogeneous scalar Helmholtz and wave equations for sources with a finite aperture. Spatially extended forcing terms can be used to model realistic sources, and the analytical expressions derived in this work can be utilized in problems originally solved with point sources. As an illustration, the classical problem of acoustic scattering from elastic objects⁷ can be revisited with incident waves generated by sources with specific spatial apertures and directivity patterns. In addition, spatially extended sources can be implemented in FD codes without any regularization, and the subsequent numerical error is only due to the discretization and the properties of the chosen FD scheme. The remainder of the manuscript is organized as follows. The Helmholtz and wave equations are briefly presented in Sec. 2. The far-field analytical solution to the Helmholtz equation for spatially extended sources is described in Sec. 3. The solution to the wave equation is then derived for a particular kind of non-localized source in Sec. 4. The effects of the source spatial aperture on the acoustic far field and on the accuracy of numerical methods based on FD schemes are discussed in Sec. 5. Concluding remarks are finally drawn in Sec. 6.

The wave equation for the acoustic pressure p(r, t) due to a source $s(r,t)$ is³ 

where $r∈ℝd$ is the position vector in $d∈{1,2,3}$ spatial dimensions, $t∈ℝ$ denotes time, and $c∈ℝ>0$ is the speed of sound. Herein, we consider fixed sources and we assume that the function s is separable in space and time,

where $b(r)$ describes the spatial aperture with support $As⊆ℝd$ centered around the reference location $rc∈As$⁠, and q(t) is the temporal envelope or transmitted waveform. The function $b(r)$ is a non-negative integrable function normalized so as

As a consequence, when the support A_s is reduced to the point $rc$⁠, the aperture $b(r)$ tends, in the distributional sense, to the Dirac generalized function $δ(r−rc)$ (Ref. 8, Sec. 3.3). After applying the temporal Fourier transform, the acoustic pressure in the frequency domain f, $P(r,f)$⁠, satisfies the Helmholtz equation

where $k=2πf/c$ is the wave number.

The solution to the Helmholtz Eq. (4) is provided by the integral³ 

where $GdD(r|rs,f)$ denotes the d-dimensional Green's function, i.e., the solution to the Helmholtz Eq. (4) for $b(r)=δ(r−rs)$⁠. In the far-field, i.e., for $k‖r−rs‖2≫1$⁠, using the Fraunhofer approximation,³ the Green's function can be approximated as

where $k′=kr̂′$ is the wavenumber vector in the direction of $r̂′=(r−rc)/‖r−rc‖2$ and $rs′=rs−rc$⁠. The expressions of the Green's function term $GdD(r|rc,f)$ for $d∈{1,2,3}$ are found in Table 1. Therefore, for far-field propagation, Eq. (5) can be approximated as

In Eq. (7), the integral is the spatial Fourier transform of the aperture in the direction of $k′$⁠, i.e., the transmitter beampattern

while the product $GdD(r|rc,f)Q(f)$ is the sound field induced by a point source with temporal envelope q(t) located at the center $rc$ of the aperture. Indeed, when the support A_s is shrunk to $rc$⁠, the aperture is described by the Dirac distribution $b(rs′)=δ(rs′)$⁠, so that $B(k′)=1$⁠. Concisely, the far-field pressure $P(r,f)$ at r due to a source centered at $rc$⁠, with a waveform spectrum Q(f) and a beampattern $B(k′)$⁠, is given by the product

Sections 3.1 and 3.2 describe the waveform spectra Q(f) and beampatterns $B(k′)$ for specific source functions typically employed in underwater acoustics.

We consider two waveforms, namely the linear frequency modulated pulse (LFMP) and the Gaussian-weighted polynomial pulse (GWPP), which include the continuous wave and the Ricker wavelet, respectively, as special cases. Hereafter, we denote the pulse center as t_c and the central frequency as f_c.

An LMFP is here defined as the waveform

where τ_q is the pulse duration and $Δf$ the frequency bandwidth. A continuous wave is a special case of an LFMP with $Δf=0$⁠. The spectrum of the LFMP (10) can be written explicitly as (Ref. 9, integral 3.462.2, p. 365)

where $σ1=1/τq2−jπΔf/τq$ and $σ2=σ1*$ with the asterisk denoting complex conjugation.

A GWPP is here defined by the expression

where $n∈ℕ0$ is the degree of the polynomial with coefficients $βi∈ℝ$⁠. Note that the GWPP in Eq. (12) can be interpreted as a weighted truncated Taylor series and can therefore approximate an infinitely differentiable waveform.

The spectrum of the GWPP is (Ref. 9, integral 3.462.4, p. 365)

where H_i denotes the Hermite polynomial of order i. The Ricker wavelet q(t) is a special case of GWPP with parameters n = 2, $β0=1, β1=0$⁠, and $β2=−(2πfc)2/2$⁠, resulting in the frequency spectrum

Table 2 summarizes the waveforms considered in this work and the corresponding spectra.

We focus on typical aperture shapes encountered in underwater acoustic applications, i.e., the rectangular and the elliptical aperture, which are generalized in 3D as the parallelepipedal and the ellipsoidal or cylindrical sources, respectively. We include also the theoretically interesting case of a Gaussian aperture function.

The aperture of a parallelepipedal source is defined as

where $αl∈ℝ+, l=1,…,d$ are the sides and “$rect$” is the rectangular function. Its beampattern is given by the expression

with $sinc$ the unnormalized sine cardinal function.

The aperture of an ellipsoidal source is defined as

where $αl∈ℝ+, l=1,…,d$ are the semiaxes and Γ is the Gamma function. Its beampattern is provided by the integral

where $z$ is a scaled position vector with components $rsl′/αl, l=1,…,d$⁠. The term $rect(‖z‖22/2)$ is radial and, according to Ref. 10 (Sec. B.5), the beampattern (18) is given by the expression

where the term $k′∘α$ is the element-wise product between $k′$ and $α=(α1,…,αd)∈ℝd$⁠, and $Jν$ is the vth-order Bessel function of the first kind. The aperture of an elliptic cylindrical source is defined as

where $αl∈ℝ+, l=1,…,2$ are the semiaxes and $α3∈ℝ+$ is the cylinder's height. Its beampattern is given by

The beampattern of the Gaussian aperture

where $αl∈ℝ+, l=1,…,d$ are the standard deviations, is finally given by the expression (Ref. 9, integral 3.462.2, p. 365)

Table 3 summarizes the aperture/beampattern pairs.

The far-field sound pressure in the frequency domain $P(r,f)$ is given by Eq. (9). Depending on the application of interest, the corresponding terms for the Green's function term $GdD(r|rc,f)$⁠, the waveform Q(f) and the beampattern $B(k′)$ can be inserted from the Tables 1, 2, and 3, respectively.

The pressure signal $p(r,t)$ induced in the far-field of a spatially extended source is obtained through the inverse temporal Fourier transform of Eq. (9),

While closed-form expressions for the wavefield induced by realistic sources are generally unavailable, an explicit formula for the sound pressure in the time domain can be derived for the theoretical case of a Gaussian aperture with a Ricker waveform. Introduce the retarded time t_r,

and the characteristic frequency $f̃c$⁠,

where $r̂l′$ is the lth component of the unit vector $r̂′$⁠.

For d = 1, inserting $G1D$⁠, Q, and B from Tables 1, 2, and 3 into Eq. (24) yields (Ref. 9, integral 3.462.2, p. 365)

For d = 2, the sound pressure in the time domain is (Ref. 9, integral 3.462.1, p. 365)

where $ℜ$ indicates the real part and

with $Kν(e)$ the vth-order exponentially scaled modified Bessel function of the second kind. Finally, for d = 3, the sound pressure in the time domain is (Ref. 9, integral 3.462.2, p. 365)

By virtue of Eq. (3), the beampattern $B(k′)$ tends to 1 in the limit $f→0$⁠. Furthermore, since $b(r)$ is a non-negative integrable function, $B(k′)$ goes to 0 as $|f|→∞$ according to the Riemann-Lebesgue lemma (Ref. 10, proposition 2.2.17). Therefore, since $P(r,f)=Prc(r,f)B(k′)$⁠, the beampattern $B(k′)$ behaves as a low-pass frequency filter of the ponctual response $Prc(r,f)=GdD(r|rc,f)Q(f)$⁠. In particular, as a result of the scaling property of the Fourier transform, the broader is the source aperture, the narrower is the beampattern, and consequently the wider is the range of suppressed frequencies. This low-pass filtering effect is exemplified in Fig. 1 for a 1D Gaussian aperture with standard deviation α₁.

The analytical expressions for the far-field acoustic pressure due to spatially extended sources, as derived in this study, can serve as benchmarks for evaluating the accuracy of numerical algorithms. As an illustration, we consider a one-dimensional problem with a point source located at $rc=0.11$ mm, and transmitting a Ricker temporal waveform with parameters $c=1500$ m/s, $tc=0.2$ ms, and $fc=15$ kHz. The analytical solution p(r, t) is given by Eq. (25) with $α1=0$⁠. The numerical solution $pΔ(r,t)$ is calculated as follows: the wave equation is first reformulated as a first-order hyperbolic system and is discretized on an N-node grid with uniform spacing Δ: for $i=1,…,N, ri∈[rmin,rmax), ri=rmin+(i−1)Δ, N=⌊1+(rmax−rmin)/Δ⌋$⁠; then, the spatial derivative is computed by the standard fourth-order FD scheme and time integration is performed via the classical fourth-order Runge-Kutta method; the point source is finally implemented via the continuous approximant⁴ 

which does not vanish only at the end-nodes of the cell containing r_c. The $L1$-norm error $E1(t)$ at instant t is used to quantify the accuracy of the algorithm

In Fig. 2(a), the numerical waveform calculated at $t=tmax=2.5$ ms with $Δ=c/(20fc)$ is compared to the analytical solution. Spurious oscillations clearly contaminate the numerical result. Figure 2(b) illustrates the behaviour of $E1(tmax)$ for different values of $c/(fcΔ)$⁠. The blue curve is obtained by including all the computational nodes in the sum in Eq. (30), whereas only points r_i in the range $ri∈(3,4)$ m are employed for the red curve. In this second case, the error due to the spurious oscillations is not accounted for in $E1(tmax)$⁠. Nevertheless, although the algorithm is formally fourth-order accurate in both space and time, $E1$ goes to zero as $Δ2$ at best. This result can be explained by noting that $δΔ(r−rc)$ is only a second-order accurate approximation of $δ(r−rc)$⁠. More specifically, as highlighted by⁶ 

for any infinitely differentiable compactly supported test function $ϕ$⁠. As a consequence, the main contribution to $E1$ is not due to the FD discretization but instead to the approximation $δΔ$ of the Dirac distribution δ. To solve this issue, higher-order discretizations of δ have been proposed,^4–6 with the works on particle methods leading the research field.¹¹ The present work provides an alternative approach for validating numerical schemes. Spatially extended Gaussian sources can indeed be implemented in FD-based methods without any regularization. To clarify this point, the one-dimensional wave equation is reconsidered with an extended source centered at r_c and with Gaussian aperture of width $α1=c/(5fc)$⁠. The analytical expression of p(r, t) is given by Eq. (25). Two numerical solutions are computed: in both simulations, the classical fourth-order Runge-Kutta method is again employed for time integration; the spatial derivative is now calculated with the standard second- and fourth-order centered FD schemes. Figure 3 depicts the error $E1$ for different values of Δ and clearly demonstrates that $E1$ behaves as desired according to the formal order of the FD scheme.

Explicit analytical expressions for the acoustic far-field induced by spatially extended sources have been derived. In particular, a general far-field solution to the non-homogeneous Helmholtz equation has been developed; far field solutions to the wave equation have been provided for the theoretical case of sources with Gaussian spatial aperture and Ricker temporal envelope. These formulas can be useful in numerous acoustic applications as well as in computational physics, especially for evaluating the accuracy of numerical methods based on a direct discretization of the Helmholtz and wave equations.

This work was performed under Project SAC000C04 of the STO-CMRE Programme of Work, funded by the NATO Allied Command Transformation.