Yukawa Potential, Panharmonic Measure and Brownian Motion

In [25] a Walk On Spheres (WOS) algorithm for Monte Carlo simulation of the solutions of the Yukawa and the Helmholtz PDE's was developed by using the so-called Duffin correspondence. In this paper we investigate the foundations behind the algorithm for the case of the Yukawa PDE. We study the panharmonic measure that is a generalization of the harmonic measure for the Yukawa PDE. We show that there are natural stochastic definitions for the panharmonic measure in terms of the Brownian motion and that the harmonic and the panharmonic measures are all mutually equivalent. Furthermore, we calculate their Radon--Nikodym derivatives explicitly for some balls, which is a key result behind the WOS algorithm.


Introduction and Preliminaries
The harmonic measure is a fundamental tool in geometric function theory, and it has interesting applications in the study of bounded analytic functions, quasiconformal mappings and potential theory. For example, the harmonic measure has proven very useful in the study of quasidisks and related topics (see, e.g., [1][2][3]). Results involving the harmonic measure have been given by numerous authors since the 1930s (see [4] and references therein). In this paper, we consider the panharmonic measure, which is a natural counterpart of the classical harmonic measure, whereby the harmonic functions related are replaced with the smooth solutions to the Yukawa equation: Equation (1) first arose from the work of the Japanese physicist Hideki Yukawa in particle physics. Here, u : D → R is a two-times differentiable function and D ⊂ R n , n ≥ 2 is a domain. The Yukawa equation was first studied in order to describe the nuclear potential of a point charge. This model led to the concept of the Yukawa potential (also called a screened Coulomb potential), which satisfies an equation of the type given by Equation (1). The Yukawa equation also arises from certain problems related to optics (see [5]). Clearly, when µ = 0, we have the Laplace equation, and, indeed, the results given in this paper reduce to the classical results.
Using the terminology of Duffin [6,7], we call a function u : D → R panharmonic, or µ -panharmonic, in a domain D if its second derivatives are continuous and if it satisfies the Yukawa equation (Equation (1)) for all x ∈ D . The function u is called panharmonic at x 0 ∈ D if there is a neighborhood of x 0 where u is panharmonic.
In Definition 1 of the panharmonic measure below, and in all that follows, we always assume that n ≥ 2, although some results are also true in the dimension n = 1 . For Definition 1, we need the notions of smallness and the regularity of a domain. These are best given by using the stochastic characterization via the Brownian motion. We refer to any of the classical textbooks [8][9][10] for further details.
We recall that the n -dimensional Brownian motion W = (W(t); t ≥ 0) starting from the point x ∈ R n is the time-homogeneous Markov process with the Markov semigroup given by that is, 1 2 ∆ is the generator of the Markov semigroup of the Brownian motion. A domain D ⊂ R n is regular if the Brownian motion does not dwell on its boundary; more where P x is the probability measure under which P x [W(0) = x] = 1 , and is the first hitting time of the Brownian motion in the set D c . We call a regular domain D (Wiener) small if a Brownian motion starting inside D eventually will leave the domain; that is, D is small if For example, all bounded domains are small. All half-spaces are also small.
The panharmonic, or µ -panharmonic, measure is a generalization of the harmonic measure: Definition 1. Let D ⊂ R n be a small regular domain, and let µ 2 ≥ 0 . The µ -panharmonic measure on a boundary ∂D with a pole at x ∈ D is the measure H x µ (D; ·) such that any bounded µ -panharmonic function u onD admits the representation The existence and uniqueness of a panharmonic measure is established by Theorem 1 and Corollary 2 later. Indeed, by Theorem 1, all bounded solutions to the Dirichlet problem ∆u − µ 2 u = 0 on a small regular domain with continuous and bounded boundary data are given by the panharmonic measure as in Equation (2). By Corollary 2, if µ 2 > 0 , then the assumption that the domain is small can be removed; that is, all bounded solutions on a regular domain are of the form given by Equation (2) if the boundary data is bounded and continuous. Of course, it is well known that there are unbounded solutions to the Laplace equation that do not admit the harmonic measure representation. The same is true for the Yukawa equation. We refer to Evans [11] for more details on the solutions of the Laplace equation.
We note that if we replace the "killing parameter" µ 2 in the Yukawa equation (Equation (1)) with a "creation parameter" λ < 0, we obtain another important partial differential equation (PDE), the Helmholtz equation. In principle, the stochastic approaches taken in this paper can be applied to the solutions of the Helmholtz equation if the domain D is small enough compared to the parameter λ . For details, we refer to Chung and Zhao [8]. If we replace µ 2 by a (positive) function, we obtain the Schrödinger equation. Again, the stochastic approaches taken in this paper can be applied, in principle, to the Schrödinger equation, but the results may not be mathematically very tractable. Again, we refer to Chung and Zhao [8] for details. The rest of the paper is organized as follows: In Section 2, we show three different connections between the panharmonic measures and the Brownian motion. The first two (Theorem 1 and Corollary 1) are essentially well known. The third (Corollary 2) is new. In Section 3, we show that the panharmonic measures and the harmonic measures are all mutually equivalent (Theorem 3) and provide some corollaries; namely, we provide a domination principle for the Dirichlet problem related to the Yukawa equation (Corollary 3) and analogs of theorems of Riesz-Riesz, Makarov and Dahlberg for the panharmonic measures (Corollary 4). In Section 4, we consider the panharmonic measures on balls and prove an analogue of the Gauss mean value theorem, or the average property, for the panharmonic functions (Theorem 4), and as a corollary, we obtain the Liouville theorem for panharmonic functions (Corollary 5). Finally, in Section 5, we discuss extensions to the Schrödinger and the Helmholtz PDEs and the walk-on-spheres (WOS) simulation of PDEs.

Yukawa Equation and Brownian Motion
We first recall the celebrated connection between the harmonic measure and the Brownian motion first noticed by Kakutani [12] in the 1940s; the harmonic measure is the hitting measure: Theorem 1 below is a variant of the Kakutani connection (Equation (3)). A key component of the variant is the following disintegration of the harmonic measure at the time the associated Brownian motion hits the boundary ∂D : Lemma 1. Let D ⊂ R n be a regular domain, and let x ∈ D . Then is the harmonic kernel.
Proof. First, we show the existence of the regular conditional distribution: For this, we note that the random vector (W(τ D ), τ D ) can be considered as a function from a space of continuous functions that are the Brownian trajectories equipped with the metric For Brownian trajectories, the metric d is almost surely finite because of the independent increments of the Brownian motion and the Borel-Cantelli lemma. Additionally, with the metric d , the space of Brownian paths is a Polish space. Now, by Theorem A1.2 of [13], Polish spaces are Borel spaces. Consequently, for any fixed x ∈ D , by Theorems 6.3 and 6.4 of [13], the probability kernel (Equation (5)) exists and is measurable with respect to t . Consequently, the harmonic kernel is measurable with respect to t .
Second, we show that the distribution of the hitting time τ D is absolutely continuous with respect to the Lebesgue measure. Let ε > 0 be small enough so that . Now, the distribution of τ B is absolutely continuous; see, for example, the section of Bessel processes in Borodin and Salminen [14]. Additionally, because of the rotation symmetry of the Brownian motion, τ B and τ D − τ B are independent. Hence, by disintegration and independence, we obtain that Thus, the distribution of τ D is absolutely continuous when the distribution of τ B is absolutely continuous.
Third, we show that the formula given by Equation (4) holds. By disintegrating and conditioning, and by using the continuity of the distribution of τ D , we obtain that The claim follows now from the Kakutani connection (Equation (3)).
The following Theorem 1 is a version of the Kakutani theorem [12] for the Yukawa equation. In some sense, it is a special case of the Kakutani connection to the Schrödinger equation studied extensively by Chung and Zhao [8]. However, it seems that this version with an unbounded and non-small domain D does not appear in any classical texts. Theorem 1. Let D ⊂ R n be a regular domain, and let f : ∂D → R be bounded and continuous.
is a solution to the Yukawa-Dirichlet problem: (ii) Moreover, if u is bounded and D is small, then Equation (6) is the only solution to the Yukawa-Dirichlet problem. (iii) As a consequence, the panharmonic measure admits the representation where h x (D; ·, ·) is the harmonic kernel defined in Equation (4).
Proof. The first and the second claim of Theorem 1 follow from the classical Kakutani theorem (cf., e.g., [15] Sections 4.4. and 4.6). Indeed, we note that the difficulties involving the Schrödinger equation in [15], Section 4.6 vanish, because To show the third claim, we condition on {τ D = t} and use the law of total probability:

Remark 1.
Unfortunately, even for very simple D , the harmonic kernel (Equation (4)) is quite difficult to find. The same is true for the regular conditional distribution (Equation (5)). For smooth boundaries ∂D , one can try the following approach: If ∂D is smooth, then the harmonic kernel h x (D; dy, t) is absolutely continuous with respect to the Lebesgue measure dy . Indeed, define p : Then p is the Brownian transition kernel: and, as a result of [16], Theorem 1, the harmonic kernel can be written as where n y is the inward normal at y ∈ ∂D and p(D; ·, ·) is the transition density of a Brownian motion that is killed when it hits the boundary ∂D , which can be written as as a result of [10], Equation (3), on page 34.
Consequently, for C 3 boundaries, the harmonic measure admits a Poisson kernel representation, and therefore, as a result of the representation given by Equation (7), the panharmonic measure also admits a Poisson kernel representation: Theorem 1 gives an interpretation of the panharmonic measure in terms of exponentially discounted Brownian motion. We give a second interpretation in terms of exponentially killed Brownian motion. Indeed, exponential discounting is closely related to exponential killing. The exponentially killed Brownian motion W µ is where † is a coffin state (by convention f ( †) = 0 for all functions f ) and Y µ is an independent exponential random variable with mean 2/µ 2 ; that is, Then we have the following representation of the panharmonic measure: Corollary 1. Let D ⊂ R n be a regular domain. Then the panharmonic measure admits the representation Proof. Let f : ∂D → R be bounded. Then, by Theorem 1 and the independence of W and Y µ , Because f was arbitrary, the claim follows.
The two representations, Theorem 1 and Corollary 1, for the panharmonic measures are, at least in spirit, classical. Now we give a third representation for the panharmonic measure in terms of an escaping Brownian motion. This representation is apparently new in spirit. The representation is due to the following Duffin correspondence Theorem 2. The functionū defined by Equation (11) is harmonic onD if and only if u is µ -panharmonic on D .
Proof. We first show that D is regular if and only ifD is regular. LetW = (W,W) be (n + 1) -dimensional Brownian motion. Denote We note that for {τ =x} to happen,x must be an endpoint of the interval I . Then, by the independence of W andW , because I is clearly regular. This shows thatD is regular if and only if D is regular.
We then show that u satisfies the Laplace equation if and only ifū satisfies the Yukawa equation; this is straightforward calculus: LetW be a 1-dimensional standard Brownian motion that is independent of W . ThenW = (W,W) is a (n + 1) -dimensional standard Brownian motion. Now the idea of how to use the Duffin correspondence is clear. We start the Brownian particlē W and count the boundary data on the side of the cylinderD = D × I , if the Brownian motion does not escape the cylinder from the bottom or from the top. In that case we count zero in the boundary; whence the name escaping Brownian motion.
Here we have chosen I = (− π 2µ , π 2µ ) in the Duffin correspondence. Consequently, all bounded solutions to the Yukawa-Dirichlet problem on a regular domain with µ 2 > 0 and continuous and bounded boundary data are given by the panharmonic measure.
Finally, we note that for a regular domain D , the domainD is regular and small. (12) is exceptionally well suited for calculations of the panharmonic measures on upper half-spaces H n + = {x ∈ R n ; x n > 0} . Indeed, Duffin ([6], Theorem 5) used it to calculate the Poisson kernel representation for panharmonic measures in the dimension n = 2 . Similar calculations can also be carried out for the general case, n ≥ 2 .

Equivalence of Harmonic and Panharmonic Measures
The probabilistic interpretation provided by Corollary 1 implies that the harmonic measure and the panharmonic measures are equivalent. Indeed, the harmonic measure counts the Brownian particles on the boundary, and the panharmonic measures count the killed Brownian particles on the boundary. However the killing happens with independent exponential random variables. Thus, if the Brownian motion can reach the boundary with positive probability, so can the killed Brownian motion, and vice versa. Additionally, it does not matter, as far as the equivalence is concerned, what the starting point is of the Brownian motion, killed or not.
Theorem 3 below makes the heuristics above precise. As corollaries of Theorem 3, we obtain a domination principle for the Dirichlet problem related to the Yukawa equation (Corollary 3) and analogs of theorems of Riesz-Riesz, Makarov and Dahlberg for the panharmonic measures (Corollary 4).
The same arguments that give the existence of the regular conditional law (Equation (5)) in the proof of Lemma 1 also give the existence and regular measurability of the following conditional Radon-Nikodym derivative: Theorem 3. Let D be a regular domain. Then all the panharmonic measures H x µ (D; ·) , µ ≥ 0, x ∈ D are mutually equivalent. The Radon-Nikodym derivative of H x µ (D; ·) with respect to H x (D; ·) is the function Z x µ (D; ·) given by Equation (13). Moreover Z x µ (D; y) is strictly decreasing in µ , and 0 < Z x µ (D; y) ≤ 1 .

Remark 3.
By Corollary 1, the Radon-Nikodym derivative Z x µ (D; ·) in Equation (13) can be interpreted as the probability that a Brownian motion killed with intensity µ 2 /2, and that would exit the domain D at y ∈ ∂D , survives to the boundary ∂D : where Y µ is an exponentially distributed random variable with mean 2/µ 2 that is independent of the Brownian motion W . To see that Z x µ (D; ·) is the Radon-Nikodym derivative, we note that, by the representation given by Equation (7) and the Kakutani connection (Equation (3)), Finally, the fact that 0 < Z x µ (D; ·) ≤ 1 is clear from the representation given by Equation (13). The fact that Z z µ (D; ·) is strictly decreasing follows immediately from the representation given by Equation (14).
From Theorem 3, we obtain immediately the following domination principle for the Dirichlet problem related to panharmonic functions: Corollary 3. Let D be a regular domain, and let u µ ≥ 0 be µ -panharmonic and u ν ≥ 0 be ν -panharmonic, respectively, on D with µ ≤ ν . Then, u ν ≤ u µ on ∂D implies u ν ≤ u µ on D .
Because domains with a rectifiable boundary are regular, we obtain immediately from Theorem 3 the following analogs of the theorems of F. Riesz and M. Riesz, Makarov and Dalhberg (see [17][18][19], respectively).

The Average Property for Panharmonic Measures and Functions
By using the representation given by Equation (7), one can calculate the panharmonic measures if one can calculate the corresponding harmonic kernels; or, equivalently, one can calculate the panharmonic measures if one can calculate the corresponding harmonic measures and the Radon-Nikodym derivatives given by Equation (13).
The harmonic kernels for balls are calculated in [16]. We do not present the general formula here. Instead, we confine ourselves to the case in which the center of the ball and the pole of the panharmonic measure coincide, and give the Gauss mean value theorem, or the average property, for panharmonic measures. As a corollary, we have the Liouville theorem for the panharmonic measures.
Let D ⊂ R n be a regular domain. For the harmonic measure, the Gauss mean value theorem states that a function u : D → R is harmonic if and only if for all balls B n (x, r) ⊂ D we have the average property: u(x) = y∈∂B n (x,r) u(y) σ n (r; dy), where σ n (r; dy) = Γ(n/2) 2π n/2 r 1−n dy is the uniform probability measure on the sphere ∂B n (x, r) .
For the panharmonic measures, the situation is similar to the harmonic measure: the only difference is that the uniform probability measure has to be replaced by a uniform sub-probability measure that depends on the killing parameter µ and the radius of the ball r . Indeed, we denote where ν = (n − 2)/2 , and is the modified Bessel function of the first kind of order ν . for all open balls B n (x, r) ⊂ D . Equivalently, H x µ (B n (x, r); dy) = ψ n (µr) σ n (r; dy).
Proof of Theorem 4. We note that we may assume x = 0 .
Denote by τ n r the first hitting time of the Brownian motion W on the boundary ∂B n (0, r) ; that is, τ n r is identical in law to the first hitting time of the Bessel process with index ν = (n − 2)/2 reaching the level r when it starts from zero.
From the rotation symmetry of the Brownian motion, it follows that the hitting place is uniformly distributed on ∂B n (0, r) for all hitting times t . Consequently, by Theorem 1 and the independence of the hitting time τ n r and place W(τ n r ) , (µr) ν 2 ν Γ(ν + 1)I ν (µr) .
The claim follows from this.

Remark 5.
The Radon-Nikodym derivative, or the "killing constant", ψ n (µ) , is rather complicated. However, some of its properties are easy to see (cf. Figure 1): The items (i)-(iv) are clear, because ψ n (µ) is the probability that an exponentially killed Brownian motion starting from the origin and with killing intensity µ 2 /2 is not killed before it hits the boundary of the unit ball. A non-probabilistic argument for (i)-(iv) is to note that ψ n (µr) = E 0 e − µ 2 2 τ n r and to use the monotone convergence. The item (v) is somewhat surprising: the higher the dimension n , the more likely it is for the killed Brownian motion to survive to the boundary of the unit ball. A possible intuitive explanation is that the higher the dimension, the more transitive the unit ball is, combined with the remarkable result by Ciesielski and Taylor [21] that the probability distribution for the total time spent in a ball by (n + 2) -dimensional Brownian motion is the same as the probability distribution of the hitting time of n -dimensional Brownian motion on the boundary of the ball. u(y) σ n (r; dy) + ψ n (µr) ∂B n (0,r) u(y) σ n (r; dy) ≤ 2ψ n (µr) u ∞ , which tends to 0 as r → ∞ by property (iii) of Remark 5. This shows that u is constant. However, for a constant u , the Yukawa equation (Equation (1)) yields 0 = µ 2 u , which implies that u ≡ 0 .

Discussion on Extensions and Simulation
The Yukawa equation (Equation (1)) is a special case of the Schrödinger equation: The Schrödinger equation and its connection to the Brownian motion has been studied, for example, by Chung and Zhao [8]. Our investigation here can be seen as a special case. For example, analogs of Theorem 1 and Corollary 1 are known for the Schrödinger equation. However, analogs of the Duffin correspondence (Equation (11)) and Corollary 2 are not known even to exist. Moreover, the results given here cannot easily be calculated for the Schrödinger equation. The problem is that the prospective Radon-Nikodym derivate of the measure associated with the solutions of the Schrödinger equation with respect to the harmonic measure takes the form where e q (t) = e − 1 2 t 0 q(W(s)) ds is the Feynman-Kac functional. Thus, we see that in order to calculate the Radon-Nikodym derivative, we need to know the joint density of the Feynman-Kac functional and the Brownian motion when the Brownian motion hits the boundary ∂D . If q is constant, that is, we have either the Yukawa equation or the Helmholtz equation, then it is enough to know the joint distribution of the hitting time and place of the Brownian motion on the boundary ∂D . These distributions are well studied (see, e.g., [14,16,[21][22][23]), but few joint distributions involving the Feynman-Kac functionals are known.
It is possible to also provide a Duffin correspondence for the Helmholtz equation. Indeed, for example, settingū (x) =ū(x,x) = u(x) cosh(λx) provides a correspondence (see [24] for details). Thus, our results extend in a straightforward manner to the Helmholtz equation (Equation (18)) for domains that are small enough with respect to the creation parameter λ that the associated Feynman-Kac functional is finite: Finally, we note that Theorem 1, Corollary 1 and Corollary 2 give three different ways to simulate the panharmonic measures. Indeed, in [24,25], the classical WOS algorithm given by Muller [26] was extended for the Yukawa PDE, and also for the Helmholtz PDE, by using the results mentioned above.
Author Contributions: Antti Rasila was responsible for the general idea of connecting the Yukawa equation to the harmonic measure and the Brownian motion through the Duffin correspondence, connections to classical analysis and potential theory, and computer experiments with Mathematica. Tommi Sottinen was responsible for the proofs and methodologies making use of stochastic analysis. Much of the paper was shaped in discussions between the authors.