Mesh‐free Monte Carlo method for electrostatic problems with floating potentials

Abstract

Numerical simulation plays a crucial role in the analysis and design of power equipment, such as lightning protection devices, which may become inefficient using traditional grid‐based methods when handling complex geometries of large problems. The authors propose a grid‐free Monte Carlo method to handle electrostatic problems of complex geometry for both the interior and exterior domains, which is governed by the Poisson equation with a floating potential boundary condition that is neither a pure Dirichlet nor a Neumann condition. The potential and gradient at any given point can be expressed in terms of integral equations, which can be estimated recursively within the walk‐on‐sphere algorithm. Numerical examples have been demonstrated, including the evaluation of the mutual capacitance matrix of multi‐conductor structures and lighting striking near real fractal trees. The proposed method shows advantages in terms of geometric flexibility and robustness, output sensitivity, and parallelism, which may become a candidate for game‐changing numerical methods and exhibit great potential applications in high‐voltage engineering.

INTRODUCTION

Electrostatic problems are subject to the solution of the Poisson equation in an open domain with specific boundary conditions. This kind of problem is often solved using a finite-dimensional approximation. In general, one needs to first build a mesh of the domain that can be used to define basis functions and then compute the solution that satisfies the Poisson equation in a weak sense. Even recent meshless finite element methods (FEMs) still require a careful arrangement of nodes within the domain to discretise the interior of the domain [1, 2]. Consequently, errors appear when one approximates the geometries and underlying physics using basis functions. Moreover, the discretisation of domains or boundaries can be the most expensive procedure or even fail when the geometry is complicated. In addition, comprehensive pre-processing is required for non-standard geometric data, such as point clouds from laser scanning for digital twin applications. For open-domain problems, the domain of interest often needs to be truncated when FEM are employed. This would introduce an error if we simply ground the outermost boundaries far from the region of interest because the potential, in fact, approaches zero at infinity; alternatively, an artificial layer can be introduced as an infinite domain, which would be time-consuming. Even within boundary element methods for which only the boundaries are discretised, one also needs to solve linear systems, which will be a challenge when the number of unknowns is large.

In addition to traditional grid-based numerical methods for partial differential equations (PDEs), novel approaches that incorporate machine learning [3] and neural networks [4], such as methods based on physics-informed neural networks [5, 6], have been proposed. In particular, grid-free Monte Carlo methods have been introduced in the field of geometry processing and rendering [7, 8], which handle the Laplace equation and the Poisson equation. According to Kakutani's principle and the mean value theorem, the walk-on-sphere (WoS) algorithm can determine the expectation of the boundary value of which a random Brownian walk starts from a given point and hits the boundary for the first time to solve a Laplace equation at the given point [9, 10]. Within Monte Carlo methods, no spatial discretisation is required, which shows great advantages when dealing with complex geometries compared to traditional mesh-based methods. Recent efforts have also been devoted to handling other PDEs of complicated physics and inverse problems [11], such as Helmholtz and Yukawa equations [12], Poisson and diffusion equations with varying coefficients for inhomogeneous materials [13], real-time fluid simulation [14], etc. Along with pioneering work [8, 13], many acceleration techniques have been investigated, such as the implementation of the graphics processing unit [15], fast methods to query the closest point within the spatial hierarchy [16], and boundary value caching for the WoS algorithm [17]. In addition, other Monte Carlo grid-free solvers and formulations have been proposed for certain classes of second-order PDEs, including the walk-on-boundary (WoB) solver [18] and the bidirectional formulation for the WoS algorithm [19].

The handling of boundary conditions is critical for Monte Carlo solvers. Similar to the ray tracing of the light transport simulation in rendering [7], random walks are terminated and absorbed at boundaries with Dirichlet boundary conditions. However, it becomes more complicated and a challenge when dealing with other types of boundary conditions. For the Neumann or Robin boundary conditions, it is often required to simulate reflecting random walks (or partially reflecting or absorbed walks) that bounce off the boundary [20–22]. Recently, an efficient and more general walk-on-star (WoSt) has been proposed to handle the Neumann condition [23], which can be viewed as a Monte Carlo estimator for the boundary integral equation of a Laplace problem. For common problems in an open domain, the classical WoS and WoSt algorithms become less efficient because of the enormously useless random walks tending to infinity and the non-recurrent behaviour of Brownian motion. By formulating the appropriate integral equation, the WoB algorithm is reported to handle Dirichlet, Neumann, Robin, and mixed boundary conditions for both the interior and exterior domains [18]. Another efficient approach to exterior problems with pure Dirichlet boundary conditions can be achieved by performing the Kelvin transformation [24, 25]. In this paper, we have demonstrated the spherical inversion of the domain in a more elegant form compared to the previous work [25], which can be generalised to other coordinate transformations. Moreover, we have dealt with a floating potential boundary condition that is neither a pure Dirichlet nor a pure Neumann condition for both the interior and exterior domains. The proposed method has advantages in terms of geometric flexibility and robustness, output sensitivity, and parallelism.

THEORY AND METHODS

Electrostatic BVPs of exterior domains

We start by discussing the common electrostatic boundary value problems (BVPs) with Dirichlet and floating potential boundaries in an open domain Ω_inf, as shown in Figure 1, which reads: 1 $${\nabla }_{{\boldsymbol{x}}^{\prime }}^{2}\varphi \left({\boldsymbol{x}}^{\prime }\right)=-f\left({\boldsymbol{x}}^{\prime }\right)\quad {\boldsymbol{x}}^{\prime }\in {\Omega }_{\text{inf}}$$ 2 $$\varphi \left({\boldsymbol{x}}^{\prime }\right)={\varphi }_{\mathrm{D}}\left({\boldsymbol{x}}^{\prime }\right)\quad {\boldsymbol{x}}^{\prime }\in {\Gamma }_{\mathrm{D}}$$ 3 $$\varphi \left({\boldsymbol{x}}^{\prime }\right)={\varphi }_{m}\hspace*{.5em}(\text{to}\hspace*{.5em}\text{be}\hspace*{.5em}\text{determined})\quad {\boldsymbol{x}}^{\prime }\in {\Omega }_{m}$$ 4 $${\int }_{\partial {\Omega }_{m}}{\hat{\boldsymbol{n}}}_{m}\left({\boldsymbol{x}}^{\prime }\right)\cdot {\nabla }_{{\boldsymbol{x}}^{\prime }}\varphi \left({\boldsymbol{x}}^{\prime }\right)\mathrm{d}{S}^{\prime }=\frac{{Q}_{m}}{{{\epsilon}}_{0}}\hspace*{.5em}m=1,\text{\ldots },M$$ 5 $$\underset{{\boldsymbol{x}}^{\prime }\to \infty }{\mathrm{lim}}\varphi \left({\boldsymbol{x}}^{\prime }\right)=0$$ where $${\nabla }_{{\boldsymbol{x}}^{\prime }}^{2}$$ and ∇_x′, respectively, denote the Laplacian and gradient (spatial derivative) with respect to the field location x′; Γ_D represents the boundaries for the Dirichlet boundary conditions with known potential; f(x′) = ρ(x′)/ϵ with ρ and ϵ = ϵ₀ denoting the charge density and permittivity of the medium (vacuum), respectively; Ω_m ⊂ Ω_inf denotes the free standing metals with unknown floating potential φ_m, whose total charge is denoted by Q_m (refer to as floating potential boundary conditions Equations (3) and (4)); $${\hat{\boldsymbol{n}}}_{m}$$ denotes the normal vector of the boundary ∂Ω_m towards out of the domain Ω_m; dS′ represents the surface area differential element of Ω_m in the primed coordinates denoted by x′; M denotes the number of the free standing metals; and $${\Omega }_{\text{inf}}\in {\mathbb{R}}^{2}$$ or $${\Omega }_{\text{inf}}\in {\mathbb{R}}^{3}$$. In fact, when the potentials φ_m for x ∈ Ω_m are determined, the Poisson equation, Equation (1), can be solved with pure Dirichlet boundary conditions in the open domain $${\Omega }_{\text{inf}}{\backslash}\left(\bigcup \nolimits_{m=1}^{M}{\Omega }_{m}\right)$$. For an open domain problem where the potential φ(x′) approaches zero when x′ → ∞, an artificial boundary is often needed to mimic the infinite domain when methods based on differential equations are adopted, such as the finite element method.

[IMAGE OMITTED. SEE PDF]

Coordinate transformation

Using the inversion transformation $$\boldsymbol{x}:{\Omega }_{\text{inf}}\to {\widetilde{\Omega }}_{\text{inv}}$$, that is, the Kelvin transformation, which reads: 6 $$\boldsymbol{x}\left({\boldsymbol{x}}^{\prime }\right)={\boldsymbol{x}}^{\prime }/\vert {\boldsymbol{x}}^{\prime }{\vert }^{2}$$ the infinite domain Ω_inf can be transformed into bounded region $${\widetilde{\Omega }}_{\text{inv}}$$. Now, the governing equations can be derived in $${\widetilde{\Omega }}_{\text{inv}}$$ (see Part A in Supporting Information S1 for the details), which reads: 7 $${\nabla }_{\boldsymbol{x}}^{2}\psi (\boldsymbol{x})=-\overset{~}{f}(\boldsymbol{x})\quad \boldsymbol{x}\in {\overset{~}{\Omega }}_{\text{inv}}$$ 8 $$\psi (\boldsymbol{x})={\overset{~}{\psi }}_{\mathrm{D}}(\boldsymbol{x})\quad \boldsymbol{x}\in {\overset{~}{\Gamma }}_{\mathrm{D}}$$ 9 $$\psi (\boldsymbol{x})={\overset{~}{\psi }}_{m}(\boldsymbol{x})=\frac{{\varphi }_{m}}{\vert \boldsymbol{x}\vert }(\text{tobedetermined})\quad \boldsymbol{x}\in \partial {\overset{~}{\Omega }}_{m}$$ where ψ(x) = φ(x)/|x| since the potential φ(x′) vanishes as x′ tends to infinity or x approaches zero and we need to remove singularity; $${\widetilde{\Gamma }}_{\mathrm{D}}$$ is the Dirichlet boundary of the transformed region; $$\partial {\widetilde{\Omega }}_{m}$$ denotes the boundary of the transformed domain $${\widetilde{\Omega }}_{m}$$ corresponding to Ω_m; $$\widetilde{f}(\boldsymbol{x})=f(\boldsymbol{x})\vert \boldsymbol{x}{\vert }^{-1}\vert \boldsymbol{g}{\vert }^{-1}$$ with |g| = |x|⁴ in $${\mathbb{R}}^{3}$$ and |g| = |x|² in $${\mathbb{R}}^{2}$$ denoting the determinant of the contra-variant metric tensor g of the transformed Cartesian coordinate {xⁱ} and {x^j} (e.g. xⁱ (i = 1, 2, 3) represent x, y and z coordinates, respectively); $${\widetilde{\psi }}_{\mathrm{D}}(\boldsymbol{x})={\varphi }_{\mathrm{D}}(\boldsymbol{x})/\vert \boldsymbol{x}\vert $$; $${\widetilde{\psi }}_{m}(\boldsymbol{x})={\varphi }_{m}/\vert \boldsymbol{x}\vert $$ for m = 1, …, M. In Equation (7), $${\nabla }_{\boldsymbol{x}}^{2}$$ denotes the Laplacian with respect to the transformed coordinates x. Here, the electric potential in the m-th metallic region φ_m in Equation (9) can be determined, which reads: 10 $${\int }_{\partial {\widetilde{\Omega }}_{m}}\vert \boldsymbol{x}\vert \left[\alpha (\boldsymbol{x})\psi (\boldsymbol{x})+\hat{\boldsymbol{n}}(\boldsymbol{x})\cdot {\nabla }_{\boldsymbol{x}}\psi (\boldsymbol{x})\right]\mathrm{d}S=\frac{{Q}_{m}}{{{\epsilon}}_{0}}$$

Equation (10) is satisfied when the total charge Q_m in the m-th metallic region Ω_m (or $${\widetilde{\Omega }}_{m}$$) are given, where $$\alpha (\boldsymbol{x})=\hat{\boldsymbol{n}}\cdot \boldsymbol{x}/\vert \boldsymbol{x}{\vert }^{2}$$ and ∇_x denotes the gradient with respect to the transformed coordinates x. Here, dS represents the surface area differential element of $${\widetilde{\Omega }}_{m}$$ in the un-primed coordinates denoted by x; and $$\hat{\boldsymbol{n}}$$ denotes the unit normal vector of the boundary $$\partial {\widetilde{\Omega }}_{m}$$. Since ψ(x) is formally determined according to Equation (9), Equation (10) can be rewritten as follows: 11 $${\int }_{\partial {\widetilde{\Omega }}_{m}}\left[{\varphi }_{m}\frac{\hat{\boldsymbol{n}}\cdot \boldsymbol{x}}{\vert \boldsymbol{x}{\vert }^{2}}+\vert \boldsymbol{x}\vert \frac{\partial \psi (\boldsymbol{x})}{\partial n}\right]\mathrm{d}S=\frac{{Q}_{m}}{{{\epsilon}}_{0}}$$ where $$\partial \psi (\boldsymbol{x})/\partial n=\hat{\boldsymbol{n}}\cdot \nabla {\psi }_{\boldsymbol{x}}(\boldsymbol{x})$$ denotes the directional derivative with respect to the direction $$\hat{\boldsymbol{n}}$$. Now, the original Poisson equation Equation (1) on φ(x′) in $${\Omega }_{\text{inf}}{\backslash}\left(\bigcup \nolimits_{m=1}^{M}{\Omega }_{m}\right)$$ has been transformed into another Poisson equation Equations (7)–(9) regarding ψ(x) in $${\widetilde{\Omega }}_{\text{inv}}{\backslash}\left(\bigcup \nolimits_{m=1}^{M}{\widetilde{\Omega }}_{m}\right)$$, with the source transformed from f(x′) into $$\widetilde{f}(\boldsymbol{x})$$ and the boundary conditions transformed from φ_D(x′) and φ_m into $${\widetilde{\psi }}_{\mathrm{D}}(\boldsymbol{x})$$ and ψ(x) = φ_m/|x|, respectively. Furthermore, the total charge constraint Equation (4) has been transformed into Equation (11) accordingly.

Direct integral equations for the potential and the gradient

The solution ψ(x) to the linear Poisson equation can be expressed via a boundary integral involving the associated Green's function G_Ω(x, y) of a domain Ω with a smooth boundary ∂Ω. According to the Green's second identity, the solution to Equation (7) at an evaluation point x in the interior of $$\Omega \subset {\widetilde{\Omega }}_{\text{inv}}$$ can be expressed as follows: 12 $$\psi (\boldsymbol{x})={\int }_{\Omega }{G}_{\Omega }(\boldsymbol{x},\boldsymbol{y})\widetilde{f}(\boldsymbol{y})\mathrm{d}\boldsymbol{y}-{\int }_{\partial \Omega }{\hat{\boldsymbol{n}}}_{\boldsymbol{y}}\cdot {\nabla }_{\boldsymbol{y}}{G}_{\Omega }(\boldsymbol{x},\boldsymbol{y})\psi (\boldsymbol{y})\mathrm{d}\boldsymbol{y}$$

When the domain Ω is deliberately chosen as a ball region centred on point x with radius a, denoted by B_a(x), the analytical expressions of G_Ω(x, y) and ∇_yG_Ω(x, y) can be obtained, as given in Part B in Supporting Information S1. Therefore, Equation (12) reads: 13 $$\psi (\boldsymbol{x})=\frac{1}{\vert \partial {B}_{a}(\boldsymbol{x})\vert }{\int }_{\partial {B}_{a}(\boldsymbol{x})}\psi (\boldsymbol{y})\mathrm{d}\boldsymbol{y}+{\int }_{{B}_{a}(\boldsymbol{x})}{G}_{B}(\boldsymbol{x},\boldsymbol{y})\widetilde{f}(\boldsymbol{y})\mathrm{d}\boldsymbol{y}$$ where |∂B_a(x)| equals 4πa² (or 2πa) denoting the surface area (perimeter) of the sphere (circle) B_a(x); $${G}_{B}(\boldsymbol{x},\boldsymbol{y})={G}_{{B}_{a}(\boldsymbol{x})}(\boldsymbol{x},\boldsymbol{y})$$ is the Green's function on the ball (circle) B_a(x). In addition, the gradient can be expressed as follows: 14 $\begin{align*}\hfill \nabla \psi (\boldsymbol{x})=& \frac{1}{\vert {B}_{a}(\boldsymbol{x})\vert }{\int }_{\partial {B}_{a}(\boldsymbol{x})}\psi (\boldsymbol{y}){\hat{\boldsymbol{n}}}_{\boldsymbol{y}}\mathrm{d}\boldsymbol{y}-\hfill \\ \hfill & {\int }_{{B}_{a}(\boldsymbol{x})}{\nabla }_{\boldsymbol{y}}{G}_{B}(\boldsymbol{x},\boldsymbol{y})\widetilde{f}(\boldsymbol{y})\mathrm{d}\boldsymbol{y}\hfill \end{align*}$ where $${\hat{\boldsymbol{n}}}_{\boldsymbol{y}}=(\boldsymbol{y}-\boldsymbol{x})/a$$ is the outward unit normal at point y ∈ ∂B_a(x); |B_a(x)| = 4πa³/3 (or πa²) denotes the volume (area) of the ball region B_a(x); ∇_yG_B(x, y) is given by Equation (S20) or Equation (S21) in Part B in Supporting Information S1.

Monte Carlo integration

In general, Equations (13) and (14) can be calculated using a Monte Carlo estimator, which approximates an integral by using random samples of the integrand. In particular, for any integrable function $$\phi :\Omega \to \mathbb{R}$$, I = $${\int }_{\Omega }$$ϕ(x)dx can be approximated by the following: 15 $${\widehat{I}}_{N}=\frac{1}{N}\sum\limits _{\ell =1}^{N}\phi \left({\boldsymbol{x}}_{\ell }\right)/p\left({\boldsymbol{x}}_{\ell }\right)\quad {\boldsymbol{x}}_{\ell }\sim p$$ where N is the number of sampling points, and x_ℓ are independent random samples drawn from any probability density p that is nonzero on the kernel (support) of ϕ(x), as shown in Figure 2a. In particular, p = 1/|Ω| for uniform sampling, where |Ω| denotes the volume of the integral domain. Importantly, a well-designed estimator will give the exact value of the integral in expectation. In fact, $${\widehat{I}}_{N}$$ is unbiased if the expectation of $${\widehat{I}}_{N}$$ satisfies $$\mathbb{E}\left({\widehat{I}}_{N}\right)=I$$ even when N = 1. Furthermore, according to the central limit theorem, an unbiased estimator is automatically consistent, with the variance or the ‘error’ reaching zero at a rate of $$\mathcal{O}(1/N)$$ [26], as shown in Figure 2b.

[IMAGE OMITTED. SEE PDF]

Walk-on-sphere algorithm

Within Monte Carlo estimator Equation (15), the direct integral equation Equation (13) for the potential can be written in an estimated form as follows: 16 $$\widehat{\psi }\left({\boldsymbol{x}}_{i}\right)=\frac{1}{\vert \partial {B}_{a}\left({\boldsymbol{x}}_{i}\right)\vert }\frac{\widehat{\psi }\left({\boldsymbol{x}}_{i+1}\right)}{p\left({\boldsymbol{x}}_{i+1}\vert {\boldsymbol{x}}_{i}\right)}+\frac{{G}_{B}\left({\boldsymbol{x}}_{i},{\boldsymbol{y}}_{i}\right)\widetilde{f}\left({\boldsymbol{y}}_{i}\right)}{p\left({\boldsymbol{y}}_{i}\vert {\boldsymbol{x}}_{i}\right)}$$ where $$\widehat{\psi }$$ is the sampled estimator of ψ; x_i+1 ∈ ∂B_a(x_i) denotes a random point on the boundary of the ball region B_a(x_i) centred on x_i; y_i is a random point inside the ball; p(x_i+1|x_i) and p(y_i|x_i) denote the condition probability of the sampling from point x_i to x_i+1 and y_i, respectively. Remarkably, the Equation (13) for ψ(x) can be calculated recursively, which is known as the so-called walk-on-sphere (WoS) algorithm [8, 9, 13]. As illustrated in Figure 3a, starting from x₀ of interest, a sequence of points x₁, x₂, …, x_i forms a random walk. When the next point x_i+1 is located within a small distance ɛ > 0 from the domain boundary, where the potential ψ(x) is known according to the Dirichlet boundary condition, the random walk can be terminated and returns $${\widetilde{\psi }}_{\mathrm{D}}\left(\bar{\boldsymbol{x}}\right)$$ at the boundary point $$\bar{\boldsymbol{x}}$$ closest to x_i. Here, $${\widetilde{\psi }}_{\mathrm{D}}$$ also represents $${\widetilde{\psi }}_{m}$$ for the floating potential boundaries. To take large steps, WoS typically uses the largest sphere centred at x_i and contained entirely in Ω (i.e. $$a=\vert {\boldsymbol{x}}_{i}-\bar{\boldsymbol{x}}\vert $$). Additionally, the source term involving $$\widetilde{f}$$ and G_B can be included in each step of the random walk using a single-point Monte Carlo estimate. Therefore, the recursive estimate of ψ(x) is as given below: 17 $$\widehat{\psi }\left({\boldsymbol{x}}_{i}\right)=\left\{\begin{array}{@{}ll@{}}\frac{1}{\vert \partial {B}_{a}\left({\boldsymbol{x}}_{i}\right)\vert }\frac{{\widetilde{\psi }}_{\mathrm{D}}\left(\bar{\boldsymbol{x}}\right)}{p\left({\boldsymbol{x}}_{i+1}\vert {\boldsymbol{x}}_{i}\right)}\quad \hfill & \vert {\boldsymbol{x}}_{i+1}-\bar{\boldsymbol{x}}\vert < \varepsilon \hfill \\ \frac{1}{\vert \partial {B}_{a}\left({\boldsymbol{x}}_{i}\right)\vert }\frac{\widehat{\psi }\left({\boldsymbol{x}}_{i+1}\right)}{p\left({\boldsymbol{x}}_{i+1}\vert {\boldsymbol{x}}_{i}\right)}+\frac{{G}_{B}\left({\boldsymbol{x}}_{i},{\boldsymbol{y}}_{i}\right)\widetilde{f}\left({\boldsymbol{y}}_{i}\right)}{p\left({\boldsymbol{y}}_{i}\vert {\boldsymbol{x}}_{i}\right)}\quad \hfill & \vert {\boldsymbol{x}}_{i+1}-\bar{\boldsymbol{x}}\vert \ge \varepsilon \hfill \end{array}\right.$$

[IMAGE OMITTED. SEE PDF]

Similarly, as illustrated in Figure 3b, the gradient Equation (14) can be estimated as follows: 18 $$\widehat{\nabla \psi }\left({\boldsymbol{x}}_{0}\right)=\frac{1}{\vert {B}_{a}\left({\boldsymbol{x}}_{0}\right)\vert }\frac{\widehat{\psi }\left({\boldsymbol{x}}_{1}\right)}{p\left({\boldsymbol{x}}_{1}\vert {\boldsymbol{x}}_{0}\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{1}}-\frac{{\nabla }_{{\boldsymbol{y}}_{\boldsymbol{0}}}{G}_{B}\left({\boldsymbol{x}}_{0},{\boldsymbol{y}}_{0}\right)\widetilde{f}\left({\boldsymbol{y}}_{0}\right)}{p\left({\boldsymbol{y}}_{0}\vert {\boldsymbol{x}}_{0}\right)}$$ where $$\widehat{\psi }\left({\boldsymbol{x}}_{1}\right)$$ can be calculated recursively according to Equation (17). Finally, the potential ψ(x₀) at the point x₀ of interest can be obtained using the following equation: 19 $$\psi \left({\boldsymbol{x}}_{0}\right)=\frac{1}{N}\sum\limits _{k=1}^{N}{\widehat{\psi }}_{k}\left({\boldsymbol{x}}_{0}\right)$$ where $${\widehat{\psi }}_{k}\left({\boldsymbol{x}}_{0}\right)$$ denotes the estimate of the k-th random walk and N denotes the number of random walks. Likewise, the corresponding gradient reads: 20 $$\nabla \psi \left({\boldsymbol{x}}_{0}\right)=\frac{1}{N}\sum\limits _{k=1}^{N}{\widehat{\nabla \psi }}_{k}\left({\boldsymbol{x}}_{0}\right)$$

Note that variance reduction techniques can be applied to estimate the gradient since differentiation would amplify the high-frequency variance [27, 28]. Towards this end, we can replace Equation (18) with the following equation: 21 $\begin{align*}\hfill \widehat{\nabla \psi }\left({\boldsymbol{x}}_{0}\right)=& \frac{1}{\vert {B}_{a}\left({\boldsymbol{x}}_{0}\right)\vert }\frac{\widehat{\psi }\left({\boldsymbol{x}}_{0}\right)-{\bar{\psi }}_{i-1}\left({\boldsymbol{x}}_{0}\right)}{p\left({\boldsymbol{x}}_{1}\vert {\boldsymbol{x}}_{0}\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{1}}\hfill \\ \hfill & -\frac{{\nabla }_{{\boldsymbol{y}}_{\boldsymbol{0}}}{G}_{B}\left({\boldsymbol{x}}_{0},{\boldsymbol{y}}_{0}\right)\left[\widetilde{f}\left({\boldsymbol{y}}_{0}\right)-{\overline{\widetilde{f}}}_{i-1}\right]}{p\left({\boldsymbol{y}}_{0}\vert {\boldsymbol{x}}_{0}\right)}\hfill \end{align*}$ where $${\bar{\psi }}_{i-1}\left({\boldsymbol{x}}_{0}\right)$$ denotes the estimate of ψ(x₀) for the first i − 1 walks and $${\overline{\widetilde{f}}}_{i-1}$$ is the average source sampled from the initial ball B_a(x₀) (see Part C in Supporting Information S1 for the details).

Sampling strategies

For Monte Carlo estimators, for example, Equations (17) and (18), the sampling strategies are critical, which determine the probability density functions (PDFs) for the sampling from x_i to x_i+1 and x_i to y_i and, in turn, affect the expressions of the estimators. Although Monte Carlo integration and estimators per se are unbiased, proper sampling methods would be critical to reducing variance and improving efficiency.

For the WoS algorithm Equations (17) and (18), the next point x_i+1 can be sampled uniformly on the surface of the ball region B_a(x_i) centred on x_i, which produces the PDF reading as p(x_i+1|x_i) = 1/|∂B_a(x_i)|. For the source term involving the Green's function G_B(x_i, y_i) which will approach infinity when y_i approaches the centre x_i, ‘importance sampling’ should be considered since the source closer to the centre x_i contributes more to Monte Carlo estimation, as suggested in Ref. [8]. $${\int }_{{B}_{a}(\boldsymbol{x})}{G}_{B}(\boldsymbol{x},\boldsymbol{y})\widetilde{f}(\boldsymbol{y})\mathrm{d}\boldsymbol{y}$$ can be transformed into the corresponding integration $${\int }_{{B}_{a}(\boldsymbol{x})}{G}_{B}^{2\mathrm{D}}(r)\widetilde{f}(\boldsymbol{y})r\mathrm{d}r\mathrm{d}\theta $$ in the local polar coordinates (r, θ), or $${\int }_{{B}_{a}(\boldsymbol{x})}{G}_{B}^{3\mathrm{D}}(r)\widetilde{f}(\boldsymbol{y}){r}^{2}\hspace*{.5em}\sin \hspace*{.5em}\theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\phi $$ in the local spherical coordinates (r, θ, ϕ), with $${G}_{B}^{2\mathrm{D}}$$ and $${G}_{B}^{3\mathrm{D}}$$ being the Green's functions in a circle and sphere, respectively (see Part B in Supporting Information S1); r = |y_i − x_i| denotes the distance from the centre x_i of the circle (ball), and $$r{G}_{B}^{2\mathrm{D}}(r)$$ or $${r}^{2}{G}_{B}^{3\mathrm{D}}(r)$$ would not be singular when r approaches zero. Thus, the ‘normalised Green's function’ can be used as a PDF, which reads p(y_i|x_i) = G_B(x_i, y_i)/|G_B(x_i)|, where $$\vert {G}_{B}\left({\boldsymbol{x}}_{i}\right)\vert ={\int }_{{B}_{a}(\boldsymbol{x})}{G}_{B}(\boldsymbol{x},\boldsymbol{y})\mathrm{d}\boldsymbol{y}$$ reads a²/4 and a²/6 for 2D and 3D, respectively. Since G_B(x_i, y_i) is symmetric, we have $${\boldsymbol{y}}_{i}={\boldsymbol{x}}_{i}+r{\hat{\boldsymbol{y}}}_{i}$$, where $${\hat{\boldsymbol{y}}}_{i}$$ denotes a unit vector that can be sampled uniformly on a unit-sphere or unit-circle. Furthermore, the distance r follows the radial PDF p_r(r), which can be sampled by rejection sampling or Ulirich's polar method [8]. In particular, the radial PDF reads $${p}_{r}^{2\mathrm{D}}(r)=\left(4r/{a}^{2}\right)\ln(a/r)$$ and $${p}_{r}^{3\mathrm{D}}(r)=6r(a-r)/{a}^{3}$$. When handling the source term in Equation (18) to estimate the gradient, we can slightly change the integral as $${\int }_{{B}_{a}(\boldsymbol{x})}\nabla {G}_{B}(\boldsymbol{x},\boldsymbol{y})\widetilde{f}(\boldsymbol{y})\mathrm{d}\boldsymbol{y}={\int }_{{B}_{a}(\boldsymbol{x})}r\nabla {G}_{B}(r)\cdot \widetilde{f}(\boldsymbol{y})/r\mathrm{d}\boldsymbol{y}$$ so that $${r}^{2}\nabla {G}_{B}^{2\mathrm{D}}(r)$$ or $${r}^{3}\nabla {G}_{B}^{3\mathrm{D}}(r)$$ would not be singular when r approaches zero; consequently, similar to the importance sampling for the potential, p(y₀|x₀) can be the normalised gradient of the Green's function of r∇G_B(r), which reads p^2D(y₀|x₀) = 3(1 − r²/a²)/(4πa) and p^3D(y₀|x₀) = (1 − r³/a³)/(3πa). Now, $${\boldsymbol{y}}_{0}={\boldsymbol{x}}_{0}+r{\hat{\boldsymbol{y}}}_{0}$$, where r follows the radial PDF $${p}_{r}^{0}(r)$$, which reads $${p}_{r}^{\text{0,}2\mathrm{D}}(r)=2{\uppi }r{p}^{2\mathrm{D}}\left({\boldsymbol{y}}_{0}\vert {\boldsymbol{x}}_{0}\right)=3r\left(1-{r}^{2}/{a}^{2}\right)/(2a)$$ and $${p}_{r}^{0,3\mathrm{D}}(r)=4{\uppi }{r}^{2}{p}^{3\mathrm{D}}\left({\boldsymbol{y}}_{0}\vert {\boldsymbol{x}}_{0}\right)=4{r}^{2}\left(1-{r}^{3}/{a}^{3}\right)/(3a)$$. In addition to importance sampling of Green's functions, importance sampling of source terms can be used, in particular, for point sources and line sources. Moreover, multiple importance sampling can be used [26].

Monte Carlo methods for floating potentials

We first consider the boundary condition of the floating potential in the untransformed domain, which is governed by Equation (4). Taking the i-th free-standing potential for example, the integral can also be approximated using Monte Carlo estimation according to Equation (15), which reads: 22 $$\frac{1}{{N}_{i}}\sum\limits _{j=1}^{{N}_{i}}\frac{{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}\cdot \nabla \varphi \left({\boldsymbol{x}}_{ij}^{\prime }\right)}{p\left({\boldsymbol{x}}_{ij}^{\prime }\right)}=\frac{{Q}_{i}}{{{\epsilon}}_{0}}\quad i=1,\text{\ldots },M$$ where $${\boldsymbol{x}}_{ij}^{\prime }$$ denotes the j-th randomly sampled point located on the i-th surface of the free-standing metal; N_i denotes the number of points sampled on the i-th surface of the free-standing metal; $${\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}$$ denotes the unit normal vector outward to the domain at the point $${\boldsymbol{x}}_{ij}^{\prime }$$; $$p\left({\boldsymbol{x}}_{ij}^{\prime }\right)$$ denotes the PDF of sampling $${\boldsymbol{x}}_{ij}^{\prime }$$ on ∂Ω_i, which varies for different sampling methods [29]. Here, $$\nabla {\varphi }_{i}\left({\boldsymbol{x}}_{ij}^{\prime }\right)$$ is the gradient at point $${\boldsymbol{x}}_{ij}^{\prime }$$, which can be obtained according to Equation (20), that is, 23 $$\nabla \varphi \left({\boldsymbol{x}}_{ij}^{\prime }\right)=\frac{1}{{N}_{ij}}\sum\limits _{k=1}^{{N}_{ij}}{\widehat{\nabla \varphi }}_{k}\left({\boldsymbol{x}}_{ij}^{\prime }\right)\approx \frac{1}{{N}_{ij}}\sum\limits _{k=1}^{{N}_{ij}}{\widehat{\nabla \varphi }}_{k}\left({\boldsymbol{x}}_{ij}^{\prime \sigma }\right)$$ where N_ij denotes the number of random walks for evaluating $$\nabla {\varphi }_{i}\left({\boldsymbol{x}}_{ij}^{\prime }\right)$$; σ denotes the thickness of the σ-layer on the surface of free-standing metals. Since $${\boldsymbol{x}}_{ij}^{\prime }\in \partial {\Omega }_{i}$$ is located on the boundary, as illustrated in Figure 4a, we use the estimate of $${\widehat{\nabla \varphi }}_{k}\left({\boldsymbol{x}}_{ij}^{\prime \sigma }\right)$$ at an offset point $${\boldsymbol{x}}_{ij}^{\prime \sigma }$$ to approximate $${\widehat{\nabla \varphi }}_{k}\left({\boldsymbol{x}}_{ij}^{\prime }\right)$$, which is given by Equation (18), that is, 24 $\begin{align*}\hfill {\widehat{\nabla \varphi }}_{k}\left({\boldsymbol{x}}_{ij}^{\prime \sigma }\right)=\frac{1}{\vert {B}_{\sigma }\left({\boldsymbol{x}}_{ij}^{\prime \sigma }\right)\vert }\frac{\widehat{\varphi }\left({\boldsymbol{x}}_{ijk}^{\prime \sigma }\right)}{p\left({\boldsymbol{x}}_{ijk}^{\prime \sigma }\vert {\boldsymbol{x}}_{ij}^{\prime \sigma }\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\prime \sigma }}+{\boldsymbol{T}}_{ijk}^{\sigma }\\ ={\eta }_{{\boldsymbol{x}}_{ij}^{\prime \sigma }}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\prime \sigma }}\widehat{\varphi }\left({\boldsymbol{x}}_{ijk}^{\prime \sigma }\right)+{\boldsymbol{T}}_{ijk}^{\sigma }\hfill \end{align*}$ if $${\boldsymbol{x}}_{ijk}^{\prime \sigma }$$ is uniformly sampled on $$\partial {B}_{\sigma }\left({\boldsymbol{x}}_{ij}^{\prime \sigma }\right)$$, where $${\eta }_{{\boldsymbol{x}}_{ij}^{\prime \sigma }}=\vert \partial {B}_{\sigma }\left({\boldsymbol{x}}_{ij}^{\prime \sigma }\right)\vert /\vert {B}_{\sigma }\left({\boldsymbol{x}}_{ij}^{\prime \sigma }\right)\vert ={N}^{\text{dim}}/\sigma $$ denotes the surface to volume ratio of the ball $${B}_{\sigma }\left({\boldsymbol{x}}_{ij}^{\prime \sigma }\right)$$ (N^dim denotes the dimension), $${\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\prime \sigma }}=\left({\boldsymbol{x}}_{ijk}^{\prime \sigma }-{\boldsymbol{x}}_{ij}^{\prime \sigma }\right)/\vert {\boldsymbol{x}}_{ijk}^{\prime \sigma }-{\boldsymbol{x}}_{ij}^{\prime \sigma }\vert $$ and $${\boldsymbol{T}}_{ijk}^{\sigma }$$ reads: 25 $${\boldsymbol{T}}_{ijk}^{\sigma }=-\frac{{\nabla }_{{\boldsymbol{y}}_{ijk}^{\prime \sigma }}{G}_{B}\left({\boldsymbol{x}}_{ij}^{\prime \sigma },{\boldsymbol{y}}_{ijk}^{\prime \sigma }\right)f\left({\boldsymbol{y}}_{ijk}^{\prime \sigma }\right)}{p\left({\boldsymbol{y}}_{ijk}^{\prime \sigma }\vert {\boldsymbol{x}}_{ij}^{\prime \sigma }\right)}$$

[IMAGE OMITTED. SEE PDF]

Here, the subscript k denotes the k-th estimate of the gradient at the point $${\boldsymbol{x}}_{ij}^{\prime \sigma }$$. Furthermore, the potential $$\widehat{\varphi }\left({\boldsymbol{x}}_{ijk}^{\prime \sigma }\right)$$ on the surface $$\partial {B}_{\sigma }\left({\boldsymbol{x}}_{ij}^{\prime \sigma }\right)$$ of the ball centred at the point $${\boldsymbol{x}}_{ij}^{\prime \sigma }$$ can be recursively estimated using the WoS algorithm Equation (17), which reads: 26 $$\widehat{\varphi }\left({\boldsymbol{x}}_{ijk,l}^{\prime \sigma }\right)=\left\{\begin{array}{@{}ll@{}}{\varphi }_{m}\quad \hfill & \vert {\boldsymbol{x}}_{ijk,l+1}^{\prime \sigma }-\bar{{\boldsymbol{x}}^{\prime }}\vert < \varepsilon \hfill \\ \widehat{\varphi }\left({\boldsymbol{x}}_{ijk,l+1}^{\prime \sigma }\right)+{S}_{ijkl}^{\sigma }\quad \hfill & \vert {\boldsymbol{x}}_{ijk,l+1}^{\prime \sigma }-\bar{{\boldsymbol{x}}^{\prime }}\vert \ge \varepsilon \hfill \end{array}\right.$$ if a uniform sampling of $${\boldsymbol{x}}_{ijk,l+1}^{\prime \sigma }$$ on $$\partial {B}_{a}\left({\boldsymbol{x}}_{ijk,l}^{\prime \sigma }\right)$$ is used; and $${S}_{ijkl}^{\sigma }$$ denotes the source contribution, which reads: 27 $${S}_{ijkl}^{\sigma }=\frac{{G}_{B}\left({\boldsymbol{x}}_{ijk,l}^{\prime \sigma }{\boldsymbol{y}}_{ijk,l}^{\prime \sigma }\right)f\left({\boldsymbol{y}}_{ijk,l}^{\prime \sigma }\right)}{p\left({\boldsymbol{y}}_{ijk,l}^{\prime \sigma }\vert {\boldsymbol{x}}_{ijk,l}^{\prime \sigma }\right)}$$

Note that the random walk path of the k-th estimate starting from the point $${\boldsymbol{x}}_{ijk,1}^{\prime \sigma }={\boldsymbol{x}}_{ijk}^{\prime \sigma }$$ would terminate at the ɛ-shell of the m-th Dirichlet boundary with floating potential denoted by φ_m (φ₀ denotes φ_D) after the $${N}_{ijk}^{\text{depth}}$$ path.

Now, combining Equation (22) with Equation (26), we can obtain: 28 $\begin{align*}\hfill \frac{{Q}_{i}}{{{\epsilon}}_{0}}=& \sum\limits _{j=1}^{{N}_{i}}\sum\limits _{k=1}^{{N}_{ij}}\frac{1}{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}^{\prime }\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}\cdot \left[{\eta }_{{\boldsymbol{x}}_{ij}^{\prime \sigma }}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\prime \sigma }}\widehat{\varphi }\left({\boldsymbol{x}}_{ijk}^{\prime \sigma }\right)+{\boldsymbol{T}}_{ijk}^{\sigma }\right]\hfill \\ \hfill & +\sum\limits _{j=1}^{{N}_{i}}\sum\limits _{k=1}^{{N}_{ij}}\frac{1}{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}^{\prime }\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}\cdot {\boldsymbol{T}}_{ijk}^{\sigma }\hfill \\ \hfill & +\sum\limits _{j=1}^{{N}_{i}}\sum\limits _{m=0}^{M}\sum\limits _{k=1}^{{N}_{ij,m}}\frac{1}{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}^{\prime }\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}\cdot {\eta }_{{\boldsymbol{x}}_{ij}^{\prime \sigma }}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\prime \sigma }}{\varphi }_{m}\hfill \\ \hfill & +\sum\limits _{j=1}^{{N}_{i}}\sum\limits _{k=1}^{{N}_{ij}}\sum\limits _{l=1}^{{N}_{ijk}^{\text{depth}}}\frac{1}{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}^{\prime }\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}\cdot {\eta }_{{\boldsymbol{x}}_{ij}^{\prime \sigma }}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\prime \sigma }}{S}_{ijkl}^{\sigma }\hfill \end{align*}$ where N_ij,m denotes the number of random walks that starts from $${\boldsymbol{x}}_{ij}^{\prime \sigma }$$ and terminates on ∂Ω_m, as illustrated in Figure 4b,c. Note that $${\sum }_{m=0}^{M}{N}_{ij,m}={N}_{ij}$$. Finally, one can obtain a linear system for the unknown floating potentials, which reads: 29 $$\left(\begin{array}{@{}ccccc@{}}\hfill {C}_{11}\hfill & \hfill {\cdots}\hfill & \hfill {C}_{1m}\hfill & \hfill {\cdots}\hfill & \hfill {C}_{1M}\hfill \\ \hfill {\vdots}\hfill & \hfill {\ddots}\hfill & \hfill {\vdots}\hfill & \hfill {\ddots}\hfill & \hfill {\vdots}\hfill \\ \hfill {C}_{i1}\hfill & \hfill {\cdots}\hfill & \hfill {C}_{im}\hfill & \hfill {\cdots}\hfill & \hfill {C}_{iM}\hfill \\ \hfill {\vdots}\hfill & \hfill {\ddots}\hfill & \hfill {\vdots}\hfill & \hfill {\ddots}\hfill & \hfill {\vdots}\hfill \\ \hfill {C}_{M1}\hfill & \hfill {\cdots}\hfill & \hfill {C}_{Mm}\hfill & \hfill {\cdots}\hfill & \hfill {C}_{MM}\hfill \end{array}\right)\left(\begin{array}{@{}c@{}}\hfill {\varphi }_{1}\hfill \\ \hfill {\vdots}\hfill \\ \hfill {\varphi }_{i}\hfill \\ \hfill {\vdots}\hfill \\ \hfill {\varphi }_{M}\hfill \end{array}\right)=\left(\begin{array}{@{}c@{}}\hfill {\overline{Q}}_{1}\hfill \\ \hfill {\vdots}\hfill \\ \hfill {\overline{Q}}_{i}\hfill \\ \hfill {\vdots}\hfill \\ \hfill {\overline{Q}}_{M}\hfill \end{array}\right)$$ where C_im and $${\overline{Q}}_{i}$$ can be, respectively, calculated using the following equation: 30 $${C}_{im}=\sum\limits _{j=1}^{{N}_{i}}\sum\limits _{k=1}^{{N}_{ij,m}}\frac{{\eta }_{{\boldsymbol{x}}_{ij}^{\prime \sigma }}}{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}^{\prime }\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}\cdot {\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\prime \sigma }}$$ 31 $\begin{align*}\hfill {\overline{Q}}_{i}=& \frac{{Q}_{i}}{{{\epsilon}}_{0}}-\sum\limits _{j=1}^{{N}_{i}}\sum\limits _{k=1}^{{N}_{ij}}\frac{1}{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}^{\prime }\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}\cdot {\boldsymbol{T}}_{ijk}^{\sigma }+\sum\limits _{j=1}^{{N}_{i}}\sum\limits _{k=1}^{{N}_{ij,0}}\frac{{\eta }_{{\boldsymbol{x}}_{ij}^{\prime \sigma }}{\varphi }_{\mathrm{D}}}{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}^{\prime }\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}\cdot {\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\prime \sigma }}+\sum\limits _{j=1}^{{N}_{i}}\sum\limits _{k=1}^{{N}_{ij}}\sum\limits _{l=1}^{{N}_{ijk}^{\text{depth}}}\frac{{\eta }_{{\boldsymbol{x}}_{ij}^{\prime \sigma }}{S}_{ijkl}^{\sigma }}{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}^{\prime }\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}\cdot {\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\prime \sigma }}\hfill \end{align*}$

Once the floating potentials φ_m are determined, the potential and the gradient at given points can be readily calculated using the WoS algorithm.

It should be noted that the variance of the Monte Carlo approximation of the integration Equation (22) for the floating potential boundary condition depends on the sample strategies of $${\boldsymbol{x}}_{ij}^{\prime }$$ on the boundary ∂Ω_i. In general, we can cast a ray in a random direction from an arbitrary $${\boldsymbol{x}}_{i0}^{\prime }$$ within the domain Ω_i [29], resulting in the PDF as $$p\left({\boldsymbol{x}}_{ij}^{\prime }\right)=p\left({\boldsymbol{x}}_{ij}^{\prime }\vert {\boldsymbol{x}}_{i0}^{\prime }\right)=-{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}\cdot {\nabla }_{{\boldsymbol{x}}_{ij}^{\prime }}{G}_{0}\left({\boldsymbol{x}}_{i0}^{\prime },{\boldsymbol{x}}_{ij}^{\prime }\right)$$, where G₀(x, y) denotes the Green's function in free space. Here, $${\boldsymbol{x}}_{i0}^{\prime }$$ is usually chosen as the centre of the domain. When Ω_i is non-convex, projecting a ray from x_i might have multiple intersections with its boundary ∂Ω_i, and only one is chosen as the next point; accordingly, the PDF should be multiplied by 1/N_int, where N_int denotes the number of intersections [18].

Regarding the boundaries of the floating potential in the transformed domain given by Equation (11), we can calculate the floating potential in a similar way. Now, the elements of the linear system $$\widetilde{\boldsymbol{C}}\cdot \widetilde{\boldsymbol{\psi }}=\widetilde{\boldsymbol{Q}}$$ read as follows: 32 $${\widetilde{C}}_{im}={\delta }_{im}\sum\limits _{j=1}^{{N}_{i}}\frac{{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}}\cdot {\boldsymbol{x}}_{ij}}{{N}_{i}p\left({\boldsymbol{x}}_{ij}\right)\vert {\boldsymbol{x}}_{ij}{\vert }^{2}}+\sum\limits _{j=1}^{{N}_{i}}\sum\limits _{k=1}^{{N}_{ij,m}}\frac{{\eta }_{{\boldsymbol{x}}_{ij}^{\sigma }}\vert {\boldsymbol{x}}_{ij}\vert {\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}}\cdot {\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\sigma }}}{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}\right)}$$ and 33 $\begin{align*}\hfill {\widetilde{Q}}_{i}=& \frac{{Q}_{i}}{{{\epsilon}}_{0}}-\sum\limits _{j=1}^{{N}_{i}}\sum\limits _{k=1}^{{N}_{ij}}\frac{\vert {\boldsymbol{x}}_{ij}\vert }{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}}\cdot {\widetilde{\boldsymbol{T}}}_{ijk}^{\sigma }\hfill \\ \hfill & +\sum\limits _{j=1}^{{N}_{i}}\sum\limits _{k=1}^{{N}_{ij,0}}\frac{{\eta }_{{\boldsymbol{x}}_{ij}^{\sigma }}{\widetilde{\psi }}_{\mathrm{D}}\vert {\boldsymbol{x}}_{ij}\vert }{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}}\cdot {\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\sigma }}\hfill \\ \hfill & +\sum\limits _{j=1}^{{N}_{i}}\sum\limits _{k=1}^{{N}_{ij}}\sum\limits _{l=1}^{{N}_{ijk}^{\text{depth}}}\frac{{\eta }_{{\boldsymbol{x}}_{ij}^{\sigma }}{\widetilde{S}}_{ijkl}^{\sigma }\vert {\boldsymbol{x}}_{ij}\vert }{{N}_{i}{N}_{ij}p\left({\boldsymbol{x}}_{ij}\right)}{\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}}\cdot {\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ijk}^{\sigma }}\hfill \end{align*}$ where 34 $${\widetilde{\boldsymbol{T}}}_{ijk}^{\sigma }=-\frac{{\nabla }_{{\boldsymbol{y}}_{ijk}^{\sigma }}{G}_{B}\left({\boldsymbol{x}}_{ij}^{\sigma },{\boldsymbol{y}}_{ijk}^{\sigma }\right)\widetilde{f}\left({\boldsymbol{y}}_{ijk}^{\sigma }\right)}{p\left({\boldsymbol{y}}_{ijk}^{\sigma }\vert {\boldsymbol{x}}_{ij}^{\sigma }\right)}$$ 35 $${\widetilde{S}}_{ijkl}^{\sigma }=\frac{{G}_{B}\left({\boldsymbol{x}}_{ijk,l}^{\sigma },{\boldsymbol{y}}_{ijk,l}^{\sigma }\right)\widetilde{f}\left({\boldsymbol{y}}_{ijk,l}^{\sigma }\right)}{p\left({\boldsymbol{y}}_{ijk,l}^{\sigma }\vert {\boldsymbol{x}}_{ijk,l}^{\sigma }\right)}$$

Finally, $${\widetilde{\psi }}_{m}$$ can be solved. Therefore, the potential and/or the gradient of interest can be readily calculated within WoS.

Complexity analysis

We now discuss the space and time complexity of the Monte Carlo method and the WoS algorithm. The procedure includes (i) closest point query of N_element geometric elements; (ii) single random walk that converges to the ɛ-layer with a thickness of ɛ to the Dirichlet boundary; and (iii) repetition of N_walk random walks. The time complexity $$\mathcal{T}$$, respectively, reads: (i) $${\mathcal{T}}_{\text{CPQ}}=\mathcal{O}\left(\log {N}_{\text{element}}\right)$$ if the bounding volume hierarchy is adopted [30]; (ii) $${\mathcal{T}}_{\varepsilon }=\mathcal{O}(1/\log \varepsilon )$$ [13]; (iii) $${\mathcal{T}}_{{N}_{\text{walk}}}=\mathcal{O}\left({N}_{\text{walk}}\right)$$. Within the Monte Carlo method, it is not necessary to store the intermediate results for each iteration (step) and walk. Therefore, the space complexities are $${\mathcal{S}}_{\text{geom}}=\mathcal{O}\left({N}_{\text{element}}\right)$$ for the geometry and $${\mathcal{S}}_{\text{field}}=\mathcal{O}\left({N}_{\text{field}}\right)$$ for the field (where N_field is the number of interested field points) due to the progressive and pointwise (parallelism) nature of the Monte Carlo method. In addition, when handling floating potentials within the proposed Monte Carlo method, the time complexity becomes $$\mathcal{T}=\mathcal{O}\left(M\times {N}_{i}\right)$$, where M is the number of electrodes and N_i is the number of samples on the surface of each electrode; and the space complexity remains unchanged.

RESULTS AND DISCUSSION

Impact of the thicknesses of the ɛ-layer and σ-layer

As indicated in Equation (23), the gradient $$\nabla \varphi \left({\boldsymbol{x}}_{ij}^{\prime }\right)$$ (or the normal derivative) at the point $${\boldsymbol{x}}_{ij}^{\prime }$$ on the boundary ∂Ω_i is approximated by the gradient $$\nabla \varphi \left({\boldsymbol{x}}_{ij}^{\prime \sigma }\right)$$ at the offset point $${\boldsymbol{x}}_{ij}^{\prime \sigma }$$, which is located at the distance of σ in the normal direction $${\hat{\boldsymbol{n}}}_{{\boldsymbol{x}}_{ij}^{\prime }}$$. This approximation takes into account the evaluation of potential $$\varphi \left({\boldsymbol{x}}_{ij}^{\prime }\right)$$, which can be influenced by the thickness ɛ of the absorbing ɛ-layer, according to the WoS algorithm, as indicated in Equation (16). Subsequently, we examine the impact of parameters ɛ and σ. Note that when the number of random walks N is sufficiently large, for example, N = 10⁷, the sampling error resulting from the Monte Carlo nature of the method could become insignificant in comparison to the errors originated from the ɛ-layer and the σ-layer. The first example is a simplified model of an overhead transmission line above a flat ground, as shown in Figure 5a, for which the infinite domain of interest (shaded region) can be transformed into a bounded region within the Kevin transform, producing an eccentric region to be considered. Figure 5b plots the average number of steps with respect to the ratio ɛ/r₀ between the thickness ɛ of the ɛ-layer and the curvature (the radius r₀) of the inner wire when random walks end within the ɛ-layer. It is evident that for various curvatures or radii r₀, the thinner the ɛ-layer, the more steps will be required to grab the known value on the Dirichlet boundary. Furthermore, the WoS algorithm reaches the Dirichlet boundary in $$\mathcal{O}\left[\log \left({r}_{0}/\varepsilon \right)\right]$$, on average. Figure 5c plots the relative error $${E}_{\mathrm{r}}$$ of the electric potential and of the electric field E with respect to the number of random walks N. It is evident that the relative error $${E}_{\mathrm{r}}$$ is quite small, which decreases dramatically as the number of random walks N increases. Moreover, for a limited maximum number of random walks, variant reduction methods can also be adopted to improve accuracy. Figure 5d depicts the relative error $${{E}_{\mathrm{r}}}_{\varphi }$$ of the potential with respect to the normalised thickness ɛ/r₀ of the ɛ-layer. The relative error is shown to be very small even when ɛ/r₀ is considerably large, for example, ɛ/r₀ = 10⁻², indicating that the Monte Carlo method is accurate and feasible for electrostatic problems. As the thickness of the ɛ-layer decreases, the relative errors or derivations (variances) of the solution will approach zero. Figure 5e demonstrates the impact of the thicknesses of the σ-layer on the gradient. Due to the inaccuracy in evaluating the electric potential inside the ɛ-layer, the gradient on the Dirichlet boundary may not be accurate when the thickness of σ-layer is comparable to the thickness of the ɛ-layer [18]. When σ is small, the random walk paths starting from the offset point are more likely to adhere very closely to the boundary; therefore, more random walks may be required to obtain accurate results of the electric field (gradient). Moreover, instead of using the gradient approximation at an offset point as indicated by Equation (23), an exact expression of the boundary derivative may be used, for which the WoB method is incorporated [18].

[IMAGE OMITTED. SEE PDF]

Calculation of capacitance matrix

Having demonstrated the feasibility and precision of the Monte Carlo WoS method, we discuss the floating potential condition by addressing a typical problem of a multi-conductor system, as shown in the inset of Figure 6a. In the model, there are two arrays of electrodes for which the non-excited electrodes satisfy the floating potential condition. When one calculates the mutual capacitance matrix within conventional approaches, such as FEMs or boundary element methods, geometry discretisation is necessary for the electrodes and/or the surrounding space; in addition, each electrode is defined as a port or terminal on which excitation is applied; consequently, the capacitance matrix is subsequently obtained through a steady source sweep [31]. Therefore, multiple simulations of solving linear systems are performed, in fact, which may be time-consuming when the number of electrodes is large and the geometry is complicated. In contrast to conventional methods, as indicated in Equation (30), no linear system is solved to render the capacitance matrix. Furthermore, the precision of the capacitance matrix is determined by the number of random walks N_im that start from ∂Ω_i and end in ∂Ω_m and the number of sampling points N_i in ∂Ω_i, indicating that the precision of C_im can be improved by increasing the number of random walks, as illustrated in Figure 6b. As a reference, the mutual capacitance matrix calculated within the conventional FEM simulation is visualised in Figure 6c. We stress that the convergence of random walks is independent of the size and complexity of the geometry but depends only on the thickness of the ɛ-layer, implying that the proposed method is suitable for calculating the equivalent circuit parameters when the geometry is complicated.

[IMAGE OMITTED. SEE PDF]

Lightning striking near a real tree

Finally, we consider the scenario of a lightning strike near a real tree, as shown in Figure 7a. Real lightning processes can be difficult to model accurately; in practice, a fractal-based random model is often adopted to characterise lightning physics [32]. As for real trees, when fine geometries such as leaves and branches are considered within conventional grid-based numerical methods, tremendous meshes and the number of unknowns will be introduced, which in turn makes the linear system unsolvable. Within the proposed method, no mesh is necessary and no linear system is solved, implying that very complicated geometry, such as real trees, can be handled. Furthermore, as indicated in Equations (19) and (20), the evaluation at each point of interest is independent, implying that one does not need to calculate everywhere as usual and only regions or points of interest are necessary; in addition, progressive results can be rendered as the number of samples increases. Figure 7d (Figure 7g) and Figure 7e (Figure 7h) show the distribution of the electric potential (field) when the number of samples is N = 10³ and N = 10⁶, respectively. We note that the potential exhibits a notable variance when N is small; and as the number of samples increases, the results approach the reference (see Figure 7c). We point out that the branches and leaves of the tree may not be modelled as pure metal but rather as dielectrics; it is possible to handle the boundary conditions at dielectric interfaces within Monte Carlo methods, which is currently under development and will be present elsewhere in future work.

[IMAGE OMITTED. SEE PDF]

We note that the computational time for the Monte Carlo method can be longer than that for conventional methods. However, this is not about comparing whether oranges or apples are better. As noted in Ref. [8], discretisation or meshing of complex geometry can be the most expensive procedure of the whole simulation. Therefore, the mesh-free Monte Carlo method may be feasible. Furthermore, variance reduction techniques can be used to obtain accurate results with a small number of samples, as discussed in Part C in Supporting Information S1. Moreover, it is possible to combine the Monte Carlo method with conventional methods, which may become a Swiss-knife tool for electrostatic problems with complex geometries.

CONCLUSION AND PERSPECTIVES

In summary, we have introduced a grid-free Monte Carlo method to handle electrostatic problems with Dirichlet and floating boundary conditions within the WoS algorithm. By applying Green's theorem, the potential at a point of interest can be expressed as an integral or weighted average of the potential on the sphere centred at the given point, which may be evaluated in a recursive manner when the integral is calculated within the Monte Carlo method. According to the Kakutani principle or the Feynman–Kac formula, the solution to a Laplace equation is the expectation of the boundary value for which a random Brownian walk starts from the given point and hits the boundary for the first time, which can be obtained using the WoS algorithm. For the Poisson equation, the source contribution can also be included in each random walk. When floating potential boundary conditions are involved, meshing is still not necessary and only a small linear system is solved. Moreover, when calculating the mutual capacitance matrix of a multi-conductor system, one does not need to solve for swept terminals multiple times, and progressive results towards real values can be obtained as the sampling increases. In addition, only solutions of points or regions of interest need to be evaluated without having to perform a global solution within conventional methods, showing the output sensitivity. Although only Dirichlet boundary conditions along with floating boundary conditions are investigated in this paper, we expect Neumann and Robin boundary conditions to be handled by incorporating other state-of-the-art approaches such as the WoSt [23] and WoB [18] algorithms. Theoretically, provided that an integral equation is derived for the potential or an equivalent source, the Monte Carlo method can be used to handle complex media, such as inhomogeneous media, when the material variations are not large [13], or even with an abrupt change [33, 34]. However, handling a general complex material such as anisotropic media can be a challenge, since the Green's function is complex, which would be an interesting research topic in future work.

ACKNOWLEDGEMENTS

This work was financially supported by the State Grid Corporation of China under grant no. 5500-202155584A-0-5-SF.

CONFLICT OF INTEREST STATEMENT

The authors declare no potential conflict of interests.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available on request from the corresponding author.

References

Seidl, A., Schmidt, T.: Implementation of meshless FEM for engineering applications. Proc. World Acad. Sci. Eng. Technol. 23, 233–237 (2007)

Schaback, R.: Direct discretizations with applications to meshless methods for PDEs. Dolomites Res. Notes Approx. 6(12), 37–50 (2013)

Samaniego, E., et al.: An energy approach to the solution of partial differential equations in computational mechanics via machine learning: concepts, implementation and applications. Comput. Methods Appl. Mech. Eng. 362, 112790 (2020). https://doi.org/10.1016/j.cma.2019.112790

Zang, Y., et al.: Weak adversarial networks for high‐dimensional partial differential equations. J. Comput. Phys. 411, 109409 (2020). https://doi.org/10.1016/j.jcp.2020.109409

Zeng, S., Cai, Y., Zou, Q.: Deep neural networks based temporal‐difference methods for high‐dimensional parabolic partial differential equations. J. Comput. Phys. 468, 111503 (2022). https://doi.org/10.1016/j.jcp.2022.111503

Wang, Y., Zhong, L.: NAS‐PINN: neural architecture search‐guided physics‐informed neural network for solving PDEs. J. Comput. Phys. 496, 112603 (2024) https://doi.org/10.1016/j.jcp.2023.112603

Shirley, P., Wang, C., Zimmerman, K.: Monte Carlo techniques for direct lighting calculations. ACM Trans. Graph. 15(1), 1–36 (1996). https://doi.org/10.1145/226150.226151

Sawhney, R., Crane, K.: Monte Carlo geometry processing: a grid‐free approach to PDE‐based methods on volumetric domains. ACM Trans. Graph. 39(4), 123:1–123:18 (2020). https://doi.org/10.1145/3386569.3392374

Muller, M.E.: Some continuous Monte Carlo methods for the Dirichlet problem. Ann. Math. Stat. 27(3), 569–589 (1956). https://doi.org/10.1214/aoms/1177728169

Sabelfeld, K.K., Simonov, N.A.: Random walks on boundary for solving PDEs, 1st ed. De Gruyter, Berlin, Boston (1994)

Yılmazer, E.F., Vicini, D., Jakob, W.: Solving inverse PDE problems using grid‐free Monte Carlo estimators. arXiv:2208.02114 (2022)

Yang, X., Rasila, A., Sottinen, T.: Walk on spheres algorithm for Helmholtz and Yukawa equations via Duffin correspondence. Methodol. Comput. Appl. Probab. 19(2), 589–602 (2017). https://doi.org/10.1007/s11009‐016‐9504‐9

Sawhney, R., et al.: Grid‐free Monte Carlo for PDEs with spatially varying coefficients. ACM Trans. Graph. 41(4), 53:1–53:17 (2022). https://doi.org/10.1145/3528223.3530134

Rioux‐Lavoie, D., et al.: A Monte Carlo method for fluid simulation. ACM Trans. Graph. 41(6), 240:1–240:16 (2022). https://doi.org/10.1145/3550454.3555450

Mossberg, L.: GPU‐accelerated Monte Carlo geometry processing for gradient‐domain methods. Dissertation, Linköping University (2021)

Pharr, M., Jakob, W., Humphreys, G.: Physically based rendering: from theory to implementation. MIT Press, Cambridge (2023)

Miller, B., et al.: Boundary value caching for walk on spheres. ACM Trans. Graph. 42(4), 82:1–82:11 (2023). https://doi.org/10.1145/3592400

Sugimoto, R., et al.: A practical walk‐on‐boundary method for boundary value problems. ACM Trans. Graph. 42(4), 81:1–81:16 (2023). https://doi.org/10.1145/3592109

Qi, Y., et al.: A bidirectional formulation for walk on spheres. Comput. Graph. Forum 41(4), 51–62 (2022). https://doi.org/10.1111/cgf.14586

Makarov, N.R.: Monte Carlo methods for solving boundary value problems of the second and third kinds. Russ. J. Numer. Anal. Math. Model. 13(2), 117–132 (1998). https://doi.org/10.1515/rnam.1998.13.2.117

Sipin, A.S.: Random walk on balls for a boundary value problems with mixed boundary conditions. AIP Conf. Proc. 2849(1), 100001 (2023). https://doi.org/10.1063/5.0162520

Simonov, N.A.: Walk‐on‐spheres algorithm for solving third boundary value problem. Appl. Math. Lett. 64, 156–161 (2017). https://doi.org/10.1016/j.aml.2016.09.008

Sawhney, R., et al.: Walk on stars: a grid‐free Monte Carlo method for PDEs with Neumann boundary conditions. ACM Trans. Graph. 42(4), 80:1–80:20 (2023). https://doi.org/10.1145/3592398

Bogdan, K., Żak, T.: On Kelvin transformation. J. Theor. Probab. 19(1), 89–120 (2006). https://doi.org/10.1007/s10959‐006‐0003‐8

Nabizadeh, M.S., Ramamoorthi, R., Chern, A.: Kelvin transformations for simulations on infinite domains. ACM Trans. Graph. 40(4), 97:1–97:15 (2021). https://doi.org/10.1145/3476576.3476658

Veach, E., Guibas, L.J.: Optimally combining sampling techniques for Monte Carlo rendering. In: Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Techniques, pp. 419–428. Association for Computing Machinery (1995)

Alanko, S., Avellaneda, M.: Reducing variance in the numerical solution of BSDEs. Compt. Rendus Math. 351(3), 135–138 (2013). https://doi.org/10.1016/j.crma.2013.02.010

Titsias, M., Shi, J.: Double control variates for gradient estimation in discrete latent variable models. In: Proceedings of the 25th International Conference on Artificial Intelligence and Statistics, vol. 151, pp. 6134–6151. PMLR (2022)

Barill, G., et al.: Fast winding numbers for soups and clouds. ACM Trans. Graph. 37(4), 43:1–43:12 (2018). https://doi.org/10.1145/3197517.3201337

Lauterbach, C., et al.: Fast BVH construction on GPUs. In: Computer Graphics Forum, vol. 28(2), pp. 375–384. Wiley Online Library (2009). https://doi.org/10.1111/j.1467‐8659.2009.01377.x

Zhu, Q., et al.: Four‐plate capacitive power transfer system with different grounding connections. Int. J. Electr. Power Energy Syst. 115, 105494 (2020). https://doi.org/10.1016/j.ijepes.2019.105494

Dellera, L., Garbagnati, E.: Lightning stroke simulation by means of the leader progression model. I. Description of the model and evaluation of exposure of free‐standing structures. IEEE Trans. Power Deliv. 5(4), 2009–2022 (1990). https://doi.org/10.1109/61.103696

Ermakov, S.M., Sipin, A.S.: The “walk in hemispheres” process and its applications to solving boundary value problems. Vestn. St. Petersbg. Univ. Math. 42(3), 155–163 (2009). https://doi.org/10.3103/s1063454109030029

Simonov, N.A.: Random walks for solving boundary‐value problems with flux conditions. In: International Conference on Numerical Methods and Applications, pp. 181–188. Springer (2006)