Dynamical universality classes of simple growth and lattice gas models

TL;DR

This study uses large-scale simulations of the 2D octahedron model to analyze surface growth universality classes, confirming theoretical conjectures and revealing finite size effects in KPZ and EW models.

Contribution

It provides high-precision autocorrelation data, compares different stochastic dynamics, and supports the KPZ ansatz and Kallabis-Krug conjecture in 2+1 dimensions.

Findings

Confirmed the KPZ growth exponent as 0.2414(2).

Compared stochastic dynamics and their effects on correlations.

Identified finite size corrections occurring before steady state.

Abstract

Large scale, dynamical simulations have been performed for the two dimensional octahedron model, describing the Kardar-Parisi-Zhang (KPZ) for nonlinear, or the Edwards-Wilkinson (EW) class for linear surface growth. The autocorrelation functions of the heights and the dimer lattice gas variables are determined with high precision. Parallel random-sequential (RS) and two-sub-lattice stochastic dynamics (SCA) have been compared. The latter causes a constant correlation in the long time limit, but after subtracting it one can find the same height functions as in case of RS. On the other hand the ordered update alters the dynamics of the lattice gas variables, by increasing (decreasing) the memory effects for nonlinear (linear) models with respect to random-sequential. Additionally, we support the KPZ ansatz and the Kallabis-Krug conjecture in $2+1$ dimensions and provide a precise growth exponent value $\beta=0.2414(2)$. We show the emergence of finite size corrections, which occur long before the steady state roughness is reached.

Tables

Table 1. Table 1: Autocorrelation limits for KPZ with SCA dynamics for different deposition rates p 𝑝 p and q = 0 𝑞 0 q=0 , as functions of the waiting time s 𝑠 s . Fit errors are shown, which are below the given number of digits in case of the slopes.

$s/\text{SI}MCS$	$o_h$		$o_s$
	$p=\text{num}.5$	$p=\text{num}.95$	$p=\text{num}.5$	$p=\text{num}.95$
30	\num0.00320(3)	\num0.055398(8)	\num0.012871	\num0.221623
100	\num0.00331(5)	\num0.05520(3)	\num0.014286	\num0.219827
500	\num0.0031(2)	\num0.05457(8)	\num0.013944	\num0.218547
1000	\num0.0035(3)	\num0.0548(2)	\num0.013903	\num0.218330

Table 2. Table 2: Summary of KPZ and EW autocorrelation λ C subscript 𝜆 𝐶 \lambda_{C} and aging b 𝑏 b exponents, assuming z KPZ = \num ​ 1.611 ​ ( 3 ) subscript 𝑧 KPZ \num 1.611 3 z_{\mathrm{KPZ}}=\num{1.611}(3) and z EW = 2 subscript 𝑧 EW 2 z_{\mathrm{EW}}=2 , respectively. There are no independent estimates for b 𝑏 b under SCA dynamics. Values for the EW case provided without error margin correspond to the analytical solution  [ 48 , 49 ] .

		${\lambda }_{C,h}$	${\lambda }_{C,s}$	$b_h$	$b_s$
KPZ	RS	\num1.98(5)	3.8(2)	\num-.4828(4)	\num0.76(2)
	SCA	\num2.01(2)	\num1.25(2)
EW	RS	2	$\text{num}1.4\left(4\right)$	0	1.1(2)
	SCA	2	$\approx \text{num}4$

Table 3. Table 3: Summary of fitting of β eff subscript 𝛽 eff \beta_{\mathrm{eff}} using the KPZ ansatz in the most general form ( 32 ) and using ( 33 ) with even powers, more suitable to describe corrections of the roughness-scaling up to maximum orders N = 6 𝑁 6 N=6 . Error margins given are uncertainties of the fit parameters, not actual error estimates for β 𝛽 \beta . The reduced sums of residuals χ red subscript 𝜒 red \chi_{\mathrm{red}} are given to judge the quality of the fits. For each dataset we underlined the value of χ red subscript 𝜒 red \chi_{\mathrm{red}} , where the quality does not improve significantly by increasing N 𝑁 N .

	random-sequential (RS)				stochastic cellular automaton (SCA)
	DTrDB, TC=1,1		DTrDB, TC=2,2		$p=\text{num}.5$		$p=\text{num}.75$		$p=\text{num}.95$
$N$	$\beta$	${\chi }_{\mathrm{red}}$	$\beta$	${\chi }_{\mathrm{red}}$	$\beta$	${\chi }_{\mathrm{red}}$	$\beta$	${\chi }_{\mathrm{red}}$	$\beta$	${\chi }_{\mathrm{red}}$
equation (32)
$1$	$\text{num}0.2430\left(\text{num}1\right)$	\num10.04	$\text{num}0.2433\left(\text{num}1\right)$	\num9.24	$\text{num}0.2420\left(\text{num}1\right)$	\num7.03	$\text{num}0.2419\left(\text{num}1\right)$	\num8.63	$\text{num}0.2418\left(\text{num}1\right)$	\num3.15
$2$	$\text{num}0.2396\left(\text{num}1\right)$	\num1.82	$\text{num}0.2396\left(\text{num}1\right)$	\num1.29	$\text{num}0.2403\left(\text{num}1\right)$	\num1.62	$\text{num}0.2400\left(\text{num}1\right)$	\num1.31	$\text{num}0.2403\left(\text{num}1\right)$	\num0.51
$3$	$\text{num}0.2411\left(\text{num}1\right)$	\num0.79	$\text{num}0.2408\left(\text{num}1\right)$	\num0.55	$\text{num}0.2414\left(\text{num}1\right)$	\num0.64	$\text{num}0.2404\left(\text{num}1\right)$	\num1.18	$\text{num}0.2401\left(\text{num}1\right)$	\num0.49
$4$	$\text{num}0.2421\left(\text{num}2\right)$	\num0.67	$\text{num}0.2421\left(\text{num}1\right)$	\num0.32	$\text{num}0.2414\left(\text{num}1\right)$	\num0.64	$\text{num}0.2385\left(\text{num}1\right)$	\num0.63	$\text{num}0.2403\left(\text{num}2\right)$	\num0.48
$5$	$\text{num}0.2386\left(\text{num}3\right)$	\num0.47	$\text{num}0.2409\left(\text{num}2\right)$	\num0.28	$\text{num}0.2388\left(\text{num}3\right)$	\num0.51	$\text{num}0.2394\left(\text{num}3\right)$	\num0.61	$\text{num}0.2399\left(\text{num}5\right)$	\num0.48
$6$	$\text{num}0.2377\left(\text{num}8\right)$	\num0.47	$\text{num}0.2403\left(\text{num}5\right)$	\num0.28	$\text{num}0.2333\left(\text{num}9\right)$	\num0.45	$\text{num}0.2420\left(\text{num}7\right)$	\num0.59	$\text{num}0.2485\left(\text{num}8\right)$	\num0.40
equation (33)
$1$	$\text{num}0.2421\left(\text{num}1\right)$	\num7.28	$\text{num}0.2422\left(\text{num}1\right)$	\num6.27	$\text{num}0.2416\left(\text{num}1\right)$	\num5.22	$\text{num}0.2415\left(\text{num}1\right)$	\num6.54	$\text{num}0.2412\left(\text{num}1\right)$	\num1.76
$2$	$\text{num}0.2414\left(\text{num}1\right)$	\num1.16	$\text{num}0.2414\left(\text{num}1\right)$	\num1.10	$\text{num}0.2411\left(\text{num}1\right)$	\num0.65	$\text{num}0.2410\left(\text{num}1\right)$	\num1.72	$\text{num}0.2409\left(\text{num}1\right)$	\num0.81
$3$	$\text{num}0.2412\left(\text{num}1\right)$	\num0.65	$\text{num}0.2412\left(\text{num}1\right)$	\num0.33	$\text{num}0.2411\left(\text{num}1\right)$	\num0.64	$\text{num}0.2409\left(\text{num}1\right)$	\num1.51	$\text{num}0.2408\left(\text{num}1\right)$	\num0.52
$4$	$\text{num}0.2413\left(\text{num}1\right)$	\num0.53	$\text{num}0.2412\left(\text{num}1\right)$	\num0.29	$\text{num}0.2412\left(\text{num}1\right)$	\num0.59	$\text{num}0.2407\left(\text{num}1\right)$	\num0.86	$\text{num}0.2407\left(\text{num}1\right)$	\num0.50
$5$	$\text{num}0.2412\left(\text{num}1\right)$	\num0.51	$\text{num}0.2412\left(\text{num}1\right)$	\num0.28	$\text{num}0.2412\left(\text{num}1\right)$	\num0.59	$\text{num}0.2406\left(\text{num}1\right)$	\num0.68	$\text{num}0.2405\left(\text{num}1\right)$	\num0.41
$6$	$\text{num}0.2412\left(\text{num}1\right)$	\num0.51	$\text{num}0.2412\left(\text{num}1\right)$	\num0.28	$\text{num}0.2411\left(\text{num}1\right)$	\num0.57	$\text{num}0.2403\left(\text{num}1\right)$	\num0.50	$\text{num}0.2407\left(\text{num}1\right)$	\num0.38

Equations

$${\partial }_{t}h(\mathbf{r},t)=\sigma {\nabla }^{2}h(\mathbf{r},t)+\lambda (\nabla h(\mathbf{r},t){)}^2+\eta (\mathbf{r},t),$$

$$W^{2}(L,t)=\frac{1}{L^{2}}\,\sum_{i,j}^{L}\,h^{2}_{i,j}(t)-\Bigl{(}\frac{1}{L}\,\sum_{i,j}^{L}\,h_{i,j}(t)\Bigr{)}^{2}\ .$$

$$W(L,t)\propto L^{\alpha }f(t/L^z),$$

$$f(u)\propto\left\{\begin{array}[]{lcl}u^{\beta}&{\rm if}&u\ll 1\\ {\rm const.}&{\rm if}&u\gg 1\end{array}\right.$$

$$z=\alpha/\beta\$$

$${\Phi }_L^n[{\phi }_L]$$

$$S_L[{\phi }_L]$$

$$Q_L[{\phi }_L]$$

$${\Psi }_L[W^2(L)]$$

$${\Psi }_L[h_L(\mathbf{r})]$$

$$C(t,s)$$

$$\propto s^{-b}(t/s{)}^{-{\lambda }_C/z},$$

$$C_h(t,s)$$

$$=⟨h^2(s,\mathbf{r})⟩-{⟨h(s,\mathbf{r})⟩}^2=W^2(L\to \infty ,s)\propto s^{-b_h}\cdot f_C(1).$$

$$b_h=-2\beta ,$$

$$C_s(t,s)$$

$$=⟨n(t,\mathbf{r})n(s,\mathbf{r})⟩-{\overline{n}}^2$$

$$=s^{-b_s}{f}_C^{\prime }\left(\frac{t}{s}\right),$$

$$\theta =(D+4)/z-2.$$

$$R(t,s)={\frac{\delta ⟨\varphi (t)⟩}{\delta j(s)}|}_{j=0}=s^{-1-a}{f}_R\left(\frac{t}{s}\right)$$

$$\left(\begin{array}[]{cc}0&1\\ 0&1\end{array}\right)\overset{p}{\underset{q}{\rightleftharpoons}}\left(\begin{array}[]{cc}1&0\\ 1&0\end{array}\right)\quad,$$

$$h_{i,j}=\sum _{l=1}^i[2{\sigma }_x(l,1)-1]+\sum _{k=1}^j[2{\sigma }_y(i,k)-1].$$

$${\beta }_{\mathrm{eff}}\left(\frac{t_i-t_{i/2}}{2}\right)=\frac{\ln W(L\to \infty ,t_i)-\ln W(L\to \infty ,t_{i/2})}{\ln (t_i)-\ln (t_{i/2})}.$$

$$t_{i+1}=(t_i+10){\mathrm{e}}^m,$$

$$C_{\mathrm{h}}^{\mathrm{EW}}=c_0\ln \left(\frac{t+s}{t-s}\right),$$

$$f_{\mathrm{SCA}}(t,p)=\tilde{t}(p)=t\cdot p{\mathrm{e}}^p.$$

$$R(t)=R_{\infty }+a_R/{\tilde{t}}^{2\beta }+b_R/{\tilde{t}}^{4\beta },$$

$$h(t)$$

$${\beta }_{\mathrm{eff}}$$

$$2=\alpha +z=\alpha \left(1+1/\beta \right).$$

$$f(t)=o+c\cdot t^{-x},$$

$$f_{C,\text{lx@glossaries@gls@link}\mathrm{acronym}\mathrm{sca}\text{leavevmode}\mathrm{SCA}}(t/s,p)\propto c_1\cdot (t/s{)}^{-{\lambda }_{C,\mathrm{s}}/z}+c_2\cdot (t/s{)}^{-{\lambda }_{C,\mathrm{s},\text{lx@glossaries@gls@link}\mathrm{acronym}\mathrm{sca}\text{leavevmode}\mathrm{SCA}}/z}$$

$$C_q(t,s)=⟨h_{-q}(t)h_q(s)⟩.$$

$$T\chi (t,s;r)=-{\partial }_r^{2}C(t,s;r).$$

$$1+a=b+2/z$$

Full text

Keywords: driven lattice gas, surface growth, autocorrelation, Kardar–Parisi–Zhang class, stochastic cellular automaton, Edwards–Wilkinson class

1 Introduction

Nonequilibrium systems are known to exhibit dynamical scaling, when the correlation length diverges as $\xi\propto t^{1/z}$, characterized by the exponent $z$. Simplest models are driven lattice gases (DLG) [1], which in certain cases can be mapped onto surface growth [2, 3]. Therefore, understanding DLGs, which is far from being trivial due to the broken time reversal symmetry [4], and is possible mostly by numerical simulations only, sheds some light on the corresponding interface phenomena [5]. The simplest example is the asymmetric simple exclusion process (ASEP) of particles [6], in which particles and holes can be mapped onto binary surface slopes [7, 8] and the corresponding continuum model can be described by the Kardar–Parisi–Zhang (KPZ) equation [9]

[TABLE]

where the scalar field $h(\mathbf{r},t)$ is the height, progressing in the $D$ dimensional space relative to its mean position, that moves linearly with time $t$. This equation was inspired in part by the stochastic Burgers equation [10] and can describe the dynamics of simple growth processes in the thermodynamic limit [11], randomly stirred fluids [12], directed polymers in random media [13], dissipative transport [14, 15], and the magnetic flux lines in superconductors [16]. In case of surface growth $\sigma$ represents a surface tension, competing with the nonlinear up–down anisotropy of strength $\lambda$ and a zero mean valued Gaussian white noise $\eta$. This field exhibits the covariance $\langle\eta(\mathbf{r},t)\eta(\mathbf{x^{\prime}},t^{\prime})\rangle=2\Gamma\delta^{D}(\mathbf{r-r^{\prime}})(t-t^{\prime})$. The $\lambda=0$, linear equation describes the Edwards–Wilkinson (EW) [17] surface growth, an exactly solvable equilibrium system.

Several discrete models obeying these equations have been studied [7, 18, 2]. The morphology of a surface of linear size $L$ is usually described by the squared interface width

[TABLE]

In the absence of any characteristic length simple growth processes are expected to be scale-invariant [19]

[TABLE]

with the universal scaling function $f(u)$:

[TABLE]

Here $\alpha$ is the roughness exponent in the stationary regime, when the correlation length $\xi$ has grown to exceeded $L$, and $\beta$ is the growth exponent, describing the intermediate time behavior. The dynamical exponent $z$ can be expressed as the ratio of the growth exponents:

[TABLE]

Apart from the exponents, the shapes of the rescaled width and height distributions of the interfaces $\Psi_{L}(\varphi_{L})$ were shown to be universal in KPZ models in both the steady state [20] and the growth regime [21]. Here, $\varphi_{L}$ denotes the interface observable in question, $W^{2}$ or $h$, in a system of linear size $L$. In fact many people define the universality classes by these quantities, which can be obtained exactly in one dimension for various surface geometries, like flat [22] or curved [23, 24, 22] interfaces. The non-rescaled probability distributions are denoted by $P_{L}(\varphi_{L})$ and their moments are defined via the distribution averages as:

[TABLE]

for the surface height. Note, that $\left\langle h_{L}\right\rangle\equiv\Phi_{L}^{0}[h_{L}]\equiv 0$ in the co-moving frame of the surface.

While many systems are described by a single dynamical length scale, aging ones are best characterized by two-time quantities, such as the dynamical correlation and response functions [25]. In the aging regime: $s\gg\tau_{\rm m}$ and $t-s\gg\tau_{\rm m}$, where $\tau_{\rm m}$ is a microscopic time scale and $s$ is the start time, when the snapshot is taken, one expects the following law for the autocorrelation function

[TABLE]

here $\langle\rangle$ denotes averaging over both lattice sites and independent samples; $\lambda_{C}$ is the autocorrelation and $b$ is the aging exponent. The function $\phi$ denotes the measured quantity, which can be the particle density of the lattice gas or the surface height $h(t,\mathbf{r})$. In the latter case,

[TABLE]

This implies the relation

[TABLE]

which must be satisfied in the $L\to\infty$ and $s\to\infty$ limit. We have also calculated the auto-correlation of the slope (lattice gas occupancy variables) $n(t,\mathbf{r})$ as:

[TABLE]

where $\overline{n}=\num{0.5}$ is the conserved average occupancy of sites. However, $C_{s}(t,s)$ decays much faster than the height auto-correlator and obtaining reasonable signal/noise ratio requires much higher statistics.

A dynamic, perturbative renormaliztion group (RG) analysis of the KPZ equation [26, 27] suggested that the short and the long time scaling behavior of the height correlation function are identical and deduced a scaling relation for the exponent of $C_{h}(t\gg s,s\to 0)\propto(t/s)^{-\theta}$ as:

[TABLE]

Since $\theta=\lambda_{C}/z+2\beta$, due to (14) the relation (16) holds exactly in the solvable $1+1$ dimensional case. In $D\geq 2$ dimensions perturbative RG can’t access the strong coupling KPZ fixed point [28], thus the validity of this law should be tested by precise exponent estimates.

A conjecture based on a purely geometric argument, advanced by Kallabis and Krug [29], which can also be deduced from the scaling relation (16) claims $D=\lambda_{C}$ in any dimensions. In [30] we provided marginal agreement for this in $D=2$ and Ref. [31] also suggested it using solid-on-solid models, although clear power laws were not reached within the times studied. Here we provide stronger numerical evidence in case of different lattice gas dynamics.

In aging systems a similar scaling form is expected for the autoresponse function of the field $\phi$:

[TABLE]

where $j$ is the external conjugate to $\phi$ and $a$ denotes the so-called aging exponent $a$. The universal scaling function exhibits the asymptotic behavior $f_{R}(t/s)\sim(t/s)^{-\lambda_{R}/z}$ with the autoresponse exponent $\lambda_{R}$. In equilibrium $\lambda_{C}=\lambda_{R}$ and $a=b$ due to the fluctuation-dissipation (FD) symmetry [32]. In nonequilibrium systems these exponents can be completely independent. Therefore, we shall determine them one-by-one and investigate if some extended FD relation may occur among them. This has been done, using our very recent aging response exponents [33], determined to test the validity of a logarithmic extension of local scale-invariance (LSI) [32] proposed in [34] and work on such extensions in other models has been continued recently [35, 36].

Throughout the present study we compare results obtained using two common updating schemes for lattice models: random-sequential (RS) and stochastic cellular automaton (SCA) (checkerboard) updates. We find constant as well as non-trivial corrections to the dynamical correlation functions produced by the SCA and also observe differences in the corrections to scaling.

This paper is structured as follows. The investigated model and the simulation algorithms are introduced in Sect. 2. Roughness growth results are presented in Sect. 3.1, while autocorrelation and aging date can be found in Sect. 3.2. We conclude the paper with a discussion of the main implications of our results in Sect. 4.

2 Models and simulation algorithms

Discrete models set up for KPZ have been studied a lot in the past decades [7, 18, 2]. A mapping between KPZ surface growth in two dimensions and driven lattice gases has been advanced in [37, 38] an extension of the "rooftop" model of [7, 8]. We called it octahedron model, characterized by binary slope variables $\sigma_{x/y}$ at the edges connecting top vertexes of octahedra [38] representing atoms. The $\sigma_{x/y}$ take the values [math] or $1$ to encode down or up slopes, respectively. Thus deposition or removal of octahedra corresponds to a stochastic cellular automaton, with the simple Kawasaki update rules [38]

[TABLE]

where $p$ and $q$ denote the acceptance probabilities. Projecting the edges onto a plane yields a square lattice of slopes, which can then be considered as occupancy variables. This maps the octahedron model onto self-reconstructing dimers following an oriented migration along the bisection of the $x$ and $y$ directions of the surface (see figure 1 or figure 11 of the supplemental material for a 3D depiction). In this picture the surface heights must be defined relative to a reference point $h_{1,1}=0$ and can be reconstructed from the slope variables as

[TABLE]

Discrete surface and DLG models usually apply random sequential dynamics. On the other hand in certain cases synchronous, so called SCA-like site updating can prove to be useful, especially for simulations on parallel computers. This study is based on massively parallel simulations on graphics cards (GPUs). Synchronous updating in case of one-dimensional ASEP models has already been investigated [39, 40, 41]. One-point quantities in the bulk, like particle current or surface growth have been shown to exhibit the same behavior as in case of RS. However, n-point correlation functions may be different.

Here we extend the parallel two-sublattice scheme developed for ASEP [40] to the two dimensional dimer model as shown on figure 1, and compare the dynamical scaling results with those of the RS dynamics.

While the latter produces uncorrelated deposition and removal processes, SCA dynamics attempts updates in a checkerboard pattern, which are thusly correlated. Because of this, blocks of sites to be updated can be visited in a sequential order within a SCA sub-lattice step, allowing for very efficient implementation [42], matching perfectly parallel processors of GPU architectures [43].

Performing RS simulations on GPUs is less straight-forward, because unwanted correlations may be introduced [30]. In order to eliminate these and to achieve results as close to really sequential simulations as possible, we apply a new DD scheme, with two layers of DD to match the GPU architecture. At level one, in a double tiling decomposition, the origin is moved randomly after each sweep of the lattice (DTr). At level two, these tiles are subdivided further, with a logical dead border (DB) scheme. Collective update attempts are preformed inside these cells, excluding one lattice site-wide borders around each. Here, the decomposition origin is moved randomly after each collective update attempt. This scheme will be referred to as DTrDB in the following. Details of the new implementation are documented elsewhere [44].

In order to estimate the asymptotic values of different exponents for $t\to\infty$, local slope analyses of the scaling laws were performed [5]. For example in case of the interface width growth we used

[TABLE]

In our studies the simulation time, measured in Monte Carlo steps (MCS), between two measurements was increased exponentially

[TABLE]

using $m=\num{.01}$ and $t_{0}=0$. A flat initial state is realized by a zig-zag pattern with $W^{2}(L,t_{0})=\num{.25}$. The simulations are subject to periodic boundary conditions.

Statistical uncertainties are provided as $1\sigma$–standard errors, defined as $\Delta_{1\sigma}x=\sqrt{\langle x^{2}\rangle-\langle x\rangle^{2}/(N-1)}$. Throughout this study we used the implementation of the Levenberg–Marquadt algorithm [45, 46] in the gnuplot software [47] for non-linear least squares fitting.

For $p=q>0$ the octahedron model describes the surface growth of the EW equation [38] in $2+1$ dimensions. In this case the autocorrelation function of heights has been derived [48, 49]:

[TABLE]

where $c_{0}$ is a model-dependent constant. This function approaches [math] for $t\gg s$ as a power law (PL) with the exponent $\lambda_{C,\text{h,\lx@glossaries@gls@link{acronym}{ew}{\leavevmode EW}}}/z_{\text{\lx@glossaries@gls@link{acronym}{ew}{\leavevmode EW}}}=1$, where $z_{\text{\lx@glossaries@gls@link{acronym}{ew}{\leavevmode EW}}}=2$. In Sect. 3.2.2 we shall reproduce this result numerically as a test of our simulations.

3 Results

Extensive dynamical simulations were performed using both RS and SCA updating schemes. To avoid finite-size effects we considered large systems with lateral sizes of $L=2^{16}$.

In SCA simulations the deposition probability must be $p<1$ in order to allow stochastic noise. We investigated three cases: $p=\num{.5},\num{.75}$ and $\num{.95}$ in depth. While KPZ runs were performed without removals: $q=0$, in the EW growth we applied $p=q=1$ for RS and $p=q=\num{.5}$ for SCA.

The roughness scaling of the interface width is analyzed in section 3.1. This is followed by autocorrelation and aging studies of the height as well as lattice-gas variables in section 3.2.

3.1 Roughness scaling

To compare numerical results coming from different updates we determined experimentally a scaling function $f(t,p)$, that provides collapses of $W(L,f(t,p))$ for different dynamics. In case or RS dynamics this function is linear $f_{\mathrm{RS}}(t,p)\propto p$. Since the $p\to 0$ limit of the SCA corresponds to RS updating, we tried to extend the linear form analytically. A smaller survey study of SCA for a larger number of different $p\leq\num{.95}$ values was used to obtain this function numerically and resulted in the following nonlinear extension:

[TABLE]

The speedup with respect to linear function of RS can be understood as follows. A dimer, that was moved at a given time step becomes the target of another update at the next sub-lattice step in the $p\to 1$, $q=0$ case. This is more effective than random sequential updating. Therefore, the roughness growth is faster under SCA than under RS dynamics. One can test this function by observing a reasonably collapse on figure 2(a) for different $p$-s of SCA as compared to the RS results.

Figure 2(b) shows the effective scaling exponents $\beta_{\mathrm{eff}}$, as defined in $\eqref{eq:beff}$, for SCA and RS simulations as the function of the rescaled time variable. Most notably, the $\beta_{\mathrm{eff}}$ exponent results exhibit slightly shifted plateaus for almost two decades in time, but the difference lies well within the error-margin of our best published result $\beta_{\mathrm{eff}}=\num{.2415}(15)$ [50].

Random-sequential

The pronounced plateau visible for DTrDB, suggests that corrections at these late times are small, thus the $\beta_{\mathrm{eff}}$ here should be close to the asymptotic value for $\beta$. This leads to the estimate $\beta=\num{.2414}(2)$, where the error margin is about the size of the $1\sigma$-error bars attached to the effective exponents at late times.

Stochastic cellular automaton

Like in the RS case, there are almost two decades long plateaus in the effective exponents, depending on $p$ for $L=2^{16}$. The plateau value differences are beyond the statistical fluctuations and the DTrDB result. The deviation from the RS result shrinks as we decrease $p$, i.e. as we introduce more and more randomness. This is plausible, but smaller $p$ also means less effective simulations.

The plots also show a break down of $\beta_{\mathrm{eff}}$ at late times. Such behavior can be attributed to the onset of the steady state, which is not apparent from existing finite size scaling studies [50], where $\xi\propto t^{z}\sim L$ appears to be reached about one decade later than the left end of the displayed plot.

Most importantly, the $\beta_{\mathrm{eff}}$ curve does not show this cut-off in the plateau in case of our largest sized $L=2^{17}$ data, but matches perfectly the RS result. It only shows noise related oscillations within the $1\sigma$-error margin. This indicates that the cut-off is related to finite sizes that will be investigated further in the following section.

3.1.1 Distribution of interface heights in the growth regime

In order to get information about the shapes of the distribution of the SCA interface heights we calculated their lowest moments for $L=2^{16}$ and for a smaller $L=2^{12}$, to distinguish finite time from finite size-corrections. Figure 3 shows the evolution of the cumulant ratios $S[h]$ and $Q[h]$, defined by (7) and (8). The curves approach their growth regime asymptotic values, but move away again at late times. These values: $S_{\infty}$ and $Q_{\infty}$ can be determined by performing a fit of the form:

[TABLE]

where $\beta$ is the growth exponent, which is motivated by the KPZ ansatz discussed in the next section. We use $\widetilde{t}$, defined in (23), so the timescales match between RS and SCA run across various deposition probabilities. In (34), $R$ is a placeholder for $S[h]$ or $Q[h]$, in the interval: $200\leq\widetilde{t}\leq\SI{200000}{MCS}$, which excludes early time oscillations as well as the cut-off at late times, coming from $\xi\to L$. This yields $S_{\infty}=\num{-0.427}(2)$ and $Q_{\infty}=\num{.352}(3)$ for the growth regime, in agreement with literature values [51, 52, 53, 54] of the KPZ universality class. The sign of $S$ depends on the choice $p\gtrless q$ in the simulations, corresponding to $\operatorname{sign}(\lambda)$ of the KPZ equation (1).

Panel (b) of figure 3 shows the deviations from these asymptotic values. The error estimates given above originate from this representation: The error is assumed to be on the order of the closest approach of the numerical data to the asymptotic value.

After the closest approach to the asymptotic values in the growth regime, $S(\widetilde{t})$ and $Q(\widetilde{t})$, both, move in the direction to their respective values in the steady state: $S_{\infty}\approx\num{.26}$ and $Q_{\infty}\approx\num{.13}$ [42, 55, 56, 57]. The shape of the distribution of surface heights changing in this way is an indication of finite-size effects becoming relevant at $\widetilde{t}_{\mathrm{fs}}\approx\SI{3e5}{MCS}$. This coincides with the time at which the cut-off 111 The cut-off was observed at $\widetilde{t}^{\prime}\approx\SI{1.7e5}{MCS}$ at two different times contributing in the calculation of $\beta_{\mathrm{eff}}(\widetilde{t}^{\prime})$: $\widetilde{t}_{1}\approx\SI{2e4}{MCS}$ and $\widetilde{t}_{2}\approx\SI{3.2e5}{MCS}\gtrsim t_{\mathrm{fs}}$. in $\beta_{\mathrm{eff}}(\widetilde{t})$ was observed in SCA runs for $L=2^{16}$ (see figure 2(b)). Hence it becomes clear, that this change in $\beta_{\mathrm{eff}}$ is caused by finite-size effects. In figure 3(b) we also plotted $S[h](\widetilde{t})$ of a smaller system ($L=2^{12}$), for which the steady state is reached after the relaxation time $\tau=L^{z}\approx\SI{6e5}{MCS}$, but finite-size corrections are evidently relevant long before then.

3.1.2 KPZ ansatz for the growth regime

Analytical and numerical investigations of KPZ models in $1+1$ dimensions found that finite-time corrections to $h(t)$ took the form $\propto t^{-\beta}$ for the interface height [58, 24, 59, 60, 61]:

[TABLE]

with non-universal parameters $c_{n}$ and $N$ 222 This assumes that $\xi$ and $\zeta$ are independent, but there is no guaranty for that. For ballistic deposition in $D=1$, a strange correction exponent, close to $1/2$ was observed (see e.g. [60]), while in our recent RSOS model simulations [57] we also found correction exponent $\beta$ in case of $N>1$ levels.. Moreover, good agreement between the numerics and experiments has been found  [62, 63, 31]. In the $2+1$ dimensional restricted solid-on-solid model (RSOS) model, the dominant corrections to the roughness growth were found to be of order $~{}\propto t^{-4\beta}$ [54], which motivates the inclusion of higher orders in these forms in higher dimensions [51, 54, 52, 64]. Ideally, such a model would fit the data well as soon as all relevant orders are included. Adding more terms should not improve the fit quality further. However, adding more free parameters in this way can result in overfitting of a noisy data, if not convergence-problems.

Figure 4 shows fitting results using (33) on the previously introduced datasets. It is immediately apparent, that (33) with $N=1$ does not describe the presented data, $n=2$-terms are required, as in case of the RSOS model. In case of RS simulations, the ansatz appears to fit reasonably well early times: $\widetilde{t}\geq 100$ as well. The SCA runs on the other hand show strong oscillations at early times, caused by the synchronous updates, and are not described well by the KPZ ansatz here. Still, late times before the finite size cutoff becomes effective, (in the interval $\num{1e-5}\leq 1/\widetilde{t}\leq\num{3e-3}$) can be fitted well by (33), suggesting universality of the corrections. This becomes true in the $p\to 0$ limit, as in case of figure 2.

The spread of $\beta$ values for larger $N$ provides an estimate for overfitting and may serve as an error estimate for a small confidence interval of $1\sigma$. For simulations with RS dynamics, this yields $\beta=\num{0.2414}(2)$, which, remarkably, is identical to the result based on the average of the late-time plateaus of $\beta_{\mathrm{eff}}$. Fit parameters for $N\leq 6$ can be found in table 3 of the supplemental material.

3.1.3 The KPZ dynamical exponent

The dynamical exponent $z=\alpha/\beta$ of the KPZ class is related to the roughness exponent $\alpha$ by the Galilean symmetry [12]:

[TABLE]

Inserting the our estimate $\beta=\num{.2414}(2)$ into this equation yields $\alpha=\num{.3889}(3)$ and $z=\num{1.611}(3)$. The latter is used to calculate autocorrelation exponents in the next section. It should be noted, that the above value for $\alpha$, while in agreement with earlier numerical estimates [42, 50, 65, 64], marginally disagrees with the currently most accepted one $\alpha=\num{.3869}(4)$ [66]. Combining this roughness exponent results with our own estimate for $\beta$ violates equation (26) by about $\num{2.5}\sigma$. A slight violation of the Galilean invariance, which was proposed for discrete systems [67], may explain this disagreement. If this is the case the correct dynamical exponent would be $z=\num{1.603}(3)$. However, the validity of the Galilean invariance is still widely accepted in literature [68, 65, 18], for this reason we use our $\alpha$ and $z$ estimates, obtained using  (26), for consistency.

3.2 Autocorrelation

3.2.1 Autocorrelation of interface heights in the KPZ

case

Aging

The autocorrelation results of the interface heights under RS dynamics are summarized in figure 5. A near-perfect collapse of the $C_{h}(t,s)$ functions could be achieved by using $b$ from the relation (14).

Autocorrelation exponent obtained by RS dynamics

We calculated effective exponents of $\lambda_{C,h}/z$ and its $t/s\to\infty$ behavior, by an analysis shown in the right inset of figure 5. In order to read-off the appropriate correction to scaling we linearised the left tail of the curves by plotting them on the $\sqrt{s/t}$ scale. This leads to the extrapolations $\lambda_{C,\mathrm{h}}/z=\num{1.254}(9)$ depending on $s$.

To clarify the situation we attempted a different type of local slope analysis presented in the left inset of figure 5, using tail effective exponents, where each $\lambda_{C,\mathrm{h}}/z$ value was determined as the exponent of a PL-fit to $C_{h}(t^{\prime},s)$ for $t^{\prime}\geq t$. These can be expected to converge more monotonically to the asymptotic value as before, because the left tail data of $C_{h}(t,s)$ are included in the procedure for all $t_{\mathrm{min}}$ with an increasing weight as $t_{\mathrm{min}}$ increases. Indeed, the curves of different $s$ values in figure 5(b) behave more linearly with some additional oscillations. However, all curves seem to fluctuate around a common mean, which is not the case for the local slope analysis. A single linear fit for the combination of all curves, yields an averaged extrapolation of $\widetilde{\lambda}_{C,\mathrm{h}}/z=\num{1.23}(3)$ in a marginal agreement with our previous result for this $\lambda_{C,h}/z=1.21(1)$ [30] and with the value obtained in [64, 31] for intermediate times.

The present larger error margin takes the uncertainty due to the actually unknown corrections into account. This problem is illustrated in the comparison between effective exponents and tail effective exponent, where the former show a smaller apparent extrapolation error. In the following we use the simpler extrapolation method based on effective exponents, but estimate the error from their direct fluctuations rather than the uncertainty of the extrapolation fit.

Using our $z$ value the corresponding autocorrelation exponent is $\lambda_{C,\mathrm{heights}}=\num{1.98}(5)$. These results also hold in simulations with more coarse DD, where one would expect an observable difference if any artificial correlations were present.

SCA autocorrelation functions and aging

SCA updates are spatially correlated, therefore they introduce a contribution to the autocorrelation function, which depends on the update probability $p<1$. If we want to model cellular automaton like systems this is not a problem, but for describing the KPZ equation this is artificial. Figure 6 compares the autocorrelation functions of height variables at $p=\num{.95}$ and $p=\num{.5}$. The most apparent property is the finite asymptotic value (figure 6(a)). This is the consequence of frozen regions, arising in ordered domains, which are difficult to randomize by the SCA dynamics. In the dimer model updates can happen at the boundaries only, besides this alternating domains are also stable in case of SCA, they flip-flop at even-odd sublattice steps, when $p\to 1$.

We applied an iterative fitting procedure to determine the functional behavior as follows. As a first approximation the $C_{\mathrm{h}}(t\to\infty,s,p)=o(p)$ limit was determined using a linear extrapolation from the function’s right tail. Subtracting the appropriate value from each curve revealed a PL approach to this constant. To obtain refined $o(p)$ values, the exponent $x$ was read off from the data, allowing a subsequent fit for the tail in the form:

[TABLE]

with free parameters $o$ and $c$. The corrected exponents converged as $x^{\prime}\to\lambda_{C}/z$, after subtracting the refined $o(p)$ values. These iterations yielded self-consistent estimates for $o(p)$ and the autocorrelation exponent of the SCA. This procedure is more prone to statistical error for small $t/s$, because $C_{\mathrm{h}}(t,s)$ is farther away from the asymptotic behavior in this case, allowing noise in the tail to influence the extrapolated value more strongly. Table 1 lists the calculated $o(p)$ limits (including those for the lattice-gas variables, see Section 3.2.3).

The limiting value turned out to depend exponentially on $p$. Note, that similar $e^{p}$ dependence has been found in $f_{\mathrm{SCA}}(t,s)$ relating SCA and RS timescales 333 These limits could also be determined from the small survey study presented in figure 2(a), comprising much smaller sample sizes than the results presented in detail in the following. This data suggests an exponential dependence $o(p)\propto\exp(\nu p)$ with a similar, or possibly the same, value for the parameter $\nu$ for both slopes and heights. However, these autocorrelation measurements used the same waiting time $s$, without taking into account the $p$-dependent time-scale. Thus the actual waiting times $\widetilde{s}$ decrease with $p$, which makes the fit performed on the $o(p)$ across these runs unsuitable to determine a reliable value for $\nu$. .

Figure 6(b) shows the corrected $C_{\mathrm{h}}(t,s)$ functions, after subtracting the limiting $o(p)$ values. A nearly perfect data collapse could be achieved using the aging exponent $b_{\mathrm{h}}$, coming from the RS simulations. Even more, the corrected SCA and the displayed RS autocorrelation functions show identical behavior.

Autocorrelation exponent: SCA

Local slope analyses of the corrected autocorrelation functions are displayed in figures 6(c) and 6(d) at $p=\num{.95}$ and $p=\num{.5}$, respectively. Assuming a rescaling of the abscissa: $\sqrt{s/t}$, allows one to observe a linear behavior of the effective exponents for intermediate times. In case of $p=\num{.95}$ the $o$-values could not be determined precisely enough for $s>100$, thus we considered extrapolations at $s=30,100$ only in a weighted average of the results. This yielded: $\lambda_{C,\mathrm{h}}/z=\num{1.26}(1)$ and so $\lambda_{C,\mathrm{h}}=\num{2.01}(2)$. These values are in good agreement with those obtained from a local slope analysis of RS calculations for small $s$.

The effective exponents for $p=\num{.5}$ show a slightly decreasing tendency with $s$ in figure 6(d), moving towards the RS estimate $\widetilde{\lambda}_{C,\mathrm{h}}/z=\num{1.23}(3)$. However, we can’t consider the extrapolated values for $s=500$ and $s=1000$ more precise, than those at $s=30,100$, because the determination of the $o(s)$ constant becomes more uncertain at higher times, increasing the possible error of the exponent estimates.

3.2.2 Autocorrelation of interface heights in the EW

case

Since the autocorrelation function in the EW case is known exactly (22), we can verify our simulations by a comparison with it. Indeed, the expected form could be reproduced by our RS implementation. A more interesting result is, that the SCA simulations also fit it perfectly. The finite saturation value, caused by correlated updates, observed in the KPZ case is not present here.

The agreement with the analytical form is exemplified in the inset of figure 7. A small deviation at very early times can be observed here, as well as in the RS results and should be related to the initial conditions of the simulation with respect to those of the analytical calculations. The application of a fit with  (22) results in $c_{0}\simeq 0.152$ for different waiting times $s$. Using the consistency relation (14) for $s\to t$ we can expect the same value, which was derived for the octahedron model in [38] for the $s\to\infty$ limit.

These numerical results do not only show the correctness of the SCA and RS implementations of the roughening kinetics, but provide an example, where the correlations introduced by SCA do not affect the dynamical behavior.

3.2.3 Autocorrelation of lattice-gas variables in the KPZ case

Next we show results for the lattice-gas variables corresponding to the binary slope values of heights of the KPZ growth presented earlier (see figure 8) using RS dynamics. Here again, the $C_{s}(t,s)$ functions of different waiting times collapse almost perfectly with the value: $b_{\mathrm{s}}=\num{0.76}(2)$. In a previous paper [30] we reported: $b_{\mathrm{s}}=0.70(1)$, which were obtained by a smaller sized analysis.

Autocorrelation exponent: RS

Since the density autocorrelation functions decay much more rapidly than those of the heights, the signal-to-noise ratio in the present sample is insufficient for a reliable extrapolation based on the effective exponents. A weighted average of direct PL fits for $4\leq t/s\leq 90$ yielded $\lambda_{C,\mathrm{s}}/z=\num{2.312}(2)$. However, the effective exponents show curvature as $t/s\to\infty$ and suggest an asymptotic value $\lambda_{C,\mathrm{s,eff}}/z=\num{2.39}(2)$. In [30] we obtained $\lambda_{C,\mathrm{s,eff}}/z=\num{2.35}(2)$, coming from $s=30$, $L=2^{13}$ sized CPU simulations.

SCA density autocorrelation functions

Similarly to the case of interface heights the $C_{\mathrm{s}}(t,s)$ functions approach finite values asymptotically, as shown in figure 9(a) as the consequence of the SCA dynamics. The computed values of $o(p,s)$ are listed in table 1.

Figure 9(b) shows the corrected functions $C_{\mathrm{s}}(t,s)-o(p)$ in comparison with our RS result. Data collapse for $p=0.5$ and $p=0.95$ could be achieved using the common aging exponent value $b_{\mathrm{s}}=0.76$, obtained from our previous RS calculations. This indicates, that the density correlation behavior is not changed by the application of SCA updates as in case of the height variables.

SCA density autocorrelation exponent:

In contrast with the interface height results the density correlation exponent of the corrected $C_{\mathrm{s}}(t,s)$ exhibits a more complex behavior. The dataset for $p=\num{.95}$ clearly exhibits a different exponent than what we observed in case of the RS simulations. We show effective exponents fitting for $s=30$ and $100$, where the signal-to-noise ratio is better on figure 9(c). A direct linear fit extrapolating to $s/t\to 0$ yields an estimate of $\lambda_{C,\mathrm{s,\lx@glossaries@gls@link{acronym}{sca}{\leavevmode SCA}}}/z=\num{.75(2)}$.

At $p=\num{.5}$ we can find a crossover form the RS to a different, SCA asymptotic behavior in figure 9(d). A linear extrapolation for the tail of this crossover curve results in $\lambda_{C,\mathrm{s,\lx@glossaries@gls@link{acronym}{sca}{\leavevmode SCA}}}/z$, in good agreement with the $p=\num{.95}$. This leads to the following numerical form for the tail of the autocorrelation function under SCA dynamics:

[TABLE]

3.2.4 Autocorrelation of lattice-gas variables in the EW

case

RS autocorrelation functions

In case of RS simulations, the tail of $C_{\mathrm{s}}(t,s)$ does not decay with a simple PL as can be observed in the figure 10. The pronounced curvature in the log-log plot suggests a slower than PL decay at first glance. However, the effective exponents (inset) suggest a PL with an asymptotic exponent $\lambda_{C,\mathrm{s}}^{\text{\lx@glossaries@gls@link{acronym}{ew}{\leavevmode EW}}}/z_{\text{\lx@glossaries@gls@link{acronym}{ew}{\leavevmode EW}}}=\num{0.7}(2)$, following a cross-over from an early-time regime.

SCA autocorrelation functions

While the SCA dynamics seems to reproduce the expected autocorrelation function of the surface heights after the removal of the constant, the evolution of the underlying lattice gas is different. In case of $p=q=\num{.5}$ the density autocorrelation exhibits a PL tail, characterized by $\lambda_{C,\mathrm{s}}^{\text{\lx@glossaries@gls@link{acronym}{sca}{\leavevmode SCA}},\num{.5}}/z_{\text{\lx@glossaries@gls@link{acronym}{ew}{\leavevmode EW}}}\approx 2$ (see inset of figure 10).

In the $p\to 0$ , $q\to 0$ limit the SCA crosses over to an effective RS dynamics, because we avoid the correlated updates of the lattices. This is indeed the case here, evidenced by the $p=q=1/32$ results (see figure 10). Following a rescale of time $t=p\cdot t$ one can find a good collapse with the RS results. Therefore, the update dynamics seems to affect the scaling behavior of the density autocorrelation function.

Aging

The aging exponent obtained from the presented simulations is $b_{\mathrm{s}}^{\text{\lx@glossaries@gls@link{acronym}{ew}{\leavevmode EW}}}=\num{1.1}(2)$. This value holds for both RS and SCA dynamics, but breaks down for very small values of $s$.

4 Discussion and conclusions

We performed extensive simulations of the octahedron model by random-sequential (RS) and stochastic cellular automaton (SCA) dynamics. Precise estimates were obtained for the dynamical behavior: exponents as well as probability distributions of the Kardar–Parisi–Zhang (KPZ) and Edwards–Wilkinson (EW) universality classes. The main advance of this work, in the long story of KPZ research, is the influence of correlated SCA dynamics on the universal properties of these models. Furthermore, we determined the aging properties of the underlying driven lattice gasess (DLGs) of the octahedron model.

By determining moments of the probability distributions we could study finite size effects and arrived at the conclusion that the corrections related to this become relevant much before the occurrence of the steady state. Our surface growth simulations support the validity of the KPZ ansatz hypothesis in $(2+1)D$ and yield a growth exponent $\beta=\num{.2414}(2)$, from which $\alpha=\num{.3889}(3)$ and $z=\num{1.611}(3)$ can be deduced. The growth exponent value lies within the error margins of [57, 52, 50], but not within those of the early landmark result of Forrest and Tang $\beta=\num{.240}(1)$ [69]. However, the small simulation cells used then demanded shorter simulations times, which could have lead to a smaller estimate due finite-time corrections. Under SCA dynamics marginally lower growth exponents were observed for deposition probabilities $p>\num{.5}$ and additional corrections to scaling caused the failure of the KPZ ansatz at early times.

Our estimate of the roughness exponent does not agree with the direct estimate $\alpha=\num{.3869}(4)$, obtained recently through a finite-size scaling analysis of the restricted solid-on-solid model (RSOS) model by Pagnani and Parisi [66], which was based on SCA simulations with $p=\num{0.5}$. Numerical differences between SCA and RS dynamics might be a cause of this. However, since our estimate was derived using (26), a slight violation of the Galilean invariance, which was proposed for discrete systems [67], may also explain this disagreement.

Both our RS and SCA simulations reproduced the expected autocorrelation behavior of interface heights in the EW universality class. In the KPZ case correlated updates resulted in $C_{\mathrm{h}}(t,s)$ to approach a finite value asymptotically. However, after the subtraction of this constant we found the same universal power law (PL) tails for both types of site-selection dynamics.

In case of the underlying lattice-gas variables, we found the relevance of the SCA dynamics for the asymptotic autocorrelation decay exponents, but the aging exponent seems to be insensitive for this. Interestingly, in case of the non-linear (KPZ) model the SCA dynamics slows the decay of the autocorrelations, while in the linear (EW) model this results in a shorter memory of the dimer model. This is the consequence of the effectivity of the ordered SCA updates, which enhances the build up (KPZ) or distortion (EW) of homogeneous areas, correlated for long times.

Our estimates for the autocorrelation exponents of the KPZ class are summarized in table 2. We provided numerical results for $C(t,s)$ in the KPZ case with unprecedented accuracy, drawn from timescales up to $t/s=1000$ due to the high signal-to-noise ratios we could achieve by these parallel algorithms implemented on GPUs. These simulations can be help to test predictions of theories like local scale-invariance with logarithmic corrections [34].

The KPZ autocorrelation exponent in $(1+1)$ dimensions was derived analytically $\lambda^{\mathrm{1d}}_{C,\mathrm{h}}=1$ [26, 27]. Later Kallabis and Krug conjectured, that in higher dimensions $\lambda_{C,\mathrm{h}}=D$ [29] applies, but rigourous proof is still missing. Our estimates for $\lambda_{C,\mathrm{h}}$, summarized in table 2, support this hypothesis within error margin both for RS and SCA dynamics.

We have tested the validity of the relation (16) by Krech with our numerical data. The value for the short time dynamical exponent: $\theta=\lambda_{C}/z+2\beta=1.23(2)+0.2414(2)\simeq 1.71(3)$ agrees well with: $(D+4)/z-2=6/1.611(3)-2\simeq 1.72(1)$, therefore we can support the validity of the relation obtained by a perturbative RG analysis [26, 27].

A possible continuation of this work could be the study of the height correlations in the momentum space:

[TABLE]

In particular one should be able to test, whether $C_{q}(t,s)$ decays in an exponential or in a stretched exponential way as predicted in references [70, 71, 72, 73]. Note, that we have already succesfully used the extension of the dimer model to determine the power spectrum density: $S(k,t)=\left\langle h_{-q}(t)h_{q}(t)\right\rangle$ in case of Kuramoto-Sivashinsky type of systems [74].

We can also compare the present estimates with our recently published values for the autoresponse $\lambda_{R}=\num{2.00}(6)$ and the corresponding aging exponent $a=0.24(2)$ [33]. $\lambda_{C}\simeq\lambda_{R}$ seems to hold within error margins. In $(1+1)$ dimensions an exceptional fluctuation-dissipation relation (FDR) exists [75, 12]:

[TABLE]

This implies the exponent relations $\lambda_{C}=\lambda_{R}$ and

[TABLE]

confirmed by simulations [76]. Our $(2+1)D$ results support the first one, but the latter is not satisfied by our numerics:

[TABLE]

This calls for the existence of a possible FDR in higher dimensions. For example the genaralized form

[TABLE]

is satisfied by the exponents within error limits in $D=1,2$ both. Confirmation of this assumption should be a target of further research. An intermediate step in this direction could also occur as an inequality, like one found in the KPZ steady state [77, 78].

The autocorrelation and aging exponents which we found for the driven lattice gas of slopes differ from another two dimensional extension of the totally asymmetric exclusion process (TASEP) described in [79], where $\lambda_{C,\mathrm{s}}/z=1$ and $b_{\mathrm{s}}=1$ are reported.

Finally we point out that the SCA simulations are more efficient because they they allow for optimal memory access patters in contrast to the random accesses required for the RS ones. Technical details of our implementations are published elsewhere [43, 44]. The extension of these algorithms for other surface models, like those with conservation laws [74, 80] or in higher dimensions [81] is straightforward. However, the efficiency of RS implementations, using the approach employed here, decreases with the number of dimensions due to the volume of local cells increasing. SCA simulations do not suffer from this problem and are thus more suitable for higher dimensional problems.

The code used in this work can be found at https://github.com/jkelling/CudaKpz.

Acknowledgments

We are grateful for the useful comments from Malte Henkel and Timothy Halpin-Healy and thank Herbert Spohn, Giorgio Parisi and Uwe Täuber for helpful discussions. Support from the Hungarian research fund OTKA (Grant No. K109577), the Initiative and Networking Fund of the Helmholtz Association via the W2/W3 Programme (W2/W3-026) and the International Helmholtz Research School NanoNet (VH-KO-606) is acknowledged. We gratefully acknowledge computational resources provided by the HZDR computing center, NIIF Hungary and the Center for Information Services and High Performance Computing (ZIH) at TU Dresden.

Appendix Supplement A Octahedron Model

Figure 11 shows a depiction of the update process in the octahedron model after [38].

Appendix Supplement B KPZ ansatz for the growth regime

In the main manuscript we test how the KPZ ansatz [58, 24, 60, 59, 51, 52, 64] hypothesis describes the late-time corrections to the growth exponent $\beta$. Its most general form gives

[TABLE]

A list of best fit parameters for various versions with $N\leq 6$ is provided in table 3, which also lists the reduced sums of residuals $\chi_{\mathrm{red}}$ to quantify agreement between the fitting model and the data.

Where the KPZ ansatz does indeed apply, fits of the more general models (32) should not show increased agreement with the data. The table shows them to be less consistent with respect to the resulting estimates for $\beta$. They provide the best description of the data with only the term $\propto\widetilde{t}^{-2\beta}$, but one or two additional odd terms present ($N=2,3$). The best fits resulting from models (33) are consistently better than those of (32), across all datasets, which justifies discarding the latter class of models and thereby supports the KPZ ansatz hypothesis for the roughness growth.

Appendix Supplement C Distribution of interface heights in the growth regime

In the main manuscript, the finite time corrections to the cumulant ratio of the distribution of surface heights are shown to be described by the KPZ ansatz, yielding the form:

[TABLE]

a fit can also be made, resulting in the exponents $c_{S}\approx\num{-.54}$ and $c_{Q}\approx\num{-.60}$ for $S[h](t)$ and $Q[h](t)$, respectively. The value $c_{R}\approx\num{-.6}$ also holds for a number of other dimensionless cumulant rations, not displayed here. The values for $S_{\infty}$ and $Q_{\infty}$ resulting from this extrapolation stay well within the margins of error given in the main manuscript.