The relief is a random surface that evolves in time: a stationary Gaussian process with values in the distributions on the torus, the solution of a damped stochastic wave equation driven by a cylindrical Brownian motion.

Setting

Let $\mathbb{T}^2=\mathbb{R}^2/\mathbb{Z}^2$ and let $\mathcal{H}=L^2_0(\mathbb{T}^2)$ be the space of square-integrable functions with mean zero, with $\langle f,g\rangle=\int_{\mathbb{T}^2}fg\,dx$. Let $\mathbb{Z}^2_+=\{k\in\mathbb{Z}^2:k_2>0\}\cup\{k\in\mathbb{Z}^2:k_2=0,\ k_1>0\}$ and, for $k\in\mathbb{Z}^2\setminus\{0\}$, $$e_k(x)=\begin{cases}\sqrt2\cos(2\pi k\cdot x), & k\in\mathbb{Z}^2_+,\\ \sqrt2\sin(2\pi k\cdot x), & -k\in\mathbb{Z}^2_+.\end{cases}$$ Then $(e_k)_{k\neq0}$ is an orthonormal basis of $\mathcal{H}$ and $-\Delta e_k=4\pi^2|k|^2e_k$.

Noise and coefficients

Let $W$ be a cylindrical Wiener process on $\mathcal{H}$, that is, $\langle W_t,\varphi\rangle=\sum_{k\neq0}\langle e_k,\varphi\rangle\,W^k_t$ with independent standard Brownian motions $(W^k)_{k\neq0}$. For $H\in(0,1)$ and $\omega_0>0$ let $$\Gamma=\omega_0\Big(-\frac{\Delta}{4\pi^2}\Big)^{1/4},\qquad B=2\,\Gamma^{3/2}\Big(-\frac{\Delta}{4\pi^2}\Big)^{-(1+H)/2}.$$ Then $\Gamma e_k=\gamma_ke_k$ and $Be_k=\beta_ke_k$ with $\gamma_k=\omega_0|k|^{1/2}$ and $\beta_k=2\gamma_k^{3/2}|k|^{-(1+H)}$.

Equation

The field $X$ solves the damped stochastic wave equation $\ddot X=-\Gamma^2X-2\Gamma\dot X+B\dot W$ in the weak sense: for every $\varphi\in C^\infty(\mathbb{T}^2)$ with mean zero and every $t\ge0$, $$\begin{aligned}\langle X_t,\varphi\rangle&=\langle X_0,\varphi\rangle+\int_0^t\langle\dot X_s,\varphi\rangle\,ds,\\ \langle\dot X_t,\varphi\rangle&=\langle\dot X_0,\varphi\rangle-\int_0^t\big(\langle X_s,\Gamma^2\varphi\rangle+2\langle\dot X_s,\Gamma\varphi\rangle\big)\,ds\\ &\quad+\sum_{k\neq0}\beta_k\int_0^t\langle e_k,\varphi\rangle\,dW^k_s.\end{aligned}$$ The initial value $(X_0,\dot X_0)$ is independent of $W$ and centred Gaussian with independent coordinates, $\operatorname{Var}\langle X_0,e_k\rangle=|k|^{-2(1+H)}$ and $\operatorname{Var}\langle\dot X_0,e_k\rangle=\gamma_k^2|k|^{-2(1+H)}$.

Law

The solution is stationary and centred Gaussian. With two-sided Brownian motions $W^k$, $$\langle X_t,\varphi\rangle=\sum_{k\neq0}\beta_k\int_{-\infty}^t(t-r)\,e^{-\gamma_k(t-r)}\,\langle e_k,\varphi\rangle\,dW^k_r.$$ For $s,t\in\mathbb{R}$, $k,l\neq0$ and $x,y\in\mathbb{T}^2$, $$\mathbb{E}\big[\langle X_s,e_k\rangle\langle X_t,e_l\rangle\big]=\delta_{kl}\,|k|^{-2(1+H)}\big(1+\gamma_k|t-s|\big)\,e^{-\gamma_k|t-s|},$$ $$\mathbb{E}\big[X_t(x)\,X_t(y)\big]=\sum_{k\neq0}|k|^{-2(1+H)}\cos\big(2\pi k\cdot(x-y)\big).$$ $X_t$ is Hölder continuous of every order below $H$, and $t\mapsto\langle X_t,\varphi\rangle$ is of class $C^1$. $\dot X_t\in\mathcal{H}$ if and only if $H>1/2$; otherwise $\dot X_t$ is a distribution. $B$ is Hilbert–Schmidt on $\mathcal{H}$ if and only if $H>3/4$. As $H\to0$ the law of $X_t$ tends to the Gaussian free field.

References

How the picture is computed →

← The model

Scale

The torus is drawn with side $64$; heights are in the same unit. Time $t$ is in seconds, and $\omega_0$ is set in $\mathrm{s}^{-1}$ (default $0.018$). The code keeps $\omega_0=0.018\,\mathrm{s}^{-1}$ in $\Gamma$ and $B$ and runs its clock $c=\omega_0/0.018$ times as fast. This gives the same law: $t\mapsto X_{ct}$ solves the equation with $\omega_0$ replaced by $c\,\omega_0$, since $B$ is proportional to $\Gamma^{3/2}$. For $\omega_0=0$ the picture stands still.

Truncation

Only the coordinates $\langle X_t,e_k\rangle$ and $\langle\dot X_t,e_k\rangle$ with $k\in\Lambda=\{k\in\mathbb{Z}^2:2\le|k|\le32\}$ are computed ($|\Lambda|=3200$). Since $\Gamma$ and $B$ are diagonal in $(e_k)$, these coordinates are independent processes, and the computed field is exactly the orthogonal projection $P_\Lambda X_t$ onto $\operatorname{span}\{e_k:k\in\Lambda\}$. The coordinates with $|k|<2$ are left out: at $H=1/2$ they carry $61\%$ of the variance of $P_{\{|k|\le32\}}X_t$, and one wave across the torus would dominate the picture. The bound $|k|\le32$ keeps the shortest wavelength at $8$ grid cells.

Time

Over a step $\Delta=1/30\,\mathrm{s}$ each pair $\big(\langle X_t,e_k\rangle,\langle\dot X_t,e_k\rangle\big)$ moves by its exact Gaussian transition, with mean $e^{A_k\Delta}$ applied to the pair and covariance $\int_0^\Delta e^{A_ku}b_kb_k^{\top}e^{A_k^{\top}u}\,du$, where $$A_k=\begin{pmatrix}0&1\\-\gamma_k^2&-2\gamma_k\end{pmatrix},\qquad b_k=\begin{pmatrix}0\\ \beta_k\end{pmatrix}.$$ There is no discretisation error. The initial value is drawn from the stationary law; a change of $H$ rescales the state so that it stays stationary.

Space

Every $0.12\,\mathrm{s}$, $P_\Lambda X_t$ is evaluated on a $256\times256$ grid by an inverse FFT. Between two evaluations the grid values are interpolated linearly in time, and between grid points bilinearly.

Relief

Let $\sigma_\Lambda^2=\mathbb{E}\big[(P_\Lambda X_t)(x)^2\big]=\sum_{k\in\Lambda}|k|^{-2(1+H)}$ and $Z_t(x)=(P_\Lambda X_t)(x)/\sigma_\Lambda$. The drawn height is $F(Z_t(x))$ with $F(z)=L\big(\rho(G(z)-G(-3))\big)$, $$\begin{aligned}G(z)&=6\,\frac{e^{gz}-1}{e^{3g}-1},\\ \rho(u)&=\epsilon\log\big(1+e^{u/\epsilon}\big),\\ L(y)&=y+\ell T\Big(1-\frac{y}{T}\Big)_+^2,\end{aligned}$$ where $T=\rho\big(G(3)-G(-3)\big)$, $g$ is the setting sharp ($G(z)=2z$ for $g=0$), $\epsilon=0.35$ and $\ell=0.18$. $G$ makes the peaks steeper and the valleys flatter, $\rho$ is a smooth floor at $0$, and $L$ lifts the valleys and leaves the height $T$ unchanged. $F$ is increasing, so the peaks and valleys of the picture are those of $P_\Lambda X_t$.

Camera

At horizontal distance $d$ from the viewer the drawn height is multiplied by $w(d)=0.45+0.55\,S\big((d-6)/30\big)$, where $S(u)=3u^2-2u^3$ for $u\in[0,1]$, $S(u)=0$ for $u<0$ and $S(u)=1$ for $u>1$.

Code

In js/field.js, for $k\in\Lambda\cap\mathbb{Z}^2_+$, cRe and cIm hold $\langle X_t,e_k\rangle$ and $\langle X_t,e_{-k}\rangle$, and vRe and vIm the same coordinates of $\dot X_t$. The grid holds $P_\Lambda X_t/\sqrt2$; the factor cancels in $Z_t$.