Non-localized impulsive sources are ubiquitous in underwater acoustic applications. However, analytical expressions of their acoustic field are usually not available. In this work, far-field analytical solutions of the non-homogeneous scalar Helmholtz and wave equations are developed for a class of spatially extended impulsive sources. The derived expressions can serve as benchmarks to verify the accuracy of numerical solvers.

## 1. Introduction

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^{3} 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,^{2} numerical methods based on a direct discretization of a PDE, for instance, finite-difference (FD) algorithms,^{1} 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.^{5} To cope with this issue, spatially extended sources that do not require any regularization may be employed.^{1}

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^{7} 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.

## 2. Problem statement

The wave equation for the acoustic pressure *p*(r, t) due to a source $s(r,t)$ is^{3}

where $r\u2208\mathbb{R}d$ is the position vector in $d\u2208{1,2,3}$ spatial dimensions, $t\u2208\mathbb{R}$ denotes time, and $c\u2208\mathbb{R}>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\u2286\mathbb{R}d$ centered around the reference location $rc\u2208As$, 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 $\delta (r\u2212rc)$ (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\pi f/c$ is the wave number.

## 3. Acoustic pressure in the frequency domain

where $GdD(r|rs,f)$ denotes the *d*-dimensional Green's function, i.e., the solution to the Helmholtz Eq. (4) for $b(r)=\delta (r\u2212rs)$. In the far-field, i.e., for $k\Vert r\u2212rs\Vert 2\u226b1$, using the Fraunhofer approximation,^{3} the Green's function can be approximated as

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

Space dimension . | $GdD(r|rc,f)$ . |
---|---|

d = 1 | $jejk\Vert r\u2212rc\Vert 22k$ |

d = 2 | $+j4H0(1)(|k|\Vert r\u2212rc\Vert 2),k>0$ |

$\u2212j4H0(2)(|k|\Vert r\u2212rc\Vert 2),k<0$ | |

$\u2243sgn(k)j42\pi |k|\Vert r\u2212rc\Vert 2ejk\Vert r\u2212rc\Vert 2\u2212sgn(k)j\pi /4$ | |

d = 3 | $ejk\Vert r\u2212rc\Vert 24\pi \Vert r\u2212rc\Vert 2$ |

Space dimension . | $GdD(r|rc,f)$ . |
---|---|

d = 1 | $jejk\Vert r\u2212rc\Vert 22k$ |

d = 2 | $+j4H0(1)(|k|\Vert r\u2212rc\Vert 2),k>0$ |

$\u2212j4H0(2)(|k|\Vert r\u2212rc\Vert 2),k<0$ | |

$\u2243sgn(k)j42\pi |k|\Vert r\u2212rc\Vert 2ejk\Vert r\u2212rc\Vert 2\u2212sgn(k)j\pi /4$ | |

d = 3 | $ejk\Vert r\u2212rc\Vert 24\pi \Vert r\u2212rc\Vert 2$ |

In Eq. (7), the integral is the spatial Fourier transform of the aperture in the direction of $k\u2032$, 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\u2032)=\delta (rs\u2032)$, so that $B(k\u2032)=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\u2032)$, is given by the product

### 3.1 Waveforms

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}#### 3.1.1 LFMP

An LMFP is here defined as the waveform

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

where $\sigma 1=1/\tau q2\u2212j\pi \Delta f/\tau q$ and $\sigma 2=\sigma 1*$ with the asterisk denoting complex conjugation.

#### 3.1.2 GWPP

A GWPP is here defined by the expression

where $n\u2208\mathbb{N}0$ is the degree of the polynomial with coefficients $\beta i\u2208\mathbb{R}$. 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, $\beta 0=1,\u2009\beta 1=0$, and $\beta 2=\u2212(2\pi fc)2/2$, resulting in the frequency spectrum

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

Pulse . | q(t)
. | Q(f)
. |
---|---|---|

LFMP | $e\u2212(t\u2212tc)2/\tau q2\u2009cos\u2009(2\pi fc(t\u2212tc)+\pi \Delta f\tau q(t\u2212tc)2)$ | $\u2211i=12\pi 4\sigma ie[(\u22121)i2\pi 2fc/\sigma i+j2\pi tc]f\u2212\pi 2(f2+fc2)/\sigma i$ |

GWPP | $e\u2212\pi 2fc2(t\u2212tc)2\u2211i=0n\beta i(t\u2212tc)i$ | $2\pi ej2\pi ftc\u2212(f/fc)2\u2211i=0nj\u2212i\beta i(2\pi fc)\u2212(i+1)Hi(\u2212ffc)$ |

Pulse . | q(t)
. | Q(f)
. |
---|---|---|

LFMP | $e\u2212(t\u2212tc)2/\tau q2\u2009cos\u2009(2\pi fc(t\u2212tc)+\pi \Delta f\tau q(t\u2212tc)2)$ | $\u2211i=12\pi 4\sigma ie[(\u22121)i2\pi 2fc/\sigma i+j2\pi tc]f\u2212\pi 2(f2+fc2)/\sigma i$ |

GWPP | $e\u2212\pi 2fc2(t\u2212tc)2\u2211i=0n\beta i(t\u2212tc)i$ | $2\pi ej2\pi ftc\u2212(f/fc)2\u2211i=0nj\u2212i\beta i(2\pi fc)\u2212(i+1)Hi(\u2212ffc)$ |

### 3.2 Beampatterns

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.

#### 3.2.1 Parallelepipedal aperture

The aperture of a parallelepipedal source is defined as

where $\alpha l\u2208\mathbb{R}+,\u2009l=1,\u2026,d$ are the sides and “$rect$” is the rectangular function. Its beampattern is given by the expression

with $sinc$ the unnormalized sine cardinal function.

#### 3.2.2 Ellipsoidal and elliptic cylindrical apertures

The aperture of an ellipsoidal source is defined as

where $\alpha l\u2208\mathbb{R}+,\u2009l=1,\u2026,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\u2032/\alpha l,\u2009l=1,\u2026,d$. The term $rect(\Vert z\Vert 22/2)$ is radial and, according to Ref. 10 (Sec. B.5), the beampattern (18) is given by the expression

where the term $k\u2032\u2218\alpha $ is the element-wise product between $k\u2032$ and $\alpha =(\alpha 1,\u2026,\alpha d)\u2208\mathbb{R}d$, and $J\nu $ is the *v*th-order Bessel function of the first kind. The aperture of an elliptic cylindrical source is defined as

where $\alpha l\u2208\mathbb{R}+,\u2009l=1,\u2026,2$ are the semiaxes and $\alpha 3\u2208\mathbb{R}+$ is the cylinder's height. Its beampattern is given by

#### 3.2.3 Gaussian aperture

The beampattern of the Gaussian aperture

where $\alpha l\u2208\mathbb{R}+,\u2009l=1,\u2026,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.

Aperture . | $b(rs\u2032)$ . | $B(k\u2032)$ . |
---|---|---|

Parallelepipedal | $\u220fl=1d1\alpha lrect(rsl\u2032\alpha l)$ | $\u220fl=1dsinc(kl\u2032\alpha l2)$ |

Ellipsoidal | $\Gamma (1+d2)(\u220fl=1d1\pi 1/2\alpha l)rect(12\u2211l=1drsl\u20322\alpha l2)$ | $2d/2\Gamma (1+d2)Jd/2(\Vert k\u2032\u2218\alpha \Vert 2)\Vert k\u2032\u2218\alpha \Vert 2d/2$ |

Cylindrical | $1\pi (\u220fl=131\alpha l)rect(12\u2211l=12rsl\u20322\alpha l2)rect(rs3\u2032\alpha 3)$ | $2\u2009J1(\Vert k\u2032\u2218\alpha \Vert 2)\Vert k\u2032\u2218\alpha \Vert 2\u2009sinc(k3\u2032\alpha 32)$ |

Gaussian | $1(2\pi )d\u220fl=1d1\alpha le\u2212rsl\u20322\u200a/\u200a(2\alpha l2)$ | $\u220fl=1de\u2212(kl\u2032\alpha l)2/2$ |

Aperture . | $b(rs\u2032)$ . | $B(k\u2032)$ . |
---|---|---|

Parallelepipedal | $\u220fl=1d1\alpha lrect(rsl\u2032\alpha l)$ | $\u220fl=1dsinc(kl\u2032\alpha l2)$ |

Ellipsoidal | $\Gamma (1+d2)(\u220fl=1d1\pi 1/2\alpha l)rect(12\u2211l=1drsl\u20322\alpha l2)$ | $2d/2\Gamma (1+d2)Jd/2(\Vert k\u2032\u2218\alpha \Vert 2)\Vert k\u2032\u2218\alpha \Vert 2d/2$ |

Cylindrical | $1\pi (\u220fl=131\alpha l)rect(12\u2211l=12rsl\u20322\alpha l2)rect(rs3\u2032\alpha 3)$ | $2\u2009J1(\Vert k\u2032\u2218\alpha \Vert 2)\Vert k\u2032\u2218\alpha \Vert 2\u2009sinc(k3\u2032\alpha 32)$ |

Gaussian | $1(2\pi )d\u220fl=1d1\alpha le\u2212rsl\u20322\u200a/\u200a(2\alpha l2)$ | $\u220fl=1de\u2212(kl\u2032\alpha l)2/2$ |

### 3.3 Acoustic far-field in the frequency domain

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\u2032)$ can be inserted from the Tables 1, 2, and 3, respectively.

## 4. Acoustic pressure in the time domain

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\u0303c$,

where $r\u0302l\u2032$ is the *l*th component of the unit vector $r\u0302\u2032$.

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 $\u211c$ indicates the real part and

with $K\nu (e)$ the *v*th-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)

## 5. Applications

### 5.1 Influence of the source spatial aperture on the far-field acoustic pressure

By virtue of Eq. (3), the beampattern $B(k\u2032)$ tends to 1 in the limit $f\u21920$. Furthermore, since $b(r)$ is a non-negative integrable function, $B(k\u2032)$ goes to 0 as $|f|\u2192\u221e$ according to the Riemann-Lebesgue lemma (Ref. 10, proposition 2.2.17). Therefore, since $P(r,f)=Prc(r,f)B(k\u2032)$, the beampattern $B(k\u2032)$ 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 *α*_{1}.

### 5.2 Source representation in FD algorithms

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 $\alpha 1=0$. The numerical solution $p\Delta (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,\u2026,N,\u2009ri\u2208[rmin,rmax),\u2009ri=rmin+(i\u22121)\Delta ,\u2009N=\u230a1+(rmax\u2212rmin)/\Delta \u230b$; 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^{4}

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 $\Delta =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\Delta )$. 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\u2208(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 $\Delta 2$ at best. This result can be explained by noting that $\delta \Delta (r\u2212rc)$ is only a second-order accurate approximation of $\delta (r\u2212rc)$. More specifically, as highlighted by

^{6}

for any infinitely differentiable compactly supported test function $\varphi $. As a consequence, the main contribution to $E1$ is not due to the FD discretization but instead to the approximation $\delta \Delta $ 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.^{11} 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 $\alpha 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.

## 6. Concluding remarks

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.

## Acknowledgment

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