{'fraction_non_alphanumeric': 0.04463584069405589, 'fraction_numerical': 0.02087197057806557, 'mean_word_length': 4.396776929601357, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large scale spatio-temporal behaviour in surface growth

This paper presents new findings concerning the dynamics of the slow height variations in surfaces produced by the two-dimensional isotropic Kuramoto-Sivashinsky equation with an additional nonlinear term. In addition to the disordered patterns of specific size evident at small scales, slow height variations of scale-free character become increasingly evident when the system size is increased. The surface spectrum at small wave numbers has a power-law shape with a lower cut-off due to the finite system size. The temporal properties of these long-range height variations are investigated by analysing the time series of surface roughness fluctuations. The resulting power-spectral densities can be expressed as a sum of white noise and a generalized Lorentzian whose cut-off frequency varies with system size. The dependence of this lower cut-off frequency on the smallest wave number connects spatial and temporal properties and gives new insight into the surface evolution on large scales.


INTRODUCTION
Detailed understanding of surface growth physics has fueled several advances in science and technology, including more accurate dating in archeology [1], better integrated circuit technology [2], as well as production of novel materials [3]. In many cases of scientific and technological interest, the evolution of growing surfaces can be described by so-called continuum models that consist of nonlinear partial differential equations and often display rich and interesting dynamics [4]. Even though not all of this dynamics is currently accessible experimentally, it still is worthwhile to investigate, especially in the view of rapid experimental [5] and theoretical [6] progress.
The object of this study is a continuum surface growth model described by the two-dimensional generalized Kuramoto-Sivashisky equation with a single independent parameter α, considered in [7], that produces chaotically evolving disordered spatial patterns. Equations of this type (with and without added noise) have been successfully used as models for amorphous solid surface growth [8][9][10] and nano-scale pattern formation induced by ion beam sputtering (IBS) [11][12][13][14][15]. Eq. (1) in two spatial dimensions describes the evolution of a (2+1)-dimensional interface, i.e., a surface whose height h(r, t) is defined as a function on a twodimensional plane r ∈ R 2 that is growing in the direction h perpendicular to that plane as time t goes by. Numerical studies of Eq. (1) in one dimension have also been performed by Muñoz-Garcia et. al. [16], and a good correspondence to the IBS experiments has been found [17].
Eq. (1) has the celebrated Kuramoto-Sivashinsky (KS) equation [18][19][20][21] as its special case when parameter α = 0: The latter equation stands as a paradigmatic model for chaotic spatially extended systems and has been used to study the connections between chaotic dynamics at small scales and apparent stochastic behaviour at large scales [22][23][24]. Various generalizations and modifications of the KS equation (2) with local and non-local damping terms, anisotropy, and noise have been used to study pattern formation due to ion-beam erosion [25][26][27][28][29][30]. Eq. (2) itself in one-and two-dimensional cases has been a subject of active research for about three decades, and its scaling properties have even been an object of controversy. It has been suggested by Yakhot [31] and subsequently confirmed and reiterated by different authors (see e.g., [22,23,32]) that the large-scale behaviour of the deterministic KS equation (2) in the one-dimensional case can be described by a stochastic equation where η represents random uncorrelated Gaussian noise. Equation (3) has become known as the Kardar-Parisi-Zhang (KPZ) equation [33]. It was originally proposed as a continuum model for surface growth due to ballistic deposition [34], since it showed the same dynamic scaling behaviour [35]. However, the correspondence between KS and KPZ in two-dimensions has led to disagreements by the same authors [23,32,36,37], because of the lack of conclusive analytical results. More recent results [24,38] tend to support the conjecture that KPZ and KS equations belong to the same universality class, although the numerical results for the deterministic (noiseless) KS are not conclusive due to the extremely long transient effects. Recent numerical results for the two-dimensional KS [7] with much longer simulation times show the same scaling properties of the saturated surface roughness as obtained by Manneville and Chaté [39] for the two-dimentional KPZ, thus, further supporting the argument that the two-dimensional KS and KPZ equations belong to the same universality class. However, [7] has found different scaling behaviour in the less researched generalized KS case (1) with α > 0.
The purpose of this paper is to demonstrate in a wider parameter range the validity of the scaling relations for the surface roughness reported in the previous work [7], and to investigate the dynamics of the scale-free lowwavenumber spatial variations that these relations imply.
The article is structured as follows. After commenting on numerical methods and parameters used in the investigation, Sec. I demonstrates the initial transient kinetics of the surface roughness and eventual cross-over to the saturation regime observed in the investigated parameter range. Sec. II presents morphologies of surfaces produced by the generalized KS equation (1) in the saturated regime and presents results on the finite-size scaling of the saturated surface roughness. The analysis of temporal behaviour of the surface roughness for the KS case is demonstrated in greater detail in Sec. III, and the results of the same analysis for the generalized-KS case are then given in Sec. IV. We conclude with a summary of our findings and propose several directions for future work in Sec. V.  (1) with several values of parameter α. Each simulation starts from a random uncorrelated surface profile with initial roughness w(0) = 10 −4 . The left panel shows w(t) on linear scale for 3 realizations with each α. The right panel highlights the initial exponential increase of w(t) by using a semi-logarithmic scale. The inset demonstrates the saturation roughness wsat and the cross-over time t×.

I. SURFACE ROUGHNESS: TRANSIENT KINETICS AND SATURATION
In this investigation, the equation (1) is solved numerically for different values of α using the finite difference method with periodic boundary conditions, the time step ∆t = 0.005, and spatial discretization step ∆x = 0.71086127010534 . Such a seemingly bizarre number for the discretization step ∆x is actually a good approximation of the value that is needed in order for the system with periodic boundary conditions to be able to contain hexagonal patterns that appear in some other versions of the generalized KS equation (see for example [25,30]). The equation is solved for system sizes L ranging from about 36 to about 711 (i.e., on the N × N lattices with N from 50 to 1000, where L = N ∆x). Different methods of numerical solution for (1) are presented and compared in [40].
In the KS case (α = 0), the evolving surfaces reach the regime where the dynamics is chaotic, but statistically stationary. This type of behaviour also persists for α > 0, at least up to α = 5. However, for larger values of α, this stationary chaotic behaviour gives way to non-stationary effects that prevent the saturation in the surface evolution. Indeed, in the limiting case when α → ∞ in (1), by rescaling h, one arrives at the conserved Kuramoto-Sivashinsky equation [41], which produces a non-stationary regime with ever increasing surface roughness due to the uninterrupted coarsening of the surface patterns. Thus, by increasing α, there must be a route from the stationary chaotic evolution to non-stationary coarsening behaviour. Nonetheless, the long-time behaviour at the intermediate α values seems to be quite complicated and has not been studied in detail so far. Even though negative α values do not follow from surface growth or erosion models (see, e.g., [9,15]), it is desirable to understand the dynamics in that range for completeness. However, for α < −0.14 we find that large local gradients in the surface emerge and grow. They eventually exceed the numerical capacity of the simulation, thus making the required long-time calculations unstable. Therefore, this work focuses on the evolution of surfaces produced by (1) in a moderate range, −0.12 ≤ α ≤ 5, of parameter values where the long time behaviour is stationary.
One of the most important quantities characterizing a surface [42,43] is the surface roughness w(t), also called the surface width: The scaling properties of this quantity are often used to characterize and classify various surface growth models into various universality classes [22,35,38,43]. The roughness of an evolving surface changes with time. In the range of parameter values considered here, the kinetics of w(t) due to the surface evolution according to (1) seems to follow a distinct pattern (see, e.g., [7,9,10,41]): starting from a random surface with some small initial roughness w(t = 0) 1, the roughness begins to grow at an exponential rate, but at some time t × 100 this growth slows down and, later on, crosses over to a stationary regime where it oscillates about some average (saturation) value w sat . This transient behaviour is shown in Fig. 1) for several parameter values.
The value of saturated surface roughness can be defined as follows: Here t 0 t × is a time at which all initial transient effects have decayed and are virtually undetectable, i.e., the time at which the stationary regime has been reached. In practice, the total observation time T has to be much larger than the typical time scale in the kinetics of w(t).
In the investigation presented here, the saturation values for surface roughness w sat are calculated using t 0 = 2·10 4 and T = 8 · 10 4 . Note that these times are significantly larger than those recently achieved by Muñoz-Garcia et. al. in the numerical investigation of an equation equivalent to (1) in the one-dimensional case [16,17]. There, although the 'interrupted coarsening' is observed, the saturated stationary regime appears not to have been fully reached.
The surface roughness w(t) represents the integral effect of all modes contributing to the surface morphology. Therefore, in this work, the time series of chaotic fluctuations of w(t) in the stationary regime are used to investigate the long-time dynamics of surfaces, in particular, the temporal behaviour of the large-scale height variations observed in Ref. [7].

II. SURFACE MORPHOLOGIES AND THE SCALING OF ROUGHNESS
The surface profiles produced by (1) in the stationary regime have a disordered cellular structure [7,9] (c.f. Fig. 2 and Fig. 3). Since (1) is isotropic, and, consequently, the resulting profiles have no distinct direction on the r-plane, the surface morphologies are investigated by averaging the surface height autocorrelation function over all directions at a distance r = |r|:  (7) changes from monotonically decreasing at α = −0.12 (corresponding to 'flaky' surface profiles, with 'flakes' of widely varying size) to having a short flat region at α = 0 (corresponding to a profile with 'cells' of similar size), and to a function with at least one distinct peak at α > 0 whose distance increases with α (corresponding to the surface 'cells' becoming almost round 'humps' whose size increases with α).
Another thing that can be noticed in Fig. 2 is that the normalized correlation function C(r)/w 2 decays slowly for α = 0 and faster for increasing α. Also, perhaps surprisingly, the autocorrelation function for α < 0 decays faster than for α = 0. These are the first indications of the influence of parameter α on long-range height correlations.
Simulations show that the resulting saturated surface roughness (6) increases with the system size. This indicates that the surface profiles of larger systems contain additional spatial Fourier components of smaller wave number k, since the structure on smaller scales remains virtually unchanged [7].
Large-scale height variations in surfaces produced by (1) become more distinct as the system size is chosen to be many times larger than the typical cell size (see Fig. 3).

A. Scaling of roughness
The isotropic surface power-spectral density (PSD), S(k), defined as the absolute square of the Fourier transform of the surface profile integrated over all directions of the wave vectors k, can be obtained from the isotropic surface autocorrelation function C(r) [7,44]: Here J 0 (kr) is the Bessel function of the 1st kind: An example of numerically calculated surface PSD S(k) using (7)-(9) for α = 1 is shown in Fig. 4. Here, one can see a distinct peak that corresponds to an average size of a hump-shaped cells in the surface pattern (see Figs. 2 and 3) and a power-law trend for small wave numbers.
The integral of the PSD S(k) (8) over all wave numbers k equals the variance of the surface profile which is the square of surface roughness: Since the surfaces in numerical simulations are represented on a discrete (N × N ) lattice of finite size L with a discretization step ∆x, wave numbers that can fit into the system are k n = n∆k with n = 1, . . . , N and ∆k = 2π/L. For large enough systems with N 1, according to (10), the square of the surface roughness can then be expressed as: where If the discretization step ∆x is kept constant (implying k max = const), and the surface patterns at different system sizes L (up to the smallest wave number k min ∝ L −1 ) remain statistically the same (as in Fig. 4), then, by increasing the system size L, the calculated dependence w 2 (L) should yield, according to (12), the shape of the surface PSD S(k) for small wave numbers k → 0. This is useful, since, for larger systems, the direct calculation of the two-dimensional autocorrelation function (7) and surface spectrum (8) can take a very long computation time. In [7], an assumption was made that the PSD S(k) (8) of surfaces produced by (1) has a power-law shape for small wave numbers (below some value k s ): By substituting (12) into (11), one gets three qualitatively distinct scaling behaviours w 2 (L) for L > 2π k −1 s , depending on the value of spectral exponent γ in (12): For asymptotically large systems L → ∞, (13) would become corresponding to asymptotically constant roughness for γ < 1, logarithmically increasing square of the surface roughness for γ = 1, and power-law scaling for γ > 1.
It has been shown in [7] that the assumption (12) of a power-law surface PSD with (0 < γ ≤ 1) at small wave numbers is indeed valid for surfaces produced by (1) with parameter values 0 ≤ α ≤ 1, since the relations (13) fit the numerically calculated surface roughness exceptionally well.
Investigations of a broader parameter range, −0.12 ≤ α ≤ 5, presented in this paper, show that the same assumption (12) also holds for other parameter values. Fig. 5 shows the calculated square of the surface roughness w 2 dependence on the system size L = N ∆x. In order to fit the results with different α values in the same plot, the numerical results and their fits for each α have been divided by the corresponding w 2 values at N = 250. At α = 0 the resulting spectral exponent γ = 1 gives the logarithmic dependence w 2 (N ) (see (13)) which is a straight line in the log-linear scale. This scaling is the same as found by Manneville and Chaté for the twodimensional KPZ equation [39].
As the parameter increases from α = 0 to α = 5, the γ values are found to decrease from γ = 1 to γ ≈ 0.55 (see Fig. 5). This corresponds to slower-than-linear growth of w 2 with ln N . Hence, for large systems w 2 approaches a finite value. Perhaps unexpectedly, for α < 0, the exponent γ has also been observed to become smaller than 1. Therefore, we conclude that the scaling properties of the generalized KS equation (1)

III. ANALYSIS OF ROUGHNESS DYNAMICS IN THE KURAMOTO-SIVASHINSKY CASE
Model equation (1) produces disordered spatial patterns that evolve in time. As shown in Sec. II, with an increase of the system size L new long range height variations appear in the resulting surface profiles in addition to the small scale patterns. The apparent scalefree character of these slow height variations is very different from the cellular patterns on small scales which have a characteristic length (the average size of a 'cell' or 'hump'). Also, the spatial properties of both, the small scale patterns and the large scale height variations, depend strongly on the value of parameter α in (1). This section investigates the corresponding dynamics of these surfaces.
In order to understand the complex spatio-temporal behaviour of (1), we investigate the dynamics of surfaces it produces by analysing the numerically obtained time series of the surface roughness w(t) which contains the collective behaviour of all modes. The time series of w(t) are investigated in the time interval t ∈ [2 · 10 4 , 10 5 ) with sampling time τ sample = 1 (i.e., sampled every 200 time steps ∆t = 0.005), that is, 8 · 10 4 values in total for every realization. The results are averaged over 5 to 10 realizations (differing in the initial surface profile) for every parameter α value. For the range of parameter values explored here, the surface evolution can be considered stationary and ergodic, since the statistical properties of w(t) (average, standard deviation, skewness, autocorrelation function) seem to vary little from realization to realization. Moreover, their values calculated in large enough subintervals of the total time interval differ only slightly from each other.
In this section, the analysis of w(t) is presented in more detail for parameter value α = 0, that is, the Kuramoto-Sivashinsky case (2). The same analysis performed on other parameter values is discussed in Sec. IV. A. Occurrence of slow modes Fig. 6 shows a representative sample of a surface roughness w(t) time series for α = 0 and system sizes varying from N = 125 to N = 1000. Even though the roughness dynamics is dominated by white noise for small systems (N = 125), additional slow modes appear as the system size is increased. For relatively large systems (N = 1000), the time series in question is similar to a signal produced by a random walk. This transition can be visualized even more clearly by using the recurrence plot technique [45][46][47][48] (see Fig. 7). There, a time series s(t) is depicted by plotting a matrix R ti,tj . In the plot, the axes represent the discrete time t i and t j . A black dot (R ti,tj = 1) is put at a point (t i , t j ) if the values of the time series s(t) at these times coincide (recur) to a given accuracy . The pixel remains white otherwise (value R ti,tj = 0), that is: where Θ(x) is the Heaviside step function. Each of the recurrence plots in Fig. 7 is made for a single realization of w(t) in the time interval t ∈ [8·10 4 , 10 5 ), i.e., one fourth of the total length of the time series is investigated. The slow fluctuations of w(t) that appear when the system size is increased can be attributed to the low wave number spatial modes that occur in larger systems. By investigating the scaling properties of these fluctuations, connections between spatial and temporal properties of the corresponding large scale height variations can be made.

B. Autocorrelation functions
The character of the slow fluctuations that appear in the time series (TS) of the surface roughness w(t) resulting from (1) for large systems (see Figs. 6 and 7) is captured by their autocorrelation functions, wherew = w(t) t is the average value of w(t) in the stationary regime. Fig. 8 shows the autocorrelation functions obtained from TS of w(t) with α = 0 for four different system sizes increasing by the factor of 2: N = 125, 250, 500, 1000 . In the top panel of Fig. 8, the normalized (i.e., divided by the variance σ 2 w = (w(t) −w) 2 t ≡ A(0)) autocorrelation functions are displayed in the log-linear scale. In this plot, one can immediately recognize the way in which the characteristic time scales in w(t) grow with N . For instance, by defining some characteristic correlation time τ corr as, for example, the lag τ at which the autocorrelation function decays to the 10% (dashed horizontal line) of its initial value at τ = 0, i.e., one can see that it increases by about the same factor (corresponding to almost constant shifts along a logarithmic scale of τ axis) as the system size N increases by a factor of 2. This indicates that the characteristic time τ corr grows as a power law of N : Further insight into the dynamics can be gained by looking at the same autocorrelation functions in a semilogarithmic plot, as displayed on the bottom panel of Fig. 8. Plotted this way, the autocorrelation functions A(τ ) appear almost as straight lines (with an additional kink at very small τ ) indicating that the their shape should be approximately exponential: C. Power spectra and characteristic frequencies In order to obtain more quantitative results, it is essential to look at the shape of the corresponding power spectra of w(t). As stated by the Wiener-Khinchin theorem [49], the power spectral density (PSD) W (f ) of a signal can be obtained by Fourier transforming its autocorrelation function (16): By substituting the exponentially decaying autocorrelation function A(τ ) ∝ e −λ|τ | (as in (19)) into (20), the PSD W (f ) of a Lorentzian shape is obtained: where f 0 = λ/2π is the characteristic frequency that signifies the cross-over between different behaviours of W (f ), namely: Thus, f 0 represents the lowest frequency (or the lowest decay rate λ ∝ f 0 ) that affects the dynamics of w(t). The above considerations suggest that f 0 must correspond to the lowest wave number, k min ∝ L −1 , of a spatial mode occurring in the system of size L. Plotted in semi-logarithmic scale (bottom panel of Fig. 8), the autocorrelation functions A(τ ) appear as almost straight lines corresponding to the approximate exponential decay (19) whose PSD is a Lorentzian (21). Nevertheless, there are deviations from this trend at very short lag times τ . These deviations correspond to additional fluctuations with a very short correlation timea white noise whose PSD is a constant. Therefore, the resulting PSD of w(t) can be fitted by a Lorentzian plus a constant: where A, B and f 0 are fit parameters. The PSDs obtained from the autocorrelation functions of the surface roughness at α = 0 for different system sizes N = L/∆x are shown in the top panel of Fig. 9 together with their fits by (23). A closer analysis shows that a function with a generalized Lorentzian plus a constant B, with β = 1.8 fits the calculated PSDs even better (see the bottom panel of Fig. 9). The cross-over frequency f 0 obtained as a fit parameter represents the lowest frequency (corresponding to the longest time scale) in the kinetics of w(t). In Fig. 9, it is clearly visible that f 0 decreases as the system size N is increased. Since the lowest wave number k min of the spatial modes occurring in the system is inversely proportional to the system size, k min ∝ N −1 , the f 0 (N ) dependence connects the spatial and the temporal scales. Indeed, by defining the some critical wave number k 0 as one can obtain a dispersion relation f 0 (k 0 ) -a connection between the lowest wave number in the system and its corresponding frequency. The resulting f 0 dependence on k 0 ∆x/(2π) = N −1 for α = 0 is shown in Fig. 10 in the double-logarithmic scale. Plotted this way, the results appear to lie on a straight line, meaning that the relation is approximately a power-law f 0 ∝ k ξ 0 with the exponent ξ ≈ 1.89, as the fit shows (c.f. Fig. 10).

IV. DYNAMICS OF ROUGHNESS FOR OTHER PARAMETER VALUES
The fluctuations of w(t) change character as parameter α is varied. This can already be seen from their time series (Fig. 11). This section presents some of the results on spatio-temporal properties of surfaces evolving according to (1) with parameter values α = 0 in order to point out the similarities and differences from the α = 0 case presented in Sec. III.

A. Spatio-temporal properties
The same type of analysis, as presented in Sec. III for parameter α = 0, has also been performed for other parameter values.
As in the α = 0 case, for α = 0, the occurrence of slow modes can also be observed as the system size increases. However, since the character of low wave number spatial variations depends on α, as shown in Sec. II, their temporal properties also differ.
The PSDs of w(t) for −0.12 ≤ α ≤ 5 can be fitted very well (see Fig. 15) by a generalized Lorentzian with an added constant (24) at different system sizes N (except for some cases discussed in the following subsection). The exponent β in the fit (24) increases monotonically from β ≈ 1.7 for α = −0.12 to β ≈ 3 for α = 5. From these fits at different system sizes N , the relations between the lowest frequencies f 0 in the dynamics and lowest wave numbers of spatial variations k 0 ∝ N −1 are obtained (Fig. 12), as is done in Sec. III for α = 0. Fig. 12 reveals how the spatio-temporal behaviour of evolving surfaces depend on parameter α.
The relations f 0 (k 0 ) shown in Fig. 12 indicate that for small k 0 , the power-law behaviour f 0 ∼ k ξ 0 observed in Fig. 10 for α = 0, also persists for α = 0 with exponent ξ decreasing with increasing α: from ξ ≈ 2.2 for α = −0.12 to ξ ≈ 0.9 for α = 2. However, for α = 0, this powerlaw behaviour flattens out at larger values of k 0 . For α = 5 (not shown in Fig. 12), the possible power-law behaviour is more difficult to determine, since the curve f 0 (k 0 ) appears flat almost through the whole range of k 0 , except for only two points with smallest k 0 -way less than enough to make conclusions.
One can interpret f 0 at some k 0 /2π = l −1 0 as the approximate rate of processes at the length scale l 0 , or n 0 ≡ l 0 /∆x in lattice units. Then the results displayed in Fig. 12 imply that at smaller scales -say, n 0 < 200 (k 0 ∆x/2π > 0.005 in Fig. 12) -the rate is monotonically decreasing with α. On the other hand, for larger scales, this does not hold any more. For example, for α = 0, 0.5, 1 the relation of between f 0 and α reverses (becomes monotonically increasing) already at n 0 > 300. For large enough scales, f 0 should become monotonically increasing with α for all values, at least in −0.12 ≤ α ≤ 2, if the power-law trends f 0 (k 0 ) ∝ k ξ 0 shown as straight lines in Fig. 12 continue for even larger systems, N > 1000.

B. Fits by two generalized Lorentzians
The fits of the PSDs W (f ) by a generalized Lorentzian plus a constant (24) seem to be suitable for most cases investigated for −0.12 ≤ α ≤ 5 with system sizes 100 ≤ N ≤ 1000. However, for α = 0.5 and α = 1, and system sizes N ≥ 700, some larger deviations from the fits can be observed. For example, Fig. 13 displays the apparent occurrence of a second hump in the PSD for α = 1 at N = 1000 which renders the fit (24) less suitable, although at smaller N it works very well (dotted lines in Fig. 13). In these cases, however, the sum of two generalized Lorentzians and a constant with the same exponent β, fits the PSD almost perfectly (orange long-dashed line in Fig. 13 and red dashed line in Fig. 14).
As can be seen in Fig. 14, the characteristic frequencies f 0 1 and f 0 2 of the two-generalized-Lorentzian fit (26) have the frequency f 0 of the original single-generalized-Lorentzian fit (24) between them, i.e., f 0 1 < f 0 < f 0 2 . Moreover, the frequency f 0 seems to follow the power-law trend (blue diamonds and dash-dotted line in Fig. 12), even if the fit is not that good as for smaller N values. Fig. 15 displays the PSDs with their fits (24) and (26) for the whole parameter α range investigated at system size N = 1000.

V. SUMMARY AND OUTLOOK
The results presented in this paper give some new insights into the complex spatio-temporal behaviour of surfaces produced by the two-dimensional generalized  Kuramoto-Sivashinsky equation (1) and might be interesting to a broader circle of researchers working in the field of continuum systems with complex nonlinear dynamics. The scaling properties (13) of the saturated surface roughness indicate that additional large scale height variations of scale free character appear when the system size N is increased. The dynamics of these slow height variations can be investigated by analysing the time series of the fluctuating surface roughness w(t) where the occurrence of slow modes with increasing system size can also be observed (see Fig. 6 and Fig. 7). This analysis shows that the resulting power-spectral densities (PSDs) can be expressed as the sum of a generalized Lorentzian and a constant, (24), or, in some cases, as two generalized Lorentzians (26), as shown in Figs. 9, 13 and 15.
The characteristic frequency f 0 obtained as a fit parameter corresponds to the smallest rate (largest time scale) that plays a role in the surface evolution. It can be attributed to the spatial mode of lowest wave number k 0 (which is inversely proportional to the system size) that can appear in the system. The dependence of this characteristic frequency on the system size gives the 'dispersion relation' f 0 (k 0 ) that connects spatial and temporal scales of surface dynamics. These relations have the power-law f 0 ∼ k ξ 0 character (see Fig. 12) for large systems (small k 0 ), thus, suggesting that the underlying temporal behaviour is scale free. Also, the exponent ξ is found to decrease with increasing value of parameter α. These results indicate, among other things, that although the characteristic time scale of dynamics on smaller scales decreases very strongly with increasing α, on large enough scales, this relation is reversed, i.e., the evolution on large scales is slower for smaller α.
The findings presented in this paper also raise some interesting questions for further research. For example, it is apparent from Fig. 11 and from the values of the PSD exponent β that the character of surface roughness dynamics depends quite strongly on parameter α. The question arises how temporal properties on various scales change with α and what are the statistical properties of the apparent bursts observed for larger values of α.
The Lorentzian shape, W (f ) ∼ (f 2 0 + f 2 ) −1 , of the PSD and relation f 0 ∼ k ξ 0 with ξ ≈ 2 for α ≈ 0 also suggests a possible analogy between the large-scale fluctuations of surface roughness and a diffusive process with the probability density Fourier transformed in space and time [50],P (k 0 , f ) ∝ k 2 0 (k 2 0 D) 2 + f 2 , where D is the diffusion constant independent of the system size. This correspondence becomes apparent when f 0 = Dk 2 0 is substituted in (21). Thus, perhaps the slow kinetics of the surface roughness might even be reproduced by a random walk of a particle in some external potential which is implied by the fact that the process w(t) is bounded and, consequently, k 0 does not go to zero for systems of finite size. For larger α values where the corresponding PSD exponent β ≈ 3 and ξ < 2 this process would then correspond to anomalous diffusion. Moreover, the fact that, for some parameter values, one more generalized Lorentzian has to be added to the in order to fit the calculated PSD for large systems (see Figs. 13 and 14) suggests the emergence of one more time scale, or perhaps, the whole interval of time scales. Any conclusive answers about both, the exact character and the occurrence mechanism, of this regime require more data obtained from simulations on even larger systems.