Application of Lippmann-Schwinger Scattering Series with an Adaptive Preconditioner

1. Introduction

Modeling high-frequency acoustic wave scattering in heterogeneous media is important for applications like medical ultrasound, acoustic imaging, and materials characterization. For example, intravascular and small-animal ultrasound systems routinely operate above 20–30 MHz, where individual red blood cells (RBCs) and other microstructures strongly influence the backscattered signal [1]. In this non-Rayleigh regime (e.g., 73.2 MHz in our study), the exact shape of a scatterer – such as the biconcave geometry of an RBC – can significantly alter the acoustic field. In order to predict these fields, one must solve the Helmholtz equation for the pressure field in a fluid with embedded objects. Direct time-domain simulation (followed by Fourier analysis) is one approach, but it can be prohibitively expensive and inaccurate at high frequencies due to long runtimes and dispersion error. Similarly, naive direct solvers of the frequency-domain equation require large memory. As a result, fast frequency-domain methods have attracted attention. In particular, iterative Born-series solvers can efficiently handle large Helmholtz problems by summing successive scattering interactions. The CBS method is a preconditioned Born-series formulation that guarantees convergence even for strong scatterers [2]. Recent work has shown that a first-order, matrix-free Born-series solver (using FFTs and a special split-preconditioner) can compute pressure and velocity fields efficiently with minimal memory.

In this paper, we apply such an iterative solver to study acoustic scattering by two distinct rigid scatterers in 2D: (1) a circular disc and (2) an idealized red blood cell (RBC)–shaped biconcave plate. Key physical parameters are chosen to mimic strong scattering: the scatterer has an internal sound speed of 1650 m/s while the outside fluid has the speed of a typical fluid. We focus on monochromatic plane-wave excitation at frequencies around 73.2 MHz, so that the nondimensional size parameter ka is on the order of unity to several (where a is radius).

Our numerical approach is a frequency-domain Helmholtz solver based on the Convergent Born Series. The algorithm iteratively updates the pressure field using FFT-based convolution and a “universal split” preconditioner, avoiding any large matrix inversions. We also compare it to a reference solution (RS) to verify accuracy. Both methods operate in a matrix-free fashion: at each step the scattered field is computed via fast Fourier transforms without constructing the full system matrix. This framework is capable of handling the strong density/sound-speed contrast we impose.

Frequency dependence: As expected, higher frequency (shorter wavelength) produces more concentric interference fringes in the pressure field around the scatterer. The number of rings visible in the spatial pressure map increases with frequency, and the central beam width narrows. This is consistent with wave optics: finer features emerge when ka grows. Shape effects: The biconcave (RBC) scatterer yields a noticeably narrower central pressure lobe than the simple circular disk, even though the total number of rings is similar. This agrees with prior theoretical work showing that RBC geometry becomes important above 20 MHz. In other words, the disk and RBC of the same effective size scatter differently due to shape. Solver performance: Both the RS and CBS solvers produced nearly identical pressure fields at convergence, validating the implementation. The CBS-based solver converged reliably for our strong-contrast case. However, the number of iterations required grew rapidly at higher frequencies, reflecting the known fact that larger wavenumbers (and stronger contrasts) slow down Born-series convergence.

Overall, the FFT-based solver proved efficient and accurate for this 2D scattering problem. Despite these promising results, there are limitations to note. Our simulations assume an ideal, two-dimensional free-field (non-viscous) environment and rigid (sound-hard) obstacles. In reality, ultrasound waves in tissues would experience absorption and 3D propagation, and actual RBCs may deform or have viscoelastic properties. Future work will extend this framework to three-dimensional geometries, include absorption and more realistic boundary conditions, and consider ensembles of scatterers (e.g., blood suspensions). Experimental validation (for example using micrometer-scale scatterers in water) would also help confirm the predicted scattering patterns. Overall, this study provides a detailed computational exploration of how scatterer shape and frequency influence the acoustic pressure field and demonstrates the capability of modern Born-series solvers in handling high-frequency, high-contrast wave problems.

2. Materials and Methods

2.1. PA Model

In the PA effect, optical absorption produces a localized temperature rise in the medium. This temperature increase causes the medium to undergo thermal expansion. However, the surrounding medium restricts this expansion, resulting in the generation of a thermoelastic pressure rise. This pressure rise subsequently propagates through the medium as an acoustic wave can be expressed as [3,4,5,6],

$$\tilde{\mathrm{\Theta }}\left(\mathbf{r}\right)=\beta K▵T^0,$$

(1)

where $\tilde{\mathrm{\Theta }}\left(\mathbf{r}\right)$ denotes the thermoelastic pressure rise, $\beta$ is the volumetric thermal expansion coefficient of the medium. K (RBC, $3\times {10}^9$ Pa) [3] represents the bulk modulus, and $▵T^0$ is the temperature increase caused by optical energy absorption. When the conditions of thermal and stress confinement hold, the resulting PA pressure in a homogeneous and non-viscous acoustic medium can be formulated with thermal diffusion time by the following equation [3],

$$C_p\{v_s^2{\nabla }^2-{\partial }_t^2)\mathrm{\Theta }(\mathbf{r},t)=-\beta {\partial }_t\mathcal{H}(\mathbf{r},t),$$

(2)

where $C_p$ is the specific heat capacity at constant pressure, $v_s$ is the speed of sound inside the source, and heat deposition function $\mathcal{H}(\mathbf{r},t)$ per unit time and volume. To obtain the time-independent form, we assume a time-harmonic dependence of the acoustic field and the heat source. Therefore, the temporal variation of the field can be written as [4],

$$\begin{array}{cc}\hfill C_p\{{\nabla }^2+k_s^2\}\tilde{\mathrm{\Theta }}\left(\mathbf{r}\right)& =i\omega \mu \beta I_0,\text{within the source}\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill \{{\nabla }^2+k_f^2\}\tilde{\mathrm{\Theta }}\left(\mathbf{r}\right)& =0,\text{in the surrounding medium}\hfill \end{array}$$

(3b)

where $k_s$ and $k_f$ are the wave number within the source and ambient medium. The subscripts s and f state the source (RBC) and the ambient (plasma) medium. $\omega$ is the angular frequency of the incident laser source and $I_0$ incidendent beam intensity with optical absorption coefficient $\mu$.

2.2. Reference Solution (RS) Method

The RS method provides an analytical benchmark for validating numerical approaches. In this method the governing PA wave equation is solved analytically for a simple geometry by applying appropriate boundary conditions at the interface between the absorbing object and the surrounding fluid. For a cylindrically symmetric configuration, the acoustic field can be expressed in term of cylindrical wave function. The analytical solution is obtained by solving the Helmholtz equation inside and outside the particle while enforcing the continuity of acoustic pressure and normal component of the particle velocity at the boundary. As a result, the acoustic fields are represented by Bessel functions inside the particle and Hankel functions in the surrounding medium. Following this approach, the frequency domain solutions for the PA field inside and outside the cylindrical symmetry are given by [7],

$$\tilde{\mathrm{\Theta }}\left(\mathbf{r}\right)=\left\{\begin{array}{cc}{\displaystyle \frac{i{v}_s\mu \beta I_0}{C_p{k}_s}}[1+{\displaystyle \frac{\widehat{\rho }\widehat{c}{J}_0\left(k_{s}r\right)H_1^{\left(1\right)}\left(k_{f}a\right)}{J_1\left(k_{s}a\right)H_0^{\left(1\right)}\left(k_{f}a\right)-\widehat{\rho }\widehat{c}{J}_0\left(k_{s}a\right)H_1^{\left(1\right)}\left(k_{f}a\right)}}],\hfill & \mathrm{Inside},\hfill \\ {\displaystyle \frac{i{v}_s\mu \beta I_0}{C_p{k}_s}}\left[{\displaystyle \frac{J_1\left(k_{s}a\right)H_0^{\left(1\right)}\left(k_{f}r\right)}{J_1\left(k_{s}a\right)H_0^{\left(1\right)}\left(k_{f}a\right)-\widehat{\rho }\widehat{c}{J}_0\left(k_{s}a\right)H_1^{\left(1\right)}\left(k_{f}a\right)}}\right],\hfill & \mathrm{Outside}.\hfill \end{array}\right.$$

(4)

where $\widehat{\rho }={\rho }_s/{\rho }_f$ and $\widehat{c}=v_s/v_f$, The notations J and $H^{\left(1\right)}$ represent the Bessel and Hankel function of first kind. These expressions serve as a RS method for comparision with numerical results obtained from Lippmann-Schwinger scattering series (CBS) based methods. The detailed derivation are shown in Appendix A.

2.3. Evans–Fung RBC Model Equation

The geometric shape of a RBC is commonly described by the Evans-Fung model, which represents the chracteristic biconcave structure of the cell. In the present study, the two-dimensional cross-section of the RBC is modeled. The thickness profile of the RBc along the radial direction is given by [8],

$$z\left(r^{\prime }\right)=\pm R^{\prime }\sqrt{1-{\left({\displaystyle \frac{r^{\prime }}{R^{\prime }}}\right)}^2}\left(c_0+c_1{\left({\displaystyle \frac{r^{\prime }}{R^{\prime }}}\right)}^2+c_2{\left({\displaystyle \frac{r^{\prime }}{R^{\prime }}}\right)}^4\right)$$

(5)

where $r^{\prime }$ denotes the radial coordinate measured from the cell center and $z\left(r^{\prime }\right)$ is the half thickness of the RBC. The parameter $R^{\prime }=D/2$ corresponds to the cell radius, and D is the RBC diameter. The coeffients $c_0$, $c_1$, and $c_2$ determine the biconcave profile and are obtained from the geometric parameters of the cell. These parameters are defined as $X_1={\displaystyle \frac{d^2}{4R^{\prime 2}}}$, $X_2={\displaystyle \frac{d^4}{16R^{\prime 4}}}$, $X_3={\displaystyle \frac{h}{2R^{\prime }\sqrt{1-X_1}}}$, where d denotes the diameter of the central dimple and h represents the maximum height of the RBC rim as shown in Figure 1. The coefficients are then expressed as, $c_0={\displaystyle \frac{t^{\prime }}{2R^{\prime }}}$, $c_1={\displaystyle \frac{(X_3-c_0)-X_2{c}_2}{X_1}}$, $c_2=-{\displaystyle \frac{-c_0{X}_1-(X_1-2)(X_3-c_0)}{4X_1^2-5X_1-X_2+2X_1{X}_2}}$. For the numerical simulations, the geometric parameters of the RBC are taken as $D=7.65\mu$m, $d=0.7D$, $h=2.84\mu$m, $t=1.44\mu$m, where t denotes the minimum thickness at the center of the cell. $R^{\prime }={\displaystyle \frac{D}{2}}$. These parameters produce the typical biconcave morphology observed in healthy human RBCs.

2.4. Lippmann-Schwinger Scattering Series

The time-independent PA wave pressure field $\tilde{\mathrm{\Theta }}\left(\mathbf{r}\right)$ at point $\mathbf{r}$ with opto-acoustic source $Q\left(\mathbf{r}\right)$ satisfies the Helmholtz equation [9],

$$\begin{array}{cc}\hfill \{{\nabla }^2+(k_f^2+i\widehat{\mho })+\Phi \left(\mathbf{r}\right)\}\tilde{\mathrm{\Theta }}\left(\mathbf{r}\right)& =-Q\left(\mathbf{r}\right),\end{array}$$

(6)

where, source term $Q\left(\mathbf{r}\right)=i\omega \mu \beta I_0/C_p$, sacttering potential $\Phi \left(\mathbf{r}\right)=(k_s^2-k_f^2-i\widehat{\mho })$ and $\widehat{\mho }$ is the real positive number. After inserting the source term and the potential function, the opto-acoustic field inside and outside the source region satisfies the following governing equation,

$$\left\{{\nabla }^2+(k_f^2+i\widehat{\mho })\right\}\tilde{\mathrm{\Theta }}\left(\mathbf{r}\right)=\left\{\begin{array}{cc}i\omega \mu \beta I_0/C_p-(k_s^2-k_f^2-i\widehat{\mho })\tilde{\mathrm{\Theta }}\left(\mathbf{r}\right),\hfill & \mathrm{Inside},\hfill \\ i\widehat{\mho }\tilde{\mathrm{\Theta }}\left(\mathbf{r}\right),\hfill & \mathrm{Outside}.\hfill \end{array}\right.$$

(7)

The solution in the form of volume integral over all the space (${\mathbb{R}}^3$ or ${\mathbb{R}}^2$) of the PA wavefield $\tilde{\mathrm{\Theta }}\left(\mathbf{r}\right)$ expressed as (Morse 1953),

$$\begin{array}{c}\hfill \tilde{\mathrm{\Theta }}\left(\mathbf{r}\right)=\int g(\mathbf{r},{\mathbf{r}}_0)[\Phi \left({\mathbf{r}}_0\right)\tilde{\mathrm{\Theta }}\left({\mathbf{r}}_0\right)+Q\left({\mathbf{r}}_0\right)]d^3{\mathbf{r}}_0.\end{array}$$

(8)

The Green’s function satisfies the inhomogeneous Helmholtz equation with Dirac delta source $\delta$, which represents the reponse of the medium to a unit point excitation located at ${\mathbf{r}}_{\mathbf{0}}$.

$$\begin{array}{c}\hfill \{{\nabla }^2+(k_f^2+i\widehat{\mho })\}g(\mathbf{r},{\mathbf{r}}_0)=-\delta (\mathbf{r}-{\mathbf{r}}_0).\end{array}$$

(9)

Again, the differential Equation (6) can be transformed into an equvalent integral equation of the Lippmann-Schwinger type [10],

$$\begin{array}{c}\hfill \tilde{\mathrm{\Theta }}\left(\mathbf{r}\right)={\tilde{\mathrm{\Theta }}}^{\left(0\right)}\left(\mathbf{r}\right)+{\int }_{\Xi }{g}^{\left(0\right)}(\mathbf{r},{\mathbf{r}}_0)\Phi \left({\mathbf{r}}_0\right)\tilde{\mathrm{\Theta }}\left({\mathbf{r}}_0\right)d^3{\mathbf{r}}_0.\end{array}$$

(10)

where $g^{\left(0\right)}$ and ${\tilde{\mathrm{\Theta }}}^{\left(0\right)}$ are the Green’s function and the PA wave field in the background medium over the space $\Xi$ where $\Phi \left({\mathbf{r}}_0\right)\ne 0$. For a three dimensional, the Green’s function corresponds to a spherical wave that decays inversely with distance from the source. In a two dimensional domain, the solution is expressed in terms of the Hankel function of the first kind, which represents cylindrical wave propagation. In contrast, the one dimensional form corresponds to a pair of propagating plane waves traveling away from the source location. In In three different dimensions, the free space Helmholtz Green’s function can be written as [11],

$$g(\mathbf{r},{\mathbf{r}}_0)=\left\{\begin{array}{cc}{\displaystyle \frac{e^{i{k}_f|\mathbf{r}-{\mathbf{r}}_0|}}{4\pi |\mathbf{r}-{\mathbf{r}}_0|}},\hfill & \mathrm{in}{\mathbb{R}}^3\hfill \\ {\displaystyle \frac{i}{4}}{H}_0^{\left(1\right)}(k_f|\mathbf{r}-{\mathbf{r}}_0\left|\right),\hfill & \mathrm{in}{\mathbb{R}}^2\hfill \\ {\displaystyle \frac{i}{2k_f}}{e}^{i{k}_f|\mathbf{r}-{\mathbf{r}}_0|},\hfill & \mathrm{in}{\mathbb{R}}^1\hfill \end{array}\right.$$

(11)

For the numerical implementation Equation (10) is discretized on a Cartesian computational grid, which involves convolution sums, and get the following matrix equation,

$$\begin{array}{c}\hfill \tilde{\mathrm{\Theta }}=G^0\Phi \tilde{\mathrm{\Theta }}+G^{0}Q,\end{array}$$

(12)

After discretization, the scattering potential is represented by a diagonal matrix $\Phi$. For simplicity, $\tilde{\mathrm{\Theta }}$ is used to denote both the continuous field and its discretized counterpart. With these conditions, Equation (12) can be expressed in matrix form as,

$$\begin{array}{c}\hfill (I-G^0\Phi )\tilde{\mathrm{\Theta }}=G^{0}Q,\end{array}$$

(13)

where I the identity matrix and $G^{0}Q$ the incident field. The resulting linear system may then be solved using standard numerical techniques, for instance by performing matrix inversion,

$$\begin{array}{c}\hfill \tilde{\mathrm{\Theta }}={(I-G^0\Phi )}^{-1}{G}^{0}Q,\end{array}$$

(14)

where $G^0={\mathcal{F}}^{-1}\tilde{g}\left(\mathbf{p}\right)\mathcal{F}$, $\mathcal{F}$ and ${\mathcal{F}}^{-1}$ are the forward and inverse Fourier transform operators, respectively. The spectral representation of the free-space Green’s function can be obtained by applying the Fourier transform to the Helmholtz equation. In the wavenumber domain, the Green’s function takes a simple algebric form that depends on the magnitude of the spectral variable $\mathbf{p}$. The Fourier domain Green’s function is expressed as [2],

$$\begin{array}{c}\hfill \tilde{g}\left(\mathbf{p}\right)={\displaystyle \frac{1}{\left(\right|\mathbf{p}{|}^2-k_f^2-i\widehat{\mho })}},\end{array}$$

(15)

2.5. Lippmann-Schwinger Scattering Series with Spatially Adaptive Preconditioner

To improve the convergence behavior of the Lippmann-Schwinger scattering formulation, the governing integral equation is multiplied by a suitable preconditioning operator $\Gamma$. This leads to a modified form of the scattering equation in which the field variable can be expressed in terms of an iteration operator. After rearranging the equation, the solution can be written in an operator form involving the background Green’s operator, the scattering potential, and the source term.

$$\begin{array}{c}\hfill \Gamma \tilde{\mathrm{\Theta }}=\Gamma G^0\Phi \tilde{\mathrm{\Theta }}+\Gamma G^{0}Q.\end{array}$$

(16)

By defining an iteration operator M that contains the effect of the Green’s operator and the scattering potential, the PA field can be expressed as a recursive relation.

$$\begin{array}{c}\hfill \tilde{\mathrm{\Theta }}=M\tilde{\mathrm{\Theta }}+\Gamma G^{0}Q,\end{array}$$

(17)

where $M=\Gamma G^0\Phi -\Gamma +1$. Expanding this recursive formulation leads to an infinite Neumann-type scattering series, where the total field is represented as the summation of successive scattering contributions. Each term in the series corresponds to higher order interactions between the wave field and the scattering medium.

$$\begin{array}{c}\hfill \tilde{\mathrm{\Theta }}=\sum _{n=0}^{\infty }{M}^n\Gamma G^{0}Q,\end{array}$$

(18)

The convergence of this series strongly depends on the choice of the preconditioner. Osnabrugge preconditioner provides stable convergence for a wide range of scattering problems. A spatially adaptive preconditioning stategy is employed. In this approach, the preconditioner varies locally according to the spatial distribution of the scattering potential. This adaptive formulation improves numerical stability and accelerates the convergence of the iterative solution, particularly for heterogeneous geometries such as the RBC model considered in this study. The above equation converges only if $M<1$. Moreover, converging solutions for all strong contrast can be obtained if the preconditioner chose $\Gamma ={\displaystyle \frac{i\Phi }{(\widehat{\mho }+\kappa |\Phi \left|\right)}}$ and $\widehat{\mho }\ge Max|k_s^2-k_f^2|$ where $\widehat{\mho }=0.8k_f^2$ [2].

Algorithm 1: Lippmann-Schwinger scattering series with spatially adaptive preconditioner

Require:  Square computational domain of size $(N_y\times N_x)$, resolution $(dx=dy)$, $N_{\mathrm{mid}}$, opto-acoustic parameters ($I_0,\mu ,\beta ,C_p,v_s,v_f,k_s,k_f,\omega ,a,R,f)$, tolerance, ${10}^{-4}$, Evans-Fung parameters ($t,d,D,R^{\prime },h$) and $N_{\mathrm{throw}}$. Ensure:  PA field at stedy state ${\tilde{\mathrm{\Theta }}}^j$. 1: Compute auxiliary variables ($X_1,X_2,X_3,c_1,c_2$). 2: Define bounding box. 3: for$i=1$ to $N_{\mathrm{throw}}$ do 4:     Generate random point: 5:     $x\leftarrow -L_x+2L_x·\mathrm{rand}\left(\right)$ 6:     $y\leftarrow -L_y+2L_y·\mathrm{rand}\left(\right)$ 7:     if $\left|x\right|\le R$ then 8:         $d_1\leftarrow R^{\prime }{\left(1-x^2/R^{\prime 2}\right)}^{1/2}$ 9:         $d_2\leftarrow c_0+(c_1{x}^2/R^{\prime 2})+(c_2{x}^4/R^{\prime 4})$ 10:         if $\left|y\right|\le d_1·d_2$ then 11:            $c_{\mathrm{in}}\leftarrow c_{\mathrm{in}}+1$ 12:            Store $(x,y)$ 13:         end if 14:     end if 15: end for 16: for each sampled point $(x_t,y_t)$ do 17:     $j\leftarrow \mathrm{round}(x_t/dx)+N_{\mathrm{mid}}$ 18:     $i\leftarrow \mathrm{round}(y_t/dx)+N_{\mathrm{mid}}$ 19:     if $i,j$ inside domain then 20:         $RBC\mathcal{M}ask(i,j)\leftarrow 1$ 21:     end if 22: end for 23: Define: $\widehat{\mho },\Phi ,Q$, $\Gamma ={\displaystyle \frac{i\Phi }{(\widehat{\mho }+\kappa |\Phi \left|\right)}}$, $\tilde{g}\left(\mathbf{p}\right)$. 24: Initial field: ${\tilde{\mathrm{\Theta }}}^0\left(\mathbf{r}\right)=\Gamma {\mathcal{F}}^{-1}\left[\tilde{g}\left(\mathbf{p}\right)\mathcal{F}Q\left(\mathbf{r}\right)\right])$. 25: for$l=1$ to max iterations do 26:     Compute update: ${\tilde{\mathrm{\Theta }}}^j\left(\mathbf{r}\right)={\tilde{\mathrm{\Theta }}}^{j-1}\left(\mathbf{r}\right)-\Gamma \left({\tilde{\mathrm{\Theta }}}^{j-1}\left(\mathbf{r}\right)-{\mathcal{F}}^{-1}\left[\tilde{g}\left(\mathbf{p}\right)\mathcal{F}[\Phi \left(\mathbf{r}\right){\tilde{\mathrm{\Theta }}}^{j-1}\left(\mathbf{r}\right)+Q\left(\mathbf{r}\right)]\right]\right)$. 27:     if error < tolerance then 28:         break 29:     else 30:         ${\tilde{\mathrm{\Theta }}}^{j-1}\leftarrow \mathcal{M}as{k}_{abl}\times {\tilde{\mathrm{\Theta }}}^j$ 31:     end if 32: end for 33: Calculate steady state PA pressure field and analysis.

2.6. Numerical Experiment

The numerical experiments are carried out in four main stages. First, a RS is obtained for a circular symmetric disc. Then, the Lippmann-Schwinger scattering equation is solved iteratively. Similarly, for an irregular biconcave (RBC) geometry, both the analytical formulation and the corresponding scattering series are considered.

Figure 1. The two-dimensional cross-sectional view of the RBC is illustrated using the Evans-Fung model, and the associated geometric parameters are indicated. By characteristic biconcave discoid shape the RBC exhibits a concave profile at the center and a thicker rim near the edges, resulting in a non-uniform geometry. A Monte Carlo random method is employed to generate rondom points within boundary box and only those points that satisfy the Evans-Fung condition are retained.

Figure 1. The two-dimensional cross-sectional view of the RBC is illustrated using the Evans-Fung model, and the associated geometric parameters are indicated. By characteristic biconcave discoid shape the RBC exhibits a concave profile at the center and a thicker rim near the edges, resulting in a non-uniform geometry. A Monte Carlo random method is employed to generate rondom points within boundary box and only those points that satisfy the Evans-Fung condition are retained.

Figure 2. Shematic diagram of the computational domain for implementing the Lippmann-Schwinger scattering series with spatially adaptive preconditioner algorithm. The black region denotes the ABL whereas the yellow region indicates the source. The computational domain consistis of a $2048x2048$ grid.

Figure 2. Shematic diagram of the computational domain for implementing the Lippmann-Schwinger scattering series with spatially adaptive preconditioner algorithm. The black region denotes the ABL whereas the yellow region indicates the source. The computational domain consistis of a $2048x2048$ grid.

All simulations are performed within a square computational domain of size $2048x2048$. An ABL of 100 grid points is applied arround the domain to minimize wave wrapping effects. The spatial resolution is set to $dx=100$ nm. To further reduce artificial reflections, an exponentially decaying attenuation function is applied within the ABL region, where the parameter $\widehat{\mho }$ controls the decay near the boundaries.

$$\mathcal{M}ask{\left(\mathbf{r}\right)}_{abl}=e^{-\sqrt{\widehat{\mho }}\left|\mathbf{r}\right|},$$

(19)

For the circular disc, two different radii ($a_1=75dx$ and $a_2=90dx$) are considered, and the pressure field as well as energy density ($|\widehat{\mathrm{\Theta }}{|}^2/2{\rho }_f{v}_f^2$) are studied under the different contrast conditions. A similar analysis is carried out for the RBC geometry and its does not depend on radii. In all cases, the density inside and outside the scatterer is assumed to be identical (${\rho }_s={\rho }_f=1000$ kg/$m^3$). The background fluid sound speed is taken as $v_f=1500$ m/s.

To investigate strong scattering effects, the sound speed inside the scatterer is varied, and two contrast values are considered: $v_s=1650$ m/s and $v_s=1350$ m/s. These choices allow us to examine the robustness of the numerical solver under high contrast conditions. All opto-thermo-mechanical parameters are assumed to be unity for simplicity. For far field and near field pressure analysis, the point sensor location is placed at a distance of $350dx$ from the center of the source. All numerical simulations are implemented in matrix laboratory (MATLAB). The computations are performed on a system equpped with an NVIDIA GeForce RTX 4060 GPU (302 CUDA cores, 8 GB RAM) running Ubuntu 24.04.2 LTS, using MATLAB version R2023b. MATLAB’s built in multithreading capabilities provide additional computational speed up compared to purely sequential execution.

The PA field data are collected at four different frequencies, $7.32$ MHz, $73.2$ MHz, 146 MHz and 366 MHz. Before starting the iterations, the pressure field is initialized.

Figure 3. (a) Circular symmetry and (b) biconcave shape. A point sensor is placed at a fixed distance from the center. The blue region represents the surrounding fluid medium.

Figure 3. (a) Circular symmetry and (b) biconcave shape. A point sensor is placed at a fixed distance from the center. The blue region represents the surrounding fluid medium.

$$\begin{array}{c}\hfill {\tilde{\mathrm{\Theta }}}^0\left(\mathbf{r}\right)=\Gamma {\mathcal{F}}^{-1}\left[\tilde{g}\left(\mathbf{p}\right)\mathcal{F}Q\left(\mathbf{r}\right)\right]).\end{array}$$

(20)

The solution is then computed iteratively over successive iterations, and at each steps, the ABL function is applied to suppress boundary reflections [2].

$$\begin{array}{cc}\hfill {\tilde{\mathrm{\Theta }}}^j\left(\mathbf{r}\right)& ={\tilde{\mathrm{\Theta }}}^{j-1}\left(\mathbf{r}\right)-\Gamma \left({\tilde{\mathrm{\Theta }}}^{j-1}\left(\mathbf{r}\right)-{\mathcal{F}}^{-1}\left[\tilde{g}\left(\mathbf{p}\right)\mathcal{F}[\Phi \left(\mathbf{r}\right){\tilde{\mathrm{\Theta }}}^{j-1}\left(\mathbf{r}\right)+Q\left(\mathbf{r}\right)]\right]\right),\hfill \end{array}$$

(21)

At every iteration, the error is evaluated along the center line of the computational domain.

$$\begin{array}{c}\hfill error={\displaystyle \frac{{\sum }_{m=1}^{N_x}|{\tilde{\mathrm{\Theta }}}^j-{\tilde{\mathrm{\Theta }}}^{j-1}|}{{\sum }_{m=1}^{N_x}\left|{\tilde{\mathrm{\Theta }}}^{j-1}\right|}}.\end{array}$$

(22)

The system is considered to have reached a steady state when the error falls below a prescribed tolerance. The overall procedure is summarized in Algorithm 1.

3. Simulation Results

The Figure 4 and Figure 5 present the variation of the real part of the pressure amplitude along the x-axis. First, for the circular disc ($75dx$), the pressure amplitude is shown using both RS and CBS methods along the lines $N_x/2+1$, $N_x/2-20$, and $N_x/2+20$. In this case, the study is conducted at a frequency $f=73.2$ MHz, with the sound speed inside the source set to $v_s=1650$ m/s, while the outside density and sound speed are kept constant as before. The maximum oscillation occurs at the center of the source, and the oscillation gradually decreases towards the sides.

Figure 4a along the center line, the carrier wave shows better agreement, especially at the central maxima. Moving left and right from the center line, a slight mismatch can be observed. The black dots represent the BS method, and the red line represent the CBS method. For the biconcave case, the carrier wave is also shown in Figure 5, where the real amplitude becomes slightly lower compared to the disc. Although the number of oscillations remains the same, the width of the central maxima apperas slightly narrower.

The Figure 6 and Figure 7 whow the spatial pressure distribution and energy density distribution for both the disc and RBC cases. To better visualize the distribution map, scalling has been applied up to 250 grid points on both sides along the x and y directions from the center. From Figure 6a–h, for the same radius and contrast, increasing the frequency from left to right results in an incraese in the number of rings, making the concentric circular patterns more clearly visible. As the frequency increases, the number of iterations required for convergence also increases (28, 92, 155, 310). At lower frequencies, the pressure and the energy distributions appears relatively flat and smooth, whereas at high frequencies, clear concentric circular patterns are observed. The center region defines the maximum intensity.

Similarly, in Figure 6i–p, for a radius 90dx and sound speed, $v_s=1650$ m/s, the number of iterations required for convergence are 28, 91, 153, and 315, respectively. The same observations have also been made for the RBC case, and these results are presented in the Figure 7.

4. Discussion

In all cases the pressure field has a strong central maxima at the scatterer center, with oscillatory rings of decreasing amplitude outward. Both solver (RS and CBS) agree closely on the pressure profile (black dots vs. red line), indicating numerical consistency. The central focus is expected, plane wave impinging on a symmetric obstacle produces constructive interference along the symmetry axis. This pattern alligns with general scattering theory (Huygens-Fresnel).

Increasing the frequency (reducing wavelength) produces more visible interference rings and narrower main lobe. Physically, higher frequency waves resolve more of the scatterer’s geometry, leading to finer wavefront structure. In our results, the number of concentric rings grows (as shown from Figure 6a–h left to right), which is consistent with diffraction theory. Importantly, higher frequency also increases the number of iterations from 28 at the lowest frequency to over 300 at the highest. This behaviour is consistent with known properties of the CBS solver-higher contrast (higher k) slows convergence. In other words, as the effective wavenumber increases, each iteration propagates the solution over distances, requiring more itarations to reach convergence.

The RBC shape produces a slightly narrower central peak compared to the flat disc. Both methods yield excellent agreement in the computed fields at the center, giving confidence in the numerical implementation. The energy density distribution follows the same spatial pattern as the pressure.

5. Conclusions

The pressure field from a disk scatterer exhibits a centered high-amplitude region with concentric interference rings. These rings become more numerous and pronounced as frequency increases, reflecting the shorter wavelength; concurrently, the solver iteration count rises dramatically due to the higher effective wavenumber (consistent with Born-series convergence theory ). The RBC (biconcave) scatterer produces a qualitatively similar interference pattern, but with a modestly narrower central lobe. This confirms that the RBC’s shape influences the scattering, as predicted by prior RBC models. The two numerical methods (standard versus CBS-preconditioned) agree very well on the results, illustrating that the CBS approach is robust even under strong scattering conditions.

Acoustic intensity (energy density) maps mirror the pressure distribution, peaking at the obstacle center and forming the same ring structure, indicating conservation of energy in the scattering process. These results contribute a detailed comparison of solver methods and scatterer geometries in high-frequency acoustics. The effective validation of the CBS solver supports its use for large-scale wave problems. In practical terms, understanding how frequency and shape affect the field pattern can inform the design of ultrasonic imaging systems and targets (e.g., modeling how red blood cells interact with ultrasound).

Subsequent research could extend this framework to three-dimensional scatterers and include absorption or complex boundary conditions. For instance, modeling a suspension of RBCs or adding layered media would bring the simulation closer to biomedical scenarios. Further optimization of the solver (e.g., GPU implementation or alternative preconditioners) could improve efficiency for very high frequencies. Finally, experimental validation of these predicted patterns (especially the narrower RBC lobe) would be valuable for confirming the models.

Appendix A.1. Analytical Solution to the Wave Equation for Cylindrical Symmetry in Frequency Domain

Applying after appropriate boundary condition, the parametric form of Helmholtz wave equation inside and outside the cylinder,

$${\psi }_s\left(\mathbf{r}\right)=U{J}_0\left(k_{s}r\right)+{\displaystyle \frac{A}{k_s^2}}$$

(A1)

and

$${\psi }_f\left(\mathbf{r}\right)=S{H}_0^1\left(k_{f}r\right)$$

(A2)

Now applying the conservation equations of photoacoustics, at the point r. Both equation will satisfy at $r=a$.

$${\psi }_s{{|}_{r=a}={\psi }_f|}_{r=a}$$

and

$${\displaystyle \frac{\nabla {\psi }_s{|}_{r=a}}{{\rho }_s}}={\displaystyle \frac{\nabla {\psi }_f{|}_{r=a}}{{\rho }_f}}$$

.

now,

$${\psi }_s\left(a\right)={\psi }_f\left(a\right)$$

$$U{J}_0\left(k_{s}a\right)+{\displaystyle \frac{A}{k_s^2}}=S{H}_0^1\left(k_{f}a\right)$$

(A3)

and

$$\nabla {\psi }_s{|}_a=U{\displaystyle \frac{d}{dr}}[J_0\left(k_{s}r\right)+{\displaystyle \frac{A}{k_s^2}}]=U{k}_s{\displaystyle \frac{d}{d\left(k_{s}r\right)}}{J}_0^1\left(k_{s}r\right)$$

Now applying the properties of Bessel’s function first kind zeroth order,

$$\left[\therefore {\displaystyle \frac{d}{dr}}{J}_0^1\left(r\right)=-J_0^1\left(r\right)\right]$$

$$\nabla {\psi }_s{{|}_a=-U{J}_1^1\left(k_{s}r\right)|}_a=-U{J}_1^1\left(k_{s}a\right)$$

and now,

$$\nabla {\psi }_f{|}_a=S{k}_f{\displaystyle \frac{d}{d\left(k_{f}r\right)}}{H}_0^1\left(k_{f}r\right)=-S{k}_f{H}_0^1\left(k_{f}a\right)$$

since

$$[{\displaystyle \frac{d}{dz}}{H}_0^1\left(z\right)=0-H_1^1\left(z\right)]$$

substitute these value into second boundary condition,

$${\displaystyle \frac{U{k}_s}{{\rho }_s}}{J}_1^1\left(k_{s}a\right)={\displaystyle \frac{S{k}_f}{{\rho }_f}}{H}_1^1\left(k_{f}a\right)$$

and

$$\therefore {\displaystyle \frac{{\rho }_s}{{\rho }_f}}=\widehat{\rho }\mathrm{and}{\displaystyle \frac{k_s}{k_f}}=\widehat{c}$$

$$\therefore S=U{\displaystyle \frac{J_1^1\left(k_{s}a\right)}{\widehat{\rho }\widehat{c}{H}_1^1\left(k_{f}a\right)}}$$

from A3,

$$U{J}_0\left(k_{s}a\right)+{\displaystyle \frac{A}{k_s^2}}=U{\displaystyle \frac{J_1\left(k_{s}a\right)}{\widehat{\rho }\widehat{c}{H}_1^1\left(k_{f}a\right)}}{H}_0^1\left(k_{f}a\right)$$

$$U{J}_0\left(k_{s}a\right)\widehat{\rho }\widehat{c}{H}_1^1\left(k_{f}a\right)+{\displaystyle \frac{A\widehat{\rho }\widehat{c}{H}_1^1\left(k_{f}a\right)}{k_s^2}}=U{J}_1\left(k_{s}a\right)H_0^1\left(k_{f}a\right)$$

$${\displaystyle \frac{A\widehat{\rho }\widehat{c}{H}_1^1\left(k_{f}a\right)}{k_s^2}}=U[J_1\left(k_{s}a\right)H_0^1\left(k_{f}a\right)-J_0\left(k_{s}a\right)\widehat{\rho }\widehat{c}{H}_1^1\left(k_{f}a\right)]$$

$$U={\displaystyle \frac{A\widehat{\rho }\widehat{c}}{k_s^2}}{\displaystyle \frac{H_1^1\left(k_{f}a\right)}{J_1\left(k_{s}a\right)H_0^1\left(k_{f}a\right)-\widehat{\rho }\widehat{c}{J}_0\left(k_{s}a\right)H_1^1\left(k_{f}a\right)}}$$

From (A1) we get, the pressure field inside the cylinder,

$${\psi }_s={\displaystyle \frac{A\widehat{\rho }\widehat{c}}{k_s^2}}{\displaystyle \frac{H_1^1\left(k_{f}a\right)}{J_1\left(k_{s}a\right)H_0^1\left(k_{f}a\right)-\widehat{\rho }\widehat{c}{J}_0\left(k_{s}a\right)H_1^1\left(k_{f}a\right)}}+{\displaystyle \frac{A}{k_s^2}}$$

$${\psi }_s={\displaystyle \frac{A}{k_s^2}}[1+{\displaystyle \frac{\widehat{\rho }\widehat{c}{J}_0\left(k_{s}r\right)H_1^1\left(k_{f}a\right)}{J_1\left(k_{s}a\right)H_0^1\left(k_{f}a\right)-\widehat{\rho }\widehat{c}{J}_0\left(k_{s}a\right)H_1^1\left(k_{f}a\right)}}]$$

now for outside the cylinder,

$$\therefore S=U{\displaystyle \frac{J_1\left(k_{s}a\right)}{\widehat{\rho }\widehat{c}{H}_1^1\left(k_{f}a\right)}}={\displaystyle \frac{A\widehat{\rho }\widehat{c}}{k_s^2}}{\displaystyle \frac{H_1^1\left(k_{f}a\right)}{[J_1\left(k_{s}a\right)H_0^1\left(k_{f}a\right)-\widehat{\rho }\widehat{c}{J}_0\left(k_{s}a\right)H_1^1\left(k_{f}a\right)]}}{\displaystyle \frac{J_1\left(k_{s}a\right)}{\widehat{\rho }\widehat{c}{H}_1^1\left(k_{f}a\right)}}$$

$$S={\displaystyle \frac{A}{k_s^2}}{\displaystyle \frac{J_1\left(k_{s}a\right)}{[J_1\left(k_{s}a\right)H_0^1\left(k_{f}a\right)-\widehat{\rho }\widehat{c}{J}_0\left(k_{s}a\right)H_1^1\left(k_{f}a\right)]H_1^1\left(k_{f}a\right)}}$$

now put in (A2)

$${\psi }_f={\displaystyle \frac{A}{k_s^2}}{\displaystyle \frac{J_1\left(k_{s}a\right)H_0^1\left(k_{f}r\right)}{J_1\left(k_{s}a\right)H_0^1\left(k_{f}a\right)-\widehat{\rho }\widehat{c}{J}_0\left(k_{s}a\right)H_1^1\left(k_{f}a\right)}}$$

now,

$$\left[{\displaystyle \frac{A}{k_s^2}}={\displaystyle \frac{i\mu \beta I_0{v}_s{k}_s}{C_p{k}_s^2}}={\displaystyle \frac{i\mu beta{I}_0{v}_s}{C_p{k}_s}}={\displaystyle \frac{i\mu beta{I}_0{v}_{s}a}{C_p{k}_{s}a}}\right]$$

$$\therefore A={\displaystyle \frac{i\mu \beta I_0\omega }{C_p}},{\displaystyle \frac{\omega }{v_s}}=k_s,{\displaystyle \frac{{\rho }_s}{{\rho }_f}}=\widehat{\rho },{\displaystyle \frac{v_s}{v_f}}={\displaystyle \frac{k_f}{k_s}}=\widehat{c}$$

The PA field at a point $\mathbf{r}$ (inside the source) for an infinite cylinder of radius a becomes Diebold 1991,

$${\psi }_s={\displaystyle \frac{A}{k_s^2}}\left[1+{\displaystyle \frac{\widehat{\rho }\widehat{c}{J}_0\left(k_{s}r\right)H_1^1\left(k_{f}a\right)}{J_1\left(k_{s}a\right)H_0^1\left(k_{f}a\right)-\widehat{\rho }\widehat{c}{J}_0\left(k_{s}a\right)H_1^1\left(k_{f}a\right)}}\right].$$

(A4)

The PA field at a point $\mathbf{r}$ (outside the source) for an infinite cylinder of radius a becomes [7],

$${\psi }_{ex}={\psi }_f={\displaystyle \frac{A}{k_s^2}}\left[{\displaystyle \frac{J_1\left(k_{s}a\right)H_0^1\left(k_{f}r\right)}{J_1\left(k_{s}a\right)H_0^1\left(k_{f}a\right)-\widehat{\rho }\widehat{c}{J}_0\left(k_{s}a\right)H_1^1\left(k_{f}a\right)}}\right].$$

(A5)

Here, $A=i\mu \beta I_o{v}_s/C_p$; J and $H^1$ are the Bessel function and the Hankel function of first kind, respectively. The subscripts 0 and 1 denote the orders of each function. The dimensionless quantities are defined as, $\widehat{\rho }={\rho }_s/{\rho }_f$ and $\widehat{c}=v_s/v_f$. A schematic diagram is presented in Figure 1(a). The subscript ex denotes the exact method. Equation (A5) can also be used to yield the PA field generated by a solid circle/disk if we restrict ourselves in two dimensions. Equation (A5) in the far field can be approximated as,