Langevin updates and lattice numerics

Intuition Beginner

A Langevin update is a recipe for sampling a probability distribution by simulating a noisy ball rolling downhill on a potential landscape. The ball's position represents a field configuration. The landscape is the action, and lower-action configurations are deeper valleys. You let the ball roll a small step downhill, then give it a small random kick, then repeat. After a long time, the position of the ball is a sample from the Gibbs distribution that weights configurations by $e^{-S}/Z$.

This is the whole game of stochastic field theory on a lattice. You replace the continuous space of field configurations with a finite grid, you write down the action as a function of the values on grid points, and you simulate the noisy-ball recipe on the resulting finite-dimensional landscape. The random kicks are calibrated against the step size so that the equilibrium is the Euclidean path-integral measure of the field theory.

Why does this matter? The path-integral measure on field configurations is hard to make sense of directly; the noisy-ball construction gives you a concrete computer-implementable way to draw samples from it. Stochastic quantisation, due to Parisi and Wu in 1981, takes this seriously as a quantisation procedure: the fictitious time of the noisy-ball simulation is an extra dimension, and the original physics lives at the equilibrium of the noisy dynamics. Gauge theories benefit because the equilibrium measure is gauge-invariant; the noisy ball does not need to be told which gauge to pick.

Visual Beginner

A schematic showing a small two-dimensional lattice on the left with field values at each site, and on the right a noisy trajectory in stochastic time $\tau$ that visits many configurations and relaxes toward the equilibrium distribution. Arrows on the trajectory indicate that each step is a small downhill move on the action landscape plus a small random kick.

The picture captures the essence: a lattice fixes the space of field configurations, the action assigns a number to each configuration, and the noisy ball samples that number via $e^{-S}/Z$ as the long-time equilibrium of a simulated stochastic flow.

Worked example Beginner

Sample the harmonic oscillator on a one-site lattice using the Langevin recipe and check the variance of the samples.

Step 1. The configuration is a single number $\varphi \in \mathbb{R}$, the action is $S(\varphi )=\frac{1}{2}{\omega }^2{\varphi }^2$ with $\omega =1$, and the temperature is $T=1$. The Gibbs distribution is $e^{-S(\varphi )/T}/Z=e^{-{\varphi }^2/2}/\sqrt{2\pi }$, a standard normal. We expect the equilibrium variance of $\varphi$ to equal $T/{\omega }^2=1$.

Step 2. Pick a step size $\Delta \tau =0.01$. Start from ${\varphi }_0=0$. At each step, compute the downhill direction as the negative gradient of the action: $-S^{\prime }(\varphi )=-{\omega }^2\varphi =-\varphi$. Draw a random kick ${\xi }_n$ from a standard normal distribution.

Step 3. The Euler-Maruyama update reads ${\varphi }_{n+1}={\varphi }_n-\Delta \tau \cdot {\varphi }_n+\sqrt{2\Delta \tau }\cdot {\xi }_n$. With $\Delta \tau =0.01$ this becomes ${\varphi }_{n+1}=0.99{\varphi }_n+0.1414\cdot {\xi }_n$. Each step moves $\varphi$ slightly toward the origin and adds a small random kick.

Step 4. Run the update for $N=10000$ steps, starting from ${\varphi }_0=0$. The first few hundred samples are biased by the initial condition. Discard the first $1000$ steps as burn-in. Collect the remaining $9000$ samples ${\varphi }_{1001},{\varphi }_{1002},\dots ,{\varphi }_{10000}$.

Step 5. Compute the sample mean by adding up all $9000$ values of ${\varphi }_n$ and dividing by $9000$; the answer is close to $0$. Compute the sample variance by adding up all $9000$ values of $({\varphi }_n-\bar{\varphi }{)}^2$ and dividing by $9000$; the answer is close to $1$. The variance is close to the theoretical value $T/{\omega }^2=1$. The Gibbs distribution has been recovered as the equilibrium of the noisy dynamics.

What this tells us: the Langevin recipe samples the Gibbs distribution as the long-time equilibrium of a noisy dynamics. The simulation step is short and explicit; the only inputs are the gradient of the action and a random-number generator. The same recipe extends to a lattice with many sites, where the gradient at each site is the local-action derivative and the random kicks are independent across sites.

Check your understanding Beginner

Formal definition Intermediate+

Fix a finite lattice $\Lambda =(a\mathbb{Z}/L\mathbb{Z}{)}^d$ in $d$ spatial dimensions with lattice spacing $a$ and side length $L$, so $|\Lambda |=(L/a{)}^d$ is the number of sites. A real scalar field configuration is a function $\varphi :\Lambda \to \mathbb{R}$, which we identify with a vector in ${\mathbb{R}}^{|\Lambda |}$.

The lattice action is the function $$ S[\phi] = \sum_{x \in \Lambda} \left[ \tfrac{1}{2} |\nabla \phi(x)|^2 + V(\phi(x)) \right], $$ where the lattice gradient is $(\nabla \varphi {)}_{\mu }(x)=(\varphi (x+a{\hat{e}}_{\mu })-\varphi (x))/a$ and $V:\mathbb{R}\to \mathbb{R}$ is a fixed scalar potential — say $V(\varphi )=\frac{1}{2}{m}^2{\varphi }^2+\frac{\lambda }{4!}{\varphi }^4$ for a massive ${\varphi }^4$ theory. The action defines the Euclidean Gibbs measure $$ d\mu_S(\phi) = \frac{1}{Z} e^{-S[\phi]} \prod_{x \in \Lambda} d\phi(x), \qquad Z = \int_{\mathbb{R}^{|\Lambda|}} e^{-S[\phi]} \prod_{x \in \Lambda} d\phi(x). $$ This is a probability measure on ${\mathbb{R}}^{|\Lambda |}$, and the expectation ${\mathbb{E}}_{{\mu }_S}[f]=Z^{-1}\int f(\varphi )e^{-S[\varphi ]}d\varphi$ recovers the lattice path-integral expectation value of an observable $f$.

The Langevin equation with stochastic time $\tau \ge 0$ is the system of stochastic differential equations $$ \frac{d \phi(x, \tau)}{d \tau} = -\frac{\delta S}{\delta \phi(x)}[\phi(\cdot, \tau)] + \eta(x, \tau), \qquad x \in \Lambda, $$ where $\eta (x,\tau )$ is a Gaussian white-noise field on $\Lambda \times {\mathbb{R}}_{+}$ with mean zero and covariance $$ \mathbb{E}[\eta(x, \tau), \eta(y, \tau')] = 2, \delta_{xy}, \delta(\tau - \tau'). $$ The drift is the negative gradient of the action; the noise is independent across sites and white in stochastic time. The variance factor $2$ is the conventional normalisation, chosen so that the equilibrium temperature equals $1$ in our units.

The Euler-Maruyama discretisation with step size $\Delta \tau >0$ replaces the continuous SDE by the explicit update rule $$ \phi_{n+1}(x) = \phi_n(x) - \Delta\tau \cdot \frac{\partial S}{\partial \phi(x)}[\phi_n] + \sqrt{2 \Delta\tau} \cdot \xi_n(x), $$ where ${\xi }_n(x)\sim \mathcal{N}(0,1)$ are i.i.d. standard normal across sites $x\in \Lambda$ and stochastic-time steps $n\in \mathbb{N}$. After $N$ steps the discrete trajectory ${\varphi }_N\in {\mathbb{R}}^{|\Lambda |}$ is a random variable whose distribution converges, in a sense made precise in the next section, to the Gibbs measure ${\mu }_S$.

A lattice observable is a measurable function $f:{\mathbb{R}}^{|\Lambda |}\to \mathbb{R}$ — typically a polynomial in the field values such as $f(\varphi )=\varphi (x_0)\varphi (x_1)$ for the two-point function. The Langevin estimator of the expectation $⟨f{⟩}_{{\mu }_S}$ is the trajectory average $$ \hat{f}N = \frac{1}{N - N_0} \sum{n = N_0 + 1}^{N} f(\phi_n), $$ where $N_0$ is a burn-in cutoff chosen to discard the transient relaxation of ${\varphi }_n$ from its initial condition. Under coercivity hypotheses on $V$ (e.g. $V(\varphi )\to \infty$ as $|\varphi |\to \infty$), the ergodic theorem for Langevin SDEs guarantees ${\hat{f}}_N\to ⟨f{⟩}_{{\mu }_S}$ almost surely as $N\to \infty$ followed by $\Delta \tau \to 0$.

Counterexamples to common slips

• The noise normalisation matters. The variance factor $2$ in $\mathbb{E}[\eta {\eta }^{\prime }]=2\delta \delta$ is set by the choice of "temperature" $T=1$ in the Gibbs measure $e^{-S/T}/Z$; rescaling the noise to $\eta \to \sqrt{2T}\eta$ shifts the equilibrium to $e^{-S/T}/Z$ with a different temperature. Mixing conventions across units is a common source of factor-of-two errors.
• The drift sign matters. With $+\delta S/\delta \varphi$ instead of $-\delta S/\delta \varphi$, the flow goes uphill on the action landscape; the equilibrium does not exist (or, with periodic boundary conditions in $\tau$, would be $e^{+S}/Z$, which is non-normalisable for actions unbounded above).
• The Euler-Maruyama scheme has order-$\Delta \tau$ bias in the stationary distribution. The samples come from a perturbed measure $e^{-S-\Delta \tau \cdot R}/Z$ with $R$ a correction depending on the second derivatives of $S$. The bias vanishes as $\Delta \tau \to 0$, or one switches to an exact-sampling variant (hybrid Monte Carlo) that adds a Metropolis accept/reject.

Key theorem with proof Intermediate+

Theorem (Langevin equilibrium identification; Chatterjee Lecture 11 ^[Chatterjee], Parisi-Wu 1981 ^[Parisi-Wu]). Let $\Lambda$ be a finite lattice, $S:{\mathbb{R}}^{|\Lambda |}\to \mathbb{R}$ a $C^2$ action satisfying the coercivity condition $S(\varphi )\to \infty$ as $|\varphi |\to \infty$ and $\int e^{-S}d\varphi <\infty$. Let $\varphi (\tau )\in {\mathbb{R}}^{|\Lambda |}$ solve the Langevin SDE* $$ d\phi(\tau) = -\nabla S(\phi(\tau)), d\tau + \sqrt{2}, dW(\tau), $$ *with $W$ a standard $|\Lambda |$-dimensional Brownian motion. Then the unique stationary distribution of $\varphi (\tau )$ is the Gibbs measure* $$ d\mu_S(\phi) = \frac{1}{Z} e^{-S(\phi)} d\phi, \qquad Z = \int_{\mathbb{R}^{|\Lambda|}} e^{-S(\phi)} d\phi. $$ *For every continuous bounded observable $f:{\mathbb{R}}^{|\Lambda |}\to \mathbb{R}$ and every initial condition $\varphi (0)$, $\mathbb{E}[f(\varphi (\tau ))]\to {\mathbb{E}}_{{\mu }_S}[f]$ as $\tau \to \infty$, with the convergence exponentially fast in $\tau$ when ${\nabla }^{2}S\ge \rho I$ for some $\rho >0$ (Bakry-Émery curvature-dimension criterion CD($\rho ,\infty$)).

Proof. The Langevin SDE is a special case of the Itô SDE $d\varphi =b(\varphi )d\tau +\sigma (\varphi )dW$ with drift $b(\varphi )=-\nabla S(\varphi )$ and constant diffusion matrix $\sigma =\sqrt{2}I$. The Itô formula gives, for any $C^2$ test function $f$, $$ df(\phi(\tau)) = \nabla f(\phi(\tau)) \cdot d\phi(\tau) + \tfrac{1}{2} \sum_{x, y} (\sigma \sigma^T){xy} , \tfrac{\partial^2 f}{\partial \phi(x) \partial \phi(y)} d\tau. $$ Taking the expectation, the stochastic-integral term drops, and substituting $\sigma {\sigma }^T=2I$ and $b=-\nabla S$ gives the generator $L$ of the process: $$ (Lf)(\phi) = -\nabla S(\phi) \cdot \nabla f(\phi) + \Delta f(\phi), $$ where $\Delta ={\sum }_x{\partial }^2/\partial \varphi (x{)}^2$ is the Laplacian on ${\mathbb{R}}^{|\Lambda |}$. The Fokker-Planck equation for the time-dependent density $p(\varphi ,\tau )$ of $\varphi (\tau )$ is $$ \partial\tau p = L^* p = \nabla \cdot (p \nabla S) + \Delta p, $$ where $L^{\ast }$ is the formal $L^2$-adjoint of $L$.

A stationary density $p_{\mathrm{eq}}$ satisfies $L^{\ast }{p}_{\mathrm{eq}}=0$. Substitute $p_{\mathrm{eq}}(\varphi )=e^{-S(\varphi )}/Z$: $$ L^* p_{\mathrm{eq}} = \nabla \cdot (e^{-S}/Z \cdot \nabla S) + \Delta(e^{-S}/Z). $$ Compute the second term using $\nabla e^{-S}=-e^{-S}\nabla S$ and $\Delta e^{-S}=-\nabla \cdot (e^{-S}\nabla S)=-e^{-S}\Delta S+e^{-S}|\nabla S{|}^2$. So $$ \Delta(e^{-S}/Z) = (1/Z)(- e^{-S} \Delta S + e^{-S} |\nabla S|^2), $$ and $$ \nabla \cdot (e^{-S}/Z \cdot \nabla S) = (1/Z)(- e^{-S} |\nabla S|^2 + e^{-S} \Delta S). $$ The two terms cancel: $L^{\ast }(e^{-S}/Z)=0$. So the Gibbs density is a stationary solution of the Fokker-Planck equation.

For uniqueness, recall from 08.10.02 that the Fokker-Planck operator $L^{\ast }=\Delta +\nabla \cdot (\cdot \nabla S)$ on $L^2({\mathbb{R}}^{|\Lambda |})$ has a one-dimensional kernel spanned by $e^{-S}/Z$ whenever $S$ is coercive, by the spectral gap of the symmetric form on weighted $L^2$. The Bakry-Émery criterion ${\nabla }^{2}S\ge \rho I$ for $\rho >0$ identifies $\rho$ as a lower bound for the spectral gap, giving exponential convergence ${\mathrm{Ent}}_{{\mu }_S}(p(\tau )/p_{\mathrm{eq}})\le e^{-2\rho \tau }{\mathrm{Ent}}_{{\mu }_S}(p(0)/p_{\mathrm{eq}})$ of the relative entropy.

To conclude: for every continuous bounded $f$, write $\mathbb{E}[f(\varphi (\tau ))]-{\mathbb{E}}_{{\mu }_S}[f]=\int f\cdot (p(\varphi ,\tau )-p_{\mathrm{eq}}(\varphi ))d\varphi$. The Csiszár-Kullback-Pinsker inequality bounds the $L^1$ distance $\Vert p(\tau )-p_{\mathrm{eq}}{\Vert }_{L^1}$ by the relative entropy, which decays exponentially. So $\mathbb{E}[f(\varphi (\tau ))]\to {\mathbb{E}}_{{\mu }_S}[f]$ at the rate $e^{-\rho \tau }$. $\square$

Bridge. The Langevin equilibrium identification builds toward every numerical implementation of Euclidean field theory: the same identification on a finite lattice ${\mathbb{R}}^{|\Lambda |}$ produces the path-integral measure of a regulated quantum field theory as the long-time limit of a stochastic dynamics, and this is exactly the foundational reason stochastic quantisation is a viable construction route. The central insight is exactly that the Fokker-Planck generator $L^{\ast }$ of the Langevin SDE is symmetric with respect to the inner product on weighted $L^2({\mu }_S)$, putting these together with the spectral-gap argument from 08.10.02 to give a unique exponentially-attracting stationary measure. This is exactly the same detailed-balance / reversibility identification that appears again in 08.10.07 when Wightman functions of a Euclidean field theory are recovered from Osterwalder-Schrader reflection positivity. The foundational reason the cancellation $L^{\ast }(e^{-S}/Z)=0$ works is the algebraic identity $\nabla \cdot (e^{-S}\nabla S)=-\Delta e^{-S}$, which generalises to the statement that the Fokker-Planck generator is exactly $L^{\ast }=-\nabla \cdot (e^{-S}\nabla (e^S\cdot ))$ in the gradient-flow form — the central insight that identifies Langevin dynamics with gradient descent in the Wasserstein metric, the bridge between stochastic processes and optimisation theory. Putting these together with the Euler-Maruyama discretisation gives a computer-implementable sampler whose bias is order $\Delta \tau$, and the bridge to exact sampling is the hybrid Monte Carlo correction discussed below.

Exercises Intermediate+

Advanced results Master

Theorem (Parisi-Wu stochastic quantisation; Parisi-Wu 1981 ^[Parisi-Wu], Damgaard-Hüffel 1987 ^{[Damgaard-Huffel]}). Let $S[\varphi ]$ be a lattice action on ${\mathbb{R}}^{|\Lambda |}$. The Langevin SDE $d\varphi =-\nabla S(\varphi )d\tau +\sqrt{2}dW$ in fictitious stochastic time $\tau$ has the Euclidean Gibbs measure $e^{-S[\varphi ]}/Z$ as its unique stationary distribution, and the long-time average of any observable along a Langevin trajectory equals the path-integral expectation value of the lattice quantum field theory. The construction extends to lattice gauge theory: on the configuration manifold $G^{|E|}$ of group-valued link variables (with $G$ a compact Lie group), the Lie-algebra-valued Langevin SDE with the appropriate Haar-measure drift correction has the Wilson measure as its equilibrium, without conventional gauge fixing.

The proof is the finite-dimensional Fokker-Planck argument of the Intermediate-tier theorem, applied first to ${\mathbb{R}}^{|\Lambda |}$ for scalar fields and then to the compact group manifold $G^{|E|}$ for gauge fields. On the group manifold, the Fokker-Planck operator picks up curvature-of-Haar-measure corrections, and the Stratonovich-to-Itô conversion gives a drift that respects the Haar volume; the Wilson measure $e^{-S_W}/Z\cdot dU$ then satisfies $L^{\ast }(e^{-S_W}/Z\cdot dU)=0$ by the same algebraic-identity computation, with the group Laplacian replacing the Euclidean Laplacian. Damgaard-Hüffel's 1987 review consolidates the framework and the perturbative agreement with the conventional Faddeev-Popov gauge-fixed path integral order-by-order in coupling.

Theorem (Zwanziger stochastic gauge fixing; Zwanziger 1981 ^[Zwanziger]). Augment the gauge-theory Langevin equation with an additional drift along gauge orbits, parametrised by an arbitrary gauge-fixing function $G[U]$. The augmented SDE $$ \frac{dU(x, \mu, \tau)}{d\tau} = \big( -i \nabla_a S_W - i \nabla_a G[U] \big) T^a U + i\eta_a T^a U $$ has the same gauge-invariant Wilson equilibrium measure as the unmodified SDE. The extra drift selects a representative on each gauge orbit without changing the expectation value of any gauge-invariant observable, providing a Lorentz-covariant computational alternative that does not require a Faddeev-Popov ghost determinant and does not suffer the Gribov ambiguity.

The argument: the Wilson action is gauge-invariant, so the gauge-orbit drift $-{\nabla }_{a}G$ has no projection onto the gauge-invariant equilibrium density. The Fokker-Planck operator splits as $L^{\ast }=L_{\mathrm{gauge}-\mathrm{inv}}^{\ast }+L_{\mathrm{gauge}-\mathrm{orbit}}^{\ast }$, and the gauge-invariant part is what was identified in the unmodified theorem; the gauge-orbit part annihilates any function depending only on gauge-invariant observables. So expectation values of Wilson loops and other gauge-invariant operators are unaffected by $G$. The Lorentz-covariant choice $G[U]=\frac{1}{2\xi }(\partial \cdot A{)}^2$ for $U=e^{iaA}$ in the small-$a$ continuum limit recovers the conventional covariant gauge with no ghost insertion.

Theorem (Drummond-Horgan-Pendleton lattice-Langevin algorithm; 1983 ^{[Drummond-Horgan-Pendleton]}). The discrete-time stochastic Langevin update for $\mathrm{SU}(N)$ lattice gauge theory, with a single-step generator that preserves the Haar measure to leading order in $\Delta \tau$, samples the Wilson measure with $\Delta \tau$-bias of order $\Delta \tau$ and converges in distribution to the continuous-time Langevin SDE as $\Delta \tau \to 0$. Numerical simulations on $\mathrm{SU}(2)$ and $\mathrm{SU}(3)$ recover the conventional Monte Carlo plaquette expectation values within statistical errors at $\Delta \tau$ small enough for the bias to be subleading.

This was the first systematic Langevin implementation of lattice gauge theory and remains a reference benchmark. The update is the group-manifold analogue of Euler-Maruyama: starting from $U_n$, propose $U_{n+1}=\exp (i{T}^a(-\Delta \tau {\nabla }_a{S}_W+\sqrt{2\Delta \tau }{\xi }_n^a))\cdot U_n$ with ${\xi }_n^a$ standard Gaussian. The exponential map keeps the trajectory on the group, and the Stratonovich-style discretisation produces the correct Itô correction.

Theorem (Duane-Kennedy-Pendleton-Roweth hybrid Monte Carlo; 1987 ^[DKPR]). Augment the configuration space ${\mathbb{R}}^{|\Lambda |}$ with a fictitious momentum $\pi \in {\mathbb{R}}^{|\Lambda |}$ and the augmented Hamiltonian $H(\varphi ,\pi )=S(\varphi )+\frac{1}{2}|\pi {|}^2$. The HMC update consists of (i) Gibbs sampling of $\pi \sim \mathcal{N}(0,I)$, (ii) leapfrog integration of Hamilton's equations for $H$ over $K$ steps with size $\Delta \tau$, (iii) a Metropolis accept/reject step with acceptance probability $\min (1,e^{-\Delta H})$. The augmented Markov chain has $e^{-H(\varphi ,\pi )}/Z$ as its exact stationary distribution; marginalising over $\pi$ recovers $e^{-S(\varphi )}/Z$ exactly, with no order-$\Delta \tau$ bias.

HMC has become the dominant sampling algorithm for lattice QCD because it bypasses the order-$\Delta \tau$ bias of plain Langevin updates while maintaining the favourable scaling of stochastic methods over local Metropolis updates. The combination — exact target, no rejection-of-individual-coordinates inefficiency, long-trajectory exploration — is the central reason essentially all production lattice-QCD calculations since the early 1990s use HMC or its descendants (Rational HMC, mass preconditioning, multi-time-scale integration).

Theorem (Kennedy-Bowler pseudofermion Langevin; 1985 ^{[Kennedy-Bowler]}). Lattice gauge theory with dynamical fermions has an action $S_W[U]+\log \det M[U]$, with $M[U]$ the lattice Dirac operator. The fermion determinant is replaced by a Gaussian pseudofermion integral $\det M^{†}M=\int e^{-{\varphi }^{†}(M^{†}M{)}^{-1}\varphi }d{\varphi }^{†}d\varphi$. The Langevin update on $(U,{\varphi }^{†},\varphi )$ with the joint action samples the dynamical-fermion measure, with $(M^{†}M{)}^{-1}\varphi$ computed by iterative inversion at each step.

This was the workhorse algorithm for lattice QCD with dynamical fermions in the second half of the 1980s, before HMC overtook it. The Langevin-style updates allowed the first systematic studies of unquenched lattice QCD; the pseudofermion trick is what made it computationally feasible to handle the non-local fermion determinant.

Theorem (complex Langevin; Parisi 1983 ^{[Parisi-1983]}, Aarts-Seiler-Stamatescu 2010 ^[ASS-2010]). For an action with a complex or signed contribution (e.g., finite chemical potential, real-time path integral), naively the integrand $e^{-S}/Z$ is not a probability density, and conventional Langevin sampling does not apply. Complex Langevin extends the framework by complexifying the configuration space to ${\mathbb{C}}^{|\Lambda |}$ and integrating the Langevin SDE with complex drift. When the drift on the complexified manifold is holomorphic with sufficient decay at infinity, the resulting expectation values agree with the analytic continuation of the original path integral.

The convergence of complex Langevin in the general case is controversial. Aarts-Seiler-Stamatescu 2010 ^[ASS-2010] identifies the holomorphic-drift criterion and the requirement of exponentially-decaying measures on the complexified configuration space; specific failure modes include drift singularities, slow decay at infinity, and "incorrect convergence" to wrong-but-stable distributions. The framework remains an active research topic in lattice field theory at finite chemical potential and real-time.

Theorem (connection to Markov-chain Monte Carlo theory). The discretised Langevin update is a Markov chain on ${\mathbb{R}}^{|\Lambda |}$ with transition kernel $$ P_{\Delta\tau}(\phi, d\phi') = \mathcal{N}(\phi - \Delta\tau \nabla S(\phi),, 2 \Delta\tau \cdot I)(d\phi'). $$ The chain does not satisfy detailed balance with respect to the Gibbs measure at finite $\Delta \tau$ — the bias is order $\Delta \tau$ — but it does in the $\Delta \tau \to 0$ limit. The Metropolis-Hastings correction (Roberts-Tweedie 1996; the "MALA" algorithm: Metropolis-adjusted Langevin algorithm) accepts proposals with probability $\min (1,q({\varphi }^{\prime },\varphi )/q(\varphi ,{\varphi }^{\prime })\cdot e^{-S({\varphi }^{\prime })+S(\varphi )})$, restoring exact detailed balance. The HMC scheme of Duane-Kennedy-Pendleton-Roweth 1987 is the higher-order symplectic generalisation.

The connection identifies Langevin sampling as a particular case of MCMC theory, with the additional structure of continuous-time relaxation. The unifying perspective gives the modern toolkit of Wasserstein-gradient-flow methods (Jordan-Kinderlehrer-Otto 1998), the entropy-method convergence analysis (Bakry-Émery, Otto-Villani), and the high-dimensional concentration estimates used in machine-learning Langevin samplers (Langevin dynamics MCMC, score-based generative models, diffusion models).

Synthesis. Langevin updates and lattice numerics put together the foundational reason stochastic quantisation is a viable construction route for Euclidean field theory: a finite-dimensional Itô SDE on ${\mathbb{R}}^{|\Lambda |}$ with drift the negative gradient of the action and white noise has the Euclidean Gibbs measure as its unique exponentially-attracting stationary distribution, and this is exactly the algebraic content of the Fokker-Planck identity $L^{\ast }(e^{-S}/Z)=0$ derived in 08.10.02. The central insight is exactly that the Langevin-Fokker-Planck pair encodes a reversible Markov dynamics whose detailed-balance distribution is the path-integral measure, putting these together with Euler-Maruyama discretisation and Metropolis accept/reject corrections to produce a computer-implementable sampler with controlled bias. The bridge is the recognition that this framework generalises in three directions: to gauge theories on compact group manifolds (Drummond-Horgan-Pendleton 1983), where the Haar measure replaces Lebesgue measure and the gauge-invariant Wilson measure emerges as the equilibrium without conventional gauge fixing; to fermionic theories (Kennedy-Bowler 1985), where pseudofermion auxiliary fields convert the determinant into an integral that the Langevin sampler can handle; and to signed-measure regimes (Parisi 1983), where complexification of the configuration space gives a partial — and controversial — extension to actions with imaginary parts.

The bridge between Langevin sampling and the Wightman framework of 08.10.07 is Osterwalder-Schrader: the Schwinger functions sampled by Langevin updates on the Euclidean side reconstruct the Wightman functions on the Minkowski side by analytic continuation, identifying stochastic-quantisation expectation values with relativistic vacuum expectations of operator products. This is exactly the same reflection-positivity / detailed-balance identification that appears again in the Glimm-Jaffe constructive program, where the Langevin / Markov-chain framework supplies the existence of the measures that the Osterwalder-Schrader axioms require.

The synthesis runs in the opposite direction as well. Hybrid Monte Carlo (Duane-Kennedy-Pendleton-Roweth 1987) trades the order-$\Delta \tau$ bias of plain Langevin for exact sampling at the cost of a Metropolis accept/reject; the central insight is that augmenting the configuration space with a fictitious momentum produces a Hamiltonian flow whose volume-preserving symplectic integrator gives unbiased proposals. This is dual to the Langevin gradient-flow picture in the same way Hamiltonian dynamics is dual to gradient descent: the leapfrog integrator generalises Euler-Maruyama by adding ballistic exploration along the level sets of the augmented Hamiltonian. Putting these together with mass preconditioning (Hasenbusch 2001), rational HMC (Clark-Kennedy 2007), and multi-time-scale integration produces the modern production toolkit for lattice QCD, which has computed the proton mass, the QCD spectrum, hadronic structure functions, and electroweak matrix elements to percent-level precision from first principles. The foundational reason the toolkit works is that Langevin-style stochastic dynamics on the lattice path integral converges to the right equilibrium; the central insight is that the equilibrium is gauge-invariant from the outset, so no Faddeev-Popov machinery is needed to extract physics. The bridge between lattice numerics and machine learning is the recognition that score-based generative models, denoising diffusion models, and Langevin-dynamics MCMC samplers in modern probabilistic ML are all instances of the Parisi-Wu construction with a learned drift; the central insight that identifies stochastic quantisation with score-based generation is that both produce samples from a target measure as the long-time equilibrium of a noisy gradient flow.

Full proof set Master

Proposition (detailed balance for the Langevin SDE). The continuous-time Langevin SDE $d\varphi =-\nabla S(\varphi )d\tau +\sqrt{2}dW$ on ${\mathbb{R}}^{|\Lambda |}$ satisfies detailed balance with respect to the Gibbs measure ${\mu }_S(d\varphi )=e^{-S(\varphi )}/Zd\varphi$. Equivalently, the generator $L=\Delta -\nabla S\cdot \nabla$ is symmetric on $L^2({\mathbb{R}}^{|\Lambda |},{\mu }_S)$.

Proof. For $f,g\in C_c^{\infty }({\mathbb{R}}^{|\Lambda |})$, compute the bilinear form $⟨Lf,g{⟩}_{{\mu }_S}=\int (Lf)(\varphi )g(\varphi )e^{-S(\varphi )}/Zd\varphi$. Substituting $Lf=\Delta f-\nabla S\cdot \nabla f$ and integrating the Laplacian by parts in ${\mathbb{R}}^{|\Lambda |}$: $$ \int (\Delta f)(\phi), g(\phi), e^{-S(\phi)}, d\phi = -\int \nabla f \cdot \nabla(g, e^{-S}), d\phi = -\int \nabla f \cdot \nabla g, e^{-S}, d\phi + \int \nabla f \cdot \nabla S, g, e^{-S}, d\phi. $$ The boundary term at infinity vanishes by the compact support of $f,g$. Adding the explicit $-\nabla S\cdot \nabla f\cdot g{e}^{-S}$ term from $L$: $$ \langle Lf, g \rangle_{\mu_S} = -\frac{1}{Z} \int \nabla f \cdot \nabla g, e^{-S}, d\phi. $$ The right-hand side is symmetric under $f\leftrightarrow g$, so $⟨Lf,g{⟩}_{{\mu }_S}=⟨f,Lg{⟩}_{{\mu }_S}$. Symmetry of $L$ on $L^2({\mu }_S)$ is equivalent to detailed balance of the SDE: the transition density $p_{\tau }(\varphi ,{\varphi }^{\prime })$ satisfies $p_{\tau }(\varphi ,{\varphi }^{\prime })e^{-S(\varphi )}=p_{\tau }({\varphi }^{\prime },\varphi )e^{-S({\varphi }^{\prime })}$ for every $\tau \ge 0$ and almost every $\varphi ,{\varphi }^{\prime }$. $\square$

Proposition (exponential ergodicity under Bakry-Émery). If the action $S$ satisfies ${\nabla }^{2}S(\varphi )\ge \rho I$ pointwise for some $\rho >0$ (in the sense of symmetric matrices), the Langevin SDE converges exponentially to the Gibbs measure in relative entropy: $$ H(p(\cdot, \tau), |, \mu_S) \leq e^{-2 \rho \tau}, H(p(\cdot, 0), |, \mu_S), $$ where $H(p\Vert {\mu }_S)=\int p\log (p\cdot Z{e}^S)d\varphi$ is the Kullback-Leibler divergence.

Proof. Differentiate the relative entropy along the Fokker-Planck flow ${\partial }_{\tau }p=\Delta p+\nabla \cdot (p\nabla S)$: $$ \frac{d}{d\tau} H(p, |, \mu_S) = \int \partial_\tau p \cdot \log(p \cdot e^S \cdot Z), d\phi = -\int \frac{|\nabla(p, e^S)|^2}{p, e^S}, e^{-S}, d\phi = -I(p, |, \mu_S), $$ where $I(p\Vert {\mu }_S)$ is the Fisher information of $p$ relative to ${\mu }_S$. The Bakry-Émery curvature-dimension criterion ${\nabla }^{2}S\ge \rho I$ implies the log-Sobolev inequality $H(p\Vert {\mu }_S)\le \frac{1}{2\rho }I(p\Vert {\mu }_S)$ for all probability densities $p$. So $$ \frac{d}{d\tau} H(p | \mu_S) \leq -2\rho, H(p | \mu_S), $$ and Gronwall's inequality gives $H(p(\tau )\Vert {\mu }_S)\le e^{-2\rho \tau }H(p(0)\Vert {\mu }_S)$. The Csiszár-Kullback-Pinsker inequality bounds $\Vert p(\tau )-{\mu }_S{\Vert }_{L^1}^2\le 2H(p(\tau )\Vert {\mu }_S)\to 0$, giving convergence in total variation. $\square$

Proposition (Euler-Maruyama weak convergence). Let $b:{\mathbb{R}}^{|\Lambda |}\to {\mathbb{R}}^{|\Lambda |}$ be globally Lipschitz and bounded, $\sigma >0$ constant. For every $T>0$ and every $C^4$ bounded $f:{\mathbb{R}}^{|\Lambda |}\to \mathbb{R}$, $$ \big| \mathbb{E}[f(\phi_N^{\Delta\tau})] - \mathbb{E}[f(\phi(T))] \big| \leq C(f, T, b), \Delta\tau, $$ where $\varphi (\tau )$ solves $d\varphi =b(\varphi )d\tau +\sigma dW$ and ${\varphi }_n^{\Delta \tau }$ is the Euler-Maruyama discretisation with $N\Delta \tau =T$.

Proof. Define $u(\varphi ,s)=\mathbb{E}[f(\varphi (T))\mid \varphi (T-s)=\varphi ]$, which satisfies the backward Kolmogorov equation ${\partial }_{s}u=Lu$ with $L=b\cdot \nabla +\frac{{\sigma }^2}{2}\Delta$ and terminal condition $u(\varphi ,0)=f(\varphi )$. Standard regularity for parabolic PDEs gives $u\in C^4$ uniformly bounded on $[0,T]$ with bounds depending on $f,T,b$.

The error decomposes telescopically: $$ \mathbb{E}[f(\phi_N^{\Delta\tau})] - u(\phi_0, T) = \sum_{n=0}^{N-1} \big( \mathbb{E}[u(\phi_{n+1}^{\Delta\tau}, T - (n+1)\Delta\tau)] - \mathbb{E}[u(\phi_n^{\Delta\tau}, T - n \Delta\tau)] \big). $$ Each summand is the single-step error at step $n$. Taylor-expand $u({\varphi }_{n+1}^{\Delta \tau },T-(n+1)\Delta \tau )$ around $({\varphi }_n^{\Delta \tau },T-n\Delta \tau )$ to fourth order. With ${\varphi }_{n+1}^{\Delta \tau }={\varphi }_n^{\Delta \tau }+b({\varphi }_n^{\Delta \tau })\Delta \tau +\sigma \sqrt{\Delta \tau }{\xi }_n$, ${\xi }_n\sim \mathcal{N}(0,I)$: $$ u(\phi_{n+1}^{\Delta\tau}, s - \Delta\tau) = u(\phi_n^{\Delta\tau}, s) + (\partial_s u + b \cdot \nabla u + \tfrac{\sigma^2}{2} \Delta u)(-\Delta\tau) + O(\Delta\tau^2) \cdot \xi_n^{\otimes \geq 3}. $$ Take the expectation: the linear-in-${\xi }_n$ term has mean zero, the quadratic-in-${\xi }_n$ term gives the ${\sigma }^2/2\cdot \Delta u$ contribution that combines with $b\cdot \nabla u$ and ${\partial }_{s}u$ to recover $Lu-{\partial }_{s}u=0$ at leading order, and the remainder is $O(\Delta {\tau }^2)$ from the cubic and quartic Taylor terms. Summing $N$ steps gives total error $O(N\Delta {\tau }^2)=O(T\Delta \tau )$. $\square$

Proposition (HMC detailed balance). The hybrid Monte Carlo update (Gaussian momentum refresh, leapfrog integration of Hamiltonian $H(\varphi ,\pi )=S(\varphi )+\frac{1}{2}|\pi {|}^2$ for $K$ steps with step size $\Delta \tau$, Metropolis accept/reject with probability $\min (1,e^{-\Delta H})$) satisfies detailed balance with respect to the augmented measure ${\mu }_H(d\varphi ,d\pi )=e^{-H(\varphi ,\pi )}/Z^{\prime }d\varphi d\pi$, hence samples its marginal $e^{-S(\varphi )}/Zd\varphi$ exactly.

Proof. The augmented update consists of three steps: (a) sample $\pi \sim \mathcal{N}(0,I)$ — this is a Gibbs step on $\pi$ conditional on $\varphi$, which preserves ${\mu }_H$ by direct verification: the conditional law of $\pi$ given $\varphi$ under ${\mu }_H$ is exactly $\mathcal{N}(0,I)$. (b) Integrate Hamilton's equations $\dot{\varphi }=\pi$, $\dot{\pi }=-\nabla S(\varphi )$ with the leapfrog scheme: $$ \pi_{1/2} = \pi - \tfrac{\Delta\tau}{2} \nabla S(\phi),\quad \phi' = \phi + \Delta\tau, \pi_{1/2},\quad \pi' = \pi_{1/2} - \tfrac{\Delta\tau}{2} \nabla S(\phi'), $$ iterated $K$ times. The leapfrog map $(\varphi ,\pi )\mapsto ({\varphi }^{\prime },{\pi }^{\prime })$ is symplectic (preserves the canonical symplectic form $d\varphi \wedge d\pi$) and time-reversible ($({\varphi }^{\prime },{\pi }^{\prime })\mapsto (\varphi ,\pi )$ under $\pi \to -\pi$). (c) Accept the proposal $({\varphi }^{\prime },{\pi }^{\prime })$ with probability $\min (1,e^{-\Delta H})$ where $\Delta H=H({\varphi }^{\prime },{\pi }^{\prime })-H(\varphi ,\pi )$ is the integrator energy drift.

Detailed balance: for any test sets $A,B\subset {\mathbb{R}}^{|\Lambda |}\times {\mathbb{R}}^{|\Lambda |}$, $$ \int_A \mu_H(d\phi, d\pi) \cdot P((\phi, \pi), B) = \int_B \mu_H(d\phi', d\pi') \cdot P((\phi', \pi'), A), $$ where $P$ is the HMC kernel. Symplecticity of leapfrog gives unit Jacobian, time-reversibility gives the matching backward proposal, and the Metropolis accept/reject ratio is calibrated to make the two sides equal. The technical content is the standard Metropolis-Hastings construction with the proposal density determined by the leapfrog map. Marginalising over $\pi$ leaves $e^{-S(\varphi )}/Zd\varphi$, the target measure. $\square$

Proposition (lattice harmonic-oscillator equilibrium recovered). Let $\Lambda$ be a 1D lattice of $N$ sites with periodic boundary conditions, $\varphi \in {\mathbb{R}}^N$, and the action $S[\varphi ]=\frac{1}{2}{\omega }^2{\sum }_x\varphi (x{)}^2$. The Langevin SDE has equilibrium variance $\mathbb{E}[\varphi (x{)}^2]=1/{\omega }^2$ at every site, recovering the exact path-integral result $T/{\omega }^2=1/{\omega }^2$ (with $T=1$).

Proof. The action is a sum of independent harmonic oscillators at each site, so the equilibrium measure factors as ${\prod }_x{e}^{-{\omega }^2\varphi (x{)}^2/2}/\sqrt{2\pi /{\omega }^2}$. Each factor is a Gaussian with variance $1/{\omega }^2$. So $\mathbb{E}[\varphi (x{)}^2]=1/{\omega }^2$ at every site, independent of $x$. The lattice Langevin samples reproduce this value to within Monte Carlo statistical error. $\square$

Proposition (Gaussian-free-field two-point function from Langevin samples). On the 2D torus lattice $\Lambda =(\mathbb{Z}/L\mathbb{Z}{)}^2$ with action $S[\varphi ]=\frac{1}{2}{\sum }_x[(\nabla \varphi {)}^2+m^2{\varphi }^2]$, the Langevin equilibrium two-point function is the inverse lattice Laplacian: $$ \langle \phi(0), \phi(x) \rangle_{\mu_S} = (-\Delta_\Lambda + m^2)^{-1}(0, x), $$ recovering the Gaussian-free-field propagator of 08.06.01.

Proof. The action is a positive-definite quadratic form $S[\varphi ]=\frac{1}{2}{\varphi }^{T}H\varphi$ with $H=-{\Delta }_{\Lambda }+m^2$ the lattice Helmholtz operator. The Gibbs measure is the multivariate Gaussian $\mathcal{N}(0,H^{-1})$ on ${\mathbb{R}}^{L^2}$, so the equilibrium covariance is $\mathbb{E}[\varphi (x)\varphi (y)]=H^{-1}(x,y)=(-{\Delta }_{\Lambda }+m^2{)}^{-1}(x,y)$. In Fourier transform on the dual lattice $\hat{\Lambda }$ with momenta $k_{\mu }=2\pi n_{\mu }/L$, $n_{\mu }\in \{0,1,\dots ,L-1\}$: $$ \langle \phi(0) \phi(x) \rangle = \frac{1}{L^2} \sum_{k \in \hat{\Lambda}} \frac{e^{i k \cdot x}}{\hat{k}^2 + m^2}, \quad \hat{k}^2 = 4 \sum_\mu \sin^2(k_\mu/2). $$ Langevin samples ${\varphi }_n$ give the empirical estimator $\hat{C}(0,x)=(N-N_0{)}^{-1}{\sum }_n{\varphi }_n(0){\varphi }_n(x)\to ⟨\varphi (0)\varphi (x)⟩$ by ergodicity, with statistical error scaling as $1/\sqrt{N-N_0}$ corrected by the autocorrelation time of the Langevin chain. $\square$

Proposition (HMC dominance over plain Langevin in lattice QCD). Hybrid Monte Carlo achieves $O(1)$ acceptance with $\Delta \tau$ chosen so that the integrator energy error per trajectory is $O(1)$, requiring $\Delta \tau =O(N^{-1/4})$ in $N$ degrees of freedom. Plain Langevin requires $\Delta \tau =O(N^{-1/3})$ to keep the bias subleading. So HMC has favourable computational complexity $O(N^{5/4})$ per independent sample versus $O(N^{4/3})$ for plain Langevin, leading to a speedup factor of $N^{1/12}$ that grows with system size.

Proof sketch. Acceptance of HMC requires $|\Delta H|=O(1)$ over a trajectory. The leapfrog integrator has local energy error $O(\Delta {\tau }^3)$ per step (symplectic order-2 integrator), and $K=T/\Delta \tau$ steps in a trajectory give total energy error $O(K\Delta {\tau }^3)=O(T\Delta {\tau }^2)$. For $N$ degrees of freedom the energy fluctuates as $\sqrt{N}\Delta {\tau }^2$, so acceptance requires $\Delta \tau =O(N^{-1/4})$. Plain Langevin has order-$\Delta \tau$ stationary-distribution bias, with the constant scaling as $\sqrt{N}$, so the bias requires $\Delta \tau =O(N^{-1/2})$ for $O(1)$ accuracy on global observables (or $O(N^{-1/3})$ when restricted to local averages). The arithmetic cost per step is $O(N)$ in both schemes, giving total cost $N\cdot 1/\Delta \tau \cdot T_{\mathrm{auto}}$ where $T_{\mathrm{auto}}$ is the autocorrelation time. The favourable HMC scaling in $\Delta \tau$ dominates the autocorrelation overhead in practice. Creutz 1988 ^{[Damgaard-Huffel]} and Kennedy et al. 1990s subsequent work confirm this scaling in lattice QCD simulations. $\square$

Connections Master

• Fokker-Planck equation and equilibrium distribution 08.10.02. The Langevin SDE $d\varphi =-\nabla Sd\tau +\sqrt{2}dW$ is the Itô-form prerequisite for which 08.10.02 supplies the equivalent Fokker-Planck description: $L^{\ast }{p}_{\mathrm{eq}}=0$ for $p_{\mathrm{eq}}=e^{-S}/Z$, with the spectral-gap argument identifying $e^{-S}/Z$ as the unique stationary solution. The present unit specialises that infinite-dimensional Fokker-Planck framework to a finite lattice ${\mathbb{R}}^{|\Lambda |}$ and discretises the SDE for computational implementation.

• Path integral formulation of statistical mechanics 08.07.01. The path-integral measure $d{\mu }_S(\varphi )=e^{-S}/Zd\varphi$ on lattice field configurations is what the Langevin updates sample. The present unit gives a concrete computer-implementable algorithm for drawing samples from the measure that 08.07.01 introduces as the formal foundation of Euclidean field theory.

• Gaussian field theory and free boson 08.06.01. The Gaussian-free-field measure with covariance $C=(-{\Delta }_{\Lambda }+m^2{)}^{-1}$ is the Langevin equilibrium for the quadratic action $S[\varphi ]=\frac{1}{2}{\varphi }^T(-{\Delta }_{\Lambda }+m^2)\varphi$. Langevin samples produce the propagator $⟨\varphi (0)\varphi (x)⟩=C(0,x)$ numerically, with the worked-example on the 2D torus serving as the calibration benchmark.

• Wick rotation 08.09.01. The fictitious stochastic time $\tau$ of the Langevin SDE is unrelated to the physical Wick-rotated time $t_E$ of the Euclidean field theory: $\tau$ is an extra dimension over which the dynamics relaxes to equilibrium, while $t_E$ is one of the lattice spatial dimensions in the Euclidean signature obtained from 08.09.01. The Wick-rotated action is what enters the Langevin drift; the stochastic time is the extra parameter Parisi-Wu introduce.

• Wilson's lattice gauge theory 08.08.01. The Drummond-Horgan-Pendleton extension of Langevin updates to $\mathrm{SU}(N)$-valued link variables is the gauge-theoretic application of the present framework. The Wilson plaquette action $S_W[U]$ supplies the drift; the gauge-invariant equilibrium measure is the Wilson measure of 08.08.01, sampled without conventional Faddeev-Popov gauge fixing.

• Wightman axioms 08.10.07. Schwinger functions sampled by lattice Langevin updates reconstruct Wightman functions on the Minkowski side by Osterwalder-Schrader analytic continuation. The present stochastic-sampling framework is the constructive-QFT input that Glimm-Jaffe and successors use to verify the Wightman axioms in two and three spacetime dimensions; the Langevin equilibrium is precisely the measure whose Schwinger functions Osterwalder-Schrader reconstruction turns into a Wightman theory.

Historical & philosophical context Master

The Langevin equation originated in Paul Langevin's 1908 note Sur la théorie du mouvement brownien (C. R. Acad. Sci. Paris 146, 530) as a phenomenological model for Brownian motion, with deterministic friction and a stochastic force. The mathematical foundations were assembled in the following decades: Norbert Wiener's 1923 construction of the Wiener measure on continuous paths, Kiyosi Itô's 1944 stochastic integral ^[Itô-1944] and 1951 chain rule ^[Itô-1951] (Proc. Imp. Acad. Tokyo 20, 519 and Nagoya Math. J. 3, 55), and Andrey Kolmogorov's 1931 axiomatisation of continuous Markov processes via the forward (Fokker-Planck) and backward equations. The discretised Euler-Maruyama scheme is due to Gisiro Maruyama's 1955 paper Continuous Markov processes and stochastic equations (Rend. Circ. Mat. Palermo 4, 48) ^{[Maruyama-1955]}, establishing weak convergence of the discrete scheme to the continuous SDE.

The application to quantum field theory is the work of Giorgio Parisi and Yong-Shi Wu, Perturbation theory without gauge fixing (Sci. Sin. 24, 483 (1981)) ^[Parisi-Wu]. Their proposal: the Euclidean path-integral measure $e^{-S[\varphi ]}/Z$ on field configurations can be obtained as the equilibrium distribution of a Langevin process in a fictitious "fifth time" $\tau$, with the configuration field $\varphi (x,\tau )$ evolving according to ${\partial }_{\tau }\varphi =-\delta S/\delta \varphi +\eta$. The proposal had three immediate appeals: it gave a mathematically clean construction of the Euclidean measure as the long-time limit of a stochastic process (sidestepping the measure-theoretic complications of defining an "infinite-dimensional Gaussian times Boltzmann factor"); it bypassed conventional gauge fixing for gauge theories, because the equilibrium measure is automatically gauge-invariant; and it suggested a computational route based on time-stepping the SDE on a lattice. Daniel Zwanziger's 1981 paper Covariant quantization of gauge fields without Gribov ambiguity (Nucl. Phys. B 192, 259) ^[Zwanziger] introduced stochastic gauge fixing as a Lorentz-covariant alternative that does not require Faddeev-Popov ghosts and does not suffer the Gribov copies of conventional gauge fixing.

The numerical lattice implementation was carried out by I. T. Drummond, R. R. Horgan, and B. J. Pendleton in Improved gauge fixing and Langevin algorithms (Nucl. Phys. B 220, 119 (1983)) ^{[Drummond-Horgan-Pendleton]}, who implemented discretised Langevin updates on $\mathrm{SU}(2)$ and $\mathrm{SU}(3)$ link variables and compared against conventional Metropolis Monte Carlo. A. D. Kennedy and K. C. Bowler's 1985 paper Updating fermions with Langevin algorithm (Nucl. Phys. B 245, 217) ^{[Kennedy-Bowler]} extended the framework to dynamical fermions via pseudofermion auxiliary fields, making lattice QCD with sea quarks computationally feasible for the first time. The 1987 hybrid Monte Carlo paper by Simon Duane, Anthony Kennedy, Brian Pendleton, and Duncan Roweth (Phys. Lett. B 195, 216) ^[DKPR] combined the Langevin-style stochastic exploration with a Metropolis accept/reject to eliminate the order-$\Delta \tau$ bias, and HMC has been the standard sampling algorithm for lattice QCD ever since. The canonical review of the stochastic-quantisation framework is Damgaard and Hüffel's 1987 Phys. Rep. 152, 227 article ^{[Damgaard-Huffel]}.

Parisi's 1983 paper On complex probabilities (Phys. Lett. B 131, 393) ^{[Parisi-1983]} proposed complex Langevin as an extension to actions with imaginary parts — the sign-problem regime at finite chemical potential or in real-time path integrals. Convergence in the general case proved controversial; Gert Aarts, Erhard Seiler, and Ion-Olimpiu Stamatescu's 2010 paper The complex Langevin method: when can it be trusted? (Phys. Rev. D 81, 054508) ^[ASS-2010] identified the holomorphic-drift criterion and exhibited specific failure modes. The framework remains active in finite-density lattice QCD and real-time field theory research.

Bibliography Master

@article{ParisiWu1981,
  author  = {Parisi, Giorgio and Wu, Yong-Shi},
  title   = {Perturbation theory without gauge fixing},
  journal = {Scientia Sinica},
  volume  = {24},
  year    = {1981},
  pages   = {483--496}
}

@article{DamgaardHuffel1987,
  author  = {Damgaard, Poul H. and H{\"u}ffel, Helmuth},
  title   = {Stochastic quantization},
  journal = {Physics Reports},
  volume  = {152},
  number  = {5-6},
  year    = {1987},
  pages   = {227--398}
}

@article{Zwanziger1981,
  author  = {Zwanziger, Daniel},
  title   = {Covariant quantization of gauge fields without {G}ribov ambiguity},
  journal = {Nuclear Physics B},
  volume  = {192},
  year    = {1981},
  pages   = {259--269}
}

@article{DrummondHorganPendleton1983,
  author  = {Drummond, I. T. and Horgan, R. R. and Pendleton, B. J.},
  title   = {Improved gauge fixing and {L}angevin algorithms},
  journal = {Nuclear Physics B},
  volume  = {220},
  year    = {1983},
  pages   = {119--136}
}

@article{KennedyBowler1985,
  author  = {Kennedy, A. D. and Bowler, K. C.},
  title   = {Updating fermions with {L}angevin algorithm},
  journal = {Nuclear Physics B},
  volume  = {245},
  year    = {1985},
  pages   = {217--232}
}

@article{Parisi1983Complex,
  author  = {Parisi, Giorgio},
  title   = {On complex probabilities},
  journal = {Physics Letters B},
  volume  = {131},
  year    = {1983},
  pages   = {393--395}
}

@article{DKPR1987,
  author  = {Duane, Simon and Kennedy, A. D. and Pendleton, Brian J. and Roweth, Duncan},
  title   = {Hybrid {M}onte {C}arlo},
  journal = {Physics Letters B},
  volume  = {195},
  year    = {1987},
  pages   = {216--222}
}

@article{Ito1944,
  author  = {It{\^o}, Kiyosi},
  title   = {Stochastic integral},
  journal = {Proceedings of the Imperial Academy of Tokyo},
  volume  = {20},
  year    = {1944},
  pages   = {519--524}
}

@article{Ito1951,
  author  = {It{\^o}, Kiyosi},
  title   = {On a formula concerning stochastic differentials},
  journal = {Nagoya Mathematical Journal},
  volume  = {3},
  year    = {1951},
  pages   = {55--65}
}

@article{Maruyama1955,
  author  = {Maruyama, Gisiro},
  title   = {Continuous {M}arkov processes and stochastic equations},
  journal = {Rendiconti del Circolo Matematico di Palermo},
  volume  = {4},
  year    = {1955},
  pages   = {48--90}
}

@book{KloedenPlaten1992,
  author    = {Kloeden, Peter E. and Platen, Eckhard},
  title     = {Numerical Solution of Stochastic Differential Equations},
  publisher = {Springer-Verlag},
  series    = {Stochastic Modelling and Applied Probability},
  volume    = {23},
  year      = {1992}
}

@article{AartsSeilerStamatescu2010,
  author  = {Aarts, Gert and Seiler, Erhard and Stamatescu, Ion-Olimpiu},
  title   = {The complex {L}angevin method: when can it be trusted?},
  journal = {Physical Review D},
  volume  = {81},
  year    = {2010},
  pages   = {054508}
}

@article{Langevin1908,
  author  = {Langevin, Paul},
  title   = {Sur la th{\'e}orie du mouvement brownien},
  journal = {Comptes Rendus de l'Acad{\'e}mie des Sciences Paris},
  volume  = {146},
  year    = {1908},
  pages   = {530--533}
}

@article{Kolmogorov1931,
  author  = {Kolmogorov, Andrey N.},
  title   = {{\"U}ber die analytischen {M}ethoden in der {W}ahrscheinlichkeitsrechnung},
  journal = {Mathematische Annalen},
  volume  = {104},
  year    = {1931},
  pages   = {415--458}
}

@article{JordanKinderlehrerOtto1998,
  author  = {Jordan, Richard and Kinderlehrer, David and Otto, Felix},
  title   = {The variational formulation of the {F}okker-{P}lanck equation},
  journal = {SIAM Journal on Mathematical Analysis},
  volume  = {29},
  year    = {1998},
  pages   = {1--17}
}

@article{BakryEmery1985,
  author  = {Bakry, Dominique and {\'E}mery, Michel},
  title   = {Diffusions hypercontractives},
  journal = {S{\'e}minaire de Probabilit{\'e}s XIX, Lecture Notes in Math.},
  volume  = {1123},
  year    = {1985},
  pages   = {177--206}
}