Introduction: Toward the Nonequilibrium World¶
In previous lectures, we established a complex field-theoretic framework for describing the dynamics of near-equilibrium systems. Starting from the Onsager-Machlup functional, we introduced the more powerful Janssen-De Dominicis (J-D) action and, based on it, systematically derived the fluctuation-dissipation theorem (FDT), a cornerstone of near-equilibrium statistical physics. Subsequently, Lecture 32 broadened our perspective to far-from-equilibrium, showing through the Jarzynski equality and the Crooks fluctuation theorem that, during arbitrary finite-time protocols, nonequilibrium work and equilibrium free-energy differences obey deep equality relations. While these theories are powerful, their focus is mainly on how systems respond to external perturbations or are driven between two equilibrium states.
This lecture marks another crucial shift in the course perspective: from transient nonequilibrium processes to systems that are continuously driven and reach nonequilibrium steady states (NESS). The central question here is no longer relaxation alone, but whether such driven-dissipative systems can exhibit phase-transition-like behavior akin to equilibrium systems. These nonequilibrium phase transitions are not governed by free-energy minimization; rather, they arise from a complex dynamical balance between drive and dissipation and often display new critical phenomena.
The real world abounds with such systems. A living cell is the most intuitive example of a system maintained in a nonequilibrium steady state (NESS): it maintains a highly ordered, far-from-equilibrium state through continuous exchange of energy and matter. The core feature of many nonequilibrium phase transitions is precisely the transition from a fluctuating, active phase to an irreversible, quiescent absorbing state.
A typical example is epidemic spreading: when an infectious disease spreads in a population, the system is in an active phase with constantly changing numbers of infected individuals; once the last infected individual recovers or is isolated, the system enters the absorbing zero-infection state, from which infection cannot spontaneously reappear. This transition from outbreak to extinction depends on whether the basic reproduction number \(R_0\) exceeds the threshold 1. A similar paradigm can be used to describe forest fire spreading, where the active phase is burning flames and the absorbing state is complete extinction of flames; this transition is controlled by whether tree density exceeds a percolation threshold. This transition from active to permanent "silence" can even find its counterpart in the microscopic chemical world, for example on catalytic surfaces, where persistent chemical reactions form an active phase, whereas irreversible coverage of active sites by a poison halts reactions permanently—an absorbing, poisoned state. These examples clearly demonstrate that, from macroscopic ecosystems to microscopic reactions, transitions to absorbing states are extremely widespread and important nonequilibrium phenomena in nature.
In the study of such nonequilibrium phase transitions, a core concept is the absorbing state. An absorbing state refers to a state (or set of states) that, once entered, cannot be left. This concept has already been encountered in previous linear birth-death processes, where the state of zero population is a typical absorbing state—once extinct, the population remains zero without external intervention. The presence of absorbing states fundamentally breaks detailed balance, injecting irreversibility. The transition from a fluctuating active phase to a frozen absorbing phase defines a distinct and important class of nonequilibrium phase transitions.
To study these phenomena concretely, we need a minimal model. In the field of nonequilibrium phase transitions, Directed Percolation (DP) plays the role that the Ising model does in equilibrium phase transition theory—a foundational cornerstone. In its most concise form, it captures the continuous transition from an active, continuously fluctuating phase to a quiescent, static absorbing phase. The goal of this lecture is to deeply analyze the DP model, first through a simplified mean-field approach to understand its basic physical picture. However, as emphasized in previous lectures, mean-field theory ignores crucial fluctuations near criticality. Therefore, to build a more complete theory, it is necessary to revisit and develop the mathematical tools for handling discrete stochastic processes—the Master Equation and its solution methods—as a solid foundation for subsequent field-theoretic models that can accurately describe fluctuations.
1. Mean-Field Theory of Directed Percolation¶
Directed Percolation (DP) is a paradigmatic model for how fluctuation processes with intrinsic directionality propagate through a medium. Unlike ordinary percolation (such as water permeating porous media), the "directed" constraint emphasizes the irreversibility of the process, typically corresponding to the unidirectional flow of time. The model precisely captures the essence of a class of nonequilibrium phase transitions, namely the critical behavior of a system transitioning from a self-sustaining, continuously active active phase to an absorbing phase where dynamics eventually come to a complete halt. In the real world, such phenomena are common: for example, epidemic spreading, where the active phase is the spread of disease in a population, and the absorbing phase is the "disease-free" state with zero infections; forest fire spreading follows the same pattern, with burning flames as the active phase and complete flame extinction as the absorbing phase. Therefore, understanding the DP model is the theoretical starting point for studying such nonequilibrium critical phenomena.
To capture the essence of absorbing-state phase transitions mathematically, we need to construct a kinetic model that captures particle creation, annihilation, and interactions. The goal of this section is to first establish the microscopic rules of the Directed Percolation (DP) model, and then, through the mean-field approximation as a simplification method, ignore random fluctuations and spatial details to derive a deterministic macroscopic kinetic equation. By analyzing the steady-state solutions of this equation, we reveal how the phase transition from "active" to "quiescent" occurs in the DP model.
1.1 Microscopic Reaction Mechanisms¶
To capture the essence of an absorbing-state transition, we set up reactions encoding particle creation, annihilation, and interactions. Consider a particle system defined on a lattice, where particles (denoted as A) undergo a series of random birth-death processes. These processes constitute the microscopic kinetic rules of the model and can be summarized as three basic chemical reactions:
1. Spontaneous Decay: \(A \xrightarrow{\alpha_0} \varnothing\)
A particle A spontaneously disappears or dies at rate \(\alpha_0\). Here \(\varnothing\) represents an empty site or zero-particle state. This is the decay mechanism driving the system toward "quiescence".
2. Offspring Creation (Branching): \(A \xrightarrow{\alpha_2} 2A\)
A particle A produces a new particle at rate \(\alpha_2\), increasing the total particle count by one. This represents the growth mechanism of the system, the source maintaining system "activity".
3. Coagulation: \(2A \xrightarrow{\beta_1} A\)
Two adjacent particles A merge into one particle at rate \(\beta_1\). This process becomes important at high densities, simulating growth-inhibiting effects due to resource competition, stabilizing the particle count and preventing infinite growth.
This simple set of rules contains the three core elements: life (reproduction), death (decay), and interaction (coagulation). The competition and balance among them will collectively determine the macroscopic fate of the entire system.
1.2 Deterministic Rate Equation¶
To transition from these microscopic, random rules to a description of the system's macroscopic behavior, we first adopt the mean-field approximation. The core idea of this approximation is:
Neglect spatial correlations: Assume the system is "well-mixed", so particles have equal probability of appearing at any position.
Neglect random fluctuations: Replace discrete, random particle numbers with a continuous macroscopic average—the particle density \(n(t)\)—to describe the entire system's state.
Under this approximation, the temporal evolution of particle density follows a deterministic ordinary differential equation, the rate equation:
This equation is actually a direct application of the law of mass action in chemical reaction kinetics. Each term has a clear physical origin:
Linear term \((\alpha_2 - \alpha_0)n(t)\): represents the net effect of single-particle processes.
\(+\alpha_2 n(t)\): the rate of increase in total particle number. It is proportional to the existing particle density \(n(t)\) that can reproduce.
\(-\alpha_0 n(t)\): the rate of decrease in total particle number. It is proportional to the existing particle density \(n(t)\) that can decay.
Nonlinear term \(-\beta_1 n^2(t)\): represents the effect of two-particle processes.
Under mean-field approximation, the probability of two particles meeting is proportional to the square of particle density \(n^2(t)\). This term is nonlinear and is crucial for preventing explosive particle growth, embodying the inhibitory effect at high densities.
Through this averaging procedure, the equation transforms discrete particle birth-death events into the deterministic evolution of a continuous variable \(n(t)\). Although this simplification ignores the critical role of fluctuations, it can reveal different macroscopic phases that the system may exhibit at minimal cost.
1.3 Steady-State Analysis and the Absorbing-State Transition¶
The long-term behavior of the system is determined by its steady states, i.e., states where particle density no longer changes with time, so \(\frac{d n}{dt} = 0\). By solving the algebraic equation \((\alpha_2 - \alpha_0)n - \beta_1 n^2 = 0\), we can find all fixed points (steady-state solutions).
The key to the analysis is introducing a control parameter \(\Delta = \alpha_2 - \alpha_0\). This parameter has a clear physical meaning: it represents the net growth rate of a single particle in the absence of interactions (i.e., at very low density, \(n \to 0\)). Depending on the sign of \(\Delta\), the system exhibits two radically different macroscopic behaviors:
1. Absorbing Phase (\(\Delta < 0\))
When the decay rate \(\alpha_0\) exceeds the reproduction rate \(\alpha_2\), the net growth rate is negative. At this point, the only physical steady-state solution of the rate equation is \(n_{ss} = 0\).
Physical meaning: This means that any initial particle population will eventually vanish, and the system inevitably enters a "quiescent" state with no particles. This \(n = 0\) state is exactly an absorbing state, because once entered, reproductive processes cannot initiate, and the system can never escape.
2. Active Phase (\(\Delta > 0\))
When the reproduction rate \(\alpha_2\) exceeds the decay rate \(\alpha_0\), the net growth rate is positive. At this point, aside from \(n = 0\) (which is now unstable), a new, stable nonzero steady-state solution emerges: $$ n_{act} = \frac{\Delta}{\beta_1} = \frac{\alpha_2 - \alpha_0}{\beta_1} $$
Physical meaning: In this phase, the system can maintain a finite, dynamically balanced particle density. The net growth of particles is balanced by the coagulation process at high densities, forming a persistent, continuously fluctuating "living" state.
At \(\Delta = 0\), the system is exactly at the critical point from active phase to absorbing phase. This is a continuous phase transition, where the steady-state density of the active phase \(n_{act}\) plays the role of an order parameter in this transition, growing continuously from zero. This transition from a dynamically active state to a static absorbing state is called an absorbing-state phase transition.
2. Critical Phenomena in Mean Field¶
After determining that the system has two macroscopic phases, the next step is to deeply study the behavior near the phase transition critical point. In physics, the dynamics and scaling behavior near critical points often reveal universal laws of the system, which do not depend on the microscopic details of the model. This section will analyze the dynamical solutions of the mean-field rate equation, revealing the universal phenomenon of critical slowing down, introducing critical exponents to quantitatively characterize critical behavior, and finally proposing the important conjecture of the Directed Percolation universality class.
2.1 Dynamical Relaxation and Critical Slowing Down¶
In addition to steady-state properties, the dynamical process of how the system approaches the steady state also contains important information about phase transitions. The above nonlinear rate equation is a Bernoulli equation that can be solved exactly. Given the initial density \(n(t=0) = n_0\), its solution is:
When \(t \to \infty\), the exponential term disappears, and the fraction simplifies to \(\frac{n_0 n_{act}}{n_0} = n_{act}\), correctly approaching the steady-state density of the active phase.
From this solution, we can analyze the system's relaxation behavior:
Exponential Relaxation (\(\Delta > 0\))
Within the active phase, the system approaches the steady-state density \(n_{act}\) exponentially. From the denominator of the solution, we can see that its characteristic relaxation time is:
This relaxation time scale measures how fast the system recovers from a perturbed state to the steady state.
Critical Slowing Down¶
When the system approaches the critical point from the active phase, i.e., \(\Delta \to 0^+\), the relaxation time \(\tau_{RELAX} \to \infty\). This means that near the critical point, the system's evolution becomes extremely slow, and the "memory" of the initial state persists for a very long time.
Physical Origin: This is one of the universal characteristics of all continuous phase transitions. Physically, this arises because the "force" driving the system back to the steady state (embodied by \(\Delta\) here) disappears at the critical point, causing the system's response to perturbations to become infinitely slow.
Power-Law Relaxation at the Critical Point (\(\Delta = 0\))
Exactly at the critical point, the dynamical behavior undergoes a qualitative change. Exponential relaxation is replaced by a slower algebraic relaxation (or power-law relaxation). At this point, the rate equation simplifies to \(\frac{\partial n(t)}{\partial t} = -\beta_1 n^2(t)\), with solution:
The system density decays with time in a \(t^{-1}\) power-law form, which is much slower than any exponential decay.
2.2 Critical Exponents and Universality Class¶
To characterize critical behavior in a parameter-independent way, we introduce critical exponents:
- Order-parameter exponent \(\beta\): \(n_{\mathrm{act}} \sim |\Delta|^{\beta}\). Mean field gives \(n_{\mathrm{act}} = \Delta/\beta_1\), hence \(\beta=1\).
- Relaxation-time exponent \(\nu z\): \(\tau_{\mathrm{RELAX}} \sim |\Delta|^{-\nu z}\). Mean field gives \(\tau_{\mathrm{RELAX}} = 1/\Delta\), hence \(\nu z = 1\).
Adding diffusion yields a reaction-diffusion equation $$ \frac{\partial n(\mathbf{x}, t)}{\partial t} = (\alpha_2 - \alpha_0) n(\mathbf{x}, t) - \beta_1 n^2(\mathbf{x}, t) + D \, \nabla^2 n(\mathbf{x}, t). $$ where \(D\) is the diffusion coefficient and \(n(\mathbf{x}, t)\) is the density field depending on space and time. Diffusion introduces a correlation length \(\xi\) describing spatial range of fluctuations; near criticality, \(\xi \sim |\Delta|^{-\nu}\). Time and length scales are linked by the dynamic exponent \(z\): \(\tau \sim \xi^z\). Mean-field analysis yields the DP exponents $$ \beta = 1, \quad \nu = \tfrac{1}{2}, \quad z = 2. $$
This motivates a deep conjecture: all nonequilibrium transitions with a single absorbing state, a scalar order parameter, and short-range interactions belong to the DP universality class. That is, regardless of microscopic details (rates, lattice structure, etc.), critical behavior near the transition is governed by the same exponents as DP.
However, precise simulations and experiments show that the true exponents (e.g., in 1+1D, \(\beta \approx 0.276\)) differ markedly from mean-field predictions-revealing mean field's limitation: it neglects crucial fluctuations at criticality. To go beyond mean field, we need a field theory that treats stochastic fluctuations properly, and the starting point is the Master Equation.
3. Mathematical Tools for Discrete Stochastic Processes¶
Mean-field theory yields a macroscopic, deterministic picture of DP, but its main flaw is the neglect of fluctuations-central to critical phenomena. For a more fundamental, accurate theory, we return to direct descriptions of discrete particle birth-death processes, whose language is the Master Equation.
Directly tackling master equations with nonlinear interactions (e.g., DP coagulation) is difficult. As a "detour," we first analyze a simple, exactly solvable process-the linear death process-to develop a toolbox for solving master equations.
3.1 Framework Recap and Challenge¶
For continuous variables, we built a tiered framework (see Lecture 31):
- Microscopic: Langevin equation, e.g., \(dx = A(x)\, dt + C(x)\, dW\), describing single stochastic trajectories.
- Mesoscopic: Fokker-Planck equation for the deterministic evolution of the probability density \(P(\mathbf{x}, t)\).
- Global: path-integral representation, a functional integral over histories with an action \(S\) encoding the dynamics.
Our challenge is to construct an analogous framework for discrete variables (e.g., particle number \(n \in \{0,1,2,\dots\}\)). For such systems, the starting point is the Master Equation, not the Langevin equation.
3.2 Linear Death Process and the Master Equation¶
To develop the mathematical tools for solving master equations, we begin with the simplest nontrivial stochastic process: the linear death process. This process corresponds to spontaneous decay: \(A \xrightarrow{\lambda} \varnothing\), where each particle independently disappears at a constant per-capita rate \(\lambda\).
The Master Equation describes the temporal evolution of \(P_n(t)\), the probability of having exactly \(n\) particles at time \(t\). It is a balance equation for probability flow, with the following form:
The physical meaning of this equation embodies the "gain-loss" balance idea established in Lecture 7:
Gain Term: \(+\lambda(n+1)P_{n+1}(t)\)
The only way for the system to enter state \(n\) is by transitioning from state \(n+1\). When the system has \(n+1\) particles, the total death rate is \(\lambda(n+1)\). Therefore, the total probability flow from state \(n+1\) into state \(n\) is proportional to the product of the source state's probability \(P_{n+1}(t)\) and the transition rate.
Loss Term: \(-\lambda n P_n(t)\)
If the system is in state \(n\), it will leave this state when any of its \(n\) particles decays. The total departure rate is \(\lambda n\). Therefore, the total probability flow out of state \(n\) is proportional to the product of the current state's probability \(P_n(t)\) and this rate.
This equation is actually an infinite-dimensional, coupled system of linear ordinary differential equations, and solving it directly is quite cumbersome.
3.3 Generating-Function Method: From Infinite to Single¶
To solve this problem, we introduce a powerful mathematical tool already used in Lecture 4—the generating function \(G(x, t)\). It encodes the infinite sequence of probabilities \(P_n(t)\) as a single function of an auxiliary variable \(x\):
The advantage of this function is that it converts discrete operations on different values of \(n\) in the Master Equation into continuous differentiation on the generating function. Multiplying the Master Equation by \(x^n\) and summing over all \(n\) yields a single first-order partial differential equation for the generating function:
Left side: \(\sum_n \frac{dP_n}{dt} x^n = \frac{\partial}{\partial t} \sum_n P_n x^n = \frac{\partial G}{\partial t}\).
Loss term: \(\sum_n (-\lambda n P_n) x^n = -\lambda x \frac{\partial}{\partial x} \sum_n P_n x^n = -\lambda x \frac{\partial G}{\partial x}\).
Gain term: \(\sum_n \lambda (n+1) P_{n+1} x^n = \lambda \frac{\partial}{\partial x} \sum_{n+1} P_{n+1} x^{n+1} = \lambda \frac{\partial G}{\partial x}\).
By introducing the generating function, we have successfully transformed an infinite-dimensional system of ordinary differential equations into a single partial differential equation, paving the way for analytical solution.
4. Spectral Method and Physical Interpretation¶
The previous section successfully transformed the Master Equation of the linear death process into a partial differential equation for the generating function. Now, we introduce a deeper and more powerful solution method—the spectral method. Its core idea is to borrow operator theory from quantum mechanics, transforming the problem of solving a partial differential equation into finding the eigenvalues and eigenfunctions of an operator. This method not only gives the final solution but also reveals the profound algebraic structure underlying the system dynamics.
4.1 Operator Form and Quantum Analogy¶
We can rewrite the partial differential equation for the generating function in operator form. Define a linear operator \(\mathcal{H}\):
Then the equation can be written in a form very similar to the (imaginary-time) Schrödinger equation:
The formal solution of this equation is \(G(t) = e^{-\mathcal{H} t} G(0)\). To obtain the explicit solution, we need to solve the eigenvalue problem for the operator \(\mathcal{H}\):
where \(\lambda_j\) are eigenvalues and \(\varphi_j(x)\) are the corresponding eigenfunctions.
Physical meaning: Although \(\mathcal{H}\) formally resembles a Hamiltonian, it describes the evolution of classical probabilities, not quantum wave functions. The value of this analogy lies in our ability to borrow mature algebraic tools from quantum mechanics to solve classical stochastic processes.
4.2 Ladder Operators and the Spectrum¶
To solve this eigenvalue problem algebraically, we introduce a pair of ladder operators, adopting the same approach used for the quantum harmonic oscillator:
Annihilation operator: \(a := \frac{\partial}{\partial x}\)
Creation operator: \(a^{\dagger} := (x-1)\)
These two operators satisfy a commutation relation similar to bosonic operators:
Using these two operators, the Hamiltonian operator \(\mathcal{H}\) can be expressed concisely as:
By analyzing the commutation relations between the operator \(\mathcal{H}\) and the ladder operators, we can show that \(a^{\dagger}\) raises an eigenvalue by \(\lambda\), while \(a\) lowers it by \(\lambda\). Because probabilities must remain bounded, the spectrum must be bounded from below; hence there must exist a ground state \(\varphi_0\) annihilated by \(a\):
\(a \varphi_0(x) = \frac{\partial \varphi_0(x)}{\partial x} = 0 \implies \varphi_0(x) = \text{const}\). We can take \(\varphi_0(x) = 1\).
The corresponding ground-state eigenvalue is \(\lambda_0 = 0\).
Starting from the ground state, the entire spectrum can be constructed by repeatedly applying the creation operator \(a^{\dagger}\):
Eigenvalues: \(\lambda_j = j \lambda\), for \(j \in \{0, 1, 2, \ldots\}\)
Eigenfunctions: \(\varphi_j(x) \propto (a^{\dagger})^j \varphi_0 \propto (x - 1)^j\)
4.3 Full Time-Dependent Solution and Interpretation¶
With the complete eigenvalue spectrum, we can expand the generating function \(G(x, t)\) in this eigenfunction basis:
where the coefficients \(G_{j0}\) are determined by the initial condition. Assuming the system has exactly \(N\) particles at \(t = 0\), i.e., \(P_n(0) = \delta_{n,N}\), the corresponding initial generating function is \(G(x, 0) = x^N\).
To find the coefficients \(G_{j0}\), we cleverly use the binomial theorem to expand \(x^N\) about \((x-1)\):
By comparing coefficients, we immediately obtain \(G_{j0} = \binom{N}{j}\) for \(j \le N\), and \(G_{j0} = 0\) for \(j > N\).
Substituting the coefficients back into the time-dependent solution and again using the binomial theorem for summation, we obtain the final compact form of the generating function:
The final step is to extract the probability \(P_n(t)\) we actually care about from the generating function. According to the definition of the generating function, \(P_n(t)\) is the coefficient of \(x^n\) in the expansion of \(G(x,t)\). Expanding the above expression again using the binomial theorem:
From this, we can directly read off the probability distribution:
This is a binomial distribution.
This exact solution has extremely elegant and intuitive physical meaning. At time \(t\), each of the initial \(N\) particles has probability \(p = e^{-\lambda t}\) of surviving (this is the solution for a single particle's survival probability) and probability \(1 - p = 1 - e^{-\lambda t}\) of having decayed. Since the decay of each particle is an independent event (this is the fundamental assumption of a "linear" process), the probability of having exactly \(n\) particles survive at time \(t\) equals the probability of \(n\) "successes" (survivals) in \(N\) independent Bernoulli trials—this is exactly the definition of a binomial distribution. This exact solution perfectly validates the power of spectral methods as an abstract mathematical tool and profoundly reveals how the independence of microscopic processes manifests as binomial statistics at the macroscopic probability distribution level.
5. Code Practice¶
5.1 Code Practice I: Visualizing a Directed-Percolation Process (Lattice Model vs Mean Field)¶
To build an intuitive understanding of Directed Percolation (DP), we first use a 2D Lattice Model to directly observe the dynamic process of "percolation". This stochastic simulation shows how "active" particles propagate and die out in space. Subsequently, we will show how mean-field theory describes the macroscopic average behavior of this complex process at a higher level of abstraction by ignoring spatial details.
The code below simulates a simplified version of the Directed Percolation process (also known as the Contact Process). On a two-dimensional grid, each "active" (red) cell attempts to activate its neighbors, while itself has a certain probability of becoming "inactive" (black). We will compare the system evolution under two different spreading probabilities: one that eventually dies out (absorbing phase), and another that can sustain propagation (active phase).
# --- Use dark background for plots ---
import matplotlib.pyplot as plt
plt.style.use('dark_background')
import numpy as np
import matplotlib.animation as animation
from matplotlib.animation import PillowWriter
# --- Model parameters ---
L = 60 # lattice size (1D or 2D projected as needed)
T = 120 # time steps
p_spread_low = 0.20 # absorbing phase (below threshold)
p_spread_high = 0.35 # active phase (above threshold)
def dp_step(state, p):
# simple 1D nearest-neighbor directed spreading to next time layer
new_state = np.zeros_like(state)
for i in range(1, len(state)-1):
if state[i-1] or state[i] or state[i+1]:
if np.random.rand() < p:
new_state[i] = 1
return new_state
def simulate_dp(p):
grid = np.zeros((T, L), dtype=int)
mid = L // 2
grid[0, mid] = 1
for t in range(1, T):
grid[t] = dp_step(grid[t-1], p)
return grid
grid_low = simulate_dp(p_spread_low)
grid_high = simulate_dp(p_spread_high)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
im1 = ax1.imshow(grid_low, aspect='auto', cmap='magma', origin='lower')
ax1.set_title('Absorbing Phase (Low Spread Probability)')
im2 = ax2.imshow(grid_high, aspect='auto', cmap='magma', origin='lower')
ax2.set_title('Active Phase (High Spread Probability)')
plt.tight_layout()
ani = animation.ArtistAnimation(fig, [], interval=100)
plt.close(fig)
ani.save('dp_lattice_simulation.gif', writer=PillowWriter(fps=6))
Left panel (Absorbing Phase): Lower spreading probability. Although the initial active cluster (bright yellow) expands slightly in the early stages, the death rate exceeds the spreading rate. Eventually, all active sites vanish, and the system enters a completely black, static absorbing state.
Right panel (Active Phase): Higher spreading probability. The initial active cluster successfully spreads outward, forming a persistent, continuously changing dynamic pattern. Despite local birth and death, overall the system remains "alive"—this is the active phase.
This simulation demonstrates that the system has a critical spreading probability; crossing this threshold causes a qualitative change in the system's long-term behavior.
5.2 Code Practice II: Mean-Field Dynamics and Phase Flow¶
Now we return to mean-field theory. It ignores all spatial structure and random details in the above simulation, focusing only on the system's overall average density \(n(t)\). The following animation visualizes the evolution of the average density in phase space, showing the dynamic laws at the macroscopic level.
# --- Use dark background for plots ---
import matplotlib.pyplot as plt
plt.style.use('dark_background')
import numpy as np
import matplotlib.animation as animation
from matplotlib.animation import PillowWriter
# --- Parameter setup ---
alpha0 = 0.20
alpha2_values = [0.10, 0.20, 0.30]
beta1 = 0.10
dt = 0.02
T = 30.0
steps = int(T/dt)
def mean_field_traj(alpha2, n0):
n = np.zeros(steps)
n[0] = n0
for k in range(steps-1):
n[k+1] = n[k] + dt*((alpha2 - alpha0)*n[k] - beta1*n[k]*n[k])
n[k+1] = max(n[k+1], 0)
return n
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_title('Mean-Field DP Dynamics')
ax.set_xlabel('Time')
ax.set_ylabel('Density n')
for a2 in alpha2_values:
n = mean_field_traj(a2, n0=0.2)
ax.plot(np.linspace(0, T, steps), n, label=f'alpha2={a2:.2f}')
ax.legend()
plt.tight_layout()
plt.close(fig)
ani2 = animation.ArtistAnimation(plt.figure(), [], interval=100)
ani2.save('dp_mean_field_dynamics.gif', writer=PillowWriter(fps=6))
This phase-space flow animation is an abstract summary of the macroscopic behavior in the 2D lattice simulation. It demonstrates the system's fate as predicted by mean-field theory:
Active Phase (green): The system has a stable nonzero density fixed point, and the dynamical flow "pushes" the system toward this state.
Absorbing Phase (red): The system's only stable fixed point is zero density, and all initial states will "flow" toward extinction.
Critical Point (blue): As the system evolves near zero, the restoring force (i.e., \(\partial_t n\)) becomes extremely small, leading to critical slowing down.
5.3 Code Practice III: Dynamic Probability Distribution of the Linear Death Process¶
For stochastic processes, we care more about the evolution of the entire probability distribution. The following code will simulate the linear death process and dynamically display how its probability distribution \(P(n, t)\) evolves from an initial certain state (delta distribution) to the final absorbing state (all particles dead).
This animation will intuitively display the solution of the Master Equation and compare it in real time with the theoretical binomial distribution solution we derived.
# --- Use dark background for plots ---
import matplotlib.pyplot as plt
plt.style.use('dark_background')
import numpy as np
from matplotlib.animation import FuncAnimation, PillowWriter
from math import comb
# --- Simulation Parameters ---
N = 50
lambda_rate = 0.2
num_trials = 20000
# --- Run full Gillespie simulations to get all trajectories ---
def simulate_death(N, lambda_rate, t):
# For linear death, survival is independent with p=exp(-lambda*t)
p = np.exp(-lambda_rate * t)
# Sample survivor count from Binomial(N, p)
return np.random.binomial(N, p)
# Pre-calculate data for animation
times = np.linspace(0, 5.0, 100)
fig, ax = plt.subplots(figsize=(9, 6))
bar_container = ax.bar(range(N+1), np.zeros(N+1), alpha=0.7, label='Simulation (hist)')
line, = ax.plot(range(N+1), np.zeros(N+1), 'r-', lw=2, label='Theory (binomial)')
time_text = ax.text(0.95, 0.90, '', transform=ax.transAxes, ha='right')
def update_death_process(frame):
current_time = times[frame]
# Simulation histogram
samples = [simulate_death(N, lambda_rate, current_time) for _ in range(num_trials)]
counts_sim, _ = np.histogram(samples, bins=np.arange(-0.5, N+2) - 0.5)
freq_sim = counts_sim / num_trials
# Update bars
for count, rect in zip(freq_sim, bar_container.patches):
rect.set_height(count)
# Theory (binomial)
n_values = np.arange(0, N+1)
p_survival = np.exp(-lambda_rate * current_time)
prob_theory = comb(N, n_values) * (p_survival**n_values) * ((1 - p_survival)**(N - n_values))
line.set_ydata(prob_theory)
time_text.set_text(f'Time t = {current_time:.2f}')
max_prob = max(np.max(freq_sim), np.max(prob_theory)) if freq_sim.size > 0 else 0.1
ax.set_ylim(0, max_prob * 1.15)
return bar_container.patches + [line, time_text]
# --- Setup plot aesthetics ---
ax.set_xlabel('Number of Remaining Particles n', fontsize=12)
ax.set_ylabel('Probability P(n, t)', fontsize=12)
ax.set_title(f'Probability Distribution Evolution for Linear Death Process (N0={N})', fontsize=16)
ax.legend()
ax.set_xlim(-1, N + 1)
# --- Generate and save the animation ---
ani_death = FuncAnimation(fig, update_death_process, frames=len(times), interval=100, blit=True)
ani_death.save('linear_death_distribution.gif', writer=PillowWriter(fps=10))
plt.show()
This animation displays the core process described by the Master Equation—the evolution of the probability distribution.
Initial State: At \(t=0\), the probability distribution is a sharp peak (delta function), with all probability concentrated at the initial particle number \(N=50\).
Evolution Process: As time progresses, this probability peak begins to shift left (mean particle number decreases) and gradually widens (uncertainty increases due to randomness). One can see how probability "flows" from high-particle-number states to low-particle-number states.
Agreement Between Theory and Simulation: At every step of the animation, the histogram (cyan) from tens of thousands of random simulations perfectly matches the theoretical binomial curve (red) we derived using the spectral method. This strongly demonstrates the correctness of the theory.
Convergence to Absorbing State: When time is sufficiently long, the entire probability distribution ultimately concentrates completely at the absorbing state \(n=0\), corresponding to the complete extinction of the population.
Through this dynamic visualization, we are no longer observing a single random trajectory but rather viewing the ensemble behavior of the entire stochastic process—this is exactly what the Master Equation and spectral methods aim to describe.
Summary: Toward Nonlinear Stochastic Field Theory¶
This lecture has deeply explored two core themes in nonequilibrium statistical physics: physical models and mathematical methods. By combining theory with code practice, we have deepened our understanding of these concepts.
1) Physics of Absorbing-State Phase Transitions: Through mean-field analysis of the Directed Percolation (DP) model, we have revealed the basic characteristics of an important class of nonequilibrium phase transitions—absorbing-state phase transitions. The key physical phenomena include the existence of a continuous transition point from active to absorbing phase, and the Critical Slowing Down and Power-Law Dynamics exhibited near the critical point. More importantly, the proposal of the DP universality conjecture suggests that there are profound universal laws behind such phase transitions, motivating physicists to seek more precise theories beyond mean-field theory that can describe fluctuations.
2) Mathematical Methods of Master Equations: To handle discrete stochastic processes, we introduced the Master Equation as a fundamental tool. For simple linear master equations like the linear death process, this lecture demonstrated how to obtain exact solutions through Generating Functions and Spectral Methods (particularly borrowing operator algebra from quantum mechanics). This method is not only elegant, but its results have clear physical interpretation, perfectly connecting abstract mathematical derivation with underlying random events.
However, this lecture has also revealed future challenges. The reason why the spectral method was successfully applied to the linear death process lies in the fact that its Hamiltonian operator \(\mathcal{H}\) is linear (first-order in \(\frac{\partial}{\partial x}\)). For systems like DP that contain nonlinear terms (such as \(2A \to A\) corresponding to \(n^2\) terms), the master equation will lead to a more complex operator containing higher-order derivatives (such as \(\frac{\partial^2}{\partial x^2}\)) in generating function space, making spectral methods usually difficult to apply directly.
For continuous variable systems, a mature pathway has been established from Langevin equations to Fokker-Planck equations to path integrals. For discrete variable systems, starting from the master equation, although it can be approximated as a Fokker-Planck equation through Kramers-Moyal expansion and then construct a path integral, this approximate pathway loses the discreteness information of variables and is not always applicable to systems with strong noise or far from equilibrium.
A more fundamental and powerful approach is: to directly construct a path integral representation from the master equation for discrete stochastic processes. This is precisely the topic of the next lecture. We will introduce a technique called Doi-Peliti path integral, which can directly transform discrete master equations into an effective action of continuous field theory. This framework can not only naturally handle nonlinear interactions, but also systematically incorporate fluctuation effects, serving as the theoretical foundation for using modern field theory tools such as Renormalization Group to study Directed Percolation and other nonequilibrium critical phenomena.