International Journal of Modern Physics C, c World Scientific Publishing Company

ANOMALIES IN SIMULATIONS OF NEAREST NEIGHBOR BALLISTIC DEPOSITION

Raissa M. D’Souza Department of Physics, and Lab for Computer Science, MIT Cambridge, MA, 02139, U.S.A. E-mail: raissa@mit.edu

Received (received date) Revised (revised date)

Ballistic Deposition (BD) is a prototypical model for interface growth and for dynamic scaling behavior in non-equilibrium systems. BD is typically investigated with computer simulations where randomness is replaced by use of deterministic Pseudo Random Number Generators (PRNGs). In this study of BD, several results discrepant with the prevailing paradigm, were observed. First, the value of the roughness exponent, χ, obtained is below the value for a random walk (i.e. χ < 1/2). The value χ = 1/2 is predicted by the KPZ equation, and many models of growth obtain this exponent. Second, height fluctuations of the growing interface appear not to satisfy simple scaling. Third, a decrease in the surface roughness is observed in a conjectured steady state regime. Computer implementations of BD may be responsible for the discrepancies. A coupling between the BD algorithm and a PRNG algorithm is identified, and statistically discrepant results are obtained for an implementation with a different PRNG.

1. Overview

Understanding growth patterns, both of clusters and solidification fronts, has become an increasingly interesting problem, relevant to non-equilibrium processes in general.^{1,2,3} Dynamic scaling characterizes many of these processes, and Ballistic Deposition^{4} (BD) is a prototypical model for this class of system, and for interface growth in general. In the BD model free particles, following ballistic trajectories, encounter the active growth interface of the substrate at which point they aggregate. The resulting growth patterns are compact clusters with a rough surface, which may be an accurate model of thin film growth. The active growth interface exhibits dynamic scaling behavior.

2 Anomalies in Simulations of Nearest Neighbor Ballistic Deposition

BD has been investigated for more than a decade and the findings have defined pathways for countless successive studies of growth processes.^{5,6,7} However it is disconcerting that a true consensus has never been reached on the value of the roughness exponent, χ (defined in detail below). Several anomalies were observed in the implementation of BD reported in this study, namely a growth surface not described by simple scaling, a roughness exponent below the theoretically predicted value, and a decrease in the surface roughness in a conjectured steady state regime. Coupling to PRNGs is observed, which may resolve the anomalies and discrepancies in past work, and challenges the implicit assumption that randomness can be replaced by deterministic Pseudo Random Number Generators (PRNGs).

The organization of the manuscript is as follows. It begins with a brief overview of the BD algorithm and of existing theories. The implementation details and results are then discussed. Pathological situations which may cause long crossover times or a change in exponents are investigated. For example, long range spatial and/or temporal correlations in the pseudorandom sequence of numbers employed are sought for, but no evidence of such two-point correlations is found. Yet a dynamics in the value of the surface roughness is observed in a conjectured steady state regime, and a different dynamics is observed for an implementation with a different PRNG. We thus conclude that the BD model couples sensitively to, as yet undetermined, non-randomness of pseudorandom sequences, even sequences which pass all standard statistical tests.

2. Review of past work

The surface configuration for BD is completely described by the height, h(x, t), along each position, x, of an underlying substrate, with t denoting the time duration of growth. Throughout this paper growth on a one-dimensional substrate (of length L) is considered, for which exact theoretical predictions exist. The surface evolves as follows. A column along the substrate is chosen at random and a particle is added to the surface of that column at the height:

h(x, t′ + 1) = max[h(x − 1, t′), h(x, t′) + 1, h(x + 1, t′)].  (1)

(The deposited particle occupies the highest empty site with one or more occupied nearest neighbor sites.) Here t′ is the number of individual deposition events and is proportional to t (t = t′/L). Beginning with an initially flat substrate, the width of the active growth interface, ξ(L, t), increases from zero to an asymptotic value which depends on the underlying, finite size, length scale, L. A measure of ξ(L, t) is the standard deviation of the surface heights, {h(x, t)},

ξ^{2}(L, t) = \frac{1}{L} ∑_{x=1}^{L}(h(x, t) − \overline{h(t)})^{2},  (2)

where \overline{h(t)} is the mean height of the surface at time t.

It was pointed out by Family and Vicsek^{8} that the scaling forms for the growth and saturation of the width of the interface can be described by a scaling ansatz, similar to that applicable to critical systems;

ξ(L, t) = L^{χ}f(t/L^{z}),  (3)

where f(x) ∼ x^{β} for x ≪ 1 and f(x) = const for x ≫ 1. For short times, the width of the interface should increase as ξ ∼ t^{β}. In the asymptotic regime the width of the interface should scale as ξ ∼ L^{χ}.

An analytical theory by Kardar, Parisi and Zhang (KPZ)^{9} describes the evolution of fluctuations on growing surfaces using a symmetry motivated differential equation:

\frac{∂h}{∂t} = ν∇^{2}h + \frac{λ}{2}(∇h)^{2} + η(x, t).  (4)

Here, ν is related to surface relaxations, (λ/2)(∇h)^{2} is introduced to account for lateral coarsening, and η(x, t) is white noise. In one dimension, the values of the scaling exponents can be obtained exactly from KPZ theory: χ = 1/2 and β = 1/3. Note that ξ(L, t) ∼ L^{1/2} is equivalent to a random walk. While these exponents agree very well with simulations of several growth models (e.g. Restricted Solid on Solid (RSOS)^{10}, discussed further in the conclusions), it is puzzling that the reported values of the roughness exponent, χ, obtained for the BD model have all been less than the theoretically predicted value. Moreover the reported values of χ have a substantial range (greater than 15% of the lowest value), (for a summary see Table 1).

Table 1. A table of scaling exponents for two growth models, BD and RSOS, as determined by selected numerical investigations. Error bars are included when available.

Scaling Exponents from Selected Numerical Simulations

Model  Reference  β  χ
BD  FV85^{8}  0.30 ± 0.02  0.42 ± 0.03
MRSB86^{11} 0.331 ± 0.006 0.47
M93^{5}  0.33  0.45
HHZ95^{7}  0.31  not reported
This study 0.31 ± 0.02  χ_{loc} 0.42 ± 0.02
χ_{glob} 0.455 ± 0.015
RSOS  KK89^{10}  0.332 ± 0.005 0.50
KPZ Theory KPZ86^{9}  1/3  1/2

3. Implementation and Results

In this study of BD, initially flat substrates of length varying from L = 127 to L = 2047 are considered. At each update a pseudo random number (PRN) is generated indicating which site along the substrate will have the next event. A particle is added to that column at a height described by Equation 1. The source for PRNs is a C library subroutine “random()”^{12}, which is a non-linear additive feedback PRNG, initialized with a 256 bit seed, giving a repeat period of 2^{64} numbers. The lowest 20 to 24 bits of each returned number were shifted off, leaving only the highest (and least correlated) bits. Many PRNGs were investigated before this variant of random() was chosen, as it performed best in initial tests. All of the simulations were carried out on a desktop workstation, with the runs on the shortest length substrates requiring a few hours, and the runs on the longest, a few days.

The dynamic scaling exponent, β, is determined using a plot of the width of the interface, ξ(L, t), versus the time into the simulation, for all of the substrate lengths L. Note one unit of time is L particle depositions. Consistent with previous studies (see Table 1) the results are β = 0.29 ± 0.01 for times 3 ≤ t < 100, and β = 0.31 ± 0.02 for times 100 ≤ t < 2000. We chose to report the value for the longer times, as the scaling interval is of greater absolute size.

Determination of the asymptotic roughness exponent χ is not as straightforward. The primary complication is that ξ ∼ L^{χ} only in the regime where the correlation length has reached the full system size (t ≫ L^{z}), which can be computationally prohibitive for large substrate lengths. The conservative estimate of time for “relaxation”, τ = 10L^{z= \frac{5}{3}} is employed (which exceeds the expected scaling of z = 3/2). This time corresponds to an average surface height h̄ = 20L^{\frac{5}{3}} .

Once in this regime, the local roughness of the surface is investigated by studying “windows” of length l < L (i.e ξ(l, t) for l < L). Finite size effects (imposed by periodic boundary conditions) cause a rounding of these curves, restricting the scaling regime to lengths l < L/2, where self-affine scaling suggests ξ(l, t) ∼ l^{χ}. The BD local surface roughness is compared to that of a random walk. We calculate the expectation value of the local surface roughness of a random walk over the ensemble of all possible walks of total length L:

< ξ(l, t)_{rwalk}^{2} >=< \frac{1}{l} ∑_{x=1}^{l} (h(x) − \overline{h_{l}})^{2} > .  (5)

Here h(x + 1) = h(x) + δ_{x}, δ_{x} = ±1, ∑_{x=1}^{(L−1)} δ_{x} = 0 (which implements periodic boundary conditions), and the correlations between δ′_{x}s are assumed uniform. This yields the scaling relation

< ξ(l, t)_{rwalk}^{2} >∼ l (1 − \frac{l}{2L}) ,  (6)

for l, L ≫ 1. The second term in Equation 6 represents a finite size correction to scaling. Comparison with our results for BD is shown in Figure 1, where the prefactor for the random walk in Equation 6 was selected so as to agree with the BD simulations at the shortest lengths. The two curves are very distinct; Both appear to exhibit scaling behavior over more than one decade, however the scaling exponents differ. If the δ′_{x}s are chosen in agreement with the empirically obtained step height distribution of the BD simulations, no ad hoc prefactor is needed and similar results are obtained.

From curves of ξ(l, t) in the asymptotic regime, three distinct scaling behaviors can be identified. For 3 ≤ l < 20 the relation ξ(l, t) ∼ l^{χ_{loc_{0}}} is obtained, with χ_{loc_{0}} = 0.35 ± 0.01. For 30 < l < 400 the relation ξ(l, t) ∼ l^{χ_{loc}} is obtained, with χ_{loc} = 0.42 ± 0.02. Looking only at the longest length from each substrate the relation ξ(L, t) ∼ L^{χ_{glob}} is obtained, with χ_{glob} = 0.455 ± 0.015. The data for all substrate lengths is shown in Figure 2, along with a comparison to the equivalent data for a random walk with periodic boundary conditions. For the random walk, only one scaling relation can be extracted ξ(l, t) = l^{1/2}, for both l ≤ L/2 and l = L.

In an attempt to formulate a consistent scaling picture, we collapse the data from the five different length substrates, using Equation 3 (see Figure 3). For

Fig. 1. The width of the growth interface plotted as a function of increasing window size, for the BD simulations and a theoretical calculation of a random walk with periodic boundary conditions.

Fig. 2. The width of the growth interface plotted as a function of increasing window size for the longest times simulated. The data from all five substrate lengths are included. Dependent upon the range of lengths we examine, three apparent scaling exponents may be obtained, as indicated by the solid lines. The insert is the corresponding data for a random walk with periodic boundary conditions.

ξ(L, t) collapse is achieved by plotting ξ(L, t)/L^{χ_{glob}} vs t/L^{z}, with the exponents χ_{glob} = 0.45 ± 0.02 and z = 1.45 ± 0.03. For ξ(l, t), with l = 0.1L, collapse is achieved by plotting ξ(l, t)/l^{χ_{loc}} vs t/l^{z}, with the exponents χ_{loc} = 0.40 ± 0.02 and z = 1.41 ± 0.03.

Fig. 3. The scaling function, f( \frac{t}{l^{z}}), determined for data sampled at l = 0.1L and at l = L. For the data sampled at l = 0.1L collapse is achieved with the exponents χ_{loc} = 0.40 ± 0.02 and z = 1.41 ± 0.03. For data sampled at l = L collapse is achieved with χ_{glob} = 0.45 ± 0.02 and z = 1.45 ± 0.03

Exponents obtained by the collapse are consistent with those obtained by a linear fit and with those reported in previous studies (see Table 1). The measurement of χ_{loc} agrees with the FV85^{8} measurement of χ. The measurement of χ_{glob} is consistent with the M93^{5} and the MRSB86^{11} measurements of χ. However, these results are puzzling on three counts. First χ_{loc} ≠ χ_{glob} – the surface is not self affine. Second, all of the exponents obtained are below the value predicted by the KPZ theory. Third, the characteristic identity χ + z = 2 is not obeyed. These are not only violations of the KPZ scaling, but also most of its extensions (e.g. Medina, et al.^{13}). Thus we are lead to ask, what kind of continuum model could explain a decrease in the values of the scaling exponents (i.e. a hypo-rough surface)?

4. Possible resolutions

One possibility is that the model suffers from a long crossover regime. For example the Wolf-Villain model of growth has a length dependence in the adjacent step height distribution, which persists for long times (t ∼ 10^{4}), and long lengths (L ∼ 256)^{14}. The step height distribution for the BD model is a well behaved quantity, which reaches a steady state, length independent, value within the completion of a few monolayers. With reference to the work on crossover, it has been suggested that the structure factor is the only way of accurately determining exponents.^{15} The structure factor is related to the Fourier transform of the height-height correlation function

S(k, t) =< ĥ(k, t)ĥ(−k, t)) >,  (7)

with ĥ(k, t) = L^{−d/2} ∑_{x}[h(x, t) − \overline{h(t)}]e^{ik·x}. Crossover behavior is manifested by a change in the scaling exponent of S(k, t) vs k. There is no evidence of crossover in S(k, t) for our data, including the longest lengths and longest times simulated. The exponent obtained is consistent with χ_{glob} = 0.45.

A second possibility is an intrinsic width correction,^{16} as introduced by Kertesz and Wolf to the Eden Model. This correction accounts for voids in the bulk, and allows for a clear scaling regime to be obtained for the Eden model in 2+1 and 3+1 dimensions.^{17} Introducing such a correction to our data destroys the scaling regime.

A third possibility is correlated noise. There has been extensive past theoretical and numerical work into growth models with positively correlated noise (Medina, et al., is the first study).^{13,5,7} For noise with long-range correlations in space and/or time, the noise term in Equation 3 is of the form < η(x, t) >= 0, with

< η(k, ω)η(k′, ω) >= 2D(k, ω)δ^{d−1}(k + k′)δ(ω + ω′),  (8)

where the noise spectrum D(k, ω) has power-law singularities of the form D(k, ω) ∼ |k|^{−2ρ}ω^{−2θ}. For uncorrelated noise, the noise spectrum is a constant.

As shown in the original study,^{13} for spatially correlated noise (i.e ρ > 0, θ = 0), χ increases and z decreases, preserving the identity of Galilean invariance, χ+z = 2. An extension to anti-correlated noise (i.e. (ρ, θ) < 0) may show a decrease in the values of the scaling exponents and a breaking of that identity. It is clear that extreme forms of anti-correlated noise will result in a flat interface; Consider sequential updating of a BD growth algorithm, this generates a perfectly flat interface at each timestep. Hence it is reasonable that correlations in the noise may change the roughness of the surface.

We have conducted extensive tests to identify potential long-range spatial and/or temporal correlations in the PRNs generated by the function call random(). The PRNs produced pass every mathematical and physically motivated test, including those directly relevant to the growth algorithm. The Fourier transform of the noise in each column, and along the columns of the substrate, produce flat, white-noise spectra. The number of calls to each column are Poisson distributed (as discussed below). Waiting times between successive calls to each column are exponentially distributed. The simulated space was divided into checkerboard sublattices, but no discrepancy between events on even or odd numbered sites was found. The autocorrelation function for every property of the generated surface tested is a simple decaying exponential. There is no bias for relevant successive events to be in the same column or left or right neighboring columns. Other researchers have recently pointed out some physical models which manifest the pathologies of certain PRNGs.^{18} The PRNs produced by random() pass even those sensitive tests.

Nonetheless, a systematic decrease in the width of the BD growth interface is observed in the “asymptotic” regime (Figure 4). For the longer lengths simulated, the decrease occurs even before the strictly defined asymptotic regime. In order to bound the decrease outside of statistical error, 200 independent runs on substrates of lengths L = 127 and L = 511 are studied. Note that all the data reported so far were obtained with 20 independent runs on each substrate length, and all error bars reported are statistical. A plot of ξ(l, t) vs l, during the four greatest times, on the L = 127 substrate, is shown in Figure 4. As the time increases, ξ(l, t) decreases systematically. The error bars should be noted; The value of ξ(l, t) for the longest time is over four standard errors away from the value for the shortest time. This effect may be due to anti-correlation in the system (where system refers to the BD algorithm coupled to the PRNG algorithm). The effect of anti-correlation in the PRNG on a growth algorithm is illustrated by a simple Poisson process, the Random Deposition (RD) model. Implementing RD also serves as a check of the BD computer code.

Fig. 4. The width of the growth interface versus increasing window size, plotted for subsequent time steps into the BD simulation. It can be seen that as time increases, the width of the interface decreases systematically. The inset plot is the width of the interface for the largest window size, l = L, plotted with increasing time. Note the error bars included on all points.

In the RD model the site filled by deposition is the next available height in the active column (not even nearest neighbor interactions exist). No boundary conditions are needed and hence there is no dependence on L. The surface height values should theoretically be Poisson distributed: (ξ ∼ t^{1/2}, for all t). Figure 5 is of ξ(t) vs t for an RD algorithm implemented by altering only one line of our BD code (that line describing Equation 1). Two different PRNGs were used: random() (described earlier), and rand(), a 16 bit version of the standard C-library subroutine.^{19} At the close of one repeat cycle, rand() has sampled all numbers evenly, yielding a flat surface for RD. The results using random() agree with theory out to the longest times simulated, yet a systematic decrease is observed using rand(), as it is for the BD model using random(). As mentioned, the repeat period of random() is on the order of 10^{19}, almost ten orders of magnitude greater than the total number of calls made to it.

Fig. 5. The width of the growth interface versus the number of function calls for simulations of RD using two different PRNGs. The solid line corresponds to the result for ideal random numbers, ξ ∝ t^{1/2}. The cycle length of rand() is apparent since the interface width decreases to zero at the end of a cycle.

To quantify that the decrease of the surface roughness is a result of evolving the system with random() the full BD simulations were run, using the identical code, but a different PRNG, ran2()^{20} (which combines two distinct types of PRNGs in an attempt to eliminate correlations inherent to each one separately). The systematic decrease was not observed, but several unsystematic fluctuations in the value of ξ(l, t), outside of statistical error bounds, were observed. Selecting individual time samples for each length substrate, it is possible to construct a scaling curve, with the scaling exponents χ_{loc} = 0.45 ± 0.01, and χ_{glob} = 0.51 ± 0.02. Note that these are statistical error bounds, which do not address the issue of fluctuations. Constructing the curve with statistical outliers, changes the values of the exponents. We do not assert that these exponents are the “true exponents”, but instead wish to focus on the fact that different PRNGs yield different results for an identical system, implemented with the identical code.

As such, the more revealing comparison is between average values obtained by the two different PRNGs. For example, the average asymptotic values of ξ(L, t) differ by more than three standard errors. A “t-test” comparing these average values fails at the 99% confidence level^{21}. Statistically discrepant results for average values and a statistically distinct dynamics for the interface width fluctuations show a breakdown of basic sampling assumptions, and moreover, that the observed dynamics is not inherent to the BD model. Assuming the dynamics does not reside in the PRNGs, it must reside in a coupling between the BD and PRNG algorithms. This is in line with the observation that the short time scaling exponent (β) is not effected, but the asymptotic scaling exponent (χ) is; The effects of coupling may take some time accumulate. A detailed statistical analysis quantifying the breakdown of sampling assumptions and ergodic exploration of phase space is reported elsewhere.^{22}

5. Discussion and Conclusions

It is reasonable that the BD algorithm is more sensitive to correlations in PRNs than standard Monte Carlo (MC) algorithms. In standard MC, comparison to the Boltzmann probability causes rejection of some PRNs produced. In BD, all PRNs are used (in the sequence produced). It should be noted that in restricted models of growth (where physical constraints cause rejection of PRNs), such as RSOS, the theoretically predicted scaling exponents are recovered with great precision in numerical simulations. It may be interesting to implement the BD algorithm with random rejection of PRNs. In addition, we have been considering the use of massive physical simulations (run on a special purpose cellular automata machine^{23}) which have a vast amount of initial state, as a source of randomness.

BD is a sensitive physical test for correlations present in pseudorandom sequences, and it would be desirable to identify the exact nature of the correlations detected. Discrepancies in the reported values for the roughness exponent χ, and anomalies found in this study, may be attributed to distinct couplings between the BD and the PRNG algorithms. Disagreement in results generated by two different PRNGs is strong evidence for this and, more importantly, shows the dynamics observed in the asymptotic regime are not inherent to the BD model itself. Results from previous studies of BD have not indicated which PRNG was used, in addition many of these past simulations utilized power of two substrate lengths (L = 2^{n}), a system size for which PRNGs manifest their greatest pathologies.^{24} As of yet we are not able to identify the universality class of the BD model. A non-self affine growth surface may be due to unidentified crossover behavior, but it appears that coupling to the PRNG algorithms becomes dominant before the steady state regime is achieved.

Acknowledgements

The author would like to thank Y. Bar-Yam, M. Kardar, and M. Smith for countless, helpful discussions. This work was supported by NSF grants, numbers DMR-93-03667 and DMS-95-96217, and by DARPA contract number DABT63-95-C-0130.

References

1. F. Family and T. Vicsek, Dynamics of Fractal Surfaces (World Scientific, Singapore, 1991).

2. A. Barabasi and U. E. Stanley, Fractal Concepts in Surface Growth (Press Syndicate of the Univ of Cambridge, New York, 1995).

3. C. Godrèche, (ed.) Solids Far From Equilibrium (Cambridge University Press, Cambridge, 1992).

4. M. J. Vold, Jour. Colloid Sci. 14, 168 (1959).

5. P. Meakin, Physics Reports 235, 189 (1993).

6. J. Villain, Journal de Physique I 1, 19 (1991).

7. T. Halpin-Healy and Y. Zhang, Physics Reports 254(4-6), 215 (1995).

8. F. Family and T. Vicsek, Jour. of Phys. A 18, L75 (1985).

9. M. Kardar, G. Parisi, and Y. Zhang, Phys. Rev. Let. 56(9), 889 (1986).

10. J. Kim and J. Kosterlitz, Phys. Rev. Let. 62(19), 2289 (1989).

11. P. Meakin, P. Ramanlal, L. Sander, and R. Ball, Phys. Rev. A 34(4), 5103 (1986).

12. GNU C Library, Sun Release 4.1 (Original code is copyright Regents of the University of California, 1983).

13. E. Medina, T. Hwa, M. Kardar, and Y. Zhang, Phys. Rev. A 39(6), 3053 (1989).

14. M. Schroeder, J. Shore, M. Siegert, D. Wolf, and M. Plischke, Europhys. Lett. 24(7), 563 (1993).

15. M. Siegert, Phys. Rev. E 53(4), 3209 (1996).

16. J.G. Zabolitzky and D. Stauffer, Phys. Rev. A 34(2), 1523 (1986).

17. J. Kertesz and D. Wolf, J. Phys. A 21, 747 (1988).

18. I. Vattulainen, T. Ala-Nissila, and K. Kankaala, Phys. Rev. Let. 73, 2513 (1994).

19. GNU C Library, Sun Release 4.1

20. W.H. Press etal., Numerical Recipes in C: the art of scientific computation, 2nd ed. (Cambridge University Press, 1992).

21. W. Feller, Introduction to Probability and Statistics (John Wiley & Sons, Inc., New York, 1979).

22. R.M. D’Souza, Y. Bar-Yam, and M. Kardar, Manuscript in progress.

23. N. Margolus, in A. Lawniczak, and R. Kapral, (eds.) Pattern Formation and Lattice-Gas Automata (American Mathematical Society, 1994).

24. D. E. Knuth, The Art of Computer Programming (Addison-Wesley, Reading, 1973).