Walk-on-Cubes Monte Carlo Simulation for nonisotropic fractional Laplace, Helmholtz, and Yukawa equations
July 12, 2026
Abstract
We study the nonisotropic fractional analogs of Laplace, Helmholtz and Yukawa equations. We provide a Duffin correspondence for the Yukawa equation and a Feynman–Kac reconstruction for the Helmholtz equation. The foundation of our analysis is the fact that the nonisotropic fractional Laplace equation is related to a symmetric \(\alpha\)-stable Lévy process with independent identically distributed components. By using this relation we provide a Walk-on-Cubes algorithm that simulates the solutions of Helmholtz and Yukawa equations.
1 Introduction↩︎
Let \(\mathbf{x}=(x_1,\ldots,x_d)\in\mathbb{R}^d\). Let \(\alpha\in(0,2)\) and let \(A_\mathbf{x}^\alpha\) be the nonisotropic fractional Laplacian: \[A^\alpha_\mathbf{x}= \sum_{i=1}^d -(-\Delta_{x_i})^{\alpha/2},\] where \(-(-\Delta_{x_i})^{\alpha/2}\) is the one-dimensional fractional Laplacian defined by the singular integral \[\label{eq:flaplace} (-\Delta_{x_i})^{\alpha/2}f(x_i) = C_\alpha\int_{-\infty}^\infty \frac{f(x_i)-f(y)}{|x_i-y|^{1+\alpha}}\, \mathrm{d}y.\tag{1}\] Here the positive one-dimensional normalization is \[\label{eq:calpha-definition} C_\alpha :=\frac{2^\alpha \Gamma((1+\alpha)/2)}{\pi^{1/2}|\Gamma(-\alpha/2)|} =\frac{\alpha 2^{\alpha-1}\Gamma((1+\alpha)/2)}{\pi^{1/2}\Gamma(1-\alpha/2)}.\tag{2}\]
Let \(D\subset \mathbb{R}^d\) be a bounded domain, and let \(\lambda\in\mathbb{R}\). We consider the Dirichlet problem of finding for a given \(g\colon D^c\to\mathbb{R}\) a solution \(u\colon D\to \mathbb{R}\) such that \[\begin{align} \tag{3} A^\alpha_\mathbf{x}u(\mathbf{x}) &=& \lambda u(\mathbf{x}) \quad if\mathbf{x}\in D, \\ \tag{4} u(\mathbf{x}) &=& g(\mathbf{x}) \quad if\mathbf{x}\in D^c. \end{align}\] If \(\lambda<0\) 3 –4 is called the fractional Helmholtz equation, if \(\lambda>0\) it is called the fractional Yukawa equation, and for \(\lambda=0\) it is called the fractional Laplace equation.
Recently Kyprianou et al. [1] studied the isotropic fractional Laplace equation and provided the Walk-on-Spheres algorithm for its simulation. Our work is related to theirs, with two major differences: first, we consider the nonisotropic case corresponding to independent identically distributed component Lévy process as opposed to the isotropic Lévy process. Second, we also consider fractional Helmholtz and Yukawa equations.
To our knowledge, the nonisotropic operator \(A^\alpha_\mathbf{x}\) is relatively unstudied in the literature. One related study is Dybiec and Szczepaniec [2].
The rest of the paper is organized as follows: In Section 2 we provide the Duffin correspondence relating the fractional Laplace equation and the fractional Yukawa equation. In Section 3 we provide the theoretical background for our Monte Carlo simulation, including the obstruction to the same spatial Duffin lifting in the Helmholtz regime. Finally, in Section 4 we introduce our Walk-on-Cubes algorithm and provide a simulation.
2 Duffin Correspondence↩︎
Duffin [3] provided a correspondence that transforms the classical Dirichlet problem of the Yukawa equation into a Dirichlet problem of the Laplace equation. This was later studied and extended to the Helmholtz equation in [4]–[6].
Since the operator \(A^\alpha_\mathbf{x}\) has a “coordinate-sum” form, we can provide the Duffin correspondence for it in a similar way as for the classical Laplacian. The key ingredient for the Duffin correspondence is to find the non-zero solution \(f=f_\lambda\) to the one-dimensional eigenvalue problem \[\label{eq:eigenfunction} (-\Delta)^{\alpha/2}_yf(y) = \lambda f(y) \quadfor\quad y\in\mathbb{R}.\tag{5}\] For the Yukawa regime (\(\lambda > 0\)), 5 has the explicit pointwise solution \(f(y)=\cos(\lambda^{1/\alpha}y)\). For the Helmholtz regime (\(\lambda < 0\)), the heavy-tailed polynomial decay of the fractional operator rules out bounded analytical solutions, so the standard Duffin lifting is not available, as will be shown in Section 3.3.
Analytical Derivation of the One-Dimensional Eigenfunction↩︎
A standard reference for the Fourier transformation of the fractional Laplacian is Proposition 3.3 of Di Nezza, Palatucci, and Valdinoci [7].
Proposition 1 (Fourier symbol of the one-dimensional fractional Laplacian). Let \(0<\alpha<2\) and let \(\varphi\) be sufficiently regular, for instance \(\varphi\in\mathcal{S}(\mathbb{R})\). With the Fourier transform convention \[\mathcal{F}\varphi(\xi)=\int_{\mathbb{R}}\mathrm{e}^{-\mathrm{i}y\xi}\varphi(y)\,\mathrm{d}y,\] Proposition 3.3 of [7], applied with \(s=\alpha/2\), gives \[\mathcal{F}\!\left[(-\Delta)^{\alpha/2}\varphi\right](\xi) =|\xi|^\alpha \mathcal{F}\varphi(\xi).\] Equivalently, \[(-\Delta)^{\alpha/2}\varphi =\mathcal{F}^{-1}\!\left(|\xi|^\alpha\mathcal{F}\varphi(\xi)\right).\]
Applying this proposition to the Fourier basis \(u(y)=e^{\mathrm{i}\omega y}\), whose Fourier transform is concentrated at the frequency \(\xi=\omega\), gives \[(-\Delta)^{\alpha/2}_y e^{\mathrm{i}\omega y}=|\omega|^\alpha e^{\mathrm{i}\omega y}.\] For completeness, we give the direct derivation below, although the general idea is the same as in [7]. Thus the one-dimensional eigenvalue equation 5 is solved by plane waves with frequencies satisfying \(|\omega|^\alpha=\lambda\). In this direct computation we use the normalization constant from 2 ; equivalently, \[C_\alpha =\left(\int_{\mathbb{R}}\frac{1-\cos z}{|z|^{1+\alpha}}\,\mathrm{d}z\right)^{-1} =\frac{\Gamma(1+\alpha)\sin(\pi\alpha/2)}{\pi}.\] This is the standard normalization that makes the singular-integral definition have Fourier symbol \(|\xi|^\alpha\); see [8] and [7].
Let \(f(y)=e^{i\omega y}\), where \(\omega\in\mathbb{R}\) is fixed and \(y\in\mathbb{R}\) is the spatial variable. Substituting this into 5 gives \[(-\Delta)^{\alpha/2}_y e^{i\omega y} = C_\alpha\int_{-\infty}^{\infty} \frac{e^{i\omega y}-e^{i\omega x}}{|y-x|^{1+\alpha}}\,\mathrm{d}x.\] With the change of variables \(v=x-y\), we obtain \[(-\Delta)^{\alpha/2}_y e^{i\omega y} = C_\alpha e^{i\omega y}\int_{-\infty}^{\infty} \frac{1-e^{i\omega v}}{|v|^{1+\alpha}}\,\mathrm{d}v.\] Using \(e^{i\omega v}=\cos(\omega v)+i\sin(\omega v)\), we may split the integral into real and imaginary parts. And the denominator is even and \(\sin(\omega v)\) is odd so the principal value integral of a function is zero. The imaginary part by symmetry. Hence \[(-\Delta)^{\alpha/2}_y e^{i\omega y} = C_\alpha e^{i\omega y}\int_{-\infty}^{\infty} \frac{1-\cos(\omega v)}{|v|^{1+\alpha}}\,\mathrm{d}v.\] Near \(v=0\), the numerator behaves like \(\omega^2 v^2/2\), so the integrand behaves like \(|v|^{1-\alpha}\) and is locally integrable for \(\alpha\in(0,2)\). Thus the principal value integral agrees with the improper integral. For \(\omega\ne0\), set \(u=|\omega|v\); the case \(\omega=0\) is immediate. Then \[\int_{-\infty}^{\infty}\frac{1-\cos(\omega v)}{|v|^{1+\alpha}}\,\mathrm{d}v =|\omega|^\alpha \int_{-\infty}^{\infty}\frac{1-\cos u}{|u|^{1+\alpha}}\,\mathrm{d}u.\] The cancellation with \(C_\alpha\) is explicit: \[\begin{align} C_\alpha \int_{-\infty}^{\infty}\frac{1-\cos(\omega v)}{|v|^{1+\alpha}}\,\mathrm{d}v &= |\omega|^\alpha \frac{\Gamma(1+\alpha)\sin(\pi\alpha/2)}{\pi} \int_{-\infty}^{\infty}\frac{1-\cos u}{|u|^{1+\alpha}}\,\mathrm{d}u\\ &= |\omega|^\alpha \frac{\Gamma(1+\alpha)\sin(\pi\alpha/2)}{\pi} \cdot \frac{\pi}{\Gamma(1+\alpha)\sin(\pi\alpha/2)}\\ &=|\omega|^\alpha. \end{align}\] Therefore, we obtain \[(-\Delta)^{\alpha/2}_y e^{i\omega y}=|\omega|^\alpha e^{i\omega y}.\] Consequently, the solutions of 5 with \(\lambda>0\) are generated by frequencies satisfying \(|\omega|^\alpha=\lambda\), that is, \[\omega=\pm \lambda^{1/\alpha}.\] The real eigenspace is therefore spanned by \[f_\lambda(y)=A\cos\!\left(\lambda^{1/\alpha}y\right)+B\sin\!\left(\lambda^{1/\alpha}y\right). \label{eq:report3-real-eigenspace}\tag{6}\] The following is the core of the Duffin correspondence.
Theorem 1 (Duffin Correspondence). Let \(D\subset \mathbb{R}^d\) be a bounded domain. Let \(\lambda>0\). Let \(f\colon\mathbb{R}\to\mathbb{R}\) be a non-zero eigenfunction given by 5 . Define \[\mathcal{A}_{\mathbf{x},y}^{\alpha}:=A_{\mathbf{x}}^{\alpha}-(-\Delta_y)^{\alpha/2}.\] Set \(v(\mathbf{x},y) = u(\mathbf{x})f(y)\). Then \(v\) solves \[\label{eq:duffin95correspondence-1} \mathcal{A}_{\mathbf{x},y}^{\alpha}v(\mathbf{x},y) = 0 \quadon\quad D\times\mathbb{R}\qquad{(1)}\] if and only if \(u\) solves \[\label{eq:duffin95correspondence-2} A^{\alpha}_\mathbf{x}u(\mathbf{x}) = \lambda u(\mathbf{x}) \quadon\quad D.\qquad{(2)}\]
Proof. Suppose \(u\) satisfies ?? . Then, since \(f\) satisfies 5 , we have \[\begin{align} \mathcal{A}_{\mathbf{x},y}^{\alpha}v(\mathbf{x},y) &=& \sum_{i=1}^{d} -(-\Delta_{x_i})^{\alpha/2}v(\mathbf{x},y) - (-\Delta_y)^{\alpha/2}v(\mathbf{x},y) \\ &=& f(y)\sum_{i=1}^{d} -(-\Delta_{x_i})^{\alpha/2}u(\mathbf{x}) - u(\mathbf{x})(-\Delta_y)^{\alpha/2}f(y) \\ &=& f(y)A^{\alpha}_\mathbf{x}u(\mathbf{x}) -u(\mathbf{x})\lambda f(y) \\ &=& f(y)\left[\lambda u(\mathbf{x}) - \lambda u(\mathbf{x})\right] \\ &=& 0. \end{align}\] So, \(v\) satisfies ?? .
Suppose now that \(v\) satisfies ?? . Then we have by the calculations above that \[\begin{align} 0 &=& \mathcal{A}_{\mathbf{x},y}^{\alpha}v(\mathbf{x},y) \\ &=& f(y)\left[A^{\alpha}_\mathbf{x}u(\mathbf{x}) -\lambda u(\mathbf{x})\right]. \end{align}\] Since this identity holds for all \(y\in\mathbb{R}\) and \(f\not\equiv0\), there exists \(y_0\) such that \(f(y_0)\neq0\); hence \(u\) must satisfy ?? . ◻
In particular, the Duffin correspondence provides us a way to solve the Yukawa fractional Dirichlet problem 3 –4 as the Laplace problem \[\label{eq:fractional-laplace} \mathcal{A}_{\mathbf{x},y}^{\alpha}v(\mathbf{x},y) = 0 \quadon\quad D\times\mathbb{R},and \tag{7}\] \[\label{eq:fractional-laplace-boundary} v(\mathbf{x},y) = g(\mathbf{x})f(y) \quadon the exterior of\quad D\times\mathbb{R}.\tag{8}\]
Recall the Kakutani connection to the fractional Laplace equation [9]. Consider the fractional Laplace Dirichlet problem \[A^\alpha_\mathbf{x}u(\mathbf{x}) = 0onD\] \[u(\mathbf{x}) = g(\mathbf{x})onD^c\] Let \(\mathbf{X}=(X^1,\ldots,X^d)\) be an i.i.d. component symmetric \(\alpha\)-stable Lévy process, i.e., each component \(X^j\), \(j=1,\ldots, d\) have the same law given by the characteristic function \[\mathbb{E}\left[\mathrm{e}^{\mathrm{i}\theta X^1_t}\right] = \mathrm{e}^{-|\theta|^\alpha t}.\] Then we have a stochastic representation of the solution \(u\): \[\label{eq:fractional-kakutani} u(\mathbf{x}) = \mathbb{E}^\mathbf{x}\left[g(\mathbf{X}_\tau)\right]\tag{9}\] where \(\tau =\inf\{t>0; \mathbf{X}_t\in D^c\}\). The connection 9 provides a method for solving the fractional Laplace Dirichlet problem stochastically. One simply simulates a large number of paths \(\mathbf{X}(i)\), \(i=1,\ldots, N\), of the \(\alpha\)-stable Lévy process \(\mathbf{X}\) starting from \(\mathbf{x}\) and the approximate solution to \(u\) is \[\label{eq:mc0} \hat{u}(\mathbf{x}) = \frac{1}{N}\sum_{i=1}^N g(\mathbf{X}_{\tau(i)}(i)),\tag{10}\] where \(\tau(i)\) is the first time the sample \(\mathbf{X}(i)\) enters the domain \(D^c\).
Remark 1. To use the formula 10 one needs to simulate whole trajectories with fine time mesh of \(\mathbf{X}\), and this is computationally heavy and leads to accumulation of errors. To overcome the problem we provide the Walk-on-Cubes algorithm for the simulation. This algorithm is based on the well-known Walk-on-Spheres algorithm for the classical Laplace and Brownian motion case, see Muller [10] generalizations of which where studied in [11].
3 Theoretical Background↩︎
3.1 Operator Definition and Unified Framework in the Pure Fractional Laplace Regime↩︎
We consider the exterior Dirichlet problem \[\label{eq:unified-framework-pde} A_{\mathbf{x}}^{\alpha}u(\mathbf{x})-\lambda u(\mathbf{x})=0,\qquad \mathbf{x}\in D\subset\mathbb{R}^d,\tag{11}\] with boundary condition \[u(\mathbf{x})=g(\mathbf{x}),\qquad \mathbf{x}\in D^c,\] where \[A_{\mathbf{x}}^{\alpha}=\sum_{i=1}^d-(-\Delta_{x_i})^{\alpha/2},\qquad 0<\alpha<2.\] The operator is Cartesian separable: \[A_{\mathbf{x}}^{\alpha}=\sum_{i=1}^d A_{x_i}^{\alpha},\] hence each coordinate is governed by a one-dimensional fractional generator. This additive decomposition reflects the nonisotropic structure and leads to a coordinatewise jump mechanism. The same decomposition also explains why the Walk-on-Cubes (WoC) method is natural for this model. At each state \(\mathbf{z}\in D\), WoC selects the maximal axis-aligned cube \(\mathbf{z}+rQ\) with \(Q=[-1,1]^d\) and \[r=\operatorname{dist}_{\infty}(\mathbf{z},\partial D).\]
Because the \(d\) coordinates evolve independently and the local geometry is described by the coordinate-wise bounds \(|z_i-z_{0,i}|<r\) for \(i=1,\ldots,d\), the \(L_{\infty}\) metric is the natural metric here. Together with the self-similarity of symmetric \(\alpha\)-stable motions, this gives an exact rescaling rule from the unit-cube exit law to any local cube, which underlies the WoC construction. When \(\lambda=0\), 11 reduces to the pure fractional Laplace equation \[A_{\mathbf{x}}^{\alpha}u(\mathbf{x})=0\quad\text{in }D.\]
From the probabilistic viewpoint, the associated process \(\mathbf{X}_t=(X_t^1,\ldots,X_t^d)\) is a drift-free symmetric pure-jump Lévy motion with independent components, which describes unbiased anomalous diffusion. In particular, there is no deterministic drift, killing, or damping term. If \[\tau=\inf\{t>0:\,\mathbf{X}_t\notin D\},\] then the solution admits the exit distribution representation \[u(\mathbf{x})=\mathbb{E}_{\mathbf{x}}\!\left[g(\mathbf{X}_{\tau})\right].\] This representation underlies the Monte Carlo estimator used in the WoC algorithm.
3.2 Duffin Extension and Oscillatory Screening for the Fractional Yukawa Regime↩︎
We now turn to the Yukawa case \(\lambda>0\) in 11 . Set \[\kappa:=\lambda^{1/\alpha},\qquad f_\lambda(w):=\cos(\kappa w).\] Since the one-dimensional pointwise eigen relation \[(-\Delta_w)^{\alpha/2}f_\lambda(w)=\lambda f_\lambda(w),\qquad w\in\mathbb{R},\] holds, we define the lifted field \[U(\mathbf{x},w):=u(\mathbf{x})f_\lambda(w).\] Introduce the potential free lifted operator \[\mathcal{A}_{\mathbf{x},w}^{\alpha}:=A_{\mathbf{x}}^{\alpha}-(-\Delta_w)^{\alpha/2}.\] Then the Yukawa equation in \(D\subset\mathbb{R}^d\) is equivalent to the pure fractional Laplace equation in the augmented space: \[\mathcal{A}_{\mathbf{x},w}^{\alpha}U(\mathbf{x},w)=0,\qquad (\mathbf{x},w)\in D\times\mathbb{R},\] with exterior data \[U(\mathbf{x},w)=g(\mathbf{x})\cos(\lambda^{1/\alpha}w),\qquad (\mathbf{x},w)\in D^c\times\mathbb{R}.\] Let \(\widetilde{\mathbf{X}}_t=(\mathbf{X}_t,W_t)\), where \(\mathbf{X}_t\) is the \(d\)-dimensional nonisotropic symmetric \(\alpha\)-stable motion generated by \(A_{\mathbf{x}}^{\alpha}\) and \(W_t\) is an independent one-dimensional symmetric \(\alpha\)-stable process. With \[\tau_D:=\inf\{t>0:\mathbf{X}_t\notin D\},\] the Duffin representation becomes \[u(\mathbf{x})=\mathbb{E}_{(\mathbf{x},0)}\!\left[g(\mathbf{X}_{\tau_D})\cos\!\big(\lambda^{1/\alpha}W_{\tau_D}\big)\right].\] Hence the Monte Carlo estimator uses i.i.d. samples of the oscillatory payoff \[\xi_k:=g(\mathbf{X}^{(k)}_{\tau_D})\cos\!\big(\lambda^{1/\alpha}W^{(k)}_{\tau_D}\big),\] and, under \(\mathbb{E}|\xi_1|<\infty\), the strong law yields \[\frac{1}{N}\sum_{k=1}^N\xi_k\xrightarrow[N\to\infty]{\mathrm{a.s.}}u(\mathbf{x}).\] The screening mechanism is oscillatory rather than dissipative. For fixed \(t>0\), \[\mathbb{E}\!\left[\cos\!\big(\lambda^{1/\alpha}W_t\big)\right] =\Re\,\mathbb{E}\!\left[\exp\!\big(i\lambda^{1/\alpha}W_t\big)\right] =e^{-\lambda t},\] so phase mixing leads to cancellation at the expectation level. From a pathwise viewpoint, trajectories started deep inside \(D\) typically require many WoC updates before they exit, and the heavy-tailed jump structure gives a broad right tail for \(\tau_D\); accordingly, the scale of \(W_{\tau_D}\) (of order \(\tau_D^{1/\alpha}\)) becomes large. The factor \(\cos(\lambda^{1/\alpha}W_{\tau_D})\in[-1,1]\) then oscillates strongly across samples, which leads to substantial cancellation in empirical averages. At the level of the computed solution, this appears as an interior region that is nearly flat and close to zero, together with a sharp boundary layer where \(\tau_D\) is small and the boundary data are less screened.
3.3 Failure of the Duffin Extension in the Helmholtz Regime and a Feynman–Kac Reconstruction↩︎
The Helmholtz regime corresponds to \(\lambda<0\) in 11 . Writing \[\kappa:=-\lambda>0,\] the equation becomes \[A_{\mathbf{x}}^{\alpha}u(\mathbf{x})+\kappa u(\mathbf{x})=0\quad\text{in }D.\] A direct reuse of the Duffin space extension paradigm leads to a genuine functional analytic obstruction. Formally, one would seek an auxiliary profile \(f\) solving \[(-\Delta_w)^{\alpha/2}f(w)=-\kappa f(w),\] whose natural candidates are exponentially growing/decaying modes. However, for the fractional operator with heavy-tailed kernel, \[\label{eq:helmholtz-kernel-pv} (-\Delta_w)^{\alpha/2}f(w)=c_\alpha\int_{\mathbb{R}}\frac{f(w)-f(w+z)}{|z|^{1+\alpha}}\,\mathrm{d}z,\tag{12}\] the tail \(|z|^{-(1+\alpha)}\) is only polynomially decaying. If \(f(w)=e^{\beta w}\) with \(\beta\neq0\), then \[\frac{f(w)-f(w+z)}{|z|^{1+\alpha}} =e^{\beta w}\frac{1-e^{\beta z}}{|z|^{1+\alpha}},\] then for \(\beta>0\) the tail \(z\to+\infty\) grows exponentially, while for \(\beta<0\) the tail \(z\to-\infty\) grows exponentially. Hence the integral in 12 diverges at infinity. Therefore the naive lifted space Duffin correspondence is not well-posed for \(\lambda<0\), and a simple spatial augmentation is not available.
The appropriate replacement is a time enhanced Feynman–Kac reconstruction. Let \(\mathbf{X}_t=(X_t^1,\ldots,X_t^d)\) be the nonisotropic symmetric \(\alpha\)-stable process generated by \(A_{\mathbf{x}}^{\alpha}\), and let \[\tau_D:=\inf\{t>0:\mathbf{X}_t\notin D\}.\] Then the Helmholtz solution is represented by \[u(\mathbf{x})=\mathbb{E}_{\mathbf{x}}\!\left[e^{-\lambda\tau_D}g(\mathbf{X}_{\tau_D})\right] =\mathbb{E}_{\mathbf{x}}\!\left[e^{\kappa\tau_D}g(\mathbf{X}_{\tau_D})\right].\] Thus the Helmholtz weight grows exponentially rather than decays. This is the source of the instability in Monte Carlo simulation. In WoC form, each local step from a cube of radius \(r\) inherits the space–time self-similarity of the \(\alpha\)-stable dynamics: \[\Delta\tau=r^{\alpha}\tau_{\mathrm{std}},\] where \(\tau_{\mathrm{std}}\) is the first exit time from the unit cube. Along one trajectory, \[\tau_{\mathrm{total}}=\sum_{j=0}^{N-1} r_j^{\alpha}\tau_{\mathrm{std}}^{(j)}.\] Here \(\tau_{\mathrm{total}}\) is a discrete approximation of the continuous stopping time \(\tau_D\), obtained from the space–time self-similarity of the \(\alpha\)-stable process. Under appropriate regularity conditions on the boundary and the fractional Poisson kernel, and for a consistent WoC refinement (\(\epsilon\to0\), \(N\to\infty\)), one expects \[\tau_{\mathrm{total}}^{(\epsilon,N)}\xrightarrow[\epsilon\to0,\,N\to\infty]{\mathcal{D}}\tau_D.\] Hence the path payoff is reweighted as \[\Xi=g(\mathbf{X}_{\tau_D})\exp\!\big(-\lambda\tau_{\mathrm{total}}\big),\] and the Monte Carlo estimator is the empirical mean of i.i.d. copies of \(\Xi\). Let \(\lambda_1(D)>0\) denote the principal Dirichlet eigenvalue of \(-A_{\mathbf{x}}^{\alpha}\) on \(D\). The single particle Feynman–Kac representation is meaningful only in the subcritical gauge regime [12] \[-\lambda<\lambda_1(D),\] which ensures finiteness of the relevant exponential moments. When \(-\lambda\ge\lambda_1(D)\), the variance, or even the expectation, can blow up, and one typically needs branching random walk mechanisms rather than a pure single trajectory reweighting scheme.
3.4 Existence of the Principal Eigenvalue and Its WoC Simulation↩︎
We now show that the principal Dirichlet eigenvalue \(\lambda_1(D)\) exists and explain why the Walk-on-Cubes exit time simulation can be used to estimate it. This is the quantity that controls the admissible range of the Helmholtz parameter in the Feynman–Kac reconstruction.
Definition 1 (Rectilinear stable process and killed semigroup). Let \(\mathbf{X}=(X^1,\ldots,X^d)\) be the rectilinear symmetric \(\alpha\)-stable Lévy process, i.e., the coordinates are independent one-dimensional symmetric \(\alpha\)-stable Lévy processes. Equivalently, by the Lévy–Khintchine representation [13], \[\mathbb{E}_0\!\left[e^{\mathrm{i}\langle \xi,X_t\rangle}\right] =\exp\!\left(-t\Psi(\xi)\right), \qquad \Psi(\xi)=\sum_{i=1}^d|\xi_i|^\alpha .\] The generator is \(A^\alpha=-L^\alpha\), where \[L^\alpha:=\sum_{i=1}^d(-\Delta_{x_i})^{\alpha/2}.\] For a bounded open set \(D\subset\mathbb{R}^d\), define \[\tau_D:=\inf\{t>0:X_t\notin D\},\] and the killed semigroup on \(D\) by \[T_t^D f(x):=\mathbb{E}_x\!\left[f(X_t)\mathbf{1}_{\{\tau_D>t\}}\right].\] The terminology “rectilinear stable process” and the Dirichlet transition density theory used below follow Chen, Hu, and Zhao [14].
The corresponding Dirichlet form is the coordinate axis form \[\label{eq:rectilinear-dirichlet-form} \mathcal{E}_D(f,f) =\frac{C_\alpha}{2}\sum_{i=1}^d \int_{\mathbb{R}^d}\int_{\mathbb{R}} \frac{\left(f(x+he_i)-f(x)\right)^2}{|h|^{1+\alpha}}\,\mathrm{d}h\,\mathrm{d}x, \qquad f=0\;on D^c,\tag{13}\] where \(e_i\) is the \(i\)th coordinate vector. Formula 13 is the quadratic form version of the singular integral generator 1 ; see also the equivalent definitions of the fractional Laplacian in [8].
Definition 2 (Principal Dirichlet eigenvalue). The principal Dirichlet eigenvalue of \(L^\alpha=-A^\alpha\) on \(D\) is \[\lambda_1(D):= \inf\left\{\mathcal{E}_D(f,f): f\in\mathcal{F}_D,\;\|f\|_{L^2(D)}=1\right\},\] where \(\mathcal{F}_D\) is the closure of smooth functions compactly supported in \(D\) under the energy norm associated with 13 . This is the Rayleigh–Ritz characterization of the bottom of the Dirichlet spectrum; compare the spectral treatment of fractional Dirichlet eigenvalues in Kwaśnicki [15].
Proposition 2 (Existence and survival tail characterization of \(\lambda_1(D)\)). Assume that \(D\) is bounded and belongs to the geometric class for which the killed rectilinear stable process is irreducible and has a strictly positive Dirichlet transition density in the sense of [14]. Then \(T_t^D\) is a compact, self-adjoint, positivity improving operator on \(L^2(D)\) for every \(t>0\). Consequently, there are eigenvalues \[\label{eq:discrete-dirichlet-spectrum} 0<\lambda_1(D)<\lambda_2(D)\le\lambda_3(D)\le\ldots,\qquad \lambda_n(D)\to\infty,\qquad{(3)}\] and an orthonormal basis \(\{\phi_n\}_{n\ge1}\) of \(L^2(D)\) such that \[\label{eq:killed-semigroup-eigen-exp} T_t^D\phi_n=e^{-t\lambda_n(D)}\phi_n.\qquad{(4)}\] Moreover, \(\lambda_1(D)\) is simple and \(\phi_1\) can be chosen strictly positive on \(D\). For every point \(x\) for which the transition density expansion is valid, \[\label{eq:survival-spectral-expansion} S_x(t):=\mathbb{P}_x(\tau_D>t) =\sum_{n=1}^{\infty}e^{-t\lambda_n(D)}\phi_n(x)\int_D\phi_n(y)\,\mathrm{d}y,\qquad{(5)}\] and therefore \[\label{eq:survival-principal-tail} \mathbb{P}_x(\tau_D>t)=C_D(x)e^{-\lambda_1(D)t}(1+o(1)), \qquad C_D(x)>0.\qquad{(6)}\] In particular, \[\label{eq:lambda1-log-survival-limit} \lambda_1(D) =-\lim_{t\to\infty}\frac{1}{t}\log \mathbb{P}_x(\tau_D>t).\qquad{(7)}\]
Proof. The transition density \(p_D(t,x,y)\) and its strict positivity under the stated irreducibility assumptions are supplied by [14]. Thus \[T_t^D f(x)=\int_D p_D(t,x,y)f(y)\,\mathrm{d}y.\] Let \(p(t,x,y)\) denote the free transition density of the rectilinear stable process on the whole space \(\mathbb{R}^d\), i.e. before imposing the killing at \(\partial D\). Because the coordinates are independent, \[p(t,x,y)=\prod_{i=1}^d p_t^{(1D)}(y_i-x_i),\] where the one dimensional symmetric \(\alpha\)-stable density is \[p_t^{(1D)}(z) =\frac{1}{2\pi}\int_{\mathbb{R}}e^{\mathrm{i}\xi z}e^{-t|\xi|^\alpha}\,\mathrm{d}\xi =\frac{1}{\pi}\int_0^\infty e^{-t\xi^\alpha}\cos(\xi z)\,\mathrm{d}\xi.\] Since \(p_D(t,x,y)\le p(t,x,y)\), the stable scaling in [16] gives \[p(2t,x,x)\le C t^{-d/\alpha}.\] Using symmetry and the semigroup property, \[\iint_{D\times D}p_D(t,x,y)^2\,\mathrm{d}x\,\mathrm{d}y =\int_D p_D(2t,x,x)\,\mathrm{d}x \le |D|Ct^{-d/\alpha}<\infty.\] Thus \(T_t^D\) is Hilbert–Schmidt, hence compact on \(L^2(D)\). The spectral theorem for compact self-adjoint operators [17] gives ?? –?? . Because \(T_t^D\) is positivity improving, the Jentzsch–Krein–Rutman theorem for compact positive operators [18] gives simplicity and positivity of the leading eigenfunction.
Finally, \[S_x(t)=T_t^D1(x)=\int_Dp_D(t,x,y)\,\mathrm{d}y.\] Expanding \(1\) in the eigenbasis gives ?? . Since \(\phi_1>0\), the coefficient \[C_D(x):=\phi_1(x)\int_D\phi_1(y)\,\mathrm{d}y\] is strictly positive. The term with \(n=1\) dominates the expansion as \(t\to\infty\), which yields ?? and the logarithmic limit ?? . ◻
The probabilistic meaning of ?? is that the principal eigenvalue is the exponential decay rate of the survival probability. This is exactly the quantity that can be observed from simulated exit times; no eigenfunction needs to be computed.
Proposition 3 (WoC estimator for the principal eigenvalue). Let \(Q=[-1,1]^d\), and suppose the precomputed jump pool contains independent samples of the unit-cube exit pair \[\left(Y_{\sigma_Q},\sigma_Q\right), \qquad \sigma_Q:=\inf\{t>0:Y_t\notin Q\},\] where \(Y_t\) is the rectilinear stable process started at the origin. If the current WoC position is \(z\in D\) and \[r(z):=\operatorname{dist}_{\infty}(z,D^c),\] then the stable self-similarity and the strong Markov property imply the exact local rescaling law \[\label{eq:woc-local-rescaling} \left(X_{\tau_{z+rQ}},\tau_{z+rQ}\right) \stackrel{d}{=} \left(z+rY_{\sigma_Q},\,r^\alpha\sigma_Q\right), \qquad X_0=z.\qquad{(8)}\] Here \(\tau_{z+rQ}:=\inf\{t>0:X_t\notin z+rQ\}\) is the first exit time from the local cube. Consequently, a WoC path has accumulated time \[\tau_{\epsilon} =\sum_{j=0}^{N_\epsilon-1} r_j^\alpha\sigma_Q^{(j)},\] where \(r_j=r(Z_j)\) and the path stops when \(r_j\le\epsilon\) or when the sampled exit point leaves \(D\). If the unit-cube exit law is sampled exactly and \(\epsilon\downarrow0\), then \(\tau_{\epsilon}\) converges in distribution to \(\tau_D\). Hence independent WoC paths give the empirical survival curve \[\widehat S_{x,N}^{(\epsilon)}(t) =\frac{1}{N}\sum_{k=1}^N \mathbf{1}_{\{\tau_{\epsilon}^{(k)}>t\}},\] and, for each fixed fitting grid, the strong law of large numbers gives \[\widehat S_{x,N}^{(\epsilon)}(t) \xrightarrow[N\to\infty]{\mathrm{a.s.}} \mathbb{P}_x(\tau_{\epsilon}>t) \xrightarrow[\epsilon\downarrow0]{} \mathbb{P}_x(\tau_D>t).\] Combining this with ?? , the slope of the late-time log survival curve estimates \(\lambda_1(D)\): \[\log \widehat S_{x,N}^{(\epsilon)}(t) \approx -\lambda_1(D)t+b_x.\]
The rescaling ?? is the reason WoC is appropriate here. The local domain used by the algorithm is an \(L_\infty\) cube, and the rectilinear process has independent coordinate jumps and the scaling \(X_{ct}\stackrel{d}{=}c^{1/\alpha}X_t\); see the stable process references [13], [16]. The construction is the cube analogue of the walk-on-spheres principle of Muller [10], and it plays the same role for fractional exit problems as the isotropic fractional walk-on-spheres algorithm in [1].
For implementation, choose late-time grid points \(t_1<\cdots<t_m\) with nonzero empirical survival counts, and compute \[(\widehat\lambda_1,\widehat b_x) =\arg\min_{\ell,b} \sum_{q=1}^m \left(\log\widehat S_{x,N}^{(\epsilon)}(t_q)+\ell t_q-b\right)^2.\] The first component of the minimizer is the estimate \(\widehat\lambda_1\). The fitting window should be late enough that the first eigenmode dominates in ?? , but not so late that \(\widehat S\) is dominated by a small number of surviving samples.
Finally, this eigenvalue estimate is not merely diagnostic. If \(\lambda<0\) and \(\kappa=-\lambda>0\), the Helmholtz Feynman–Kac weight is \(\mathrm{e}^{\kappa\tau_D}\). From the tail identity \[\mathbb{E}_x\left[\mathrm{e}^{a\tau_D}\right] =1+a\int_0^\infty e^{at}\mathbb{P}_x(\tau_D>t)\,\mathrm{d}t,\qquad a>0,\] together with ?? , one obtains \[\mathbb{E}_x[e^{a\tau_D}]<\infty \quad\Longleftrightarrow\quad a<\lambda_1(D).\] This is the spectral form of the gaugeability threshold for Schrödinger-type perturbations of Markov processes [9], [12]. Thus the expectation in the Helmholtz representation is finite when \[\kappa=-\lambda<\lambda_1(D),\] and, for bounded boundary data, the Monte Carlo variance is finite under the stronger \(L^2\) condition \[2\kappa=2|\lambda|<\lambda_1(D).\] This is why the WoC survival profiler can be used before the Helmholtz solve: it estimates the largest safe exponential growth rate allowed by the domain.
4 Walk-on-Cubes↩︎
Based on the regime analysis in Section 3, we formulate a unified Walk-on-Cubes (WoC) Monte Carlo solver for the Dirichlet problem \(A_{\mathbf{x}}^{\alpha}u-\lambda u=0\) on \(D\subset\mathbb{R}^d\). The implementation uses a common cube-exit kernel together with regime-dependent weighting: for \(\lambda=0\), coordinatewise \(\alpha\)-stable pure-jump scaling determines the spatial transport; for \(\lambda>0\), a Duffin extension adds an auxiliary coordinate \(w\) and the boundary payoff includes the oscillatory factor; for \(\lambda<0\), a Feynman–Kac time enhancement accumulates survival time through \(\Delta\tau=r^{\alpha}\tau_{\mathrm{std}}\) and applies the exponential reweighting \(\exp(-\lambda\tau_{\mathrm{total}})\) under the subcritical gauge constraint. In this way, the Laplace, Yukawa, and Helmholtz cases are treated within a single implementation framework that is well suited to parallel computation on modern GPUs.
4.1 Walk-on-Cubes Algorithm↩︎
We consider the numerical solution of the nonisotropic fractional-order partial differential equation with a potential term: \[A_{\mathbf{x}}^\alpha u(\mathbf{x}) - \lambda u(\mathbf{x}) = 0, \quad \mathbf{x}\in D \subset \mathbb{R}^d.\] subject to the Dirichlet exterior boundary condition \(u(\mathbf{x}) = g(\mathbf{x})\) for \(\mathbf{x}\in D^c\). Here, \(\alpha \in (0, 2)\), and \(A_{\mathbf{x}}^\alpha\) is the nonisotropic fractional Laplacian, defined as the sum of one-dimensional fractional operators: \[A_{\mathbf{x}}^\alpha = \sum_{i=1}^d -(-\Delta_{x_i})^{\alpha/2}.\]
4.2 Mathematical Formulation of the Algorithm↩︎
To implement this problem efficiently on a parallel GPU architecture, we separate the regime-dependent mathematical input from the online path update used in the WoC solver. The role of the potential parameter \(\lambda\) depends on the regime. For \(\lambda>0\), the Duffin extension introduces an auxiliary coordinate \(w\) and the lifted field satisfies \[\mathcal{A}_{\mathbf{x},w}^{\alpha}U(\mathbf{x},w)=0,\quad (\mathbf{x},w)\in D\times\mathbb{R},\] where \(\mathcal{A}_{\mathbf{x},w}^{\alpha}:=A_{\mathbf{x}}^{\alpha}-(-\Delta_w)^{\alpha/2}\) and \(U(\mathbf{x},w)=u(\mathbf{x})f_\lambda(w)\). For \(\lambda<0\), no spatial lifting is used; instead, the solver accumulates survival time and applies the exponential reweighting \(\exp(-\lambda\tau_{\mathrm{total}})\). For \(\lambda=0\), neither correction is needed. The main implementation issue is that a direct simulation of fractional microsteps inside each CUDA thread produces severe warp divergence: the number of microsteps required before leaving a local cube is random and can vary substantially from one thread to another. To reduce this cost, we move this variable length part of the computation to an offline preprocessing stage and precompute a Jump Pool \(\mathcal{P}\) of standardized exits from the unit \(L_\infty\) cube \(Q=[-1,1]^d\). The online solver then uses only a distance query, a pool lookup, and a local rescaling at each WoC update. The infinitesimal generator of the nonisotropic operator \(A_{\mathbf{x}}^\alpha\) corresponds to a \(d\)-dimensional Lévy process with independent symmetric \(\alpha\)-stable components. The standard exit vectors are sampled empirically by simulating discrete paths from the origin. By the self-similarity of the \(\alpha\)-stable process (\(\mathbf{X}_{ct} \stackrel{d}{=} c^{1/\alpha} \mathbf{X}_t\)), the position update at step \(n\) for a time step \(\Delta t\) is given by \[\mathbf{X}^{(n+1)} = \mathbf{X}^{(n)} + (\Delta t)^{1/\alpha} \boldsymbol{\xi}^{(n)},\] where \(\boldsymbol{\xi}^{(n)} = (\xi_1^{(n)}, \dots, \xi_d^{(n)})\) consists of independent standard symmetric \(\alpha\)-stable random variables, in agreement with the nonisotropic sum formulation.
To generate each component \(\xi_i\), we use the Chambers-Mallows-Stuck (CMS) method [19]. We sample a uniform angle \(V\) and an independent standard exponential variable \(W\): \[V \sim
\mathcal{U}\left(-\frac{\pi}{2}, \frac{\pi}{2}\right), \quad W \sim \text{Exp}(1).\] The exact \(\alpha\)-stable increment is computed by \[\xi_i = \frac{\sin(\alpha V)}{(\cos
V)^{1/\alpha}} \left( \frac{\cos(V - \alpha V)}{W} \right)^{\frac{1-\alpha}{\alpha}}.\] A path halts at step \(N\) when it breaches the boundary (\(\mathbf{X}^{(N)} \notin Q\)). This
terminal state yields one standard space–time sample \((\Delta \mathbf{X}_{\mathrm{std}},\tau_{\mathrm{std}})\), and the resulting pool is \[\mathcal{P} = \left\{ \left(\Delta
\mathbf{X}_{\mathrm{std}}^{(k)},\tau_{\mathrm{std}}^{(k)}\right) \right\}_{k=1}^{N_{\mathrm{pool}}}.\] After the jump pool has been constructed, the physical domain \(D\) is discretized on a uniform grid, and
independent curandState objects are assigned to the grid points so that the random number streams are uncorrelated. The online WoC solver is summarized in Algorithm 1. In the algorithm, \(\mathbf{Z}^{(n)}\) denotes the full path state; its spatial component is \(\mathbf{x}^{(n)}\). Thus \(\mathbf{Z}^{(n)}=(\mathbf{x}^{(n)},w^{(n)})\) in the Yukawa
regime \(\lambda>0\), while \(\mathbf{Z}^{(n)}=\mathbf{x}^{(n)}\) in the Laplace and Helmholtz regimes.
The regime-dependent terminal payoff is \[\label{eq:woc-payoff} \text{Payoff} = \begin{cases} g(\mathbf{x}_{\mathrm{exit}}), & \lambda = 0,\\ g(\mathbf{x}_{\mathrm{exit}}) \cdot f_\lambda(w_{\mathrm{exit}}), & \lambda > 0,\\ g(\mathbf{x}_{\mathrm{exit}}) \cdot \exp\!\big(-\lambda\tau_{\mathrm{total}}\big), & \lambda < 0. \end{cases}\tag{14}\] In the Helmholtz case, this uses the Feynman–Kac exponential weight, which is amplifying because \(\lambda<0\), and it is valid in the subcritical gauge regime [12] so that the relevant exponential moments remain finite.
5 Numerical Considerations and Examples↩︎
Variance Blow-up from Boundary Conditions↩︎
Care is required when defining the exterior Dirichlet boundary condition \(g(\mathbf{x})\). If \(g\) grows too rapidly with distance, the numerical solution may exhibit large irregular fluctuations. This is caused by the heavy-tailed nature of the fractional process. The domain exterior Poisson kernel decays according to a power law (\(\propto |\mathbf{x}|^{-d-\alpha}\)) [20], so particles still have a non-negligible probability of making very large jumps into the far exterior. Coupling these heavy-tailed exit locations with an unbounded rapidly growing \(g(\mathbf{x})\) can make the empirical variance of the Monte Carlo estimator infinite or extremely large, thereby destroying the smoothness of the computed expectation surface.
To validate the robustness of the generalized solver, we first report two Yukawa Green function manufactured benchmarks, then a one-dimensional Helmholtz benchmark with an exact Fourier mode solution, and finally four additional test cases spanning multiply connected, non-convex, and curvilinear domains in both 2D and 3D.
Green Function Manufactured Benchmark↩︎
The Green function benchmarks use an exact solution generated from the full-space resolvent of the positive operator \[L^\alpha:=-A^\alpha=\sum_{i=1}^d(-\Delta_{x_i})^{\alpha/2}.\] For \(\lambda>0\), using the Fourier multiplier characterization of the fractional Laplacian [7], the full-space resolvent kernel is \[G_{\lambda,\alpha}(\mathbf{x}) =\mathcal{F}^{-1}\!\left[ \frac{1}{\lambda+\sum_{i=1}^d|\xi_i|^\alpha} \right](\mathbf{x}).\] Here \(\mathcal{F}^{-1}\) denotes the inverse Fourier transform. Equivalently, using the \(q\)-resolvent/potential operator representation for Lévy semigroups [13] and the rectilinear transition density product formula [14], if \(p_t\) denotes the transition density of the rectilinear stable process generated by \(A^\alpha\), then \[G_{\lambda,\alpha}(\mathbf{x})=\int_0^\infty e^{-\lambda t}p_t(\mathbf{x})\,\mathrm{d}t.\] For the broader anisotropic coordinate-sum setting, potential kernels and Green function estimates are treated by Bogdan and Sztonyk [21].
By construction, \[(L^\alpha+\lambda)G_{\lambda,\alpha}=\delta_0 \quadin \mathcal{D}'(\mathbb{R}^d),\] because the Fourier transform of the left-hand side is \[\left(\lambda+\sum_{i=1}^d|\xi_i|^\alpha\right) \frac{1}{\lambda+\sum_{i=1}^d|\xi_i|^\alpha}=1.\] Now place the pole at a point \(\mathbf{a}\notin \overline{D}\) and set \[u_{\mathrm{exact}}(\mathbf{x}):=G_{\lambda,\alpha}(\mathbf{x}-\mathbf{a}).\] Then \((L^\alpha+\lambda)u_{\mathrm{exact}}=0\) in \(D\), or equivalently \[A^\alpha u_{\mathrm{exact}}=\lambda u_{\mathrm{exact}}\quadin D.\] Thus this gives an exact Yukawa regime benchmark for 3 with \(\lambda>0\).
The nonlocal boundary condition is imposed on the whole exterior \(D^c\), not only on \(\partial D\); this is the standard exterior Dirichlet formulation for fractional Laplace problems [22], [23]. Therefore we prescribe \[g(\mathbf{x})=u_{\mathrm{exact}}(\mathbf{x})=G_{\lambda,\alpha}(\mathbf{x}-\mathbf{a}),\qquad \mathbf{x}\in D^c.\]
For one-dimensional tests the exact value can be evaluated by the Fourier integral \[G_{\lambda,\alpha}(\mathbf{x}-\mathbf{a}) =\frac{1}{\pi}\int_0^\infty \frac{\cos(\xi(x_1-a_1))}{\lambda+\xi^\alpha}\,\mathrm{d}\xi, \qquad \mathbf{x}=(x_1),\quad \mathbf{a}=(a_1).\] For the coordinate-sum operator in several dimensions, the time resolvent formula is convenient because the free transition density factorizes coordinatewise: \[p_t(\mathbf{x})=\prod_{i=1}^d p_t^{(1D)}(x_i).\] Here the one-dimensional factor is the symmetric \(\alpha\)-stable transition density defined by the Fourier representation \[p_t^{(1D)}(z) =\frac{1}{2\pi}\int_{\mathbb{R}} e^{i\xi z}e^{-t|\xi|^\alpha}\,\mathrm{d}\xi =\frac{1}{\pi}\int_0^\infty e^{-t\xi^\alpha}\cos(\xi z)\,\mathrm{d}\xi.\] In practice, the pole is kept a positive distance outside \(D\) so that the singularity of \(G_{\lambda,\alpha}\) is never sampled inside the computational domain. Moderate values such as \(\alpha>1\) and \(\lambda\) of order \(10^{-1}\)–\(1\) keep the Green function well behaved and make the benchmark a clean test of the WoC implementation.
Case 11: One-Dimensional Yukawa Green Comparison↩︎
The first Green function comparison case validates the positive Yukawa branch \(\lambda>0\) for the one-dimensional operator \[A^\alpha=-(-\Delta)^{\alpha/2}.\] The computational domain is \(D=(-1,1)\), with \(\mathbf{x}=(x_1)\), and the benchmark equation is \[A^\alpha u=\lambda u, \qquadequivalently\qquad \big((-\Delta)^{\alpha/2}+\lambda\big)u=0 \quadin D.\] The exact solution is generated by a one-dimensional fractional Yukawa Green function with a pole placed outside the interval: \[\begin{align} u_{\mathrm{exact}}(\mathbf{x}) &=G_{\lambda,\alpha}(\mathbf{x}-\mathbf{a}) =G_{\lambda,\alpha}(x_1-a_1),\\ G_{\lambda,\alpha}(r) &=\frac{1}{\pi}\int_0^\infty \frac{\cos(\xi r)}{\lambda+\xi^\alpha}\,\mathrm{d}\xi, \qquad \mathbf{a}=(a_1),\quad a_1=2.5. \end{align}\] Since \(\mathbf{a}\notin[-1,1]\), the singularity never enters the computational domain and the distributional resolvent identity gives \[\big((-\Delta_{x_1})^{\alpha/2}+\lambda\big)G_{\lambda,\alpha}(x_1-a_1)=0, \qquad \mathbf{x}=(x_1)\in(-1,1).\] The exterior Dirichlet data are therefore fixed by the same full-space manufactured solution, \[g(\mathbf{x})=G_{\lambda,\alpha}(x_1-2.5), \qquad \mathbf{x}=(x_1)\in D^c.\]
| Parameter | Value |
|---|---|
| Domain | \((-1,1)\) |
| Fractional order | \(\alpha=1.5\) |
| Yukawa parameter | \(\lambda=0.1\) |
| Auxiliary frequency | \(\omega=\lambda^{1/\alpha}\approx0.215443\) |
| Pole location | \(\mathbf{a}=(2.5)\) |
| Iteration shots | \(1{,}000{,}000\) |
| Jump pool size | \(1{,}000{,}000\) |
| Grid size | \(64\) |
| Random seed | \(20260505\) |
| Time step | \(10^{-4}\) |
| Exit tolerance | \(10^{-5}\) |
| Maximum steps | \(20000\) |
| Quadrature cutoff and panels | \(300.0,\;6000\) |
| Green function table size and range | \(2049,\;|r|\le40\) |
Figures 2–3 show that the Monte Carlo curve follows the analytic Yukawa Green profile, and the pointwise error remains on the same scale as the sampling uncertainty.
Case 12: Two-Dimensional Separable Yukawa Green Comparison↩︎
The second Green function comparison checks the coordinate-sum operator \[A_{\mathbf{x}}^\alpha=-(-\Delta_{x_1})^{\alpha/2}-(-\Delta_{x_2})^{\alpha/2}\] on the square \(D=(-1,1)^2\), where \(\mathbf{x}=(x_1,x_2)\). The equation is \[A_{\mathbf{x}}^\alpha u=\lambda u, \qquadequivalently\qquad \big((-\Delta_{x_1})^{\alpha/2}+(-\Delta_{x_2})^{\alpha/2}+\lambda\big)u=0 \quadin D.\] Here the manufactured solution is separable: \[u_{\mathrm{exact}}(\mathbf{x}) =G_{\mu,\alpha}(x_1-a_1)\,G_{\nu,\alpha}(x_2-a_2), \qquad \mu+\nu=\lambda.\] In the reported run \[\lambda=0.1,\qquad \mu=\nu=0.05,\qquad \mathbf{a}=(a_1,a_2)=(2.5,2.25).\] This construction is matched to the nonisotropic coordinate-sum operator. It is not the radial full-space Green kernel of the isotropic fractional Laplacian; rather, it uses the one-dimensional resolvent identity in each coordinate. Since both poles lie outside the corresponding coordinate intervals, \[\begin{align} (-\Delta_{x_1})^{\alpha/2}G_{\mu,\alpha}(x_1-a_1)=-\mu G_{\mu,\alpha}(x_1-a_1), \\ (-\Delta_{x_2})^{\alpha/2}G_{\nu,\alpha}(x_2-a_2)=-\nu G_{\nu,\alpha}(x_2-a_2) \end{align}\] for \(\mathbf{x}\in D\). Therefore \[\left((-\Delta_{x_1})^{\alpha/2}+(-\Delta_{x_2})^{\alpha/2}\right)u_{\mathrm{exact}} =-(\mu+\nu)u_{\mathrm{exact}} =-\lambda u_{\mathrm{exact}},\] or equivalently \(A_{\mathbf{x}}^\alpha u_{\mathrm{exact}}=\lambda u_{\mathrm{exact}}\) in \(D\). The exterior data are again prescribed on the whole complement: \[g(\mathbf{x})=u_{\mathrm{exact}}(\mathbf{x}), \qquad \mathbf{x}\in D^c.\]
| Parameter | Value |
|---|---|
| Domain | \((-1,1)^2\) |
| Fractional order | \(\alpha=1.5\) |
| Yukawa parameter | \(\lambda=0.1\) |
| Separable parameters | \(\mu=\nu=0.05\) |
| Auxiliary frequency | \(\omega=\lambda^{1/\alpha}\approx0.215443\) |
| Pole location | \(\mathbf{a}=(2.5,2.25)\) |
| Iteration shots | \(1{,}000{,}000\) |
| Jump pool size | \(1{,}000{,}000\) |
| Grid size | \(64\times64\) |
| Random seed | \(20260505\) |
| Time step | \(10^{-4}\) |
| Exit tolerance | \(10^{-5}\) |
| Maximum steps | \(20000\) |
| Quadrature cutoff and panels | \(300.0,\;6000\) |
| Green function table size and range | \(2049,\;|r|\le60\) |
Figures 4–5 show the expected agreement between the WoC estimator and the separable manufactured solution. The absolute error plot stays at the scale predicted by the Monte Carlo standard error, which supports the use of the Green function construction as a benchmark for the Yukawa regime of the nonisotropic operator.
5.1 One-Dimensional Helmholtz Walk-on-Intervals Validation↩︎
We first validate the one-dimensional solver against the exact Fourier mode derived in 6 . For the generator \[A^\alpha=-(-\Delta)^{\alpha/2},\] the manufactured exact solution \[u_{\mathrm{exact}}(\mathbf{x})=\cos(kx_1),\qquad \mathbf{x}=(x_1),\] satisfies \[A^\alpha u_{\mathrm{exact}}(\mathbf{x})=-k^\alpha u_{\mathrm{exact}}(\mathbf{x}),\] so the corresponding Helmholtz parameter is \(\lambda=-k^\alpha<0\). The computational domain is \(D=(-1,1)\), and the exterior Dirichlet data are taken from the exact solution: \[g(\mathbf{x})=u_{\mathrm{exact}}(\mathbf{x})=\cos(kx_1),\qquad \mathbf{x}\in D^c.\]
Before the production runs, the solver performs the same two stability checks used in the main Helmholtz code path. For the exterior data, the spatial \(L^2\) criterion checks the polynomial growth of the exterior boundary data against the fractional exit tail: if \(|g(x)|\le C(1+|x|)^p\) at infinity, then the one-dimensional condition \(p<\alpha/2\) is sufficient for the sampled exterior contribution to have a finite second moment. The detected growth rate was \[p\approx0.401810<\alpha/2=0.750000,\] so the spatial \(L^2\) criterion passed. For the time gauge, the estimated principal eigenvalue and amplification rate were \[\lambda_1\approx 1.611357,\qquad -\lambda=0.353553,\] which places the test in the subcritical regime \(-\lambda<\lambda_1\).
| Parameter | Value |
|---|---|
| Dimension | \(1\) |
| Domain | \([-1,1]\) |
| \(L\) | \(1.0\) |
| Boundary condition | \(u(\x)=\cos(kx_1)\) |
| Exact solution | \(u(\x)=\cos(kx_1)\) |
| Fractional order | \(\alpha=1.5\) |
| Iteration shots | \(1{,}000{,}000\) |
| Jump pool size | \(50{,}000\) |
| Time step | \(3\times 10^{-4}\) |
| \(\epsilon\) | \(10^{-5}\) |
| Maximum steps | \(20{,}000\) |
| Grid size | \(64\) |
The lowest-ratio run corresponds to 000_helmholtz_1m_solution_d1; its shared configuration is listed in Table 3. The parameters are \[\frac{|\lambda|}{\lambda_1}=0.15,\qquad
\lambda=-0.2417036195774674,\qquad
k=0.38802118360634363,\] with random seed 20269505. The observed errors were \[\begin{align}
\|u_{\mathrm{MC}}-u_{\mathrm{exact}}\|_{\infty}=0.001297267461,
\\
\|u_{\mathrm{MC}}-u_{\mathrm{exact}}\|_{2}=0.0005009832629540645,
\end{align}\] with mean standard error \(0.0003333879742063492\) and zero maximum step hits.
The highest-ratio run corresponds to 003_helmholtz_1m_solution_d1, with \[\frac{|\lambda|}{\lambda_1}=0.65,\qquad
\lambda=-1.0473823515023588,\qquad
k=1.0313438926461413,\] and random seed 20269508. Its errors were \[\begin{align}
\|u_{\mathrm{MC}}-u_{\mathrm{exact}}\|_{\infty}=0.052417647075,
\\
\|u_{\mathrm{MC}}-u_{\mathrm{exact}}\|_{2}=0.014328343996343449,
\end{align}\] with mean standard error \(0.008218872176079367\) and zero maximum step hits.
Comparing the two runs shows the expected deterioration as \(|\lambda|/\lambda_1\) increases toward the critical threshold. The low-ratio case remains highly accurate, while the higher-ratio case exhibits visibly larger bias and variance, consistent with the amplifying Feynman–Kac weight in the Helmholtz regime.
Systematic Dataset Diagnostics↩︎
To complement the individual benchmark plots, we also ran systematic experiments, denoted Experiments 1–3. These experiments cover Monte Carlo convergence, jump pool convergence, and principal eigenvalue survival tail validation.
| Test | Quantity varied | Main numerical conclusion |
|---|---|---|
| Exp1 | \(N_{\mathrm{shots}}=10^3,\ldots,10^6\) for the 1D Yukawa Green benchmark | Mean \(L^2\) error decreases from \(7.234846\times10^{-3}\) to \(2.973546\times10^{-4}\); the fitted log–log slope is \(-0.464887\), close to the expected Monte Carlo rate. |
| Exp2 | Jump pool size \(N_{\mathrm{pool}}=10^3,\ldots,10^6\) for 1D Yukawa, 2D Yukawa, and 1D Helmholtz tests | The finest pool errors are \(2.556641\times10^{-4}\), \(2.708984\times10^{-4}\), and \(1.237600\times10^{-3}\), respectively, showing stabilization once the pool is sufficiently large. |
| Exp3 | Survival tail regression for intervals, rectangles, and hypercubes | The median relative error in \(\lambda_1\) is \(1.434713\times10^{-2}\), the maximum relative error is \(5.962481\times10^{-2}\), and the median tail regression \(R^2\) is \(0.999956\). |
The Monte Carlo convergence test in Figure 8 supports the expected square-root sampling trend. The fitted slope is slightly shallower than \(-1/2\), which is consistent with finite-grid effects, tabulated Green function error, and jump pool discretization being present together with sampling noise. The jump pool experiment in the same figure shows that very small pools introduce a visible additional error, especially in the Green function tests, while larger pools bring the error down to the sampling level.
The eigenvalue survival tail validation in Figure 9 checks the numerical mechanism used to estimate the principal threshold in the Helmholtz regime. The tests include interval scaling in one dimension, additivity on product rectangles, and dimension scaling on hypercubes. The high tail regression \(R^2\) values and small relative errors indicate that the survival profiler gives a reliable estimate of the exponential decay rate controlling the Feynman–Kac admissibility threshold.
For the following four Laplace test cases, the common numerical parameters are \[\begin{align} \alpha&=1.5, & \lambda&=0, & N_{\mathrm{shots}}&=10^8\;\text{per grid point},\\ \Delta t_{\mathrm{micro}}&=10^{-4}, & \epsilon&=10^{-5}. \end{align}\]
Case 1: 2D Doubly Connected Domain↩︎
This case investigates a non-trivial topological domain containing an internal obstacle.
Domain: \(D = [-1, 1]^2 \setminus [-0.5, 0.5]^2\) (A square with a square hole).
Boundary Condition: \(g(\mathbf{x}) = \sin(\pi x_1) \cos(\pi x_2)\) for \(\mathbf{x}=(x_1,x_2) \in D^c\).
Case 2: 3D Asymmetric Non-convex Domain↩︎
This case tests the algorithm’s capability to handle sharp re-entrant corners in 3D without suffering from meshing singularities typical in Finite Element Methods (FEM).
Domain: \(D = [-1, 1]^3 \setminus [0, 1]^3\) (A cube missing one octant).
Boundary Condition: \(g(\mathbf{x}) = x_1 + x_2 - x_3\) for \(\mathbf{x}=(x_1,x_2,x_3) \in D^c\).
Case 3: 2D Unit Disk↩︎
Simulating the WoC algorithm on curvilinear domains requires careful formulation of the maximal inscribed \(L_\infty\) cube to maintain \(O(1)\) jump efficiency.
Domain: \(D = \{ \mathbf{x}=(x_1,x_2) \mid x_1^2 + x_2^2 < 1 \}\) (Unit Disk).
Boundary Condition: \(g(\mathbf{x}) = \cos(x_1) + \sin(x_2)\) for \(\mathbf{x}\in D^c\).
Case 4: 3D Spherical Shell↩︎
This case extends the geometric distance derivation to 3D and tests a bounded exponential boundary condition.
Domain: \(D = \{ \mathbf{x}=(x_1,x_2,x_3) \mid 0.5 < |\mathbf{x}| < 1.0 \}\) (Spherical Shell).
Boundary Condition: \(g(\mathbf{x}) = \exp(-|\mathbf{x}|^2)\) for \(\mathbf{x}\in D^c\).
