Electrostatic potential calculation by application of walk-on-spheres method

The paper is devoted to mixed boundary-value problem solving for Laplace equation with the use of walk-on-spheres algorithm. The problem under study is reduced to finding a solution of integral equation with the kernel nonzero only at some sphere in the domain considered. Ulam-Neumann scheme is applied for integral equation solving; the appropriate Markov chain is introduced. The required solution value at a certain point of the domain is approximated by the expected value of special statistics defined on Markov paths. The algorithm presented guarantees the average Markov trajectory length to be finite and allows one to take into account boundary conditions on required solution derivative and to avoid Markov paths ending in the neighborhood of the boundaries where solution values are not given. The method is applied for calculation of electric potential in the injector of linear accelerator. The purpose of the work is to verify the applicability and effectiveness of walk-on-spheres method for mixed boundary-value problem solving with complicated boundary form and thus to demonstrate the suitability of Monte Carlo methods for electromagnetic fields simulation in beam forming systems.
The numerical experiments performed confirm the simplicity and convenience of this method application for the problem considered.


Introduction
The problems of charged beam dynamics modeling imply electromagnetic field calculation in beam propagation domain. For example, for magnetic field determination in quadrupole focusing channel it is necessary to solve boundary-value problem for Laplace equation for static magnetic potential. Space charge field potential satisfies Poisson equation, and general solution finding involves Laplace equation solving for the appropriate boundary conditions. The questions of electromagnetic fields modeling and simulation in different focusing and accelerating systems, and corresponding beam dynamics investigation and optimization, are widely discussed. Electromagnetic fields analytical representation is often used in scientific research, in particular, [Ovsyannikov, 2012;Ovsyannikov and Altsybeyev, 2013;Rubtsova, 2014a,b;Rubtsova, 2016a;Rubtsova and Ovsyannikov, 2018;Vinogradova, Starikova and Varayun ′ , 2017;Altsybeyev et al., 2018]. Analytical field description makes possible beam dynamics optimization theory developing on the basis of analytical methods, including trajectory ensemble control methods [Ovsyannikov, 1997;Ovsyannikov et al., 2006;Ovsyannikov, 2012;Ovsyannikov and Altsybeyev, 2013;Rubtsova and Suddenko, 2012;Rubtsova, 2016b;Rubtsova and Ovsyannikov, 2018;Balabanov, Mizintseva, Ovsyannikov, 2018;Nikolskii and Belyaevskikh, 2017;Nikolskii and Belyaevskikh, 2018]. The problems of electromagnetic fields calculation in electrophysical systems may be reduced to solving boundary-value problems for partial derivative equations. A significant number of authors use finitedifference methods for solution calculation [Svistunov and Kozynchenko, 2004;Kozynchenko and Svistunov, 2006;Kozynchenko and Ovsyannikov, 2009;Altsybeyev and Ponomarev, 2015;Ma-yu-shan and Altsy-beyev, 2016]. Sometimes it is possible to reduce partial derivative equation to ordinary differential equation and to use corresponding methods [Golovkina, 2017]. V.A. Kozynchenko suggested the approximate analytical expressions for internal beam field [Kozynchenko and Kozynchenko, 2014a,b]. The approaches mentioned are applied for computer simulation and optimization of the processes in accelerating and focusing systems. The appropriate algorithms, software and results are discussed, in particular, in the articles [Bondarev et al., 2001;Kozynchenko, 2014;Altsybeyev and Ponomarev, 2015]. This paper deals with electric field calculation in the injection system of linear accelerator with due account of accelerating structure geometry. Electrostatic potential is described to be a solution of mixed boundaryvalue problem for Laplace equation.
In some cases (complicated boundary form, large dimension and so on) deterministic methods of such problems solving are cumbersome, and it is proposed to use Monte Carlo methods. Monte Carlo methods find an increasing application for solving applied problems, and the extensive literature is devoted to this subject, in particular, [Sobol, 1973;Sabelfeld, 1989;Mikhailov and Vojtishek, 2006;Ermakov and Sipin, 2014;Vladimirova, 2014;Vladimirova, Ovsyannikov and Rubtsova, 2015;Vladimirova and Rubtsova, 2016].
In this paper electrostatic potential in injection system is calculated with the use of walk-on-spheres algorithm. Boundary-value problem for Laplace equation is reduced to inhomogeneous integral equation solving. Markov chain is constructed on the basis of this equation; kernel function is considered to be transition probability density. The special random variable is defined on Markov trajectories for unknown function estimation at a certain point of the domain. The essence of the method is to simulate this statistics and to approximate the solution value by the mean value of the estimator. This algorithm and its analogs are simple for realization and successfully used by a number of authors [Ermakov and Mikhailov, 1976;Ermakov, 1975;Sabelfeld, 1989;Mikhailov and Vojtishek, 2006;Ermakov and Sipin, 2014]. In particular, this method is applied to Dirichlet problems solving for the Helmholtz equation [Ermakov and Mikhailov, 1976;Mikhailov and Vojtishek, 2006]. In the monograph [Ermakov and Sipin, 2014] Monte-Carlo methods are used for calculation of the solutions of the first boundary-value problem for Laplace and Poisson equations. Specifically, random walks on spheres, hemispheres, ellipsoids are used. Walk-on-spheres algorithm is not very sensitive to the dimension of the problem and region boundary form. The method provides the possibility of solution calculating at certain point of search area using no informa-tion about the solution at other points. In the process of solution calculation its derivatives and error estimation computing may be performed with negligible runtime increasing. As most of Monte Carlo methods, walk-on-spheres method is perfectly suitable for parallel computing: the statistical trials may be distributed between the processes. The main drawback of the algorithm is its slow convergence (∼N −1/2 ), where N is the number of trials. To compensate this limitation, walk-on-spheres method may be successfully combined with other methods of computational mathematics [Ermakov, 1975]. The features of our research are as follows. Firstly, we apply walk-on-spheres algorithm in specific area, namely, electromagnetic fields simulation in accelerating and focusing systems. Secondly, we solve mixed boundary-value problem for Laplace equation and therefore propose algorithm modification allowing us to take into account boundary conditions on solution derivative and avoid Markov paths ending in the neighborhood of the boundaries where solution values are not given. Thirdly, we consider the region with complicated border shape to make sure that it is not a serious difficulty for Monte Carlo method unlike deterministic methods. We set a goal to verify the applicability and effectiveness of walk-on-spheres method for mixed boundaryvalue problem solving with complicated boundary form and thus to confirm the suitability of Monte Carlo methods for electromagnetic fields simulation in beam forming systems.
2 Physical Statement of the Problem Linear accelerator injector presents cylindrical tube of radius R u and length L u . The tube with symmetry axis Oz is equipped with electrode system (Fig. 1). Electric field is created by the electrodes A, B, C, D with potentials u A , u B , u C , u D correspondingly. Potential values at injector ends are u 0 and u L . The potential between the electrodes is supposed to be distributed linearly.
The problem is to determine electric potential in beam propagation region.
3 Mathematical Statement of the Problem Due to device axial symmetry the problem is reduced to potential determination in the domain G with the boundaries Γ 0 , Γ b , Γ e , Γ L (Fig. 2).
The required potential u(r, z) (where r, z are cylindrical coordinates) is the solution of mixed boundaryvalue problem: The total boundary of the domain G is Γ = Γ 0 ∪ Γ 1 . We use the following approach to solving the problem (1)-(3). This problem is reduced to finding a solution of special integral equation with the use of Markov chains. Now let us present the integral equation of appropriate type and discuss the solving method.

Inhomogeneous Integral Equation
Consider the integral equation Let us introduce the notation for integral operators: The solution of the equation (4) may be presented by Neumann series [Ermakov and Mikhailov, 1976] If ∥K∥ < 1 , the series (5) converges in norm and the solution of the equation (4) exists.
5 Neumann-Ulam Scheme for Integral Equation Solving by Monte Carlo Method Let us introduce Markov chain with initial probability distribution density p(x) and transition probability density p(x, y), corresponding the transition from the state x to the state y. The densities mentioned satisfy the conditions: where g(x) is the probability of path ending.

Markov Paths Modeling Algorithm
The simple algorithm is as follows [Ermakov, 1975]. Let us introduce the trajectories ω τ of random length τ in the domain G: following the rules: 1. initial state x 0 ∈ G is determined using initial density p(x) simulation; 2. the path may end in the state x µ , µ = 0, 1, 2 . . ., with a probability g(x µ ) (in this case τ = µ), or continue with probability 1 − g(x µ ); 3. if the path continues, the next state x µ+1 is simulated in accordance with transition density corresponds to Markov chain introduced. The shortened form of this expression is p 0 p 01 . . . p τ −1 τ g τ .

Unbiased Estimator of Scalar Product
Let us consider the scalar product The estimator of scalar product (h, u) is given by the statistics: Here It is proved [Ermakov, 1975] that the following conditions (A) the convergence of the series are necessary and sufficient to fulfill the equality The result (8) implies that random variable (7) is unbiased estimator of the functional (h, u).
On the basis of law of large numbers we obtain approximate equality where N is the number of independent trajectories of Markov chain. Consequently, in view of (8) The result (9) provides the way of required solution calculation.

Integral Equation for Dirichlet Problem
Consider Dirichlet problem for the Laplace equation: If the function u(x) satisfying the equation (10) is defined on sphere S R with radius R and center x 0 , mean value theorem allows us to affirm [Ermakov, 1975]: where |S R | is the sphere area value.
Using (11) and (12) one can obtain integral equation (4) for unknown function u(x) with the kernel where The point x ∈ G is the center of sphere S R . The sphere is supposed to be located inside the domain G or can contain boundary points (for example, the sphere may be tangent to the boundary Γ). The kernel (13) is nonzero only at sphere S R and takes there constant value 1/|S R |. The term (14) is nonzero only at boundary Γ of the domain G. So the Dirichlet problem (10), (11) is reduced to special integral equation (4).
7 Walk-on-spheres Algorithm for Dirichlet Problem Let us apply Neumann-Ulam scheme to integral equation (4). The aim is to find a solution of Dirichlet problem (10)-(11) at a certain point, i.e. to obtain u(x 0 ). The results (8), (9) are extended to special case [Sobol, 1973]: Due to the properties of kernel (13) one can consider it to be transition density for Markov chain. Let us assume that p(x, y) = k(x, y), p(x) = h(x). The estimator (7) takes the form: The algorithm is as follows. Consider the state x µ . There are two cases: Case 1. x µ ∈ Γ. In view of (13) k(x µ , y) = 0 ∀y; the equation (6) gives g(x µ ) = 1; x µ is the final state of the path; τ = µ; in view of (14) f (x τ ) = φ(x τ ); the statistics (16) takes the value ξ(ω τ ) = f τ = φ(x τ ). Case 2. x µ is the inner point of the domain G. The state x µ is not the final state of the path: g(x) = 0 with due account of (13), (6). To determine the value of statistics (16) one should pass to the next state. Note that f (x µ ) = 0 in view of (14). We can obtain nonzero value of the estimator (16) only by reaching the boundary. The transition to the next state is to be performed as follows. The point x µ is considered to be the center of some sphere belonging G or containing boundary points. In view of (13) the transition density p(x µ , x µ+1 ) is nonzero only on sphere S R ; the density is constant on this sphere. Consequently, the state x µ+1 is the realization of random variable uniformly distributed on the sphere S R . Now we can repeat the consideration of two cases for x µ+1 and so on. The process continues until the boundary is reached. Following this algorithm, we obtain N trajectories starting at point x 0 , N values ξ(ω τj ), j = 1, N of estimator (16) and receive in view of (8), (9), (14),(15). It should be noted that radius R depends on x µ . To diminish running time, one can choose R as large as possible. However, the following problem arises. The event of reaching of boundary Γ by the finite Markov path has the zero probability since the norm of operator with kernel (13) is equal 1 [Ermakov, 1975].

Biased Estimator of Integral Equation Solution
To avoid this difficulty, let us assume the kernel function in equation (4) to be zero in the layer Γ(ε) of small thickness ε covering the boundary Γ from inside the domain G (Fig. 3): The known function in equation (4) is now introduced as follows: Let us denote the unknown function by u ε (x). Now Markov paths end when they reach the layer Γ(ε). In other words, the domain boundary is supposed to be "thick" with thickness ε. The event of reaching the layer Γ(ε) by the finite Markov path has positive probability. The unknown function values in the layer Γ(ε) may be taken from the nearest points of the boundary Γ: Here x * τ is boundary point nearest to x τ . At Fig. 3 the point x 3 ∈ Γ(ε) is the final state of Markov path; x * 3 ∈ Γ is the nearest boundary point. In this case the norm of integral operator K with kernel (17) is less than 1, the corresponding series (5) converges in norm and the solution of the equation (4) exists although the estimator of solution value turns out to be biased.
Instead of unbiased estimator ξ(ω τ ) of the solution value u(x 0 ) of integral equation (4), (13), (14) we obtain the ε-biased estimator ξ ε (ω τ ) taking the value φ(x * τ ) on the finite Markov trajectory. As pointed out in the monograph [Ermakov and Mikhailov, 1976], where C is a constant . Walk-on-spheres algorithm is now as follows. Consider x 0 ∈ G\Γ(ε). The transition is performed to the point x 1 presenting the realization of random variable uniformly distributed on the sphere S R with the center x 0 . Radius R is the maximal radius providing S R ⊂ G ∪ Γ. If x 1 ∈ Γ(ε) then τ = 1, the estimator value is ξ ε (ω τ ) = φ(x * 1 ). If x 1 / ∈ Γ(ε) we build the new sphere with center x 1 of maximal allowable radius, and so on (Fig. 3) until the layer Γ(ε) is reached. We repeat the process to obtain N Markov trajectories. The average path length is finite [Ermakov and Mikhailov, 1976].

Probabilistic Error of the Method
In most cases it is easy to estimate the value of dispersion D(ξ ε ) during the calculation of expected value M (ξ ε ): As known [Sobol, 1973], the probabilistic error r N is .
10 Mixed Boundary-Value Problem for Laplace Equation Consider the Laplace equation with the following conditions on the boundary Γ of the domain G: where Γ 1 ∪ Γ 0 = Γ. In view of condition (19) the following formula is true near the boundary Γ 0 : u(x + Hn 0 (x)) = u(x − Hn 0 (x)) + o(2H). (20) Here n 0 (x) is the unit outward-pointing normal to the boundary Γ 0 in the point x ∈ Γ 0 , H is absolute value of increment. The values of required solution u(x) at the boundary Γ 0 are unknown. So the algorithm of ε-biased estimation of the solution value must be modified. The idea of modification is as follows. Let Γ 0 (2ε 0 ) be 2ε 0 -layer surrounding the boundary Γ 0 on both sides (from inside and outside the domain G). If some point x ν of Markov path falls in the layer Γ 0 (2ε 0 ) one should direct the trajectory from the border into the domain G using the condition (19). The way to achieve this goal is to increase the sphere radius to obtain the point x ν outside the domain G ∪ Γ 0 (2ε 0 ). In this case we can reflect the point x ν symmetrically to the border Γ 0 into the domain G\Γ 0 (2ε 0 ) using (20). The resulting point x ν allows us to continue (or to finish) the trajectory. Markov path must end in ε-layer Γ 1 (ε) covering the boundary Γ 1 from inside the domain G. Let us present the version of implementation of this idea. Consider the state x µ of Markov chain; suppose that x µ ∈ G\{Γ 1 (ε) ∪ Γ 0 (2ε 0 )}. Let ε 0 = ε/4 to provide 2H ≤ ε in equality (20). The algorithm is as follows. 1. Let d 0 be the minimal distance from the point x µ to the points of the boundary Γ 0 : d 0 = ϱ(x µ , Γ 0 ). Analogically, let d 1 = ϱ(x µ , Γ 1 ). We introduce the sphere S R (x µ ) with the center x µ and radius R = min{d 0 + 2ε 0 , d 1 }. 2. Next point x µ+1 is obtained as a realization of random variable uniformly distributed on the sphere S R (x µ ). 3. Several cases arise.

Numerical Results
Let us present the results of the described algorithm application to solving problem (1)-(3) for electrostatic potential in the injection system. (See the schemes on Fig. 1, 2).
In the case considered Γ 1 = Γ b ∪ Γ e ∪ Γ L . The device parameters are as follows.
The length of tube L u = 40 cm, its radius is R u = 1 cm. The potentials of electrodes are: u A = −30 kV, u B = −80 kV, u C = −30 kV, u D = −80 kV.
The potential values at device input and exit are u 0 = −40 kV and u L = 0 correspondingly. The numerical solution of the problem (1)-(3) is obtained with the use of walk-on-spheres method. Modified algorithm of ε-biased estimation of solution value (presented in section 10) is applied. The initial points of Markov trajectories were chosen at the nodes of the uniform grid introduced in the domain G. The grid steps along the axes z and r are h z = 1 cm, h r = 0.1 cm while ε = h r /10. The number of Markov paths N = 3000. Fig. 5 presents the graph of potential distribution in the region: r ∈ [0.1; 1], z ∈ [0; 40] (the values are given in centimeters).

Conclusion
The paper deals with the problem of electrostatic potential determining in injection system with due account of geometry of accelerating-focusing structure. The problem under study is formulated as a mixed Figure 5. Potential distribution obtained with the use of walk-onspheres method boundary-value problem for Laplace equation. The potential numerical calculation is performed with the use of Monte Carlo method, specifically, walk-on-spheres algorithm. Our goal is to try the method for determining electromagnetic fields in beam forming systems, and this investigation is successfully performed using injecting system as an example. Modified algorithm of ε-biased estimation of the solution value of mixed boundary-value problem is proposed. The appropriate software is developed and numerical experiments are performed confirming the simplicity and convenience of walk-on-spheres method application to the problem considered. The complicated border shape does not present the serious difficulty for Monte Carlo method. We see the prospect of this investigation in application of parallel computing to reduce runtime.