# An Introduction to Computational
Fluctuating Hydrodynamics

Alejandro L. Garcia  San Jose State University  John B. Bell  Lawrence
Berkeley National Laboratory  Andrew Nonaka  Lawrence Berkeley National
Laboratory  Ishan Srivastava  Lawrence Berkeley National Laboratory  Daniel
Ladiges  Lawrence Berkeley National Laboratory  Changho Kim  University of
California Merced

These notes are an introduction to fluctuating hydrodynamics (FHD) and the
formulation of numerical schemes for the resulting stochastic partial
differential equations (PDEs). Fluctuating hydrodynamics was originally
introduced by Landau and Lifshitz [1] as a way to put thermal fluctuations
into a continuum framework by including a stochastic forcing to each
dissipative transport process (e.g., heat flux). While FHD has been useful in
modeling transport and fluid dynamics at the mesoscopic scale, theoretical
calculations have been feasible only with simplifying assumptions. As such
there is great interest in numerical schemes for Computational Fluctuating
Hydrodynamics (CFHD). There are a variety of algorithms (e.g., spectral,
finite element, lattice Boltzmann) but in this introduction we focus on finite
volume schemes. Accompanying these notes is a demonstration program in Python
available on GitHub.

### Stochastic Heat Equation

Our first example is the stochastic heat equation (SHE) for the transport of
internal energy by thermal diffusion. This section presents a short
derivation, for details see [2]. We start with the continuity equation for
energy

$\frac{\partial}{\partial t}(\rho e)=-\nabla\cdot\boldsymbol{Q}$ (1)

where $\rho$ is the (constant) mass density and $e$ is the specific internal
energy. The heat flux, $\boldsymbol{Q}$, can be separated into deterministic
and stochastic contributions (denoted by $\overline{\boldsymbol{Q}}$ and
$\widetilde{\boldsymbol{Q}}$, respectively) and is written in Onsager form as

$\boldsymbol{Q}=\overline{\boldsymbol{Q}}+\widetilde{\boldsymbol{Q}}=L\boldsymbol{X}+\widetilde{\boldsymbol{Q}}$
(2)

where $L$ is the Onsager coefficient and $\boldsymbol{X}$ is the thermodynamic
force. Fluctuation-dissipation gives that the stochastic flux has zero mean
and correlation

$\langle\widetilde{\boldsymbol{Q}}(\mathbf{r},t)\widetilde{\boldsymbol{Q}}(\mathbf{r}^{\prime},t^{\prime})\rangle=2k_{B}L~{}\delta(\mathbf{r}-\mathbf{r}^{\prime})~{}\delta(t-t^{\prime})$
(3)

where $k_{B}$ is the Boltzmann constant and $\delta(\cdot)$ denotes the Dirac
delta function. Readers familiar with the Langevin equation for Brownian
motion [3] will notice the similarities between that stochastic ordinary
differential equation (ODE) and the above stochastic PDE.

From non-equilibrium thermodynamics [2, 4], the deterministic total rate of
entropy change in a region $\Omega$ is

$\frac{d\overline{S}}{dt}=\int_{\Omega}\frac{d\overline{s}}{dt}\;d\mathbf{r}=\int_{\Omega}\boldsymbol{X}\cdot\boldsymbol{J}\;d\mathbf{r}-\left[\frac{\overline{\boldsymbol{Q}}}{T}\right]_{\partial\Omega}$
(4)

where $\boldsymbol{J}$ is the thermodynamic flux. The first term is the
internal entropy production and the second is the change due to heat flow at
the boundary. Using the Gibbs equation, $d\overline{u}=\rho
d\overline{e}=Td\overline{s}$, and Eq. (1) gives

$\int_{\Omega}\frac{d\overline{s}}{dt}\;d\mathbf{r}=\int_{\Omega}\frac{\rho}{T}\frac{d\overline{e}}{dt}\;d\mathbf{r}=-\int_{\Omega}\frac{1}{T}\nabla\cdot\overline{\boldsymbol{Q}}\;d\mathbf{r}=\int_{\Omega}\frac{1}{T}\nabla\cdot(\lambda\nabla
T)\;d\mathbf{r}.$ (5)

The last equality uses the phenomenological Fourier law for heat flow,
$\overline{\boldsymbol{Q}}=-\lambda\nabla T$, where $\lambda$ is the thermal
conductivity. Applying Green’s first identity gives

$\displaystyle\int_{\Omega}\frac{1}{T}\nabla\cdot(\lambda\nabla
T)\;d\mathbf{r}$
$\displaystyle=\left[-\frac{\overline{\boldsymbol{Q}}}{T}\right]_{\partial\Omega}-\int_{\Omega}\nabla\left(\frac{1}{T}\right)\cdot(\lambda\nabla
T)\;d\mathbf{r}.$ (6)

The entropy expressions (4) and (5) can be equated, and so the chain rule
gives

$\displaystyle\int_{\Omega}\boldsymbol{X}\cdot\boldsymbol{J}\;d\mathbf{r}$
$\displaystyle=-\int_{\Omega}\nabla\left(\frac{1}{T}\right)\cdot(\lambda\nabla
T)\;d\mathbf{r}=\int_{\Omega}\nabla\left(\frac{1}{T}\right)\cdot\left(\lambda
T^{2}\;\nabla\frac{1}{T}\right)\;d\mathbf{r}.$ (7)

From this equality,
$\boldsymbol{X}\cdot\boldsymbol{J}=\boldsymbol{X}\cdot(L\boldsymbol{X})$, and
Eq. (2), we can identity the thermodynamic force $\boldsymbol{X}$, the
thermodynamics flux $\boldsymbol{J}$, and the Onsager coefficient $L$ to be
[4, 5]

$\boldsymbol{X}=\nabla\left(\frac{1}{T}\right)\;\mathrm{,}\;\;\;\boldsymbol{J}=\lambda
T^{2}\;\nabla\left(\frac{1}{T}\right)=\overline{\boldsymbol{Q}}\;,\;\;\;\mathrm{and}\;\;\;L=\lambda
T^{2},$ (8)

respectively. The noise correlation is then given by111Formally, the
temperature in the noise should be the deterministic value but in CFHD the
instantaneous temperature is typically used.

$\langle\widetilde{\boldsymbol{Q}}(\mathbf{r},t)\widetilde{\boldsymbol{Q}}(\mathbf{r}^{\prime},t^{\prime})\rangle=2k_{B}\lambda
T^{2}~{}\delta(\mathbf{r}-\mathbf{r}^{\prime})~{}\delta(t-t^{\prime}).$ (9)

Finally, by writing $e=c_{V}T$, where $c_{V}$ is the specific heat at constant
volume, the SHE is given as

$\rho c_{V}\frac{\partial T}{\partial t}=\nabla\cdot(\lambda\nabla
T+\sqrt{2k_{B}\lambda T^{2}}~{}\widetilde{\boldsymbol{Z}})$ (10)

where $\widetilde{\boldsymbol{Z}}$ is a Gaussian white noise vector field
uncorrelated in space and time, that is,

$\langle\widetilde{\boldsymbol{Z}}_{s}(\mathbf{r},t)\widetilde{\boldsymbol{Z}}_{s^{\prime}}(\mathbf{r}^{\prime},t^{\prime})\rangle=~{}\delta(\mathbf{r}-\mathbf{r}^{\prime})~{}\delta(t-t^{\prime})~{}\delta^{\mathrm{Kr}}_{s,s^{\prime}}$
(11)

where $s\in\\{x,y,z\\}$ is a Cartesian coordinate direction and
$\delta^{\mathrm{Kr}}$ is Kronecker delta.

For simplicity we focus on the one-dimensional case (e.g., thin rod) and write
the SHE in compact form as

$\displaystyle\partial_{t}T$
$\displaystyle=\partial_{x}\left(\kappa\partial_{x}T+\alpha T\tilde{Z}\right)$
(12) $\displaystyle=\kappa\partial_{x}^{2}\,T+\alpha\partial_{x}\,T\tilde{Z}$
(13)

where $\kappa=\lambda/\rho c_{V}$ and $\alpha=\sqrt{2k_{B}\lambda}/\rho c_{V}$
are taken to be constants. In the one-dimensional case the noise variance is

$\langle\tilde{Z}(x,t)\tilde{Z}(x^{\prime},t^{\prime})\rangle=\frac{1}{\mathcal{A}}~{}\delta(x-x^{\prime})~{}\delta(t-t^{\prime})$
(14)

where $\mathcal{A}$ is the cross-sectional area of the system.

### CFHD Schemes for SHE

Equation (13) can be numerically integrated in time to produce sample
trajectories for $T(x,t)$. We discretize space and time as $x_{i}=i\Delta x$
and $t^{n}=n\Delta t$ so temperature is $T_{i}^{n}=T(x_{i},t^{n})$. We use
forward difference in time

$\partial_{t}f\Rightarrow\frac{f_{i}^{n+1}-f_{i}^{n}}{\Delta t}$ (15)

and centered differences in space

$\partial_{x}f\Rightarrow\frac{f_{i+{\frac{1}{2}}}^{n}-f_{i-{\frac{1}{2}}}^{n}}{\Delta
x}\;,\qquad\partial_{x}^{2}f\Rightarrow\frac{f_{i+1}^{n}-2f_{i}^{n}+f_{i-1}^{n}}{\Delta
x^{2}}\;.$ (16)

Note that for the spatial index the integer values are at cell centers while
half-integer values are at cell faces.

With this we have the forward Euler (FE) scheme

$\displaystyle T_{i}^{n+1}=T_{i}^{n}$ $\displaystyle+\frac{\kappa\Delta
t}{\Delta x^{2}}\left(T_{i+1}^{n}-2T_{i}^{n}+T_{i-1}^{n}\right)$
$\displaystyle+\frac{\alpha\Delta t}{\Delta
x}\left(T_{i+{\frac{1}{2}}}^{n}\tilde{Z}_{i+{\frac{1}{2}}}^{n}-T_{i-{\frac{1}{2}}}^{n}\tilde{Z}_{i-{\frac{1}{2}}}^{n}\right).$
(17)

For the temperature at a face we can use the arithmetic average,
$T_{i+{\frac{1}{2}}}^{n}=(T_{i+1}^{n}+T_{i}^{n})/2$. Similarly, the standard
Predictor-Corrector (PC) scheme is

$\displaystyle T_{i}^{*}=T_{i}^{n}$ $\displaystyle+\frac{\kappa\Delta
t}{\Delta x^{2}}\left(T_{i+1}^{n}-2T_{i}^{n}+T_{i-1}^{n}\right)$
$\displaystyle+\frac{\alpha\Delta t}{\Delta
x}\left(T_{i+{\frac{1}{2}}}^{n}\tilde{Z}_{i+{\frac{1}{2}}}^{n}-T_{i-{\frac{1}{2}}}^{n}\tilde{Z}_{i-{\frac{1}{2}}}^{n}\right)$
(18)

and

$\displaystyle T_{i}^{n+1}=\frac{1}{2}\Big{[}T_{i}^{n}+T_{i}^{*}$
$\displaystyle+\frac{\kappa\Delta t}{\Delta
x^{2}}\left(T_{i+1}^{*}-2T_{i}^{*}+T_{i-1}^{*}\right)$
$\displaystyle+\frac{\alpha\Delta t}{\Delta
x}\left(T_{i+{\frac{1}{2}}}^{*}\tilde{Z}_{i+{\frac{1}{2}}}^{n}-T_{i-{\frac{1}{2}}}^{*}\tilde{Z}_{i-{\frac{1}{2}}}^{n}\right)\Big{]}$
(19)

where $T_{i}^{*}$, computed in the predictor step, is used in the corrector
step.222Random numbers used in predictor and corrector steps may be the same
or different; see [6].

The discretized noise has variance

$\langle\tilde{Z}_{i+{\frac{1}{2}}}^{n}\tilde{Z}_{j+{\frac{1}{2}}}^{m}\rangle=\frac{1}{\mathcal{A}}~{}\frac{\delta^{\mathrm{Kr}}_{i,j}}{\Delta
x}~{}\frac{\delta^{\mathrm{Kr}}_{n,m}}{\Delta
t}=\frac{\delta^{\mathrm{Kr}}_{i,j}}{\Delta
V}~{}\frac{\delta^{\mathrm{Kr}}_{n,m}}{\Delta t}$ (20)

where $\Delta V$ is the volume of a grid point cell. Numerically this noise is
generated as

$\tilde{Z}_{i+{\frac{1}{2}}}^{n}=\frac{1}{\sqrt{\Delta V\Delta
t}}~{}\mathfrak{N}_{i+{\frac{1}{2}}}^{n}$ (21)

where $\mathfrak{N}$ are independent, normal (Gaussian) distributed psuedo-
random numbers.

### CFHD example in Python

The Python program, StochasticHeat, illustrates the numerical methods
described in the previous section for the stochastic heat equation (SHE). The
program is listed in the appendix and the most up-to-date version can be
downloaded as a Jupyter Notebook from github.com/AlejGarcia/IntroFHD. The
simulation computes trajectories for the one-dimensional SHE and analyzes
their statistical properties (e.g., variance of temperature). The reader is
encouraged to download and run the program since results are obtained in about
5 minutes on a laptop.

At the top of the program there are various global options and parameters:

  * •

Periodic or Dirichlet boundary conditions

  * •

Thermodynamic equilibrium or constant temperature gradient

  * •

Initialize with or without temperature perturbations

  * •

Run simulation with or without thermal fluctuations

  * •

Use Forward Euler (FE) or Predictor-Corrector (PC) scheme

  * •

Number of grid cells, $N$

  * •

Number of time steps, $M$

Other parameters (e.g., system length, $\ell$) can also be changed from their
default values within the program. Simulation results in this section are from
runs using $N=32$ grid cells, running for either $M=2\times 10^{6}$
(equilibrium cases) or $M=2\times 10^{7}$ (temperature gradient case) time
steps. Default values were used for the other parameters (see program listing
in Appendix).

The program computes time-averaged quantities, for example the average
temperature in a cell is

$\langle
T_{i}\rangle=\frac{1}{M_{\mathrm{s}}}~{}\sum_{n=n_{0}+1}^{n_{0}+M_{\mathrm{s}}}T_{i}^{n}$
(22)

where $M_{\mathrm{s}}$ is the number of statistical samples; the system
“relaxes” for the first $n_{0}$ time steps before sampling begins. At
thermodynamic equilibrium the variance and static correlation of temperature
fluctuations ($\delta T_{i}^{n}=T_{i}^{n}-\langle T_{i}\rangle$) are predicted
by statistical mechanics to be [7]

$\langle\delta T_{i}^{2}\rangle=\frac{k_{B}\langle T_{i}\rangle^{2}}{\rho
c_{V}\Delta V}\qquad\mathrm{and}\qquad\langle\delta T_{i}\delta
T_{j}\rangle=\langle\delta T_{i}^{2}\rangle~{}\delta^{\mathrm{Kr}}_{i,j}$ (23)

with the average temperature, $\langle T_{i}\rangle$, equal to the equilibrium
temperature, $T_{\mathrm{eq}}$. Figures 1 and 2 show that the simulation
results for Dirichlet boundaries (i.e., fixed temperature at the boundaries)
are in good agreement with theory, especially for the PC scheme.

Figure 1: Temperature variance $\langle\delta T_{i}^{2}\rangle$ versus $x_{i}$
at equilibrium with Dirichlet boundary conditions for (left) FE scheme;
(right) PC scheme. Theory line given by (23).

Figure 2: Temperature correlation $\langle\delta T_{i}\delta T_{j}\rangle$
versus $x_{i}$ at equilibrium with Dirichlet boundary conditions with
$x_{j}=\ell/4$ for (left) FE scheme; (right) PC scheme. Theory line given by
(23).

Figure 3: Temperature correlation $\langle\delta T_{i}\delta T_{j}\rangle$
versus $x_{i}$ for periodic boundary conditions using the PC scheme. Results
are plotted both with (left) and without (right) the point at $x_{j}=\ell/4$.

Figure 3 shows that temperature correlations are slightly different for
periodic boundary conditions. This is because, due to conservation of energy,
$\sum_{i}\rho c_{V}\langle\delta T_{i}\rangle=0$, the correlation is shifted
by an additive constant such that $\sum_{i}\langle\delta T_{i}\delta
T_{j}\rangle=0$. This is a common feature in systems with conserved variables,
such as density fluctuations in closed systems.

The assessment of numerical schemes in CFHD is best done by measuring the
static structure factor for fluctuations. For the SHE, we define the
unnormalized structure factor,
$S_{k}=\langle\hat{T}_{k}\hat{T}^{*}_{k}\rangle$, where

$\hat{T}_{k}=\sum_{j=0}^{N-1}T_{j}\exp(-2\pi\mathrm{i}\,jk/N)$ (24)

is the discrete Fourier transform with wavenumber index $k=0,1,\ldots,N-1$,
and $\hat{T}^{*}_{k}$ is its complex conjugate transpose.333In CFHD the
structure factor is often normalized to unity using the equilibrium variance.
At thermodynamic equilibrium the discretized structure factor is

$S_{k}=\frac{k_{B}T_{\mathrm{eq}}^{2}}{\rho c_{V}}~{}N.$ (25)

Figure 4 shows that the FE scheme has significant error for large $k$ while
the PC scheme is accurate for all $k$. An analysis of discretization errors
for various numerical schemes is presented in [6]; for FE the error in $S_{k}$
is $O(\kappa\Delta t\,k^{2})$ while for PC it is $O(\kappa^{2}\Delta
t^{2}k^{4})$.

Figure 4: Static structure factor, $S_{k}$, versus $k$ at equilibrium for
(left) FE scheme; (right) PC scheme. Theory line given by (25).

Finally we consider a non-equilibrium scenario, specifically a linear
temperature gradient imposed by Dirichlet boundary conditions. In this case
there is a weak, long-range correlation of temperature fluctuations [8]:

$\langle\delta T_{i}\delta T_{j}\rangle=\frac{k_{B}T_{\mathrm{eq}}^{2}}{\rho
c_{V}\Delta V}\delta^{\mathrm{Kr}}_{i,j}+\frac{k_{B}(\nabla T)^{2}}{\rho
c_{V}\mathcal{A}\ell}\times\left\\{\begin{array}[]{cc}x_{i}(\ell-
x_{j})&\mathrm{if}~{}x_{i}<x_{j},\\\ x_{j}(\ell-
x_{i})&\mathrm{otherwise}.\end{array}\right.$ (26)

Figure 5 shows that the simulation results are in good agreement with this
theoretical result.

Figure 5: Temperature correlation $\langle\delta T_{i}\delta T_{j}\rangle$
versus $x_{i}$ for linear temperature gradient imposed by Dirichlet boundary
conditions using the PC scheme. Results are plotted with (left) and without
(right) the point at $x_{j}=\ell/4$; theory line given by (26).

### Stochastic Species Diffusion

In fluctuating hydrodynamics there are several stochastic diffusion equations
that are closely related to the stochastic heat equation. For example, the
diffusion of species mass is described by a similar stochastic PDE with the
main difference being the functional form of the stochastic flux. As above we
start with

$\partial_{t}(\rho c)=-\nabla\cdot\boldsymbol{F}$ (27)

where $c$ is the species mass fraction and the species flux is

$\boldsymbol{F}=\overline{\boldsymbol{F}}+\widetilde{\boldsymbol{F}}=\mathcal{L}\boldsymbol{\mathcal{X}}+\widetilde{\boldsymbol{F}}.$
(28)

Assuming that the system is isothermal (i.e., no Soret effect), irreversible
thermodynamics tells us that the thermodynamic force is [4, 5]

$\boldsymbol{\mathcal{X}}=-\nabla\left(\frac{\mu}{T}\right)=-\frac{\nabla\mu}{T}$
(29)

where $\mu$ is the chemical potential. For ideal solutions

$\mu=\mu_{0}(T)+\frac{k_{B}T}{m}\ln
c\qquad\mathrm{so}\qquad\nabla\mu=\frac{k_{B}T}{mc}~{}\nabla c$ (30)

where $m$ is the particle mass. Since the phenomenological law for species
flow is Fick’s law, $\overline{\boldsymbol{F}}=-\rho D\nabla c$, then the
Onsager coefficient is

$\mathcal{L}=\frac{m\rho D\,c}{k_{B}}\;.$ (31)

The noise correlation is

$\langle\widetilde{\boldsymbol{F}}(\mathbf{r},t)\widetilde{\boldsymbol{F}}(\mathbf{r}^{\prime},t^{\prime})\rangle=2m\rho
D\,c~{}\delta(\mathbf{r}-\mathbf{r}^{\prime})~{}\delta(t-t^{\prime})$ (32)

so the stochastic species diffusion equation is

$\partial_{t}\,c=\nabla\cdot\left(D\nabla
c+\sqrt{\frac{2D\,c}{n_{0}}}~{}\widetilde{\boldsymbol{Z}}\right)$ (33)

where $n_{0}=\rho/m$ is the number density of the pure state ($c=1$). This may
also be written in terms of species number density, $n=c\,n_{0}$, as

$\partial_{t}\,n=\nabla\cdot(D\nabla
n+\sqrt{2D\,n}~{}\widetilde{\boldsymbol{Z}}).$ (34)

This stochastic PDE is also known as the Dean-Kawasaki equation. Recall that
the stochastic heat equation is

$\partial_{t}\,T=\nabla\cdot\left(\kappa\nabla T+\sqrt{\frac{2\kappa
k_{B}\,T^{2}}{\rho c_{V}}}~{}\widetilde{\boldsymbol{Z}}\right).$ (35)

Note that the deterministic parts of the two equations above are functionally
similar while the stochastic parts are distinctly different.

The stochastic species diffusion equation can also be derived by coarse-
graining the discrete random walk model for diffusion. Additional examples of
discrete models that have associated stochastic diffusion equations are: the
excluded random walk model; the “train model” for the diffusion of momentum
[9, 10]; and the “Knudsen chain” model for gas effusion [11].

### A Menagerie of Fluctuating Hydrodynamic Equations

The general form of the fluctuating hydrodynamic equations for a set of state
variables $\boldsymbol{U}$ is

$\partial_{t}\boldsymbol{U}=-\nabla\cdot(\boldsymbol{\mathcal{F}}_{H}(\boldsymbol{U})+\boldsymbol{\mathcal{F}}_{D}(\boldsymbol{U})+\boldsymbol{\mathcal{N}}(\boldsymbol{U})\widetilde{\boldsymbol{Z}})$
(36)

where $\boldsymbol{\mathcal{F}}_{H}$ and $\boldsymbol{\mathcal{F}}_{D}$ are
the hyperbolic and diffusive fluxes with $\boldsymbol{\mathcal{N}}$ being the
noise amplitude. For example, for the stochastic species diffusion equation
$\boldsymbol{U}=n$, $\boldsymbol{\mathcal{F}}_{H}=0$,
$\boldsymbol{\mathcal{F}}_{D}=-D\nabla n$, and
$\boldsymbol{\mathcal{N}}=\sqrt{2Dn}$. In general, the fluctuation-dissipation
theorem relates $\boldsymbol{\mathcal{F}}_{D}$ to $\boldsymbol{\mathcal{N}}$
and tells us that the stochastic term is independent of the hyperbolic term.

This brief introduction to CFHD has focused on stochastic diffusion equations
due to the pedagogical benefit of their simplicity. We close with an outline
of other fluctuating hydrodynamic equations including references describing
finite volume schemes for solving them.

  * •

Stochastic Burger’s equation [12]:

$\partial_{t}\,c=-\partial_{x}\,(ac(1-c)-D\,\partial_{x}c+\sqrt{2Dc(1-c)}\,\tilde{Z}).$
(37)

Note that for $c\ll 1$ this resembles the stochastic species diffusion
equation with the addition of a hyperbolic flux, $\boldsymbol{F}_{H}=ac(1-c)$.

  * •

Stochastic “train model” equations [9, 10]:

$\displaystyle\partial_{t}\,\rho$
$\displaystyle=-\partial_{x}(-D\partial_{x}\rho+\sqrt{2mD\rho}\,\tilde{Z}),$
(38) $\displaystyle\partial_{t}\,(\rho u)$
$\displaystyle=-\partial_{x}(-D\partial_{x}(\rho u)+\sqrt{2mD\rho
u^{2}}\,\tilde{Z}).$ (39)

These equations are a simple model for the transport of momentum density,
$\rho u$, in a fluid.

  * •

Reaction-diffusion systems [13]:
The reaction-diffusion system has $N_{s}$ species diffusing and undergoing
$N_{r}$ reactions. By denoting the number density of species $s$ by $n_{s}$,
the equations of fluctuating reaction-diffusion can be written formally as the
stochastic PDEs

$\frac{\partial n_{s}}{\partial t}=\nabla\cdot\left(D_{s}\nabla
n_{s}+\sqrt{2D_{s}n_{s}}\boldsymbol{\mathcal{Z}}_{s}^{(D)}\right)+\sum_{r=1}^{N_{r}}\nu_{sr}\left(a_{r}({\bf
n})+\sqrt{a_{r}({\bf n})}\mathcal{Z}_{r}^{(R)}\right)$ (40)

where $D_{s}$ is the diffusion coefficient of species $s$, $a_{r}({\bf n})$ is
the propensity function indicating the rate of reaction $r$, and $\nu_{sr}$ is
the stoichiometric coefficient of species $s$ in reaction $r$.

  * •

Compressible Navier-Stokes equations (single species fluid) [14, 15]:

$\displaystyle\partial_{t}\,\rho=-\nabla\cdot\left(\rho\boldsymbol{u}\right),$
(41)
$\displaystyle\partial_{t}\,(\rho\boldsymbol{u})=-\nabla\cdot\left[(\rho\boldsymbol{u}\otimes\boldsymbol{u}+p\mathbb{I})+\overline{\Pi}+\widetilde{\Pi}\right],$
(42) $\displaystyle\partial_{t}\,(\rho E)=-\nabla\cdot\left[\left(\rho
E+p\right)\boldsymbol{u}+\overline{\boldsymbol{Q}}+\widetilde{\boldsymbol{Q}}+\left(\overline{\Pi}+\widetilde{\Pi}\right)\cdot\boldsymbol{u}\right],$
(43)

where $\boldsymbol{u}$, $p$, and $\Pi=\overline{\Pi}+\widetilde{\Pi}$ are
fluid velocity, pressure, and stress tensor; total specific energy is
$E=e+\frac{1}{2}u^{2}$. Note that (43) reduces to the stochastic heat equation
when $\boldsymbol{u}=0$. See references for explicit expression of the
deterministic and stochastic fluxes.

  * •

Compressible Navier-Stokes equations (multi-species fluid) [15, 16]:
Same as above but replace (41) with

$\partial_{t}\,\rho_{k}=-\nabla\cdot(\rho_{k}\boldsymbol{u}+\overline{\boldsymbol{F}}_{k}+\widetilde{\boldsymbol{F}}_{k})$
(44)

where $\rho_{k}$ is the mass density for species $k$ and
$\rho=\sum_{k}\rho_{k}$. For ideal mixtures the diffusive flux,
$\boldsymbol{F}_{k}$, is similar to that of the species diffusion equation
described in the previous section.

  * •

Hydrodynamics plus chemical reactions [17, 18]:
Add source terms to the right-hand side of (44) for the deterministic and
stochastic rates of reactions, $\overline{\Omega}_{k}$ and
$\widetilde{\Omega}_{k}$. These can be formulated from the chemical Langevin
equation [19].

  * •

Incompressible / low Mach hydrodynamic equations [20, 21, 22, 23]:
The multi-species methodology can be extended to model isothermal mixtures of
miscible incompressible liquids. Incompressibility of the fluids leads to a
constrained evolution, similar to the incompressible Navier-Stokes equations,
given by

$\displaystyle\partial_{t}\,(\rho\boldsymbol{u})=-\nabla\cdot\left[\rho\boldsymbol{u}\otimes\boldsymbol{u}+\overline{\Pi}+\widetilde{\Pi}\right]-\nabla\pi,$
(45)
$\displaystyle\partial_{t}\,\rho_{k}=-\nabla\cdot(\rho_{k}\boldsymbol{u}+\overline{\boldsymbol{F}}_{k}+\widetilde{\boldsymbol{F}}_{k}),$
(46)
$\displaystyle\nabla\cdot\boldsymbol{u}=-\nabla\cdot\left(\sum_{k}\frac{\overline{\boldsymbol{F}}_{k}+\widetilde{\boldsymbol{F}}_{k}}{\bar{\rho}_{k}}\right).$
(47)

Here, $\pi$ is a perturbational pressure and $\bar{\rho}_{k}$ is the pure
component density for species $k$. The resulting system retains in influence
of fluctuations on mixing but is considerable more computationally efficient,
particularly for liquids with a high sound speed.

  * •

Stochastic Poisson-Nernst-Planck equations [24, 25]:
The incompressible version of CFHD for multi-species fluids can model
electrolyte solutions by including charged species (ions). The chemical
potential becomes the electrochemical potential with the electrical mobility
given by the Nernst–Einstein relation and the diffusive flux by the Nernst-
Planck equation. For scales comparable to the Debye length, the electric field
is obtained by solving the Poisson equation. At scales larger than the Debye
length, the fluid is electro-neutral. Electro-neutrality can be imposed as a
constraint by solving a modified elliptic equation to compute the electric
potential.

  * •

Multi-phase fluids [26, 27, 28, 29]:
By adding non-local gradient terms to the free energy, one can incorporate
interfacial tension in CFHD model to treat multiphase systems. This can be
done for both single component systems (e.g., liquid/vapor) and for
multicomponent systems (e.g., oil/water). For multicomponent systems, the low
Mach number version of the methodology has also be derived. These types of
systems introduce additional higher order operators into the momentum
equations and, for multicomponent systems, the species transport equations.
See the references for specific forms of these equations.

Finally, github.com/AMReX-FHD is a repository of CFHD codes written using the
AMReX framework.

## References

  * [1] L. D. Landau and E. M. Lifshitz.  Fluid Mechanics, Course of Theoretical Physics, Vol. 6.  Pergamon Press, 1959.
  * [2] J. M. Ortiz de Zarate and J. V. Sengers.  Hydrodynamic Fluctuations in Fluids and Fluid Mixtures.  Elsevier Science, 2007.
  * [3] Crispin W. Gardiner.  Handbook of Stochastic Methods.  Springer, Berlin, 1985.
  * [4] S. R. DeGroot and P. Mazur.  Non-Equilibrium Thermodynamics.  North-Holland Publishing Company, Amsterdam, 1963.
  * [5] A. L. Garcia.  Essentials of Modern Thermodynamics.  Amazon, 2022.
  * [6] Aleksandar Donev, Eric Vanden-Eijnden, Alejandro Garcia, and John Bell.  On the accuracy of finite-volume schemes for fluctuating hydrodynamics.  Communications in Applied Mathematics and Computational Science, 5(2):149–197, 2010.
  * [7] Raj Kumar Pathria.  Statistical Mechanics.  Elsevier, 2016.
  * [8] Alejandro L. Garcia, M. Malek Mansour, George C. Lie, and Enrico Cementi.  Numerical integration of the fluctuating hydrodynamic equations.  Journal of Statistical Physics, 47:209–228, 1987.
  * [9] Alejandro Garcia, F. Baras, and M. Malek Mansour.  A simple model for nonequilibrium fluctuations in a fluid.  American Journal of Physics, 64:1488, 1996.
  * [10] Francis J. Alexander, Alejandro L. Garcia, and Daniel M. Tartakovsky.  Algorithm refinement for stochastic partial differential equations: II. Correlated systems.  Journal of Computational Physics, 207(2):769–787, 2005.
  * [11] Alejandro L. Garcia.  Thermal fluctuations in a Knudsen flow system.  Physics Letters A, 119(8):379–382, 1987.
  * [12] John B. Bell, Jasmine Foo, and Alejandro L. Garcia.  Algorithm refinement for the stochastic Burgers’ equation.  Journal of Computational Physics, 223(1):451–468, 2007.
  * [13] Changho Kim, Andy Nonaka, John B. Bell, Alejandro L. Garcia, and Aleksandar Donev.  Stochastic simulation of reaction-diffusion systems: A fluctuating-hydrodynamics approach.  The Journal of Chemical Physics, 146:124110, 2017.
  * [14] John B. Bell, Alejandro L. Garcia, and Sarah A. Williams.  Numerical methods for the stochastic Landau-Lifshitz Navier-Stokes equations.  Physical Review E, 76:016708, 2007.
  * [15] Ishan Srivastava, Daniel R. Ladiges, Andy J. Nonaka, Alejandro L. Garcia, and John B. Bell.  Staggered scheme for the compressible fluctuating hydrodynamics of multispecies fluid mixtures.  Physical Review E, 107:015305, 2023.
  * [16] Kaushik Balakrishnan, Alejandro L. Garcia, Aleksandar Donev, and John B. Bell.  Fluctuating hydrodynamics of multispecies nonreactive mixtures.  Physical Review E, 89:013017, 2014.
  * [17] Amit Kumar Bhattacharjee, Kaushik Balakrishnan, Alejandro L. Garcia, John B. Bell, and Aleksandar Donev.  Fluctuating hydrodynamics of multi-species reactive mixtures.  The Journal of Chemical Physics, 142:224107, 2015.
  * [18] Changho Kim, Andy Nonaka, John B. Bell, Alejandro L. Garcia, and Aleksandar Donev.  Fluctuating hydrodynamics of reactive liquid mixtures.  The Journal of Chemical Physics, 149:084113, 2018.
  * [19] Daniel T. Gillespie.  The chemical Langevin equation.  The Journal of Chemical Physics, 113(1):297–306, 2000.
  * [20] Aleksandar Donev, Andy Nonaka, Yifei Sun, Thomas Fai, Alejandro Garcia, and John Bell.  Low Mach number fluctuating hydrodynamics of diffusively mixing fluids.  Communications in Applied Mathematics and Computational Science, 9(1):47–105, 2014.
  * [21] Aleksandar Donev, Andy Nonaka, Amit Kumar Bhattacharjee, Alejandro L. Garcia, and John B. Bell.  Low Mach number fluctuating hydrodynamics of multispecies liquid mixtures.  Physics of Fluids, 27(3), 2015.
  * [22] Andrew Nonaka, Yifei Sun, John Bell, and Aleksandar Donev.  Low Mach number fluctuating hydrodynamics of binary liquid mixtures.  Communications in Applied Mathematics and Computational Science, 10(2):163–204, 2015.
  * [23] Florencio Balboa, John B. Bell, Rafael Delgado-Buscalioni, Aleksandar Donev, Thomas G. Fai, Boyce E. Griffith, and Charles S. Peskin.  Staggered schemes for fluctuating hydrodynamics.  Multiscale Modeling & Simulation, 10(4):1369–1408, 2012.
  * [24] Jean-Philippe Péraud, Andy Nonaka, Anuj Chaudhri, John B. Bell, Aleksandar Donev, and Alejandro L. Garcia.  Low Mach number fluctuating hydrodynamics for electrolytes.  Physical Review Fluids, 1:074103, 2016.
  * [25] Aleksandar Donev, Andrew J. Nonaka, Changho Kim, Alejandro L. Garcia, and John B. Bell.  Fluctuating hydrodynamics of electrolytes at electroneutral scales.  Physical Review Fluids, 4:043701, 2019.
  * [26] Barry Z. Shang, Nikolaos K. Voulgarakis, and Jhih-Wei Chu.  Fluctuating hydrodynamics for multiscale simulation of inhomogeneous fluids: Mapping all-atom molecular dynamics to capillary waves.  The Journal of Chemical Physics, 135:044111, 2011.
  * [27] Anuj Chaudhri, John B. Bell, Alejandro L. Garcia, and Aleksandar Donev.  Modeling multiphase flow using fluctuating hydrodynamics.  Physical Review E, 90:033014, 2014.
  * [28] Katherine Klymko, Andrew Nonaka, John B. Bell, Sean P. Carney, and Alejandro L. Garcia.  Low Mach number fluctuating hydrodynamics model for ionic liquids.  Physical Review Fluids, 5:093701, 2020.
  * [29] Bryn Barker, John B. Bell, and Alejandro L. Garcia.  Fluctuating hydrodynamics and the Rayleigh-Plateau instability.  Proceedings of the National Academy of Sciences, 120(30):e2306088120, 2023.

## APPENDIX: Python program StochHeat

    # StochasticHeat - Program to calculate trajectories for the stochastic heat equation

    # Set up configuration options and special features
    import numpy as np                 # Use numerical python library
    import matplotlib.pyplot as plt    # Use plotting python library
    # Display plots in notebook window
    get_ipython().run_line_magic(’matplotlib’, ’inline’)

    # In[ ]:

    #* Set option flags and key parameters for the simulation
    BC_FLAG = 0                   # Boundary condition flag (0: Periodic; 1: Dirichlet)
    NONEQ_FLAG = 0                # Non-equilibrium flag (0: Thermodynamic equilibrium; 1: Temperature gradient)
    if( BC_FLAG == 0 and NONEQ_FLAG == 1 ):
        print(’ERROR: Periodic boundary conditions are only for thermodynamic equilibrium’)
        NONEQ_FLAG = 0            # Reset the flag to equilibrium
    PERTURB_FLAG = 1              # Start from perturbed initial condition (0: No; 1: Yes)
    STOCH_FLAG = 1                # Thermal fluctuations? (0: No, deterministic; 1: Yes, stochastic)
    SCHEME_FLAG = 1               # Numerical scheme (0: Forward Euler; 1: Predictor-corrector)

    Ncells = 32                   # Number of cells (including boundary cells ’0’ and ’Ncells-1’ )
    MegaSteps = 2                 # Number of simulation timesteps in millions
    # Recommend at least 2 million steps for equilibrium systems, 20 million for non-equilibrium systems

    # In[ ]:

    #* Set physical parameters for the system (iron bar)
    kB = 1.38e-23       # Boltzmann constant (J/K)
    mAtom = 9.27e-26    # Mass of iron atom (kg)
    rho = 7870.         # Mass density of iron (kg/m^3)
    c_V = 450.          # Specific heat capacity of iron (J/(kg K))
    ThCond = 70.        # Thermal conductivity of iron (W/(m K))
    Length = 1.0e-8     # System length (m)
    Area = (5.0e-9)**2  # System cross-sectional area (m^2)

    # In[ ]:

    #* Set numerical parameters (time step, grid spacing, etc.).
    if( BC_FLAG == 0 ):
        cellLo = 0;  cellHi = Ncells-1;     # For periodic BC, range of cells = { 0, 1, ..., Ncells-2 }
        faceLo = 0;  faceHi = Ncells-1;     # Face ii is left side of cell ii
    else:
        cellLo = 1;  cellHi = Ncells-1;     # For Dirichlet BC, range of cells = { 1, 2, ..., Ncells-2 }
        faceLo = 1;  faceHi = Ncells-1;     # Face ii is left side of cell ii

    dx = Length/(Ncells-1)     # Grid size (m)
    dV = Area*dx               # Volume of grid cell (m^3)
    xx = np.zeros(Ncells)
    for ii in range(Ncells):
        xx[ii] = ii * dx       # Cell center positions (m)

    kappa = ThCond / (rho * c_V )                # Coefficient in deterministic heat equation
    # Coefficient in stochastic heat equation
    if( STOCH_FLAG == 1 ):
        alpha = np.sqrt( 2 * kB * kappa / (rho * c_V) )
    else:
        alpha = 0.  # Set alpha = 0 to turn off thermal fluctuations

    stabilityFactor = 0.1                         # Numerical stability if stabilityFactor < 1.
    dt = stabilityFactor * dx**2 /( 2. * kappa )  # Timestep (s)

    # In[ ]:

    #* Set initial conditions for temperature
    Tref = 300.      # Reference temperature (K)
    if( NONEQ_FLAG == 0 ):
        T_Left = Tref; T_Right = Tref
    else:
        Tdiff = 400.   # Temperature difference across the system
        T_Left = Tref - Tdiff/2.; T_Right = Tref + Tdiff/2.

    # Standard deviation of temperature in a cell at the reference temperature
    Tref_SD = np.sqrt(kB*Tref**2 / (rho * c_V * dV))

    # Set initial temperature and its deterministic steady-state value
    T = np.zeros(Ncells); T0 = np.zeros(Ncells)
    for ii in range(cellLo,cellHi):
        T0[ii] = T_Left + (T_Right-T_Left) * ii/(Ncells-1)      # Linear profile
        T[ii] = T0[ii]
        if( PERTURB_FLAG == 1 ):
            T[ii] += Tref_SD*np.random.randn()  # Add random perturbation

    if( BC_FLAG == 0 ):
        T0[cellHi] = T0[0]       # Copy first cell into last cell for periodic BCs
        T[cellHi] = T[0]
    else:
        T0[0] = T_Left;   T[0] = T_Left      # Set values in first and last cells
        T0[-1] = T_Right; T[-1] = T_Right

    # In[ ]:

    #* Summarize system parameters and initial state
    print( ’System is iron bar with about ’, int(rho*Length*Area/mAtom), ’ atoms’ )
    print( ’System length = ’,Length*1.0e9,’ (nm); volume = ’,
          (Length*Area)*1.0e27,’ (nm**3)’ )
    print( ’Number of timesteps = ’, MegaSteps,’ million’ )
    print( ’Time step = ’, dt*1.e15, ’ (fs)’ )
    print( ’Number of cells = ’, Ncells )
    if( BC_FLAG == 0 ):
        print(’** PERIODIC boundary conditions **’)
    else:
        print(’++ DIRICHLET boundary conditions ++’)

    if( NONEQ_FLAG == 0 ):
        print(’Equilibrium temperature = ’,Tref,’ Kelvin’)
    else:
        print(’Fixed temperatures (K) = ’,T_Left,’ left; ’,T_Right,’ right’)
    print( ’At T = ’,Tref,’K, standard deviation sqrt(< dT^2 >) = ’, Tref_SD,’K’ )
    if( SCHEME_FLAG == 0 ):
        print(’-- Forward Euler Scheme --’)
    else:
        print(’== Predictor-Corrector Scheme ==’)

    # Plot temperature versus position
    plt.plot(xx, T,’g*’,xx, T0,’b--’,xx[0],T[0],’sk’,xx[-1],T[-1],’sk’)
    plt.xlabel(’Position (m)’); plt.ylabel(’Temperature (K)’)
    plt.title(’Initial state (stars), steady state (dashed), boundary cells (squares)’)
    plt.show()

    # In[ ]:

    #* Initialize sampling plus misc. coefficients and arrays
    Nsamp = 0                   # Count number of statistical samples
    sumT  = np.zeros(Ncells)    # Running sum of T_i ; used to compute mean
    sumT2 = np.zeros(Ncells)    # Running sum of (T_i)**2 ; used to compute variance
    sumTT = np.zeros(Ncells)    # Running sum for correlation T_i * T_iCorr
    sumSk = np.zeros(Ncells)    # Running sum for structure factor S(k)

    iCorr = np.int_(Ncells/4)   # Grid point used for correlation

    # Coefficients used in the main loop calculations
    coeffDetFE = kappa * dt / dx**2
    coeffStoFE = alpha * dt / dx
    coeffZnoise = 1. / np.sqrt( dt * dV )

    # Arrays used in the main loop calculations
    Determ = np.zeros(Ncells); Stoch = np.zeros(Ncells)
    Znoise = np.zeros(Ncells); Tface = np.zeros(Ncells)
    PreT = np.copy(T)             # Used by Predictor-Corrector scheme
    Spectrum = np.zeros(Ncells)   # Used to compute structure factor S(k)

    # In[ ]:

    #* Loop over the desired number of time steps.
    NstepInner = 10                                       # Number of time steps (inner loop)
    NstepOuter = np.int_(MegaSteps*1000000/NstepInner)    # Number of time steps (outer loop)
    NskipOuter = NstepOuter/10        # Number of outer steps to skip before sampling begins
    NdiagOuter = NstepOuter/20        # Number of outer steps between diagnostic outputs

    for iOuter in range(NstepOuter):     # Outer loop

        # Print diagnostics
        if (iOuter % NdiagOuter) == 0 :
            print( ’Finished ’,np.int_(100*iOuter/NstepOuter),’ percent of the time steps’)

        for iInner in range(NstepInner):       # Inner loop

            # Deterministic update for temperature
            for ii in range(cellLo,cellHi):
                Determ[ii] = coeffDetFE * ( T[ ii+1 ] + T[ ii-1 ] - 2*T[ ii ] )
            if( BC_FLAG == 0 ):
                Determ[0] = coeffDetFE * ( T[ 1 ] + T[ cellHi-1 ] - 2*T[ 0 ] )  # Periodic BC

            # Generate random noise Z
            Znoise = coeffZnoise * np.random.normal(size=Ncells)
            if( BC_FLAG == 0 ):
                Znoise[-1] = Znoise[0]   # Periodic BC

            # Tface[ i ] is average between T of cells i-1 and i; value on the left face of cell i
            for ii in range(cellLo,cellHi):
                Tface[ii] = 0.5 * (T[ ii-1 ] + T[ ii ])
            if( BC_FLAG == 0 ):
                Tface[cellHi] = Tface[0]   # Periodic BC
            else:
                Tface[cellHi] = 0.5 * (T[ cellHi-1 ] + T[ cellHi ])   # Dirichlet BC

            # Stochastic update for temperature
            for ii in range(cellLo,cellHi):
                Stoch[ii] = coeffStoFE * ( Tface[ii+1]*Znoise[ii+1] - Tface[ii]*Znoise[ii] )

            if( SCHEME_FLAG == 0 ):
                # Forward Euler scheme
                for ii in range(cellLo,cellHi):
                    T[ii] += Determ[ii] + Stoch[ii]   # Total update for temperature
            else:
                # Predictor-Corrector scheme
                for ii in range(cellLo,cellHi):
                    PreT[ii] = T[ii] + Determ[ii] + Stoch[ii]     # Predictor step
                if( BC_FLAG == 0 ):
                    PreT[cellHi] = PreT[0]     # Periodic BC
                # Corrector step
                for ii in range(cellLo,cellHi):
                    Determ[ii] = coeffDetFE * ( PreT[ ii+1 ] + PreT[ ii-1 ] - 2*PreT[ ii ] )
                if( BC_FLAG == 0 ):
                    Determ[0] = coeffDetFE * ( PreT[ 1 ] + PreT[ cellHi-1 ] - 2*PreT[ 0 ] )  # Periodic BC
                for ii in range(cellLo,cellHi):
                    Tface[ii] = 0.5 * (PreT[ ii-1 ] + PreT[ ii ])
                if( BC_FLAG == 0 ):
                    Tface[cellHi] = Tface[0]   # Periodic BC
                else:
                    Tface[cellHi] = 0.5 * (PreT[ cellHi-1 ] + PreT[ cellHi ])
                for ii in range(cellLo,cellHi):
                    Stoch[ii] = coeffStoFE * ( Tface[ii+1]*Znoise[ii+1] - Tface[ii]*Znoise[ii] )
                for ii in range(cellLo,cellHi):
                    T[ii] = 0.5*( T[ii] + PreT[ii] + Determ[ii] + Stoch[ii] )

            if( BC_FLAG == 0 ):
                T[cellHi] = T[0]     # Periodic BC

        # End of Inner loop

      # Take statistical sample
        if( iOuter > NskipOuter ):
            Nsamp += 1
            for ii in range(cellLo,cellHi):
                sumT[ii] += T[ii]            # Running sum for temperature average
                sumT2[ii] += T[ii]**2        # Running sum for temperature variance
                sumTT[ii] += T[ii]*T[iCorr]  # Running sum for temperature correlation

            if( NONEQ_FLAG == 0 ):
                # Take Fourier transform and record sampled spectrum in running sum
                Spectrum[cellLo:cellHi] = np.abs( np.fft.fft( T[cellLo:cellHi] ) )**2
                for ii in range(cellLo,cellHi):
                    sumSk[ii] += Spectrum[ii]

    # End of Outer loop

    # In[ ]:

    #* Calculate average, variance, and correlation
    aveT = np.zeros(Ncells)     # Average < T_i >
    varT = np.zeros(Ncells)     # Variance < (T_i - <T_i>)**2 >
    corrT = np.zeros(Ncells)    # Correlation < T_i T_iCorr >

    for ii in range(cellLo,cellHi):
        aveT[ii] = sumT[ii] / Nsamp
        varT[ii] = sumT2[ii]/Nsamp - aveT[ii]**2
    for ii in range(cellLo,cellHi):
        corrT[ii] = sumTT[ii]/Nsamp - aveT[ii]*aveT[iCorr]

    # In[ ]:

    #* Plot temperature versus x
    plt.plot(xx[cellLo:cellHi], T[cellLo:cellHi],’g*’,
             xx[cellLo:cellHi],aveT[cellLo:cellHi],’ro’,xx[cellLo:cellHi],T0[cellLo:cellHi],’b--’)
    plt.xlabel(’Position (m)’); plt.ylabel(’Temperature (K)’)
    plt.title(’Final state (stars), Average (circles), Deterministic (line)’)
    plt.show()

    # In[ ]:

    #* Plot varicance of temperature versus x; compare with theory
    varT_Th = np.zeros(Ncells)   # Theoretical value
    varT_Th[cellLo:cellHi] = kB*aveT[cellLo:cellHi]**2 / (rho * c_V * dV)
    # Conservation of energy correction for periodic boundary case
    if( BC_FLAG == 0 ):
        varT_Th *= (1 - 1/(Ncells-1))

    plt.plot(xx[cellLo:cellHi], varT[cellLo:cellHi],’ro’,
             xx[cellLo:cellHi],varT_Th[cellLo:cellHi],’b--’)
    plt.xlabel(’Position (m)’); plt.ylabel(’Variance < dT**2 >’)
    plt.title(’Average (circles), Equilibrium theory (dashed)’)
    plt.show()

    # In[ ]:

    #* Plot temperature cell correlation versus x
    varTheory = kB*aveT[iCorr]**2 / (rho * c_V * dV)  # Temperature variance in cell iCorr
    if( BC_FLAG == 0 ):
        corrT_Th = - varTheory * 1/(Ncells-1) * np.ones(Ncells)  # Conservation of energy correction
        corrT_Th[iCorr] = varTheory * (1 - 1/(Ncells-1))         # for periodic boundary case
    else:
        corrT_Th = np.zeros(Ncells)                              # Kronecker delta correlation
        corrT_Th[iCorr] = varTheory                              # for Dirichlet boundary case

    plt.plot(xx[cellLo:cellHi], corrT[cellLo:cellHi],’ro’,
             xx[cellLo:cellHi],corrT_Th[cellLo:cellHi],’b--’)
    plt.xlabel(’Position (m)’); plt.ylabel(’Correlation < dT(x) dT(x*) >’)
    plt.title(’Average (circles), Theory (dashed)’)
    plt.show()
    print(’Correlation for x* = ’, xx[iCorr],’ (m)’)

    # In[ ]:

    #* Plot temperature cell correlation versus x WITHOUT the equilibrium delta function
    corrTd = np.array(corrT)
    corrTd[iCorr] = np.nan        # Delete data point at iCorr using NaN

    # Calculate theoretical values for < dT(x) dT(x’) > where x’ = xx[iCorr]
    if( NONEQ_FLAG == 0 ):
        corrTd_Th = np.array(corrT_Th)
        corrTd_Th[iCorr] = np.nan        # Delete data point at iCorr using NaN
    else:
        corrTd_Th = np.zeros(Ncells)
        for ii in range(Ncells):
            # Non-equilibrium correlation; see J. Stat. Phys. 47 209 (1987)
            corrTd_Th[ii] = kB*(T_Right-T_Left)**2 / (rho * c_V * Area * Length**3)
            if( ii < iCorr ):
                corrTd_Th[ii] *= xx[ii]*(Length - xx[iCorr])
            else:
                corrTd_Th[ii] *= xx[iCorr]*(Length - xx[ii])

    plt.plot(xx[cellLo:cellHi], corrTd[cellLo:cellHi],’ro’,
             xx[cellLo:cellHi],corrTd_Th[cellLo:cellHi],’b--’)
    plt.xlabel(’Position (m)’); plt.ylabel(’< dT(x) dT(x*) > [delta-function removed]’)
    plt.title(’Average (circles); Theory (dashed)’)
    plt.show()

    # In[ ]:

    #* Compute the measured structure factor and compare with linear theory
    if( NONEQ_FLAG == 0 ):     # Only compute S(k) for equilibrium
        NN = len(range(cellLo,cellHi))
        Sk = np.zeros(NN)                           # Power spectrum S(k), AKA the structure factor
        Sk[0:NN] = sumSk[cellLo:cellHi] / Nsamp     # Average S(k)
        Nyq = np.int_(np.floor(NN/2))               # Nyquist frequency is kSpect[NNNy+1]
        ki = np.arange(1,Nyq)                       # Wavenumber index, skipping the zero index
        Sk_eq = varTheory*NN * np.ones(Nyq+1)       # Equilibrium structure factor, which is constant

        # Plot structure factor, skipping the k=0 wavenumber index
        plt.plot( ki, Sk[1:Nyq],’ro’, ki, Sk_eq[1:Nyq],’b--’,)
        plt.xlabel(’Wavenumber index’); plt.ylabel(’Structure factor S(k)’)
        plt.title(’Average (circles); Theory (dashed)’)
        plt.show()


# Some simplifications in the proof of the Sims conjecture

László Pyber<EMAIL_ADDRESS>Alfréd Rényi Institute of Mathematics, Hungarian
Academy of Sciences, Reáltanoda utca 13-15, H-1053, Budapest, Hungary  Gareth
Tracey<EMAIL_ADDRESS>Alfréd Rényi Institute of Mathematics, Hungarian
Academy of Sciences, Reáltanoda utca 13-15, H-1053, Budapest, Hungary

###### Abstract

We prove an elementary lemma concerning primitive amalgams and use it to
greatly simplify the proof of the Sims conjecture in the case of almost simple
groups.

##  1 Introduction ††The work of the authors on the project leading to this
application has received funding from the European Research Council (ERC)
under the European Unions Horizon 2020 research and innovation programme
(grant agreement No. 741420).

Over the last forty years or so, the Classification of Finite Simple Groups
(henceforth abbreviated to CFSG) has been greatly utilized to resolve a number
of open questions and conjectures in finite group theory. One of the most
well-known of these is the _Sims conjecture_ , first proposed by Charles Sims
in [8]. The conjecture states that if $G$ is a primitive permutation group,
and $h>1$ is the length of a non-trivial orbit of a point stabilizer $H$ in
$G$, then $H$ has order bounded above by a function of $h$.

This conjecture was proved by P.J. Cameron, C.E. Praeger, J. Saxl, and G.M.
Seitz in [1]. Therefore, it is now a theorem, and reads precisely as follows.

######  Theorem 1.1.

There exists a function $f:\mathbb{Z}\rightarrow\mathbb{R}$ such that whenever
$G$ is a primitive permutation group, and $h>1$ is the length of a non-trivial
orbit of a point stabilizer $H$ in $G$, then $|H|\leq f(h)$.

In this note, we present a simplified proof of Theorem 1.1. Our proof does
still require the CFSG, but the more involved methods used in [1, page 502 and
Section 4] are not required.

As a by-product of our approach, we also make progress on a more general
version of the Sims conjecture, which we will refer to as the _non-geometric
Goldschmidt-Sims conjecture_ , following Goldschmidt’s work in [5].

To state this, we first need to introduce the language of _primitive
amalgams_. A triple of finite groups ($H,M,K$), usually written $H>M<K$, is
called a _primitive amalgam_ if $M$ is a subgroup of both $H$ and $K$ and:

  1. (i)

Whenever $A$ is a normal subgroup of $H$ contained in $M$, we have
$N_{K}(A)=M$; and

  2. (ii)

whenever $B$ is a normal subgroup of $K$ contained in $M$, we have
$N_{H}(B)=M$.

The non-geometric Goldschmidt-Sims conjecture reads as follows:

######  The non-geometric Goldschmidt-Sims conjecture.

There exists a function $f:\mathbb{Z}\times\mathbb{Z}\rightarrow\mathbb{R}$
such that whenever $H>M<K$ is a primitive amalgam with $|H:M|=h$ and
$|K:M|=k$, then $|M|\leq f(h,k)$.

We remark that the conjecture as stated above is referred to as the ‘baby
Goldschmidt-Sims problem’ by Fan in [4]. The general ‘Goldschmidt-Sims
problem’ seeks to classify all primitive amalgams, usually under certain
additional conditions (for example, the case where $M$ is a $p$-group and
$|H:M|$, $|K:M|$ are both prime and distinct from $p$ is dealt with in [4]).
We also remark that if ($H,M,K$) is a primitive amalgam, then there is an
associated graph $\Gamma=\Gamma(H,M,K)$ which encodes many properties of the
amalgam. The ‘geometric Goldschmidt-Sims conjecture’ is then a statement about
this graph. See [3, 5] for more details.

Notice that if $G$ is a primitive permutation group, and $G_{\alpha}$ and
$G_{\beta}$ are point stabilizers in $G$, then $G_{\alpha}>G_{\alpha}\cap
G_{\beta}<G_{\beta}$ is a primitive amalgam. Thus, the non-geometric
Goldschmidt-Sims conjecture really is a generalization of the Sims conjecture.

A number of papers (see [2, 3, 6, 10, 12, 13, 11]) have addressed and/or made
progress on the non-geometric Goldschmidt-Sims conjecture. In particular a
result of of Wielandt (published by Knapp in [6, Theorem 2.1]), extending an
earlier result of Thompson [10, Main Theorem], implies the so-called
_Thompson–Wielandt theorem_ (see Theorem 2.1 below) which asserts that if
$H>M<K$ is a primitive amalgam, then there exists a prime $p$ such that
$|H:O_{p}(H)|$ can be bounded in terms of $|H:M|$ and $|K:M|$. We remark that
[6, Theorem 2.1] is only stated in the case where $H$ and $K$ are maximal
subgroups in a finite group $G$, and $M:=H\cap K$ is core-free in $G$. It can
be readily seen, however, that the proof works in the more general setting of
primitive amalgams stated in Theorem 2.1 below.

Our contribution to the non-geometric Goldschmidt-Sims conjecture reads as
follows.

######  Theorem 1.2.

There exists a function $f:\mathbb{Z}\times\mathbb{Z}\rightarrow\mathbb{R}$
such that whenever $H>M<K$ is a primitive amalgam with $|H:M|=h$ and
$|K:M|=k$, and either $H$ or $K$ is abelian by $(h,k)$-bounded, then $|H|\leq
f(h,k)$ and $|K|\leq f(h,k)$.

We remark that in Theorem 1.2, and throughout the paper, the terminology _$X$
is $(n_{1},\ldots,n_{r})$-bounded_, for natural numbers
$X,n_{1},\ldots,n_{r}$, means that $X$ can be bounded by a function of
$n_{1},\ldots,n_{r}$. All other notation used in the paper is standard, though
for a finite simple group $X$ of Lie type, we will use
$\mathop{\mathrm{Inndiag}}(X)$ to denote the subgroup of
$\mathop{\mathrm{Aut}}(X)$ generated by the group $\mathop{\mathrm{Inn}}(X)$
of inner automorphisms of $X$, together with the set of diagonal automorphisms
of $X$. That is,
$\mathop{\mathrm{Inndiag}}(X)=\langle\mathop{\mathrm{Inn}}(X),\alpha\text{ :
}\alpha\text{ a diagonal automorphism of }X\rangle$.

##  2 Proofs of Theorems 1.1 and 1.2

We begin preparations toward the proof of Theorem 1.2, and the simplified
proof of Theorem 1.1, with the Thompson–Wielandt theorem [6, 10] mentioned in
the introduction.

######  Theorem 2.1.

Let $H>M<K$ be a primitive amalgam, with $|H:M|=h$ and $|K:M|=k$. Then there
exists a prime $p$ such that $p$, $|H:O_{p}(H)|$, and $|K:O_{p}(K)|$ are all
$(h,k)$-bounded.

We remark that the proof of the Thompson–Wielandt theorem does not require the
CFSG.

We now discuss the proof of the Sims conjecture in [1]. The proof has five
steps:

  1. (1)

Reduction to the case where the primitive group $G$ in question is almost
simple with socle $X$ a finite group of Lie type.

  2. (2)

Proof that it suffices to bound $|O_{p}(H)\cap\mathop{\mathrm{Inndiag}}(X)|$
in terms of $h$, where $p$ is as in the Thompson–Wielandt theorem.

  3. (3)

Proof that $|H|$ is $h$-bounded if $X$ has characteristic $p$.

  4. (4)

Proof that $|H|$ is $h$-bounded if $X$ is a classical group in $p^{\prime}$
characteristic.

  5. (5)

Proof that $|H|$ is $h$-bounded if $X$ is a finite group of Lie type in
$p^{\prime}$ characteristic, and of $h$-bounded (Lie) rank.

Steps (1), (3), and (4) make use of classical results in group and
representation theory, and the theory of algebraic groups. We will not repeat
the details here, and instead refer the reader to [1, Section 1 and first
three paragraphs of Section 2] for Step (1); the two paragraphs in [1, Section
3] for Step (3); and [1, Section 4a] for Step (4).

Steps (2) and (5) are much more involved, and rely on delicate results of
Seitz [7] from the theory of algebraic groups (see page 502 and Section 4 of
[1] for the details). For this reason, our aim in this note is to replace the
proofs in Steps (2) and (5) with much simpler arguments.

Our next result, part (i) of which is Theorem 1.2, will allow us to do so. We
note that within its proof, we frequently use the standard fact that if $A$,
$B$, and $C$ are subgroups of a group $G$ with $C$ a finite index subgroup of
$B$, then $|A\cap B:A\cap C|\leq|B:C|$. We will also use the notation $J(X)$
for the _Thompson subgroup_ of the finite group $X$. That is, $J(X):=\langle
A<X\text{ : }A\text{ an abelian subgroup of }X\text{ of maximal
order}\rangle$. We remark that the Thompson subgroup is usually only defined
for $p$-groups, but our proofs here require the more general version.

######  Theorem 2.2.

Let $H>M<K$ be a primitive amalgam with $|H:M|=h$ and $|K:M|=k$.

  1. (i)

Suppose that $H$ is abelian by $(h,k)$-bounded. Then $|H|$, and hence $|K|$,
is $(h,k)$-bounded.

  2. (ii)

Suppose that $H$ and $K$ are subgroups of a finite group $G$, and that $I$ is
a normal subgroup of $G$. If $H\cap I$ and $K\cap I$ are abelian by
$(h,k)$-bounded, then the exponents of $H\cap I$ and $K\cap I$ are
$(h,k)$-bounded.

###### Proof.

Adopt the notation of the statement of the theorem, and set $N:=H$, $L:=K$, in
case (i), and $N:=H\cap I$, $L:=K\cap I$, in case (ii).

Let $A$ be an abelian subgroup of $N$ of maximal order. Then $|N:A|$ is
$(h,k)$-bounded, since $N$ is abelian by $(h,k)$-bounded. Suppose now that
$K\cap A$ is core-free in $K$. Then $K$ is a transitive permutation group of
degree $|K:K\cap A|=|K:K\cap N||K\cap N:K\cap A|\leq|K:K\cap N||N:A|$. If we
are in case (i), then $N=H$ and $|K|$ can be bounded in terms of $|K:K\cap
H||H:A|\leq k|H:A|$. Thus, $|K|$, and hence $|H|$, is $(h,k)$-bounded. If we
are in case (ii), then $\mathop{\mathrm{core}}_{K}(L\cap A)=1$. Thus, $K$ is a
permutation group of degree $|K:L\cap A|$. The groups induced by the normal
subgroup $L$ on each of its orbits are all permutation groups of degree
$|L:L\cap A|$, and are all isomorphic to each other. Hence, the exponent of
$L$ is $|L:L\cap A|$-bounded. But $|L:L\cap A|=|L:L\cap N||L\cap N:L\cap
A|=|I\cap K:I\cap M||L\cap N:L\cap A|\leq|K:M||N:A|=k|N:A|$. Thus, $L$ has
$(h,k)$-bounded exponent. Since $|N:N\cap L|=|I\cap H:I\cap M|\leq h$, it
follows that $N$ has $(h,k)$-bounded exponent.

Thus, we may assume that $C:=\mathop{\mathrm{core}}_{K}(K\cap A)$ is non-
trivial. Since $A$ is abelian and $C\leq A$, we have $A\leq N_{H}(C)$, and
since the amalgam is primitive, we have $N_{H}(C)=M$. Thus, $A\leq M$. Hence,
$J(N)$ is contained in $M\leq H$. It follows that $J(H)=J(M)$ in case (i), and
$J(N)=J(M\cap I)$ in case (ii). The same argument with $H$ replaced by $K$
(and $N$ replaced by $L$ in case (ii)), however, also shows that $J(K)=J(M)$
in case (i), and $J(L)=J(M\cap I)$ in case (ii). But then we have, in either
case, that $J(N)=J(L)$ is normal in both $H$ and $K$ \- a contradiction. This
completes the proof. ∎

We remark that the proof of Theorem 2.2 does not require the Thompson–Wielandt
theorem. Also, part (i) above proves Theorem 1.2.

We can now replace Steps (2) and (5) in the proof of the Sims conjecture with
the following easy consequences of Theorem 2.2(i). We begin with Step (2):

######  Corollary 2.3.

Let $G$ be an almost simple primitive permutation group, with point stabilizer
$H$, and let $h>1$ be the length of a non-trivial $h$-orbit. Suppose that the
socle $X$ of $G$ is a finite group of Lie type, and let $p$ be the prime from
the Thompson–Wielandt theorem. If $|O_{p}(H)\cap\mathop{\mathrm{Inndiag}}(X)|$
is $h$-bounded, then $|H|$ is $h$-bounded.

###### Proof.

Since $\mathop{\mathrm{Aut}}(X)/\mathop{\mathrm{Inndiag}}(X)$ has a normal
subgroup of order at most $6$ with cyclic quotient, the Thompson–Wielandt
theorem and the hypothesis of the corollary guarantee that $H$ is cyclic by
$h$-bounded. Theorem 2.2(i) then yields the result. ∎

Next, Step (5):

######  Corollary 2.4.

Let $G$ be an almost simple primitive permutation group, with point stabilizer
$H$, and let $h>1$ be the length of a non-trivial $h$-orbit. Let $p$ be the
prime from the Thompson–Wielandt theorem, and suppose that the socle $X$ of
$G$ is a finite group of Lie type in $p^{\prime}$ characteristic, and that $X$
has $h$-bounded (Lie) rank. Then $|H|$ is $h$-bounded.

###### Proof.

Write $h=|H:H\cap H^{g}|$, for an element $g$ of $G$, and set $K:=H^{g}$,
$M:=H\cap K$, and $I:=\mathop{\mathrm{Inndiag}}(G)$. Also, we may write
$X=O^{p^{\prime}}(\overline{X}_{\sigma})$ for a simple algebraic group
$\overline{X}$ of adjoint type and a Steinberg endomorphism $\sigma$ of
$\overline{X}$. Since $p$ is coprime to the defining characteristic of $X$,
[9, II, Theorem 5.16] implies that $O_{p}(H)\cap I$ normalizes of a maximal
torus $\overline{T}$ in $\overline{X}$. Then $(O_{p}(H)\cap I)/(O_{p}(H)\cap
I\cap\overline{T})$ is a subgroup of the Weyl group $W$ of $X$. Since the Lie
rank of $X$ is $h$-bounded, the order of the group $W$ is $h$-bounded. We
deduce that (a) $O_{p}(H)\cap I$ is abelian by $h$-bounded; (b) $O_{p}(H)\cap
I$ has $h$-bounded subgroup rank; and (c) $O_{p}(H)\cap I$ has bounded derived
length.

It follows from (a) and the Thompson–Wielandt theorem that $H\cap I$ and
$K\cap I=(H\cap I)^{g}$ are abelian by $h$-bounded. Since $H>M<K$ is a
primitive amalgam, Theorem 2.2(ii) then implies that $H\cap I$ and $K\cap I$
have $h$-bounded exponents. This, together with (b) and (c) above, implies
that $|O_{p}(H)\cap I|$ and $|O_{p}(K)\cap I|$ are $h$-bounded. The result now
follows from Corollary 2.3. ∎

This completes our simplified proof of the Sims conjecture.

## References

  * [1] P.J. Cameron, C.E. Praeger, J. Saxl, and G.M. Seitz. On the Sims conjecture and distance transitive graphs. Bull. London Math. Soc. 15 (1983), 499–506.
  * [2] P. S. Fan. Amalgams of finite groups and the Goldschmidt-Sims conjecture. Proceedings of the Rutgers group theory year, 1983-1984 (ed. M. Aschbacher et al. CUP, Cambridge, 1984), 161–166.
  * [3] P. S. Fan. The Thompson-Wielandt theorem. Proc. Amer. Math. Soc. 97 (1986), 590–592.
  * [4] P. S. Fan. Amalgams of prime index. J. Algebra. 98 (1986), 375–421.
  * [5] D. Goldschmidt. Automorphisms of trivalent graphs. Ann. of Math. 111 (1980), 377–406.
  * [6] W. Knapp. On the point stabilizer in a primitive permutation group. Math. Z. 133 (1973), 137–168.
  * [7] G. M. Seitz. The root subgroups for maximal tori in finite groups of Lie type. Pacific J. Math. 106 (1983), 153–244.
  * [8] C.C. Sims. Graphs and finite permutation groups. Math. Z. 95 (1967), 76–86.
  * [9] T. A. Springer and R. Steinberg. Conjugacy classes. Seminar on algebraic groups and related finite groups (eds A. Borel et al). Lecture Notes in Mathematics 131 (Springer, Berlin, 1970), 168–266.
  * [10] J.G. Thompson. Bounds for the orders of maximal subgroups. J. Algebra 14 (1970), 135–138.
  * [11] J. Van Bon. Thompson-Wielandt Theorems revisited. Bull. London Math. Soc. 35 (2003), 30–36.
  * [12] R. Weiss, On a theorem of Goldschmidt. Ann. of Math. 126 (1987), 429–438.
  * [13] H. Wielandt. Subnormal subgroups and permutation groups. Lecture notes (ed. F. Demana, W. McWorter and S. Seghal, Ohio State University, Columbus, Ohio, 1971).