Rayleigh–Sommerfeld ASM Overview

• Rayleigh–Sommerfeld ASM is a hybrid method that computes the initial planar field using a Rayleigh–Sommerfeld integral before applying FFT-based angular-spectrum propagation.
• The approach accelerates simulation time significantly in tFUS planning while maintaining focal metrics close to full-wave solvers in homogeneous media.
• It trades some model fidelity, particularly in refractive and scattering effects beyond the transfer plane, for tremendous computational efficiency.

Rayleigh–Sommerfeld Angular Spectrum Method (RS-ASM, also written RSASM) is a hybrid steady-state wave-propagation scheme in which the initial planar field is obtained from a Rayleigh–Sommerfeld surface integral and subsequent propagation is carried out by angular-spectrum evolution in the spatial-frequency domain. In the 2025 transcranial focused ultrasound study that explicitly formulates RS-ASM, the method is positioned as an accelerated alternative to full-wave simulation for rapid targeting, with the transfer plane computed from a closed-form Rayleigh–Sommerfeld integral and the deeper field propagated by FFT-based angular-spectrum steps in a homogeneous medium (Gao et al., 11 Jul 2025).

1. Definition and mathematical formulation

In the scalar approximation with time-harmonic $e^{i\omega t}$ convention, the first Rayleigh–Sommerfeld integral for the pressure at an observation point $(x,y,z)$ produced by an insonified source patch $S$ in the plane $z=0$ is reduced, under the rigid-piston (Dirichlet) or uniform-velocity (Neumann) source assumption, to the scalar form

$$p(x,y,z) \;=\; \iint_{S} p_0(x',y')\, \frac{1 - i\,k\,R}{R^2}\, e^{\,i k R} \;dS', \qquad R=\sqrt{(x-x')^2+(y-y')^2+z^2}, \qquad k=\frac{2\pi f}{c}.$$

Here $p_0(x',y')$ is the complex source pressure on the transducer surface element $dS'$ (Gao et al., 11 Jul 2025).

The angular-spectrum stage is then introduced on a planar transfer surface. Denoting the 2D Fourier transform in $(x,y)$ by

$P(k_x,k_y,0)=\mathcal{F}_{x,y}\{p(x,y,0)\},$

the field at depth $z$ is propagated according to

$(x,y,z)$0

followed by inverse transformation: $(x,y,z)$1 The defining feature of RS-ASM is therefore not the angular-spectrum propagator alone, but the specific use of the Rayleigh–Sommerfeld integral to construct the initial planar field before spectral propagation (Gao et al., 11 Jul 2025).

Within this formulation, RS-ASM is distinct from a full-wave solver. No time-stepping is needed; RS-ASM is a steady-state method with a single Fourier-domain pass. This makes it particularly suitable when the field can be decomposed into a source-to-plane stage and a plane-to-target stage without retaining full transient dynamics.

2. Discretization and algorithmic workflow

The implementation described for transcranial focused ultrasound uses an isotropic grid with $(x,y,z)$2 mm, stated as approximately $(x,y,z)$3 points per wavelength at $(x,y,z)$4 kHz. The 2D DFT/IDFT is performed on an $(x,y,z)$5 grid, zero-padded if necessary to avoid wrap-around (Gao et al., 11 Jul 2025).

The Rayleigh–Sommerfeld integral is discretized as a double sum over the transducer surface elements: $(x,y,z)$6 with

$(x,y,z)$7

This produces the field on the planar transfer surface. The subsequent propagation step is purely spectral: compute $(x,y,z)$8, form $(x,y,z)$9, $S$0, and $S$1, multiply by the transfer function $S$2, and recover the target-plane field by $S$3 (Gao et al., 11 Jul 2025).

Operationally, the workflow can be summarized in two stages. The first stage is source-to-plane projection by the Rayleigh–Sommerfeld kernel. The second stage is plane-to-target propagation by the angular spectrum. This decomposition is central to the method’s computational profile: the near-source geometry is handled explicitly through the distance kernel $S$4, while deeper propagation in a homogeneous region is treated by spectral phase accumulation.

3. Assumptions, approximations, and model fidelity

The principal approximation in RS-ASM is that the medium between the transfer plane and the target plane is treated as homogeneous water with constant $S$5, so the angular-spectrum propagation neglects refractive effects inside the skull (Gao et al., 11 Jul 2025). Skull heterogeneity is captured only by the field computed on the transfer plane. In the hybrid kW-ASM pipeline this field comes from a k-Wave simulation; in RS-ASM it comes from the closed-form Rayleigh–Sommerfeld integral.

This leads to a specific fidelity trade-off. Surface-only modeling omits multiple scattering and mode conversion in bone, and the homogeneous-propagation assumption omits attenuation and refractive heterogeneity beyond what is already encoded on the transfer surface. The paper also notes that the Rayleigh–Sommerfeld integral is $S$6 and can be expensive for very large transducer or field grids if not carefully optimized (Gao et al., 11 Jul 2025).

A common misconception is to treat RS-ASM as interchangeable with full-wave propagation. The formulation does not support that equivalence. Full-wave k-Wave retains volumetric wave physics, whereas RS-ASM explicitly trades some model fidelity for speed. The paper’s own summary describes this as a method that “trades some model fidelity (homogeneous propagation) for tremendous speed,” even while reporting small focal deviations and low normalized pressure error under the tested conditions (Gao et al., 11 Jul 2025).

4. Reported performance in MRI-derived synthetic CT tFUS planning

In the reported application, MRI-derived synthetic CT was combined with full-wave and accelerated simulation methods, including kW-ASM and RS-ASM, for CT-free transcranial focused ultrasound targeting. Across five skull models, both full-wave and hybrid pipelines using sCT demonstrated sub-millimeter targeting deviation, focal shape consistency with FWHM approximately $S$7–$S$8 mm, and less than $S$9 normalized pressure error compared to the CT-based gold standard (Gao et al., 11 Jul 2025).

The detailed focal metrics reported for comparison against the k-Wave-CT reference are as follows.

Metric	k-Wave-CT	RS-ASM-sCT
Axial focal length (FWHM in $z=0$0)	$z=0$1 mm	$z=0$2 mm
Longitudinal (axial) focal shift	—	$z=0$3 mm
Lateral focal width (FWHM in $z=0$4/$z=0$5)	$z=0$6 mm	$z=0$7 mm
Transverse (lateral) error	—	$z=0$8 mm
Normalized pressure error in focal region	—	$z=0$9

The runtime summaries in the same source require careful reading because two different summaries are presented. The abstract states that the kW-ASM and RS-ASM pipelines reduced simulation time from approximately $$p(x,y,z) \;=\; \iint_{S} p_0(x',y')\, \frac{1 - i\,k\,R}{R^2}\, e^{\,i k R} \;dS', \qquad R=\sqrt{(x-x')^2+(y-y')^2+z^2}, \qquad k=\frac{2\pi f}{c}.$$0 s to $$p(x,y,z) \;=\; \iint_{S} p_0(x',y')\, \frac{1 - i\,k\,R}{R^2}\, e^{\,i k R} \;dS', \qquad R=\sqrt{(x-x')^2+(y-y')^2+z^2}, \qquad k=\frac{2\pi f}{c}.$$1 s and $$p(x,y,z) \;=\; \iint_{S} p_0(x',y')\, \frac{1 - i\,k\,R}{R^2}\, e^{\,i k R} \;dS', \qquad R=\sqrt{(x-x')^2+(y-y')^2+z^2}, \qquad k=\frac{2\pi f}{c}.$$2 s respectively, corresponding to approximately $$p(x,y,z) \;=\; \iint_{S} p_0(x',y')\, \frac{1 - i\,k\,R}{R^2}\, e^{\,i k R} \;dS', \qquad R=\sqrt{(x-x')^2+(y-y')^2+z^2}, \qquad k=\frac{2\pi f}{c}.$$3 and $$p(x,y,z) \;=\; \iint_{S} p_0(x',y')\, \frac{1 - i\,k\,R}{R^2}\, e^{\,i k R} \;dS', \qquad R=\sqrt{(x-x')^2+(y-y')^2+z^2}, \qquad k=\frac{2\pi f}{c}.$$4 time savings. The detailed runtime comparison reports mean $$p(x,y,z) \;=\; \iint_{S} p_0(x',y')\, \frac{1 - i\,k\,R}{R^2}\, e^{\,i k R} \;dS', \qquad R=\sqrt{(x-x')^2+(y-y')^2+z^2}, \qquad k=\frac{2\pi f}{c}.$$5 SD over five skull models, deep and shallow focus, as approximately $$p(x,y,z) \;=\; \iint_{S} p_0(x',y')\, \frac{1 - i\,k\,R}{R^2}\, e^{\,i k R} \;dS', \qquad R=\sqrt{(x-x')^2+(y-y')^2+z^2}, \qquad k=\frac{2\pi f}{c}.$$6 s for full-wave k-Wave, $$p(x,y,z) \;=\; \iint_{S} p_0(x',y')\, \frac{1 - i\,k\,R}{R^2}\, e^{\,i k R} \;dS', \qquad R=\sqrt{(x-x')^2+(y-y')^2+z^2}, \qquad k=\frac{2\pi f}{c}.$$7 s for kW-ASM, and $$p(x,y,z) \;=\; \iint_{S} p_0(x',y')\, \frac{1 - i\,k\,R}{R^2}\, e^{\,i k R} \;dS', \qquad R=\sqrt{(x-x')^2+(y-y')^2+z^2}, \qquad k=\frac{2\pi f}{c}.$$8 s for RS-ASM (Gao et al., 11 Jul 2025).

Taken together, these results place RS-ASM in the category of accelerated planning methods that retain focal metrics close to the full-wave reference while substantially reducing wall-clock time. The paper’s stated conclusion is that MRI-derived sCT combined with rapid simulation techniques enables fast, accurate, and radiation-free tFUS planning (Gao et al., 11 Jul 2025).

5. Relation to discrete Rayleigh–Sommerfeld sampling theory

A broader theoretical context for RS-based propagation is provided by a sampling theorem for computational diffraction. That result begins from the Helmholtz equation

$$p(x,y,z) \;=\; \iint_{S} p_0(x',y')\, \frac{1 - i\,k\,R}{R^2}\, e^{\,i k R} \;dS', \qquad R=\sqrt{(x-x')^2+(y-y')^2+z^2}, \qquad k=\frac{2\pi f}{c}.$$9

with boundary condition $p_0(x',y')$0, and separates the field into homogeneous (propagating) and inhomogeneous (evanescent) angular-spectrum components. The inhomogeneous part $p_0(x',y')$1 satisfies

$p_0(x',y')$2

so for $p_0(x',y')$3 it is negligible and one may set $p_0(x',y')$4 (Merthe, 2013).

Under the equivalent assumption that the spatial-frequency content of $p_0(x',y')$5 is band-limited to direction cosines $p_0(x',y')$6, corresponding to numerical aperture $p_0(x',y')$7, the boundary field may be sampled on a lattice of spacing

$p_0(x',y')$8

which exactly captures all propagating plane-wave components. In that setting the continuous Rayleigh–Sommerfeld integral is transformed into a discrete sum

$p_0(x',y')$9

with a closed-form discrete kernel involving ordinary exponentials and Bessel functions (Merthe, 2013).

The same work states that the discrete formula is exact for any field composed purely of homogeneous plane waves, while neglect of the evanescent part introduces an error $dS'$0 for $dS'$1. It also states that the discrete formula has the form of a discrete convolution in $dS'$2, so one may apply 2D FFTs for all observation points at once in $dS'$3 time rather than an $dS'$4 brute-force double integral (Merthe, 2013).

This suggests a theoretical affinity between RS-ASM and sampled-diffraction reformulations of Rayleigh–Sommerfeld propagation. In both cases, the practical objective is to replace direct evaluation of highly oscillatory integrals by discrete operations better aligned with digital computation. The sampling-theorem result, however, is a statement about exact representation of propagating plane-wave content under explicit band-limiting assumptions, whereas RS-ASM in tFUS is presented as a hybrid source-to-plane plus spectral-propagation workflow.

6. Comparative methods, numerical caveats, and extensions

A separate line of work proposes sinc-series approximations for computing Rayleigh–Sommerfeld and Fresnel diffraction integrals and uses that framework to critique FFT-based angular-spectrum implementations. In that treatment, methods based on the fast Fourier transform, such as the angular spectrum method and its variants, approximate the optical fields in the source and observation planes using Fourier series, introduce artificial periodic boundary conditions, violate the preservation of bandwidth property, and show limited accuracy that decreases for longer propagation distances (Cubillos et al., 2021).

More specifically, that work states that Fourier-series approximation on a finite square imposes artificial periodicity on both the source field and the propagated field, that high-frequency modes travel transversely and reappear at the opposite boundary, and that the Fourier series of a compactly supported function is not bandlimited, so ASM violates bandwidth preservation and suffers aliasing. Its numerical results for Gaussian-beam propagation and circular-aperture diffraction report that ASM error grows with propagation distance, whereas the sinc method’s error remains unchanged and decays super-algebraically with $dS'$5 (Cubillos et al., 2021).

These observations are relevant to RS-ASM because its second stage is an angular-spectrum propagation performed by FFT2 and IFFT2. A plausible implication is that the windowing, padding, and grid-design choices in RS-ASM are not merely implementation details but can materially affect numerical fidelity in long-range or wide-angle regimes. This implication is consistent with the explicit instruction in the RS-ASM implementation to use zero-padding if necessary to avoid wrap-around (Gao et al., 11 Jul 2025).

The 2025 RS-ASM paper also lists several acceleration strategies and extensions. The Rayleigh–Sommerfeld integral can be vectorized and evaluated via graphics-card GPU kernels or multi-threaded CPU loops. The angular-spectrum steps rely on FFT2/IFFT2 and may use optimized MKL/FFTW on CPU or cuFFT on GPU. Potential extensions include split-step or multi-segment angular-spectrum methods with known inhomogeneous phase screens for skull regions, extension to broadband pulses by superposition across frequency bins, real-time implementation by precomputing and caching RS kernels, and integration into treatment navigation suites with on-the-fly adjustment of focal depth and transducer steering (Gao et al., 11 Jul 2025).

In this sense, RS-ASM occupies an intermediate position between full-wave solvers and purely spectral free-space propagators. It preserves the explicit source-to-plane Rayleigh–Sommerfeld geometry while exploiting the efficiency of homogeneous angular-spectrum propagation, and its practical value depends on how acceptable that trade-off is for the intended physical regime and error tolerance.