The complex Langevin (CL) simulation method described in Section 6.4 is a versatile tool for bypassing the sign problem that arises in sampling field theories with non-positive definite weights, i.e. theories with a complex Hamiltonian H[w]. In this appendix we discuss the theoretical basis for the method.
The complex Langevin technique was devised independently by Parisi (1983) and Klauder (1983) for evaluating averages such as
(D.1)
$$<G\left(x\right)>=\frac{{\int}_{-\infty}^{\infty}\text{}dx\text{}\text{}\text{}G\left(x\right)\text{exp}[-H\left(x\right)]}{{\int}_{-\infty}^{\infty}\text{}dx\text{}\text{}\text{exp}[-H\left(x\right)]}$$
where the integration path is along the real axis for the variable
x, but the Hamiltonian
H(x) is a complex (not strictly real) function of
x. We shall begin by discussing the case where
x is a scalar, but then generalize to the more important situation where
x is replaced by an
M-vector so that the integrals in eqn (
D.1) are
M-dimensional integrals taken along the real axis.
A convenient way to rewrite eqn (D.1) is in the form
(D.2)
$$<G\left(x\right)>=\int \text{}dx\text{}\text{}G\left(x\right){P}_{c}\left(x\right)$$
where
P _{c} (x) is a so-called “complex probability weight” defined by
(D.3)
$${P}_{c}\left(x\right)=\frac{\text{exp}[-H\left(x\right)]}{\int \text{}dx\text{}\text{}\text{}\text{}\text{exp}[-H\left(x\right)]}$$
and it is understood that the path of integration is the real axis. In spite of its name,
P _{c}(
x) is not a true probability density because it is not positive semidefinite for
H(x) complex. This also implies that eqn (
D.2) cannot be directly tackled by Monte Carlo importance sampling of
P _{c}(
x) (Landau and Binder,
2000).
^{119} The basic idea behind the CL technique is to assume that one can find a
real, non-negative probability density
P(x, y) so that eqn (
D.2) can be reexpressed as
(D.4)
$$<G\left(x\right)>=\int \text{}\text{}dx\text{}\text{}\int \text{}\text{}\text{}dy\text{}\text{}\text{}\text{}G(x+iy)P(x,y)$$
Equation (
D.4) amounts to the assumption that the
line integral in eqn (
D.2) along the real axis can be exactly rewritten as an
area integral over the entire
(p.401)
complex plane of
z =
x +
iy. If such a probability density
P(x, y) exists, so that eqns (
D.2) and (
D.4) are equivalent, then eqn (
D.4) can be approximately evaluated with the importance sampling formula
(D.5)
$$<G\left(x\right)>\approx \frac{1}{{N}_{C}}\sum _{j=1}^{{N}_{C}}G\left({z}^{j}\right)$$
where
z ^{j} =
x ^{j} +
iy ^{j} for
j = 1,2,3, …
N _{C} are a set of random points in the complex plane selected from the distribution
P(x, y). The sign problem discussed in the context of eqn (
6.129) in Section
6.3 would thereby be avoided, because no complex phase factor appears in eqn (
D.5) multiplying the observable
G to be averaged.
For such a strategy to be realized, we require two things:
(D.6)
$$\begin{array}{l}{P}_{c}\left(x\right)=\int \text{}dy\text{}\text{}\text{}P(x-iy,y)\\ \phantom{\rule[-0.0ex]{2.8em}{0.0ex}}=\int \text{}\text{}\text{}d{x}^{\prime}\text{}\text{}\text{}\int \text{}\text{}d{y}^{\prime}\text{}\text{}\text{}\text{}\delta (x-{x}^{\prime}-i{y}^{\prime})P({x}^{\prime},{y}^{\prime})\end{array}$$
where we have assumed that
G(x), P _{c}(
x), and
P(
x, y) are analytic functions of their
x arguments. The second line of this expression will prove especially useful in the following. Necessary and sufficient conditions for the existence of
P(
x, y) have recently been identified (Salcedo,
1997; Weingarten,
2002). These conditions lead us to expect that most, if not all, physically realistic polymer field theory models will possess a real, non-negative distribution
P satisfying eqn (
D.6).
The complex Langevin (CL) scheme is a stochastic dynamics that, if convergent to a steady state, provides a method for sampling the distribution P(x, y) and verifying that it exists. The method amounts to writing a Langevin equation analogous to eqn (6.144), but generalizing it to trajectories z(t) = x(t) + iy(t) in the complex plane according to (Parisi, 1983; Klauder, 1983)
(D.7)
$$\begin{array}{l}\frac{d}{dt}x\left(t\right)=-\lambda \text{}\text{}\text{}\text{Re}\left[\frac{dH}{dz\left(t\right)}\right]+\eta \left(t\right)\\ \frac{d}{dt}y\left(t\right)=-\lambda \text{}\text{}\text{}\text{Im}\left[\frac{dH}{dz\left(t\right)}\right]\end{array}$$
In these equations, Re and Im denote the operations of taking the real and imaginary parts of a complex function and
dH(z)/dz is the complex derivative
(p.402)
for an analytic Hamiltonian
H(z) (Ahlfors,
1979).
^{120} The random force η
(t) is a
real, Gaussian white noise defined by (van Kampen,
1981; Kloeden and Platen,
1992)
(D.8)
$$\begin{array}{l}\phantom{\rule[-0.0ex]{2.4em}{0.0ex}}<\eta \left(t\right)>=0\\ <\eta \left(t\right)\eta \left({t}^{\prime}\right)>=2\lambda \text{}\text{}\delta (t-{t}^{\prime})\end{array}$$
The “kinetic coefficient” λ appearing in eqns (
D.7) and (
D.8) must be real and positive, although its value is arbitrary and can be absorbed into the time variable. Here we keep it explicit for reasons that will become apparent.
There are several notable features of the CL eqns (D.7). The first is the asymmetry with respect to the addition of the random force – the force is added only to the equation for the real component x(t) of the complex trajectory z(t). This asymmetry is necessary to preserve the broken symmetry of the original model in the complex x−y plane; namely, the fact that the integral in eqn (D.2) is taken along the real axis. The noise covariance in eqn (D.8) is consistent with the usual fluctuation-dissipation theorem for Brownian dynamics (van Kampen, 1981; McQuarrie, 1976), which states that the noise strength should be twice the dissipative coefficient λ appearing in front of the force terms in a Langevin equation. Another important feature of eqns (D.7) is that with the random force η(t) removed, the equations constitute a relaxational dynamics towards a saddle point z* of the model satisfying
(D.9)
$${\frac{dH\left(z\right)}{dz}|}_{z=z*}=0$$
Indeed, without the random force, the CL eqns (
D.7) reduce to the relaxation scheme eqn (
5.106) presented in Chapter
5 for the numerical computation of saddle points.
The physical content of eqns (D.7) should now be clear. In the absence of the noise, the CL equations evolve deterministically towards a nearby saddle point. However, with the random force present, the second of the two equations drives the stochastic sampling path to a value of y that is approximately consistent with the local constant phase condition^{121}
(D.10)
$$\text{Im}\frac{dH}{dz}=\frac{\partial}{\partial x}{H}_{I}(x,y)=0$$
The second equation in (
D.7) thus attempts to maintain the dynamic trajectory
z(t) on a locally constant phase path by adjusting the imaginary component
y(t). In contrast, the first equation stochastically drives the trajectory along the path through the action of the random force on the real component
x(t). As a result, if the Langevin dynamics converge to a stationary distribution
P(x, y)
(p.403)
in the complex plane, we expect
P(x, y) to have
maximum intensity centered around a constant phase ascent path passing through one or more saddle points of the model. The beauty of the technique is that it is
fully adaptive – namely, the dominant saddle point and constant phase path need not be determined in advance of running a CL simulation!
Our next task is to use eqns (D.7) to derive a Fokker–Planck equation (van Kampen, 1981) for the time-dependent probability distribution P(x, y, t) implied by the CL stochastic dynamics. The steady state solution of this equation, if it exists, is the real probability density P(x, y). Integrating both sides of eqns (D.7) from t to t + Δt leads to
(D.11)
$$\begin{array}{l}\begin{array}{c}\Delta \text{}x\equiv x(t+\Delta \text{}t)-x\left(t\right)=\lambda \text{}\text{}{\int}_{t}^{t+\Delta \text{}t}\text{}\text{}\text{}ds\text{}\text{}\text{}{F}_{R}(x\left(s\right),y\left(s\right))+\mu \end{array}\\ \Delta y\equiv y(t+\Delta \text{}t)-y\left(t\right)=\lambda \text{}\text{}{\int}_{t}^{t+\Delta \text{}t}\text{}\text{}\text{}ds\text{}\text{}\text{}{F}_{I}(x\left(s\right),y\left(s\right))\end{array}$$
where
F(
z) ≡ −
dH (
z)/
dz is the complex force. The quantity
$\mu \equiv {\int}_{t}^{t+\Delta \text{}t}\text{}\text{}\text{}ds\text{}\text{}\eta \left(s\right)$ is a new Gaussian random force acting over the timestep with mean and variance that follow immediately from eqn (
D.8):
(D.12)
$$\begin{array}{l}\phantom{\rule[-0.0ex]{.5em}{0.0ex}}<\mu >=0\\ <{\mu}^{2}>={\int}_{t}^{t+\Delta \text{}t}\text{}\text{}\text{}ds\text{}{\int}_{t}^{t+\Delta \text{}t}\text{}\text{}\text{}d{s}^{\prime}\text{}\text{}\text{}<\eta \left(s\right)\eta \left({s}^{\prime}\right)>=2\lambda \Delta t\end{array}$$
It is important to note that μ is characteristically
O((Δ
t)
^{1/2}). Assuming continuity of
dH/dz, eqns (
D.11) can be approximated by
(D.13)
$$\begin{array}{l}\Delta \text{}x=\lambda \Delta t\text{}\text{}\text{}{F}_{R}+\mu \\ \Delta \text{}y=\lambda \Delta t\text{}\text{}\text{}{F}_{I}\end{array}$$
with errors that are
O((Δ
t)
^{2}). Using these equations, it is straightforward to show that the first two moments of the random variables Δ
x and Δ
y, averaged over all realizations of the Gaussian force μ, are given to
O(Δ
t) by
(D.14)
$$\begin{array}{l}\phantom{\rule[-0.0ex]{1.5em}{0.0ex}}<\Delta \text{}x>=\lambda \Delta t\text{}\text{}\text{}{F}_{R},\text{}\text{}\text{}\text{}<\Delta y>=\lambda \Delta t\text{}\text{}{F}_{I}\\ <{\left(\Delta \text{}x\right)}^{2}>=2\lambda \Delta t,\text{}\text{}\text{}\text{}<{\left(\Delta y\right)}^{2}>=<\Delta \text{}x\Delta \text{}y>=0\end{array}$$
These results can now be used to derive a Fokker–Planck equation for the probability density P(x, y, t). The starting point is a Chapman–Kolmogorov (CK) equation strictly analogous to eqn (2.58) for the continuous Gaussian chain. Defining a two-component state vector according to x = (x, y)^{T}, the CK equation can be written
(D.15)
$$P(\mathbf{x},t+\Delta t)=\int d\left(\Delta \mathbf{x}\right)\Phi (\Delta \mathbf{x};\mathbf{x}-\Delta \mathbf{x})P(\mathbf{x}-\Delta \mathbf{x},t)$$
where Φ(Δ
x; x) is the transition probability density for a displacement Δ
x in the complex plane, starting at the point
x, over a time interval of Δ
t. This function
(p.404)
is normalized so that ∫
d(Δ
x) Φ = 1 and its first two moments are summarized by eqn (
D.14). Following the procedure outlined in Section
2.4, eqn (
D.15) can be converted to a Fokker–Planck equation by expanding the left-hand side in powers of Δ
t to
O(Δ
t) and expanding the right-hand side in powers of Δ
x to
O((Δ
x)
^{2}) =
O(Δ
t). This leads to
(D.16)
$$\begin{array}{l}\Delta t\frac{\partial}{\partial t}P(\mathbf{x},t)=-{\nabla}_{\mathbf{x}}\cdot \left[<\Delta \mathbf{x}>P(\mathbf{x},t)\right]\\ \phantom{\rule[-0.0ex]{7.0em}{0.0ex}}+\frac{1}{2!}{\nabla}_{\mathbf{x}}{\nabla}_{\mathbf{x}}:\left[<\Delta \mathbf{x}\Delta \mathbf{x}>P(\mathbf{x},t)\right]+O\left({\left(\Delta t\right)}^{2}\right)\end{array}$$
Finally, substituting eqn (
D.14) for the moments of Φ produces the desired Fokker–Planck equation for the CL process
(D.17)
$$\begin{array}{l}\frac{\partial}{\partial t}P(x,y,t)=-\lambda \frac{\partial}{\partial t}\left[{F}_{R}(x,y)P(x,y,t)\right]-\lambda \frac{\partial}{\partial y}\left[{F}_{I}(x,y)P(x,y,t)\right]\\ \phantom{\rule[-0.0ex]{7.0em}{0.0ex}}+\lambda \frac{{\partial}^{2}}{\partial {x}^{2}}P(x,y,t)\end{array}$$
In spite of its linearity, this Fokker–Planck equation apparently has no closed form solution for an arbitrary force
F(z), even in the steady state limit where
P(x, y, t) →
P(x, y).
Our final task is to prove that if a steady state solution P(x, y) of the above equation exists, then averages computed with this solution using eqn (D.4) are equivalent to averages computed with eqn (D.2) using the complex weight P _{c}(x) (Schoenmaker, 1987; Lee, 1994). This can be shown by combining eqns (D.6) and (D.17). Specifically, applying the operation ∫ dx′ ∫ dy′ δ(x − x′ − iy′) to both sides of the Fokker–Planck equation written for P(x′, y′, t) leads to
(D.18)
$$\frac{\partial}{\partial t}{P}_{c}(x,t)={T}_{1}(x,t)+{T}_{2}(x,t)+{T}_{3}(x,t)$$
where
P _{c} (
x, t) ≡ ∫
dy P(
x − iy, y, t) and
T _{j}(
x, t) is the function obtained by applying the indicated operation to the
j th term on the right-hand side of eqn (
D.17).
T _{1} can be manipulated as follows:
(D.19)
$$\begin{array}{l}{\begin{array}{c}T\end{array}}_{1}(x,t)=\lambda \int d{x}^{\prime}\int d{y}^{\prime}\delta (x-{x}^{\prime}-i{y}^{\prime})\frac{\partial}{\partial {x}^{\prime}}\left[\text{Re}\left(\frac{dH\left({x}^{\prime}+i{y}^{\prime}\right)}{d({x}^{\prime}+i{y}^{\prime})}\right)P({x}^{\prime},{y}^{\prime},t)\right]\\ \phantom{\rule[-0.0ex]{4.0em}{0.0ex}}=\lambda \int d{y}^{\prime}\frac{\partial}{\partial x}\left[\text{Re}\left(\frac{dH\left(x\right)}{dx}\right)P(x-i{y}^{\prime},{y}^{\prime},t)\right]\\ \phantom{\rule[-0.0ex]{4.0em}{0.0ex}}=\lambda \frac{\partial}{\partial x}\left[\text{Re}\left(\frac{dH\left(x\right)}{dx}\right){P}_{c}(x,t)\right]\end{array}$$
Similarly, the second term can be written
(p.405)
(D.20)
$$\begin{array}{l}{\begin{array}{c}T\end{array}}_{2}(x,t)=\lambda \int d{x}^{\prime}\int d{y}^{\prime}\delta (x-{x}^{\prime}-i{y}^{\prime})\frac{\partial}{\partial {y}^{\prime}}\left[\text{Im}\left(\frac{dH\left({x}^{\prime}+i{y}^{\prime}\right)}{d({x}^{\prime}+i{y}^{\prime})}\right)P({x}^{\prime},{y}^{\prime},t)\right]\\ \phantom{\rule[-0.0ex]{4.0em}{0.0ex}}=i\lambda \int d{x}^{\prime}\int d{y}^{\prime}\delta (x-{x}^{\prime}-i{y}^{\prime})\frac{\partial}{\partial x}\left[\text{Im}\left(\frac{dH\left(x\right)}{dx}\right)P({x}^{\prime},{y}^{\prime},t)\right]\\ \phantom{\rule[-0.0ex]{4.0em}{0.0ex}}=i\lambda \frac{\partial}{\partial x}\left[\text{Im}\left(\frac{dH\left(x\right)}{dx}\right){P}_{c}(x,t)\right]\end{array}$$
Finally, the last term is
(D.21)
$$\begin{array}{l}{\begin{array}{c}T\end{array}}_{3}(x,t)=\lambda \int d{x}^{\prime}\int d{y}^{\prime}\delta (x-{x}^{\prime}-i{y}^{\prime})\frac{{\partial}^{2}}{{\left(\partial {x}^{\prime}\right)}^{2}}P({x}^{\prime},{y}^{\prime},t)\\ \phantom{\rule[-0.0ex]{4.0em}{0.0ex}}=\lambda \int d{y}^{\prime}\frac{{\partial}^{2}}{\partial {x}^{2}}P(x,i{y}^{\prime},{y}^{\prime},t)\\ \phantom{\rule[-0.0ex]{4.0em}{0.0ex}}=\lambda \frac{{\partial}^{2}}{\partial {x}^{2}}{P}_{c}(x,t)\end{array}$$
Combining these results, we see that the function
P _{c}(
x, t) satisfies the following
complex Fokker–Planck (FP) equation:
(D.22)
$$\frac{\partial}{\partial t}{P}_{c}(x,t)=\frac{\partial}{\partial x}\lambda [\frac{\partial}{\partial x}+\frac{dH\left(x\right)}{dx}]{P}_{c}(x,t)$$
This equation has a complex steady state solution
P _{c}(
x)∞ exp[−
H(x)] corresponding to eqn (
D.3), but also a second “spurious” steady state (Lee,
1994)
(D.23)
$${P}_{spur}\left(x\right)\propto \text{exp}[-H\left(x\right)]{\int}^{x}dy\text{}\text{}\text{}\text{}\text{exp}\left[H\left(y\right)\right]$$
This spurious solution is usually not relevant, because it leads to expectation values that are incompatible with the most common choices of boundary conditions, e.g. periodic. Thus, if the real FP eqn (
D.17)
converges to a steady state P(x, y), then the associated function
P _{c}(
x) = ∫
dy P(
x − iy, y) corresponds to the desired complex probability given by eqn (
D.3). When this condition is met, a “time average” computed with eqn (
D.5) along a complex Langevin trajectory will converge to the ensemble average (
D.2) in the limit of
N _{C} → ∞.
Unfortunately, an analytical proof that eqn (D.17) has a steady state solution is not in hand. Nevertheless, there is a simple test to ensure that a numerical CL simulation is working properly (Gausterer and Lee, 1993; Lee, 1994). If the expectation values for analytic observables G(x) become time independent over the course of a CL simulation, then it can be proven that these values are correct and in agreement with eqn (D.2). In practice, we have not encountered convergence problems in CL simulations for any of the models considered in Chapter 4.
The complex Langevin equations given in (D.7) constitute the “standard” CL approach. However, there are several extensions of the formalism that are potentially useful in conducting field-theoretic simulations. The first of these is
(p.406)
a generalization to include a noise source acting on both the real and imaginary parts of the field:^{122}
(D.24)
$$\begin{array}{c}\frac{d}{dt}x\left(t\right)=-\lambda \text{}\text{}\text{}\text{Re}\left[\frac{dH}{dz\left(t\right)}\right]+{\eta}_{R}\left(t\right)\\ \frac{d}{dt}y\left(t\right)=-\lambda \text{}\text{}\text{}\text{Im}\left[\frac{dH}{dz\left(t\right)}\right]+{\eta}_{I}\left(t\right)\end{array}$$
The random forces η
_{R}(
t) and η
_{I}(
t) can be taken to be real, Gaussian white noise processes with vanishing mean values, <η
R(
t)> = <η
I(
t)> = 0. The covariance of the noise can be further generalized from eqn (
D.8) to
(D.25)
$$\begin{array}{l}<{\eta}_{R}\left(t\right){\eta}_{R}\left({t}^{\prime}\right)>=2(\lambda +\epsilon )\delta (t-{t}^{\prime})\\ \phantom{\rule[-0.0ex]{.5em}{0.0ex}}<{\eta}_{I}\left(t\right){\eta}_{I}\left({t}^{\prime}\right)>=2\epsilon \text{}\text{}\delta (t-{t}^{\prime})\\ \phantom{\rule[-0.0ex]{0.3em}{0.0ex}}<{\eta}_{R}\left(t\right){\eta}_{I}\left({t}^{\prime}\right)>=0\end{array}$$
where λ > 0 and ε ≥ 0 are real parameters that determine the relative strengths of the two noise components. Equations (
D.24) and (
D.25) evidently reduce to the “standard” CL eqns (
D.7)–(
D.8) in the special case of ε = 0. Since the parameter λ can be absorbed into the time scale of the stochastic process, ε is effectively a free parameter that can be adjusted to optimize the performance of a CL simulation.
The role of ε can be established by deriving the Fokker–Planck equation corresponding to these generalized CL equations. By repeating the steps leading to eqn (D.17) it is straightforward to show that eqns (D.24)–(D.25) are consistent with the Fokker–Planck equation
(D.26)
$$\begin{array}{c}\frac{\partial}{\partial t}P(x,y,t)=-\lambda \frac{\partial}{\partial x}\left[{F}_{R}(x,y)P(x,y,t)\right]-\lambda \frac{\partial}{\partial y}\left[{F}_{I}(x,y)P(x,y,t)\right]\\ +(\lambda +\epsilon )\frac{{\partial}^{2}}{\partial {x}^{2}}P(x,y,t)+\epsilon \frac{{\partial}^{2}}{\partial {y}^{2}}P(x,y,t)\end{array}$$
Moreover, it can be shown that if this generalized Fokker–Planck equation has a steady state
P(
x, y), then the associated function
P _{c}(
x) = ∫
dy P(
x – iy, y) reduces to eqn (
D.3). Thus, the generalized CL equations are a suitable alternative to the standard CL approach. For non-zero ε, however, eqn (
D.26) shows that the generalized scheme contains an extra dispersive term
$\epsilon {\partial}_{y}^{2}P$ that will tend to
smooth the steady state distribution in the variable
y.
To illustrate the role of finite ε, we return to the “toy” model of eqn (5.7) for the case of p = 1, i.e. H(x) = ix + x ^{2}/2. The steady state solution of eqn (D.26) for this simple quadratic Hamiltonian is
(D.27)
$$P(x,y)\sim \text{exp}(-\frac{\lambda {x}^{2}}{2(\lambda +\epsilon )}-\frac{\lambda {(y+1)}^{2}}{2\epsilon})$$
(p.407)
apart from a normalization constant. The distribution function has a Gaussian ridge of maximum probability centered on the line
y = − 1 in the complex plane that passes through the saddle point,
z* = −
i. This line is, of course, the constant-phase ascent path for the model. The width of the ridge normal to this line is
O(ε
^{1/2}), while the decay in probability density away from the saddle point along the line
y = −
1 is much slower for ε ≪ λ. For ε → 0 +, the conventional CL theory is recovered and eqn (
D.27) reduces to the singular distribution
(D.28)
$$P(x,y)\sim \delta (y+1)\text{exp}(-{x}^{2}/2)$$
At least for this simple model, it is clear that the generalized CL theory with ε > 0 produces a smoother steady state distribution function that is centered on the constant-phase ascent path.
The CL theory can be immediately extended to the multi-dimensional case of model Hamiltonians H(x), where x = (x _{1}, x _{2}, …, x _{M})^{T} is a real M-vector and H is complex valued. In particular, the generalized CL eqns (D.24)–(D.25) can be written in the multi-variate case as
(D.29)
$$\begin{array}{c}\frac{d}{dt}\mathbf{x}\left(t\right)=-\lambda \text{}\text{}\text{}\text{Re}\left[\frac{\partial H}{\partial \mathbf{z}\left(t\right)}\right]+{\eta}_{R}\left(t\right)\\ \frac{d}{dt}\mathbf{y}\left(t\right)=-\lambda \text{}\text{}\text{}\text{Im}\left[\frac{\partial H}{\partial \mathbf{z}\left(t\right)}\right]+{\eta}_{I}\left(t\right)\end{array}$$
where
z =
x +
iy is a complex
M-vector. The Gaussian random forces η
_{R}(
t) and η
_{I} (
t) are
M-vectors with vanishing mean and covariance matrices given by
(D.30)
$$\begin{array}{l}<{\eta}_{R}\left(t\right){\eta}_{R}\left({t}^{\prime}\right)>=2(\lambda +\epsilon )\text{}\text{}\mathbf{1}\delta (t-{t}^{\prime})\\ \phantom{\rule[-0.0ex]{.5em}{0.0ex}}<{\eta}_{I}\left(t\right){\eta}_{I}\left({t}^{\prime}\right)>=2\epsilon \text{}\text{}\mathbf{1}\delta (t-{t}^{\prime})\\ \phantom{\rule[-0.0ex]{0.3em}{0.0ex}}<{\eta}_{R}\left(t\right){\eta}_{I}\left({t}^{\prime}\right)>=0\end{array}$$
where
1 is the
M ×
M unit tensor. By setting ε = 0, the above scheme reduces to the standard multi-variate CL theory.
Finally, in the multi-dimensional case, the parameters λ and ε can be replaced by a positive definite, M × M “kinetic coeffficient” matrix ε and a positive semidefinite, M × M “noise” matrix ε. This generalization amounts to the scheme
(D.31)
$$\begin{array}{c}\frac{d}{dt}\mathbf{x}\left(t\right)=-\lambda \text{}\text{}\cdot \text{}\text{Re}\left[\frac{\partial H}{\partial \mathbf{z}\left(t\right)}\right]+{\eta}_{R}\left(t\right)\\ \frac{d}{dt}\mathbf{y}\left(t\right)=-\lambda \text{}\cdot \text{}\text{Im}\left[\frac{\partial H}{\partial \mathbf{z}\left(t\right)}\right]+{\eta}_{I}\left(t\right)\end{array}$$
with
(D.32)
$$\begin{array}{l}<{\eta}_{R}\left(t\right){\eta}_{R}\left({t}^{\prime}\right)>=2(\lambda +\epsilon )\text{}\text{}\delta (t-{t}^{\prime})\\ <{\eta}_{I}\left(t\right){\eta}_{I}\left({t}^{\prime}\right)>=2\epsilon \text{}\text{}\delta (t-{t}^{\prime})\\ <{\eta}_{R}\left(t\right){\eta}_{I}\left({t}^{\prime}\right)>=0\end{array}$$
(p.408)
Again it can be proven that if the above CL equations converge to a steady state
P(
x, y), the steady state solution is consistent with the proper complex probability weight
P _{c}(
x) ∼ exp[−
H(
x)]. By adjusting the form of the matrices λ and ε it may be possible to achieve better performance in numerical simulations than the simplest choice of λ = λ
1 and ε = ε
1. A potentially important class of rank-
M matrices is those that are translationally invariant on a
d-dimensional collocation grid and can be diagonalized by a discrete Fourier transform. By employing such matrices, efficient pseudo-spectral numerical solutions of the CL equations can be achieved.