# Branching random walk solutions to the Wigner equation

The stochastic solutions to the Wigner equation, which explain the nonlocal oscillatory integral operator Θ_V with an anti-symmetric kernel as the generator of two branches of jump processes, are analyzed. All existing branching random walk solutions are formulated based on the Hahn-Jordan decomposition Θ_V=Θ^+_V-Θ^-_V, i.e., treating Θ_V as the difference of two positive operators Θ^±_V, each of which characterizes the transition of states for one branch of particles. Despite the fact that the first moments of such models solve the Wigner equation, we prove that the bounds of corresponding variances grow exponentially in time with the rate depending on the upper bound of Θ^±_V, instead of Θ_V. In other words, the decay of high-frequency components is totally ignored, resulting in a severe numerical sign problem. To fully utilize such decay property, we have recourse to the stationary phase approximation for Θ_V, which captures essential contributions from the stationary phase points as well as the near-cancelation of positive and negative weights. The resulting branching random walk solutions are then proved to asymptotically solve the Wigner equation, but gain a substantial reduction in variances, thereby ameliorating the sign problem. Numerical experiments in 4-D phase space validate our theoretical findings.

## Authors

• 7 publications
• 4 publications
11/25/2020

### Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media (extended version)

In this article, we present new random walk methods to solve flow and tr...
04/07/2022

### Symmetric cooperative motion in one dimension

We explore the relationship between recursive distributional equations a...
07/19/2021

### Global random walk solvers for reactive transport and biodegradation processes in heterogeneous porous media

Flow and multicomponent reactive transport in saturated/unsaturated poro...
04/28/2021

### A Gaussian approximation theorem for Lévy processes

Without higher moment assumptions, this note establishes the decay of th...
06/02/2021

### Random walk approximation for irreversible drift-diffusion process on manifold: ergodicity, unconditional stability and convergence

Irreversible drift-diffusion processes are very common in biochemical re...
10/19/2012

### Markov Random Walk Representations with Continuous Distributions

Representations based on random walks can exploit discrete data distribu...
02/02/2022

### Well-posedness and dynamics of solutions to the generalized KdV with low power nonlinearity

We consider two types of the generalized Korteweg - de Vries equation, w...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

We are intended to discuss the probabilistic interpretation of the backward Wigner equation [1, 2, 3], arising from the recently developed particle-based simulation of the Wigner quantum dynamics [4, 5, 6, 3]. The backward Wigner equation is a partial integro-differential equation defined in phase space with an “initial” condition .

 ∂∂tφ(x,k,t)+ℏkm⋅∇xφ(x,k,t)=ΘV[φ](x,k,t),  0≤t≤T, (1.1) φ(x,k,T)=φT(x,k), (1.2)

Here is the dual Wigner function, is the mass, represents the reduced Planck constant and the pseudo-differential operator (PDO) reads

 ΘV[φ](x,k,t)=1iℏ(2π)n∫Rn×Rnei(k−k′)⋅yDV(x,y,t)φ(x,k′,t)dydk′, (1.3)

with (i.e., the central difference of the external potential ). Obviously, is anti-symmetric in -variable,

 DV(x,y,t)=−DV(x,−y,t). (1.4)

It is well known that , a nonlocal operator with an anti-symmetric symbol, actually characterizes a deformation of the classical Poisson bracket [7] and exactly reflects the nonlocal nature of quantum mechanics [8, 9, 10, 11].

The subsequent analysis will be based on two equivalent representations of the PDO. The first form is the kernel representation:

 ΘV[φ](x,k,t) =∫RnVW(x,k−k′,t)φ(x,k′,t)dk′, (1.5)

with the real-valued kernel function (termed the Wigner kernel)

 (1.6)

Here

denotes the Fourier transform of the potential function

in -variable, and

 ψ(k,t)=1iℏ(2π)neik⋅(x−z(x))Fx→kV(k,t). (1.7)

It is realized that the kernel is anti-symmetric in -variable

 VW(x,k−k′,t)=−VW(x,k′−k,t). (1.8)

due to the anti-symmetry of (see Eq. (1.4)).

The second form is the oscillatory integral representation:

 ΘV[φ](x,k,t)=∫Rneiz(x)⋅k′ψ(k′,t)(φ(x,k−k′2,t)−φ(x,k+k′2,t))dk′, (1.9)

which facilitates the derivation of its asymptotic expansion (see Theorem 2).

In order to extend to a bounded operator from to itself, say, there exists a uniform upper bound such that

 ∥ΘV[φ](t)∥2≤KV∥φ(t)∥2, (1.10)

we make the following assumptions for a finite time interval .

• and is localized in -space for any , with the minimal compact support denoted by .

• Suppose either of the following conditions holds:

• ;

• ,

and there exist a radial function , such that in and holds for sufficiently large and given constants and ;

Here is short for norm, say, .

The prototypes for the latter condition in (A2) arise from quantum molecular systems and fractional diffusion problems [12, 13]. When the potential is of the Coulomb type , it is easy to verify that and , so that the symbol functions may have singularities at and . Therefore, we need to focus on the weakly singular convolution [14], instead of solely treating it in the classical symbol class .

Now we turn to the probabilistic perspective. The starting point of the stochastic solution is to cast Eq. (1.1) into its equivalent integral formulation by adding a term on both sides of Eq. (1.1) [3],

 φ(x,k,t)=(1−G(T−t))φT(x(T−t),k)+∫TtdG(t′−t)×∫Rn(−VW(x(t′−t),k′,t′)γ0+δ(k′))φ(x(t′−t),k−k′,t′)dk′, (1.11)

the derivation of which will be put in Section 2. The constant parameter

turns out to be the intensity of an exponential distribution as follows,

 G(t′−t)=1−e−γ0(t′−t),dG(t′−t)=γ0e−γ0(t′−t),t′≥t. (1.12)

The main problem is how to resolve the negative values of kernel . In constrast to nonlocal operators with nonnegative and symmetric kernels [12, 15, 13], the existing stochastic approach is based on the unique Hahn-Jordan decomposition (HJD) [16]:

 ΘV[φ](x,k,t) =Θ+V[φ](x,k,t)−Θ−V[φ](x,k,t), (1.13) Θ±V[φ](x,k,t) =∫RnV±W(x,k−k′,t)φ(x,k′,t)dk′, (1.14) V±W(x,k,t) =max{±VW(x,k,t),0}, (1.15)

so that become positive semi-definite kernels. Moreover, we assume that there exists a uniform normalizing bound for ,

 γ0≥˘ξ=max0≤t≤Tmaxx∈X∫RnV±W(x,k,t)1{k∈2K}dk. (1.16)

It follows that the probabilistic interpretation is to seek a branching random walk model (BRW) such that its first moment satisfies the renewal-type Wigner (W) equation (1.11) [17, 5, 3], dubbed WBRW-HJD hereafter. Such model describes a mass distribution of a random cloud starting at

and frozen at random states and exhibiting both random motion and random growth. The random variable is a family history

, a denumerable random sequence corresponding to a unique family tree [18], and is the Borel extension of cylinder sets on . The particles in the family history move according to the following five rules.

• (Markov property) The motion of each particle is described by a right continuous Markov process.

• (Memoryless life-length) The particle at dies in the age time interval

with probability

.

• (Frozen state) The particle at is frozen at the state when its life-length .

• (Branching property) The particles at , carrying a weight , dies at age at state when , and produces at most five new offsprings at states , , endowed with updated weights , , respectively.

• (Independence) The only interaction between the particles is that the birth time and state of offsprings coincide with the death time and state of their parent.

We are able to define a probability measure on the measurable space and thus the stochastic process based on a specific setting of the transition kernels and particle weights in the fourth rule (vide post). Roughly speaking, WBRW-HJD can be categorized into the weighted-particle (wp) [3] and signed-particle (sp) [17, 5, 19] implementations, denoted by and associated with the probability laws and , respectively. It has been shown in [3] that (taking as an example),

 ΠwQXwt=φ(x,k,t) (1.17)

holds on some kind of probability space , where means the expectation of with respect to the probability measure

. However, to the best of our knowledge, the related variance estimation has not been established. To this end, our first contribution is to estimate the variance of WBRW-HJD, as stated in Theorem

1.

###### Theorem 1 (Variance of WBRW-HJD)

Suppose (A1) and (A2) are satisfied and let . Then the variances of and satisfy

 ∥ΠwQ(Xwt−φ(t))2∥1 ≤(1+γ1γ0(T−t))e2max(KV,˘ξ2γ0)(T−t)∥φT∥22−∥φ(t)∥22, (1.18) ∥ΠsQ(Xst−φ(t))2∥1 ≤(1+γ1γ0(T−t))e2˘ξ(T−t)∥φT∥22−∥φ(t)∥22. (1.19)

Two key observations are readily seen from Theorem 1. One is the exponential rate for spWBRW-HJD is , that depends on the volume of the support and thus cannot be improved. This poses a huge challenge for high dimensional problems since usually depends on exponentially. By contrast, the rate for wpWBRW-HJD can be reduced by increasing , and the optimal exponential rate is . Definitely, is usually far less than , implied by Eqs. (1.10) and (1.16). In this sense, the latter outperforms the former. The other is the large exponential rates and , introduced by HJD (1.13), lead to a rapid growth of variance. Such phenomenon is called “numerical sign problem” [20] as the Hahn-Jordan decomposition of a signed measure totally ignores the near-cancellation of positive and negative weights.

Our second contribution is to formulate a new class of BRW solutions, dubbed WBRW-SPA, to diminish the variance growth. The motivation comes from the stationary phase method, a useful technique in microlocal analysis [21], which makes full use of the essential contribution from the localized parts (see Theorem 2). As a consequence, the upper bounds in Eqs. (1.18) and (1.19) can be significantly reduced especially in the region where the module is sufficiently large (see Theorem 3).

###### Theorem 2 (Stationary phase approximation)

Suppose and the amplitude function . Then for a sufficiently large , we have a stationary phase approximation to PDO

 ΘV[φ](x,k,t) =Θλ0V[φ](x,k,t)+O(λ−n/20), (1.20) Θλ0V[φ](x,k,t) =Λ<λ0[φ](x,k,t)+Λ>λ0+[φ](x,k,t)+Λ>λ0−[φ](x,k,t), (1.21)

where

 Λ<λ0[φ](x,k,t) =∫B(λ0|z(x)|)eiz(x)⋅k′ψ(k′,t)Δk′[φ](x,k,t)dk′, Λ>λ0±[φ](x,k,t) =∫+∞λ0|z(x)|e±ir|z(x)|(2π±ir|z(x)|)n−12rn−1ψ(rσ±,t)Δrσ±[φ](x,k,t)dr,

in the sense that there exists a positive constant , which depends on and its first derivate but is independent on , such that

 ∥ΘV[φ](t)−Θλ0V[φ](t)∥2≤Cλ−n/20∥φ(t)∥L2x×H1k. (1.22)

Here the norm is and is a closed ball with radius centered at the origin, (short for represent two critical points on the

-dimensional unit spherical surface with normal vectors pointing in (or opposite to) the direction of

, which can be parameterized by

 σ±=(cosϑ±1,sinϑ±1cosϑ±2,…,sinϑ±1⋯sinϑ±n−2cosϑ±n−1,sinϑ±1⋯sinϑ±n−1),

with

 ϑ±i=arccot(±zi/√z2i+1+⋯+z2n)∈[0,π],  i=1,2,…,n−2,ϑ±n−1=2arccot(±(zn−1+√z2n−1+z2n)/zn)∈[0,2π), (1.23)

and is the central difference operator

 Δk′[φ](x,k,t)=φ(x,k−k′2,t)−φ(x,k+k′2,t). (1.24)

Intuitively speaking, the parameter serves as a filter to decompose PDO into a low-frequency component and a high-frequency one, the leading terms of which are , and use the resulting nonlocal operator to directly formulate WBRW-SPA, instead of as adopted in WBRW-HJD. Specifically, we still use HJD to deal with and tackle by another two branches of particles, yielding two stochastic processes: the “wp” implementation and the “sp” implementation , associated with the probability measures and , respectively. In order to estimate the effect of low-frequency parts, we further need the following assumptions.

• and is localized in -space for any with the minimal compact support denoted by ;

• For the positive constant in Eq. (1.16) there exist positive constants and such that

 α∗˘ξ=max0≤t≤Tmaxx∈X∫RnV±W(x,k,t)1{|2k|<λ0/|z(x)|}dk. (1.25)

The assumption (A4) indicates that the normalizing bound for can be diminished when is restricted in a smaller domain, which holds if is large enough. For instance, , it requires the displacement is sufficiently large. Accordingly, we are able to show that the first moment of WBRW-SPA turns out to be an asymptotic approximation to the solution of Eq. (1.1). We also study its deviation from the dual Wigner function by estimating the second moment (also termed “variance” hereafter) and find that, in contrast to Eqs. (1.18) and (1.19), the exponential growth rate in the upper bound is suppressed, so that a moderate increase of variance can be achieved.

###### Theorem 3 (Wbrw-Spa)

Suppose (A2)-(A4) are satisfied and let . Then for a sufficient large , there exist a weighted-particle branching random walk model and a signed-particle one on the probability spaces and , respectively, such that

 ΠwQYwt=ΠsQYst=φ(x,k,t)+O(λ−n/20), (1.26)

and their variances satisfy

 ∥ΠwQ(Ywt−φ(t))2∥1 ≲(1+4γ2γ0(T−t))e2max(KV,α∗˘ξ2γ0)(T−t)∥φT∥22−∥φ(t)∥22, (1.27) ∥ΠsQ(Yst−φ(t))2∥1 ≲(1+2(KV+γ2γ0)(T−t))e2α∗˘ξ(T−t)∥φT∥22−∥φ(t)∥22. (1.28)

The rest is organized as follows. Section 2 briefly reviews the basic of the Wigner equation. Section 3 derives the -boundedness and the stationary phase approximation to PDO. WBRW-HJD and WBRW-SPA are analyzed in Sections 4 and 5, respectively. In Section 6, a typical numerical experiment is performed to verify our theoretical analysis. This paper is concluded in Section 7.

## 2 The Wigner equation

The Wigner equation, introduced by Wigner in his pioneering work [1], provides a fundamental phase space description of quantum mechanics, and quantum behavior is completely characterized by the nonlocal pseudo-differential operator defined in Eq. (1.3). Mathematically speaking, it is a partial integro-differential equation defined in phase space

 ∂∂tf(x,k,t)+ℏkm⋅∇xf(x,k,t)=ΘV[f](x,k,t),  0≤t≤T,f(x,k,0)=f0(x,k), (2.1)

with an initial value . The weak formulation of the Wigner equation is of great importance since any quantum observable can be expressed by its Weyl symbol averaged by the Wigner function [8], namely, with

 ⟨f(t),φ(t)⟩=∫Rn×Rnf(x,k,t)φ(x,k,t)dxdk. (2.2)

Thus it motivates to study the dual system [2] and derive the adjoint equation of Eq. (2.1) under a non-degenerate inner product:

 (2.3)

where is a fixed time instant and is a test function with a compact support in . Using the anti-symmetry (1.8) of the Wigner kernel, we have

 ⟨ΘV[f],φ⟩=−⟨f,ΘV[φ]⟩ (2.4)

and integration by parts directly leads to

 ⟨∂f∂t+ℏkm⋅∇xf−ΘV[f],φ⟩T=⟨∂f∂t,φ⟩T+⟨ℏkm⋅∇xf,φ⟩T−⟨ΘV[f],φ⟩T=⟨fT,φT⟩−⟨f0,φ0⟩−⟨f,∂φ∂t+ℏkm⋅∇xφ−ΘV[φ]⟩T,

where and are short for and , respectively. Therefore the adjoint correspondence, i.e., the backward Wigner equation (1.1), is immediately derived by setting

 ⟨φT,fT⟩=⟨φ0,f0⟩. (2.5)

Formally, Eq. (2.5) allows us to evaluate the quantum mechanical observable only by the “initial” data [3].

The backward Wigner equation (1.1) can be cast into a renewal-type equation by adding a term on both sides,

 ∂∂tφ(x,k,t)+ℏkm⋅∇xφ(x,k,t)−γ0⋅φ(x,k,t)=ΘV[φ](x,k,t)−γ0⋅φ(x,k,t), (2.6)

with being a prescribed constant (see Eq. (1.16)), and the mild solution reads

 φ(x,k,t)=e(T−t)AφT(x,k)−∫Tte(t′−t)A(ΘV[φ](x,k,t′)−γ0⋅φ(x,k,t′))dt′,

by the variation-of-constant formula [22]. Here is short for the semigroup generated by , and its action on a given function can be now readily performed. For instance, we have

 e(t′−t)Ag(x,k,t′)=e−γ0(t′−t)g(x(t′−t),k,t′),t′≥t, (2.7)

where gives the forward-in-time trajectory of with a positive time increment . That is, the backward renewal-type equation Eq. (1.11) is thus verified.

## 3 L2-boundedness and stationary phase approximation

Before proceeding to the probabilistic aspect, we first need to establish the -boundedness of under the assumptions (A1) and (A2). For , Eq. (1.10) is readily verified by Young’s convolution inequality, whereas the -boundedness for weakly singular kernels is obtained by the Hardy-Littlewood-Sobolev theorem [23]. After that, we present the stationary phase approximation and detail its remainder estimate.

Suppose and the second condition of hold, then for due to the Hölder’s inequality:

 ∥φ(t)∥p=∥φ(t)⋅1X⋅1K∥p≤∥1X⋅1K∥2p2−p∥φ(t)∥2<∞ (3.1)

as has a finite measure. Next we introduce a smooth cut-off function :

 χε,R(r)={1,r∈[ε,2R],0,r∈[0,ε/2)∪(3R,+∞), (3.2)

and let , . Here is introduced to remove the singularity at and is chosen sufficient large to ensure , as stated in assumption . Then it is readily verified that the truncated operator has the following estimate

 ∥ΘεV[φ](t)∥L2k≤2n+1∥∫Rne2i(k−k′)⋅z(x)(ψε+ψ∞)(2(k−k′),t)φ(x,k′,t)dk′∥L2k≤2n+1∥∫Rnψε(2(k−k′),t)(e−2ik′⋅z(x)φ(x,k′,t))dk′∥L2k+2n+1∥∫Rnψ∞(2(k−k′),t)(e−2ik′⋅z(x)φ(x,k′,t))dk′∥L2k. (3.3)

The first term is bounded from to itself as is locally integrable, say,

 ∥∫Rnψε(2(k−k′),t)(e−2ik′⋅z(x)φ(x,k′,t))dk′∥L2k≤∥Ψ⋅χε,R∥L1k⋅∥φ(t)∥L2k, (3.4)

and the bound is independent of . The second term is also bounded from to , with , owing to the Hardy-Littlewood-Sobelev theorem,

 ∥∫Rnψ∞(2(k−k′),t)(e−2ik′⋅z(x)φ(x,k′,t))dk′∥L2k≤∥∫Rnχ2R,+∞(2|k−k′|)⋅|φ(x,k′,t)|2n−α|k−k′|n−αdk′∥L2k≤Cp∥φ(t)∥Lpk≤~Cp∥φ(t)∥L2k. (3.5)

Now let in Eq. (3.3). By combining Eq. (3.1), we obtain that there exists a uniform such that

 ∥ΘV[φ](t)∥2≤KV∥φ(t)∥2. (3.6)

A remarkable feature of the oscillatory integral operator is the decay property as the integrand becomes more and more oscillating. As stated by Hörmander’s theorem [21, 24], when is sufficiently smooth and compactly supported, it has a sharp estimate for a sufficiently large

 ∥ΘV[φ](x,k,t)∥L2k≤C|z(x)|−n/2∥φ(t)∥L2k. (3.7)

The physical meaning of Eq. (3.7) is also clear. When we consider the two-body interacting potential like , turns out to be the spatial displacement between two bodies, so that the estimate (3.7) characterizes the decay rate of quantum interaction as the distance increases. A similar result like Eq. (3.7) can also be found in our framework, and the decay property will be fully utilized by the stationary phase method as presented in Theorem 2, which is definitely ignored by HJD (1.13).

###### Proof (Proof of Theorem 2)

It starts by splitting the wavevector into its modulus and orientation parts with the modulus and the orientation , where denotes the -dimensional spherical surface, and focusing on the high-frequency component

 ΘV[φ]−Λ<λ0[φ]=∫+∞λ0|z(x)|rn−1dr∫Sn−1dσ eir|z(x)|z′⋅σψ(rσ,t)Δrσ[φ](x,k,t), (3.8)

where represents the orientation of , and denotes the induced Lebesgue measure on . After choosing the equatorial plane normal to , the unit sphere can be decomposed into an upper hemisphere and a lower one satisfying . Accordingly, the inner surface integral of the first kind over in Eq. (3.8) equals to the sum of those over and . Without loss of generality, it suffices to assume that , which be realized by a rotation otherwise. Let us start from the graph

 Sn−1±={σ∈Rn∣∣σn=±ϕ(σ1,…,σn−1),σi∈[−1,1],i=1,…,n−1}, (3.9)

with

 ϕ(σ1,…,σn−1)=√1−σ21−⋯−σ2n−1 (3.10)

and take the surface integral of the first kind over the upper hemisphere as an example. Now the phase function of the integrand becomes

 S(x,r,σ1,…,σn−1)=r|z(x)|z′⋅σ=r|z(x)|ϕ(σ1,…,σn−1). (3.11)

For such phase function, it can be easily verified that there is only one critical point satisfying , and the determinant of its Hessian matrix at turns out to be

 det(Hess(S