Walking on Heat Stars for Parabolic Heat Equations with Neumann Boundary Conditions

Abstract: Monte Carlo methods have proven highly effective for elliptic partial differential equations through algorithms such as Walk on Spheres and Walk on Stars, which evaluate solutions at individual points without volumetric meshing or global linear solves. Extending these methods to the transient regime has remained an open challenge: parabolic equations couple space and time through an anisotropic scaling, requiring joint sampling of spatial displacements and backward time steps whose distribution was not previously available in a unified, exact form. We present Walk on Heat Stars, a grid-free Monte Carlo solver that closes this gap by extending the boundary integral framework of Walk on Stars to the parabolic setting. Our method introduces a non-cylindrical boundary integral formulation that accommodates the time-varying domains induced by heat-ball sampling. The heat ball geometry is parameterized by a logarithmic time coordinate and a spatial direction, revealing that the double-layer kernel factorizes into independent Gamma and uniform components. This factorization enables exact directional importance sampling of the recursive next walk position, the Neumann flux contribution, and the volumetric source term. We further derive a decoupled gradient estimator that expresses spatial derivatives as weighted boundary integrals of the solution, requiring no recursion on the gradient, and adapt a heteroscedastic regression-based denoiser to the space-time domain for variance reduction. We validate our method on analytical solutions across a range of geometries and spatial frequencies, confirm convergence at the expected Monte Carlo rate, and demonstrate practical applicability on heat sink and cooling scenes with mixed or pure Neumann boundary conditions.

• The paper introduces Walk on Heat Stars, a grid-free Monte Carlo method for transient heat equations that supports mixed Dirichlet–Neumann boundaries, volumetric sources, and non-convex geometries.
• The method factorizes the heat-sphere double-layer kernel into Gamma-distributed logarithmic time and uniform spatial-direction samples, enabling unit-weight importance sampling with standard ray-casting operations.
• Experiments confirm expected Monte Carlo convergence and agreement with analytical and finite-element solutions, while revealing higher errors for high-frequency fields, narrow corridors, and gradient estimation without denoising.

Problem and contribution

This paper presents Walk on Heat Stars (WoHSt), a grid-free Monte Carlo solver for the transient heat equation with mixed Dirichlet–Neumann boundary conditions, volumetric sources, and complex geometries. The method extends the boundary integral framework of Walk on Stars (WoSt) to the parabolic setting, where space and time are coupled through anisotropic scaling: a time increment $\Delta t$ corresponds to a spatial displacement of order $\sqrt{\Delta t}$. Prior transient Monte Carlo solvers—Walk on Moving Spheres (WoMS) and the heat-ball walk of Deaconu and Herrmann—handle only pure Dirichlet problems, while the global random walk on spheres of Sabelfeld likewise excludes Neumann conditions. WoHSt is, to the authors' knowledge, the first method to simultaneously support the transient regime, Neumann boundary conditions, source terms, and non-convex geometry, matching the feature set of WoSt in the elliptic case.

The paper's central technical insight is geometric: the heat-sphere surface, parameterized by a logarithmic time coordinate $u = \log(c/s)$ and a spatial direction $\omega$, yields a double-layer kernel that factorizes exactly into independent Gamma and uniform distributions. This factorization enables unit-weight importance sampling throughout the estimator, including on truncated "heat star" boundaries where a radial trace intersects a Neumann wall before reaching the heat-sphere surface.

Non-cylindrical boundary integral formulation

The method departs from the standard cylindrical boundary integral representation of Costabel by adopting a non-cylindrical formulation for domains $\Omega_s$ that evolve in backward time. The authors derive this representation self-contained from Green's second identity, the Reynolds transport theorem, and the adjoint heat kernel identity, recovering the moving-boundary correction term $\int G\,\langle V, \vec{n}\rangle\,u$ that vanishes for stationary domains. The formulation is stated under smoothness assumptions; the more general Sobolev-space treatment of Brügger et al. is cited but not reproduced.

The heat ball $\mathcal{H}(x_0, t_0, c)$ is defined as the super-level set of the fundamental solution $\Phi$ above the threshold $\tau(c) = (4\pi c)^{-n/2}$. Its spatial cross-sections have radius $r(s) = \sqrt{2n\,s\log(c/s)}$, growing from zero, attaining a maximum $\sqrt{\Delta t}$0 at $\sqrt{\Delta t}$1, and contracting back to zero—hence a spindle shape in space-time. The Green's function $\sqrt{\Delta t}$2 vanishes on the heat-sphere boundary, eliminating the single-layer term there and yielding the non-uniform mean value property for caloric functions. The initial-condition term vanishes because the earliest slice of the heat ball degenerates to a point.

Exact importance sampling via kernel factorization

The core sampling result (Proposition 2 in the paper) establishes that the double-layer kernel $\sqrt{\Delta t}$3 on the heat-star boundary factorizes as

$\sqrt{\Delta t}$4

The proof exploits the solid-angle identity $\sqrt{\Delta t}$5 to cancel the boundary orientation, after which all geometry collapses into the Gamma density over the logarithmic diffusion time. The consequence is a two-stage walk step: draw $\sqrt{\Delta t}$6 and $\sqrt{\Delta t}$7 uniformly, then trace the radial curve $\sqrt{\Delta t}$8 to the boundary. If the trace strikes the Neumann wall first (Neumann branch), the spatial step is clamped to the visibility distance and the time step is recovered by inverting the radial trace relation; otherwise the walk lands on the heat-sphere surface. In either case the importance weight is exactly one, requiring only standard ray-casting primitives.

The same machinery extends to the Neumann flux and volumetric source terms. Flux sampling reuses the Gamma–Uniform draw to determine the spatial support and samples uniformly on the visible Neumann patch, reusing the ball–boundary intersection primitive from elliptic WoSt unchanged. Source sampling applies a heat-spherical change of variables under which the Green's function depends only on the radial coordinate $\sqrt{\Delta t}$9, the volume element factorizes into independent components, and the integral of $u = \log(c/s)$0 over the heat ball evaluates analytically to $u = \log(c/s)$1, so no explicit kernel evaluation is needed. The radial coordinate is drawn by rejection with acceptance rate $u = \log(c/s)$2.

The authors verify numerically that their sampling scheme is indistinguishable from that of Deaconu and Herrmann in dimensions $u = \log(c/s)$3 and $u = \log(c/s)$4, confirming that the classical heat-ball walk is recovered as the special case of pure Dirichlet data with no source. Their Gamma parameterization additionally unifies all dimensions under a single distribution, whereas prior work treated the logarithmic coordinate as an implementation detail.

Gradient estimation

The paper derives a decoupled spatial gradient estimator from the volumetric heat mean value property. Differentiating with respect to the query point and applying the divergence theorem on each spatial cross-section expresses $u = \log(c/s)$5 as a weighted boundary integral of $u = \log(c/s)$6 alone—no recursion on the gradient is required. A structural difference from the Laplace case arises because the mean value weight $u = \log(c/s)$7 is non-constant: its gradient contributes an additional volume term that is absent in the elliptic formula. Both surface and volume terms share a Gamma-distributed time marginal $u = \log(c/s)$8 with analytically precomputed normalization constants.

The authors are explicit that this estimator is a proof of concept. The $u = \log(c/s)$9 factor amplifies contributions from small backward-time steps where the Gamma distribution has a heavy tail, producing variance one to two orders of magnitude above the value estimator in their tests. Raw estimates are reported as too noisy for direct use; readable results require the heteroscedastic denoiser. A finite-difference alternative based on Walk on Heat Spheres with step length $\omega$0 performed worse, supporting the decoupled formulation despite its variance limitation.

Variance reduction

For denoising, the paper adapts the heteroscedastic regression technique of Bao et al. to the space-time domain: a SIREN network jointly regresses the mean and log-variance of the raw WoHSt estimates under the $\omega$1-NLL loss with $\omega$2, with the sole modification of including time $\omega$3 in the network input. Denoised solution fields closely match ground truth on configurations where raw estimates exhibit substantial variance.

Experimental evaluation

Validation uses separable analytical solutions $\omega$4 with spatial frequencies $\omega$5 on three geometries of varying boundary complexity (a tool, a gear-shaped domain, and a spot with a narrow corridor), evaluated over a $\omega$6 grid at five time instants under a fixed 500-second compute budget. Two consistent trends emerge. Error increases with spatial frequency, as expected under a fixed budget. More notably, the pure-Dirichlet Walk on Heat Spheres baseline consistently achieves lower MSE than the full WoHSt estimator—by roughly an order of magnitude in several configurations (e.g., $\omega$7 versus $\omega$8 on the ball geometry at $\omega$9)—reflecting the simpler estimator structure without Neumann or silhouette-ray queries. On the spot geometry with a narrow corridor, the WoHSt error at $\Omega_s$0 reaches approximately $\Omega_s$1, indicating that highly constrained geometries with high-frequency solutions remain challenging under realistic budgets. Convergence curves at the hardest frequency decay at the expected Monte Carlo rate across all geometries and time instants, confirming statistical well-behavedness of the importance sampling on non-convex domains. Temporal accuracy at fixed spatial points yields relative MSE between $\Omega_s$2 and $\Omega_s$3 depending on geometry and location.

Application experiments compare WoHSt against a finite element reference (FastTetWild meshing with deal.II) on a 3D heat sink with mixed boundary conditions, showing convergence to the FEM solution across time instants. A turbine housing with pure Neumann conditions and constant surface heat flux is simulated with only 16 walks per pixel, capturing the cooling process at three sensor locations; the authors note that pure-Neumann problems are well-posed here given the initial temperature constraint, so no Tikhonov regularization is required.

Limitations and open questions

Two limitations are conceded directly. First, Robin boundary conditions—which model convective heat exchange and are common in heat transfer—are not supported; extending the Brakhage–Werner trick and adaptive radius selection of elliptic Walkin' Robin to the parabolic setting remains open. Second, the gradient estimator's variance, driven by the $\Omega_s$4 singularity of the weight gradient, precludes production use without denoising; the robust derivative estimators developed for the elliptic case via Robin reformulation do not yet have parabolic counterparts. The experimental results additionally suggest that accuracy on narrow-corridor geometries at high spatial frequency degrades substantially, and the paper does not quantify the overhead of the silhouette-distance and visibility queries that the heat-star construction requires relative to pure heat-ball walks.

Conclusion

Walk on Heat Stars extends grid-free Monte Carlo PDE solving to transient heat conduction with mixed boundary conditions by combining a non-cylindrical boundary integral formulation with an exact Gamma–Uniform factorization of the heat-sphere double-layer kernel. The construction achieves unit-weight importance sampling on truncated heat-star boundaries using only standard ray-casting primitives, generalizes the classical heat-ball walk, and is validated against analytical solutions and finite element references on engineering geometries. The principal open problems are Robin boundary support and variance reduction for derivative estimation in the parabolic regime.