Abstract

The transmission of information by propagating or diffusive waves is common to many fields of engineering and physics. Such physical phenomena are governed by a Helmholtz (real wavenumber) or pseudo-Helmholtz (complex wavenumber) equation. Since these equations are linear, it would be useful to be able to use tools from signal theory in solving related problems. The aim of this paper is to derive multidimensional input/output transfer function relationships in the spatial domain for these equations in order to permit such a signal theoretic approach to problem solving. This paper presents such transfer function relationships for the spatial (not Fourier) domain within appropriate coordinate systems. It is shown that the relationships assume particularly simple and computationally useful forms once the appropriate curvilinear version of a multidimensional spatial Fourier transform is used. These results are shown for both real and complex wavenumbers. Fourier inversion of these formulas would have applications for tomographic problems in various modalities. In the case of real wavenumbers, these inversion formulas are presented in closed form, whereby an input can be calculated from a given or measured wavefield.

1. Introduction

The transmission of information over space and time is often governed by the theory of waves. Many important physical phenomena are described by a Helmholtz equation, for example, in the fields of electromagnetism, acoustics, and optics [14]. Other physical phenomena are well described by the use of damped or diffusion waves, such as the propagation of heat or photonic waves. For example, ultrasound has long been used for medical imaging [3] while the use of optical radiation via diffuse photon density waves for imaging inhomogeneities in turbid media is a newer development [58]. Similarly, photothermal tomographic imaging methods have also been used for nondestructive evaluation [912], and more recently for biomedical imaging [13, 14]. Furthermore, a combination of ultrasonic and optical techniques known as photoacoustics has also gained prominence in recent years [4, 15].

Modelled as a linear process, the propagation and scattering of waves obey the principles of superposition and homogeneity. One common method of studying linear processes is to view them as linear systems with an input and output. The input/output relationships of the system can then be characterized via the system transfer function, which implies interpreting the problem as a signal theory problem and thus enabling the application of powerful system theoretic results and concepts. This type of approach applied to propagating or diffusion wave problems has already yielded promising results [1618] and the goal of this paper is to expand the scope of this approach.

The seeds of such an approach to the solution of tomographic-type problems can be seen in [1, 8, 10, 19, 20]; however, in those papers there is no clear reference to any such system theoretic concepts. As such, a unifying framework for a signal theoretic approach to wavefields does not exist. The goal of this work is to propose a clear and cohesive signal theoretic framework for working with acoustic, thermal, photonic, or other wavefield-related problems. This would be particularly powerful, enabling the application of a vast array of results from a signal processing point of view and in particular, recent novel results from algebraic signal processing [2127] should find ground-breaking applications.

With these ideas in mind, the goal of this exposition is to derive multidimensional spatial transfer function relationships for the Helmholtz and pseudo-Helmholtz equation so that the tools of linear system theory can be used to solve related problems.

We consider the input quantity to be the time- and space-dependent inhomogeneity term (the “forcing” term in the relevant partial differential equation). The output quantity is then taken as the resulting field present at any point in space. The goal is to develop transfer function relationships in spatial multidimensions so that general results of the theory of signals and systems may be used. It is known that the relationship between the input (inhomogeneity term) and the output (resulting wavefield) is one of convolution with the impulse response (Green's function in classical terms). In the Fourier domain, the equivalent relationship is a transfer function relationship, that is to say one of multiplication. In this paper, key transfer function relationships are shown to hold for the spatial (not Fourier) domain with the use of appropriate coordinate systems. That is to say, the relationship between input and output is shown to be one of multiplication (not convolution) even in the spatial domain, hence yielding the signal theoretic construction that is the aim of the paper.

Rather than focusing on a specific application area (acoustics, medical imaging, seismic, etc.), the end goal of this paper is a unifying framework. The backbone of this framework is Fourier transforms in multidimensions which are also transcribed into curvilinear coordinate systems. The representations in this paper are embedded in a time-dependent, spatially one- to three-dimensional description.

The outline of this paper is as follows. Sections 2 to 7 establish the background preliminaries, definitions, and sign conventions that are important to the developments in the rest of the paper. The definition of multidimensional Fourier transforms is given in Section 2. Section 3 presents the mathematical forms of the types of signal to which this development applies, namely, any signal governed by a Helmholtz or pseudo-Helmholtz equation. Section 4 introduces the Green's and transfer functions for the (generalized) Helmholtz equation and points out that the spatial transfer function assumes a particularly simple form in curvilinear coordinates. Section 5 explains the sign conventions that will be employed in this paper to ensure consistency and simplicity of results. Sections 6 and 7 develop multidimensional Fourier transforms in curvilinear coordinates. Sections 8, 9, and 10 present the relevant transfer function relationships in curvilinear coordinates in 3D, 2D, and 1D, respectively, and form the core results of this paper. Section 11 summarizes the results of the previous three sections in a tabular format that is easy to refer to. Sections 12 and 13 present some applications while Section 14 concludes the paper.

2. Fourier Transforms in 𝑛-Dimensions

The theory of Fourier transforms can be extended to 𝑛 in a way that is completely analogous to the treatment of one-dimensional Fourier transforms. Let elements of 𝑛 be denoted by 𝑥=(𝑥1,𝑥2,,𝑥𝑛) or more generally 𝑟, and elements of the corresponding Fourier dual space Ω𝑛 be denoted by 𝜔=(𝜔1,𝜔2,,𝜔𝑛). Under the suitability of integration of 𝑓, then the Fourier transform in 𝑛 is defined for 𝜔in Ω𝑛 as𝐹=𝜔𝑛𝑓𝑒𝑥𝑗𝜔𝑥𝑑𝑥.(2.1) Under suitable conditions, the function 𝑓 can be recovered from the inverse transform through 𝑓=1𝑥(2𝜋)𝑛Ω𝑛𝐹𝑒𝜔𝑗𝜔𝑥𝑑𝜔.(2.2) Various other conventions are possible regarding the location of the positive and negative signs and also the factors of 2𝜋 in (2.1) and (2.2).

3. The Helmholtz and Pseudo-Helmholtz Equation

All wave fields governed by the wave equation (such as acoustic or electromagnetic waves) lead to the Helmholtz equation once a Fourier transform is used to transform the time domain to the frequency domain:2𝑢𝑟,𝜔+𝑘2𝑠𝑢𝑟,𝜔=𝑠𝑟,𝜔,(3.1) where 𝑘2𝑠=𝜔2/𝑐2𝑠. Here, 𝑠(𝑟,𝜔) is the temporal Fourier transform of the inhomogeneous time- and space-dependent source term for the wave equation. From a signal theory perspective, this is considered to be the input to the system. The variable 𝑢(𝑟,𝜔) represents a physical variable that is governed by the wave equation, for example, acoustic pressure, and is considered to be the output from a signal theory point of view. Both 𝑠(𝑟,𝜔) and 𝑢(𝑟,𝜔) are functions of position, 𝑟, and (temporal) frequency, 𝜔. The variable 𝑐𝑠 represents the speed of the wave, which for an acoustic wave would be the speed of sound and for an electromagnetic wave would be the speed of light. For wavefields governed by the wave equation, 𝑘2𝑠 is a real (and positive) quantity.

Other physical phenomena lead to a “pseudo-” Helmholtz equation upon Fourier transformation of the time variable to a (temporal) frequency variable. For example, the equation for a diffuse photon density wave (DPDW) which describes the photon density 𝑢(𝑟,𝑡) in a solid due to incident energy intensity 𝑠(𝑟,𝑡) (optical source function) is given in the time domain by [2] D2𝑢1𝑟,𝑡𝑐𝜕𝑢𝜕𝑡𝑟,𝑡𝜇𝑎𝑢𝑟,𝑡=𝑠𝑟,𝑡.(3.2)

In the above equation, 𝜇𝑎 is the optical absorption coefficient (m−1), 𝑐 is the speed of light in the turbid medium (m/s), and D is the optical diffusion coefficient (m). Taking the temporal Fourier transform (FT) of (3.2) to transform time to frequency gives the following:2𝑢𝑟,𝜔+𝑘2𝑝𝑢𝑠𝑟,𝜔=𝑟,𝜔D,(3.3) where 𝑘2𝑝=(𝜇𝑎/D)(𝑖𝜔/𝑐D) is the complex wavenumber. Equation (3.3) is a pseudo-Helmholtz equation. The term “pseudo” is used because although (3.3) has the form of a Helmholtz equation similar to (3.1), the wavenumber 𝑘2𝑝 is a complex variable, with the imaginary part indicating a decaying or damped wave. In this development the terms “Helmholtz” and “pseudo-Helmholtz” will generally be used interchangeably unless otherwise noted. Similarly, the standard heat equation is given by 2𝑢1𝑟,𝑡𝛼𝜕𝑢𝜕𝑡𝑟,𝑡=𝑠𝑟,𝑡,(3.4) where 𝑠(𝑟,𝑡) is the time- and space-dependent heat source, 𝑢(𝑟,𝑡) describes the temperature in the material as a function of time and space, and 𝛼 is the thermal diffusivity of a material. As done previously, we take the temporal FT of (3.4) to obtain2𝑢𝑟,𝜔+𝑘2𝑡𝑢𝑟,𝜔=𝑠𝑟,𝜔,(3.5) where 𝑘2𝑡=𝑖𝜔/𝛼 is the complex wavenumber. Once again (3.5) has the form of a Helmholtz equation, although it is a pseudo-Helmholtz equation to be exact since the wavenumber is complex and indicates a damped wave.

Thus, we see that the general Helmholtz form of (3.1) can be used to describe several different physical phenomena, from the propagation of light or acoustic waves to the heavily damped nature of photonic or thermal waves. The exact form of the wavenumber in each case is the best indicator as to the propagation characteristics of a wave with a real 𝑘2 indicating a propagating wave and a complex 𝑘2 indicating a damped wave.

4. Green's Function and Transfer Function for the Helmholtz Equation

Taking the full spatial Fourier transform of the Helmholtz (or pseudo-Helmholtz) equation and rearranging yields𝑈𝜔𝜔,𝜔=𝑈𝑥,𝜔𝑦,𝜔𝑧=1,𝜔𝜔2𝑥+𝜔2𝑦+𝜔2𝑧𝑘2𝑆𝜔,𝜔,(4.1) where we have used the shorthand notation (𝜔)=(𝜔𝑥,𝜔𝑦,𝜔𝑧) to denote a point in 3D spatial Fourier frequency space. The wavenumber may be real or complex. A capital letter is used to denote the full 3D spatial and temporal Fourier transform of a function, although it should be clear from the arguments which is being indicated. In cases where it may not be clear, a tilde () will be used to denote the function in spatial Fourier space. For shorthand notation, let 𝜔2𝑘=𝜔2𝑥+𝜔2𝑦+𝜔2𝑧, so that 𝜔𝑘 is the length of the spatial Fourier vector. The Fourier space has the same spatial frequency dimension as the spatial dimension of the problem. So in 2D space, a 𝜔𝑧 spatial frequency variable would not be required and so forth.

By inverse spatial Fourier transformation of (4.1), the wavefield in 𝑛 dimensions is given by𝑢=1𝑟,𝜔(2𝜋)𝑛𝑆𝜔,𝜔𝜔2𝑘𝑘2𝑒𝑖𝜔𝑟𝑑𝜔.(4.2) Using the definition of the Fourier transform of 𝑠(𝑟,𝜔), the above equation can be rewritten as𝑢=1𝑟,𝜔(2𝜋)𝑛𝑒𝑖𝜔𝑟𝜔2𝑘𝑘2𝑠𝑒𝑥,𝜔𝑖𝜔𝑥𝑑𝑥𝑑𝜔.(4.3) Let us define the spatial Green’s function in 𝑛-dimensional space as𝑔=1𝑟𝑥,𝜔=𝑔𝑟𝑥,𝜔(2𝜋)𝑛𝑒𝑖𝜔(𝑟𝑥)𝜔2𝑘𝑘2𝑑𝜔.(4.4) This is the Green's function for the (pseudo-) Helmholtz equation. By switching the order of integration, (4.3) can be rewritten such that it can be clearly interpreted as a spatial convolution of Green's function with the input source:𝑢=𝑟,𝜔𝑠𝑔𝑥,𝜔𝑟𝑥,𝜔𝑑𝑥=𝑠𝑟,𝜔𝑟𝑔𝑟,𝜔.(4.5) The notation 𝑟 has been used to denote a (multidimensional) space-only convolution. The frequency domain equivalent to (4.5) is (4.1), which can be interpreted as a multiplication in terms of the spatial transfer function as 𝑈𝑆𝜔,𝜔=𝐺𝜔,𝜔𝜔,𝜔,(4.6) where the Fourier transform of the spatial Green's function is the spatial transfer function:𝐺𝑔=1𝜔,𝜔=𝔽𝑟,𝜔𝜔2𝑘𝑘2.(4.7) The dependence of the spatial Green's function on temporal frequency is via the wavenumber 𝑘 of the wavefield (which may be real or complex) and this is emphasized by writing 𝐺(𝜔,𝜔) as a function of spatial and temporal frequencies.

4.1. Spatial Transfer Function in Curvilinear Coordinates

Equations (4.4) and (4.7) hold in any 𝑛-dimensional space with Green's function of (4.4) being the 𝑛-dimensional spatial inverse Fourier transform of the spatial transfer function of (4.7). More specifically, in 1D, 𝜔2𝑘=𝜔2𝑥, and in 2D 𝜔2𝑘=𝜔2𝑥+𝜔2𝑦. Transforming to polar coordinates in the 2D spatial frequency domain via the transformation 𝜔𝑥=𝜌cos𝜓, 𝜔𝑦=𝜌sin𝜓 gives 𝜔2𝑘=𝜔2𝑥+𝜔2𝑦=𝜌2 where 𝜌 is the radial frequency variable in polar coordinates. Similarly, in 3D, 𝜔2𝑘=𝜔2𝑥+𝜔2𝑦+𝜔2𝑧 becomes 𝜔2𝑘=𝜔2𝑥+𝜔2𝑦+𝜔2𝑧=𝜌2 where 𝜌 is the spherical radial frequency variable in spherical polar coordinates in spatial frequency space via the coordinate transformation given by 𝜔𝑥=𝜌sin𝜓𝜔cos𝜃𝜔, 𝜔𝑦=𝜌sin𝜓𝜔sin𝜃𝜔 and 𝜔𝑧=𝜌cos𝜓𝜔. Thus, the spatial transfer functions assume particularly simple, radially symmetric forms in polar coordinates for 2D space and spherical polar coordinates for 3D space. These particularly simple forms for the 2D and 3D transfer functions in polar coordinates motivates the development of the general 2D and 3D Fourier transforms in curvilinear coordinates [28, 29] so that this simple form of the transfer function may be exploited.

5. Notation and Sign Conventions

Wavenumbers are defined by the squares of their quantities and arise as a result of taking the Fourier transform of the corresponding propagation equation (be it acoustic, thermal, or otherwise), which in turn leads to a Helmholtz or pseudo-Helmholtz equation. These wavenumbers are defined so that the wave propagation equation transformed to the temporal frequency domain all have the form of a pseudo-Helmholz equation as given by (3.1) or (3.5) with a complex or real wavenumber. For the rest of the paper, we will consider that 𝑘𝑡 represents a generic complex wavenumber while 𝑘𝑠 represents a generic real wavenumber.

Rather than the squared wavenumber, the quantity of interest will prove to be the wavenumber itself, namely, 𝑘𝑠 and 𝑘𝑡, which are the square roots of the given squared wavenumber in the Helmholtz equation. Each 𝑘 can be considered as the sum of a real and an imaginary part, so that 𝑘𝑡=𝑘𝑡𝑟+𝑖𝑘𝑡𝑖 with 𝑘𝑡𝑟 denoting the real part of 𝑘𝑡 and 𝑘𝑡𝑖 denoting the imaginary part. Since the square root of any 𝑘2 can be ±𝑘, we will use the convention that for a given complex 𝑘, the required square root of the corresponding 𝑘2 is defined such that the imaginary part of 𝑘 is negative. Hence, 𝑘𝑡 is chosen as the square root of 𝑘2𝑡 such that 𝑘𝑡𝑖<0 and so forth. If this sign convention is adopted, then it was shown in [36] that the many results for complex wavenumber use the same mathematical form of travelling wave solution as for the real wavenumbers.This makes the notation and book-keeping considerably simpler. It is noted that this sign convention is the opposite of what this author adopts in [30].

5.1. Sommerfeld Radiation Condition

To aid in the selection of a causal solution, the Sommerfeld radiation condition is required. The Sommerfeld radiation condition states that the sources in the field must be sources not sinks of energy. Therefore, energy radiated from sources must scatter to infinity and cannot radiate from infinity into the field. Mathematically, a solution 𝑢(𝑥), where 𝑥 is the spatial variable, to the Helmholtz equation is considered to be radiating if it satisfieslim|𝑥||𝑥|(𝑛1)/2𝜕𝜕|𝑥|+𝑖𝑘𝑢(𝑥)=0,(5.1) where 𝑛 is the dimension of the space and 𝑘 is the wavenumber in the Helmholtz equation. This is the radiation equation based on an implied time variation of 𝑒𝑖𝜔𝑡 which is implicit in our chosen definition of the Fourier transform, (2.1). For example, with the standard definition of the Fourier transform, the transform of 𝑓(𝑡) is 𝑖𝜔𝐹(𝜔), clearly implying the 𝑒𝑖𝜔𝑡 dependence. Had a different definition of the Fourier transform been used, the implied time variation would have been 𝑒𝑖𝜔𝑡, in which case the sign of 𝑖 in (5.1) would be reversed.

6. 2D Fourier Transforms in Terms of Polar Coordinates

Given the desire to exploit the simple form of the spatial transfer function in curvilinear coordinates, we consider 2D and 3D Fourier transform in terms of curvilinear coordinates. Let us first consider the Fourier transforms in 2D in polar coordinates. In order to define this, some preliminary definitions of Hankel transforms are required first.

6.1. Hankel Transforms

The Hankel transform of order 𝑛 is defined by the integral [31]𝐹𝑛(𝜌)=𝑛(𝑓(𝑟))=0𝑓(𝑟)𝐽𝑛(𝜌𝑟)𝑟𝑑𝑟,(6.1) where 𝐽𝑛(𝑧) is the 𝑛th order Bessel function. If 𝑛>1/2, the transform is self-reciprocating and the inversion formula is given by 𝑓(𝑟)=𝑛1𝐹𝑛=(𝜌)0𝐹𝑛(𝜌)𝐽𝑛(𝜌𝑟)𝜌𝑑𝜌.(6.2)

6.2. Connection between the 2D Fourier Transform and Hankel Transform

To develop Fourier transforms in 2D in polar coordinates, both the function 𝑓(𝑟) and its Fourier transform 𝐹(𝜔) are expressed in polar coordinates. In general 𝑓(𝑟)=𝑓(𝑟,𝜃) is not radially symmetric and is a function of both 𝑟 and 𝜃 so that the 𝜃 dependence can be expanded into a Fourier series due to the 2𝜋 periodicity of the function in 𝜃:𝑓(𝑟,𝜃)=𝑛=𝑓𝑛(𝑟)𝑒𝑗𝑛𝜃,(6.3) where the Fourier coefficients 𝑓𝑛(𝑟) can be found from𝑓𝑛1(𝑟)=2𝜋02𝜋𝑓(𝑟,𝜃)𝑒𝑗𝑛𝜃𝑑𝜃.(6.4) Similarly, the 2D Fourier transform 𝐹(𝜔)=𝐹(𝜌,𝜓) can also be expanded into its own Fourier series so that𝐹(𝜌,𝜓)=𝑛=𝐹𝑛(𝜌)𝑒𝑗𝑛𝜓,(6.5) and where those Fourier coefficients are similarly found from𝐹𝑛1(𝜌)=2𝜋02𝜋𝐹(𝜌,𝜓)𝑒𝑗𝑛𝜓𝑑𝜓,(6.6)

The relationship between the Fourier coordinates in normal space 𝑓𝑛(𝑟) in (6.3) and the spatial frequency space Fourier coordinates 𝐹𝑛(𝜌) in (6.5) is desired. The details of this are given in [28], and the results are summarized here. It is emphasized that the relationship between 𝑓𝑛(𝑟) and 𝐹𝑛(𝜌) is not a Fourier transform. The relationship between them is given by𝐹𝑛(𝜌)=2𝜋𝑖𝑛0𝑓𝑛(𝑟)𝐽𝑛(𝜌𝑟)𝑟𝑑𝑟=2𝜋𝑖𝑛𝑛𝑓𝑛(𝑟),(6.7) where 𝑛 is the 𝑛th order Hankel transform. The reverse relationship is given by𝑓𝑛𝑖(𝑟)=𝑛2𝜋0𝐹𝑛(𝜌)𝐽𝑛𝑖(𝜌𝑟)𝜌𝑑𝜌=𝑛2𝜋𝑛𝐹𝑛(𝜌).(6.8) The 𝑛th term in the Fourier series for the original function will Hankel transform into the 𝑛th term of the Fourier series of the Fourier transform. However, it is an 𝑛th order Hankel transform for the 𝑛th term, namely, all the terms are not equivalently transformed. The mapping from 𝑓𝑛(𝑟) to 𝐹𝑛(𝜌) is one of 𝑛th-order Hankel transform, which in general is not a 2D Fourier transform.

The operation of taking the 2D Fourier transform of a function is thus equivalent to (1) first finding its Fourier series expansion in the angular variable and (2) then finding the 𝑛th-order Hankel transform (of the radial variable to the spatial radial variable) of the 𝑛th coefficient in the Fourier series. Since each of these operations involves integration over one variable only with the others being considered parameters vis-à-vis the integration, the order in which these operations are performed is interchangeable.

7. Spherical Hankel Transform

We introduce the spherical Hankel transform, which will form part of the 3D Fourier transform in spherical coordinates.

7.1. Definition of the Spherical Hankel Transform

The spherical Hankel transform can then be defined as [32, 33]𝐹𝑛(𝜌)=𝕊𝑛{𝕗𝑓(𝑟)}=0𝑓(𝑟)𝑗𝑛(𝜌𝑟)𝑟2𝑑𝑟.(7.1)𝕊𝑛 is used to specifically denote the spherical Hankel transform of order 𝑛. The inverse transform is given by [29]2𝑓(𝑟)=𝜋0𝐹𝑛(𝜌)𝑗𝑛(𝜌𝑟)𝜌2𝑑𝜌.(7.2) The spherical Hankel transform is particularly useful for problems involving spherical symmetry.

7.2. Spherical Harmonics

The spherical harmonics are the solution to the angular portion of Laplace's equation in spherical polar coordinates and can be shown to be orthogonal. These spherical harmonics are given by [32]𝑌𝑚𝑙(𝜓,𝜃)=(2𝑙+1)(𝑙𝑚)!𝑃4𝜋(𝑙+𝑚)!𝑚𝑙(cos𝜓)𝑒𝑖𝑚𝜃,(7.3) where 𝑌𝑚𝑙 is called a spherical harmonic function of degree 𝑙 and order 𝑚, 𝑃𝑚𝑙 is an associated Legendre function, 0𝜓𝜋 represents the colatitude and 0𝜃2𝜋 represents the longitude. With the normalization of the spherical harmonics as given in (7.3), the spherical harmonics are orthonormal so that 02𝜋𝜋0𝑌𝑚𝑙𝑌𝑚𝑙sin𝜓𝑑𝜓𝑑𝜃=𝛿𝑙𝑙𝛿𝑚𝑚.(7.4) Here 𝛿𝑖𝑗 is the kronecker delta and the overbar indicates the complex conjugate. It is important to note that there are several different normalizations of the spherical harmonics that are possible, that will differ from (7.3) and thus lead to a different version of (7.4). Several of these are nicely catalogued in [34], including which disciplines tend to use which normalization. This is important to note since any result that uses orthogonality will differ slightly depending on the choice of normalization.

The spherical harmonics form a complete set of orthonormal functions and thus form a vector space. When restricted to the surface of a sphere, functions may be expanded on the sphere into a series approximation much like a Fourier series. This is in fact a spherical harmonic series. Any square-integrable function may be expanded as 𝑓(𝜓,𝜃)=𝑙𝑙=0𝑚=𝑙𝑓𝑚𝑙𝑌𝑚𝑙(𝜓,𝜃),(7.5) where 𝑓𝑚𝑙=02𝜋𝜋0𝑓(𝜓,𝜃)𝑌𝑚𝑙(𝜓,𝜃)sin𝜓𝑑𝜓𝑑𝜃.(7.6) The coefficients 𝑓𝑚𝑙 are sometimes referred to as the spherical Fourier transform of 𝑓(𝜓,𝜃) [35].

7.3. 3D Fourier Transforms in Spherical Polar Coordinates

To find 3D Fourier transforms in spherical polar coordinates, both the function and its 3D Fourier transform are written in terms of polar coordinates in the spatial and spatial frequency domains. That is, a function in 3D space is expressed as 𝑓(𝑟)=𝑓(𝑟,𝜓𝑟,𝜃𝑟) as a function of spherical polar coordinates and its 3D Fourier transform is written as 𝐹(𝜔)=𝐹(𝜌,𝜓𝜔,𝜃𝜔) in frequency spherical polar coordinates. Certain relationships can be shown to hold. The relationship between the function and its transform are summarized here, with the relevant details omitted for brevity. The function itself can be expanded as a series in terms of the spherical harmonics as𝑓𝑟=𝑓𝑟,𝜓𝑟,𝜃𝑟=𝑙𝑙=0𝑘=𝑙𝑓𝑘𝑙(𝑟)𝑌𝑘𝑙𝜓𝑟,𝜃𝑟,(7.7) where 𝑌𝑘𝑙(𝜓𝑟,𝜃𝑟) are the spherical harmonics and the Fourier coefficients are given by𝑓𝑘𝑙(𝑟)=02𝜋𝜋0𝑓𝑟,𝜓𝑟,𝜃𝑟𝑌𝑘𝑙𝜓𝑟,𝜃𝑟sin𝜓𝑟𝑑𝜓𝑟𝑑𝜃𝑟.(7.8) The 3D Fourier transform of 𝑓 can also be written in Fourier space in spherical polar coordinates as𝐹𝜔=𝐹𝜌,𝜓𝜔,𝜃𝜔=𝑙𝑙=0𝑘=𝑙𝐹𝑘𝑙(𝜌)𝑌𝑘𝑙𝜓𝜔,𝜃𝜔,(7.9) where 𝐹𝑘𝑙(𝜌)=02𝜋𝜋0𝐹𝜌,𝜓𝜔,𝜃𝜔𝑌𝑘𝑙𝜓𝜔,𝜃𝜔sin𝜓𝜔𝑑𝜓𝜔𝑑𝜃𝜔.(7.10) We emphasize again that the relationship between 𝑓𝑘𝑙(𝑟) and 𝐹𝑘𝑙(𝜌) is not that of a Fourier transform. In fact, this relationship is given by [29]𝐹𝑘𝑙(𝜌)=4𝜋(𝑖)𝑙0𝑓𝑘𝑙(𝑟)𝑗𝑙(𝜌𝑟)𝑟2𝑑𝑟=4𝜋(𝑖)𝑙𝕊𝑙𝑓𝑘𝑙,(𝑟)(7.11) where 𝕊𝑙 denotes a spherical Hankel transform of order 𝑙. The inverse relationship is given by𝑓𝑘𝑙1(𝑟)=4𝜋(𝑖)𝑙2𝜋0𝐹𝑘𝑙(𝜌)𝑗𝑙(𝜌𝑟)𝜌21𝑑𝜌=4𝜋(𝑖)𝑙𝕊𝑙1𝐹𝑘𝑙(𝜌).(7.12)

8. Spatial Transfer Function Relationship in 3D

With this long set of preliminaries, definitions, and conventions established, we can now begin the derivation of the 3D spatial transfer function relationship between input 𝑠(𝑟,𝜔) and output 𝑢(𝑟,𝜔). The foregoing section makes use of several key results in [36] which will be presented here without proof. The results in [36] address general cases, not all of which are relevant to the work here. The versions of the theorems required herein are stated based on the sign conventions presented in Section 5 with regards to the sign conventions for the complex wavenumbers and also for the chosen definition of the Fourier transform. These affect the implied time dependence, the Sommerfeld radiation condition, and thus the physically correct solution that is chosen from several mathematical possibilities. The versions of the theorems presented here ensure that the chosen waves are outwardly propagating (Sommerfeld radiation condition) and also decay to zero at infinity for a damped wavefield, thus ensuring bounded and physically meaningful solutions.

Theorem 8.1. It is shown in [36] that the following result holds true: 𝐼=0𝜙(𝑥)𝑥2𝑘2𝑗𝑛(𝑥𝑟)𝑥2𝑘𝑑𝑥=𝜋𝑖2𝑛(2)(𝑘𝑟)𝜙(𝑘).(8.1) Here, 𝜙 is any analytic function defined on the positive real line that remains bounded as 𝑥 goes to infinity, 𝑗𝑛(𝑥) is a spherical Bessel function of order 𝑛, 𝑛(2)(𝑥) is a spherical Hankel function of order 𝑛, and 𝑘 is a wavenumber which may be real or complex, chosen such that the imaginary part of 𝑘 is negative. Given the definition of the Fourier transform that is being currently used, the presented result satisfies the Sommerfeld radiation condition, ensuring an outwardly propagating wave.

8.1. 3D Transfer Function Relationship in Spherical Polar Coordinates

From (4.7) and using the conversion to spherical polar coordinates in 3D, the transfer function for the Helmholtz equation in the Fourier domain can be written as 𝐺=1𝜔,𝜔𝜌2𝑘2,(8.2) which depends on the frequency spherical radial variable 𝜌 only.

In general, as the imaginary part of the wavenumber 𝑘 gets smaller, the transfer function becomes more frequency selective, in the sense of strongly passing a smaller bandwidth. As the imaginary part of the wavenumber 𝑘 gets larger, the transfer function becomes more low-pass in nature, passing lower frequencies more strongly than the higher frequencies, with a corresponding larger bandwidth of frequencies being passed. This is physically meaningful as the imaginary part of 𝑘 represents the damping inherent in this wave modality. Wavenumbers with larger imaginary parts are more heavily damped and that damping tends to affect the higher frequencies more—that is, the higher frequencies are attenuated and the lower ones are passed through, implying a low-pass nature to the physical system.

Recall that 𝜌 is the magnitude of the spatial frequency variables. The case for the imaginary part of 𝑘 equal to zero (in other words, a real wavenumber) is easy to visualize. For the case of a real 𝑘, the function is discontinuous at 𝜌=𝑘, meaning that only those spatial frequencies where exactly 𝜌=𝑘 are passed and the rest are attenuated. In other words, the resulting wave must have 𝜌2=𝜔2𝑥+𝜔2𝑦+𝜔2𝑧=𝑘2. While these statements are fairly mathematical in nature, they explain the nature of the resolution of various imaging modalities. For example, it is far more difficult to achieve the resolution of acoustic imaging with thermal imaging and one of the explanations for this is the “low-pass” blurring nature of the complex wavenumber of a thermal wave versus the highly selective band-pass nature of the acoustic wavenumber. This is discussed further in other papers, for example, in [18].

8.2. 3D Green Function Coefficients

Let us define a set of functions that will be referred to as the Green function coefficients. The actual Green function for the system is the full 3D inverse Fourier transform of the transfer function, which for a spherically symmetric function is equivalent to a spherical Hankel transform of order zero only [33]. With the help of Theorem 8.1, we define these Green function coefficients as𝑔3D𝑛(𝑟,𝑘)=𝕊𝑛11𝜌2𝑘2=2𝜋01𝜌2𝑘2𝑗𝑛(𝜌𝑟)𝜌2𝑑𝜌=𝑖𝑘𝑛(2)(𝑘𝑟),(8.3) where 𝑔3D𝑛(𝑟,𝑘) has been used to denote the 𝑛th order Green function coefficient, with the subscript indicating the order of the coefficient and the wavenumber 𝑘 included as a parameter of the function.

8.3. Transfer Function Relationship in Spherical Polar Coordinates

The 3D Fourier transform of the input (source) function is written in polar coordinates and expanded in terms of a spherical harmonic series as𝑆=𝜔,𝜔𝑙𝑙=0𝑚=𝑙𝑆𝑚𝑙(𝜌,𝜔)𝑌𝑚𝑙𝜓𝜔,𝜃𝜔.(8.4)

Hence, (4.1) for the output wavefield becomes 𝑈=𝜔,𝜔𝑙𝑙=0𝑚=𝑙𝑆𝑚𝑙(𝜌,𝜔)𝜌2𝑘2𝑌𝑚𝑙𝜓𝜔,𝜃𝜔.(8.5) The symbol 𝑘 is used to denote the wavenumber. The spatial inverse Fourier transform and thus the output in spatial coordinates is given by𝑢𝑟,𝜓𝑟,𝜃𝑟=,𝜔𝑙𝑙=0𝑚=𝑙(𝑖)𝑙24𝜋𝜋0𝑆𝑚𝑙(𝜌,𝜔)𝑗𝑙(𝜌𝑟)𝜌2𝑘2𝜌2𝑌𝑑𝜌𝑚𝑙𝜓𝑟,𝜃𝑟.(8.6) The temporal frequency 𝜔 is a constant as far as the spatial inverse Fourier transformation is concerned. We recognize that the quantity within the curly brackets can be evaluated with the help of Theorem 8.1 and using the definition of the Green function coefficients in (8.3) to obtain2𝜋0𝑆𝑚𝑙(𝜌,𝜔)𝑗𝑙(𝜌𝑟)𝜌2𝑘2𝜌2𝑑𝜌=𝑔3D𝑙(𝑟,𝑘)𝑆𝑚𝑙(𝑘,𝜔).(8.7) It now follows that the wavefield expression becomes𝑢=𝑟,𝜔𝑙𝑙=0𝑚=𝑙𝑖𝑙𝑔4𝜋3D𝑙(𝑟,𝑘)𝑆𝑚𝑙(𝑘,𝜔)𝑌𝑚𝑙𝜓𝑟,𝜃𝑟.(8.8) If the measured wavefield itself, namely, the left-hand side of (8.8), is expanded as in a spherical harmonic series so that 𝑢𝑟,𝜔=𝑢𝑟,𝜓𝑟,𝜃𝑟=,𝜔𝑙𝑙=0𝑚=𝑙𝑢𝑚𝑙(𝑟,𝜔)𝑌𝑚𝑙𝜓𝑟,𝜃𝑟,(8.9) then (8.8) gives us the simple input-output transfer function relationship we seek:𝑢𝑚𝑙𝑖(𝑟,𝜔)=𝑙𝑔4𝜋3D𝑙(𝑟,𝑘)𝑆𝑚𝑙(𝑘,𝜔).(8.10) Note how this relationship is in the spatial domain (not the frequency spatial domain) and gives a multiplicative transfer function relationship between the input coefficients 𝑆𝑚𝑙(𝑘,𝜔) and the resulting wavefield 𝑢𝑚𝑙(𝑟,𝜔).

Using the relationships proposed between the coefficients of the function in the spatial domain and the coefficients in the Fourier domain as given in (7.11), the general relationship between Hankel and Fourier transforms is given by 𝑆𝑚𝑙(𝜌,𝜔)=4𝜋(𝑖)𝑚𝕊𝑙𝑠𝑚𝑙(𝑟,𝜔)=4𝜋(𝑖)𝑙0𝑠𝑚𝑙(𝑟,𝜔)𝑗𝑙(𝜌𝑟)𝑟2𝑑𝑟.(8.11)

Hence (8.10) can be written as𝑢𝑚𝑙(𝑟,𝜔)=𝑔3D𝑙(𝑟,𝑘)0𝑠𝑚𝑙(𝑥,𝜔)𝑗𝑙(𝑘𝑥)𝑥2𝑑𝑥,(8.12) so that the value of the measured wavefield is related to the spherical Hankel transform of the input function, evaluated at the wavenumber of the wavefield in question. Using the definition of the spherical Hankel transforms, this becomes a very compact yet powerful expression:𝑢𝑚𝑙(𝑟,𝜔)=𝑔3D𝑙(𝑟,𝑘)𝕊𝑙𝑠𝑚𝑙(𝑘,𝜔).(8.13) The value of defining the Green function coefficients now becomes apparent in that it permits the relationship between 𝑢𝑚𝑙(𝑟,𝜔) and 𝑠𝑚𝑙(𝑘,𝜔) to have a simple multiplicative transfer function relationship. In other words, the relationship between output wave coefficients and input is one of multiplication with the transfer function coefficients instead of a complicated convolution-type relationship.

If we also assume that the input function is separable in time and space (a fairly general assumption) so that 𝑠(𝑟,𝜔)=𝑞(𝑟)𝜂(𝜔). This would be the case for a spatial inhomogeneity and a temporal input. In this case, then 𝑠𝑚𝑙(𝑟,𝜔)=𝑞𝑚𝑙(𝑟)𝜂(𝜔) and the formula for the forward problem is given by 𝑢𝑚𝑙(𝑟,𝜔)=𝑔3D𝑙(𝑟,𝑘)𝜂(𝜔)0𝑞𝑚𝑙(𝑟)𝑗𝑙(𝑘𝑟)𝑟2𝑑𝑟=𝑔3D𝑙(𝑟,𝑘)𝜂(𝜔)𝕊𝑙𝑞𝑚𝑙(𝑘).(8.14)

8.4. Discussion and Relationship with the Fourier Diffraction Theorem

Equation (8.13) is the key input-output relationship that is sought. It gives the relationship between input and output and the relationship is a transfer function (multiplication) type of relationship, even in the spatial domain where the normal relationship to be expected is one of convolution. This relationship is for the spherical harmonic expansion coefficients, not between the full functions themselves. However, because it is a direct, proportional type of relationship between the (𝑙,𝑚)th term of the input and the (𝑙,𝑚)th term of the output, it is particularly useful and simple to apply.

In particular, (8.13) states that the wavefield (output) coefficients are directly proportional to the transfer function coefficients and the proportionality factor between them is the spherical Hankel transform of the input coefficients, evaluated on the sphere 𝜌=𝑘. Note that (𝑙,𝑚)th order output coefficients are related to (𝑙,𝑚)th order transfer function coefficients in proportion to the 𝑙th order spherical Hankel transform of the (𝑙,𝑚)th order input coefficient. This is in fact a generalization of the Fourier diffraction theorem of tomography [3] which loosely states that the output wavefield is proportional to the Fourier transform of the input inhomogeneity evaluated somewhere on the 𝜌=𝑘 sphere in Fourier space. This is still true but we have expressed a more precise version of this theorem in the sense that we have removed any ambiguity regarding the “somewhere” on the 𝜌=𝑘 sphere and replaced it with a relationship (8.14) that does not depend on angular location. This was enabled by the use of the Fourier transform in curvilinear coordinates so that a full Fourier transform requires an (𝑙,𝑚) set of Fourier coefficients and spherical Hankel transforms. Furthermore, the idea of proportionality between Fourier transform of the inhomogeneity and output wavefield is also made more precise by stating that the proportionality between them is actually the coefficients of the Green function itself. In essence, the output wavefield can be seen as being the Green function coefficients, with each term weighted by the spherical Hankel transform of the input inhomogeneity.

9. 2D Transfer Function Relationship in 2D Polar Coordinates

We proceed to develop the transfer function expressions for the 2D case in polar coordinates. As for the 3D case, we begin with a necessary theorem which is stated without proof. As previously mentioned, the results in [36] address general cases, and the versions of these theorems required herein are stated here based on the sign conventions presented in Section 5 for the complex wavenumbers and also for the Fourier transform which affects the implied time dependence and thus the Sommerfeld radiation condition.

Theorem 9.1. It is shown in [36] that the following result holds true: 𝐼=0𝜙(𝜌)𝜌2𝑘2𝐽𝑛1(𝜌𝑟)𝜌𝑑𝜌=𝜋𝑖2𝐻𝑛(2)(𝑘𝑟)𝜙(𝑘).(9.1) Here, 𝜙 is an analytic function defined on the positive real line that remains bounded as 𝑥 goes to infinity, 𝐽𝑛(𝑥) is a Bessel function of order 𝑛, 𝐻𝑛(2)(𝑥) is a Hankel function of order 𝑛, and 𝑘 is a wavenumber which may be real or complex, chosen such that the imaginary part of 𝑘 is negative. Given the definition of the Fourier transform that is being currently used, the presented result satisfies the Sommerfeld radiation condition, ensuring an outwardly propagating wave.

9.1. 2D Green Function Coefficients

The overall system transfer function in 2D is given by𝐺=1𝜔,𝜔𝜌2𝑘2,(9.2) which depends on the frequency radial variable only.

As for the 3D case, we define the 2D Green function coefficients as the 𝑛th-order inverse Hankel transform of the overall system spatial transfer function. These will be shown to be the required transfer functions for working in the 2D polar formulation. From the definition of the inverse Hankel transform, these are given by𝑔2D𝑛(𝑟,𝑘)=𝑛11𝜌2𝑘2=01𝜌2𝑘2𝐽𝑛(𝜌𝑟)𝜌𝑑𝜌=𝜋𝑖2𝐻𝑛(2)(𝑘𝑟).(9.3) Recall that the actual (full) Green function for the system is the full 2D inverse Fourier transform of the transfer function, which is an inverse Hankel transform of order zero only. In (9.3), 𝑔2D𝑛(𝑟,𝑘) has been used to denote the nth order Green function coefficients with the subscript indicating the order.

9.2. Fourier Theorem for 2D Fourier Transforms in Polar Coordinates

The 2D Fourier transform of the input function is written in polar coordinates as 𝑆(𝜌,𝜓,𝜔)=𝑛=𝑆𝑛(𝜌,𝜔)𝑒𝑗𝑛𝜓,(9.4) where 𝑆𝑛(𝜌,𝜔) can be found from 𝑆𝑛1(𝜌,𝜔)=2𝜋02𝜋𝑆(𝜌,𝜓,𝜔)𝑒𝑗𝑛𝜓𝑑𝜓.(9.5) The output wavefield 𝑢(𝑟,𝜔) is similarly expanded as a series so that in the spatial domain we can write𝑢𝑟,𝜔=𝑢(𝑟,𝜃,𝜔)=𝑛=𝑢𝑛(𝑟)𝑒𝑗𝑛𝜃,(9.6) while in the spatial Fourier domain this is𝑈𝜔,𝜔=𝑈(𝜌,𝜓,𝜔)=𝑛=𝑈𝑛(𝜌,𝜔)𝑒𝑗𝑛𝜓.(9.7) The output wavefield is given by 𝑈(𝜔,𝜔)=𝐺(𝜔,𝜔)𝑆(𝜔,𝜔) and since 𝐺(𝜔,𝜔) is radially symmetric and does not require a full series, the output wavefield in the Fourier domain is given by𝑈(𝜌,𝜓,𝜔)=𝑛=𝑈𝑛(𝜌,𝜔)𝑒𝑗𝑛𝜓=𝑛=𝑆𝑛(𝜌,𝜔)𝜌2𝑘2𝑒𝑗𝑛𝜓,(9.8) or equivalently as𝑈𝑛𝑆(𝜌,𝜔)=𝑛(𝜌,𝜔)𝜌2𝑘2.(9.9) The equivalent expression to (9.8) in the spatial domain is given by𝑢=𝑟,𝜔𝑛=𝑢𝑛(𝑟)𝑒𝑗𝑛𝜃=𝑛=𝑖𝑛2𝜋0𝑆𝑛(𝜌,𝜔)𝜌2𝑘2𝐽𝑛𝑒(𝜌𝑟)𝜌𝑑𝜌𝑗𝑛𝜃,(9.10) or𝑢𝑛𝑖(𝑟)=𝑛2𝜋0𝑆𝑛(𝜌,𝜔)𝜌2𝑘2𝐽𝑛(𝜌𝑟)𝜌𝑑𝜌.(9.11) The temporal frequency 𝜔 is a constant as far as the spatial inverse Fourier transformation is concerned. The quantity within the curly brackets can be evaluated with the help of Theorem 9.1 as0𝑆𝑛(𝜌,𝜔)𝜌2𝑘2𝐽𝑛(𝜌𝑟)𝜌𝑑𝜌=𝜋𝑖2𝐻𝑛(2)(𝑘𝑟)𝑆𝑛(𝑘,𝜔)=𝑔2D𝑛(𝑟,𝑘)𝑆𝑛(𝑘,𝜔).(9.12) Hence, the output wavefield expression is given by𝑢=𝑟,𝜔𝑛=𝑢𝑛(𝑟)𝑒𝑗𝑛𝜃=𝑛=𝑖𝑛𝑔2𝜋2D𝑛(𝑟,𝑘)𝑆𝑛(𝑘,𝜔)𝑒𝑗𝑛𝜃.(9.13) The definition of forward and inverse Hankel transforms, the general relationship between Hankel and Fourier transforms is given by 𝑆𝑛(𝜌)=2𝜋𝑖𝑛𝑛𝑠𝑛(𝑟)=2𝜋𝑖𝑛0𝑠𝑛(𝑟)𝐽𝑛(𝜌𝑟)𝑟𝑑𝑟,(9.14) and may be used to simplify (9.13) to yield𝑢=𝑟,𝜔𝑛=𝑢𝑛(𝑟)𝑒𝑗𝑛𝜃=𝑛=𝑔2D𝑛(𝑟,𝑘)0𝑠𝑛(𝑥,𝜔)𝐽𝑛(𝑘𝑥)𝑥𝑑𝑥𝑒𝑗𝑛𝜃,(9.15) or more compactly as 𝑢𝑛(𝑟,𝜔)=𝑔2D𝑛(𝑟,𝑘)0𝑠𝑛(𝑥,𝜔)𝐽𝑛(𝑘𝑥)𝑥𝑑𝑥.(9.16) It is in (9.15) that the interpretation of the evaluation of the 2D Fourier transform becomes apparent, namely through the evaluation of the Bessel function at the (real or complex) wavenumbers. Using the definition of the Hankel transform, (9.16) can be written in the compact and powerful formulation of 𝑢𝑛(𝑟,𝜔)=𝑔2D𝑛(𝑟,𝑘)𝑛𝑠𝑛(𝑘,𝜔).(9.17) If it is also further assumed that the inhomogeneity function is separable in time and space (a fairly general assumption) so that 𝑠(𝑟,𝜔)=𝑞(𝑟)𝜂(𝜔), and𝑛=𝑠𝑛(𝑟,𝜔)𝑒𝑗𝑛𝜃=𝜂(𝜔)𝑛=𝑞𝑛(𝑟)𝑒𝑗𝑛𝜃,(9.18) then (9.16) can be further reduced to 𝑢𝑛(𝑟,𝜔)=𝑔2D𝑛(𝑟,𝑘)𝜂(𝜔)𝑛𝑞𝑛(𝑘).(9.19) It is noted that the value of the measured wavefield at any position 𝑟 is related to the Hankel transform of the heterogeneity function evaluated at the wavenumber of the wavefield in question, with the Green function coefficient acting as the proportionality term.

Equation (9.19) is the 2D equivalent of (8.14). The same comments can be made as those made in the discussion after the 3D version of this problem, as the relationship is identical, save for translating 3D tools such as the spherical harmonic expansions and spherical Hankel transforms into 2D tools involving Fourier series and Hankel transforms.

10. Transfer Function Relationship in 1D

We proceed to develop the equivalent relationship for the 1D case. We have seen from (4.7) that the 1D transfer function in the Fourier domain is given by 𝐺𝜔𝑥=1,𝜔𝜔2𝑥𝑘2,(10.1) which depends on only a single spatial frequency. Equation (4.2) in 1D then becomes1𝑢(𝑟,𝜔)=(2𝜋)𝑆𝜔𝑥,𝜔𝜔2𝑥𝑘2𝑒𝑖𝜔𝑥𝑟𝑑𝜔𝑥.(10.2) The temporal frequency 𝜔 is a constant as far as the spatial inverse Fourier transformation is concerned. We are thus interested in calculating integrals of the form 1𝐼=2𝜋0𝜙(𝑥)𝑥2𝑘2𝑒𝑖𝑥𝑟𝑑𝑥,(10.3) where 𝜙 is an analytic function defined on the positive real line that approaches zero as 𝑥 goes to infinity. The required theorem is given in [36] and the relevant result is presented below.

Theorem 10.1. It is shown in [36] that the following is true 12𝜋𝜙(𝜌)𝜌2𝑘2𝑒𝑖𝜌𝑟1𝑑𝜌=2𝑖𝑘𝜙(𝑘)𝑒𝑖𝑘𝑟1,𝑟>0,2𝑖𝑘𝜙(𝑘)𝑒𝑖𝑘𝑟,𝑟<0.(10.4) Here, 𝜙 is an analytic function defined on the positive real line that remains bounded as 𝑥 goes to infinity and 𝑘 is a wavenumber which may be real or complex, chosen such that the imaginary part of 𝑘 is negative. Given the definition of the Fourier transform that is being currently used, the presented result satisfies the Sommerfeld radiation condition, ensuring an outwardly propagating wave.

10.1. 1D Green's Functions via Inverse Fourier Transformation of the Spatial Transfer Functions

The Green's function for the Helmholtz equation is 1D and is given by inverse Fourier transform of the spatial transfer function as𝑔1D1(𝑟,𝑘)=2𝜋𝑒𝑖𝜔𝑥𝑟𝜔2𝑥𝑘2𝑑𝜔𝑥,(10.5) and can be evaluated with the help of Theorem 10.1 to give 𝑔1D(1𝑟,𝑘)=𝑒2𝑖𝑘𝑖𝑘𝑟1𝑟>0𝑒2𝑖𝑘𝑖𝑘𝑟=1𝑟<0𝑒2𝑖𝑘𝑖𝑘|𝑟|.(10.6)

10.2. Transfer Function Relationship in 1D

Equation (10.2) can now be evaluated with the help of Theorem 10.1 as well as the result in (10.6) as 1𝑢(𝑟,𝜔)=(2𝜋)𝑆𝜔𝑥,𝜔𝜔2𝑥𝑘2𝑒𝑖𝜔𝑥𝑟𝑑𝜔𝑥=12𝑖𝑘𝑆(𝑘,𝜔)𝑒𝑖𝑘𝑟1𝑟>02𝑖𝑘𝑆(𝑘,𝜔)𝑒𝑖𝑘𝑟𝑟<0=𝑔1D(𝑟,𝑘)𝑆(sgn(𝑟)𝑘,𝜔),(10.7) where 𝑔1𝐷(𝑟,𝑘)=(1/2𝑖𝑘)𝑒𝑖𝑘|𝑟|. This is the 1D version of the transfer function relationship, with the wavefield being directly proportional to the Fourier transform of the object function evaluated at 𝑘. This is similar to the results for the 2D and 3D cases. For a real wavenumber, a “sphere” in 1D becomes the two points on the real line at ±𝑘 and the proportionality term is the Green's function for the space in question.

As before, we assume that the inhomogeneity function is separable in time and space (a fairly general assumption) so that 𝑠(𝑟,𝜔)=𝑞(𝑟)𝜂(𝜔)𝑆(𝜔𝑥,𝜔)=𝑄(𝜔𝑥)𝜂(𝜔), where 𝑄(𝜔𝑥) is the 1D Fourier transform of 𝑞(𝑟). The relevant result now reads𝑢(𝑟,𝜔)=𝑔1D(𝑟,𝑘)𝜂(𝜔)𝑄(sgn(𝑟)𝑘).(10.8) This result is directly applicable for computational purposes.

11. Summary of Results

The relationships in the previous sections are summarized in Table 1.

12. Applications of the Transfer Function Relationships to the Wave Equation

The transfer function relationships can be used to find time domain Green's functions where the temporal Fourier integral can be easily inverted. In particular, we consider the case where the input to the standard wave equation is a Dirac-delta function at the origin2𝑢1𝑟,𝑡𝑐2𝜕2𝜕𝑡2𝑢𝑟,𝑡=𝛿𝑟𝛿(𝑡).(12.1)

This corresponds to the Helmholtz equation above with 𝑘=𝜔/𝑐. The temporal and spatial Fourier transform of 𝛿(𝑟)𝛿(𝑡) is 1, which is spherically symmetric so that the transfer function relationships given above only need the zeroth-order component and simplify considerably. In 3D, the transfer function relationship is given by𝑢(𝑟,𝜔)=𝑔3D0(𝑟,𝑘)1=𝑖𝑘0(2)(𝑘𝑟).(12.2) In 2D, the relationship is given by𝑢(𝑟,𝜔)=𝑔2D0(𝑟,𝑘)1=𝜋𝑖2𝐻0(𝑘𝑟).(12.3) In 1D, the relationship is𝑢(𝑟,𝜔)=𝑔1D1(𝑟,𝑘)1=𝑒2𝑖𝑘𝑖𝑘|𝑟|.(12.4) These three relationships can now be inverse Fourier transformed in time.

For the 3D case, we use the fact that 0(2)(𝑘𝑟)=𝑗0(𝑘𝑟)𝑖𝑦0(𝑘𝑟)=𝑖exp(𝑖𝑘𝑟)/𝑘𝑟 so that the inverse Fourier transform of (12.2) gives1𝑢(𝑟,𝑡)=2𝜋1𝑟𝑒𝑖𝜔𝑟/𝑐𝑒𝑖𝜔𝑡=1𝑑𝜔𝑟𝛿𝑟𝑡𝑐.(12.5) The 1D time response can be found from the inverse Fourier transform of  (12.4)1𝑢(𝑟,𝑡)=2𝜋𝑐𝑒2𝑖𝜔𝑖𝜔|𝑟|/𝑐𝑒𝑖𝜔𝑡=𝑐𝑑𝜔42𝐻𝑡|𝑟|𝑐=𝑐14sgn𝑡|𝑟|𝑐,(12.6) where 𝐻(𝑥) is the Heaviside unit step function.

The 2D case is a little more complicated but can nevertheless be evaluated in closed form by an inverse Fourier transform of (12.3). To do this, we write the zeroth-order Hankel function in its integral form as [37]𝐻0(2)(𝑥)=𝐽0(𝑥)𝑖𝑌0(𝑥)=2𝑖𝜋1𝑒𝑖𝑥𝜏𝜏21𝑑𝜏.(12.7) Hence from (12.7) and (12.3), the inverse Fourier transform of (12.3) gives1𝑢(𝑟,𝑡)=2𝜋1𝑒𝑖𝜔𝑟𝜏/𝑐𝜏21𝑑𝜏𝑒𝑖𝜔𝑡𝑑𝜔.(12.8) Changing the order of integration gives𝑢(𝑟,𝑡)=11𝜏2112𝜋𝑒𝑖𝜔[𝑡𝑟𝜏/𝑐]𝑑𝜔𝑑𝜏.(12.9) But 12𝜋𝑒𝑖𝜔[𝑡𝑟𝜏/𝑐]𝑑𝜔=𝛿𝑡𝑟𝜏𝑐,(12.10) so that (12.9) becomes𝑢(𝑟,𝑡)=11𝜏2𝛿1𝑡𝑟𝜏𝑐𝑑𝜏=1𝜏21𝐻(𝜏1)𝛿𝑡𝑟𝜏𝑐𝑑𝜏.(12.11) Changing variables so that 𝑥=𝑟𝜏/𝑐 gives𝑐𝑢(𝑟,𝑡)=𝑟𝐻𝑐𝑥𝑟11(𝑐𝑥/𝑟)2𝑟1𝛿(𝑡𝑥)𝑑𝑥=𝐻𝑡𝑐1𝑡2𝑟2/𝑐2,(12.12) which yields the desired result in closed form.

13. Application: Inversion Formulas for Tomographic Applications with Real Wavenumbers

In the case of a real wavenumber 𝑘, an immediate application of the above formulas is for closed-form inversion formulas which are useful for tomographic applications. These inversion formulas are applicable to applications where the wavefield (e.g., acoustic field) is measured and the goal is to reconstruct the source (input) that led to that measured wavefield (output). As we have seen, formulas above give the output (wavefield) in terms of the transform of the input evaluated at 𝑘. Thus, if the output wavefield is measured, this implies knowledge of the transform of the input. An inverse transform then leads to the input itself.

13.1. Inversion Formula in 3D

Equation (8.14) admits an inversion useful for tomographic applications. It can be written as 𝑖𝑢𝑘𝑙𝑟0,𝜔𝑘𝑠𝑙(2)𝑘𝑠𝑟0=𝜂(𝜔)0𝜙𝑘𝑙(𝑥)𝑗𝑙𝑘𝑠𝑥𝑥2𝑑𝑥,(13.1) where 𝑥 has been used as a dummy integration variable in order to avoid possible confusion and 𝑟0 has been used as the radial variable and indicates the position where a measurement of the wavefield is made. Multiplying both sides by 𝑗𝑙(𝑘𝑠𝑟)𝑘2𝑠 and integrating over 𝑘𝑠 gives0𝑖𝑢𝑘𝑙𝑟0,𝜔𝑘𝑠𝑙(2)𝑘𝑠𝑟0𝑗𝜂(𝜔)𝑙𝑘𝑠𝑟𝑘2𝑠𝑑𝑘𝑠=00𝜙𝑘𝑙(𝑥)𝑗𝑙𝑘𝑠𝑥𝑥2𝑑𝑥𝑗𝑙𝑘𝑠𝑟𝑘2𝑠𝑑𝑘𝑠=00𝜙𝑘𝑙(𝑥)𝑗𝑙𝑘𝑠𝑥𝑗𝑙𝑘𝑠𝑟𝑘2𝑠𝑑𝑘𝑠𝑥2𝑑𝑥.(13.2)

Using the orthogonality of the spherical Bessel functions, it follows that (13.2) becomes 0𝑖𝑢𝑘𝑙𝑟0,𝜔𝑘𝑠𝑙(2)𝑘𝑠𝑟0𝑗𝜂(𝜔)𝑙𝑘𝑠𝑟𝑘2𝑠𝑑𝑘𝑠=0𝜙𝑘𝑙𝜋(𝑥)2𝑥2𝛿(𝑥𝑟)𝑥2𝜋𝑑𝑥=2𝜙𝑘𝑙(𝑟).(13.3) Recalling that 𝑘𝑠=𝜔/𝑐𝑠, the inversion formula for the spatial source becomes 𝜙𝑘𝑙(𝑟)=2𝑖𝜋𝑐2𝑠0𝑢𝑘𝑙𝑟0,𝜔𝑙(2)𝜔𝑟0/𝑐𝑠𝑗𝜂(𝜔)𝑙𝜔𝑟𝑐𝑠𝜔𝑑𝜔.(13.4) If the measured wavefield is also separable in space and time so that 𝑢𝑘𝑙(𝑟0,𝜔)=𝜒𝑘𝑙(𝑟0)𝑇(𝜔), then (13.4) can be further simplified to 𝜙𝑘𝑙(𝑟)=2𝑖𝜒𝑘𝑙𝑟0𝜋𝑐2𝑠0𝑇(𝜔)𝑙(2)𝜔𝑟0/𝑐𝑠𝑗𝜂(𝜔)𝑙𝜔𝑟𝑐𝑠𝜔𝑑𝜔.(13.5) Equation (13.5) gives a simple form for finding the spatial source function at any position 𝑟 once having measured the temporal response of the wavefield at a position 𝑟0. However, since the spherical harmonic transform 𝑢𝑘𝑙(𝑟0,𝜔) of the wavefield is required, this implies that sufficient angular information about the wavefield 𝑢 must be obtained in order to permit this spherical harmonic calculation.

13.2. Inversion Formula in 2D

Equation (9.19) admits an inversion leading to the input function. We assume that the source function is separable in time and space (a fairly general assumption) so that 𝑠(𝑟,𝜔)=𝜙(𝑟)𝜂(𝜔). In this case, then 𝑠𝑘𝑙(𝑟,𝜔)=𝜙𝑘𝑙(𝑟)𝜂(𝜔) and (8.12) becomes 2𝑖𝑢𝑛𝑟0,𝜔𝜋𝐻𝑛(2)𝑘𝑠𝑟0=𝜂(𝜔)0𝜙𝑛(𝑥)𝐽𝑛𝑘𝑠𝑥𝑥𝑑𝑥,(13.6) where 𝑥 has been used as a dummy integration variable in order to avoid possible confusion and 𝑟0 indicates the radial position where the wavefield is measured. Multiplying both sides of (13.6) by 𝐽𝑛(𝑘𝑠𝑟)𝑘𝑠 and integrating over 𝑘𝑠 gives02𝑖𝑢𝑛𝑟0,𝜔𝜋𝐻𝑛(2)𝑘𝑠𝑟0𝐽𝜂(𝜔)𝑛𝑘𝑠𝑟𝑘𝑠𝑑𝑘𝑠=00𝜙𝑛(𝑥)𝐽𝑛𝑘𝑠𝑥𝑥𝑑𝑥𝐽𝑛𝑘𝑠𝑟𝑘𝑠𝑑𝑘𝑠,=0𝜙𝑛(𝑥)0𝐽𝑛𝑘𝑠𝑥𝐽𝑛𝑘𝑠𝑟𝑘𝑠𝑑𝑘𝑠𝑥𝑑𝑥.(13.7) Using the orthogonality of the spherical Bessel functions, it follows that (13.7) becomes 02𝑖𝑢𝑛𝑟0,𝜔𝜋𝐻𝑛(2)𝑘𝑠𝑟0𝐽𝜂(𝜔)𝑛𝑘𝑠𝑟𝑘𝑠𝑑𝑘𝑠=0𝜙𝑛1(𝑥)𝑥𝛿(𝑥𝑟)𝑥𝑑𝑥=𝜙𝑛(𝑟).(13.8) Recalling that 𝑘𝑠=𝜔/𝑐𝑠, the inversion formula for the inhomogeneity becomes 𝜙𝑛(𝑟)=2𝑖𝜋𝑐2𝑠0𝑢𝑛𝑟0,𝜔𝐻𝑛(2)𝜔𝑟0/𝑐𝑠𝐽𝜂(𝜔)𝑛𝜔𝑟𝑐𝑠𝜔𝑑𝜔.(13.9) If the measured wavefield is also separable in space and time so that 𝑢𝑛(𝑟0,𝜔)=𝜒𝑛(𝑟0)𝑇(𝜔), then (13.4) can be further simplified to 𝜙𝑛(𝑟)=2𝑖𝜒𝑛𝑟0𝜋𝑐2𝑠0𝑇(𝜔)𝐻𝑛(2)𝜔𝑟0/𝑐𝑠𝐽𝜂(𝜔)𝑛𝜔𝑟𝑐𝑠𝜔𝑑𝜔.(13.10) Equation (13.10) gives a simple form for finding the spatial source function at any position 𝑟 once having measured the wavefield response at a radial position 𝑟0. Comparing (13.10) and (13.5), it is noted that they have identical forms with the exception of the replacement of Hankel and Bessel functions for the 2D case with spherical Hankel and Bessel functions for the 3D case.

13.3. Inversion Formula in 1D

We can also further assume that the wavefield is also separable in time and space so that 𝑢(𝑟,𝜔)=𝜒(𝑟)𝑇(𝜔), which leads to to𝜙𝑘𝑠=𝜒(𝑟)𝑇(𝜔)𝜂(𝜔)2𝑖𝑘𝑠𝑒𝑖𝑘𝑠𝑟=𝜙(𝑦)𝑒𝑖𝑦𝑘𝑠𝜙𝑘𝑑𝑦𝑟>0,𝑠=𝜒(𝑟)𝑇(𝜔)𝜂(𝜔)2𝑖𝑘𝑠𝑒𝑖𝑘𝑠𝑟=𝜙(𝑦)𝑒𝑖𝑦𝑘𝑠𝑑𝑦𝑟<0.(13.11) Using the orthogonality of the Fourier kernel, both sides are multiplied by 𝑒±𝑖𝑥𝑘𝑠 and integrated over all 𝑘𝑠 to give𝜒(𝑟)𝑇(𝜔)𝜂(𝜔)2𝑖𝑘𝑠𝑒𝑖𝑘𝑠𝑟𝑒𝑖𝑥𝑘𝑠𝑑𝑘𝑠=𝜙(𝑦)𝑒𝑖𝑦𝑘𝑠𝑒𝑖𝑥𝑘𝑠𝑑𝑘𝑠𝑑𝑦=2𝜋𝜙(𝑦)𝛿(𝑦𝑥)𝑑𝑦𝑟>0𝜙(𝑥)=𝑖𝜒(𝑟)𝜋𝑐2𝑠𝑇(𝜔)𝑒𝜂(𝜔)𝑖𝜔((𝑟𝑥)/𝑐𝑠)𝜔𝑑𝜔𝑟>0,(13.12) and similarly𝜒(𝑟)𝑇(𝜔)𝜂(𝜔)2𝑖𝑘𝑠𝑒𝑖𝑘𝑠𝑟𝑒𝑖𝑥𝑘𝑠𝑑𝑘𝑠=𝜙(𝑦)𝑒𝑖𝑦𝑘𝑠𝑒𝑖𝑥𝑘𝑠𝑑𝑘𝑠𝑑𝑦=2𝜋𝜙(𝑦)𝛿(𝑦𝑥)𝑑𝑦𝑟<0𝜙(𝑥)=𝑖𝜒(𝑟)𝜋𝑐2𝑠𝑇(𝜔)𝑒𝜂(𝜔)𝑖𝜔((𝑥𝑟)/𝑐𝑠)𝜔𝑑𝜔𝑟<0.(13.13) Several points need to be made regarding (13.12) and (13.13). First, they both give the source function as a function of position 𝑥 along the real line. The interpretation of the variable 𝑟 is that of the position at which the measurement is made and is considered to be a fixed quantity. Both (13.12) and (13.13) are inverse Fourier transforms of (𝑇(𝜔)/𝜂(𝜔))𝜔, evaluated at (𝑥𝑟)/𝑐𝑠 or (𝑟𝑥)/𝑐𝑠, depending on whether measurements are made in transmission or reflection. Clearly |𝑥𝑟|/𝑐𝑠 is the time taken for a wave to travel the distance |𝑥𝑟|.

The 2D and 3D equivalent to (13.12) and (13.13) are (13.5) and (13.10) which also similarly involve an inverse transform of (𝑇(𝜔)/𝜂(𝜔))𝜔 but have a different kernel for the integration. The reason for the differences in the nature of the kernels from the 2D/3D cases to the 1D case is that in the 2D and 3D cases the kernels used for the Fourier transforms are the Bessel and spherical Bessel functions, respectively, which are the standing wave solutions in those dimensions. However, Green's functions in those dimensions are the Hankel and spherical Hankel functions which are the travelling wave solutions. We note that the kernels for the inversion formulas of (13.5) and (13.10) involve inverse Hankel and spherical Hankel functions, which represent the inverse of Green's function. In contrast, the travelling wave solutions in 1D are the complex exponential and those are also the kernel of 1D Fourier transform. The ratio of the two complex exponentials (one represents the Fourier kernel, the other represents the inverse of 1D Green's function) finally give another complex exponential which represents the shift |𝑥𝑟|/𝑐𝑠.

14. Summary and Conclusions

This paper presents the derivation of spatial transfer function relationships for the Helmholtz equation. The focus has been on deriving forward transfer function relationships so that once given a particular input function, the resulting output wavefield can be calculated at any point in space. These are termed “transfer function” relationships because the relationship between the input and output quantities is one of multiplication and no convolution is involved, even in the spatial domain where normally a convolution would be required. Interestingly, instead of Green's function itself, Green's function coefficients for the space are required. These Green function coefficients form the kernel of the output response. Each element of the kernel is then weighted (multiplicatively) by the relevant Fourier/Hankel/Spherical Hankel transform of the input, evaluated on the 𝑛-dimensional sphere of radius 𝑘, where 𝑘 is the wavenumber and 𝑛 is the dimension of the space. Combined together, these form the final output response. This is true for dimensions 𝑛=1, 2, 3. These simple but powerful results serve to cast the entire problem as an input-output problem with a transfer function relating input to output. In this view of the problem, the input (inhomogeneity) and output (resulting wavefield) in space are related by a simple multiplicative transfer function relationship and not via a convolution. Some applications of these results were shown in the manuscript.