# Multidimensional Wave Field Signal Theory: Transfer Function Relationships

Authors: Natalie Baddour
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 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 [1–4]. 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 [5–8]. Similarly, photothermal tomographic imaging methods have also been used for nondestructive evaluation [9–12], 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 [16–18] 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 [21–27] 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 thespatial (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. Sections2 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 inn-Dimensions

The theory of Fourier transforms can be extended toℝn in a way that is completely analogous to the treatment of one-dimensional Fourier transforms. Let elements of ℝn be denoted by x⃗=(x1,x2,…,xn) or more generally r⃗, and elements of the corresponding Fourier dual space Ωn be denoted by ω⃗=(ω1,ω2,…,ωn). Under the suitability of integration of f, then the Fourier transform in ℝn is defined for ω⃗in Ωn as(2.1)F(ω⃗)=∫Rnf(x⃗)e-jω⃗⋅x⃗dx.⃗
Under suitable conditions, the function f can be recovered from the inverse transform through
(2.2)f(x⃗)=1(2π)n∫ΩnF(ω⃗)ejω⃗⋅x⃗dω⃗.
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:(3.1)∇2u(r⃗,ω)+ks2u(r⃗,ω)=-s(r⃗,ω),
where ks2=ω2/cs2. Here, s(r⃗,ω) 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 u(r⃗,ω) 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 s(r⃗,ω) and u(r⃗,ω) are functions of position, r⃗, and (temporal) frequency, ω. The variable cs 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, ks2 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 densityu(r⃗,t) in a solid due to incident energy intensity s(r⃗,t) (optical source function) is given in the time domain by [2]
(3.2)D∇2u(r⃗,t)-1c∂∂tu(r⃗,t)-μau(r⃗,t)=-s(r⃗,t).In the above equation,μa is the optical absorption coefficient (m−1), c 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:(3.3)∇2u(r⃗,ω)+kp2u(r⃗,ω)=-s(r⃗,ω)D,
where kp2=-(μa/D)-(iω/cD) 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 kp2 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 (3.4)∇2u(r⃗,t)-1α∂∂tu(r⃗,t)=-s(r⃗,t),
where s(r⃗,t) is the time- and space-dependent heat source, u(r⃗,t) 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 obtain(3.5)∇2u(r⃗,ω)+kt2u(r⃗,ω)=-s(r⃗,ω),
where kt2=-iω/α 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 k2 indicating a propagating wave and a complex k2 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(4.1)U(ω⃗,ω)=U(ωx,ωy,ωz,ω)=1(ωx2+ωy2+ωz2-k2)S(ω⃗,ω),
where we have used the shorthand notation (ω⃗)=(ωx,ωy,ωz) 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 ωk2=ωx2+ωy2+ωz2, so that ωk 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 ωz spatial frequency variable would not be required and so forth.By inverse spatial Fourier transformation of (4.1), the wavefield in n dimensions is given by(4.2)u(r⃗,ω)=1(2π)n∫-∞∞S(ω⃗,ω)(ωk2-k2)eiω⃗⋅r⃗dω⃗.
Using the definition of the Fourier transform of s(r⃗,ω), the above equation can be rewritten as(4.3)u(r⃗,ω)=1(2π)n∫-∞∞eiω⃗⋅r⃗(ωk2-k2)∫-∞∞s(x⃗,ω)e-iω⃗⋅x⃗dx⃗dω⃗.
Let us define the spatial Green’s functionin n-dimensional space as(4.4)g(r⃗∣x⃗,ω)=g(r⃗-x⃗,ω)=1(2π)n∫-∞∞eiω⃗⋅(r⃗-x⃗)ωk2-k2dω⃗.
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)u(r⃗,ω)=∫-∞∞s(x⃗,ω)g(r⃗-x⃗,ω)dx⃗=s(r⃗,ω)*r⃗g(r⃗,ω).
The notation *r⃗ 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)U(ω⃗,ω)=G(ω⃗,ω)S(ω⃗,ω),
where the Fourier transform of the spatial Green's function is the spatial transfer function:(4.7)G(ω⃗,ω)=F[g(r⃗,ω)]=1(ωk2-k2).
The dependence of the spatial Green's function on temporal frequency is via the wavenumber k of the wavefield (which may be real or complex) and this is emphasized by writing G(ω⃗,ω) 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 n-dimensional space with Green's function of (4.4) being the n-dimensional spatial inverse Fourier transform of the spatial transfer function of (4.7). More specifically, in 1D, ωk2=ωx2, and in 2D ωk2=ωx2+ωy2. Transforming to polar coordinates in the 2D spatial frequency domain via the transformation ωx=ρcosψ, ωy=ρsinψ gives ωk2=ωx2+ωy2=ρ2 where ρ is the radial frequency variable in polar coordinates. Similarly, in 3D, ωk2=ωx2+ωy2+ωz2 becomes ωk2=ωx2+ωy2+ωz2=ρ2 where ρ is the spherical radial frequency variable in spherical polar coordinates in spatial frequency space via the coordinate transformation given by ωx=ρsinψωcosθω, ωy=ρsinψωsinθω and ωz=ρ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.


## 4.1. Spatial Transfer Function in Curvilinear Coordinates

Equations (4.4) and (4.7) hold in any n-dimensional space with Green's function of (4.4) being the n-dimensional spatial inverse Fourier transform of the spatial transfer function of (4.7). More specifically, in 1D, ωk2=ωx2, and in 2D ωk2=ωx2+ωy2. Transforming to polar coordinates in the 2D spatial frequency domain via the transformation ωx=ρcosψ, ωy=ρsinψ gives ωk2=ωx2+ωy2=ρ2 where ρ is the radial frequency variable in polar coordinates. Similarly, in 3D, ωk2=ωx2+ωy2+ωz2 becomes ωk2=ωx2+ωy2+ωz2=ρ2 where ρ is the spherical radial frequency variable in spherical polar coordinates in spatial frequency space via the coordinate transformation given by ωx=ρsinψωcosθω, ωy=ρsinψωsinθω and ωz=ρ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 kt represents a generic complex wavenumber while ks represents a generic real wavenumber.Rather than the squared wavenumber, the quantity of interest will prove to be the wavenumber itself, namely,ks and kt, which are the square roots of the given squared wavenumber in the Helmholtz equation. Each k can be considered as the sum of a real and an imaginary part, so that kt=ktr+ikti with ktr denoting the real part of kt and kti denoting the imaginary part. Since the square root of any k2 can be ±k, we will use the convention that for a given complexk, the required square root of the corresponding k2 is defined such that the imaginary part ofkis negative. Hence, kt is chosen as the square root of kt2 such that kti<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 solutionu(x), where x is the spatial variable, to the Helmholtz equation is considered to be radiating if it satisfies(5.1)lim⁡|x|→∞|x|(n-1)/2(∂∂|x|+ik)u(x)=0,
where n is the dimension of the space and k is the wavenumber in the Helmholtz equation. This is the radiation equation based on an implied time variation ofeiωt 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 f′(t) is iωF(ω), clearly implying the eiωt dependence. Had a different definition of the Fourier transform been used, the implied time variation would have been e-iωt, in which case the sign of i in (5.1) would be reversed.


## 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 solutionu(x), where x is the spatial variable, to the Helmholtz equation is considered to be radiating if it satisfies(5.1)lim⁡|x|→∞|x|(n-1)/2(∂∂|x|+ik)u(x)=0,
where n is the dimension of the space and k is the wavenumber in the Helmholtz equation. This is the radiation equation based on an implied time variation ofeiωt 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 f′(t) is iωF(ω), clearly implying the eiωt dependence. Had a different definition of the Fourier transform been used, the implied time variation would have been e-iωt, in which case the sign of i 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 ordern is defined by the integral [31](6.1)F⌢n(ρ)=Hn(f(r))=∫0∞f(r)Jn(ρr)rdr,
where Jn(z) is the nth order Bessel function. If n>-1/2, the transform is self-reciprocating and the inversion formula is given by (6.2)f(r)=Hn-1{F⌢n(ρ)}=∫0∞F⌢n(ρ)Jn(ρr)ρdρ.


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

To develop Fourier transforms in 2D in polar coordinates, both the functionf(r⃗) and its Fourier transform F(ω⃗) are expressed in polar coordinates. In general f(r⃗)=f(r,θ) is not radially symmetric and is a function of both r and θ so that the θ dependence can be expanded into a Fourier series due to the 2π periodicity of the function in θ:(6.3)f(r,θ)=∑n=-∞∞fn(r)ejnθ,
where the Fourier coefficients fn(r) can be found from(6.4)fn(r)=12π∫02πf(r,θ)e-jnθdθ.
Similarly, the 2D Fourier transform F(ω⃗)=F(ρ,ψ) can also be expanded into its own Fourier series so that(6.5)F(ρ,ψ)=∑n=-∞∞Fn(ρ)ejnψ,
and where those Fourier coefficients are similarly found from(6.6)Fn(ρ)=12π∫02πF(ρ,ψ)e-jnψdψ,The relationship between the Fourier coordinates in normal spacefn(r) in (6.3) and the spatial frequency space Fourier coordinates Fn(ρ) 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 fn(r) and Fn(ρ) is not a Fourier transform. The relationship between them is given by(6.7)Fn(ρ)=2πi-n∫0∞fn(r)Jn(ρr)rdr=2πi-nHn{fn(r)},
where ℍn is the nth order Hankel transform. The reverse relationship is given by(6.8)fn(r)=in2π∫0∞Fn(ρ)Jn(ρr)ρdρ=in2πHn{Fn(ρ)}.
The nth term in the Fourier series for the original function will Hankel transform into the nth term of the Fourier series of the Fourier transform. However, it is an nth order Hankel transform for the nth term, namely, all the terms are not equivalently transformed. The mapping from fn(r) to Fn(ρ) is one of nth-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 thenth-order Hankel transform (of the radial variable to the spatial radial variable) of the nth 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.


## 6.1. Hankel Transforms

The Hankel transform of ordern is defined by the integral [31](6.1)F⌢n(ρ)=Hn(f(r))=∫0∞f(r)Jn(ρr)rdr,
where Jn(z) is the nth order Bessel function. If n>-1/2, the transform is self-reciprocating and the inversion formula is given by (6.2)f(r)=Hn-1{F⌢n(ρ)}=∫0∞F⌢n(ρ)Jn(ρr)ρdρ.


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

To develop Fourier transforms in 2D in polar coordinates, both the functionf(r⃗) and its Fourier transform F(ω⃗) are expressed in polar coordinates. In general f(r⃗)=f(r,θ) is not radially symmetric and is a function of both r and θ so that the θ dependence can be expanded into a Fourier series due to the 2π periodicity of the function in θ:(6.3)f(r,θ)=∑n=-∞∞fn(r)ejnθ,
where the Fourier coefficients fn(r) can be found from(6.4)fn(r)=12π∫02πf(r,θ)e-jnθdθ.
Similarly, the 2D Fourier transform F(ω⃗)=F(ρ,ψ) can also be expanded into its own Fourier series so that(6.5)F(ρ,ψ)=∑n=-∞∞Fn(ρ)ejnψ,
and where those Fourier coefficients are similarly found from(6.6)Fn(ρ)=12π∫02πF(ρ,ψ)e-jnψdψ,The relationship between the Fourier coordinates in normal spacefn(r) in (6.3) and the spatial frequency space Fourier coordinates Fn(ρ) 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 fn(r) and Fn(ρ) is not a Fourier transform. The relationship between them is given by(6.7)Fn(ρ)=2πi-n∫0∞fn(r)Jn(ρr)rdr=2πi-nHn{fn(r)},
where ℍn is the nth order Hankel transform. The reverse relationship is given by(6.8)fn(r)=in2π∫0∞Fn(ρ)Jn(ρr)ρdρ=in2πHn{Fn(ρ)}.
The nth term in the Fourier series for the original function will Hankel transform into the nth term of the Fourier series of the Fourier transform. However, it is an nth order Hankel transform for the nth term, namely, all the terms are not equivalently transformed. The mapping from fn(r) to Fn(ρ) is one of nth-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 thenth-order Hankel transform (of the radial variable to the spatial radial variable) of the nth 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](7.1)F⌢n(ρ)=Sn{ff(r)}=∫0∞f(r)jn(ρr)r2dr.𝕊n is used to specifically denote the spherical Hankel transform of order n. The inverse transform is given by [29](7.2)f(r)=2π∫0∞F⌢n(ρ)jn(ρr)ρ2dρ.
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](7.3)Ylm(ψ,θ)=(2l+1)(l-m)!4π(l+m)!Plm(cos⁡ψ)eimθ,
where Ylm is called a spherical harmonic function of degree l and order m, Plm 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
(7.4)∫02π∫0πYlmYl′m′¯sin⁡ψdψdθ=δll′δmm′.
Here δij 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(7.5)f(ψ,θ)=∑l=0∞∑m=-llflmYlm(ψ,θ),
where (7.6)flm=∫02π∫0πf(ψ,θ)Ylm(ψ,θ)¯sin⁡ψdψdθ.
The coefficients flm are sometimes referred to as the spherical Fourier transform of f(ψ,θ) [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 asf(r⃗)=f(r,ψr,θr) as a function of spherical polar coordinates and its 3D Fourier transform is written as F(ω⃗)=F(ρ,ψω,θω) 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(7.7)f(r⃗)=f(r,ψr,θr)=∑l=0∞∑k=-llflk(r)Ylk(ψr,θr),
where Ylk(ψr,θr) are the spherical harmonics and the Fourier coefficients are given by(7.8)flk(r)=∫02π∫0πf(r,ψr,θr)Ylk(ψr,θr)¯sin⁡ψrdψrdθr.
The 3D Fourier transform of f can also be written in Fourier space in spherical polar coordinates as(7.9)F(ω⃗)=F(ρ,ψω,θω)=∑l=0∞∑k=-llFlk(ρ)Ylk(ψω,θω),
where (7.10)Flk(ρ)=∫02π∫0πF(ρ,ψω,θω)Ylk(ψω,θω)¯sin⁡ψωdψωdθω.
We emphasize again that the relationship between flk(r) and Flk(ρ) is not that of a Fourier transform. In fact, this relationship is given by [29](7.11)Flk(ρ)=4π(-i)l∫0∞flk(r)jl(ρr)r2dr=4π(-i)lSl{flk(r)},
where 𝕊l denotes a spherical Hankel transform of order l. The inverse relationship is given by(7.12)flk(r)=14π(i)l2π∫0∞Flk(ρ)jl(ρr)ρ2dρ=14π(i)lSl-1{Flk(ρ)}.


## 7.1. Definition of the Spherical Hankel Transform

The spherical Hankel transform can then be defined as [32, 33](7.1)F⌢n(ρ)=Sn{ff(r)}=∫0∞f(r)jn(ρr)r2dr.𝕊n is used to specifically denote the spherical Hankel transform of order n. The inverse transform is given by [29](7.2)f(r)=2π∫0∞F⌢n(ρ)jn(ρr)ρ2dρ.
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](7.3)Ylm(ψ,θ)=(2l+1)(l-m)!4π(l+m)!Plm(cos⁡ψ)eimθ,
where Ylm is called a spherical harmonic function of degree l and order m, Plm 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
(7.4)∫02π∫0πYlmYl′m′¯sin⁡ψdψdθ=δll′δmm′.
Here δij 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(7.5)f(ψ,θ)=∑l=0∞∑m=-llflmYlm(ψ,θ),
where (7.6)flm=∫02π∫0πf(ψ,θ)Ylm(ψ,θ)¯sin⁡ψdψdθ.
The coefficients flm are sometimes referred to as the spherical Fourier transform of f(ψ,θ) [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 asf(r⃗)=f(r,ψr,θr) as a function of spherical polar coordinates and its 3D Fourier transform is written as F(ω⃗)=F(ρ,ψω,θω) 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(7.7)f(r⃗)=f(r,ψr,θr)=∑l=0∞∑k=-llflk(r)Ylk(ψr,θr),
where Ylk(ψr,θr) are the spherical harmonics and the Fourier coefficients are given by(7.8)flk(r)=∫02π∫0πf(r,ψr,θr)Ylk(ψr,θr)¯sin⁡ψrdψrdθr.
The 3D Fourier transform of f can also be written in Fourier space in spherical polar coordinates as(7.9)F(ω⃗)=F(ρ,ψω,θω)=∑l=0∞∑k=-llFlk(ρ)Ylk(ψω,θω),
where (7.10)Flk(ρ)=∫02π∫0πF(ρ,ψω,θω)Ylk(ψω,θω)¯sin⁡ψωdψωdθω.
We emphasize again that the relationship between flk(r) and Flk(ρ) is not that of a Fourier transform. In fact, this relationship is given by [29](7.11)Flk(ρ)=4π(-i)l∫0∞flk(r)jl(ρr)r2dr=4π(-i)lSl{flk(r)},
where 𝕊l denotes a spherical Hankel transform of order l. The inverse relationship is given by(7.12)flk(r)=14π(i)l2π∫0∞Flk(ρ)jl(ρr)ρ2dρ=14π(i)lSl-1{Flk(ρ)}.


## 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 inputs(r⃗,ω) and output u(r⃗,ω). 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:
(8.1)I=∫0∞ϕ(x)x2-k2jn(xr)x2dx=-πik2hn(2)(kr)ϕ(k).
Here, ϕ is any analytic function defined on the positive real line that remains bounded as x goes to infinity, jn(x) is a spherical Bessel function of order n,  hn(2)(x) is a spherical Hankel function of order n,  and k is a wavenumber which may be real or complex, chosen such that the imaginary part of k 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 (8.2)G(ω⃗,ω)=1(ρ2-k2),
which depends on the frequency spherical radial variable ρ only.In general, as the imaginary part of the wavenumberk 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 k 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 k 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 k equal to zero (in other words, a real wavenumber) is easy to visualize. For the case of a real k, the function is discontinuous at ρ=k, meaning that only those spatial frequencies where exactly ρ=k are passed and the rest are attenuated. In other words, the resulting wave must have ρ2=ωx2+ωy2+ωz2=k2. 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 thefull 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(8.3)gn3D(r,k)=Sn-1{1ρ2-k2}=2π∫0∞1ρ2-k2jn(ρr)ρ2dρ=-ikhn(2)(kr),
where gn3D(r,k) has been used to denote the nth order Green function coefficient, with the subscript indicating the order of the coefficient and the wavenumber k 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(8.4)S(ω⃗,ω)=∑l=0∞∑m=-llSlm(ρ,ω)Ylm(ψω,θω).Hence, (4.1) for the output wavefield becomes (8.5)U(ω⃗,ω)=∑l=0∞∑m=-llSlm(ρ,ω)ρ2-k2Ylm(ψω,θω).
The symbol k is used to denote the wavenumber. The spatial inverse Fourier transform and thus the output in spatial coordinates is given by(8.6)u(r,ψr,θr,ω)=∑l=0∞∑m=-ll(i)l4π{2π∫0∞Slm(ρ,ω)jl(ρr)ρ2-k2ρ2dρ}Ylm(ψr,θr).
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 obtain(8.7)2π∫0∞Slm(ρ,ω)jl(ρr)ρ2-k2ρ2dρ=gl3D(r,k)Slm(k,ω).
It now follows that the wavefield expression becomes(8.8)u(r⃗,ω)=∑l=0∞∑m=-llil4πgl3D(r,k)Slm(k,ω)Ylm(ψr,θr).
If the measured wavefield itself, namely, the left-hand side of (8.8), is expanded as in a spherical harmonic series so that (8.9)u(r⃗,ω)=u(r,ψr,θr,ω)=∑l=0∞∑m=-llulm(r,ω)Ylm(ψr,θr),
then (8.8) gives us the simple input-output transfer function relationship we seek:(8.10)ulm(r,ω)=il4πgl3D(r,k)Slm(k,ω).
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 Slm(k,ω) and the resulting wavefield ulm(r,ω).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 (8.11)Slm(ρ,ω)=4π(-i)mSl{slm(r,ω)}=4π(-i)l∫0∞slm(r,ω)jl(ρr)r2dr.Hence (8.10) can be written as(8.12)ulm(r,ω)=gl3D(r,k)∫0∞slm(x,ω)jl(kx)x2dx,
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:(8.13)ulm(r,ω)=gl3D(r,k)Sl{slm(k,ω)}.
The value of defining the Green function coefficients now becomes apparent in that it permits the relationship between ulm(r,ω) and slm(k,ω) 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 thats(r⃗,ω)=q(r⃗)η(ω). This would be the case for a spatial inhomogeneity and a temporal input. In this case, then slm(r,ω)=qlm(r)η(ω) and the formula for the forward problem is given by (8.14)ulm(r,ω)=gl3D(r,k)η(ω)∫0∞qlm(r)jl(kr)r2dr=gl3D(r,k)η(ω)Sl{qlm(k)}.


### 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 (l,m)th term of the input and the (l,m)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 ρ=k. Note that (l,m)th order output coefficients are related to (l,m)th order transfer function coefficients in proportion to the lth order spherical Hankel transform of the (l,m)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 ρ=k 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 ρ=k 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 (l,m) 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.


## 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 (8.2)G(ω⃗,ω)=1(ρ2-k2),
which depends on the frequency spherical radial variable ρ only.In general, as the imaginary part of the wavenumberk 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 k 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 k 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 k equal to zero (in other words, a real wavenumber) is easy to visualize. For the case of a real k, the function is discontinuous at ρ=k, meaning that only those spatial frequencies where exactly ρ=k are passed and the rest are attenuated. In other words, the resulting wave must have ρ2=ωx2+ωy2+ωz2=k2. 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 thefull 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(8.3)gn3D(r,k)=Sn-1{1ρ2-k2}=2π∫0∞1ρ2-k2jn(ρr)ρ2dρ=-ikhn(2)(kr),
where gn3D(r,k) has been used to denote the nth order Green function coefficient, with the subscript indicating the order of the coefficient and the wavenumber k 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(8.4)S(ω⃗,ω)=∑l=0∞∑m=-llSlm(ρ,ω)Ylm(ψω,θω).Hence, (4.1) for the output wavefield becomes (8.5)U(ω⃗,ω)=∑l=0∞∑m=-llSlm(ρ,ω)ρ2-k2Ylm(ψω,θω).
The symbol k is used to denote the wavenumber. The spatial inverse Fourier transform and thus the output in spatial coordinates is given by(8.6)u(r,ψr,θr,ω)=∑l=0∞∑m=-ll(i)l4π{2π∫0∞Slm(ρ,ω)jl(ρr)ρ2-k2ρ2dρ}Ylm(ψr,θr).
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 obtain(8.7)2π∫0∞Slm(ρ,ω)jl(ρr)ρ2-k2ρ2dρ=gl3D(r,k)Slm(k,ω).
It now follows that the wavefield expression becomes(8.8)u(r⃗,ω)=∑l=0∞∑m=-llil4πgl3D(r,k)Slm(k,ω)Ylm(ψr,θr).
If the measured wavefield itself, namely, the left-hand side of (8.8), is expanded as in a spherical harmonic series so that (8.9)u(r⃗,ω)=u(r,ψr,θr,ω)=∑l=0∞∑m=-llulm(r,ω)Ylm(ψr,θr),
then (8.8) gives us the simple input-output transfer function relationship we seek:(8.10)ulm(r,ω)=il4πgl3D(r,k)Slm(k,ω).
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 Slm(k,ω) and the resulting wavefield ulm(r,ω).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 (8.11)Slm(ρ,ω)=4π(-i)mSl{slm(r,ω)}=4π(-i)l∫0∞slm(r,ω)jl(ρr)r2dr.Hence (8.10) can be written as(8.12)ulm(r,ω)=gl3D(r,k)∫0∞slm(x,ω)jl(kx)x2dx,
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:(8.13)ulm(r,ω)=gl3D(r,k)Sl{slm(k,ω)}.
The value of defining the Green function coefficients now becomes apparent in that it permits the relationship between ulm(r,ω) and slm(k,ω) 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 thats(r⃗,ω)=q(r⃗)η(ω). This would be the case for a spatial inhomogeneity and a temporal input. In this case, then slm(r,ω)=qlm(r)η(ω) and the formula for the forward problem is given by (8.14)ulm(r,ω)=gl3D(r,k)η(ω)∫0∞qlm(r)jl(kr)r2dr=gl3D(r,k)η(ω)Sl{qlm(k)}.


## 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 (l,m)th term of the input and the (l,m)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 ρ=k. Note that (l,m)th order output coefficients are related to (l,m)th order transfer function coefficients in proportion to the lth order spherical Hankel transform of the (l,m)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 ρ=k 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 ρ=k 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 (l,m) 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:
(9.1)I=∫0∞ϕ(ρ)ρ2-k2Jn(ρr)ρdρ=-πi12Hn(2)(kr)ϕ(k).
Here, ϕ is an analytic function defined on the positive real line that remains bounded as x goes to infinity, Jn(x) is a Bessel function of order n,  Hn(2)(x) is a Hankel function of order n,  and k is a wavenumber which may be real or complex, chosen such that the imaginary part of k 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(9.2)G(ω⃗,ω)=1(ρ2-k2),
which depends on the frequency radial variable only.As for the 3D case, we define the 2D Green function coefficients as thenth-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(9.3)gn2D(r,k)=Hn-1{1ρ2-k2}=∫0∞1ρ2-k2Jn(ρr)ρdρ=-πi2Hn(2)(kr).
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), gn2D(r,k) 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)S(ρ,ψ,ω)=∑n=-∞∞Sn(ρ,ω)ejnψ,
where Sn(ρ,ω) can be found from (9.5)Sn(ρ,ω)=12π∫02πS(ρ,ψ,ω)e-jnψdψ.
The output wavefield u(r⃗,ω) is similarly expanded as a series so that in the spatial domain we can write(9.6)u(r⃗,ω)=u(r,θ,ω)=∑n=-∞∞un(r)ejnθ,
while in the spatial Fourier domain this is(9.7)U(ω⃗,ω)=U(ρ,ψ,ω)=∑n=-∞∞Un(ρ,ω)ejnψ.
The output wavefield is given by U(ω⃗,ω)=G(ω⃗,ω)S(ω⃗,ω) and since G(ω⃗,ω) is radially symmetric and does not require a full series, the output wavefield in the Fourier domain is given by(9.8)U(ρ,ψ,ω)=∑n=-∞∞Un(ρ,ω)ejnψ=∑n=-∞∞Sn(ρ,ω)(ρ2-k2)ejnψ,
or equivalently as(9.9)Un(ρ,ω)=Sn(ρ,ω)(ρ2-k2).
The equivalent expression to (9.8) in the spatial domain is given by(9.10)u(r⃗,ω)=∑n=-∞∞un(r)ejnθ=∑n=-∞∞in2π{∫0∞Sn(ρ,ω)(ρ2-k2)Jn(ρr)ρdρ}ejnθ,
or(9.11)un(r)=in2π∫0∞Sn(ρ,ω)(ρ2-k2)Jn(ρr)ρdρ.
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 as(9.12)∫0∞Sn(ρ,ω)(ρ2-k2)Jn(ρr)ρdρ=-πi2Hn(2)(kr)Sn(k,ω)=gn2D(r,k)Sn(k,ω).
Hence, the output wavefield expression is given by(9.13)u(r⃗,ω)=∑n=-∞∞un(r)ejnθ=∑n=-∞∞in2πgn2D(r,k)Sn(k,ω)ejnθ.
The definition of forward and inverse Hankel transforms, the general relationship between Hankel and Fourier transforms is given by (9.14)Sn(ρ)=2πi-nHn(sn(r))=2πi-n∫0∞sn(r)Jn(ρr)rdr,
and may be used to simplify (9.13) to yield(9.15)u(r⃗,ω)=∑n=-∞∞un(r)ejnθ=∑n=-∞∞gn2D(r,k)∫0∞sn(x,ω)Jn(kx)xdxejnθ,
or more compactly as (9.16)un(r,ω)=gn2D(r,k)∫0∞sn(x,ω)Jn(kx)xdx.
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 (9.17)un(r,ω)=gn2D(r,k)Hn(sn(k,ω)).
If it is also further assumed that the inhomogeneity function is separable in time and space (a fairly general assumption) so that s(r⃗,ω)=q(r⃗)η(ω), and(9.18)∑n=-∞∞sn(r,ω)ejnθ=η(ω)∑n=-∞∞qn(r)ejnθ,
then (9.16) can be further reduced to (9.19)un(r,ω)=gn2D(r,k)η(ω)Hn{qn(k)}.
It is noted that the value of the measured wavefield at any position r 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.


## 9.1. 2D Green Function Coefficients

The overall system transfer function in 2D is given by(9.2)G(ω⃗,ω)=1(ρ2-k2),
which depends on the frequency radial variable only.As for the 3D case, we define the 2D Green function coefficients as thenth-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(9.3)gn2D(r,k)=Hn-1{1ρ2-k2}=∫0∞1ρ2-k2Jn(ρr)ρdρ=-πi2Hn(2)(kr).
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), gn2D(r,k) 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)S(ρ,ψ,ω)=∑n=-∞∞Sn(ρ,ω)ejnψ,
where Sn(ρ,ω) can be found from (9.5)Sn(ρ,ω)=12π∫02πS(ρ,ψ,ω)e-jnψdψ.
The output wavefield u(r⃗,ω) is similarly expanded as a series so that in the spatial domain we can write(9.6)u(r⃗,ω)=u(r,θ,ω)=∑n=-∞∞un(r)ejnθ,
while in the spatial Fourier domain this is(9.7)U(ω⃗,ω)=U(ρ,ψ,ω)=∑n=-∞∞Un(ρ,ω)ejnψ.
The output wavefield is given by U(ω⃗,ω)=G(ω⃗,ω)S(ω⃗,ω) and since G(ω⃗,ω) is radially symmetric and does not require a full series, the output wavefield in the Fourier domain is given by(9.8)U(ρ,ψ,ω)=∑n=-∞∞Un(ρ,ω)ejnψ=∑n=-∞∞Sn(ρ,ω)(ρ2-k2)ejnψ,
or equivalently as(9.9)Un(ρ,ω)=Sn(ρ,ω)(ρ2-k2).
The equivalent expression to (9.8) in the spatial domain is given by(9.10)u(r⃗,ω)=∑n=-∞∞un(r)ejnθ=∑n=-∞∞in2π{∫0∞Sn(ρ,ω)(ρ2-k2)Jn(ρr)ρdρ}ejnθ,
or(9.11)un(r)=in2π∫0∞Sn(ρ,ω)(ρ2-k2)Jn(ρr)ρdρ.
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 as(9.12)∫0∞Sn(ρ,ω)(ρ2-k2)Jn(ρr)ρdρ=-πi2Hn(2)(kr)Sn(k,ω)=gn2D(r,k)Sn(k,ω).
Hence, the output wavefield expression is given by(9.13)u(r⃗,ω)=∑n=-∞∞un(r)ejnθ=∑n=-∞∞in2πgn2D(r,k)Sn(k,ω)ejnθ.
The definition of forward and inverse Hankel transforms, the general relationship between Hankel and Fourier transforms is given by (9.14)Sn(ρ)=2πi-nHn(sn(r))=2πi-n∫0∞sn(r)Jn(ρr)rdr,
and may be used to simplify (9.13) to yield(9.15)u(r⃗,ω)=∑n=-∞∞un(r)ejnθ=∑n=-∞∞gn2D(r,k)∫0∞sn(x,ω)Jn(kx)xdxejnθ,
or more compactly as (9.16)un(r,ω)=gn2D(r,k)∫0∞sn(x,ω)Jn(kx)xdx.
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 (9.17)un(r,ω)=gn2D(r,k)Hn(sn(k,ω)).
If it is also further assumed that the inhomogeneity function is separable in time and space (a fairly general assumption) so that s(r⃗,ω)=q(r⃗)η(ω), and(9.18)∑n=-∞∞sn(r,ω)ejnθ=η(ω)∑n=-∞∞qn(r)ejnθ,
then (9.16) can be further reduced to (9.19)un(r,ω)=gn2D(r,k)η(ω)Hn{qn(k)}.
It is noted that the value of the measured wavefield at any position r 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 (10.1)G(ωx,ω)=1(ωx2-k2),
which depends on only a single spatial frequency. Equation (4.2) in 1D then becomes(10.2)u(r,ω)=1(2π)∫-∞∞S(ωx,ω)(ωx2-k2)eiωxrdωx.
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 (10.3)I=12π∫0∞ϕ(x)x2-k2eixrdx,
where ϕ is an analytic function defined on the positive real line that approaches zero as x 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
(10.4)12π∫-∞∞ϕ(ρ)ρ2-k2eiρrdρ={12ikϕ(-k)e-ikr,r>0,12ikϕ(k)eikr,r<0.
Here, ϕ is an analytic function defined on the positive real line that remains bounded as x goes to infinity and k is a wavenumber which may be real or complex, chosen such that the imaginary part of k 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(10.5)g1D(r,k)=12π∫-∞∞eiωxrωx2-k2dωx,
and can be evaluated with the help of Theorem 10.1 to give (10.6)g1D(r,k)={12ike-ikrr>012ikeikrr<0}=12ike-ik|r|.


### 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 (10.7)u(r,ω)=1(2π)∫-∞∞S(ωx,ω)(ωx2-k2)eiωxrdωx={12ikS(-k,ω)e-ikrr>012ikS(k,ω)eikrr<0=g1D(r,k)S(-sgn⁡(r)k,ω),
where g1D(r,k)=(1/2ik)e-ik|r|. 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 k. 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 ±k 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 thats(r,ω)=q(r)η(ω)→S(ωx,ω)=Q(ωx)η(ω), where Q(ωx) is the 1D Fourier transform of q(r). The relevant result now reads(10.8)u(r,ω)=g1D(r,k)η(ω)Q(-sgn⁡(r)k).
This result is directly applicable for computational purposes.


## 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(10.5)g1D(r,k)=12π∫-∞∞eiωxrωx2-k2dωx,
and can be evaluated with the help of Theorem 10.1 to give (10.6)g1D(r,k)={12ike-ikrr>012ikeikrr<0}=12ike-ik|r|.


## 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 (10.7)u(r,ω)=1(2π)∫-∞∞S(ωx,ω)(ωx2-k2)eiωxrdωx={12ikS(-k,ω)e-ikrr>012ikS(k,ω)eikrr<0=g1D(r,k)S(-sgn⁡(r)k,ω),
where g1D(r,k)=(1/2ik)e-ik|r|. 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 k. 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 ±k 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 thats(r,ω)=q(r)η(ω)→S(ωx,ω)=Q(ωx)η(ω), where Q(ωx) is the 1D Fourier transform of q(r). The relevant result now reads(10.8)u(r,ω)=g1D(r,k)η(ω)Q(-sgn⁡(r)k).
This result is directly applicable for computational purposes.


## 11. Summary of Results

The relationships in the previous sections are summarized in Table1.Table 1DimensionTransformsPseudo-Green's (gn) functionInput/output relationship3Spherical Hankel and spherical harmonicgn3D(r,k)=-ikhn(2)(kr)ulk(r,ω)=gl3D(r,k)η(ω)𝕊l{qlk(k)}2Hankel and Fourier Seriesgn2D(r,k)=-πi2Hn(2)(kr)un(r,ω)=gn2D(r,k)η(ω)ℍn{qn(k)}11D Fourierg1D(x)=-i2ke-ik|x|u(x,ω)=g1D(x)η(ω)Q(-sgn(x)k)


## 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 origin(12.1)∇2u(r⃗,t)-1c2∂2∂t2u(r⃗,t)=-δ(r⃗)δ(t).This corresponds to the Helmholtz equation above withk=ω/c. The temporal and spatial Fourier transform of δ(r⃗)δ(t) 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(12.2)u(r,ω)=g03D(r,k)⋅1=-ikh0(2)(kr).
In 2D, the relationship is given by(12.3)u(r,ω)=g02D(r,k)⋅1=-πi2H0(kr).
In 1D, the relationship is(12.4)u(r,ω)=g1D(r,k)⋅1=12ike-ik|r|.
These three relationships can now be inverse Fourier transformed in time.For the 3D case, we use the fact thath0(2)(kr)=j0(kr)-iy0(kr)=iexp(-ikr)/kr so that the inverse Fourier transform of (12.2) gives(12.5)u(r,t)=12π∫-∞∞1re-iωr/ceiωtdω=1rδ(t-rc).
The 1D time response can be found from the inverse Fourier transform of  (12.4)(12.6)u(r,t)=12π∫-∞∞c2iωe-iω|r|/ceiωtdω=c4[2H(t-|r|c)-1]=c4sgn⁡(t-|r|c),
where H(x) 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](12.7)H0(2)(x)=J0(x)-iY0(x)=2iπ∫1∞e-ixττ2-1dτ.
Hence from (12.7) and (12.3), the inverse Fourier transform of (12.3) gives(12.8)u(r,t)=12π∫-∞∞∫1∞e-iωrτ/cτ2-1dτeiωtdω.
Changing the order of integration gives(12.9)u(r,t)=∫1∞1τ2-112π∫-∞∞eiω[t-rτ/c]dωdτ.
But (12.10)12π∫-∞∞eiω[t-rτ/c]dω=δ(t-rτc),
so that (12.9) becomes(12.11)u(r,t)=∫1∞1τ2-1δ(t-rτc)dτ=∫-∞∞1τ2-1H(τ-1)δ(t-rτc)dτ.
Changing variables so that x=rτ/c gives(12.12)u(r,t)=cr∫-∞∞H(cxr-1)1(cx/r)2-1δ(t-x)dx=H(t-rc)1t2-r2/c2,
which yields the desired result in closed form.


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

In the case of a real wavenumberk, 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 k. 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 (13.1)iulk(r0,ω)kshl(2)(ksr0)η(ω)=∫0∞ϕlk(x)jl(ksx)x2dx,
where x has been used as a dummy integration variable in order to avoid possible confusion and r0 has been used as the radial variable and indicates the position where a measurement of the wavefield is made. Multiplying both sides by jl(ksr)ks2 and integrating over ks gives(13.2)∫0∞iulk(r0,ω)kshl(2)(ksr0)η(ω)jl(ksr)ks2dks=∫0∞∫0∞ϕlk(x)jl(ksx)x2dxjl(ksr)ks2dks=∫0∞∫0∞ϕlk(x)jl(ksx)jl(ksr)ks2dksx2dx.Using the orthogonality of the spherical Bessel functions, it follows that (13.2) becomes (13.3)∫0∞iulk(r0,ω)kshl(2)(ksr0)η(ω)jl(ksr)ks2dks=∫0∞ϕlk(x)π2x2δ(x-r)x2dx=π2ϕlk(r).
Recalling that ks=ω/cs, the inversion formula for the spatial source becomes (13.4)ϕlk(r)=2iπcs2∫0∞ulk(r0,ω)hl(2)(ωr0/cs)η(ω)jl(ωrcs)ωdω.
If the measured wavefield is also separable in space and time so that ulk(r0,ω)=χlk(r0)T(ω), then (13.4) can be further simplified to (13.5)ϕlk(r)=2iχlk(r0)πcs2∫0∞T(ω)hl(2)(ωr0/cs)η(ω)jl(ωrcs)ωdω.
Equation (13.5) gives a simple form for finding the spatial source function at any position r once having measured the temporal response of the wavefield at a position r0. However, since the spherical harmonic transform ulk(r0,ω) of the wavefield is required, this implies that sufficient angular information about the wavefield u 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 s(r⃗,ω)=ϕ(r⃗)η(ω). In this case, then slk(r,ω)=ϕlk(r)η(ω) and (8.12) becomes (13.6)2iun(r0,ω)πHn(2)(ksr0)η(ω)=∫0∞ϕn(x)Jn(ksx)xdx,
where x has been used as a dummy integration variable in order to avoid possible confusion and r0 indicates the radial position where the wavefield is measured. Multiplying both sides of (13.6) by Jn(ksr)ks and integrating over ks gives(13.7)∫0∞2iun(r0,ω)πHn(2)(ksr0)η(ω)Jn(ksr)ksdks=∫0∞∫0∞ϕn(x)Jn(ksx)xdxJn(ksr)ksdks,=∫0∞ϕn(x)∫0∞Jn(ksx)Jn(ksr)ksdksxdx.
Using the orthogonality of the spherical Bessel functions, it follows that (13.7) becomes (13.8)∫0∞2iun(r0,ω)πHn(2)(ksr0)η(ω)Jn(ksr)ksdks=∫0∞ϕn(x)1xδ(x-r)xdx=ϕn(r).
Recalling that ks=ω/cs, the inversion formula for the inhomogeneity becomes (13.9)ϕn(r)=2iπcs2∫0∞un(r0,ω)Hn(2)(ωr0/cs)η(ω)Jn(ωrcs)ωdω.
If the measured wavefield is also separable in space and time so that un(r0,ω)=χn(r0)T(ω), then (13.4) can be further simplified to (13.10)ϕn(r)=2iχn(r0)πcs2∫0∞T(ω)Hn(2)(ωr0/cs)η(ω)Jn(ωrcs)ωdω.
Equation (13.10) gives a simple form for finding the spatial source function at any position r once having measured the wavefield response at a radial position r0. 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 thatu(r,ω)=χ(r)T(ω), which leads to to(13.11)ϕ̃(-ks)=χ(r)T(ω)η(ω)2ikseiksr=∫-∞∞ϕ(y)eiyksdyr>0,ϕ̃(ks)=χ(r)T(ω)η(ω)2ikse-iksr=∫-∞∞ϕ(y)e-iyksdyr<0.
Using the orthogonality of the Fourier kernel, both sides are multiplied by e±ixks and integrated over all ks to give(13.12)χ(r)∫-∞∞T(ω)η(ω)2ikseiksre-ixksdks=∫-∞∞ϕ(y)∫-∞∞eiykse-ixksdksdy=2π∫-∞∞ϕ(y)δ(y-x)dyr>0⟹ϕ(x)=iχ(r)πcs2∫-∞∞T(ω)η(ω)eiω((r-x)/cs)ωdωr>0,
and similarly(13.13)χ(r)∫-∞∞T(ω)η(ω)2ikse-iksreixksdks=∫-∞∞ϕ(y)∫-∞∞e-iykseixksdksdy=2π∫-∞∞ϕ(y)δ(y-x)dyr<0⟹ϕ(x)=iχ(r)πcs2∫-∞∞T(ω)η(ω)eiω((x-r)/cs)ωdωr<0.
Several points need to be made regarding (13.12) and (13.13). First, they both give the source function as a function of position x along the real line. The interpretation of the variable r 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 (T(ω)/η(ω))ω, evaluated at (x-r)/cs or (r-x)/cs, depending on whether measurements are made in transmission or reflection. Clearly |x-r|/cs is the time taken for a wave to travel the distance |x-r|.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 (T(ω)/η(ω))ω 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 |x-r|/cs.


## 13.1. Inversion Formula in 3D

Equation (8.14) admits an inversion useful for tomographic applications. It can be written as (13.1)iulk(r0,ω)kshl(2)(ksr0)η(ω)=∫0∞ϕlk(x)jl(ksx)x2dx,
where x has been used as a dummy integration variable in order to avoid possible confusion and r0 has been used as the radial variable and indicates the position where a measurement of the wavefield is made. Multiplying both sides by jl(ksr)ks2 and integrating over ks gives(13.2)∫0∞iulk(r0,ω)kshl(2)(ksr0)η(ω)jl(ksr)ks2dks=∫0∞∫0∞ϕlk(x)jl(ksx)x2dxjl(ksr)ks2dks=∫0∞∫0∞ϕlk(x)jl(ksx)jl(ksr)ks2dksx2dx.Using the orthogonality of the spherical Bessel functions, it follows that (13.2) becomes (13.3)∫0∞iulk(r0,ω)kshl(2)(ksr0)η(ω)jl(ksr)ks2dks=∫0∞ϕlk(x)π2x2δ(x-r)x2dx=π2ϕlk(r).
Recalling that ks=ω/cs, the inversion formula for the spatial source becomes (13.4)ϕlk(r)=2iπcs2∫0∞ulk(r0,ω)hl(2)(ωr0/cs)η(ω)jl(ωrcs)ωdω.
If the measured wavefield is also separable in space and time so that ulk(r0,ω)=χlk(r0)T(ω), then (13.4) can be further simplified to (13.5)ϕlk(r)=2iχlk(r0)πcs2∫0∞T(ω)hl(2)(ωr0/cs)η(ω)jl(ωrcs)ωdω.
Equation (13.5) gives a simple form for finding the spatial source function at any position r once having measured the temporal response of the wavefield at a position r0. However, since the spherical harmonic transform ulk(r0,ω) of the wavefield is required, this implies that sufficient angular information about the wavefield u 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 s(r⃗,ω)=ϕ(r⃗)η(ω). In this case, then slk(r,ω)=ϕlk(r)η(ω) and (8.12) becomes (13.6)2iun(r0,ω)πHn(2)(ksr0)η(ω)=∫0∞ϕn(x)Jn(ksx)xdx,
where x has been used as a dummy integration variable in order to avoid possible confusion and r0 indicates the radial position where the wavefield is measured. Multiplying both sides of (13.6) by Jn(ksr)ks and integrating over ks gives(13.7)∫0∞2iun(r0,ω)πHn(2)(ksr0)η(ω)Jn(ksr)ksdks=∫0∞∫0∞ϕn(x)Jn(ksx)xdxJn(ksr)ksdks,=∫0∞ϕn(x)∫0∞Jn(ksx)Jn(ksr)ksdksxdx.
Using the orthogonality of the spherical Bessel functions, it follows that (13.7) becomes (13.8)∫0∞2iun(r0,ω)πHn(2)(ksr0)η(ω)Jn(ksr)ksdks=∫0∞ϕn(x)1xδ(x-r)xdx=ϕn(r).
Recalling that ks=ω/cs, the inversion formula for the inhomogeneity becomes (13.9)ϕn(r)=2iπcs2∫0∞un(r0,ω)Hn(2)(ωr0/cs)η(ω)Jn(ωrcs)ωdω.
If the measured wavefield is also separable in space and time so that un(r0,ω)=χn(r0)T(ω), then (13.4) can be further simplified to (13.10)ϕn(r)=2iχn(r0)πcs2∫0∞T(ω)Hn(2)(ωr0/cs)η(ω)Jn(ωrcs)ωdω.
Equation (13.10) gives a simple form for finding the spatial source function at any position r once having measured the wavefield response at a radial position r0. 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 thatu(r,ω)=χ(r)T(ω), which leads to to(13.11)ϕ̃(-ks)=χ(r)T(ω)η(ω)2ikseiksr=∫-∞∞ϕ(y)eiyksdyr>0,ϕ̃(ks)=χ(r)T(ω)η(ω)2ikse-iksr=∫-∞∞ϕ(y)e-iyksdyr<0.
Using the orthogonality of the Fourier kernel, both sides are multiplied by e±ixks and integrated over all ks to give(13.12)χ(r)∫-∞∞T(ω)η(ω)2ikseiksre-ixksdks=∫-∞∞ϕ(y)∫-∞∞eiykse-ixksdksdy=2π∫-∞∞ϕ(y)δ(y-x)dyr>0⟹ϕ(x)=iχ(r)πcs2∫-∞∞T(ω)η(ω)eiω((r-x)/cs)ωdωr>0,
and similarly(13.13)χ(r)∫-∞∞T(ω)η(ω)2ikse-iksreixksdks=∫-∞∞ϕ(y)∫-∞∞e-iykseixksdksdy=2π∫-∞∞ϕ(y)δ(y-x)dyr<0⟹ϕ(x)=iχ(r)πcs2∫-∞∞T(ω)η(ω)eiω((x-r)/cs)ωdωr<0.
Several points need to be made regarding (13.12) and (13.13). First, they both give the source function as a function of position x along the real line. The interpretation of the variable r 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 (T(ω)/η(ω))ω, evaluated at (x-r)/cs or (r-x)/cs, depending on whether measurements are made in transmission or reflection. Clearly |x-r|/cs is the time taken for a wave to travel the distance |x-r|.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 (T(ω)/η(ω))ω 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 |x-r|/cs.


## 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 then-dimensional sphere of radius k, where k is the wavenumber and n is the dimension of the space. Combined together, these form the final output response. This is true for dimensions n=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.



---

*Source: 478295-2011-11-24.xml*