Introduction: Another Bridge Connecting Discrete Jumps to Continuous Descriptions¶

In Lecture 36, a systematic framework—the coherent-state path integral (CSP)—was successfully constructed. This framework, built upon operator algebra as its foundation, precisely transforms master equations describing discrete particle birth-death processes into a path integral paradigm where field-theoretic tools can be applied, through the introduction of Poisson representations and creation-annihilation operators. This algebraic approach provides a solid theoretical foundation for handling stochastic chemical kinetics.

This lecture will explore another path connecting discrete and continuous descriptions from a different, more physically intuitive perspective. The core question remains unchanged: How can we systematically derive an equivalent continuous description, such as the Fokker-Planck equation or the corresponding Langevin equation, starting from a master equation describing microscopic, discrete "jump" events (e.g., particle number changing from $n$ to $n + \Delta n$)?

To achieve this goal, the key mathematical tool to be introduced in this lecture is the Kramers-Moyal expansion. It is a formal Taylor series expansion of the master equation that can transform an integro-differential master equation into an infinite-order partial differential equation. The first two terms of this expansion will naturally correspond to the drift and diffusion terms in the Fokker-Planck equation, thus providing rigorous mathematical justification for the transition from discrete master equations to continuous Langevin dynamics under the "low noise limit".

Therefore, this lecture aims to establish a powerful comparison: on one side is the coherent-state path integral (CSP) from Lecture 36, derived from operator algebra, and on the other side is the Kramers-Moyal path integral (KMPI) to be constructed in this lecture, derived from probabilistic expansions. The ultimate goal is to demonstrate how these two fundamentally different theoretical perspectives converge to the same physical description under appropriate limits. This duality between different formalisms reveals the profound unity inherent in stochastic process theory. The CSP method is extremely powerful in field-theoretic analysis (such as renormalization group), while the KMPI method provides a more direct bridge connecting physical concepts like drift and diffusion. Understanding this duality is key to mastering the content of this lecture.

1. Kramers-Moyal Expansion: From Master Equation to Partial Differential Equation¶

In previous lectures, the master equation has been established as the foundation for describing stochastic processes with discrete states (such as particle numbers). However, its integro-differential form and discrete variables present challenges for directly applying field-theoretic analytical tools. The goal of this lecture is to build a bridge from the discrete world of master equations to continuous field-theoretic descriptions. The first crucial step is to transform the master equation into a more analyzable partial differential equation form through the Kramers-Moyal expansion.

The Kramers-Moyal expansion was systematically established and developed by Dutch physicist Hendrik Anthony Kramers and later Josef L. Moyal in the 1940s (Kramers, 1940; Moyal, 1949), with its historical roots tracing back to early statistical physics research on the evolution of stochastic processes (such as Brownian motion). The physical essence of this expansion is to provide a precise mathematical framework for describing the time evolution of probability density functions for arbitrary Markov stochastic processes (including non-Gaussian, non-diffusive processes): it expresses the time rate of change of the probability distribution function as an infinite-order partial differential equation of all its moments (i.e., Kramers-Moyal coefficients), namely

\[\partial_t p(x, t) = \sum_{n=1}^{\infty} (-\partial_x)^n \left[ D^{(n)}(x, t) p(x, t) \right]\]

where the coefficients $D^{(n)}$ are directly related to the moments of stochastic increments.

In practical applications, the Kramers-Moyal expansion serves as the theoretical foundation for deriving the Fokker-Planck equation (when the expansion is truncated at second order) and is widely used to extract drift and diffusion terms of dynamics (i.e., $D^{(1)}$ and $D^{(2)}$) from experimental or simulation data, thus providing key tools for identifying underlying stochastic dynamic laws in fields such as neuroscience (analyzing neuronal firing), fluid turbulence, financial time series analysis, and single-molecule biophysics.

1.1 General Master Equation¶

Consider a system described by a single random variable $n$ (e.g., particle number). The system's state changes through a series of instantaneous "jump" events, each jump changing the particle number from $n$ to $n+\Delta n$. For a system involving only one type of jump process (with step size $\Delta n$), the time evolution of its probability distribution $p(n,t)$ is described by the following master equation:

\[ \partial_t p(n,t) = \lambda(n - \Delta n) p(n - \Delta n, t) - \lambda(n) p(n,t) \]

This equation precisely describes the conservation flow of probability in discrete state space:

Gain Term: $\lambda(n - \Delta n) p(n - \Delta n, t)$ represents the total probability flow jumping from state $n-\Delta n$ to state $n$. $p(n - \Delta n, t)$ is the probability of the system being in the source state, while $\lambda(n - \Delta n)$ is the transition rate (propensity) from that source state.
Loss Term: $-\lambda(n) p(n,t)$ represents the total probability flow jumping out of state $n$. $\lambda(n)$ is the total rate of jumps occurring in the current state $n$.

Although this form is precise, its "difference" nature (i.e., simultaneously involving $p(n,t)$ and $p(n-\Delta n, t)$) makes it difficult to handle directly.

1.2 Taylor Series Expansion¶

The core idea of the Kramers-Moyal expansion is to perform a formal Taylor series expansion of the gain term in the master equation (a function of $n-\Delta n$) around the point $n$. Expanding the function $f(n) = \lambda(n)p(n,t)$ around point $n$, we can obtain the expression for $f(n-\Delta n)$:

\[ \lambda(n - \Delta n) p(n - \Delta n, t) = \sum_{\mu=0}^{\infty} \frac{(-\Delta n)^{\mu}}{\mu!} \frac{\partial^{\mu}}{\partial n^{\mu}} \left[ \lambda(n) p(n,t) \right] \]

Substituting this infinite series back into the master equation, a crucial simplification occurs:

When $\mu = 0$, the first term of the expansion is $\frac{(-\Delta n)^0}{0!} \frac{\partial^0}{\partial n^0}[\lambda(n)p(n,t)] = \lambda(n)p(n,t)$.
This term exactly cancels with the loss term $-\lambda(n)p(n,t)$ in the master equation.

Therefore, the master equation is precisely rewritten as an infinite-order partial differential equation starting from $\mu=1$, which is the Kramers-Moyal expansion:

\[ \partial_t p(n,t) = \sum_{\mu=1}^{\infty} \frac{(-\Delta n)^{\mu}}{\mu!} \frac{\partial^{\mu}}{\partial n^{\mu}} \left[ \lambda(n) p(n,t) \right] \]

This equation is mathematically equivalent to the original master equation, but its form has undergone a fundamental transformation. It converts a discrete difference equation into a partial differential equation that treats $n$ as a continuous variable and contains derivatives of all orders, paving the way for subsequent analysis.

1.3 Pawula's Theorem and Its Physical Significance¶

This infinite-order expansion raises a crucial question: Under what circumstances can this series be truncated? Pawula's theorem provides a non-trivial answer with profound physical meaning: the Kramers-Moyal expansion either has all higher-order terms strictly zero after second order ($\mu = 2$), or it must contain infinitely many non-zero terms. Truncation at third order or higher is not allowed.

This theorem clearly divides stochastic processes into two major categories:

Pure Diffusion Processes (Fokker-Planck Processes): If the expansion stops exactly at $\mu = 2$, then the master equation is strictly equivalent to a Fokker-Planck equation. This corresponds to stochastic processes where state variables can change continuously, such as Brownian motion. In this case, the Fokker-Planck equation is an exact description of the system dynamics.
Jump Processes: If the expansion does not stop after second order, then according to Pawula's theorem, it must have infinitely many terms. This means any finite-order truncation (e.g., retaining only up to second order) can only be an approximation. This corresponds to processes with inherent discrete jump characteristics, such as molecular number changes in chemical reactions or population dynamics.

Therefore, for any real physical system driven by discrete events (such as chemical reaction networks), simplifying it to a Fokker-Planck equation is essentially an approximation. The validity of this approximation depends on a core physical assumption: the step size $\Delta n$ of a single jump is very small relative to the macroscopic state of the system (such as total particle number $n$). When the system scale is large and jump events occur frequently, the cumulative effect of numerous tiny, discrete jumps appears macroscopically like a continuous diffusion process. This is precisely the essence of the so-called "low noise limit" or "system size expansion", and also the theoretical foundation for subsequently connecting master equations with Langevin equations.

2. Path Integral Representation of Jump Processes¶

The previous section successfully transformed the master equation from a discrete difference form to a continuous, infinite-order partial differential equation through the Kramers-Moyal expansion. Although this step is mathematically equivalent, its true purpose is to pave the way for constructing a more powerful theoretical framework—path integrals. The goal of this lecture is to transform this partial differential equation into path integral form, thus providing a global, "action"-based perspective on dynamics, which echoes the field-theoretic framework established for Langevin equations in Lectures 31 and 32.

2.1 Introduction of Response Field $\tilde{n}$¶

The first step in the transformation is to rewrite the Kramers-Moyal expansion partial differential equation as an integral form through a series of mathematical techniques. The core of this process is introducing an auxiliary response field $\tilde{n}$.

First, using the properties of the Dirac $\delta$ function, write the right-hand side of the partial differential equation (a function of $n$) as an integral form. This is an identity transformation aimed at introducing an integration variable $m$:

\[ \partial_t p(n,t) = \int dm\ \delta(m-n) \sum_{\mu=1}^{\infty} \frac{(-\Delta n)^{\mu}}{\mu!} \frac{\partial^{\mu}}{\partial m^{\mu}} [\lambda(m) p(m,t)] \]

Next, using the Fourier integral representation of the $\delta$ function, this is the key step for introducing the response field:

\[ \delta(m-n) = \int \frac{d\tilde{n}}{2\pi} e^{-i\tilde{n}(n-m)} \]

Substituting this representation, the evolution of the master equation is written as a double integral over the physical field $m$ and the newly introduced response field $\tilde{n}$:

\[ \partial_t p(n,t) = \int dm \int \frac{d\tilde{n}}{2\pi} e^{-i\tilde{n}(n-m)} \sum_{\mu=1}^{\infty} \frac{(-\Delta n)^{\mu}}{\mu!} \frac{\partial^{\mu}}{\partial m^{\mu}} [\lambda(m) p(m,t)] \]

The advantage of this form is that the differential operators $\partial_m^\mu$ originally acting on the probability density term $\lambda(m)p(m,t)$ are now inside an integral. Therefore, these derivatives can be transferred away from the complex probability density term through integration by parts. Performing $\mu$ times integration by parts on $m$, each time the $\partial_m$ acts on the exponential term $e^{i\tilde{n}m}$ and produces a factor $i\tilde{n}$. After repeating $\mu$ times, we obtain:

\[ \partial_t p(n,t) = \int dm\ p(m,t)\lambda(m) \int \frac{d\tilde{n}}{2\pi} e^{-i\tilde{n}(n-m)} \sum_{\mu=1}^{\infty} \frac{(\Delta n)^{\mu}}{\mu!} (i\tilde{n})^{\mu} \]

Through this operation, the originally complex differential operators are transformed into simple algebraic products of the response field $\tilde{n}$. The role played by the response field $\tilde{n}$ here is completely consistent with the response field in the J-D formalism of Lecture 32—it is a field conjugate to the main state variable $n$.

2.2 Resummation and Path Integral Construction¶

Observing the summation over $\mu$ in the above equation, we can see it is a series with a clear mathematical structure:

\[ \sum_{\mu=1}^{\infty} \frac{(i\tilde{n}\Delta n)^{\mu}}{\mu!} = \left( \sum_{\mu=0}^{\infty} \frac{(i\tilde{n}\Delta n)^{\mu}}{\mu!} \right) - 1 = e^{i\tilde{n}\Delta n} - 1 \]

This series is precisely the Taylor expansion of the exponential function $e^{i\tilde{n}\Delta n}$ minus its first term (the $\mu = 0$ term). Therefore, the entire infinite-order Kramers-Moyal expansion, after introducing the response field, can be precisely resummed into a compact exponential form. Substituting this result back, the integral form of the master equation becomes:

\[ \partial_t p(n,t) = \int dm\ p(m,t)\lambda(m) \int \frac{d\tilde{n}}{2\pi} e^{-i\tilde{n}(n-m)} \left( e^{i\tilde{n}\Delta n} - 1 \right) \]

This equation is a precise integral rewriting of the original master equation, and it constitutes a solid starting point for constructing path integrals.

To construct the path integral, we follow the standard "time slicing" method. Within a small time step $\Delta t$, the evolution of the probability distribution can be approximated as:

\[ p(n,\tau + \Delta t) \approx p(n,\tau) + \Delta t\, \partial_\tau p(n,\tau) \]

Substituting the integral form of $\partial_\tau p$ derived above and approximating $1 + x\Delta t$ as $e^{x\Delta t}$, the single-step evolution operator can be written in exponential form. This expression describes the "infinitesimal propagator" from time $\tau$ to $\tau + \Delta t$. By decomposing the entire process from initial time $t_0$ to final time $t$ into a series of infinitesimal time steps and integrating over all intermediate states between each step, we can obtain the path integral representation of the conditional probability $p(n,t|n_0,t_0)$:

\[ p(n,t|n_0,t_0) = \int_{n(t_0)=n_0}^{n(t)=n} \mathcal{D}[n,\tilde{n}] e^{iS[n,\tilde{n}]} \]

Here, the integration measure $\mathcal{D}[n,\tilde{n}]$ represents summation over all possible historical paths $n(\tau)$ and all possible response field histories $\tilde{n}(\tau)$.

2.3 Action of Jump Processes¶

The above path integral is controlled by an action $S$, whose corresponding Lagrangian $\mathcal{L}$ is:

\[ S[n,\tilde{n}] = \int_{t_0}^{t} d\tau\ \mathcal{L}[n,\tilde{n}] = \int_{t_0}^{t} d\tau\ \left[ \tilde{n} \partial_\tau n - i\lambda(n) \left(e^{i\tilde{n}\Delta n}-1\right) \right] \]

To maintain consistency with the real action commonly used in statistical physics (corresponding to probability weight $e^{-S'}$), we can obtain through a Wick rotation (letting $\tilde{n} \to i\tilde{n}'$):

\[ S'[n,i\tilde{n}] = \int_{t_0}^{t} d\tau\ \left[ i\tilde{n} \partial_\tau n - \lambda(n) \left(e^{i\tilde{n}\Delta n}-1\right) \right] \]

This action summarizes all the dynamic information of the jump process:

Dynamic term $i\tilde{n} \partial_\tau n$: This is a "kinetic energy" or "Berry phase" term commonly present in MSRJD-type actions. It mathematically establishes the response field $\tilde{n}$ as the conjugate momentum of variable $n$, ensuring the path integral correctly describes continuous time evolution.

Hamiltonian term $-\lambda(n)(e^{i\tilde{n}\Delta n} - 1)$: This term can be regarded as the system's "Hamiltonian" $H(n, i\tilde{n})$, which generates the system's random jumps. It contains all specific information about the physical process:

$\lambda(n)$: The rate at which jumps occur.

$e^{i\tilde{n}\Delta n} - 1$: The operator describing the jump event itself, where $\Delta n$ is the step size of the jump.

At this point, a path integral representation starting from discrete master equations, via Kramers-Moyal expansion, has been completely constructed. This framework lays the foundation for subsequent analysis of system behavior under the "low noise limit" and establishing connections with Langevin equations.

3. Low Noise Limit and Fokker-Planck Action¶

The path integral action derived in the previous section, although mathematically precise, contains exponential terms $\exp(i\tilde{n}\Delta n)$ that make direct analytical calculations very difficult. To extract useful physical information, a reasonable approximation needs to be introduced. The physical basis of this approximation is precisely the "low noise limit" mentioned at the end of Section 1.

The "low noise limit" physically corresponds to situations where the system size is very large (i.e., average particle number $n \gg 1$). In this case, the step size $\Delta n$ of a single jump (usually a small integer like 1 or 2) is a tiny perturbation relative to the total particle number $n$. This physical picture provides a solid basis for systematically approximating the action mathematically. The goal of this lecture is to execute this approximation and reveal how discrete jump processes naturally transition to continuous diffusion processes described by Fokker-Planck equations or Langevin equations at macroscopic scales.

3.1 Expansion of the Action¶

Under the low noise limit, the Hamiltonian part $H(n, i\tilde{n}) = \lambda(n)(e^{i\tilde{n}\Delta n} - 1)$ in the action can be Taylor expanded according to the power of the small parameter $\Delta n$. Specifically, expand the exponential term $e^{i\tilde{n}\Delta n}$ to second order:

\[ e^{i\tilde{n}\Delta n} \approx 1 + (i\tilde{n}\Delta n) + \frac{1}{2}(i\tilde{n}\Delta n)^2 + \mathcal{O}((\Delta n)^3) \]

Substituting this approximation into the Hamiltonian:

\[ \lambda(n)\left(e^{i\tilde{n}\Delta n} - 1\right) \approx \lambda(n)\left[ (i\tilde{n}\Delta n) + \frac{1}{2}(i\tilde{n})^2 (\Delta n)^2 \right] \]

Now, substituting this expanded Hamiltonian back into the complete Lagrangian $\mathcal{L}' = i\tilde{n} \partial_\tau n - H(n, i\tilde{n})$:

\[ \mathcal{L}' \approx i\tilde{n} \partial_\tau n - \lambda(n)\left[ i\tilde{n}\Delta n + \frac{1}{2}(i\tilde{n})^2 (\Delta n)^2 \right] \]

Rearranging the terms, separating those linearly related to $i\tilde{n}$ from those quadratically related:

\[ \mathcal{L}' \approx i\tilde{n} \left( \partial_\tau n - \lambda(n) \Delta n \right) - \frac{1}{2} \lambda(n) (\Delta n)^2 (i\tilde{n})^2 \]

Performing time integration on this Lagrangian yields the approximated Fokker-Planck action:

\[ S_{\text{FP}}[n,i\tilde{n}] \approx \int_{t_0}^{t} d\tau\ \left[ i\tilde{n} \left( \partial_\tau n - \lambda(n) \Delta n \right) - \frac{1}{2} (i\tilde{n})^2 \lambda(n) (\Delta n)^2 \right] \]

Physical Significance: The significance of this approximation operation is revolutionary. The original action contained infinite-order $i\tilde{n}$ terms (hidden in the exponential function), corresponding to infinite-order derivatives in the Kramers-Moyal expansion. After truncating to second order, the resulting action is quadratic (Gaussian) in the response field $i\tilde{n}$. This form of action is extremely important in stochastic field theory—this is the MSRJD action introduced multiple times before, and it has a direct and profound correspondence with a continuous Langevin stochastic differential equation, as discussed in detail in Lecture 32.

3.2 Connection to Langevin Equation¶

A general Langevin equation describing a continuous random variable $n(\tau)$ can be written as:

\[ \partial_\tau n = V(n) + \sqrt{2D(n)}\, \eta(\tau) \]

where:

$V(n)$ is the deterministic drift term.
$D(n)$ is the state-dependent diffusion coefficient.
$\eta(\tau)$ is standard Gaussian white noise, satisfying $\langle \eta(\tau) \rangle = 0$ and $\langle \eta(\tau)\eta(\tau') \rangle = \delta(\tau - \tau')$.

As established in Lecture 32, the MSRJD action corresponding to such a Langevin equation has the following general form:

\[ S_{\text{Langevin}}[n,i\tilde{n}] = \int d\tau\ \left[ i\tilde{n} \left( \partial_\tau n - V(n) \right) - D(n) (i\tilde{n})^2 \right] \]

Now, we only need to compare term by term the Fokker-Planck action $S_{\text{FP}}$ derived in the previous subsection with this general form:

\[ \begin{aligned} S_{\text{FP}} &= \int d\tau\ \left[ i\tilde{n} \left( \partial_\tau n - \underline{\lambda(n) \Delta n} \right) - \underline{\frac{1}{2} \lambda(n) (\Delta n)^2} (i\tilde{n})^2 \right] \\ S_{\text{Langevin}} &= \int d\tau\ \left[ i\tilde{n} \left( \partial_\tau n - \underline{V(n)} \right) - \underline{D(n)} (i\tilde{n})^2 \right] \end{aligned} \]

By comparing the underlined parts, we can immediately identify the drift term $V(n)$ and diffusion term $D(n)$ in the equivalent Langevin equation:

Drift: $V(n) = \lambda(n)\Delta n$

Diffusion: $D(n) = \frac{1}{2}\lambda(n)(\Delta n)^2$

This result is extremely important, as it provides a systematic derivation method from microscopic jump processes to macroscopic continuous stochastic differential equations. It clearly reveals:

The macroscopic deterministic drift $V(n)$ originates from the average effect of microscopic jumps (rate $\lambda(n)$ multiplied by step size $\Delta n$).

The macroscopic random diffusion $D(n)$ originates from the second moment or fluctuations of microscopic jumps (proportional to the square of step size $(\Delta n)^2$).

This methodology constitutes a solid bridge connecting the discrete master equation world and the continuous Langevin equation world, enabling the use of powerful continuous field-theoretic tools to analyze complex systems that are essentially driven by discrete events.

4. Revisiting Annihilation Processes and Stochastic Logistic Growth¶

In Lecture 36, the bimolecular annihilation reaction $A+A \to \emptyset$ was used as an example to demonstrate the coherent-state path integral (CSP) and the concept of "imaginary noise". This lecture will revisit this model, but this time the purpose is to show how the Kramers-Moyal path integral (KMPI) derives a standard Langevin equation through a more direct, physically intuitive path under the low noise limit. This comparison will profoundly reveal the similarities and differences between the two theoretical frameworks. Subsequently, to demonstrate the wide applicability of this method, the code practice section will be replaced with a completely new model—stochastic logistic growth.

4.1 Comparison of Two Path Integrals: Annihilation Process $A+A \to \emptyset$¶

Briefly reviewing, for the reaction $A+A \to \emptyset$, we have derived its field-theoretic description through two different paths:

Coherent-State Path Integral (CSP):

Through operator algebraic methods, we obtained an exact action corresponding to an equivalent Langevin equation containing imaginary noise. The noise term in the action is $+(λ/2)x²(iq)²$, where this positive sign is the source of imaginary noise, enabling the continuous field theory to precisely reproduce the Poisson statistical properties of discrete jumps.

Kramers-Moyal Path Integral (KMPI):

Through K-M expansion of the master equation and taking the low noise limit, we obtained an approximate Langevin equation $\partial_t n = -λn(n-1) + \sqrt{4λn(n-1)}η(t)$. The corresponding noise term in the MSRJD action is $-D(n)(i\tilde{n})^2 = -2λn(n-1)(i\tilde{n})^2$. This negative sign corresponds to a standard real noise.

Core Differences and Connections:

CSP provides an exact but abstract mapping, at the cost of introducing non-physical "imaginary noise". It is more elegant in form and serves as the starting point for rigorous field-theoretic analysis (such as renormalization group).

The low noise limit of KMPI provides an approximate but intuitive mapping that directly leads to Langevin equations familiar to physicists, describing real fluctuations, and clearly connects drift and diffusion with microscopic processes.

These two methods converge from different perspectives, describing the behavior of the same physical system from different angles, reflecting the richness and inherent consistency of theoretical physics frameworks.

4.2 Application Example: Stochastic Logistic Growth¶

Now, we will apply the KMPI recipe to a model of great importance in ecology and population dynamics: logistic growth. This model contains two basic reactions:

Birth: $A \xrightarrow{\mu} 2A$ (one organism splits into two at rate $\mu$)
Competition/Death: $A+A \xrightarrow{\lambda/2} A$ (two organisms compete for resources, causing one to die, at rate $\lambda/2$)

This model describes the growth process of populations under resource-limited conditions.

1. Define Microscopic Events:

Birth: $\Delta n = +1$, rate $\lambda_1(n) = \mu n$.

Competition: $\Delta n = -1$, rate $\lambda_2(n) = (\lambda/2)n(n-1)$.

2. Calculate Total Drift and Diffusion:

Since the two processes occur independently, the total drift $V(n)$ and diffusion $D(n)$ are the sums of contributions from both processes:

Total Drift:

\[ V(n) = V_1(n) + V_2(n) = (\mu n)(+1) + \left(\frac{\lambda}{2}n(n-1)\right)(-1)= \mu n - \frac{\lambda}{2}n(n-1) \]

Total Diffusion: $$ D(n) = D_1(n) + D_2(n) = \frac{1}{2}(\mu n)(+1)^2 + \frac{1}{2}\left(\frac{\lambda}{2}n(n-1)\right)(-1)^2 = \frac{1}{2}\mu n + \frac{1}{4}\lambda n(n-1) $$

3. Construct Langevin Equation: Substituting $\partial_t n = V(n) + \sqrt{2D(n)}η(t)$, we obtain the Langevin equation for stochastic logistic growth: $$ \partial_t n = \left( \mu n - \frac{\lambda}{2}n(n-1) \right) + \sqrt{\mu n + \frac{\lambda}{2}n(n-1)}\, \eta(t) $$

Under the approximation $n \gg 1$, the equation simplifies to: $$ \partial_t n = \mu n \left(1 - \frac{n}{K}\right) + \sqrt{\left(\mu + \frac{\lambda}{2}n\right)n}\, \eta(t) $$ where $K = 2\mu/\lambda$ is the famous carrying capacity.

4.3 Python Simulation: Population Fluctuations and Carrying Capacity¶

The following code will simulate this stochastic logistic growth process and compare it with the deterministic logistic equation solution.

import numpy as np
import matplotlib.pyplot as plt

# --- 1. Parameter Settings ---
n0 = 10.0         # Initial population size
mu = 0.5          # Birth rate
lambda_rate = 0.01 # Competition rate
T_max = 100.0     # Total simulation time
dt = 0.01         # Time step
n_steps = int(T_max / dt)
n_traj = 5        # Number of stochastic trajectories

# Carrying capacity K
K = 2 * mu / lambda_rate

# --- 2. Time Axis ---
t = np.linspace(0, T_max, n_steps + 1)

# --- 3. Simulate Stochastic Trajectories ---
plt.figure(figsize=(12, 7))
# Set black background
plt.style.use('dark_background')

for j in range(n_traj):
    n_stochastic = np.zeros(n_steps + 1)
    n_stochastic[0] = n0
    for i in range(n_steps):
        n = n_stochastic[i]
        if n <= 0:
            n_stochastic[i+1:] = 0
            break

        # Drift term V(n)
        drift = mu * n - (lambda_rate / 2) * n * (n - 1)

        # Diffusion term D(n)
        diffusion_term_squared = mu * n + (lambda_rate / 2) * n * (n - 1)

        # Noise amplitude sqrt(2D(n))
        noise_amp = np.sqrt(diffusion_term_squared) if diffusion_term_squared > 0 else 0

        # Euler-Maruyama step
        noise = np.random.normal(0, 1)
        n_stochastic[i+1] = n + drift * dt + noise_amp * np.sqrt(dt) * noise

    plt.plot(t, n_stochastic, lw=1.5, alpha=0.8, label=f'Stochastic Trajectory {j+1}')

# --- 4. Calculate Deterministic (Mean Field) Solution ---
# Solution to dn/dt = μn(1 - n/K)
n_deterministic = K / (1 + (K/n0 - 1) * np.exp(-mu * t))
plt.plot(t, n_deterministic, 'r-', lw=2.5, label='Deterministic Logistic Eq.')

# --- 5. Plotting ---
plt.axhline(K, color='cyan', linestyle=':', lw=2, label=f'Carrying Capacity K = {K:.0f}')
plt.xlabel('Time', fontsize=14)
plt.ylabel('Population Size (n)', fontsize=14)
plt.title('Stochastic vs. Deterministic Logistic Growth', fontsize=16)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.ylim(bottom=0)
plt.xlim(left=0)
plt.savefig('LogisticGrowth.png')
plt.show()

Deterministic Logistic Curve (red solid line): This classic S-shaped curve is the solution to the mean field equation $\frac{dn}{dt} = \mu n\left(1 - \frac{n}{K}\right)$. It describes an ideal process where a population grows smoothly from its initial size and eventually saturates at the carrying capacity $K$.

Carrying Capacity (cyan dashed line): $K$ represents the maximum population size that the environment can stably support under given resources.

Stochastic Trajectories: These curves show the real dynamics of finite populations. The population size does not stabilize at the $K$ value but fluctuates incessantly around it. These fluctuations arising from the discreteness of individual birth and death events are called "demographic noise". In some cases, these fluctuations may even lead to unexpected population extinction, even when the $K$ value is high.

This example demonstrates the universality of the KMPI to Langevin equation derivation "recipe", which applies not only to annihilation processes in physical chemistry but also to population dynamics in ecology, as well as the universality and importance of fluctuations in nature.

5. Application II: Revisiting Directed Percolation—From Lattice to Field Perspective¶

In Lecture 34, the directed percolation (DP) model was introduced as a paradigm for nonequilibrium phase transitions, with its phase transition behavior initially explained through mean field theory. This lecture will return to the DP model, but with a completely different perspective. The goal here is no longer to introduce what DP is, but to use it as a "touchstone" to test the consistency of the powerful stochastic field-theoretic tools developed in Lectures 36 and 37—coherent-state path integrals (CSP) and Kramers-Moyal path integrals (KMPI).

This section starts from the microscopic reaction rules of DP and derives a macroscopic stochastic dynamic equation containing fluctuation effects through the systematic "recipe" of KMPI.

5.1 From Microscopic Rules to 0D Langevin Equation¶

First, let's briefly review the derivation process from DP's microscopic reaction rules to its Langevin equation in a "0D" system (i.e., well-mixed, no spatial dimensions). According to the course blackboard, the DP model consists of three reactions: proliferation, death, and coagulation. Through the KMPI framework and low noise expansion, we can obtain its drift term $V(n)$ and diffusion term $D(n)$:

Drift term: $V(n) = (\mu-\lambda)n - \gamma n(n-1)$
Diffusion term: $D(n) = \frac{1}{2} \left[ (\mu+\lambda)n + \gamma n(n-1) \right]$

This corresponds to the 0D Langevin equation describing the time evolution of total particle number $n(t)$: $$ \frac{d n}{d t} = V(n) + \sqrt{2D(n)}\, \eta(t) $$ This equation describes particle number fluctuations at an isolated point but fails to reflect the spatial propagation characteristics implied by "percolation".

5.2 From 0D to (2+1)D: Stochastic Partial Differential Equation (SPDE)¶

To describe a real physical system evolving in two-dimensional space, we need to generalize the above equation to a stochastic partial differential equation (SPDE), where the particle number $n$ becomes a continuous density field $n(x, y, t)$. For this purpose, we need to add two key spatial terms:

Spatial Diffusion: Particles spontaneously diffuse from high-density regions to low-density regions. This process is described by a Laplacian term $D_{\text{space}} \nabla^2 n$, where $D_{\text{space}}$ is the spatial diffusion coefficient.
Directed Propagation (Advection): To reflect the essence of "directed" percolation—i.e., the irreversible propagation of activity in a preferred direction (usually regarded as "time")—we introduce an advection term $-v_y \frac{\partial n}{\partial y}$. This is equivalent to applying a constant "wind" that forces the entire density field to flow at velocity $v_y$ along the negative $y$-axis direction (i.e., from top to bottom).

Combining these terms, we obtain the complete field equation describing (2+1)D (2 spatial dimensions + 1 time dimension) directed percolation:

\[ \frac{\partial n}{\partial t} = \underbrace{V(n)}_{\text{reaction}} + \underbrace{D_{\text{space}} \nabla^2 n}_{\text{lateral diffusion}} \underbrace{- v_y \frac{\partial n}{\partial y}}_{\text{longitudinal advection}} + \underbrace{\sqrt{2D(n)}\, \eta(\mathbf{x},t)}_{\text{spacetime noise}} \]

This SPDE is the final embodiment of this lecture's theoretical derivation in spacetime, providing a complete field-theoretic model that includes both internal reaction dynamics and spatial propagation characteristics.

5.3 Code Practice: Field-Theoretic Visualization of Directed Percolation¶

The following code is a direct numerical simulation of the above SPDE. It no longer tracks discrete particles but solves the evolution of a continuous density field $n(x,y,t)$. The initial condition is set as a thin "active" line at the top of space to observe how it "percolates" downward.

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import laplace
from matplotlib.animation import FuncAnimation, PillowWriter

# --- 1. Parameter Settings ---
grid_size = 128
# Reaction parameters
lambda_death = 1.0
gamma_coag = 0.1
mu_active = 1.15   # Active phase (slightly increased for clearer pattern)
mu_absorbing = 0.95 # Absorbing phase
# Simulation parameters
D_space = 0.2      # Spatial diffusion coefficient (horizontal)
v_y = 15.0         # !!! New: Vertical advection velocity !!!
dt = 0.01
total_time = 10
n_steps = int(total_time / dt)

def run_directed_spde_simulation(mu_rate):
    """ Numerically simulate the "directed" SPDE """
    # Initialization: Activate a thin line at the top
    n_field = np.zeros((grid_size, grid_size))
    start_row = 5
    n_field[start_row, grid_size//2 - 10 : grid_size//2 + 10] = 1.0

    history = [n_field.copy()]
    # Spatial step size
    dy = 1.0 

    for i in range(n_steps):
        # Reaction term
        reaction_drift = (mu_rate - lambda_death) * n_field - gamma_coag * n_field * (n_field - 1)

        # Noise term
        noise_strength_sq = (mu_rate + lambda_death) * n_field + gamma_coag * n_field * (n_field - 1)
        noise_strength = np.sqrt(np.maximum(0, noise_strength_sq))
        space_time_noise = np.random.normal(0, 1, (grid_size, grid_size)) * noise_strength / np.sqrt(dt)

        # Spatial diffusion term (Laplacian operator)
        laplacian = laplace(n_field, mode='wrap') # Periodic boundary conditions

        # !!! New: Advection term (using simple upwind scheme) !!!
        # np.roll(n_field, 1, axis=0) shifts entire rows down by one, simulating flow from above
        advection = -v_y * (n_field - np.roll(n_field, 1, axis=0)) / dy

        # Euler method to update field
        n_field += (reaction_drift + D_space * laplacian + advection + space_time_noise) * dt
        n_field = np.maximum(0, n_field) # Density cannot be negative

        # Save a frame every few steps
        if (i+1) % 2 == 0: 
            history.append(n_field.copy())

    return history

# --- 2. Run simulations for two scenarios ---
print("Running simulation for active phase...")
history_active = run_directed_spde_simulation(mu_active)
print("Running simulation for absorbing phase...")
history_absorbing = run_directed_spde_simulation(mu_absorbing)

# --- 3. Create animation ---
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))
fig.patch.set_facecolor('black')
vmax = np.max(history_active) * 0.7 # Unified color range

def update(frame):
    ax1.clear()
    im1 = ax1.imshow(history_active[frame], cmap='hot', vmin=0, vmax=vmax)
    ax1.set_title(f'Active Phase (μ={mu_active})', color='white')
    ax1.set_xticks([])
    ax1.set_yticks([])

    ax2.clear()
    im2 = ax2.imshow(history_absorbing[frame], cmap='hot', vmin=0, vmax=vmax)
    ax2.set_title(f'Absorbing Phase (μ={mu_absorbing})', color='white')
    ax2.set_xticks([])
    ax2.set_yticks([])

    fig.suptitle(f'Field Simulation of Directed Percolation (Time: {frame*dt*2:.2f})', color='white', fontsize=16)

    # Force aspect ratio to stretch image, better reflecting "directed" nature
    ax1.set_aspect(0.5)
    ax2.set_aspect(0.5)

    return [im1, im2]

# --- 4. Generate and save animation ---
frames_to_render = min(len(history_active), len(history_absorbing))
ani = FuncAnimation(fig, update, frames=frames_to_render, interval=40, blit=False)
ani.save("directed_percolation_field.gif", writer=PillowWriter(fps=25))
plt.show()

Directed Propagation: The mandatory downward advection term $(-v_y \cdot \partial n / \partial y)$ dominates the overall motion direction of the system, forming the classic "percolation cone" propagating from top to bottom.

Lateral Expansion: The spatial diffusion term $(D_{\text{space}} \cdot \nabla^2 n)$ causes the percolation cone to expand laterally (spatial dimensions) while propagating, forming the cone's boundaries.

Left Figure (Active Phase, $\mu > \lambda$): Within the percolation cone, the proliferation term dominates, allowing the density field to maintain and develop complex, dynamically fluctuating internal structures. Activity successfully propagates to the system's "future".

Right Figure (Absorbing Phase, $\mu < \lambda$): Although an initial percolation cone is also formed, the death rate within is too high, causing the density field to decay rapidly. Although the noise term can cause temporary local brightening, it cannot change the overall fate of extinction, and eventually the activity will completely extinguish within a limited "time".

5.4 Theoretical Unity: Cole-Hopf Transformation¶

After constructing path integrals from two completely different perspectives—probability theory (KMPI) and operator algebra (CSP)—a crucial question arises: What connections exist between them? The answer lies in a profound nonlinear field transformation—the Cole-Hopf transformation, which in the stochastic processes literature is also often called the Grassberger transformation.

The Core of the Transformation¶

The core role of the Cole-Hopf transformation in the current context is to establish precise mathematical mappings between the field variables $(n, \tilde{n})$ in the KMPI framework and the "symbols" $(x, iq+1)$ of operators in the CSP framework. According to the course blackboard, the key correspondences are:

Equivalence of Jump Operator and Annihilation Operator:

$$ e^{i\tilde{n}} \quad \longleftrightarrow \quad \hat{a} \quad (\text{whose symbol is } iq+1) $$ This relationship is the cornerstone of the entire transformation. It reveals that the operator $e^{i\tilde{n}}$ describing a single particle "jump" in KMPI is algebraically equivalent to the operator $\hat{a}$ that "annihilates" a particle in CSP.

Equivalence of Particle Number Field and Particle Number Operator: $$ n \quad \longleftrightarrow \quad \hat{a}^\dagger \hat{a} \quad (\text{whose symbol is } x(iq+1)) $$

This indicates that the basic variable in KMPI—particle number $n$—is a composite quantity in the CSP framework, namely the particle number operator $\hat{a}^\dagger \hat{a}$.

The significance of this transformation goes far beyond mathematical symbol substitution; it reveals the deep unity of theory:

Unified Physical Picture: It proves that the "probability theory" perspective (KMPI focusing on changes in particle number $n$) and the "algebraic" perspective (CSP focusing on basic events of creation $\hat{a}^\dagger$ and annihilation $\hat{a}$) are completely equivalent when describing the same physical reality.
Transforming Nonlinear to Linear: Historically, the most famous application of the Cole-Hopf transformation is to precisely transform the nonlinear Burgers' equation describing shock waves in fluid mechanics into the linear heat conduction equation. This ability to simplify complexity and reveal hidden linear structures behind complex nonlinear problems is precisely why it is so powerful in theoretical physics.

In this lecture, the Cole-Hopf transformation plays the role of a "translator". Despite different derivation paths, KMPI and CSP must converge to the same Langevin equation describing macroscopic physics under the low noise limit. This provides theoretical justification for choosing the most appropriate field-theoretic tool according to the specific situation of the problem (KMPI is more intuitive, CSP is more convenient for using advanced tools like renormalization group).

KMPI provides a path from probability theory that is physically intuitive and leads directly to Langevin equations; while CSP provides a path from operator algebra that is more elegant in form and convenient for using more advanced field-theoretic tools. The final convergence of the two, and the SPDE they jointly derive that can be intuitively visualized.

Conclusion¶

This lecture systematically introduced how to construct an equivalent continuous stochastic field-theoretic description starting from master equations describing discrete jump processes, and ultimately revealed the inherent unity of stochastic process theoretical frameworks.

Bridge from Discrete to Continuous: The Kramers-Moyal expansion was established as a crucial mathematical bridge. It provides a systematic method to transform discrete, integro-differential master equations into infinite-order partial differential equations that treat state variables as continuous.

Precise Path Integral Formulation: Based on the K-M expansion, a Kramers-Moyal path integral (KMPI) that precisely represents the original master equation was constructed. Its action $S[n,\tilde{n}]$ elegantly encodes the jump rate $\lambda(n)$ and step size $\Delta n$ in a compact exponential "Hamiltonian" $\lambda(n)(e^{i\tilde{n}\Delta n} - 1)$, preserving the complete non-Gaussian statistical properties of jump processes.

Physical Significance of the Low Noise Limit: By expanding the action under the physically reasonable low noise limit (equivalent to system size expansion), a Gaussian-form Fokker-Planck (or MSRJD) action is naturally obtained. This process not only provides theoretical justification for approximating Fokker-Planck equations from master equations, but more importantly, establishes a systematic derivation "recipe" from microscopic jump rules (rates and step sizes) to macroscopic Langevin dynamics (drift and diffusion coefficients).

Theoretical Unity: The climax of this lecture lies in revealing the profound unity of theory through the analysis of the directed percolation (DP) model. The results show that under the low noise limit, the macroscopic Langevin equation derived by the probabilistic method (KMPI) is completely consistent with the results derived by the completely different operator algebraic method (coherent-state path integral, CSP) in Lecture 36. The Cole-Hopf transformation reveals the deep duality between these two seemingly unrelated formalisms, proving that regardless of which fundamental perspective we start from, the universal macroscopic physics described is the same.

The single-species theoretical framework $(n, \tilde{n})$ developed in this lecture lays a solid foundation for studying more complex systems. The next natural step is to extend these techniques to systems containing multiple interacting species. At that point, there will no longer be just one particle number density $n$ and its response field $\tilde{n}$, but a set of fields $\{n_i, \tilde{n}_i\}$ will be introduced for each species $i$.

The structure of the action (whether KMPI or CSP form) will become richer, with coupling terms appearing in addition to terms describing each species' own dynamics, used to describe interactions between species, such as:

Predator-Prey: $A + B \to 2A$
Competition: $A + B \to A$ or $A + B \to B$
Cyclic Dominance: $A + B \to 2B$, $B + C \to 2C$, $C + A \to 2A$

This generalization will enable the study of much richer dynamic phenomena than simple decay/growth processes, such as the sustained oscillations of population numbers that appear in predator-prey models (Lotka-Volterra), or the spontaneously emerging chase-escape spiral patterns in cyclic competition models (such as "rock-paper-scissors" games). This will lead us from single-variable stochastic processes into the fascinating world of complex population dynamics and spatiotemporal pattern formation driven by multi-field coupling.