Efficient Evaluation of Gravitational Lensing Amplification Factors: A Deep Learning Framework

We adopt the standard geometric configuration for gravitational lensing under the thin-lens approximation, as illustrated schematically in Figure 1. A GW signal emitted by a source S at a true angular position β is deflected by a lens mass before reaching the observer. Due to this deflection, the observer perceives the signal coming from an apparent angular position θ. As shown in Figure 1, these angular positions correspond to physical coordinate vectors in the source plane and (the impact parameter) in the lens plane, respectively. The system is characterized by the angular diameter distances from the observer to the lens D_L, to the source D_S, and between the lens and the source D_LS.

Zoom In Zoom Out Reset image size

Following P. Christian et al. (2018), we define the normalized coordinates in the lens plane , source plane , and observer plane by the Einstein radius ξ₀,

where , , and represent the physical position vectors in the lens, source, and observer planes, respectively. The Einstein radius ξ₀ serves as the natural length scale.

The amplification factor F(ω, y), defined as the ratio of the lensed and unlensed GW amplitudes, is computed by solving the Fresnel–Kirchhoff diffraction integral in the thin-lens approximation,

where is the time delay function. And the dimensionless frequency ω, related to the physical frequency f, is given by

where is the angular frequency of the GW.

The amplification factor for a PML can be derived by solving the diffraction integral, leading to the following closed-form expression for F(ω, y)​​​​​​ (R. Takahashi & T. Nakamura 2003; P. Christian et al. 2018),

where ₁F₁(a, b, z) is the confluent hypergeometric function of the first kind and Γ is the Gamma function. The term ϕ_m(y) represents the minimum time delay and is defined by

Here x_m is the position in the lens plane where the time delay is minimized, and y is the dimensionless source-lens separation in units of the Einstein radius.

For a PML, ω in Equation (3) can be written as

where M_L is the mass of point lens and z_L is the lens redshift. This frequency is widely used as the standard parameter to distinguish the wave-optics and geometric-optics regimes. If ω ≫ 1, under this situation the geometric-optics approximation is applicable. If the dimensionless frequency w does not satisfy this condition, the full diffraction integral incorporating wave-optics effects should be evaluated. Stellar-mass black holes can serve as representative lens objects in the PML model. For a lens mass of M_L = 60 M_⊙ at a redshift z_L = 1, and considering the most sensitive frequency band (10–1000 Hz) of current ground-based GW detectors, the corresponding dimensionless frequency range computed using Equation (6) is w ∈ [0.15,  14.85]. In this regime, the geometric-optics approximation may not be sufficient. For PML with smaller masses or located at lower redshifts, the upper bound of w becomes even smaller, further emphasizing the necessity of computing the full amplification factor that includes wave-optics effects.

These expressions provides a detailed description of the amplification factor for lensed GWs and serves as the reference model for this study. The complex nature of F(ω, y), particularly its highly oscillatory behavior at large y, highlights the need for accurate numerical methods or machine learning–based approximations to efficiently model this function.

The amplification factor F(ω, y) depends critically on two dimensionless inputs: frequency ω and impact parameter y. In the low-frequency limit (ω < 1), diffraction effects dominate and F approaches unity with a smooth envelope. As ω increases, interference fringes appear and become more closely spaced, reflecting rapid phase accumulation between multiple propagation paths. Similarly, small y values (near alignment) produce strong magnification peaks at low frequency, while larger y values weaken the overall amplitude but generate high-frequency, dense oscillations as path-length differences grow.

To visualize the characteristic behavior of the amplification factor in the wave-optics regime, we present a representative case in Figure 2 with a dimensionless impact parameter y = 1.0. In this figure, interference fringes are clearly resolved. We fix the lens redshift at z_L = 0.25 and consider a lens mass of M_L = 30 M_⊙. According to Equation (6), the dimensionless frequency range will be ω ∈ [0.006, 10], mapping to a physical frequency band of f = [1.4 Hz, 2.2 kHz], which effectively covers the sensitivity window of ground-based GW detectors. As illustrated in Figure 2, the resulting amplification factor ∣F(ω, y)∣ exhibits highly nonstationary behavior: The oscillation frequency increases with ω, reflecting the rapid phase accumulation between interfering paths, while the amplitude modulation follows a slow envelope. The pattern is neither strictly periodic nor purely sinusoidal, and each extremum features a sharp peak. This multiscale structure, which combines a broad diffractive envelope with dense, high-frequency fringes, poses a severe challenge for low-order basis expansions or naive interpolation. This complexity motivates our adoption of a learning model capable of capturing both global trends and fine-grained spectral features.

Zoom In Zoom Out Reset image size

Following standard wave-optics treatments (X. Guo & Y. Lu 2020), the amplification factor for an axisymmetric PML can be written as

where J₀ is the Bessel function of the first kind. For the definitions of other variables, including the dimensionless frequency ω, the impact parameter y, and the lens potential, we refer the reader to Section 2. Substituting z = x²/2, the diffraction integral transforms to

and .

3.1.1. Integral Mean Method

We approximate the improper integral (8) by its (C, 1) Cesàro mean,

which damps late-time oscillations and converges in the mean as b → ∞.

3.1.2. Asymptotic Expansion Method

Split the domain at a finite cutoff b and write the integrand as e^iωzf(z), where

Repeated integration by parts on the tail [b, ∞) leads to the asymptotic-expansion estimator

where N is the truncation order and f ⁽ⁿ⁻¹⁾(b) denotes the (n − 1)th derivative of f at z = b.

3.1.3. Levin’s Method

Levin’s method computes oscillatory integrals on a finite interval,

where f(z) is a nonoscillatory amplitude and q(z) is a smooth phase. Choose n smooth, linearly independent basis functions and collocation nodes . The coefficients are determined by the collocation system

The collocated antiderivative yields the endpoint expression

We approximate the finite-interval integral by I(b) ≈ I_n(b). In our setting the integration variable is z = x²/2 with domain [0, ∞), so we set a = 0 and q(z) = ωz. Sending b → ∞ and assuming the endpoint term at b vanishes, we obtain

3.1.4. Zero-points Integral Method

Define the partial integral

Let j_k be the kth positive root of J₀. Under z = x²/2,

Evaluating (16) at b = z_k gives

The improper integral is

A (C, 1) Cesàro average often accelerates convergence:

3.1.5. Zero-points Asymptotic Expansion Method

With the zero-points z_k from Equation (17) and the partial integral I(z_k) from Equation (19), the zero-points asymptotic expansion (ZPAE) is

For accelerated convergence, one may apply a (C, 1) Cesàro average to the asymptotic-expansion estimator I_AE(b) in Equation (11) evaluated at b = z_k:

The wave-optics amplification factor F(ω, y) is characterized by intense oscillatory behavior, where dense interference fringes modulate a diffractive envelope. Numerically, direct evaluation of Equation (4) is hindered by high computational costs and instability arising from the cancellation of special-function terms in destructive-interference intervals. From a deep learning perspective, standard MLPs are also ill-suited for this task, as they suffer from spectral bias—a phenomenon where networks prioritize low-frequency features and fail to resolve rapid variations (N. Rahaman et al. 2018). In fact, our empirical tests also confirmed that standard MLPs completely fail to handle this highly oscillatory task. To overcome this limitation, we implement a coordinate-based neural network that integrates Fourier feature mappings (M. Tancik et al. 2020) with SIRENs (see Figure 3).

Zoom In Zoom Out Reset image size

3.2.1. Fourier Feature Mapping

Before entering the network, the low-dimensional input coordinates x = (ω, y) are projected into a higher-dimensional feature space using a Fourier feature mapping . This mapping enables the network to learn high-frequency functions more efficiently by controlling the spectral decay of the neural tangent kernel (NTK; M. Tancik et al. 2020). Specifically, we employ a Gaussian random Fourier feature mapping defined as

where is a fixed weight matrix containing m sampled frequency vectors, with entries sampled from a Gaussian distribution . This projection maps the input onto a periodic manifold, enriching the input representation with high-frequency components prior to the learnable layers.

3.2.2. Sinusoidal Representation Network

The embedded features z₀ = γ(x) serve as the input to the SIREN backbone. Unlike standard MLPs, SIREN employs a periodic sine activation function for every hidden layer. The network operation is defined recursively as

followed by a final linear transformation z_L = W_Lz_L−1 + b_L to output the real and imaginary parts of the amplification factor. Here ω₀ is a hyperparameter that dictates the frequency of the activation function.

The periodic nature of the activation function provides a strong structural consistency suited for wave physics, allowing the network to represent derivatives and fine-grained structures accurately. To ensure stable signal propagation, we adopt the specialized initialization scheme proposed by  V. Sitzmann et al. (2020). For the first hidden layer (accepting the Fourier features), weights are initialized from a uniform distribution,

where d_in is the dimensionality of the input Fourier features (i.e., 2m). For subsequent layers, weights are initialized to preserve the activation variance,

where n_h denotes the number of neurons in the hidden layers. This initialization strategy prevents the gradients from vanishing or exploding, enabling the training of deep architectures capable of resolving the complex interference patterns of F(ω, y).

We employ the NTK (A. Jacot et al. 2018) to analyze the training dynamics, characterizing the convergence rate of neural networks in the frequency domain. First, consider the spectral properties of the target function. As shown in Equation (7), the amplification factor F(ω, y) is composed of a rapidly oscillating exponential prefactor and a Bessel-weighted diffraction integral. In the wave-optics regime (large ω and y), the local oscillation frequency is determined by the gradient of the phase terms. Specifically, the characteristic bandwidth ν_target is proportional to the product of the physical parameters:

This implies that the spectral support of F(ω, y) extends to high frequencies, requiring the model to resolve fine-scale variations.

Next, consider the learning process. For a network f(x;θ) trained by gradient descent, the training dynamics in the infinite-width limit decouple into independent eigenmodes of the NTK. The error decay for a specific eigenmode is governed by its eigenvalue λ_ν,

where is the error component corresponding to the spectral frequency ν of the target function’s decomposition (i.e., the spatial frequency in the parameter space (ω, y)), η is the learning rate, and λ_ν determines the convergence rate. Note that ν here characterizes the rate of variation of the function F with respect to its inputs, distinct from the physical GW frequency ω, which serves as an input coordinate.

For standard MLPs with rectified linear unit (ReLU) activation, N. Rahaman et al. (2018) demonstrated that the NTK eigenvalues decay rapidly according to a power law, typically . This induces a “spectral bias”: For the high-frequency modes present in the wave-optics regime (large ν_target), the convergence speed becomes negligibly small. Consequently, ReLU networks act as low-pass filters, effectively smoothing out the essential interference fringes.

In contrast, SIREN employs periodic activations . The NTK spectrum of SIREN is not concentrated at zero but exhibits tunable bandwidth controlled by the initialization frequency ω₀. By setting ω₀ to match the characteristic scale of the diffraction integral, we ensure that λ_ν remains significant across the target spectrum. Structurally, the composition of sinusoidal layers synthesizes a function:

This form is structurally analogous to the oscillatory nature of Equation (7). The correspondence between the network’s spectral properties and the physical model ensures the accurate representation of fine-grained diffraction patterns.

The choice of the input domains for ω and y is dictated by both the physical behavior of the amplification factor F(ω, y) and numerical stability considerations.

We select three representative masses of the PML M_L = {10, 30, 50} M_⊙, which target the stellar-mass black hole population routinely inferred in LVK catalogs. These values are used solely to determine the physically relevant boundary limits of the dimensionless frequency ω in the LVK sensitivity band, rather than to fix the lens mass during training. Based on the optimal sensitivity band of current GW detectors (advanced LVK), we select a frequency range of f ∈ [1, 1000] Hz, corresponding to the typical GW frequencies emitted by CBCs.

Throughout, we convert this physical range to the dimensionless frequency ω using Equation (6). Taking the union of these resulting intervals yields our final training domain,

where and . The model is trained directly on F(ω, y), so the dependence on M_L is fully encoded through ω, and no mass-specific training is required. This range is critical as it covers the transition from the diffraction-dominated regime to the dense oscillatory fringes of the wave-optics limit.

It is worth noting that due to the scale invariance of the dimensionless formulation, this domain is astrophysically universal. Focusing on the prime LISA frequency band of f ∈ [10⁻⁴, 10⁻¹] Hz (P. Amaro-Seoane et al. 2017), we consider the redshifted lens mass M_Lz = (1 + z_L)M_L as the primary physical parameter. In this context, our training domain of ω ∈ [0.0067, 43.93] is astrophysically versatile, as it effectively encompasses the diffraction patterns produced by M_Lz ranging from ∼10³ M_⊙ to ∼10⁷ M_⊙. Thus our selected parameter space provides a rigorous testbed applicable to a wide range of supermassive black hole lensing scenarios.

For the impact parameter y, we set

with and . The lower bound avoids numerical instability in the analytic formula as y → 0, while the upper bound extends to the weak-lensing regime where F ≈ 1. Our primary focus for detection prospects lies in the high-magnification region y ≤ 3.

The analytic amplification factor F(ω, y) (Equation (4)) was evaluated on a dense Cartesian grid covering the selected domain. Special-function computations were carried out using the mpmath library with global precision set to 50 decimal places. The resulting ground-truth labels were validated against adaptive quadrature of the Fresnel diffraction integral and geometric-optics limits, achieving agreement within 10⁻¹⁵.

Prior to network ingestion, the input coordinates (ω, y) are first normalized to the range [−1, 1] via linear scaling:

Let x = (ω_norm,  y_norm). As detailed in the model architecture (Section 3.2), we enrich the input representation using a Gaussian Fourier feature mapping. In this implementation, we set the mapping dimension to m = 20 and the Gaussian scale to σ = 10.0. The resulting projected features are concatenated with the normalized inputs x to form the final input vector passed to the SIREN backbone.

The amplification factor exhibits highly nonstationary behavior: Smooth variations at low y and rapid, dense oscillations at high y. This behavior poses a multiscale challenge for unified modeling. Inspired by the “Frequency Principle” (Z.-Q. J. Xu et al. 2020), which states that neural networks tend to fit low-frequency components first, we adopt a domain decomposition strategy to reduce the spectral dynamic range of each subtask. Analogous to the adaptive integration intervals employed by  X. Guo & Y. Lu (2020) for numerical integration, we decompose the y-domain into four overlapping intervals defined by breakpoints {0.2, 1.0, 3.0, 6.0, 10.0}.

We employ a domain decomposition based on the interference characteristics: and cover the strong-lensing and high-contrast fringe regimes relevant to CBCs, while and handle the high-frequency saturation and weak-lensing tails. This allows us to tune the network bandwidth according to local oscillatory characteristics, preventing optimization stagnation common in single networks facing mixed-frequency tasks.

Within each interval, we employ an adaptive tensor-product grid for sampling. High-frequency intervals utilize finer grids (up to N_ω = 20,000 points) to resolve dense fringes, while smoother regions use coarser sampling (see Table 1).

Table 1. Adaptive Sampling Configuration for Each Interval . and Denote Grid Points in y and ω

Interval	y-range			Total Samples
	[0.2,  1.0)	200	5000	1.0 × 106
	[1.0,  3.0)	400	5000	2.0 × 106
	[3.0,  6.0)	600	10,000	6.0 × 106
	[6.0,  10.0]	800	20,000	1.6 × 107

Download table as:  ASCIITypeset image

This modular design extends beyond training to the inference phase. By establishing a direct mapping in the parameter space (ω, y), the system automatically routes any query to the appropriate SIREN subnetwork based on the input value of y. Unlike traditional numerical integration, where cost scales with oscillation frequency due to sampling requirements, this routing mechanism decouples inference time from ω. This theoretically ensures a constant computational complexity, a property we empirically verify in Section 5.3.2.

All models are implemented in PyTorch 2.0 and trained on NVIDIA A100 GPUs using mixed precision. We use the Adam optimizer with a plateau-based scheduler. Separate loss functions are minimized for the real and imaginary components of the amplification factor.

To ensure numerical continuity across the domain decomposition boundaries (i.e., at y = 1.0, 3.0, 6.0), we employed an overlapping training strategy. Specifically, the boundary values are explicitly included in the training datasets of both adjacent intervals (e.g., y = 1.0 is present in both and ). Since both submodels converge to the same high-precision analytical ground truth at these interfaces, the discrepancy between models at the boundaries is minimized to the order of the intrinsic validation error (∼10⁻³). This ensures that the piecewise inference remains sufficiently smooth for likelihood-based sampling methods, preventing artificial jumps in the parameter estimation posterior. Key hyperparameters, selected via grid search, are summarized in Table 2.

We quantify the model performance using the complex relative error evaluated on uniformly sampled random off-grid points (ω, y):

As summarized in Table 3, our SIREN model attains per-mille-level accuracy across the domain. The overall mean relative error is 1.10 × 10⁻³, with a median of 8.30 × 10⁻⁴. The 99th percentile error is 4.94 × 10⁻³, indicating well-controlled outliers.

Table 3. Summary of the Model’s Relative-error Performance

Metric	Overall	[0.2, 1.0]	[1.0, 3.0]	[3.0, 6.0]	[6.0, 10.0]
Mean rel. error	1.10 × 10−3	1.83 × 10−3	1.58 × 10−3	8.27 × 10−4	8.75 × 10−4
Median rel. error	8.30 × 10−4	1.47 × 10−3	1.22 × 10−3	6.94 × 10−4	7.28 × 10−4
95th percentile error	3.02 × 10−3	4.73 × 10−3	4.22 × 10−3	2.05 × 10−3	2.18 × 10−3
99th percentile error	4.94 × 10−3	7.40 × 10−3	6.26 × 10−3	2.74 × 10−3	2.99 × 10−3
Maximum error	3.23 × 10−2	2.05 × 10−2	3.23 × 10−2	4.84 × 10−3	1.26 × 10−2

Note. Evaluated on random off-grid samples, reported metrics include mean, median, 95th/99th percentiles, and maximum error across different y-intervals.

Download table as:  ASCIITypeset image

The spatial distribution of errors is visualized in the heatmaps of Figure 5. We observe distinct error morphologies corresponding to different physical regimes. At small impact parameters (y ≲ 1), the relative errors are slightly elevated. This region corresponds to the quasi-geometric-optics regime near the lens caustic, where the amplification factor exhibits a high dynamic range due to the constructive interference of two bright images. The sharp central peak (∣F∣ ≫ 1) introduces high-frequency gradients in the amplitude envelope, posing a challenge for the network to resolve simultaneously with the phase oscillations. Conversely, in the large impact parameter regime (y ≳ 3), the system transitions into the weak-lensing diffraction tail. Here the challenge shifts from amplitude dynamic range to phase frequency: The oscillation frequency scales as ωy, leading to dense interference fringes. The domain decomposition strategy reduces the errors in this high-frequency limit. The SIREN model maintains a uniform accuracy profile across these physically distinct regimes, demonstrating robustness against both the amplitude singularities of strong lensing and the oscillatory catastrophes of weak lensing.

Zoom In Zoom Out Reset image size

To provide a comprehensive view of the error behavior across the entire parameter space, particularly near the complex diffraction-interference transition, we visualize the relative-error distribution in Figure 6. The parameter space is represented with a logarithmic axis for the dimensionless frequency ω and a linear axis for the impact parameter y. As illustrated, the model demonstrates remarkable robustness, with the majority of the parameter space generally maintaining relative errors at the level. Importantly, in the critical transition regime (moderate y ∼ 1 and intermediate ω ∼ 1–10) where the amplification factor exhibits dense and rapid interference fringes, the SIREN architecture maintains high stability without introducing numerical artifacts. We observe a slight elevation in relative errors, occasionally approaching , localized strictly in the very-low-frequency regime (ω ∈ [10⁻², 10⁻¹]). This minor localized deviation is primarily attributed to the relatively sparse sampling density within this specific subinterval during our training phase. Refining the sampling strategy in this low-frequency regime will be a key direction for our future optimization work.

Zoom In Zoom Out Reset image size

Figure 6 illustrates the statistical distribution of these errors. The histogram for off-grid points (Figure 7(a)) is heavily concentrated between 10⁻⁸ and 10⁻³. The cumulative distribution function (CDF) confirms that more than 99% of test points achieve subpercent accuracy. The training-grid statistics (Figure 7(b)) closely mirror the off-grid results, indicating that the model performance on off-grid samples is consistent with the training error.

Zoom In Zoom Out Reset image size

To verify the physical fidelity of the predictions, Figure 8 presents the waveforms at two representative impact parameters: y = 0.3 (high dynamic range) and y = 1.0 (intermediate regime). At y = 0.3 (Figure 8(a)), the model accurately captures the pronounced central peak (∣F∣ ≈ 2.5) and the subsequent Airy pattern with relative errors consistently below 10⁻³. At y = 1.0 (Figure 8(b)), the diffraction pattern transitions to a regime with a smoother envelope peak (∣F∣ ≈ 1.4) but increasingly rapid oscillations. The model tracks both the envelope and the phase of the oscillations robustly, maintaining per-mille-level accuracy throughout the spectrum.

Zoom In Zoom Out Reset image size

We assess the spectral fidelity of the model using representative impact parameters y = 0.3 and y = 1.0. For a fixed y, we characterize the spectral content of the complex amplification factor by its power spectral density (PSD). We estimate the PSD using the periodogram, defined as the squared modulus of the discrete Fourier transform (DFT; D. B. Percival & A. T. Walden 1993),

where denotes the DFT applied to the discretized F(ω, y) on a uniform ω-grid, and ν represents the spectral frequency conjugate to ω (quantifying the rate of oscillation of the diffraction fringes). Since we focus on the relative spectral structure, the standard normalization factor is omitted. In practice, the series is mean-subtracted, multiplied by a Tukey window (α = 0.25) to mitigate spectral leakage, and zero-padded to the next power of two prior to the DFT. We verified that the qualitative conclusions are robust to moderate variations in the windowing and padding parameters.

Figure 9 compares the PSDs of the analytic PML solution and the model prediction. At y = 0.3 (Figure 9(a)), the principal diffraction peak and the power-law tail at high ν overlap over several orders of magnitude, indicating that the model preserves both the broadband envelope and the fine-scale oscillatory content. At y = 1.0 (Figure 9(b)), the peak alignment and the asymptotic slope are accurately recovered. While minor residuals appear at intermediate frequencies due to the dense fringe pattern, they do not deviate significantly from the ground-truth scaling. Overall, these tests confirm that the SIREN model successfully reproduces the multiscale spectral structure of F(ω, y) across the relevant domain.

Zoom In Zoom Out Reset image size

To assess the practical utility of our approach, we compare the SIREN model against the classical ZPAE method (Section 3.1.5), as well as the AE and zero-points integral (ZPI) methods, in terms of both accuracy and computational efficiency. Additionally, we include the high-precision analytical evaluation computed via the mpmath library to serve as our absolute ground-truth baseline.

5.3.1. Accuracy Robustness

We evaluate the accuracy using standard empirical parameters k = 5, ν = 7 (and m = 15 for ZPI) at a fixed frequency ω = 10.0. Figure 10 compares the relative error ∣ΔF∣/∣F∣ as the impact parameter y varies from 0.2 to 5.0.

Zoom In Zoom Out Reset image size

The ZPAE and AE methods exhibit extreme sensitivity to the alignment of their expansion centers. While they achieve vanishing residuals (below 10⁻⁸) specifically at the anchor point y = 0.2, the approximation error diverges rapidly as y deviates from this configuration. The error increases by orders of magnitude, so the methods require computationally expensive retuning to remain reliable. Meanwhile, the ZPI method, although more stable across the domain, yields a consistently high baseline error under standard truncation settings. In comparison, the SIREN model maintains a uniformly low relative error of order 10⁻³ across the entire domain. This confirms that the neural network successfully captures the global oscillatory behavior of the diffraction integral without requiring case-by-case hyperparameter adjustments.

5.3.2. Computational Efficiency

To quantify inference latency, we benchmarked the SIREN model against the Mpmath, AE, ZPI, and ZPAE implementations by sampling across representative domains of both ω and y.

As illustrated in Figure 11, the SIREN model achieves a consistent latency of ∼10⁻⁵ s per sample. This represents a substantial speedup, being approximately 2 orders of magnitude faster than the high-precision analytic Mpmath baseline (∼10⁻³ s), up to 4 orders faster than the unoptimized ZPAE and AE implementations (∼10⁻¹ s), and 5 orders faster than the segmented ZPI method (∼10⁰ s).

Zoom In Zoom Out Reset image size

Furthermore, the results in Figures 11(a) and (b) demonstrate the robustness of the neural estimator across the parameter space. The computational cost remains stable at this low latency level regardless of variations in the impact parameter y or the frequency ω. This consistent sub-millisecond efficiency is essential for supporting real-time lensing computations and large-scale parameter scans in GW data analysis.

We decompose the amplification factor as

and track its phase ϕ. The phase structure is known to affect detectability in matched filtering and to introduce parameter degeneracies. Following the standard definition in circular statistics, we adopt the shortest arc distance as the phase error measure:

Here we additionally report the phase error to highlight the phase fidelity in the wave optics regime.

Table 4 summarizes the mean, median, 95th percentile, and maximum phase errors (in milliradians) for representative impact parameters in each of the four y-intervals. Notably, the maximum phase error observed in our testing is approximately 15.36 mrad. To understand this value, we compare it against the instrumental calibration uncertainty of current GW detectors. During the LIGO O3 observing run, the systematic phase calibration error was estimated to be approximately 2°–5° (∼35–87 mrad) in the sensitive frequency band. Our model’s residual phase error remains well below this instrumental floor. Furthermore, according to waveform accuracy standards, a phase error of this magnitude corresponds to a negligible mismatch, ensuring that the neural network approximation does not introduce significant systematic bias into parameter estimation or reduce the recoverable signal-to-noise ratio (SNR).

Table 4. Phase Error Statistics for Representative y Values in Each Domain Interval (Unit: Mrad)

Interval	y	Mean	Median	95th %ile	Maximum
	[0.2, 1.0]	1.43	1.10	3.93	15.36
	[1.0, 3.0]	1.33	1.05	3.46	9.15
	[3.0, 6.0]	0.26	0.22	0.61	7.56
	[6.0, 10.0]	0.99	0.82	2.46	5.60

Download table as:  ASCIITypeset image

To validate the generalization capability of our framework beyond the PML, we extended the pipeline to the SIS model. The primary objective was to verify that the SIREN architecture captures intrinsic wave-optics diffraction behavior rather than overfitting to a specific potential. We employed the identical network architecture and training strategy used for the PML baseline, modifying only the ground-truth generation module to incorporate the SIS diffraction integral (G. Pagano et al. 2020). As shown in Figure 12, the model accurately reconstructs the spectral amplification for representative impact parameters (y = 0.3 and y = 1.0), successfully capturing both the high-amplitude diffraction peaks and the dense interference fringes while maintaining relative errors generally below the 10⁻³ level. To quantitatively verify this generalization, we evaluated the model on random test points strictly independent of the training grid. The error statistics and distributions are detailed in Appendix B. The SIS model achieves an overall mean relative error of 2.29 × 10⁻⁴, maintaining stable accuracy across the entire evaluated parameter space. This robust performance on out-of-sample data further demonstrates the reliability of the architecture. This successful transferability implies that the proposed neural representation relies on the underlying wave physics rather than specific geometric symmetries, suggesting its broader applicability to more complex lens models.

Zoom In Zoom Out Reset image size

To validate the generalization capability of our framework beyond the PML, we extended the pipeline to the SIS model. The primary objective was to verify that the SIREN architecture captures intrinsic wave-optics diffraction behavior rather than overfitting to a specific potential. We employed the identical network architecture and training strategy used for the PML baseline, modifying only the ground-truth generation module to incorporate the SIS diffraction integral (G. Pagano et al. 2020). As shown in Figure 12, the model accurately reconstructs the spectral amplification for representative impact parameters (y = 0.3 and y = 1.0). To quantitatively verify this generalization, we evaluated the model on random test points strictly independent of the training grid. The error statistics and distributions are detailed in Appendix B. The SIS model achieves an overall mean relative error of 2.29 × 10⁻⁴, maintaining stable accuracy across the entire evaluated parameter space. This robust performance on out-of-sample data further demonstrates the reliability of the architecture. This successful transferability implies that the proposed neural representation relies on the underlying wave physics rather than specific geometric symmetries, suggesting its broader applicability to more complex lens models.