The Spherical Harmonic Family of Beampatterns

: The free space solution to the wave equation in spherical coordinates is well known as a separable product of functions. Re-examination of these functions, particularly the sums of spherical Bessel and harmonic functions, reveals behaviors which can produce a range of useful beampatterns from radially symmetric sources. These functions can be modiﬁed by several key parameters which can be adjusted to produce a wide-ranging family of beampatterns, from the axicon Bessel beam to a variety of unique axial and lateral forms. We demonstrate that several special properties of the simple sum over integer orders of spherical Bessel functions, and then the sum of their product with spherical harmonic functions specifying the free space solution, lead to a family of useful beampatterns and a unique framework for designing them. Examples from a simulation of a pure tone 5 MHz ultrasound conﬁguration demonstrate strong central axis concentration, and the ability to modulate or localize the axial intensity with simple adjustment of the integer orders and other key parameters related to the weights and arguments of the spherical Bessel functions.


Introduction
In the design of coherent imaging systems using optical or ultrasound excitation, there has been a longstanding interest in the behavior of beampatterns. One active area of interest is the field of localized waves or weakly diffracting beams [1]. These include the Bessel beam [2,3], the X-wave [4], and the needle pulse [5]. Ultrasound imaging systems with phased array transducers traditionally rely on focusing and apodization strategies to achieve high spatial resolution in the transverse direction [6,7] and in this context a key consideration is the beamwidth as a function of depth, along with the minimization of sidelobes [8,9]. A useful concept in designing beampatterns is the Fourier transform relation that applies to the field on a source plane and at some farfield or focal depth [10,11]. This, combined with related approximations, enables the specification and selection of beams from practical systems [12].
In this paper, an alternative framework is examined for specifying radially symmetric beampatterns, by using interpretations of the free-space separable solution to the Helmholtz equation. In particular, the sum over integer orders of the spherical Bessel and harmonic (SBH) product terms has unique properties that produce a number of useful implementations. First, localized central axis beams can be produced, with a variety of axial profiles linked directly to the selection of the integer orders. Second, some other key parameters related to the weights and arguments of the functions, can be specified to create a family of beampatterns ranging from the classical Bessel beam produced by an optical axicon of limited aperture [2,3] to a range of axial and transverse patterns. These beampatterns can be adjusted within the framework of the sums of SBH terms and without necessary recourse to traditional treatments of focusing and Fourier transform operations. Examples are provided from theoretical treatments of a 5 MHz flat radially symmetric transducer.

General Solution
Let us reconsider the classic separable solution for monochromatic waves in a spherical coordinate system, initially assuming a source surface that is conically shaped with prescribed boundary conditions radiating into the interior free space within a cone ( Figure 1).
OR PEER REVIEW 2 Examples are provided from theoretical treatments of a 5 MHz flat radially symmetric transducer.

General Solution
Let us reconsider the classic separable solution for monochromatic waves in a spherical coordinate system, initially assuming a source surface that is conically shaped with prescribed boundary conditions radiating into the interior free space within a cone ( Figure 1).

Figure 1.
Geometry for expansion of spherical waves within a cone of polar angle , and with the surface of the cone being an active source capable of generating spatial distributions in the form of spherical Bessel functions and sums of spherical Bessel functions. The azimuthal angle is , the radial coordinate lies on the active source, and z represents the zenith direction.
In spherical coordinates, the free space Helmholtz equation is known to be separable [13][14][15][16] and can be written as: Here, ( ) and ( ) are the spherical Bessel functions, ( , ) are the spherical harmonics [17], and and are the amplitudes assigned to each of the functions. In the summation limits, n represents the integer orders of the spherical Bessel functions and m represents the finite number of azimuthal angle orders. Note that this form provides general solutions, and requires boundary conditions to be specified in any specific case [18,19]. Now to simplify these, let us restrict our examinations to the solutions that are finite at r = 0, are symmetric about , and are comprised of a limited sum of integer orders up to a maximum of N: To visualize the different terms, Figure 2 shows the first few spherical Bessel functions of increasing order, and Figure 3 displays several low order spherical harmonic functions Y as a function of . In spherical coordinates, the free space Helmholtz equation is known to be separable [13][14][15][16] and can be written as: Here, j ι (kr) and y ι (kr) are the spherical Bessel functions, Y m ι (θ, ϕ) are the spherical harmonics [17], and a nm and b nm are the amplitudes assigned to each of the functions. In the summation limits, n represents the integer orders of the spherical Bessel functions and m represents the finite number of azimuthal angle orders. Note that this form provides general solutions, and requires boundary conditions to be specified in any specific case [18,19]. Now to simplify these, let us restrict our examinations to the solutions that are finite at r = 0, are symmetric about ϕ, and are comprised of a limited sum of integer orders up to a maximum of N: To visualize the different terms, Figure 2 shows the first few spherical Bessel functions of increasing order, and Figure 3 displays several low order spherical harmonic functions Y as a function of θ. Now, consistent with the classical treatments of Green's functions and the Rayleigh-Sommerfield integral theorems [20,21], we may reinterpret Equations (1) and (2) as a design specification on boundary conditions. If a single spherical Bessel function j n (kr) can be excited as a source function across the surface of the cone shown in Figure 1, in steady state, the corresponding product of that function with Y n (θ) will be uniquely generated in the interior free space region. Furthermore, the principle of superposition applies directly as indicated in Equation (2), meaning that if we can simultaneously establish the sum of spherical Bessel functions on the surface of the cone, the resulting field as a function of θ will be composed of the summation over all orders of corresponding Y n . Here, we are not talking about summation over frequency ω (or wavenumber k), but instead summation over integer orders of spherical Bessel functions of the same frequency. Thus, a closer examination of summations is warranted. Acoustics 2022, 4 FOR PEER REVIEW 3  Now, consistent with the classical treatments of Green's functions and the Rayleigh-Sommerfield integral theorems [20,21], we may reinterpret Equations (1) and (2) as a design specification on boundary conditions. If a single spherical Bessel function can be excited as a source function across the surface of the cone shown in Figure 1, in steady state, the corresponding product of that function with will be uniquely generated in the interior free space region. Furthermore, the principle of superposition applies directly as indicated in Equation (2), meaning that if we can simultaneously establish the sum of spherical Bessel functions on the surface of the cone, the resulting field as a function of will be composed of the summation over all orders of corresponding . Here, we are not talking about summation over frequency (or wavenumber k), but instead summation over integer orders of spherical Bessel functions of the same frequency. Thus, a closer examination of summations is warranted.

Special Properties of Sums of Spherical Bessel Functions
Because of the asymptotic limit of ( ) n j x being equal to [17], effectively the long tails of any group of four successive orders of equal amplitude and phase will trend towards zero. Meanwhile the early maximum or peak of each order n is progressively delayed from the origin and decreases as √ . Therefore, the summation series has properties  Now, consistent with the classical treatments of Green's functions and the Rayleigh-Sommerfield integral theorems [20,21], we may reinterpret Equations (1) and (2) as a design specification on boundary conditions. If a single spherical Bessel function can be excited as a source function across the surface of the cone shown in Figure 1, in steady state, the corresponding product of that function with will be uniquely generated in the interior free space region. Furthermore, the principle of superposition applies directly as indicated in Equation (2), meaning that if we can simultaneously establish the sum of spherical Bessel functions on the surface of the cone, the resulting field as a function of will be composed of the summation over all orders of corresponding . Here, we are not talking about summation over frequency (or wavenumber k), but instead summation over integer orders of spherical Bessel functions of the same frequency. Thus, a closer examination of summations is warranted.

Special Properties of Sums of Spherical Bessel Functions
Because of the asymptotic limit of ( ) n j x being equal to [17], effectively the long tails of any group of four successive orders of equal amplitude and phase will trend towards zero. Meanwhile the early maximum or peak of each order n is progressively delayed from the origin and decreases as √ . Therefore, the summation series has properties . The zero order is flat and successive orders are increasingly higher order in cosine terms. These are also related to the Legendre polynomials. When θ = π/2, at the surface of a flat piston transducer, the values of Y 0 n (π/2) reduce to simple oscillations around ±0.318 as a function of increasing n with a period of 4n. The colors correspond to integer orders 0 (blue line), 1 (yellow line), 2 (green line), and 3 (orange line).

Special Properties of Sums of Spherical Bessel Functions
Because of the asymptotic limit of j n (x) being equal to [17], effectively the long tails of any group of four successive orders of equal amplitude and phase will trend towards zero. Meanwhile the early maximum or peak of each order n is progressively delayed from the origin and decreases as √ n. Therefore, the summation series has properties that, in the context of signal processing, appear to be useful as an interpolation function for a sampled, bandlimited signal. Specifically: Thus, any arbitrary (well-behaved, bandlimited) sampled function can be reconstructed approximately as the sum of spherical Bessel functions. The sum described by Equation (5) is shown in Figure 4, resembling a band limited window function over the range of 0 < x < 200.

( )
Thus, any arbitrary (well-behaved, bandlimited) sampled function can be reconstructed approximately as the sum of spherical Bessel functions. The sum described by Equation (5) is shown in Figure 4, resembling a band limited window function over the range of 0 < x < 200. Additionally, the issues of apodization and truncation of a source can be recast by the apodization and termination of the series. This representation is similar to the Neumann series of Bessel functions [22], albeit with simpler interpretation of the coefficients. An example is shown in Figure 5. The advantages of using the spherical Bessel functions as interpolation functions in this way is that each order is matched to corresponding spherical harmonic functions, enabling the full description of a wave field in three-dimensional free space as a precise analytical expression. Furthermore, in the case of the active area of a transducer surface, each integer order or group of orders define a radial region on the transducer and in the wave produced in free space. These properties will be illustrated in the subsequent sections. Additionally, the issues of apodization and truncation of a source can be recast by the apodization and termination of the series. This representation is similar to the Neumann series of Bessel functions [22], albeit with simpler interpretation of the coefficients. An example is shown in Figure 5. The advantages of using the spherical Bessel functions as interpolation functions in this way is that each order is matched to corresponding spherical harmonic functions, enabling the full description of a wave field in three-dimensional free space as a precise analytical expression. Furthermore, in the case of the active area of a transducer surface, each integer order or group of orders define a radial region on the transducer and in the wave produced in free space. These properties will be illustrated in the subsequent sections.

Modifications via Imaginary Shift in Coordinates
Alonso et al. considered a modification to the classical solution of spherical harmonic functions in terms of an offset consisting of a purely imaginary constant [23]. Including a simple offset in the argument of Equation (1) is also a valid solution to the wave equation, thus → + is substituted in all expressions, where i is the imaginary unit index and q is an arbitrary constant. This imaginary parameter can be shown to remap the angular spectrum in a geometrical transformation around the unit circle defined by . Practically speaking, this parameter q is then used as a parameter which can significantly alter the beampattern and its angular spectrum.

Modifications via Imaginary Shift in Coordinates
Alonso et al. considered a modification to the classical solution of spherical harmonic functions in terms of an offset consisting of a purely imaginary constant [23]. Including a simple offset in the argument of Equation (1) is also a valid solution to the wave equation, thus z → z + iq is substituted in all expressions, where i is the imaginary unit index and q is an arbitrary constant. This imaginary parameter can be shown to remap the angular spectrum in a geometrical transformation around the unit circle defined by θ. Practically speaking, this parameter q is then used as a parameter which can significantly alter the beampattern and its angular spectrum.

Other Simple Modifications and Properties
In addition to the use of the a n in Equation (2) for apodization as shown in Figure 5, another tactic for implementation at the transducer surface is to prescribe the a n as a modulation with respect to the integer orders n. This effectively changes the angular spectrum and the spatial location of the peak region of the beampattern in free space. In the context of Equation (1), we note that both the j n and the Y 0 n (θ = π/2) have a natural oscillation, with values repeating or approximately repeating with ascending n, every n = 4 integer orders. Thus, a modulation function such as: a n = cos 2π n N mod + φ will influence the summation, and when N mod = 4, this imposed modulation of integer orders matches the natural modulation of the j n and Y 0 n .

Methods
Calculations of the free space fields were performed using Mathematica (Version 13.0.1, Wolfram Research, Champaign, IL, USA) and using a working precision of 20 or more decimal places in all calculations.
In all the following examples, we specify conditions related to a 5 MHz medical ultrasound system, with wavelength approximately 0.3 mm and with the summation of spherical Bessel integer orders n up to 200. This upper limit kr max = 200 corresponds to an active high amplitude transducer radius of approximately 10 mm, with a residual decay beyond that radius. We also assume a flat transducer so that the cone angle of Figure 1 is set to π/2; this models the use of a conventional piston transducer, albeit with active control of the amplitude and phase of the excitation across the transducer.

Results
For the simple summation of Equation (2) up to n = 200, we examine the cases referred to in Equations (4) and (5); the unweighted sum (a n = 1) and the square root weighted sum (a n = n 1/2 ). These display a strong central axis amplitude with oscillations in the radial direction as shown in Figure 6.

Results
For the simple summation of Equation (2) up to n = 200, we examine the cases referred to in Equations (4) and (5); the unweighted sum (an = 1) and the square root weighted sum (an = n 1/2 ). These display a strong central axis amplitude with oscillations in the radial direction as shown in Figure 6.  Axial plots of these are shown in Figure 7. Transverse cuts of these two cases at 8 mm axial range are shown in Figure 8. Despite the difference in axial amplitudes (above), both of the transverse beampatterns are a close fit to the Bessel function of zero order: J0[kx], thus resembling the classical Bessel beam described by Durnin as implemented within the framework of the optical axicon [2,3].  Transverse cuts of these two cases at 8 mm axial range are shown in Figure 8. Despite the difference in axial amplitudes (above), both of the transverse beampatterns are a close fit to the Bessel function of zero order: J 0 [kx], thus resembling the classical Bessel beam described by Durnin as implemented within the framework of the optical axicon [2,3]. Next, we turn to variation in the q parameter. With all other factors held the same as in Figure 6 (top), the imaginary offset is added to the z (axial) coordinate, q = 1/2 and 1, respectively within the arguments of Equation (2). It should be noted that the addition of the imaginary q offset does change the amplitudes and phases of the spherical Bessel functions at the active transducer, introducing an enhanced source strength near the origin, as shown in Figure 9. Next, we turn to variation in the q parameter. With all other factors held the same as in Figure 6 (top), the imaginary offset is added to the z (axial) coordinate, q = 1/2 and 1, respectively within the arguments of Equation (2). It should be noted that the addition of the imaginary q offset does change the amplitudes and phases of the spherical Bessel functions at the active transducer, introducing an enhanced source strength near the origin, as shown in Figure 9.
Next, we turn to variation in the q parameter. With all other factors held the same as in Figure 6 (top), the imaginary offset is added to the z (axial) coordinate, q = 1/2 and 1, respectively within the arguments of Equation (2). It should be noted that the addition of the imaginary q offset does change the amplitudes and phases of the spherical Bessel functions at the active transducer, introducing an enhanced source strength near the origin, as shown in Figure 9. Figure 9. Beampatterns with the same parameters as Figure 6 (top) with the addition of the z offset by imaginary number q = 1/2 (top) and 1 (bottom). The q factor produces stronger source amplitudes at radii greater than 9 mm (beyond the scale shown), which contribute to the emerging pattern at axial range greater than 10 mm.
Next, we examine some alternative modifications. The localization property of the summation is illustrated in Figure 10 (top), where only the sum from integer orders 150 to 250 are activated. This translates into a localized area of axial strength from 7 to 12 mm from the transducer surface. An oscillatory term is added to the coefficients and shown in Figure 10 (bottom) using Equation (7) with Nmod = 2.5. This modulates the amplitudes of the integer orders within the sum, and thereby serves as a rough equivalent to changing the angle of incidence in conical optical axicon configurations. Figure 9. Beampatterns with the same parameters as Figure 6 (top) with the addition of the z offset by imaginary number q = 1/2 (top) and 1 (bottom). The q factor produces stronger source amplitudes at radii greater than 9 mm (beyond the scale shown), which contribute to the emerging pattern at axial range greater than 10 mm.
Next, we examine some alternative modifications. The localization property of the summation is illustrated in Figure 10 (top), where only the sum from integer orders 150 to 250 are activated. This translates into a localized area of axial strength from 7 to 12 mm from the transducer surface. An oscillatory term is added to the coefficients and shown in Figure 10 (bottom) using Equation (7) with N mod = 2.5. This modulates the amplitudes of the integer orders within the sum, and thereby serves as a rough equivalent to changing the angle of incidence in conical optical axicon configurations. Finally, a single integer order of 199 is shown in Figure 11, including the offset parameter of q = 1/2. We note an extended centralized beampattern and this has some resemblance to Figure 1 of Cizmár and Dholakia [24] which demonstrates their annular source axicon. Figure 11. 5 MHz beampattern produced with only a single spherical Bessel function of integer order n = 199 and with imaginary offset parameter q = 1/2. The transducer is at the left and the axial range extends to 120 mm.

Discussion and Conclusions
The free space sums of SBH functions in Equation (2) are of particular interest along the central axis ( = 0) and in the case of a flat piston source, along the source plane ( = 2 ⁄ ). Along the central axis, the spherical harmonic functions 0 increase approximately proportional to √ with increasing n. Thus, if in Equation (2), we set the = 1, Finally, a single integer order of 199 is shown in Figure 11, including the offset parameter of q = 1/2. We note an extended centralized beampattern and this has some resemblance to Figure 1 of Cizmár and Dholakia [24] which demonstrates their annular source axicon. Finally, a single integer order of 199 is shown in Figure 11, including the offset parameter of q = 1/2. We note an extended centralized beampattern and this has some resemblance to Figure 1 of Cizmár and Dholakia [24] which demonstrates their annular source axicon. Figure 11. 5 MHz beampattern produced with only a single spherical Bessel function of integer order n = 199 and with imaginary offset parameter q = 1/2. The transducer is at the left and the axial range extends to 120 mm.

Discussion and Conclusions
The free space sums of SBH functions in Equation (2) are of particular interest along the central axis ( = 0) and in the case of a flat piston source, along the source plane ( = 2 ⁄ ). Along the central axis, the spherical harmonic functions (0) increase approximately proportional to √ with increasing n. Thus, if in Equation (2), we set the = 1, Figure 11. 5 MHz beampattern produced with only a single spherical Bessel function of integer order n = 199 and with imaginary offset parameter q = 1/2. The transducer is at the left and the axial range extends to 120 mm.

Discussion and Conclusions
The free space sums of SBH functions in Equation (2) are of particular interest along the central axis (θ = 0) and in the case of a flat piston source, along the source plane (θ = π/2). Along the central axis, the spherical harmonic functions Y o n (0) increase approximately proportional to √ n with increasing n. Thus, if in Equation (2), we set the a n = 1, the resulting terms resemble Equation (5), and then the central axis amplitude is uniform out to the limit of kz = N max , as demonstrated in Figure 6 (lower) and Figure 7 (lower). The flat piston source corresponds to θ = π/2, and is of practical interest for the specification of the source excitation. As mentioned previously, and as can be seen at the extreme edges of Figure 3, the Y o n (π/2) will oscillate between positive, zero, negative, and zero values, repeatedly with a period of 4n. This modulates the sum of product terms and generally increases the relative weights of the source for kr > N max beyond those expected from the examples of the sum of unweighted j n . In addition, the imaginary offset q tends to produce at the source plane a significantly stronger amplitude for the outer region of kr > N max . Practically, this means that the active source area needs an extended outer region and additional apodization functions may be included to restrict the active source region. Furthermore, consideration of the angular spectrum [23] and the asymptotic forms of the SBH solutions [25] indicate that fields similar to that shown in Figure 11 will have nonnegligible propagation at angles near θ = π/2, which would require an unreasonably large source radius to approximate. Utilizing a source angle slightly less than π/2, effectively a cone surface at small angle with the z axis, will help to mitigate this. We also note that the spherical Bessel functions are related to simple Bessel functions of half-integer order [22], thus Equations (1)-(6) could be rewritten in an alternative form using this substitution.
Limitations of this study are that only steady state monochromatic examples are considered. Transient behavior with broad band imaging pulses remains to be studied, and are left for future work. Additionally, practical effects of transducer element size and spacing in phased arrays will have an effect on the resulting beampatterns. These will also require further research with specific parameters.
Overall, the SBH framework presents a number of useful options for the generation of radially symmetric beampatterns, by recognizing the special properties of the sums of spherical Bessel functions and their product with spherical harmonic functions at different angles. This framework recasts some of the traditional treatments of apodization, focusing, and design specifications around the summation properties including interpolation, localization, and modulation.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.