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

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.