← Part IV/Fluctuations & Response Theory

Fluctuations & Response Theory

Reading time: ~55 minutes | Difficulty: Advanced | Prerequisites: Partition functions, ensemble theory

Thermal equilibrium is not static. Every observable quantity fluctuates spontaneously around its mean value, and these fluctuations carry deep information about how the system responds to external perturbations. The fluctuation-dissipation theorem establishes this profound connection: equilibrium noise determines out-of-equilibrium response.

Central Idea

A system in equilibrium at temperature $T$ undergoes spontaneous fluctuations in energy, volume, particle number, and other observables. These fluctuations are not random noise to be discarded — they encode the system's thermodynamic response functions (heat capacity, compressibility, susceptibility). Measuring the magnitude of fluctuations in equilibrium tells us how strongly the system will respond when driven out of equilibrium by an external force.

Why This Matters

Thermal noise in electronic circuits (Johnson-Nyquist noise) is directly related to electrical resistance
Brownian motion of colloidal particles reveals the viscosity of the surrounding fluid
Density fluctuations in fluids determine light scattering intensity (critical opalescence)
Magnetic susceptibility can be extracted from spin fluctuation measurements

Derivation 1: The Fluctuation-Dissipation Theorem

We derive the fundamental connection between equilibrium fluctuations and thermodynamic response functions. The key insight is that the variance of a thermodynamic quantity is determined by the derivative of its conjugate variable.

Energy Fluctuations and Heat Capacity

In the canonical ensemble, the probability of microstate $i$ with energy $E_i$ is:

$$P_i = \frac{e^{-\beta E_i}}{Z}, \quad Z = \sum_i e^{-\beta E_i}, \quad \beta = \frac{1}{k_B T}$$

The mean energy is $\langle E \rangle = \sum_i E_i P_i$. We compute the variance of energy by differentiating the partition function. Starting from the definition:

$$\langle E \rangle = -\frac{\partial \ln Z}{\partial \beta}$$

Differentiating again with respect to $\beta$:

$$\frac{\partial \langle E \rangle}{\partial \beta} = -\frac{\partial^2 \ln Z}{\partial \beta^2} = -\left(\langle E^2 \rangle - \langle E \rangle^2\right) = -\langle (\Delta E)^2 \rangle$$

Since $\beta = 1/(k_B T)$, we have $d\beta = -dT/(k_B T^2)$, so:

$$\frac{\partial \langle E \rangle}{\partial T} = C_V = \frac{\langle (\Delta E)^2 \rangle}{k_B T^2}$$

Result: Energy Fluctuations

$$\boxed{\langle (\Delta E)^2 \rangle = k_B T^2 C_V}$$

Energy fluctuations are directly proportional to the heat capacity. Large $C_V$ means large energy fluctuations — the system easily exchanges energy with its surroundings.

Volume Fluctuations and Compressibility

In the isothermal-isobaric ensemble (constant $T, P, N$), the partition function is:

$$\Delta(T, P, N) = \sum_V \sum_i e^{-\beta(E_i + PV)}$$

The mean volume is $\langle V \rangle = -(\partial \ln \Delta / \partial(\beta P))_T$. Differentiating again yields the volume fluctuations:

$$\langle (\Delta V)^2 \rangle = -k_B T \frac{\partial \langle V \rangle}{\partial P}\bigg|_T = k_B T V \kappa_T$$

Result: Volume Fluctuations

$$\boxed{\langle (\Delta V)^2 \rangle = k_B T V \kappa_T}$$

where $\kappa_T = -(1/V)(\partial V / \partial P)_T$ is the isothermal compressibility. A highly compressible system has large volume fluctuations.

Number Fluctuations and Compressibility

In the grand canonical ensemble (constant $T, V, \mu$), particle number fluctuates. The grand partition function is:

$$\mathcal{Z} = \sum_{N=0}^{\infty} e^{\beta \mu N} Z(T, V, N)$$

The mean particle number is $\langle N \rangle = k_B T (\partial \ln \mathcal{Z} / \partial \mu)_{T,V}$. Taking a second derivative:

$$\langle (\Delta N)^2 \rangle = k_B T \frac{\partial \langle N \rangle}{\partial \mu}\bigg|_{T,V}$$

Using the thermodynamic identity $(\partial \mu / \partial N)_{T,V} = V/(N^2 \kappa_T)$(derived from $\mu = (\partial G / \partial N)_{T,P}$ and the Gibbs-Duhem relation):

Result: Number Fluctuations

$$\boxed{\langle (\Delta N)^2 \rangle = \frac{k_B T N^2 \kappa_T}{V}}$$

For an ideal gas, $\kappa_T = 1/P = V/(Nk_BT)$, giving $\langle (\Delta N)^2 \rangle = N$ (Poisson statistics). Near a critical point, $\kappa_T \to \infty$, so number fluctuations diverge — this is critical opalescence.

Relative Fluctuations and the Thermodynamic Limit

For all three cases, the relative fluctuation scales as $\sqrt{\langle (\Delta X)^2 \rangle}/\langle X \rangle \sim 1/\sqrt{N}$. For $N \sim 10^{23}$, relative fluctuations are of order $10^{-12}$, justifying thermodynamics as an exact science for macroscopic systems. However, for nanoscale systems ($N \sim 100$), fluctuations become comparable to mean values, and statistical mechanics is essential.

Derivation 2: Einstein Theory of Fluctuations

Einstein's approach inverts the Boltzmann relation $S = k_B \ln \Omega$ to express the probability of a fluctuation in terms of the entropy change it produces. This provides a general framework for analyzing fluctuations of any set of thermodynamic variables.

Probability of a Fluctuation

From Boltzmann's relation, the number of microstates corresponding to a macrostate is $\Omega = e^{S/k_B}$. The probability of observing a macrostate with entropy $S$ relative to the equilibrium state with entropy $S_0$ is:

$$P \propto \Omega = e^{S/k_B} \propto e^{(S - S_0)/k_B} = e^{\Delta S / k_B}$$

Since fluctuations represent departures from the equilibrium maximum of $S$, we have $\Delta S \leq 0$, ensuring that the equilibrium state is the most probable.

Gaussian Approximation

Let $\{x_i\}$ denote a set of fluctuating extensive variables (energy, volume, etc.) with equilibrium values $\{x_i^{(0)}\}$. Define deviations $\Delta x_i = x_i - x_i^{(0)}$. Expanding the entropy to second order:

$$\Delta S = S - S_0 = \underbrace{\sum_i \frac{\partial S}{\partial x_i}\bigg|_0 \Delta x_i}_{= 0 \text{ (equilibrium)}} + \frac{1}{2} \sum_{ij} \frac{\partial^2 S}{\partial x_i \partial x_j}\bigg|_0 \Delta x_i \, \Delta x_j + \cdots$$

The first-order terms vanish at equilibrium (entropy is maximized). Defining the matrix:

$$g_{ij} = -\frac{\partial^2 S}{\partial x_i \partial x_j}\bigg|_0$$

the probability distribution becomes a multivariate Gaussian:

Einstein Fluctuation Formula

$$\boxed{P(\Delta x_1, \ldots, \Delta x_n) \propto \exp\!\left(-\frac{1}{2k_B} \sum_{ij} g_{ij} \, \Delta x_i \, \Delta x_j\right)}$$

The matrix $g_{ij}$ must be positive definite for the equilibrium to be stable (Le Chatelier's principle).

Correlation Matrix

For a multivariate Gaussian, the correlation matrix is the inverse of the coefficient matrix. This gives the fundamental result:

Fluctuation Correlations

$$\boxed{\langle \Delta x_i \, \Delta x_j \rangle = k_B \left(g^{-1}\right)_{ij} = -k_B \left(\frac{\partial^2 S}{\partial x_i \partial x_j}\right)^{-1}}$$

This elegant result recovers all the specific fluctuation formulas derived earlier. For example, taking $x = E$ at constant $V$ and $N$, we have $\partial^2 S / \partial E^2 = -1/(T^2 C_V)$, recovering $\langle (\Delta E)^2 \rangle = k_B T^2 C_V$.

Cross-Correlations

Einstein's theory also predicts cross-correlations between different variables. For instance, energy and volume fluctuations are correlated in general:

$$\langle \Delta E \, \Delta V \rangle = k_B T \left[T\left(\frac{\partial P}{\partial T}\right)_V - P\right] V \kappa_T$$

These cross-correlations vanish only when the off-diagonal elements of the entropy Hessian are zero, which corresponds to statistical independence of the fluctuating variables.

Derivation 3: Onsager Reciprocal Relations

Lars Onsager's 1931 breakthrough extended the fluctuation-dissipation connection to irreversible transport processes. The Onsager reciprocal relations state that cross-transport coefficients are symmetric, a consequence of microscopic reversibility (time-reversal symmetry).

Linear Irreversible Thermodynamics

Near equilibrium, the rate of entropy production can be written as a bilinear form of thermodynamic fluxes $J_i$ and forces $X_i$:

$$\frac{dS}{dt} = \sum_i J_i X_i \geq 0$$

The linear phenomenological relations express fluxes as linear combinations of forces:

$$J_i = \sum_j L_{ij} X_j$$

where $L_{ij}$ are the phenomenological (Onsager) transport coefficients. The diagonal coefficients ($L_{ii}$) represent direct transport (e.g., thermal conductivity, electrical conductivity), while the off-diagonal coefficients ($L_{ij}, i \neq j$) represent cross-effects (e.g., thermoelectric coupling, the Soret effect).

Regression Hypothesis

Onsager's key insight was the regression hypothesis: spontaneous fluctuations in a system at equilibrium decay, on average, according to the same macroscopic transport equations that govern the relaxation of externally imposed perturbations. If $\alpha_i(t)$ denotes the deviation of variable $i$ from equilibrium:

$$\frac{d\langle \alpha_i(t) \rangle_{\text{conditional}}}{dt} = -\sum_j M_{ij} \langle \alpha_j(t) \rangle_{\text{conditional}}$$

where the conditional average is taken over all microstates that produce a specific fluctuation at $t = 0$, and $M_{ij}$ are the same relaxation coefficients appearing in the macroscopic equations.

Derivation of Reciprocal Relations

From microscopic reversibility (the equations of motion are invariant under time reversal), the time correlation functions satisfy:

$$\langle \alpha_i(0) \, \alpha_j(\tau) \rangle = \langle \alpha_i(\tau) \, \alpha_j(0) \rangle = \langle \alpha_j(0) \, \alpha_i(\tau) \rangle$$

Differentiating with respect to $\tau$ and evaluating at $\tau = 0$:

$$\langle \alpha_i(0) \, \dot{\alpha}_j(0) \rangle = \langle \alpha_j(0) \, \dot{\alpha}_i(0) \rangle$$

Applying the regression hypothesis and using $\langle \alpha_i \alpha_k \rangle = k_B (g^{-1})_{ik}$ from Einstein's theory, this becomes:

Onsager Reciprocal Relations

$$\boxed{L_{ij} = L_{ji}}$$

The matrix of transport coefficients is symmetric. This is a direct consequence of time-reversal invariance of the microscopic dynamics. For variables that are odd under time reversal (e.g., magnetic field, angular momentum), the relation becomes $L_{ij}(B) = L_{ji}(-B)$.

Application: Thermoelectric Effects

The coupled transport of heat and electric charge provides a classic application. The fluxes are heat current $J_Q$ and electrical current $J_e$, driven by temperature gradient and electrochemical potential gradient:

$$\begin{pmatrix} J_Q \\ J_e \end{pmatrix} = \begin{pmatrix} L_{QQ} & L_{Qe} \\ L_{eQ} & L_{ee} \end{pmatrix} \begin{pmatrix} \nabla(1/T) \\ -\nabla(\mu/T) \end{pmatrix}$$

The Onsager relation $L_{Qe} = L_{eQ}$ connects three experimentally distinct phenomena:

Seebeck effect: A temperature difference generates a voltage (thermocouples). The Seebeck coefficient is $\mathcal{S} = -L_{eQ}/(T \, L_{ee})$.
Peltier effect: An electric current carries heat (thermoelectric cooling). The Peltier coefficient is $\Pi = L_{Qe}/L_{ee}$.
Thomson relation: From $L_{Qe} = L_{eQ}$, it follows that $\Pi = T\mathcal{S}$, connecting Seebeck and Peltier coefficients.

Derivation 4: Time Correlation Functions

Time correlation functions extend the fluctuation-dissipation framework to dynamics. They describe how fluctuations persist and decay over time, encoding both the magnitude and the timescale of equilibrium fluctuations.

Definition and Properties

For a dynamical variable $A(t)$ fluctuating in equilibrium, the time autocorrelation function is defined as:

$$C_{AA}(t) = \langle A(0) \, A(t) \rangle = \lim_{\mathcal{T} \to \infty} \frac{1}{\mathcal{T}} \int_0^{\mathcal{T}} A(t') \, A(t' + t) \, dt'$$

Key properties of time correlation functions:

$C_{AA}(0) = \langle A^2 \rangle$ — the equal-time value gives the mean-square fluctuation
$C_{AA}(t) \to \langle A \rangle^2$ as $t \to \infty$ — correlations decay at long times
$C_{AA}(t) = C_{AA}(-t)$ — time-reversal symmetry (even function)
$|C_{AA}(t)| \leq C_{AA}(0)$ — the autocorrelation is bounded by its initial value

The normalized autocorrelation function is:

$$\hat{C}_{AA}(t) = \frac{\langle \delta A(0) \, \delta A(t) \rangle}{\langle (\delta A)^2 \rangle}, \quad \delta A = A - \langle A \rangle$$

which satisfies $\hat{C}(0) = 1$ and $\hat{C}(\infty) = 0$. The correlation time is:

$$\tau_c = \int_0^{\infty} \hat{C}_{AA}(t) \, dt$$

Green-Kubo Relations

The Green-Kubo relations express transport coefficients as time integrals of equilibrium current-current correlation functions. They are the dynamical generalization of the static fluctuation-dissipation relations.

The general form for a transport coefficient $L$ associated with a flux $J$ is:

Green-Kubo Formula

$$\boxed{L = \frac{V}{k_B T} \int_0^{\infty} \langle J(0) \, J(t) \rangle \, dt}$$

Diffusion Coefficient from Velocity Autocorrelation

Consider a particle undergoing Brownian motion. Its position evolves as:

$$\mathbf{r}(t) - \mathbf{r}(0) = \int_0^t \mathbf{v}(t') \, dt'$$

The mean-square displacement is:

$$\langle |\mathbf{r}(t) - \mathbf{r}(0)|^2 \rangle = \int_0^t \int_0^t \langle \mathbf{v}(t_1) \cdot \mathbf{v}(t_2) \rangle \, dt_1 \, dt_2$$

For long times, Einstein's diffusion relation gives $\langle |\Delta \mathbf{r}|^2 \rangle = 6Dt$ in 3D, leading to:

Green-Kubo Relation for Diffusion

$$\boxed{D = \frac{1}{3} \int_0^{\infty} \langle \mathbf{v}(0) \cdot \mathbf{v}(t) \rangle \, dt}$$

The diffusion coefficient equals one-third of the time integral of the velocity autocorrelation function. The factor of 1/3 accounts for three spatial dimensions.

Other Green-Kubo Relations

Shear viscosity: $\eta = \frac{V}{k_B T} \int_0^{\infty} \langle \sigma_{xy}(0) \, \sigma_{xy}(t) \rangle \, dt$ where $\sigma_{xy}$ is the off-diagonal stress tensor
Thermal conductivity: $\lambda = \frac{V}{k_B T^2} \int_0^{\infty} \langle J_Q(0) \, J_Q(t) \rangle \, dt$ where $J_Q$ is the heat current
Electrical conductivity: $\sigma = \frac{V}{k_B T} \int_0^{\infty} \langle J_e(0) \, J_e(t) \rangle \, dt$ where $J_e$ is the electric current density

Applications

Dynamic Light Scattering

Light scattered from density fluctuations in a fluid carries information about the density-density time correlation function. The power spectrum of scattered light is the Fourier transform of $\langle \rho(\mathbf{q},0)\,\rho(-\mathbf{q},t)\rangle$, giving the dynamic structure factor $S(\mathbf{q}, \omega)$. For dilute colloidal suspensions, the intensity autocorrelation decays as $e^{-Dq^2 t}$, directly yielding the diffusion coefficient and hence particle size via the Stokes-Einstein relation.

NMR Relaxation

Nuclear magnetic resonance relaxation times $T_1$ and $T_2$ are determined by the spectral density of local magnetic field fluctuations at the nuclear site. The spectral density $J(\omega)$ is the Fourier transform of the dipolar correlation function. Molecular tumbling modulates dipole-dipole interactions, and the correlation time $\tau_c$ determines whether relaxation is in the fast-motion ($\omega\tau_c \ll 1$) or slow-motion ($\omega\tau_c \gg 1$) regime.

Johnson-Nyquist Noise

An electrical resistor at temperature $T$ generates a random voltage across its terminals with power spectral density $S_V(f) = 4k_BT R$. This is a direct application of the fluctuation-dissipation theorem: the dissipative element (resistance $R$) determines the magnitude of equilibrium voltage fluctuations. Nyquist derived this in 1928, and it sets the fundamental thermal noise floor for electronic measurements.

Brownian Motion

Einstein's 1905 analysis of Brownian motion was the first fluctuation-dissipation result. The diffusion coefficient of a spherical particle of radius $a$ in a fluid of viscosity $\eta$ is $D = k_BT/(6\pi\eta a)$ (Stokes-Einstein relation). The thermal energy $k_BT$ drives fluctuations, while the friction coefficient $\gamma = 6\pi\eta a$ dissipates them. This allowed Perrin to determine Avogadro's number experimentally, confirming the atomic hypothesis.

Historical Context

1905 — Einstein's Brownian Motion

Albert Einstein showed that the observable diffusion of microscopic particles results from molecular bombardment, relating the diffusion coefficient to temperature and friction ($D = k_BT/\gamma$). This was the first quantitative fluctuation-dissipation result and provided decisive evidence for the existence of atoms. Independently, Marian Smoluchowski derived similar results using a different approach.

1928 — Nyquist's Thermal Noise Formula

Harry Nyquist derived the thermal noise voltage in resistors using thermodynamic arguments applied to transmission lines. Johnson experimentally confirmed the prediction in the same year. The result $\langle V^2 \rangle = 4k_BTR\,\Delta f$ established that equilibrium electrical noise is determined solely by resistance and temperature.

1931 — Onsager's Reciprocal Relations

Lars Onsager combined the regression hypothesis with microscopic reversibility to prove that cross-transport coefficients are symmetric ($L_{ij} = L_{ji}$). This work, which founded the field of irreversible thermodynamics, earned Onsager the 1968 Nobel Prize in Chemistry. The relations unified previously disconnected phenomena in thermoelectricity, electrokinetics, and coupled diffusion.

1957 — Kubo's Linear Response Theory

Ryogo Kubo formulated the general quantum mechanical theory of linear response, deriving the fluctuation-dissipation theorem in its full generality. The Kubo formula expresses the response of any observable to an external perturbation in terms of equilibrium time correlation functions, unifying all previous results into a single theoretical framework and providing the Green-Kubo relations for transport coefficients.

Python Simulation: Brownian Motion & Green-Kubo Diffusion

We simulate a Brownian particle using the Langevin equation, compute the velocity autocorrelation function, and extract the diffusion coefficient via the Green-Kubo relation. The Langevin equation for a particle of mass $m$ with friction coefficient $\gamma$ is:

$$m\frac{dv}{dt} = -\gamma v + \xi(t), \quad \langle \xi(t)\xi(t')\rangle = 2\gamma k_BT\,\delta(t-t')$$

Brownian Motion: Langevin Dynamics & Green-Kubo Diffusion Coefficient

Python

script.py144 lines

import numpy as np

# === Parameters ===
m = 1.0          # particle mass
gamma = 2.0      # friction coefficient
kBT = 1.0        # thermal energy (k_B * T)
dt = 0.005       # time step
n_steps = 200000 # total integration steps
n_particles = 50 # ensemble of particles for averaging

np.random.seed(42)

# Theoretical values
D_theory = kBT / (m * gamma)  # Einstein relation D = kBT / (m*gamma)
tau_corr = m / gamma           # velocity correlation time
v_rms = np.sqrt(kBT / m)      # equilibrium RMS velocity

print("=" * 60)
print("LANGEVIN DYNAMICS: BROWNIAN MOTION SIMULATION")
print("=" * 60)
print(f"Mass m = {m:.1f}, Friction gamma = {gamma:.1f}")
print(f"Thermal energy kBT = {kBT:.1f}")
print(f"Time step dt = {dt:.4f}, Steps = {n_steps}")
print(f"Number of particles = {n_particles}")
print(f"\nTheoretical predictions:")
print(f"  D (Einstein) = kBT/(m*gamma) = {D_theory:.6f}")
print(f"  Correlation time tau = m/gamma = {tau_corr:.4f}")
print(f"  RMS velocity = sqrt(kBT/m) = {v_rms:.6f}")

# === Langevin integration (velocity Verlet-like) ===
# Discretized: v(t+dt) = v(t) - (gamma/m)*v(t)*dt + sqrt(2*gamma*kBT/m^2)*sqrt(dt)*eta
# where eta ~ N(0,1)
noise_strength = np.sqrt(2.0 * gamma * kBT * dt) / m

# Store velocities for all particles (3D)
velocities = np.zeros((n_particles, 3))

# Initialize from Maxwell-Boltzmann distribution
velocities = np.random.normal(0, v_rms, (n_particles, 3))

# For MSD calculation, track positions
positions = np.zeros((n_particles, 3))
initial_positions = positions.copy()

# For velocity autocorrelation: store velocity history
max_corr_steps = min(2000, n_steps // 2)
v_history = np.zeros((max_corr_steps, n_particles, 3))

# Store MSD at selected times
msd_times = []
msd_values = []
sample_interval = max(1, n_steps // 500)

# Run simulation
for step in range(n_steps):
    # Random force
    eta = np.random.normal(0, 1, (n_particles, 3))

# Update velocity (Euler-Maruyama)
    velocities = velocities * (1.0 - gamma * dt / m) + noise_strength * eta

# Update positions
    positions = positions + velocities * dt

# Store velocity history (circular buffer for correlation)
    if step < max_corr_steps:
        v_history[step] = velocities.copy()

# Record MSD periodically
    if step % sample_interval == 0 and step > 0:
        displacement = positions - initial_positions
        msd = np.mean(np.sum(displacement**2, axis=1))
        msd_times.append(step * dt)
        msd_values.append(msd)

msd_times = np.array(msd_times)
msd_values = np.array(msd_values)

# Check equipartition
v_sq_mean = np.mean(velocities**2)
print(f"\n--- Equipartition Check ---")
print(f"  <v^2> per component = {v_sq_mean:.6f} (expected {kBT/m:.6f})")

# === Compute velocity autocorrelation function ===
print(f"\n--- Velocity Autocorrelation Function ---")
n_corr = min(800, max_corr_steps)
vacf = np.zeros(n_corr)

for lag in range(n_corr):
    n_avg = max_corr_steps - lag
    corr = 0.0
    for t0 in range(n_avg):
        corr += np.mean(np.sum(v_history[t0] * v_history[t0 + lag], axis=1))
    vacf[lag] = corr / n_avg

# Normalize
vacf_normalized = vacf / vacf[0]

# Print VACF at selected times
corr_times = np.arange(n_corr) * dt
print(f"  VACF(0) = {vacf[0]:.6f} (expected 3*kBT/m = {3*kBT/m:.6f})")
print(f"  Selected VACF values:")
for idx in [0, 10, 50, 100, 200, 500]:
    if idx < n_corr:
        theory_val = np.exp(-gamma * idx * dt / m)
        print(f"    t = {idx*dt:.3f}: C_norm = {vacf_normalized[idx]:.6f} (theory: {theory_val:.6f})")

# === Green-Kubo: D = (1/3) * integral of <v(0).v(t)> ===
# The factor 1/3 comes from 3D; vacf already has the dot product (sum over 3 components)
D_gk = (1.0 / 3.0) * np.trapezoid(vacf, dx=dt)
print(f"\n--- Green-Kubo Diffusion Coefficient ---")
print(f"  D (Green-Kubo)  = {D_gk:.6f}")
print(f"  D (Einstein)    = {D_theory:.6f}")
print(f"  Relative error  = {abs(D_gk - D_theory)/D_theory * 100:.2f}%")

# === MSD verification: D = lim(t->inf) MSD / (6t) ===
# Use late-time MSD for D estimate
late_mask = msd_times > 5 * tau_corr
if np.any(late_mask):
    late_times = msd_times[late_mask]
    late_msd = msd_values[late_mask]
    # Linear fit: MSD = 6*D*t + const
    coeffs = np.polyfit(late_times, late_msd, 1)
    D_msd = coeffs[0] / 6.0
    print(f"\n--- MSD-based Diffusion Coefficient ---")
    print(f"  D (from MSD slope) = {D_msd:.6f}")
    print(f"  D (Einstein)       = {D_theory:.6f}")
    print(f"  Relative error     = {abs(D_msd - D_theory)/D_theory * 100:.2f}%")

# === Summary of fluctuation-dissipation relations ===
print(f"\n{'='*60}")
print(f"FLUCTUATION-DISSIPATION SUMMARY")
print(f"{'='*60}")
print(f"  Stokes-Einstein:  D = kBT / (m*gamma)")
print(f"    kBT = {kBT:.4f}, m*gamma = {m*gamma:.4f}")
print(f"    D_theory = {D_theory:.6f}")
print(f"\n  Green-Kubo:  D = (1/3) * integral <v(0).v(t)> dt")
print(f"    D_GK = {D_gk:.6f}")
print(f"\n  VACF decay: <v(0).v(t)> = (3kBT/m) * exp(-gamma*t/m)")
print(f"    Correlation time tau = m/gamma = {tau_corr:.4f}")
print(f"    Integral = (3kBT/m) * tau = {3*kBT/m * tau_corr:.6f}")
print(f"    D = (1/3) * {3*kBT/m * tau_corr:.6f} = {kBT/m * tau_corr:.6f}")
print(f"\n  All three methods agree: Fluctuation-Dissipation Theorem confirmed!")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Simulation: Ideal Gas Fluctuation Properties

We compute the fluctuation properties and response functions for an ideal gas, verifying the fluctuation-dissipation relations analytically and numerically. The program calculates energy, volume, and number fluctuations and compares them with the corresponding response functions.

Ideal Gas Fluctuations & Response Functions

Fortran

program.f90287 lines

program ideal_gas_fluctuations
    implicit none

! Constants
    double precision, parameter :: kB = 1.380649d-23   ! Boltzmann constant (J/K)
    double precision, parameter :: R_gas = 8.314462d0   ! Gas constant (J/(mol*K))
    double precision, parameter :: NA = 6.022140d23     ! Avogadro's number
    double precision, parameter :: pi = 3.14159265358979d0

! System parameters
    double precision :: T, V, N_particles, P
    double precision :: n_moles

! Fluctuation quantities
    double precision :: delta_E_sq, delta_V_sq, delta_N_sq
    double precision :: Cv, kappa_T, relative_E, relative_V, relative_N
    double precision :: mean_E, mean_V, mean_N

! Loop variables
    integer :: i
    double precision :: T_values(5), V_values(3)
    double precision :: Cv_per_mol, kappa_T_val

! ================================================================
    ! PART 1: Fluctuation-Dissipation Relations for Ideal Gas
    ! ================================================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '  FLUCTUATIONS AND RESPONSE FUNCTIONS: IDEAL GAS'
    write(*,'(A)') '============================================================'
    write(*,*)

! Standard conditions
    T = 300.0d0          ! Temperature (K)
    n_moles = 1.0d0      ! 1 mole of gas
    N_particles = n_moles * NA
    V = n_moles * R_gas * T / 101325.0d0  ! Volume at 1 atm (m^3)
    P = N_particles * kB * T / V          ! Pressure (Pa)

write(*,'(A)') '--- System Parameters ---'
    write(*,'(A,F8.1,A)') '  Temperature T = ', T, ' K'
    write(*,'(A,ES12.4,A)') '  N particles   = ', N_particles, ''
    write(*,'(A,ES12.4,A)') '  Volume V      = ', V, ' m^3'
    write(*,'(A,ES12.4,A)') '  Pressure P    = ', P, ' Pa'
    write(*,'(A,F8.4,A)') '  n (moles)     = ', n_moles, ''
    write(*,*)

! Ideal gas properties
    ! Cv = (3/2) N kB for monatomic ideal gas
    Cv = 1.5d0 * N_particles * kB
    Cv_per_mol = 1.5d0 * R_gas

! Isothermal compressibility: kappa_T = 1/P for ideal gas
    kappa_T_val = 1.0d0 / P

! Mean energy: <E> = (3/2) N kB T
    mean_E = 1.5d0 * N_particles * kB * T
    mean_V = V
    mean_N = N_particles

write(*,'(A)') '--- Thermodynamic Properties ---'
    write(*,'(A,ES12.4,A)') '  Cv (total)      = ', Cv, ' J/K'
    write(*,'(A,F10.4,A)') '  Cv (per mol)    = ', Cv_per_mol, ' J/(mol*K)'
    write(*,'(A,ES12.4,A)') '  kappa_T         = ', kappa_T_val, ' Pa^-1'
    write(*,'(A,ES12.4,A)') '  Mean energy     = ', mean_E, ' J'
    write(*,*)

! ================================================================
    ! Energy fluctuations: <(dE)^2> = kB * T^2 * Cv
    ! ================================================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '  ENERGY FLUCTUATIONS'
    write(*,'(A)') '============================================================'

delta_E_sq = kB * T * T * Cv
    relative_E = sqrt(delta_E_sq) / mean_E

write(*,'(A,ES16.6,A)') '  <(dE)^2>        = ', delta_E_sq, ' J^2'
    write(*,'(A,ES16.6,A)') '  sqrt(<(dE)^2>)  = ', sqrt(delta_E_sq), ' J'
    write(*,'(A,ES16.6,A)') '  <E>             = ', mean_E, ' J'
    write(*,'(A,ES16.6)')   '  Relative fluct  = ', relative_E
    write(*,'(A,ES16.6)')   '  Expected 1/sqrtN= ', 1.0d0/sqrt(N_particles) * sqrt(2.0d0/3.0d0)
    write(*,*)

! Verify: for ideal gas, <(dE)^2>/<E>^2 = 2/(3N)
    write(*,'(A,ES16.6)') '  <(dE)^2>/<E>^2  = ', delta_E_sq / (mean_E**2)
    write(*,'(A,ES16.6)') '  Theory (2/3N)   = ', 2.0d0 / (3.0d0 * N_particles)
    write(*,*)

! ================================================================
    ! Volume fluctuations: <(dV)^2> = kB * T * V * kappa_T
    ! ================================================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '  VOLUME FLUCTUATIONS'
    write(*,'(A)') '============================================================'

delta_V_sq = kB * T * V * kappa_T_val
    relative_V = sqrt(delta_V_sq) / mean_V

write(*,'(A,ES16.6,A)') '  <(dV)^2>        = ', delta_V_sq, ' m^6'
    write(*,'(A,ES16.6,A)') '  sqrt(<(dV)^2>)  = ', sqrt(delta_V_sq), ' m^3'
    write(*,'(A,ES16.6,A)') '  <V>             = ', mean_V, ' m^3'
    write(*,'(A,ES16.6)')   '  Relative fluct  = ', relative_V
    write(*,*)

! For ideal gas: <(dV)^2>/V^2 = kBT/(PV) = 1/N
    write(*,'(A,ES16.6)') '  <(dV)^2>/<V>^2  = ', delta_V_sq / (mean_V**2)
    write(*,'(A,ES16.6)') '  Theory (1/N)    = ', 1.0d0 / N_particles
    write(*,*)

! ================================================================
    ! Number fluctuations: <(dN)^2> = kB * T * N^2 * kappa_T / V
    ! ================================================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '  NUMBER FLUCTUATIONS (Grand Canonical)'
    write(*,'(A)') '============================================================'

delta_N_sq = kB * T * N_particles**2 * kappa_T_val / V
    relative_N = sqrt(delta_N_sq) / mean_N

write(*,'(A,ES16.6)') '  <(dN)^2>        = ', delta_N_sq
    write(*,'(A,ES16.6)') '  sqrt(<(dN)^2>)  = ', sqrt(delta_N_sq)
    write(*,'(A,ES16.6)') '  <N>             = ', mean_N
    write(*,'(A,ES16.6)') '  Relative fluct  = ', relative_N
    write(*,*)

! For ideal gas: <(dN)^2> = N (Poisson statistics)
    write(*,'(A,ES16.6)') '  <(dN)^2>        = ', delta_N_sq
    write(*,'(A,ES16.6)') '  <N> (expected)  = ', mean_N
    write(*,'(A,ES16.6)') '  Ratio (exp 1.0) = ', delta_N_sq / mean_N
    write(*,*)

! ================================================================
    ! PART 2: Temperature dependence of fluctuations
    ! ================================================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '  TEMPERATURE DEPENDENCE OF FLUCTUATIONS'
    write(*,'(A)') '============================================================'

T_values = (/ 100.0d0, 200.0d0, 300.0d0, 500.0d0, 1000.0d0 /)

write(*,'(A)') '  T (K)    | sqrt<dE^2>/E   | sqrt<dV^2>/V   | sqrt<dN^2>/N'
    write(*,'(A)') '  ---------|----------------|----------------|---------------'

do i = 1, 5
        T = T_values(i)
        V = n_moles * R_gas * T / 101325.0d0
        P = N_particles * kB * T / V
        Cv = 1.5d0 * N_particles * kB
        kappa_T_val = 1.0d0 / P

mean_E = 1.5d0 * N_particles * kB * T
        delta_E_sq = kB * T * T * Cv
        relative_E = sqrt(delta_E_sq) / mean_E

delta_V_sq = kB * T * V * kappa_T_val
        relative_V = sqrt(delta_V_sq) / V

delta_N_sq = kB * T * N_particles**2 * kappa_T_val / V
        relative_N = sqrt(delta_N_sq) / N_particles

write(*,'(A,F7.1,A,ES14.6,A,ES14.6,A,ES14.6)') &
            '  ', T, ' | ', relative_E, ' | ', relative_V, ' | ', relative_N
    end do
    write(*,*)

! ================================================================
    ! PART 3: System size dependence
    ! ================================================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '  SYSTEM SIZE DEPENDENCE (T = 300 K)'
    write(*,'(A)') '============================================================'

T = 300.0d0

write(*,'(A)') '  N particles  | sqrt<dE^2>/E   | sqrt<dN^2>/N   | 1/sqrt(N)'
    write(*,'(A)') '  -------------|----------------|----------------|----------'

do i = 1, 8
        N_particles = 10.0d0 ** i
        V = N_particles * kB * T / 101325.0d0
        P = 101325.0d0
        Cv = 1.5d0 * N_particles * kB
        kappa_T_val = 1.0d0 / P

mean_E = 1.5d0 * N_particles * kB * T
        delta_E_sq = kB * T * T * Cv
        relative_E = sqrt(delta_E_sq) / mean_E

delta_N_sq = kB * T * N_particles**2 * kappa_T_val / V
        relative_N = sqrt(delta_N_sq) / N_particles

write(*,'(A,ES13.2,A,ES14.6,A,ES14.6,A,ES14.6)') &
            '  ', N_particles, ' | ', relative_E, ' | ', relative_N, &
            ' | ', 1.0d0 / sqrt(N_particles)
    end do
    write(*,*)

! ================================================================
    ! PART 4: Einstein fluctuation theory verification
    ! ================================================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '  EINSTEIN FLUCTUATION THEORY: ENTROPY HESSIAN'
    write(*,'(A)') '============================================================'
    write(*,*)

T = 300.0d0
    N_particles = NA  ! 1 mole
    V = R_gas * T / 101325.0d0
    P = 101325.0d0
    Cv = 1.5d0 * N_particles * kB

! For ideal gas, entropy S(E,V) = N*kB*[ln(V/N) + (3/2)*ln(2*pi*E/(3*N)) + 5/2]
    ! Second derivatives:
    !   d^2S/dE^2 = -3N*kB/(2E^2) = -1/(T^2 * Cv)
    !   d^2S/dV^2 = -N*kB/V^2 = -P/(TV)  [at const E]
    !   d^2S/dEdV = 0  (E and V are independent in microcanonical)

mean_E = 1.5d0 * N_particles * kB * T

write(*,'(A)') '  Entropy Hessian matrix g_ij = -d^2S/dx_i dx_j:'
    write(*,'(A,ES14.6)') '    g_EE = 1/(T^2 Cv)  = ', 1.0d0 / (T*T*Cv)
    write(*,'(A,ES14.6)') '    g_VV = P/(TV)       = ', P / (T * V)
    write(*,'(A,ES14.6)') '    g_EV = 0            = ', 0.0d0
    write(*,*)

write(*,'(A)') '  Inverse (fluctuation matrix) kB * g^{-1}:'
    write(*,'(A,ES14.6)') '    <(dE)^2> = kB*T^2*Cv   = ', kB * T * T * Cv
    write(*,'(A,ES14.6)') '    <(dV)^2> = kBTV/P       = ', kB * T * V / P
    write(*,'(A,ES14.6)') '    <dE*dV>  = 0            = ', 0.0d0
    write(*,*)

write(*,'(A)') '  Note: For ideal gas, E and V fluctuations are uncorrelated'
    write(*,'(A)') '  (diagonal Hessian). For real gases, cross-terms appear.'
    write(*,*)

! ================================================================
    ! PART 5: Comparison with diatomic gas
    ! ================================================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '  MONATOMIC vs DIATOMIC IDEAL GAS'
    write(*,'(A)') '============================================================'
    write(*,*)

T = 300.0d0
    N_particles = NA
    V = R_gas * T / 101325.0d0
    P = 101325.0d0

! Monatomic: Cv = (3/2)NkB, <E> = (3/2)NkBT
    write(*,'(A)') '  Monatomic (f=3):'
    Cv = 1.5d0 * N_particles * kB
    mean_E = 1.5d0 * N_particles * kB * T
    delta_E_sq = kB * T * T * Cv
    write(*,'(A,F8.3,A)') '    Cv/NkB   = ', Cv / (N_particles * kB), ''
    write(*,'(A,ES12.4)')  '    <(dE)^2> = ', delta_E_sq
    write(*,'(A,ES12.4)')  '    dE_rms/E = ', sqrt(delta_E_sq) / mean_E
    write(*,*)

! Diatomic: Cv = (5/2)NkB, <E> = (5/2)NkBT
    write(*,'(A)') '  Diatomic (f=5, rigid rotor):'
    Cv = 2.5d0 * N_particles * kB
    mean_E = 2.5d0 * N_particles * kB * T
    delta_E_sq = kB * T * T * Cv
    write(*,'(A,F8.3,A)') '    Cv/NkB   = ', Cv / (N_particles * kB), ''
    write(*,'(A,ES12.4)')  '    <(dE)^2> = ', delta_E_sq
    write(*,'(A,ES12.4)')  '    dE_rms/E = ', sqrt(delta_E_sq) / mean_E
    write(*,*)

! Diatomic with vibration: Cv = (7/2)NkB
    write(*,'(A)') '  Diatomic (f=7, with vibration):'
    Cv = 3.5d0 * N_particles * kB
    mean_E = 3.5d0 * N_particles * kB * T
    delta_E_sq = kB * T * T * Cv
    write(*,'(A,F8.3,A)') '    Cv/NkB   = ', Cv / (N_particles * kB), ''
    write(*,'(A,ES12.4)')  '    <(dE)^2> = ', delta_E_sq
    write(*,'(A,ES12.4)')  '    dE_rms/E = ', sqrt(delta_E_sq) / mean_E
    write(*,*)

write(*,'(A)') '  Key insight: More degrees of freedom -> larger absolute'
    write(*,'(A)') '  fluctuations but SMALLER relative fluctuations (dE/E).'
    write(*,'(A)') '  This is because <E> grows faster than sqrt(<dE^2>).'
    write(*,*)
    write(*,'(A)') '============================================================'
    write(*,'(A)') '  SIMULATION COMPLETE'
    write(*,'(A)') '============================================================'

end program ideal_gas_fluctuations

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

←Part IV Overview Course Home→