3.3 Integral formulation

The WBE, introduced in Chapter 2, is written here as

( ∂         )             ∫    {  (      )      (       )}   (      )
  -- + vg(k)  fw (r,k,t) =   dk ′ S r,k, k′ + Vw  r,k′ - k  fw r,k′,t +
  ∂t
                          -∫ λ (r,k)fw (r,k,t) ;                                                       (3.1)
                 λ(r,k) ≡   dk ′S (r,k′,k ).                                                           (3.2)
The first term on the RHS of (3.1) denotes the in-scattering and the second term denotes the out-scattering (at a rate λ). The Wigner potential will later emerge to be interpretable as a scattering mechanism. Although the scattering rates and the Wigner potential are assumed to be time-independent, introducing a time-dependence does not present a conceptual problem.

The Wigner potential can be decomposed as

 Vw (r,k ) = Vw+ (r,k )- Vw-(r,k);
  +
Vw  (r,k ) ≡ max (0,Vw );                                                           (3.3)
Vw- (r,k ) ≡ max (0, - Vw).
The quantity, which will emerge to be the scattering rate associated with the Wigner potential, is defined as
       ∫
γ (r) ≡   dkVw+ (r,k ) ≥ 0.                                                      (3.4)

The term γ(r)fw(r,k,t) is added to both sides of (3.1), such that

                  (                          )
                    ∂-+ vg (k )+ λ (r,k )+ γ (r)  fw (r,k, t) =
                    ∂t                                                                                    (3.5)
∫   ′{  (      ′)    + (   ′    )    -(    ′   )}   (   ′ )
  dk   S r,k,k   + Vw  r,k - k  - Vw  r,k  - k  fw  r,k ,t + γ (r)fw (r,k,t).
This is expressed in the more compact form
(                   )
  ∂                               ∫    ′ (     ′  )   (   ′ )
 ∂t + vg(k) + μ(r,k)  fw (r,k,t) =   dk Γ  r,k,k ,t fw r,k ,t ;                                    (3.6)

                      (  μ(r,k)) ≡ λ (r(,k )+ γ) (r);  (       )                                       (3.7)
                    Γ  r,k,k′,t ≡ S  r,k,k′ + Vw+  r,k ′ - k +
                                  - V- (r,k′ - k )+ γ (r) δ(k - k′) .                               (3.8)
                                     w

The partial differential equation represented in (3.6) can be transformed into an ordinary differential equation by introducing the characteristics of the Wigner function, which correspond to Newton trajectories in the physical sense. The trajectory of position (assuming parabolic bands and no electro-magnetic fields) can be parametrized with τ as

                              ∫ τ                ℏk
R (τ ;rn,kn,tn) = Rn (τ) = rn +   vg(kn) dℓ = rn +--*n(τ - tn );
                               tn                  m                                               (3.9)
            Rn (τ = tn) = rn.

The trajectory Rn is initialized by the values (rn,kn,tn) (which represent values in this context and not variables). Depending on the choice of the initialization time tn, the trajectory is either forward in time (τ > tn) or backward in time (τ < tn ). The trajectory of the wavevector is constant, since the LHS of (3.1) does not have a force term that accelerates the particle and is formally introduced as

K (τ;rn,kn, tn) = Kn (τ) = kn∀τ.                                                   (3.10)
In the following the initialization parameters are only explicitly stated where it aides clarity.

The LHS of (3.6) can represent the full derivative of fw (with respect to parameter τ). To achieve this, the following property, for an arbitrary function g, is used:

   ∫x             ∫b
d-e b g(y)dy =-de-  xg(y)dy = g (x)
dx           dx

such that

  (      ∫      )      ∫              ∫            {           }      ∫
d-- f(x)ebxg(y)dy =  dfe xb g(y)dy + g(x)ebxg(y)dyf (x) = -∂-+ g (x) f (x)e xb g(y)dy.
dx                  dx                               ∂x

The value of b (integration limit) can be chosen freely and is chosen as t0 in the following.

Both sides of (3.6) are multiplied by the term exp(                  )
 - ∫t0μ (R  (y) ,k)dy
    τ to yield

           d  (                       ∫t          )
           --- fw(R (τ;r,k,t0),k,τ )e- τ0μ(R (y),k)dy  =
∫          dτ                          ∫                                                     (3.11)
  dk ′Γ (R (τ),k,k′,τ)fw (R (τ),k′,τ)e-  t0τ μ(R (y),k)dy.
This equation can be formally integrated on the interval τ (t,t0) to finally yield the integral form of (3.6):
                        ∫
f  (r,k, t) =f   (r,k)e- tt0μ(R(y),k)dy+
 w      0    w,i∫ t   { ∫                                       ∫          }
             +    0dt′   dk ′Γ (R (τ),k,k ′,τ) f (R (t′),k′,t′)e- tt0′ μ(R(y),k)dy .                            (3.12)
                t                             w
This equation evaluates the Wigner function at the phase space point (r,k) at time t0 using the initial condition fw,i        (r,k) = fw(R                     (t;r,k, t0),k,t), which is assumed to be known at time t.

The integral equation (3.12) has a form analogous to the Chamber’s path integral [115], but with the contributions of the Wigner potential added. The introduction of the exponential term in (3.11) actually relates to an analytical summation of all out-scattering – the term μ could, alternatively, also be retained on the RHS of (3.1) and integrated otherwise, thereby suppressing the exponential term. However, the form of (3.12) is preferable as it gives a clear physical interpretation, which is useful when devising the Monte Carlo algorithm (Section 3.6).

The exponential term e- tt0μ(R(y),k)dy gives the probability of a particle to remain on its trajectory, i.e. not be scattered, from time t until time t0. Therefore, the first term of (3.12) gives the contribution of particles initialized at (R (t;r,k,t ),k)
           0 to reach point (r,k) without being scattered, whereas the second term gives the contribution of all particles scattered into the appropriate trajectory (according to in-scattering rate Γ) at time tto remain on the trajectory to reach point                                            (r,k) at time to.

3.3.1 Fredholm integral form

Fredholm integral equations of the second kind1 have the form:

               ∫
f (Q ) = fi(Q )+    dQ1K (Q, Q1)f (Q1),                                                (3.13)
where the free term fi and the kernel K are known and describe the initial/boundary conditions and the propagation of particles, respectively; Q is a multi-variable representing e.g. (r,k,t).

To express (3.12) in the form of (3.13), the integral must be augmented by the variable2 r1 to complete Q1 = (r1,k1,t1):

         f (r ,k ,t ) = f (r ,k ,t )+
          w  0  0  0    i∫ 0t0 0  0∫     ∫
                       +      dt1  dr1   dk1K (R0 (t1),k0,t0,r1,k1, t1)fw (r1,k1,t1) ;                       (3.14)
                          - ∞
                                       - ∫tt0μ(R0(y),k)dy
         fi(r0,k0,t0) ≡ f(R0 (t1),k0,t1)e  1           ∫t                                                  (3.15)
K  (r0,k0,t0,r1,k1,t1) ≡ Γ (r1,k0,k1,t1) δ(r1 - R0 (t1)) e- t10μ(R0 (y),k)dyθ(t0 - t1).                          (3.16)
This is a linear Fredholm integral equation of the second kind, which has a unique solution (if one exists) [116]. The step function θ is added to the kernel to retain the limits of the time integral. Similarly, an indicator function can be added to limit the spatial domain, but is omitted here for brevity.

A very wide variety of physical phenomena can be described by integral equations of the form (3.13) and a strong theory has evolved surrounding the solution of such Fredholm integral equations using Monte Carlo algorithms [117].

It is possible to formulate an integral representation of the WBE for the entire (global) domain. Sometimes, however, the operator of the equation under consideration is too complicated to be able to formulate a global integral representation. In such a case, local integral representations are used based on the Green’s function3 Monte Carlo algorithm [105]. The theory of Green’s function Monte Carlo algorithms keeps on developing [118, 119, 120] and is often applied to solve non-linear physical problems, e.g. [121].

3.3.2 Adjoint equation

The adjoint form of integral equations, of the form as in (3.13), is often easier to solve [106, 122].

The integral form of the WBE in (3.12) and (3.14) describes a backward-in-time equation as seen from the limits of the time integral. To obtain the corresponding forward-in-time equation the adjoint equation is required [123].

The kernel of (3.13) can be interpreted as a propagator: K(Q, Q1) describes the propagation from Q1 to Q. The adjoint equation to (3.13) solves for the function g (to which no particular physical meaning is attached at this stage) and uses the (self-)adjoint kernel K(Q1,Q ) = K       (Q,Q1 ), which exchanges the position of the integration variable:

                 ∫
g (Q1 ) = gi(Q1 )+   dQK (Q, Q1)g(Q ),                                                (3.17)
where gi is the free term.

Using the newly introduced adjoint equation it can be shown that

∫                 ∫
  dQfi (Q) g(Q) =   dQ1gi (Q1) f (Q1 ).                                              (3.18)
This result is obtained, if (3.13) is multiplied by g(Q) and integrated by dQ and (3.17) is multiplied by f(Q1) and integrated by dQ1; the two resulting equations are subtracted. It follows
  ∫                  ∫
     dQfi (Q )g (Q ) =   dQ1gi (Q1)f (Q1)iff                                            (3.19)
∫                    ∫
  dQ1fi (Q1)g (Q1) =   dQgi (Q)f (Q ).
This formulation will emerge to be very convenient for solving the computational task, described in Section 3.6.