Reconstruction of the initial function from the solution of the fractional wave equation measured in two geometric settings

1.   Introduction

Photoacoustic imaging (PAI) is a new biomedical imaging modality that integrates the advantages of each of the underlying modalities while complementing the problems of optical and ultrasound imaging. It is a hybrid technology that combines the high-contrast and spectroscopic-based specificity of optical imaging with the high spatial resolution of ultrasound imaging ^[1,2]. PAI capitalizes on photoacoustic effects to form images of biological tissues without tissue damage. The photoacoustic effect, discovered by Alexander Graham Bell in 1880, refers to the generation of acoustic waves using thermal expansion by absorbing electromagnetic waves such as light or radio waves ^[3,4].

Photoacoustic tomography (PAT) is a PAI methodology that involves irradiating a non-ionizing pulse wave within the tissue of a given object of diagnosis to obtain a photoacoustic signal in the ultrasound range (several MHz to several tens of MHz). The photoacoustic signal is an acoustic signal generated during thermal expansion that is produced by irradiating a laser on the tissue that absorbs the irradiated laser energy. The generated photoacoustic signal is received by an ultrasonic detector placed near the object. Moreover, the spatial distribution of the pulse energy absorption contains the diagnostic information; one of our goals was to obtain this distribution from the received signal. The photoacoustic signal satisfies the wave equation, and its spatial distribution is the initial function.

Regarding the measurement procedures, it is almost impossible to judge which one is best, but the use of point detectors has been studied both mathematically and experimentally. Hence, in this article, the detector was assumed to be point-shaped with a sufficiently small dimension. At the time $t$, the detector measures the average pressure above the surface $S$ where the detectors are located. At this time, it can be a reasonable assumption that this average pressure is the value of a pressure wave $p(\cdot, t)$ for the small size of the transducer. Therefore, the data collected at a position of the detector on the surface $S$ is consistent with the restriction of $p$ to the surface $S$ ^[5].

One of the mathematical problems arising in PAT is finding the initial function from the data measured on the outside of the object and the measurement data satisfying the wave equation. According to [6, Chapter 3], solutions of fractional order differential equations describe real-life situations better than corresponding integer-order differential equations. In this study, we consider the initial value problem for the fractional wave equation ^[7,8,9] as follows:

where $-(-\Delta_{ {\bf x}})^{\frac{\alpha}{2}}$ is the Riesz space fractional derivative of order $\alpha$, as defined below, and $D^{\alpha}_{t}$ is the Caputo time-fractional derivative of order $\alpha$, i.e.,

$I^\alpha,$ $\alpha\ge 0$ is the Riemann-Liouville fractional integral

and $\Gamma(\cdot)$ is the gamma function. For $\alpha = m$, $m\in\mathbb N$, the Caputo fractional derivative coincides with the standard derivative of order $m$. For a smooth function $f$ on $\mathbb{R}^n$ with compact support, the Riesz fractional derivative ^{[10,11,12,13]} of order $\alpha$, $\alpha\ge0$ is defined as follows:

where $\mathcal F$ is the Fourier transform of a function $f$ defined by

The solution of the fractional wave equation (1.1) is given as follows:

Here

is the Mittag-Leffler function (for a more detailed explanation, see ^[14,15,16]) with

from [17, Lemma 2.23]. Because $E_2(-z^2) = \cos(z)$ for $\alpha = 2$, the solution $p_{\alpha}$ of (1.1) reduces to the solution of the wave equation. Therefore, we focus on the case of $1 < \alpha < 2$ because the case of $\alpha = 2$ has been well studied in many articles ^{[5,18,19,20,21,22,23,24,25,26,27,28,29]}. In this study, we demonstrate how to reconstruct the initial function $f$ from the measured data, which is the solution $p_\alpha, 1 < \alpha < 2$ of the fractional wave equation (1.1) restricted to a surface with point-like detectors. To the best of our knowledge, such a PAT model has been studied here for the first time.

2.   Preliminaries

Here, we consider two geometries where point-like detectors are located: spherical and hyperplanar geometries. As their names imply, in each case, detectors are located on the unit sphere and hyperplane, respectively (see Figure 1). Our goal was to reconstruct the initial function $f$ from the measurement data, that is, the solution of (1.1) on two geometries.

In the spherical geometry, the solution $p_{\alpha}$ of (1) is measured on the unit sphere $\mathbb{S}^{n-1}$ in $\mathbb{R}^{n}$. Let the wave forward operator $\mathcal W_S$ be defined as $\mathcal{W}_{S}f(\boldsymbol{\theta}, t; \alpha) = p_{\alpha}(\boldsymbol{\theta}, t),$ $(\boldsymbol{\theta}, t)\in\mathbb{S}^{n-1}\times[0, \infty),$ where $f$ is an initial function of (1).

Similar to the spherical geometry, the solution $p_{\alpha}$ of (1) is measured on the hyperplane $\{ {\bf x} = (\boldsymbol{x}_{*}, x_n)\in \mathbb{R}^{n} : x_n = 0, \ \boldsymbol{x}_{*} \in \mathbb{R}^{n-1} \}.$ Similarly, let the wave forward operator $\mathcal{W}_H$ be defined as $\mathcal{W}_{H}f(\boldsymbol{u}, t; \alpha) = p_{\alpha}(\boldsymbol{u}, t)$, $(\boldsymbol{u}, t) \in \mathbb{R}^{n-1} \times [0, \infty),$ where $f$ is an initial function of (1).

For both geometries, the Mellin transform is essential to finding the initial function $f$ from measurement data. Moreover, spherical harmonics are employed in the spherical geometry. The remainder of this section is devoted to introducing the Mellin transform and spherical harmonics.

Regarding the Mellin transform, the majority of the Mellin transform is derived from [34,p. 79$\sim$90]. Let $f$ be a locally integrable function defined on $(0, \infty).$ The Mellin transform of $f$ is defined as

when the integral converges. Suppose that

where $O$ is the Big $O$ notation, $\epsilon > 0$, and $a < b$. The integral (2.1) converges absolutely and defines an analytic function in the strip $a < \mathrm{Re}(s) < b.$ Furthermore, its inverse transform is given by

Then, $f$ can be recovered from its Mellin transform $\mathcal Mf$ by using the inverse Mellin transform. The Mellin transform satisfies the property

where the convolution is defined by

Regarding the spherical harmonics, let $\boldsymbol{Y}_{lk}$ denote the spherical harmonics^[30,31] that form a complete orthonormal system in $L^2(\mathbb{S}^{n-1})$. Then, $f$ can be expanded in the spherical harmonics as

where $N(n, l) = (2l+n-2)(n+l-3)!/(l!(n-2)!)$ for $l\in\mathbb{N}$ and $N(n, 0) = 1$. Moreover, we use the spherical harmonics expansions of the measurement data $\mathcal W_{S}f(\boldsymbol{\theta}, t; \alpha)$ and the Fourier transform $\mathcal{F}f({\boldsymbol{\xi}})$ of the initial function $f$ with ${\boldsymbol{\xi}} = \lambda_{ {\boldsymbol{\xi}}} \boldsymbol{\omega}_{ {\boldsymbol{\xi}}}$, as follows:

and

3.   Inversion procedure

To recover the initial function, we assume that the point-like detectors are located on the unit sphere and hyperplane. Below, we provide a method to obtain the initial function $f$ from the solution of the fractional wave equation measured on two geometries.

3.1.   Spherical geometry

This section demonstrates how to obtain the initial function $f$ from $\mathcal{W}_{S}f.$ From (1.2), the measurement data $\mathcal{W}_{S}f$ are given as

First, we consider a relation between $(\mathcal{W}_{S}f)_{lk}$ and $(\mathcal{F}f)_{lk}$.

Lemma 1. For $f\in C^\infty(\mathbb{R}^n)$ with compact support, we have

where $J_{\nu}(\cdot)$ is the Bessel function of the first kind of order $\nu$.

Proof. Changing the variables ${\boldsymbol{\xi}}\rightarrow \lambda_{ {\boldsymbol{\xi}}} \boldsymbol{\omega}_{ {\boldsymbol{\xi}}}$ in (3.1), we write the measurement data as

where in the last line, we used the Funk-Hecke theorem [30, (3.19) in Chapter 7]:

A comparison with (2.3) completes our proof.

Now we present the main theorem:

Theorem 2. For $f\in C^\infty(\mathbb{R}^n)$ with compact support, we have

where

Proof. By changing the variables $\lambda_{ {\boldsymbol{\xi}}} \rightarrow \tilde{\lambda_{ {\boldsymbol{\xi}}}}^{-1},$ (3.2) can be represented as

where

To check that the Mellin transform of $(\mathcal{W}_{S}f)_{lk}(\cdot; \alpha)$ in (3.5) is well-defined, it suffices to check that the Mellin transforms of $F_{lk}$ and $E$ are well-defined, respectively. Let us consider the Mellin transform of $F_{lk}$. Notice that

because $J_{\nu}(\tilde{\rho}) = O(\tilde{\rho}^{\nu})$ as $\tilde{\rho} \to 0^{+}$ ^[32]. Therefore, $\mathcal M \left(F_{lk}\right)(s)$ is well-defined for $\mathrm{Re}(s) < l+n-\epsilon$ for any $\epsilon > 0$. Next, we consider the Mellin transform of $E.$ Taking the Mellin transform of $E$, we obtain the following formula (see [33,Lemma 9.1] or [15,(2.18)]): for $0 < \mathrm{Re}(s) < \alpha$

Hence, we have

where in the third equality, we applied the Euler's reflection formula $\Gamma(p)\Gamma(1-p) = \dfrac{\pi}{\sin(\pi p)}.$ Thus, the Mellin transform of (3.5) is well-defined for $0 < \mathrm{Re}(s) < \alpha$. Taking the Mellin transforms on both sides of (3.5), we have

where in the second equality, we used (3.7).

Now, taking the inverse Mellin transform of $\mathcal{M} \left[(\mathcal{W}_{S}f)_{lk}(\cdot; \alpha)\right](s)$, we can reconstruct $F_{lk}$ and $(\mathcal Ff)_{lk}$.

Corollary 3. For $f\in C^\infty(\mathbb{R}^n)$ with compact support, we reconstruct $f_{lk}$ from $(\mathcal{W}_{S}f)_{lk}$ by recovering the $F_{lk}:$

Thus far, we have considered the measurement data $\mathcal{W}_{S}f$. Note that our approach can be applied to the direction dependent measurement data from the model in ^[29] as well.

Remark 4. For $f\in C^\infty(\mathbb{R}^n)$ with compact support, let

be the direction dependent measurement data modeled as described in [29,(1.2)], where $\boldsymbol{\theta}\cdot\nabla_ {\bf x}\mathcal{W}_{S}f$ is the normal derivative of $\mathcal{W}_{S}f$ and $c_{1}$ and $c_{2}\in \mathbb{R}$ are constants. Using $(2.3)$, $(3.3)$, and the Bessel function identity $\frac{{\rm d}}{{\rm d}\lambda}\left[\lambda^{-\nu}J_{\nu}(\lambda)\right] = -\lambda^{-\nu}J_{\nu+1}(\lambda)$ (see, [32,(5.13) on pp. 133]), we obtain $g_{lk}$:

where we used $(3.6)$ and

By taking the Mellin transform on both sides of $(3.8)$, and because $\mathcal M (g_{lk})$ is well-defined for $0 < \mathrm{Re}(s) < \alpha$, we have $\mathcal M(F_{lk}):$

Moreover, using the inverse Mellin transform of $\mathcal M(F_{lk})$, we can recover $F_{lk}$, $\mathcal Ff_{lk}$, and $f$ from the Mellin transform $\mathcal M(F_{lk})$.

3.2.   Hyperplanar geometry

Similar to the previous subsection, Section 3.1, we show that $f$ can be determined from $\mathcal W_H f.$ From (1.2), the measurement data $\mathcal W_H f$ are given as follows:

for $(\boldsymbol{u}, t) \in \mathbb{R}^{n-1} \times [0, \infty).$ If $f$ is odd in $x_n$, then $\mathcal W_Hf(\boldsymbol{u}, t; \alpha) = 0$. Thus we assume that $f$ is even in $x_n$. First, we analyze the analog of the Fourier slice theorem:

Lemma 5. For $f\in C^\infty(\mathbb{R}^n)$ with compact support and that is even in $x_n$, we have

The Lemma for $\alpha = 2$ has already been studied in ^[18,23,25].

Proof. Taking the $n-1$-dimensional Fourier transform of $\mathcal W_Hf$ defined in (3.9) with respect to $\boldsymbol{u}$, we have

where in the second line, we used the evenness of $\mathcal Ff$ and $E_{\alpha}$ with respect to the last variable $\xi_{n}$, and in the last line, we changed the variables $\left\vert (\boldsymbol{\eta_{*}}, \xi_{n}) \right\vert \rightarrow \lambda.$

Now we present the main theorem:

Theorem 6. For $f\in C^\infty(\mathbb{R}^n)$ with compact support and that is even in $x_n$, we have

where

Proof. By changing the variables $\lambda\rightarrow \tilde{\lambda}^{-1},$ (3.10) can be represented as

where in the second line, we used the convolution (2.2) and (3.6). To demonstrate that the Mellin transform of $\mathcal F_{\boldsymbol{u}}(\mathcal W_Hf)$ defined is well-defined, we only need to check that the Mellin transforms of $F_{\boldsymbol{\eta_{*}}}$ are well-defined because, by Theorem 2 $\mathcal M (E)$ is well-defined for $0 < \mathrm{Re}(s) < \alpha$. Notice that

Therefore, $\mathcal M \left(F_{\boldsymbol{\eta_{*}}}\right)(s)$ is well-defined for any $s\in \mathbb C$ and the Mellin transform of $\mathcal F_{\boldsymbol{u}}(\mathcal W_Hf)(s)$ is well-defined for $0 < \mathrm{Re}(s) < 2$. Taking the Mellin transform, we have

where in the second equality, we used (3.7).

Again, taking the inverse Mellin transform of $\mathcal{M}(F_{\boldsymbol{\eta_{*}}})(s)$, we reconstruct $F_{lk}$ and $(\mathcal Ff)_{lk}$.

Corollary 7. For $f\in C^\infty(\mathbb{R}^n)$ with compact support and that is even in $x_n$, we reconstruct $\mathcal{F}f$ from $\mathcal{W}_{H}f$ by recovering the $F_{\boldsymbol{\eta_{*}}}$; accordingly, for $\boldsymbol\eta = (\boldsymbol{\eta_{*}}, \eta_n)\in\mathbb R^{n-1}\times\mathbb R$,

4.   Conclusions

Recovering the initial function $f$ from the solutions of the wave equation on some surface surrounding the object is crucial for the recently developed PAT methodology. In this study, we first investigated one mathematical problem of PAT by using the fractional wave equation in order to provide a way for to reconstruct $f$ from fractional wave equation solutions restricted to the sphere and hyperplane.

We summarize both cases as follows:

For the spherical case, we can recover $f$ from $\mathcal{W}_{S}f$ through the following steps:

$1)$ Find $(\mathcal{W}_{S}f)_{lk}$ using the spherical harmonics (see Lemma 1).

$2)$Take the Mellin transform of $(\mathcal{W}_{S}f)_{lk}.$

$3)$From Theorem 2, we obtain $\mathcal{M}(F_{lk})$ from $\mathcal{M}\left[(\mathcal{W}_{S}f)_{lk}\right]$.

$4)$ Taking the inverse Mellin transform, we recover $F_{lk}$ from the Mellin transform $\mathcal{M}(F_{lk})$ (see Corollary 3).

$5)$ Next, we find $(\mathcal{F}f)_{lk}$ from $F_{lk}$ and finally get $f$.

For the hyperplane case, we can recover $f$ from $\mathcal{W}_{H}f$ through the following steps:

$1)$ Take the $n-1$-dimensional Fourier transform of $\mathcal{W}_{H}f$ to get $\mathcal F_{\boldsymbol{u}} (\mathcal W_Hf)$ (see Lemma 5).

$2)$ Take the Mellin transform of $\mathcal{F}_{\boldsymbol{u}}(\mathcal{W}_{H}f).$

$3)$ Using Theorem 6, we find $\mathcal{M}(F_{\boldsymbol{\eta_{*}}})$ from $\mathcal{M}\left[\mathcal{F}_{\boldsymbol{u}}(\mathcal W_Hf)\right]$.

$4)$ Taking the inverse Mellin transform, we recover $F_{\boldsymbol{\eta_{*}}}$ from the Mellin transform $\mathcal M (F_{\boldsymbol{\eta_{*}}})$ (see Corollary 7).

$5)$ Next, we find $\mathcal Ff$ from $F_{\boldsymbol{\eta_{*}}}$ and finally get $f$.

Acknowledgements

The authors are thankful to the referees for their multiple suggestions that helped to improve this paper. We would like to thank Editage (www.editage.co.kr) for English language editing. This study was supported by a National Research Foundation of Korea grant (MSIP) and the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the South Korean government (NRF-2022R1C1C1003464, NRF-2020R1F1A1A01065912, NRF-2020R1A4A1018190). This work was conducted as H. Jung's Master's thesis.

Conflict of interest

This study does not have any conflicts of interest.