Wigner transport in linear electromagnetic fields

The analysis of charged quantum particles in electromagnetic fields is, among others, particularly important to nanoelectronics [1–8]. The established Wigner formulation of quantum mechanics [9] (see recent reviews [10, 11] and book [12]) defines the Wigner function by applying the Weyl transform to the density matrix [13]:

$$\begin{align} f_{\textrm{w}}\left({\textbf{p}},{\textbf{x}}\right) = \int \frac{\textrm{d}{\textbf{s}}}{\left(2\pi\hbar\right)^3}\textrm{e}^{-\frac{\textrm{i}}{\hbar}{\textbf{s}}\cdot{\textbf{p}}}\rho\left({\textbf{x}}+\frac{{\textbf{s}}}{2},{\textbf{x}}-\frac{{\textbf{s}}}{2}\right) .\end{align} \tag{ 1 }$$

The density matrix ρ of a pure state is defined from the solution ψ of the Schrödinger equation as $\rho(\textbf{x},\textbf{y}) = \psi(\textbf{x})\psi^*(\textbf{y})$ and depends on two position variables. (1) is a transformation from the position space to the phase space, i.e. $f_\textrm{w}$ is a function of the momentum $\textbf{p}$ and the position $\textbf{x}$. The evolution equation for the Wigner function is obtained by applying the Weyl transform to the von Neumann equation $\textrm{i}\hbar\frac{\partial}{\partial t}\hat\rho = [\hat H,\hat\rho]_-: = \hat H\hat\rho-\hat\rho \hat H,$ with the Hamiltonian $\hat H = \frac{1}{2m}\left(-\textrm{i}\hbar\nabla\right)^2+V(\textbf{r})$ [14]. The potential energy V defines a central quantity of the standard theory, namely the Wigner potential:

$$\begin{align} V_{\textrm{w}}\left({\textbf{p}},{\textbf{x}}\right) = \frac{1}{\left(2\pi\hbar\right)^3}\int \frac{\textrm{d}{\textbf{s}}}{\textrm{i}\hbar}\textrm{e}^{-\frac{\textrm{i}}{\hbar}{\textbf{s}}\cdot\mathbf{p}}\left[V\left({\textbf{x}}+\frac{{\textbf{s}}}{2}\right)-V\left({\textbf{x}}-\frac{{\textbf{s}}}{2}\right)\right] .\end{align} \tag{ 2 }$$

The scalar potential $\phi = V/e$, with the electron charge e, and the canonical momentum operator $-\textrm{i}\hbar\nabla$, are fundamental for this picture. The choice of the gauge is implicitly assumed, i.e. the vector potential $\textbf{A}$ is chosen to be zero. However, any other couple $\textbf{A}^{^{\prime}},\phi^{^{\prime}}$ satisfying $\textbf{A}^{^{\prime}} = \textbf{A}+\nabla\chi,\quad\phi^{^{\prime}} = \phi-\partial\chi/\partial t$ for a given function χ modifies the Hamiltonian and may lead to a very different physical picture, despite that the electromagnetic environment $\textbf{B} = \nabla\times\textbf{A}$, $\textbf{E} = -\nabla\phi-\partial\textbf{A}/\partial t$ remains independent on χ [15]. An example is related to electrons governed by an electric field $\textbf{E}$ [16] in a periodic potential. If Wannier-Stark localized states [17] are used for the description, the picture involves a discrete energy spectrum accounting for the translational crystal symmetry. If accelerated Bloch states (Houston states) [18] are used, the picture of continuous acceleration of the wave vector in the crystal band structure gives rise to a periodic electron motion, called Bloch oscillations. It has been shown that the two pictures are equivalent and related to the choice of a vector $(\textbf{A} = -\textbf{E} t;~\phi = 0)$, or a scalar potential gauge $(\textbf{A} = 0, \phi = -\textbf{E}\textbf{x})$, linked by $\chi = -\textbf{E}\textbf{x} t$ [19, 20]. For the standard Wigner picture, the zero vector potential is a convenient choice, because then the canonical momentum $\textbf{p}$ and the kinetic momentum $\textbf{P}$ coincide. This is not true anymore in the case of a magnetic field when $\textbf{P} = \textbf{p}-e\textbf{A}(\textbf{x})$. In this case, using the kinetic momentum as a phase space variable offers the advantage that the latter is a physical quantity and thus gauge-invariant [21–26]. Inspired by this fact, Stratonovich [27] generalized the Weyl transform to

$$\begin{align} f_{\textrm{w}}\left({\textbf{P}},{\textbf{x}}\right) = \int \frac{\textrm{d}{\textbf{s}}}{\left(2\pi\hbar\right)^3}\textrm{e}^{-\frac{\textrm{i}}{\hbar}{\textbf{s}}\cdot\left[{\textbf{P}}+\frac{e}{2}\int_{-1}^1\textrm{d}\tau{\textbf{A}}\left({\textbf{x}}+\frac{{\textbf{s}}\tau}{2}\right)\right]}\rho\left({\textbf{x}}+\frac{{\textbf{s}}}{2},{\textbf{x}}-\frac{{\textbf{s}}}{2}\right). \end{align} \tag{ 3 }$$

Now the transform depends on the vector potential, however, the evolution equation for the Wigner function regarding the position and the kinetic momentum depends only on the electromagnetic field $\textbf{E}$, $\textbf{B}$ [28]. Thus, the Weyl–Stratonovich transform lifts the gauge dependence, offering more physical transparency to the quantum evolution. In the case $\textbf{A} = 0$, the Weyl–Stratonovich transform equals the Weyl transform and can thus be seen as an extension. For the sake of convenience, we use $\textbf{p}$ instead of $\textbf{P}$ to refer to the kinetic momentum for the remainder of this work.

There are two ways to formulate the evolution equation depending on the physical settings. If the physical system is bounded in space, in a domain enclosed in $(-\textbf{L}/2,\textbf{L}/2)$, where $\textbf{L}$ is called coherence length, the momentum space becomes discrete, involving the integer variable $\textbf{m}$: $\textbf{P}_\textbf{m} = \textbf{m}\Delta \textbf{P},\quad \textbf{m}\in \mathbb{Z}\times\mathbb{Z},\quad \Delta \textbf{P} = 2\pi\hbar/\textbf{L}$. In the limit $\textbf{L}\to\infty$, called long coherence length limit, the momentum becomes continuous [29]. For electromagnetic fields with general spatiotemporal dependence, both formulations are very challenging from a numerical point of view. A computational experience with the treatment of multidimensional sums and integrals is missing. To gain first experiences, we look for simplified physical conditions to reduce the equation's complexity, allowing in particular for the application of analytical approaches. The fact that for a homogeneous magnetic field certain integrals vanish is helpful for choosing such conditions, while the field appears as the magnetic component of the Lorentz force in the Liouville operator of the reduced equation. This prompts considering the next term in the Taylor expansion of the magnetic field $\textbf{B}(\textbf{x})$, namely linearly dependent magnetic fields. Furthermore, in the case of linear electric fields, they complete the force term in the Liouville operator to a full Lorentz force. We can thus formulate the physical settings under consideration: We consider a transport in a two-dimensional (2D) plane with coordinates $\textbf{x} = (x,y,0)^T$. A magnetic field $\textbf{B}(y) = (0,0,B_0+B_1y)^T$ points perpendicular to the plane and depends linearly on y. The electric field $\textbf{E}(x,y) = (E_xx,E_yy,0)^T$ accelerates the electron in the plane. The obtained equation using the long coherence length limit [29] is given by

$$\begin{align} \left(\frac{\partial}{\partial t} + \frac{{\textbf{p}}}{m}\cdot \frac{\partial}{\partial {\textbf{x}}} +{\textbf{F}} \cdot\frac{\partial}{\partial {\textbf{p}}}\right)f_{\textrm{w}}\left({\textbf{p}},{\textbf{x}}\right) = \frac {B_1\hbar^2}m\frac e{12} \left( \frac{\partial^2}{\partial p_y^2}\frac{\partial}{\partial x} - \frac{\partial}{\partial p_x}\frac{\partial}{\partial p_y} \frac{\partial}{\partial y}\right)f_{\textrm{w}}\left({\textbf{p}},{\textbf{x}}\right). \end{align} \tag{ 4 }$$

We note that the Lorentz force $\textbf{F} = e[\textbf{E}(x,y)+\textbf{p}\times\textbf{B}(y)/m]$ in the Liouville operator on the left depends on the electromagnetic field. The operator corresponds to a classical motion over Newtonian trajectories, accelerated by the Lorentz force, linearly dependent on the position coordinates. The term on the right-hand side depends only on the magnetic field gradient B₁ and consistently vanishes if $B_1\to 0$. This term is responsible for the quantum character of the evolution process. Indeed, the structure of (4) resembles the standard Wigner equation. The latter consists of the forceless Liouville operator, whose interplay with the Wigner potential term gives rise to a fully quantum-coherent evolution. Indeed, the equation is equivalent to the von Neumann equation and in a pure state to the Schrödinger equation [13, 30]. However, this term is given by the convolution of the Wigner function with $V_\textrm{w}$ in (2) and thus depends linearly on $f_\textrm{w}$. The corresponding term in (4) introduces high-order mixed derivatives and hence has different numerical aspects. The numerical experience with the former equation has matured for more than three decades [31–36]. Furthermore, a peculiarity of phase space formulations of quantum mechanics is the ability to use them for further development of heuristic, physics-based models, associated with quantum phenomena and processes. Good examples are quantum particle models where particles are provided with additional attributes, such as sign or affinity, while the action of the electric potential is interpreted as scattering or as particle generation [37]. In contrast, alternative quantum theories associate physical quantities and quantum processes with formal mathematical expressions, which offer little physical insight (e.g. operator mechanics).

This work provides a numerical analysis of (4) and a particle picture with the corresponding quantum evolution. These quantum particles have a numerical origin, however, they bear the basic properties of the physical models of particles in classical mechanics. The additional particle properties carry the quantum information of the evolution.

In section 2, an iterative solution to (4) is presented. The strategy is based on transforming the equation to a Fredholm integral equation, which can be solved by a resolvent expansion. In section 3, we derive two different Monte Carlo algorithms for the evaluation of the terms in the resolvent expansion. In section 4, the key findings of this work are discussed.

We introduce two new time-dependent functions of the Newtonian trajectory, which replace the phase space variables. We use two parameterizations (backward and forward), which yield different representations of the same solution. This is followed by transforming the gauge-invariant Wigner equation (4) into an integral form, i.e. the Fredholm integral equation, by using a finite difference scheme and a resolvent expansion of the Wigner function. We first present the solution for the backward parameterization and afterward for the forward parameterization. For the latter, we define and solve the adjoint formulation of the Fredholm equation. Finally, both solutions are used to evaluate the expectation value of a physical quantity A iteratively.

2.1. Newtonian trajectories with backward and forward parameterization

The two new time-dependent functions of the Newtonian trajectory are based on the actual physical behavior of an electron governed by the Lorentz force $\textbf{F}$. The parameterization can be done backward and forward in time.

2.1.1. Backward parameterization.

Consider a particle at a time t, the position $\textbf{x}$, and the momentum $\textbf{p}$ as initial values in a force field $\textbf{F}$. From there, one can determine the position and momentum at an earlier time $t^{^{\prime}}\lt t$. They are given by the two integral equations

$$\begin{align} {\textbf{x}}\left(t^{^{\prime}};{\textbf{p}},{\textbf{x}},t\right): &= {\textbf{x}} - \int_{t^{^{\prime}}}^t\frac{{\textbf{p}}\left(\tau;{\textbf{p}},{\textbf{x}},t\right)}{m}\textrm{d}\tau, \end{align} \tag{ 5 }$$

$$\begin{align} {\textbf{p}}\left(t^{^{\prime}};{\textbf{p}},{\textbf{x}},t\right):& = {\textbf{p}}-\int_{t^{^{\prime}}}^t{\textbf{F}}\left({\textbf{p}}\left(\tau;{\textbf{p}},{\textbf{x}},t\right), {\textbf{x}}\left(\tau;{\textbf{p}},{\textbf{x}},t\right)\right)\textrm{d}\tau. \end{align} \tag{ 6 }$$

2.1.2. Forward parameterization.

In this case, the particle is initialized at $t^{^{\prime}},\textbf{p}^{^{\prime}},\textbf{x}^{^{\prime}}$. $\textbf{p}$ and $\textbf{x}$ are then evaluated at a later time $t\gt t^{^{\prime}}$ as

$$\begin{align} {\textbf{x}}^{^{\prime}}\left(t;{\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right): &= {\textbf{x}}^{^{\prime}}+ \int_{t^{^{\prime}}}^t\frac{{\textbf{p}}^{^{\prime}}\left(\tau;{\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right)}{m}\textrm{d}\tau, \end{align} \tag{ 7 }$$

$$\begin{align} {\textbf{p}}^{^{\prime}}\left(t;{\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right): &= {\textbf{p}}^{^{\prime}}+\int_{t^{^{\prime}}}^t{\textbf{F}}\left({\textbf{p}}^{^{\prime}}\left(\tau;{\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right), {\textbf{x}}^{^{\prime}}\left(\tau;{\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right)\right)\textrm{d}\tau. \end{align} \tag{ 8 }$$

For convenience, we will write $\textbf{x}(t^{^{\prime}}), \textbf{p}(t^{^{\prime}})$ and $\textbf{x}^{^{\prime}}(t), \textbf{p}^{^{\prime}}(t)$ respectively. We also will use the Liouville theorem, stating that the phase space volume remains constant along the trajectories of the system, i.e. $\int \textrm{d}\textbf{p} \textrm{d}\textbf{x} = \int \textrm{d}\textbf{p}(t^{^{\prime}}) \textrm{d}\textbf{x}(t^{^{\prime}}) = \int \textrm{d}\textbf{p}^{^{\prime}}(t) \textrm{d}\textbf{x}^{^{\prime}}(t)$.

2.2. Fredholm integral representation of the gauge-invariant Wigner equation

Next, we show how the gauge-invariant Wigner equation is transformed into an integral form, i.e. the Fredholm integral equation. For this purpose, a finite difference scheme is used to replace the derivatives.

2.2.1. Integral form.

For the transformation, the variables $\textbf{x}$ and $\textbf{p}$ in (4) are replaced by the functions (5) and (6), respectively. That way, the Liouville operator on the left-hand side can be replaced by a total derivative of time and integrated on tʹ in the limits $(t_0,t)$. By setting t₀ = 0 (i.e. the time when the initial condition $f_{\textrm{w}_0}$ is known) it is obtained

$$\begin{align} f_{\textrm{w}}\left({\textbf{p}},{\textbf{x}},t\right) = & ~\textrm{e}^{-\int\limits^t_{0}\gamma\left({\textbf{p}}\left(\tau\right),{\textbf{x}}\left(\tau\right)\right)\textrm{d}\tau} f_{{\textrm{w}}_0}\left({\textbf{p}}\left(0\right),{\textbf{x}}\left(0\right)\right) + \int_0^{t} \textrm{d} t^{^{\prime}} \textrm{e}^{-\int\limits^t_{t^{^{\prime}}}\gamma\left({\textbf{p}}\left(\tau\right),{\textbf{x}}\left(\tau\right)\right)\textrm{d}\tau} \nonumber\\ & \cdot\left[ \frac {B_1\hbar^2}m\frac e{12} \left( \frac{\partial^3}{\partial p_y^2\partial x} -\frac{\partial^3}{\partial p_x\partial p_y\partial y} \right) +\gamma\left({\textbf{p}}\left(t^{^{\prime}}\right),{\textbf{x}}\left(t^{^{\prime}}\right)\right) \right]f_{\textrm{w}}\left({\textbf{p}}\left(t^{^{\prime}}\right),{\textbf{x}}\left(t^{^{\prime}}\right),t^{^{\prime}}\right).\nonumber\\ \end{align} \tag{ 9 }$$

Here, γ is an auxiliary function, which is not presented in the differential form of the equation. Indeed, after taking the derivative with respect to t₀, the terms containing γ cancel exactly. Later, we show that the introduction of γ is convenient from a numerical point of view and also has a physical meaning in the quantum particle model under development.

By taking a closer look at (9) we can gain insights into the physical background. The linear coefficient B₁ of the magnetic field determines the quantum character of the evolution. Consider the case where B₁ = 0 and γ = 0. The equation then simplifies to $f_\textrm{w}\bigl(\textbf{p},\textbf{x},t\bigr) = f_{\textrm{w}_0}\bigl(\textbf{p}(0),\textbf{x}(0)\bigr)$. This means that the Wigner function is constant along the trajectories of the system and one can evaluate $f_\textrm{w}$ at any time t by tracing the trajectory back to t = 0, which is in accordance with Liouville's theorem. Indeed, an initial classical particle density in $\textrm{d}\textbf{x}(0)\textrm{d}\textbf{p}(0)$ evolves along the trajectories until time t without any change.

2.2.2. Finite difference scheme.

The integral equation (9) is not yet of Fredholm type as it contains derivatives of the integrand function $f_\textrm{w}$. However, they can be approached by a finite difference scheme, which replaces them with linear combinations of $f_\textrm{w}$ defined in adjacent phase space points. Here, we apply a central finite difference scheme. This leads to fifteen terms represented by the indices $\textbf{i} = (i_x,i_y),\textbf{j} = (j_x,j_y)$ and coefficients $\alpha_{{\textbf{i}}{\textbf{j}}}$, where $i_x,i_y,j_x,j_y\in\{-1,0,1\}$. We also choose γ to be a constant:

$$\begin{align} \gamma = \gamma\left({\textbf{p}}\left(t^{^{\prime}}\right),{\textbf{x}}\left(t^{^{\prime}}\right)\right): = \frac {B_1\hbar^2}m\frac e{96\left(\Delta P\right)^2\Delta X} = \textrm{constant} .\end{align} \tag{ 10 }$$

The convenience of this choice will be discussed below. With the help of integrals over $\textbf{p}$ and $\textbf{x}$, and the use of δ functions the equation obtains a mathematically formal appearance:

$$\begin{align} f_{\textrm{w}}\left({\textbf{p}},{\textbf{x}},t\right) &= f_{\textrm{i}}\left({\textbf{p}},{\textbf{x}},t\right)+\int_0^\infty \textrm{d} t^{^{\prime}}\int \textrm{d}{\textbf{p}}^{^{\prime}}\int \textrm{d}{\textbf{x}}^{^{\prime}} \mathcal{K}\left({\textbf{p}},{\textbf{x}},t,{\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right)f_{\textrm{w}}\left({\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right),\nonumber\\ f_{\textrm{i}}\left({\textbf{p}},{\textbf{x}},t\right) &= \textrm{e}^{-t\gamma}f_{{\textrm{w}}_0}\left({\textbf{p}}\left(0\right),{\textbf{x}}\left(0\right)\right),\nonumber\\ \mathcal{K}\left({\textbf{p}},{\textbf{x}},t,{\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right) &= ~ \theta\left(t-t^{^{\prime}}\right)\gamma\textrm{e}^{-\left(t-t^{^{\prime}}\right)\gamma}\sum_{{\textbf{i}},{\textbf{j}}}\alpha_{{\textbf{i}}{\textbf{j}}}\delta\left({\textbf{p}}\left(t^{^{\prime}}\right)+{\textbf{i}}\Delta P-{\textbf{p}}^{^{\prime}}, {\textbf{x}}\left(t^{^{\prime}}\right)+{\textbf{j}}\Delta X-{\textbf{x}}^{^{\prime}}\right). \end{align} \tag{ 11 }$$

The Heaviside function on time takes care of the proper upper limit t. The detailed form of the kernel $\mathcal{K}$ can be found in appendix A.

2.3. Solution of the Fredholm integral equation

In this section, we present a solution for (11) and how it can be used to evaluate the expectation value of a physical quantity A of a particle. The weak formulation of this task is given as a series of integrals. This series arises from the resolvent expansion of the Wigner function. Consequently, the solution for the physical quantity is done iteratively.

2.3.1. Weak formulation of the task.

The Wigner function is a quasi-distribution function and can be used as a probability density for quantum particles [13]. Consider an arbitrary physical quantity A, which depends on position, momentum, and time. The expectation value of A at a time T can be evaluated by

$$\begin{align} \langle A\rangle\left(T\right) = \int_0^\infty \textrm{d} t\int \textrm{d}{\textbf{p}}\int \textrm{d}{\textbf{x}} f_{\textrm{w}}\left({\textbf{p}},{\textbf{x}},t\right)A\left({\textbf{p}},{\textbf{x}},t\right)\delta\left(T-t\right). \end{align} \tag{ 12 }$$

For convenience reasons we set $A_T(\textbf{p},\textbf{x},t): = A(\textbf{p},\textbf{x},t)\delta(T-t)$. The solution of Fredholm integral equations is presented by its resolvent expansion [38], as given in appendix B. It allows to represent $\langle A\rangle(T)$ as a series

$$\begin{align} \langle A\rangle\left(T\right) = \sum_{n = 0}^\infty\int_0^\infty \textrm{d} t\int \textrm{d}{\textbf{p}}\int \textrm{d}{\textbf{x}} f_n\left({\textbf{p}},{\textbf{x}},t\right)A_T\left({\textbf{p}},{\textbf{x}},t\right) = \sum_{n = 0}^\infty\langle A\rangle_n\left(T\right). \end{align} \tag{ 13 }$$

In particular, if A is chosen to be a delta function, the series yields the expansion of the Wigner function.

2.3.2. Resolvent expansion of the Wigner function.

Given the scattering indices $(\textbf{i}_k,\textbf{j}_k)_{1\unicode{x2A7D} k\unicode{x2A7D} n}$ and the scattering times $t_1\lt t_2\lt\ldots\lt t_n$, we introduce the trajectory with scattering events for backward parameterization as

$$\begin{align} {\textbf{p}}_n\left(t^{^{\prime}}\right): = \left\{\begin{array}{ll}{\textbf{p}}_{n-1}\left(t^{^{\prime}}\right) &~\textrm{for}~t_n\lt t^{^{\prime}}\unicode{x2A7D} T\\ {\textbf{p}}\left(t^{^{\prime}};{\textbf{p}}_{n-1}\left(t_{n}\right)+{\textbf{i}}_n\Delta P,{\textbf{x}}_{n-1}\left(t_{n}\right)+{\textbf{j}}_n\Delta X,t_n\right)&~\textrm{for}~0\unicode{x2A7D} t^{^{\prime}}\unicode{x2A7D} t_n,\\ \end{array}\right.\nonumber\\ {\textbf{x}}_n\left(t^{^{\prime}}\right): = \left\{\begin{array}{ll}{\textbf{x}}_{n-1}\left(t^{^{\prime}}\right) &~\textrm{for}~t_n\lt t^{^{\prime}}\unicode{x2A7D} T\\ {\textbf{x}}\left(t^{^{\prime}};{\textbf{x}}_{n-1}\left(t_{n}\right)+{\textbf{i}}_n\Delta P,{\textbf{x}}_{n-1}\left(t_{n}\right)+{\textbf{j}}_n\Delta X,t_n\right)&~\textrm{for}~0\unicode{x2A7D} t^{^{\prime}}\unicode{x2A7D} t_n,\\ \end{array}\right. \end{align} \tag{ 14 }$$

where we use the convention $\textbf{p}_0(t^{^{\prime}}): = \textbf{p}(t^{^{\prime}};\textbf{p},\textbf{x},T)$, $\textbf{x}_0(t^{^{\prime}}): = \textbf{x}(t^{^{\prime}};\textbf{p},\textbf{x},T)$.

In accordance with (B.1) we obtain

$$\begin{align} f_n\left({\textbf{p}},{\textbf{x}},t\right) = \gamma^n\textrm{e}^{-\gamma t}\int_0^t\textrm{d} t_1\int_0^{t_1}\textrm{d} t_2\ldots\int_0^{t_{n-1}}\textrm{d} t_n\sum_{{\textbf{i}}_1,{\textbf{j}}_1}\ldots\sum_{{\textbf{i}}_n,{\textbf{j}}_n}\prod_{k = 1}^n\left(\alpha_{{\textbf{i}}_k{\textbf{j}}_k}\right)f_{{\textrm{w}}_0}\left({\textbf{p}}_n\left(0\right),{\textbf{x}}_n\left(0\right)\right),\nonumber\\ \end{align} \tag{ 15 }$$

where $f_0(\textbf{p},\textbf{x},t) = \textrm{e}^{-t\gamma}f_{\textrm{w}_0}\big(\textbf{p}(0),\textbf{x}(0)\big)$, see appendix C.

The existence of the backward Newtonian trajectories invokes a picture of a pointlike particle that evolves back in time. The delta functions, which give rise to offsets of the phase space positions, can be interpreted as scattering factors. Figure 1 schematically presents the second term in the iterative expansion of the Wigner function. The particle starts at $(\textbf{p},\textbf{x},T)$ and moves back in time in the phase space according to the Lorentz force $\textbf{F}$. When the particle reaches t₁ it is scattered, i.e. a factor $(\textbf{i}_1\Delta P,\textbf{j}_1\Delta X)$ is added. Next, it follows the trajectory again until it reaches t₂. This process is repeated until t = 0 is reached.

Zoom In Zoom Out Reset image size

2.3.3. Iterative representation of physical quantities.

To evaluate the solution of $\langle A\rangle_n(T)$, we insert the solution of $f_n,n\in\mathbb{N}$ in (15) into (13). This yields

$$\begin{align} \langle A\rangle_n\left(T\right) &= \gamma^n\textrm{e}^{-T\gamma}\int \textrm{d}{\textbf{p}}\int \textrm{d}{\textbf{x}} \int_0^T \textrm{d} t_1\int_0^{t_1} \textrm{d} t_2\ldots\int_0^{t_{n-1}} \textrm{d} t_nA\left({\textbf{p}},{\textbf{x}},T\right)\nonumber\\ &\quad\cdot\sum_{{\textbf{i}}_1,{\textbf{j}}_1}\ldots\sum_{{\textbf{i}}_n,{\textbf{j}}_n}\prod_{k = 1}^n\left(\alpha_{{\textbf{i}}_k{\textbf{j}}_k}\right)f_{{\textrm{w}}_0}\left({\textbf{p}}_n\left(0\right),{\textbf{x}}_n\left(0\right)\right). \end{align} \tag{ 16 }$$

This shows us how each element $\langle A\rangle_n(T)$ is generated. In the backward parameterization case, the trajectory of $\textbf{p}$ and $\textbf{x}$ starts at T and goes back in time, according to (14). The particle is scattered at each $(t_i)_{i\in\{1,2,\ldots,n\}}$, where $T\gt t_1\gt t_2\gt\ldots\gt t_n\gt 0$. The indices $\textbf{i}_k$ and $\textbf{j}_k$ are implicitly included in the functions $\textbf{p}_n$ and $\textbf{x}_n$. Reaching the final momentum and position at t = 0, they are used as the arguments of the initial condition of the Wigner function $f_{\textrm{w}_0}$. The integration limits of the t_i 's and consequently their orders are determined by the θ functions of the kernel.

2.4. Solution of the adjoint integral equation

In this section, a solution of the Fredholm integral equation (11) is presented where forward parameterization is used. The weak formulation of this task is given by the adjoint formulation of the Fredholm integral equation. Finally, the solution for the adjoint equation is used to derive the expectation value of a physical quantity iteratively.

2.4.1. Weak formulation of the task.

The adjoint of a Fredholm integral equation has the same kernel, but the integration is over the other set of variables:

$$\begin{align} g\left({\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right) = g_{\textrm{i}}\left({\textbf{p}}^{^{\prime}},{\textbf{x}},^{^{\prime}}t^{^{\prime}}\right)+\int_{0}^\infty \textrm{d} t\int_{-\infty}^\infty \textrm{d}{\textbf{p}}\int_{-\infty}^\infty \textrm{d}{\textbf{x}}\mathcal{K}\left({\textbf{p}},{\textbf{x}},t,{\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right)g\left({\textbf{p}},{\textbf{x}},t\right) .\end{align} \tag{ 17 }$$

The free term $g_\textrm{i}$ can be determined from the weak formulation of the task, namely to find the expecation value of a physical quantity A. The following relation follows from the exchange lemma in appendix B.2 and the Liouville theorem. By choosing $g_\textrm{i}(\textbf{p}^{^{\prime}},\textbf{x}^{^{\prime}},t^{^{\prime}}): = A_T(\textbf{p}^{^{\prime}},\textbf{x}^{^{\prime}},t^{^{\prime}})$ we can show

$$\begin{align} \langle A\rangle\left(T\right) &= \int\limits_0^\infty \textrm{d} t\int \textrm{d}{\textbf{p}}\int \textrm{d}{\textbf{x}} f_{\textrm{w}}\left({\textbf{p}},{\textbf{x}},t\right)A_T\left({\textbf{p}},{\textbf{x}},t\right) = \int\limits_0^\infty \textrm{d} t\int \textrm{d}{\textbf{p}}\int \textrm{d}{\textbf{x}} f_{\textrm{w}}\left({\textbf{p}},{\textbf{x}},t\right)g_{\textrm{i}}\left({\textbf{p}},{\textbf{x}},t\right) \nonumber\\ & = \int\limits_0^\infty \textrm{d} t\int \textrm{d}{\textbf{p}}\int \textrm{d}{\textbf{x}} f_{\textrm{i}}\left({\textbf{p}},{\textbf{x}},t\right)g\left({\textbf{p}},{\textbf{x}},t\right) = \int\limits_0^\infty \textrm{d} t\int \textrm{d}{\textbf{p}}\int \textrm{d}{\textbf{x}} \textrm{e}^{-t\gamma}f_{{\textrm{w}}_0}\left({\textbf{p}},{\textbf{x}}\right)g\left({\textbf{p}}^{^{\prime}}\left(t\right),{\textbf{x}}^{^{\prime}}\left(t\right),t\right). \end{align} \tag{ 18 }$$

Like before, we consider the resolvent expansion to evaluate $\langle A\rangle(T)$, which yields

$$\begin{align} \langle A\rangle\left(T\right) = \sum_{n = 0}^\infty\int_0^\infty \textrm{d} t\int \textrm{d}{\textbf{p}}\int \textrm{d}{\textbf{x}} \textrm{e}^{-t\gamma}f_{{\textrm{w}}_0}\left({\textbf{p}},{\textbf{x}}\right)g_n\left({\textbf{p}}^{^{\prime}}\left(t\right),{\textbf{x}}^{^{\prime}}\left(t\right),t\right) = \sum_{n = 0}^\infty\langle A\rangle_n\left(T\right). \end{align} \tag{ 19 }$$

The integration over the other set of variables $\textbf{p},\textbf{x},t$ gives rise to a transition to a forward parametrization of the arguments in the δ functions in the kernel:

$$\begin{align} &\delta\left({\textbf{p}}\left(t^{^{\prime}}\right)+{\textbf{i}}\Delta P-{\textbf{p}}^{^{\prime}},{\textbf{x}}\left(t^{^{\prime}}\right)+{\textbf{j}}\Delta X-{\textbf{x}}^{^{\prime}}\right)\nonumber\\ &\quad= \delta\left({\textbf{p}}-{\textbf{p}}^{^{\prime}}\left(t;{\textbf{p}}^{^{\prime}}-{\textbf{i}}\Delta P,{\textbf{x}}^{^{\prime}}-{\textbf{j}}\Delta X,t^{^{\prime}}\right),{\textbf{x}}-{\textbf{x}}^{^{\prime}}\left(t;{\textbf{p}}^{^{\prime}}-{\textbf{i}}\Delta P,{\textbf{x}}^{^{\prime}}-{\textbf{j}}\Delta X,t^{^{\prime}}\right)\right) . \end{align} \tag{ 20 }$$

For $\mathcal{K}$ this yields

$$\begin{align} \mathcal{K}\left({\textbf{p}},{\textbf{x}},t,{\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right)& = \theta\left(t-t^{^{\prime}}\right)\gamma\textrm{e}^{-\left(t-t^{^{\prime}}\right)\gamma}\sum_{{\textbf{i}},{\textbf{j}}}\alpha_{{\textbf{i}}{\textbf{j}}}\nonumber\\ &\quad \cdot\delta\left({\textbf{p}}-{\textbf{p}}^{^{\prime}}\left(t;{\textbf{p}}^{^{\prime}}-{\textbf{i}}\Delta P,{\textbf{x}}^{^{\prime}}-{\textbf{j}}\Delta X,t^{^{\prime}}\right),{\textbf{x}}-{\textbf{x}}^{^{\prime}}\left(t;{\textbf{p}}^{^{\prime}}-{\textbf{i}}\Delta P,{\textbf{x}}^{^{\prime}}-{\textbf{j}}\Delta X,t^{^{\prime}}\right)\right). \end{align} \tag{ 21 }$$

2.4.2. Solution for the adjoint equation.

We introduce the trajectory with scattering events for forward parameterization. Given the scattering indices $(\textbf{i}_k,\textbf{j}_k)_{1\unicode{x2A7D} k\unicode{x2A7D} n}$ and the scattering times $t_1\lt t_2\lt\ldots\lt t_n$, we use the convention $\textbf{p}_0^{^{\prime}}(t): = \textbf{p}(t;\textbf{p},\textbf{x},0)$, $\textbf{x}_0^{^{\prime}}(t): = \textbf{x}(t;\textbf{p},\textbf{x},0)$ to define

$$\begin{align} {\textbf{p}}_n^{^{\prime}}\left(t\right): = \left\{\begin{array}{ll}{\textbf{p}}_{n-1}^{^{\prime}}\left(t\right) &~\textrm{for}~0\unicode{x2A7D} t \unicode{x2A7D} t_n\\ {\textbf{p}}^{^{\prime}}\left(t;{\textbf{p}}^{^{\prime}}_{n-1}\left(t_{n}\right)-{\textbf{i}}_n\Delta P,{\textbf{x}}^{^{\prime}}_{n-1}\left(t_{n}\right)-{\textbf{j}}_n\Delta X,t_n\right)&~\textrm{for}~t_n \lt t\unicode{x2A7D} T,\\ \end{array}\right.\nonumber\\ {\textbf{x}}_n^{^{\prime}}\left(t\right): = \left\{\begin{array}{ll}{\textbf{x}}_{n-1}^{^{\prime}}\left(t\right) &~\textrm{for}~0\unicode{x2A7D} t \unicode{x2A7D} t_n\\ {\textbf{x}}^{^{\prime}}\left(t;{\textbf{p}}^{^{\prime}}_{n-1}\left(t_{n}\right)-{\textbf{i}}_n\Delta P,{\textbf{x}}^{^{\prime}}_{n-1}\left(t_{n}\right)-{\textbf{j}}_n\Delta X,t_n\right)&~\textrm{for}~t_n \lt t\unicode{x2A7D} T.\\ \end{array}\right. \end{align} \tag{ 22 }$$

A depiction of these functions can be seen in figure 2. The resolvent series for the solution is then presented by the term

$$\begin{align} g_n\left({\textbf{p}}^{^{\prime}}\left(t_1\right),{\textbf{x}}^{^{\prime}}\left(t_1\right),t_1\right) = ~\gamma^n\textrm{e}^{-\left(T-t_1\right)\gamma}\int\limits_{t_1}^T \textrm{d} t_2\ldots\int\limits_{t_{n-1}}^T \textrm{d} t_n \sum_{{\textbf{i}}_1,{\textbf{j}}_1\ldots{\textbf{i}}_n,{\textbf{j}}_n } \prod_{k = 1}^n\left(\alpha_{{\textbf{i}}_k{\textbf{j}}_k}\right)A_T\left({\textbf{p}}_n\left(T\right),{\textbf{x}}_n\left(T\right),T\right), \end{align} \tag{ 23 }$$

with $g_0(\textbf{p}^{^{\prime}}(t),\textbf{x}(t),t) = A_T(\textbf{p}_0(t),\textbf{x}_0(t),t)$.

Zoom In Zoom Out Reset image size

2.4.3. Iterative representation of physical quantities.

The series for the expectation value of a physical quantity is obtained by inserting (23) into (19). The general term is then

$$\begin{align} \langle A\rangle_n\left(T\right) = \gamma^n\textrm{e}^{-T\gamma}\! \int\! \textrm{d}{\textbf{p}}\!\int \!\textrm{d}{\textbf{x}} f_{{\textrm{w}}_0}\left({\textbf{p}},{\textbf{x}}\right)\! \int\limits_0^T\!\textrm{d} t_1\ldots\!\int\limits_{t_{n-1}}^T\!\textrm{d} t_n \sum_{{\textbf{i}}_1,{\textbf{j}}_1\ldots{\textbf{i}}_n,{\textbf{j}}_n}\prod_{k = 1}^n\left(\alpha_{{\textbf{i}}_k{\textbf{j}}_k}\right)A\left({\textbf{p}}_n^{^{\prime}}\left(T\right),{\textbf{x}}_n^{^{\prime}}\left(T\right),T\right). \end{align} \tag{ 24 }$$

Since both (24) and (16) are transformations of the general solution (12), they are indeed equivalent. Equation (24) remarkably resembles the corresponding expression for the Monte Carlo averages of an ensemble of M classical (Boltzmann) electrons, which move under the action of the Lorentz force and are scattered by, e.g. lattice vibrations (phonons) [39]. They are point-like particles with an initial distribution $f_{\textrm{w}_0}$, which initializes the starting phase space points $\textbf{p},\textbf{x}$. They determine Newtonian trajectories followed by the force particles during their free flight. The free flight is interrupted by scattering events, which, at a time t₁, update the phase-space coordinates. The latter initialize a novel piece of Newtonian trajectory for the next free flight. The evolution continues until the time T is reached and then each particle l contributes with its current value A_l (e.g. velocity, energy) to the statistical sum $\sum_l^MA_l$, which evaluates $\langle A\rangle$. The process corresponds to the scheme depicted in figure 2, which suggests a picture where pointlike quantum particles follow the same sequence of events. However, several problems need to be addressed to associate (24) with a quantum particle model. The classical initial distribution is non-negative, $f_{\textrm{w}_0}\unicode{x2A7E} 0$, while in the quantum case, $f_{\textrm{w}_0}$ could be any legitimate Wigner function and thus allows for negative values. This affects the evaluation of the physical averages, as can be seen already from the zeroth order term, which dominates if the evolution time is much smaller than the mean scattering time: In order to account for the sign, the statistical sum for the envisaged quantum particle model must be generalized to $\sum_l^Mw_lA_l$ where the quantity w_l , called weight, should carry the sign of $f_{\textrm{w}_0}$ in the point of initialization of the lth particle. Next, in the classical evolution, the scattering time (e.g. t₁) exponentially depends on the frequency of interaction with phonons, while in the quantum counterpart the sequence $t_1\lt t_2\lt\cdots$ is predetermined. This suggests looking for an analogical physical interpretation of the prefactor in (24). Finally, both classical and quantum counterparts rely on Newtonian trajectories, and hence the difference between the two kinds of evolution is due to the scattering: A fundamental difference between classical and quantum scattering is expected. These problems, formulated by heuristic considerations, are rigorously addressed next by the rules of the Monte Carlo theory for integration.

The two algorithms presented in this section differ in both, their parameterization and the distribution of the scattering times. The first one is more formal and evaluates $f_\textrm{w}$ pointwise using backward parametrization and a uniform distribution for the scattering times. The other one uses the more transparent (from a physical point of view) forward parametrization and introduces an exponential distribution of the scattering, which is a characteristic of the evolution of classical particles in the presence of scattering events. This gives rise to a quantum particle model, where the evolution of pointlike particles consists of consecutive events of free-flight over the Lorentz force-governed Newtonian trajectories, followed by scattering events.

3.1. Backward algorithm

The Monte Carlo algorithm introduced in this section allows to evaluate the terms $\langle A\rangle_n(T)$ of the resolvent expansion in (16). For this purpose, the integrals and sums are expressed as an expectation value $E[X_n]$ with a probability density $P_{X_n}$ and a random variable X_n . The terms $\langle A\rangle_n(T)$ are set to

$$\begin{align} \langle A\rangle_n\left(T\right) = E\left[X_n\right] = &\int \textrm{d}{\textbf{p}}\int \textrm{d}{\textbf{x}} \int_0^T \textrm{d} t_1\int_0^{t_1} \textrm{d} t_2\ldots\int_0^{t_{n-1}} \textrm{d} t_n \sum_{{\textbf{i}}_1,{\textbf{j}}_1}\sum_{{\textbf{i}}_2,{\textbf{j}}_2}\ldots\sum_{{\textbf{i}}_n,{\textbf{j}}_n}P_{X_n}X_n. \end{align} \tag{ 25 }$$

$P_{X_n}$ acts as a selector for the scattering indices $(\textbf{i}_k,\textbf{j}_k)_{1\unicode{x2A7D} k\unicode{x2A7D} n}$, the scattering times $t_1\lt t_2\lt\ldots\lt t_n$, and the initial points $\textbf{p},\textbf{x}$ of the trajectory. Thus, it is split into a product of three probability functions:

• For the coefficients $\alpha_{{\textbf{i}}{\textbf{j}}}$ of the kernel (11), we introduce a discrete transition probability $P_{{\textbf{i}}{\textbf{j}}}: = \frac{|\alpha_{{\textbf{i}}{\textbf{j}}}|}{|\alpha|}$, where $|\alpha|: = \sum_{{\textbf{i}}{\textbf{j}}}|\alpha_{{\textbf{i}}{\textbf{j}}}| = 41$, see (A.1). This means that the direction in which the trajectory scatters is chosen randomly, distributed proportionally to $|\alpha_{{\textbf{i}}{\textbf{j}}}|$.
• For the initial points $\textbf{p},\textbf{x}$ of the trajectory, a density function P is chosen. Both A and $f_{\textrm{w}_0}$ depend on $\textbf{p}$ and $\textbf{x}$, thus a possible choice could be $P(\textbf{p},\textbf{x})\propto |A(\textbf{p},\textbf{x})f_{\textrm{w}_0}(\textbf{p},\textbf{x})|$.
• The scattering times $t_1,\ldots,t_n$ are evenly distributed on the intervals $(0,T)$ for t₁ and on $(0,t_{i-1})$ for $t_i,i\in\{2,\ldots,n\}$. The density function of a uniform distribution is normalized by the inverse of the length of the integral, which has to be considered in X_n by the product $T\prod_{i = 1}^{n-1}t_i$.

In combination they yield $P_{X_n} = |\alpha_{{\textbf{i}}{\textbf{j}}}|/|\alpha|P(\textbf{p},\textbf{x})(T\prod_{i = 1}^{n-1}t_i)^{-1}$. Since the corresponding random variable X_n is the estimator of $\langle A\rangle_n(T)$, it is evaluated and averaged for several arguments randomly selected according to $P_{X_n}$. To satisfy (25) it is given as

$$\begin{align} X_n = \gamma^n|\alpha|^n\textrm{e}^{-T\gamma}\frac{A\left({\textbf{p}},{\textbf{x}},T\right)}{P\left({\textbf{p}},{\textbf{x}}\right)} T\prod_{i = 1}^{n-1}\left(t_i\right)\prod_{k = 1}^n\left({\textrm{sign}}\left(\alpha_{{\textbf{i}}_k{\textbf{j}}_k}\right)\right)f_{{\textrm{w}}_0}\left({\textbf{p}}_n\left(0\right),{\textbf{x}}_n\left(0\right)\right). \end{align} \tag{ 26 }$$

The expectation value of a physical quantity can be obtained by algorithm 3.1 (see also figure 3).

Zoom In Zoom Out Reset image size

Algorithm 3.1. Backward algorithm.
1 Initialization of N, $(N_n)_{n\in\{0,1,2,\ldots,N\}},n\leftarrow0$ and a variable, say $(A_n)_{n\in\{0,1,2,\ldots,N\}}\leftarrow\vec{0}$. N sets the total number of terms in the iterative expansion (13). Nn determines the number of independent numerical trajectories with n scattering events. $j\leftarrow 1$ is a counter for Nn . An represents the value of the nth term in the resolvent expansion. t0 is initialized by $t_0\leftarrow T$.
2 If $n\neq 0,$ the scattering times $(t_i)_{i\in\{1,2,\ldots,n\}}$ are chosen in order because the upper limit of every ti depends on $t_{i-1}$. Each $t_i\sim \textrm{U}(0,t_{i-1})$ is generated randomly, with the uniform distribution $\textrm{U}$ on the interval $(0,t_{i-1})$. $s\leftarrow1$, which represents all factors in An that are updated at each scattering event, and $i\leftarrow0$. $(\textbf{p},\textbf{x})\sim P(\textbf{p},\textbf{x})$ are chosen randomly, and distributed according to the chosen probability function $P(\textbf{p},\textbf{x})$. The initial values $\textbf{p}_T\leftarrow\textbf{p}$ and $\textbf{x}_T\leftarrow\textbf{x}$ are stored separately. If n = 0, jump to step 5.
3 Starting from the current $\textbf{p}$ and $\textbf{x}$ the trajectory is followed until it reaches the next scattering event at $t_{i+1}$, i.e. $\textbf{p}\leftarrow\textbf{p}(t_{i+1};\textbf{p},\textbf{x},t_i)$ and $\textbf{x}\leftarrow\textbf{x}(t_{i+1};\textbf{p},\textbf{x},t_i)$, and then $i\leftarrow i+1$.
4 In the event of scattering: Values for $(\textbf{i},\textbf{j})\sim P_{{\textbf{i}}{\textbf{j}}}$ are chosen randomly, distributed according to the values of the transition probability $P_{{\textbf{i}}{\textbf{j}}}$. Then s is updated to $s\leftarrow s\cdot \gamma t_{i-1}|\alpha| \textrm{sign}(\alpha_{{\textbf{i}}{\textbf{j}}})$. The factor $t_{i-1}$ comes from the length of the time integral. Finally, $\textbf{p}\leftarrow\textbf{p}+\textbf{i}\Delta P,\textbf{x}\leftarrow\textbf{x}+\textbf{j}\Delta X$. If i < n, jump to step 3.
5 The trajectory is followed backward in the time interval $(0,t_n)$, i.e. $\textbf{p}\leftarrow\textbf{p}(0;\textbf{p},\textbf{x},t_n)$ and $\textbf{x}\leftarrow\textbf{x}(0;\textbf{p},\textbf{x},t_n)$, where $(\textbf{p},\textbf{x})$ is equal to the phase space point $(\textbf{p}_n(0),\textbf{x}_n(0))$, see figure 1.
6   $f_{\textrm{w}_0}(\textbf{p},\textbf{x})$ is evaluated at the final position $(\textbf{p},\textbf{x}) = (\textbf{p}_n(0),\textbf{x}_n(0))$ and $A_n\leftarrow A_n+s\textrm{e}^{-T\gamma}f_{\textrm{w}_0}(\textbf{p},\textbf{x})A(\textbf{p}_T,\textbf{x}_T,T)/P(\textbf{p}_T,\textbf{x}_T)$. If $j\lt N_n$, set $j\leftarrow j+1$ and jump to step 2.
7   $n\leftarrow n+1,j\leftarrow 1$, and the algorithm jumps to step 2, unless n = N. In this case, the next step is executed.
8 Finally, return $\sum_{n = 0}^N A_n/N_n$.

3.2. Forward algorithm

Finally, a Monte Carlo algorithm is presented, where the number of scattering events is not predetermined and the scattering times are exponentially distributed. We will use forward parameterization in this case. Again, X_n and $P_{X_n}$ have to satisfy the condition $\langle A\rangle_n(T) = E[X_n]$. The arguments that are randomly chosen are the same as before. The transition probability $P_{{\textbf{i}}{\textbf{j}}}$ and the density function P remain the same. For the scattering times $t_1,\ldots,t_n$, we evaluate the joint density of the number of scattering events n happening in the interval $[0,T]$, and the consecutive scattering times $(t_i)_{i\in\{1,\ldots,n\}}$. Considering an exponential distribution, the density for a single scattering event is given by $\gamma\textrm{e}^{-\gamma t}$. The joint density is equal to the density of the first n events multiplied by the probability that the next event happens after T, which yields

$$\begin{align} p\left(\left(t_i\right)_{i\in\left\{1,\ldots,n\right\}},n\right)& = \gamma^n\prod_{i = 1}^n\left(\textrm{e}^{-\gamma\left(t_i-t_{i-1}\right)}\right)\int_T^\infty\gamma\textrm{e}^{-\gamma\left(t_{n+1}-t_n\right)}\textrm{d} t_{n+1}\nonumber\\ & = \gamma^n\textrm{e}^{-\gamma t_n}\textrm{e}^{\gamma t_n}\int_T^\infty\gamma\textrm{e}^{-\gamma t_{n+1}}\textrm{d} t_{n+1} \nonumber\\ & = \gamma^n\textrm{e}^{-\gamma T}, \end{align} \tag{ 27 }$$

assuming t₀ = 0. This conveniently coincides with the prefactor in (24).

Combining all probability functions gives $P_{X_n} = \prod_{k = 1}^n(|\alpha_{{\textbf{i}}_k{\textbf{j}}_k}|)|\alpha|^{-n}P(\textbf{p},\textbf{x})\gamma^n\textrm{e}^{-\gamma T}$. By using the condition $\langle A\rangle_n(T) = E[X_n]$ and the result of $\langle A\rangle_n(T)$ in (24), we can evaluate the random variable as

$$\begin{equation} X_n = |\alpha|^n\frac{f_{{\textrm{w}}_0}\left({\textbf{p}},{\textbf{x}}\right)}{P\left({\textbf{p}},{\textbf{x}}\right)}\prod_{k = 1}^n\left({\textrm{sign}}\left(\alpha_{{\textbf{i}}_k{\textbf{j}}_k}\right)\right)A\left({\textbf{p}}_n^{^{\prime}}\left(T\right),{\textbf{x}}_n^{^{\prime}}\left(T\right),T\right). \end{equation} \tag{ 28 }$$

The expectation value of a physical quantity can be obtained by algorithm 3.2 (see also figure 4).

Zoom In Zoom Out Reset image size

Algorithm 3.2. Forward algorithm.
1 Initialization of M and a variable, say $A\leftarrow0$. M sets the total number of the independent numerical trajectories and A represents the expectation value of the physical quantity. $j\leftarrow 1$ is a counter for M.
2   ti is initialized as $t_i\leftarrow 0$. $s\leftarrow1$ represents all factors in Xn that are updated at each scattering event. $(\textbf{p},\textbf{x})\sim P(\textbf{p},\textbf{x})$ are chosen randomly, distributed according to the chosen probability function $P(\textbf{p},\textbf{x})$. Since $(\textbf{p},\textbf{x})$ will change in the following steps, the initial values $\textbf{p}_0\leftarrow\textbf{p},\textbf{x}_0\leftarrow\textbf{x}$ are also saved as they are needed at a later step.
3 An exponentially distributed variable $t^{^{\prime}}\sim \textrm{Exp}(\gamma)$ with the constant γ is chosen by generating a uniformly distributed variable $r\sim \textrm{U}(0,1)$ and setting $t^{^{\prime}}\leftarrow -\text{ln}(r)/\gamma$. If $t_i+t^{^{\prime}}\gt T$, then we jump to step 6.
4 Starting from the current $\textbf{p}$ and $\textbf{x}$ the trajectory is followed until it reaches the next scattering event at $t_i+t^{^{\prime}}$, i.e. $\textbf{p}\leftarrow\textbf{p}^{^{\prime}}(t_i+t^{^{\prime}};\textbf{p},\textbf{x},t_i)$ and $\textbf{x}\leftarrow\textbf{x}^{^{\prime}}(t_i+t^{^{\prime}};\textbf{p},\textbf{x},t_i)$.
5 In the event of scattering: Values for $(\textbf{i},\textbf{j})\sim P_{{\textbf{i}}{\textbf{j}}}$ are chosen randomly, distributed according to the values of the transition probability $P_{{\textbf{i}}{\textbf{j}}}$, defined in section 3.1. Then, s is updated to $s\leftarrow s\cdot |\alpha|\textrm{sign}(\alpha_{{\textbf{i}}{\textbf{j}}})$. Finally, $\textbf{p}\leftarrow\textbf{p}-\textbf{i}\Delta P,\textbf{x}\leftarrow\textbf{x}-\textbf{j}\Delta X$ and $t_i\leftarrow t_i+t^{^{\prime}}$. Then jump to step 3.
6 The trajectory is followed in the time interval $(t_i,T)$, i.e. $\textbf{p}\leftarrow\textbf{p}^{^{\prime}}(T;\textbf{p},\textbf{x},t_i)$ and $\textbf{x}\leftarrow\textbf{x}^{^{\prime}}(T;\textbf{p},\textbf{x},t_i)$, where $(\textbf{p},\textbf{x})$ is equal to the phase space point $(\textbf{p}_n^{^{\prime}}(T),\textbf{x}_n^{^{\prime}}(T))$, see figure 2.
7   $A(\textbf{p},\textbf{x},T)$ is evaluated at the final position $(\textbf{p},\textbf{x}) = (\textbf{p}_n^{^{\prime}}(T),\textbf{x}_n^{^{\prime}}(T))$ and $A\leftarrow A+sA(\textbf{p},\textbf{x},T)f_{\textrm{w}_0}(\textbf{p}_0,\textbf{x}_0)/P(\textbf{p}_0,\textbf{x}_0)$. $j\leftarrow j+1$.
8 Jump to step 2, unless j = M. In this case, the next step is executed.
9 Finally, return $A/M$.

The two introduced Monte Carlo algorithms constitute an important step in understanding gauge-invariant Wigner theory using classical Boltzmann concepts. The choice of linear electromagnetic fields ensures the appearance of the same Liouville operator in both transport descriptions and thus provides a convenient reference frame for insights into the quantum evolution in terms of particles. The two algorithms are derived by the application of established Monte Carlo approaches for integrating the backward or forward form of the gauge-invariant Wigner equation. In the former case, the algorithm is more formal as the evolution proceeds backward in time. It offers computational advantages when the solution is needed locally in the phase space. Furthermore, it allows us to gradually introduce concepts used in the forward algorithm, which completes the particle picture conjectured at the end of the previous section. The quantum evolution resembles to a large extent the evolution of classical Boltzmann particles. An ensemble of particles is initialized in both cases according to the initial condition. Particles are accelerated by the Lorentz force over Newtonian trajectories and interrupted by scattering events. Comparing both algorithms reveals the proper interpretation of the distribution of the scattering times. In the backward algorithm, the scattering times were chosen uniformly distributed on the interval between the beginning of the evaluation and the previous scattering event. As a result, the scattering events tend to be unevenly distributed throughout the evolution time. The distribution density of the scattering events is inversely proportional to the length of the time intervals $(0,t_i)_{i\in\{1,\ldots,n-1\}}$, and is thus higher at 0 and lower toward T. In the forward algorithm, the scattering events are evenly distributed on the interval $[0,T]$, due to the exponential distribution. This manifests in the joint probability density, which corresponds exactly to the prefactor of the terms in the resolvent expansion. As for the weights of the statistical sum of the physical quantity, their absolute value is multiplied by $|\alpha|$ for every scattering event. This factor corresponds to the weighted amount of possible directions the particle could scatter. Also, the sign of the weights can change during the scattering, depending on the sign of the corresponding coefficient $\alpha_{{\textbf{i}}{\textbf{j}}}$ in the kernel.

These considerations can be summarized as follows: The distribution of scattering times is given by the formally introduced quantity γ, (10), which now has been provided with a physical meaning of a total out-scattering rate in a striking analogy with the classical counterpart. Similarly to the latter, γ is given by the sum of the quantities $|\alpha_{{\textbf{i}}{\textbf{j}}}|$, which corresponds to the probability for scattering from different classical mechanisms such as phonons and impurities. The difference is that the terms $\alpha_{{\textbf{i}}{\textbf{j}}}$ carry a sign, so that each scattering event can change both the absolute value of the weight and the sign, which are the main attributes of a quantum particle. Indeed, in this way scattering determines the difference between classical and quantum evolution, as discussed before. Furthermore, while in the former case, scattering is local in space, causing only a shift in momentum, quantum scattering leads to spatial shifts. These shifts depend on the finite difference scheme, however, this is irrelevant to the conceptual understanding: Similarly, considering computational approaches, different numerical schemes can be applied to find the numerical solution.

The introduction of the Newtonian trajectory enables us to transform the gauge-invariant Wigner equation to a Fredholm integral equation, where a resolvent expansion gives an iterative solution. However, this involves the approximation of the high-order derivative term leading to many terms in the kernel. This consequently increases the number of possible paths of the trajectory giving rise to the accumulation of the weight of a trajectory with the evolution. Large positive and negative weight values need to cancel each other in the statistical estimators for the physical averages. Thus, the maximum simulation time T of the simulation is limited, because the larger T, the higher the impact of the terms with a higher number of scattering events. From a computational point of view, this leads to the well-known 'sign problem' of quantum mechanics. A good example is the Taylor series of e^−x for large positive x, where large terms compensate each other to give a value smaller than unity. The problem can be addressed by using the Markovian character of the evolution of the particle ensemble, which, in particular, provides the Wigner solution f_w in the entire phase space: T can be decomposed on shorter time intervals Δt, so that the solution at the end of the nth interval $f_w(n\Delta T)$ becomes the initial condition for the n + 1th interval.

This research was funded in whole by the Austrian Science Fund (FWF): P33609-N and P37080-N. For open access purposes, the author has applied a CC BY public copyright license to any author accepted manuscript version arising from this submission.

No new data were created or analysed in this study.

The detailed form of the kernel (11) is given by

$$\begin{equation} \begin{aligned} &\mathcal{K}\left({\textbf{p}},{\textbf{x}},t,{\textbf{p}}^{^{\prime}},{\textbf{x}}^{^{\prime}},t^{^{\prime}}\right) = \theta\left(t-t^{^{\prime}}\right)\gamma\textrm{e}^{-\left(t-t^{^{\prime}}\right)\gamma}\nonumber\\ &\quad\cdot\left[ 4\delta\left(p_x\left(t^{^{\prime}}\right)-p_x^{^{\prime}}, p_y\left(t^{^{\prime}}\right) + \Delta P-p_y^{^{\prime}}, x\left(t^{^{\prime}}\right) + \Delta X-x^{^{\prime}}, y\left(t^{^{\prime}}\right)-y^{^{\prime}}\right)\right.\nonumber \\ &\quad- 8 \delta\left(p_x\left(t^{^{\prime}}\right)-p_x^{^{\prime}}, p_y\left(t^{^{\prime}}\right)-p_y^{^{\prime}}, x\left(t^{^{\prime}}\right) + \Delta X-x^{^{\prime}}, y\left(t^{^{\prime}}\right)-y^{^{\prime}}\right)\nonumber \\ &\quad+ 4\delta\left(p_x\left(t^{^{\prime}}\right)-p_x^{^{\prime}}, p_y\left(t^{^{\prime}}\right) - \Delta P-p_y^{^{\prime}}, x\left(t^{^{\prime}}\right) + \Delta X-x^{^{\prime}}, y\left(t^{^{\prime}}\right)-y^{^{\prime}}\right) \nonumber \\ &\quad- 4\delta\left(p_x\left(t^{^{\prime}}\right)-p_x^{^{\prime}}, p_y\left(t^{^{\prime}}\right) + \Delta P-p_y^{^{\prime}}, x\left(t^{^{\prime}}\right) - \Delta X-x^{^{\prime}}, y\left(t^{^{\prime}}\right)-y^{^{\prime}}\right) \nonumber\\ &\quad+ 8 \delta\left(p_x(t^{^{\prime}}\right)-p_x^{^{\prime}}, p_y(t^{^{\prime}})-p_y^{^{\prime}}, x(t^{^{\prime}}) - \Delta X-x^{^{\prime}}, y(t^{^{\prime}})-y^{^{\prime}})\nonumber\\ &\quad- 4\delta(p_x(t^{^{\prime}})-p_x^{^{\prime}}, p_y(t^{^{\prime}}) - \Delta P-p_y^{^{\prime}}, x(t^{^{\prime}}) - \Delta X-x^{^{\prime}}, y(t^{^{\prime}})-y^{^{\prime}}) \nonumber\\ &\quad-\delta(p_x(t^{^{\prime}}) + \Delta P-p_x^{^{\prime}}, p_y(t^{^{\prime}}) + \Delta P-p_y^{^{\prime}}, x(t^{^{\prime}})-x^{^{\prime}}, y(t^{^{\prime}}) + \Delta X-y^{^{\prime}})\nonumber \\ &\quad+ \delta(p_x(t^{^{\prime}}) + \Delta P-p_x^{^{\prime}}, p_y(t^{^{\prime}}) - \Delta P-p_y^{^{\prime}}, x(t^{^{\prime}})-x^{^{\prime}}, y(t^{^{\prime}}) + \Delta X-y^{^{\prime}}) \nonumber\\ &\quad+\delta(p_x(t^{^{\prime}}) - \Delta P-p_x^{^{\prime}}, p_y(t^{^{\prime}}) + \Delta P-p_y^{^{\prime}}, x(t^{^{\prime}})-x^{^{\prime}}, y(t^{^{\prime}}) + \Delta X-y^{^{\prime}}) \nonumber\\ &\quad- \delta(p_x(t^{^{\prime}}) - \Delta P-p_x^{^{\prime}}, p_y(t^{^{\prime}}) - \Delta P-p_y^{^{\prime}}, x(t^{^{\prime}})-x^{^{\prime}}, y(t^{^{\prime}}) + \Delta X-y^{^{\prime}}) \nonumber\\ &\quad+ \delta(p_x(t^{^{\prime}}) + \Delta P-p_x^{^{\prime}}, p_y(t^{^{\prime}}) + \Delta P-p_y^{^{\prime}}, x(t^{^{\prime}})-x^{^{\prime}}, y(t^{^{\prime}}) - \Delta X-y^{^{\prime}}) \nonumber\\ &\quad- \delta(p_x(t^{^{\prime}}) + \Delta P-p_x^{^{\prime}}, p_y(t^{^{\prime}}) - \Delta P-p_y^{^{\prime}}, x(t^{^{\prime}})-x^{^{\prime}}, y(t^{^{\prime}}) - \Delta X-y^{^{\prime}}) \nonumber\\ &\quad- \delta(p_x(t^{^{\prime}}) - \Delta P-p_x^{^{\prime}}, p_y(t^{^{\prime}}) + \Delta P-p_y^{^{\prime}}, x(t^{^{\prime}})-x^{^{\prime}}, y(t^{^{\prime}}) - \Delta X-y^{^{\prime}})\nonumber \\ &\quad+ \delta(p_x(t^{^{\prime}}) - \Delta P-p_x^{^{\prime}}, p_y(t^{^{\prime}}) - \Delta P-p_y^{^{\prime}}, x(t^{^{\prime}})-x^{^{\prime}}, y(t^{^{\prime}}) - \Delta X-y^{^{\prime}})\nonumber\\ &\quad\left. +\delta\bigl({\textbf{p}}(t^{^{\prime}})-{\textbf{p}}^{^{\prime}},{\textbf{x}}(t^{^{\prime}})-{\textbf{x}}^{^{\prime}}\bigr) \right]. \end{aligned} \end{equation} \tag{ A.1 }$$

B.1. Fredholm integral equation of the second kind

The solution for a general Fredholm equation of the second kind, $f(s) = f_\textrm{i}(s)+\int_a^b\mathcal{K}(s,s^{^{\prime}})f(s^{^{\prime}})ds^{^{\prime}}$ is given by

$$\begin{align} f\left(s\right) &= \sum_{n = 0}^\infty f_n\left(s\right), \nonumber\\ f_0\left(s\right): &=\, f_{\textrm{i}}\left(s\right),\nonumber \\ f_n\left(s\right): &= \int_a^b\ldots\int_a^b\int_a^b\mathcal{K}\left(s,t_1\right)\mathcal{K}\left(t_1,t_2\right)\cdots\mathcal{K}\left(t_{n-1},t_n\right)f_{\textrm{i}}\left(t_n\right)\textrm{d} t_1\textrm{d} t_2\ldots \textrm{d} t_n, \end{align} \tag{ B.1 }$$

provided that the series converges. To derive the solution with forward parameterization, the adjoint integral equation of the Fredholm equation is used. In general, given a Fredholm equation problem with $\mathcal{K}:[a,b]\times[a,b]\rightarrow\mathbb{R}, f_\textrm{i}:[a,b]\rightarrow\mathbb{R}, a,b\in[-\infty,\infty], a\lt b$, like in (B.1), let $g_\textrm{i}:[a,b]\rightarrow\mathbb{R}$. Then, the adjoint equation is defined as

$$\begin{align} g\left(s^{^{\prime}}\right) = g_{\textrm{i}}\left(s^{^{\prime}}\right)+\int_a^b\mathcal{K}\left(s,s^{^{\prime}}\right)g\left(s\right)\textrm{d} s, \end{align} \tag{ B.2 }$$

with a not specified initial function $g_\textrm{i}$. The solution for g is also given by a resolvent expansion, but the order of the variables in $\mathcal{K}$ is reversed, i.e.

$$\begin{align} g\left(s\right)& = \sum_{n = 0}^\infty g_n\left(s\right),\nonumber\\ g_0\left(s\right)&: = g_{\textrm{i}}\left(s\right),\nonumber\\ g_n\left(s\right)&: = \int_a^b\ldots\int_a^b\int_a^b\mathcal{K}\left(t_1,s\right)\mathcal{K}\left(t_2,t_1\right)\cdots\mathcal{K}\left(t_{n},t_{n-1}\right)g_{\textrm{i}}\left(t_n\right)\textrm{d} t_1\textrm{d} t_2\ldots \textrm{d} t_n. \end{align} \tag{ B.3 }$$

B.2. Exchange Lemma

Let $f:[a,b]\rightarrow\mathbb{R}$ be the solution of a Fredholm equation, where $\mathcal{K}:[a,b]\times[a,b]\rightarrow\mathbb{R}, f_\textrm{i}:[a,b]\rightarrow\mathbb{R}$, and let $g:[a,b]\rightarrow\mathbb{R}$ be the solution of the adjoint equation, where $g_\textrm{i}:[a,b]\rightarrow\mathbb{R}$. Then, there holds

$$\begin{align} \int_a^bf_{\textrm{i}}\left(s\right)g\left(s\right)\textrm{d} s = \int_a^bf\left(s\right)g_i\left(s\right)\textrm{d} s. \end{align} \tag{ B.4 }$$

Proof. 

$$\begin{align*} \int_a^bf_{\textrm{i}}\left(s\right)g\left(s\right)\textrm{d} s& = \int_a^bg\left(s\right)\left[f\left(s\right)-\int_a^b\mathcal{K}\left(s,s^{^{\prime}}\right)f\left(s^{^{\prime}}\right)ds^{^{\prime}}\right]\textrm{d} s\\ & = \int_a^bg\left(s\right)f\left(s\right)\textrm{d} s-\int_a^b\int_a^b\mathcal{K}\left(s,s^{^{\prime}}\right)f\left(s^{^{\prime}}\right)g\left(s\right)\textrm{d} s^{^{\prime}}\textrm{d} s\\ & = \int_a^bf\left(s^{^{\prime}}\right)\left[g\left(s^{^{\prime}}\right)-\int_a^b\mathcal{K}\left(s,s^{^{\prime}}\right)g\left(s\right)\textrm{d} s\right]\textrm{d} s^{^{\prime}}\\ & = \int_a^bf\left(s^{^{\prime}}\right)g_{\textrm{i}}\left(s^{^{\prime}}\right)\textrm{d} s^{^{\prime}}. \end{align*}$$

 □

This fact can be used to express the solution f by the adjoint solution g at a given point $s\in[a,b]$. In particular, if we set $g_\textrm{i}(s^{^{\prime}}): = \delta(s-s^{^{\prime}})$, we can show that

$$\begin{align} f\left(s\right) = \int_a^bf\left(s^{^{\prime}}\right)\delta\left(s-s^{^{\prime}}\right)\textrm{d} s^{^{\prime}} = \int_a^bf\left(s^{^{\prime}}\right)g_{\textrm{i}}\left(s^{^{\prime}}\right)\textrm{d} s^{^{\prime}} = \int_a^bf_{\textrm{i}}\left(s^{^{\prime}}\right)g\left(s^{^{\prime}}\right)\textrm{d} s^{^{\prime}}. \end{align} \tag{ B.5 }$$

The solution of (11) is given by the series $\sum_{n = 0}^\infty f_n(\textbf{p},\textbf{x},t)$, where

$$\begin{align} & f_0\left({\textbf{p}},{\textbf{x}},t\right) = ~\textrm{e}^{-t\gamma}f_{{\textrm{w}}_0}\left({\textbf{p}}_0\left(0\right),{\textbf{x}}_0\left(0\right)\right)\nonumber\\ &f_n\left({\textbf{p}},{\textbf{x}},t\right) = ~\gamma^n\textrm{e}^{-\gamma t}\int_0^tdt_1\int_0^{t_1}dt_2\ldots\int_0^{t_{n-1}}dt_n \nonumber\\ & \quad \times \sum_{{\textbf{i}}_1,{\textbf{j}}_1}\sum_{{\textbf{i}}_2,{\textbf{j}}_2}\ldots\sum_{{\textbf{i}}_n,{\textbf{j}}_n}\prod_{k = 1}^n\left(\alpha_{{\textbf{i}}_k{\textbf{j}}_k}\right)f_{{\textrm{w}}_0}\left({\textbf{p}}_n\left(0\right),{\textbf{x}}_n\left(0\right)\right). \end{align}$$

Proof. First, we show

$$\begin{align*} &\int d{\textbf{p}}_{1}\int d{\textbf{x}}_{1}\ldots\int d{\textbf{p}}_{n}\int d{\textbf{x}}_{n} \delta\left({\textbf{p}}\left(t_1\right)+{\textbf{i}}_1\Delta P-{\textbf{p}}_1,{\textbf{x}}\left(t_1\right)+{\textbf{j}}_1\Delta X-{\textbf{x}}_1\right)\\ &\quad\cdot\prod_{k = 1}^{n-1}\delta\left({\textbf{p}}\left(t_{k+1};{\textbf{p}}_k,{\textbf{x}}_k,t_k\right)+{\textbf{i}}_{k+1}\Delta P-{\textbf{p}}_{k+1},{\textbf{x}}\left(t_{k+1};{\textbf{p}}_k,{\textbf{x}}_k,t_k\right)+{\textbf{j}}_{k+1}\Delta X-{\textbf{x}}_{k+1}\right)\\ &\quad \cdot f_{{\textrm{w}}_0}\left({\textbf{p}}\left(0\right),{\textbf{x}}\left(0\right)\right) = f_{{\textrm{w}}_0}\left({\textbf{p}}_{n}\left(0\right),{\textbf{x}}_{n}\left(0\right)\right), \end{align*}$$

by induction, assuming $(\textbf{p},\textbf{x},t)$ as initial point and $t\gt t_1\gt t_2\gt\ldots\gt t_{n+1}$.

Base Case:

$$\begin{align*} f_{{\textrm{w}}_0}\left({\textbf{p}}\left(0\right),{\textbf{x}}\left(0\right)\right) = f_{{\textrm{w}}_0}\left({\textbf{p}}_0\left(0\right),{\textbf{x}}_0\left(0\right)\right). \end{align*}$$

Induction Step:

$$\begin{align*} &\int d{\textbf{p}}_{1}\int d{\textbf{x}}_{1}\ldots\int d{\textbf{p}}_{n+1}\int d{\textbf{x}}_{n+1} \delta\left({\textbf{p}}\left(t_1\right)+{\textbf{i}}_1\Delta P-{\textbf{p}}_1,{\textbf{x}}\left(t_1\right)+{\textbf{j}}_1\Delta X-{\textbf{x}}_1\right)\\ &\quad\cdot\prod_{k = 1}^{n}\delta\left({\textbf{p}}\left(t_{k+1};{\textbf{p}}_k,{\textbf{x}}_k,t_k\right)+{\textbf{i}}_{k+1}\Delta P-{\textbf{p}}_{k+1},{\textbf{x}}\left(t_{k+1};{\textbf{p}}_k,{\textbf{x}}_k,t_k\right)+{\textbf{j}}_{k+1}\Delta X-{\textbf{x}}_{k+1}\right)f_{{\textrm{w}}_0}\left({\textbf{p}}\left(0\right),{\textbf{x}}\left(0\right)\right)\\ & = \int d{\textbf{p}}_{n+1}\int d{\textbf{x}}_{n+1}\delta\left({\textbf{p}}\left(t_{n+1}\right)+{\textbf{i}}\Delta P-{\textbf{p}}_{n+1},{\textbf{x}}\left(t_{n+1}\right)+{\textbf{j}}\Delta X-{\textbf{p}}_{n+1}\right)f_{{\textrm{w}}_0}\left({\textbf{p}}_n\left(0\right),{\textbf{x}}_n\left(0\right)\right)\\ & = f_{{\textrm{w}}_0}\left({\textbf{p}}\left(0;{\textbf{p}}_n\left(t_{n+1}\right)+{\textbf{i}}\Delta P,{\textbf{x}}_n\left(t_{n+1}\right)+{\textbf{j}}\Delta X,t_{n+1}\right), {\textbf{x}}\left(0;{\textbf{p}}_n\left(t_{n+1}\right)+{\textbf{i}}\Delta P,{\textbf{x}}_n\left(t_{n+1}\right)+{\textbf{j}}\Delta X,t_{n+1}\right)\right)\\ & = f_{{\textrm{w}}_0}\left({\textbf{p}}_{n+1}\left(0\right),{\textbf{x}}_{n+1}(0\right)\big), \end{align*}$$

where we have used (14). Using this result in the resolvent expansion (B.1) yields the desired result. □