# Singularity spectra of rough growing surfaces from wavelet analysis

## 1 Introduction

Inspired by the great technological importance of epitaxial crystal growth, the past decade has raised much theoretical research in the subject of kinetic roughening of surfaces during growth. The investigation of this effect, which is undesirable in practical applications, promises deep insight into statistical physics far from thermal equilibrium, see e.g.  for an overview. We focus on a full-diffusion Monte-Carlo model of homoepitaxial growth of a hypothetical material with simple cubic lattice structure under solid on solid conditions, i.e. the effects of overhangs and displacements are being neglected. Then, the crystal can be described by a two-dimensional array of integers which denote the height $`f(\stackrel{}{x})`$ of the surface. On each site, new particles are deposited with a rate $`r_a`$. Particles on the surface are hopping to nearest neighbour sites with Arrhenius rates $`\nu _0\mathrm{exp}((E_b+nE_n)/(k_bT))`$, where $`E_b`$ and $`E_n`$ are the binding energies of a particle to the substrate and to its $`n`$ nearest neighbours. $`\nu _0`$ is the attempt frequency, and $`k_bT`$ has its usual meaning. In contrast to earlier investigations of similar models , we permit the desorption of particles from the surface with rates $`\nu _0\mathrm{exp}((E_d+nE_n)/(k_bT))`$, where $`E_d>E_b`$.

The aim of this publication is twofold: We will first discuss the advantages of the wavelet analysis compared to the structure function (SF) approach, which has to date solely been used in the investigation of multiaffine surfaces. Then, we will apply this formalism to investigate the influence of desorption on kinetic roughening. We conclude with some remarks on the relevance of universality classes for our results.

## 2 Scaling concepts

The standard approach of dynamic scaling  assumes that the statistical properties of a growing surface before saturation remain invariant under a simultanous transformation of spatial extension $`\stackrel{}{x}`$, heigth $`f(\stackrel{}{x})`$ and time $`t`$,

$$\stackrel{}{x}b\stackrel{}{x}=\stackrel{}{x}^{};fb^\alpha f=f^{};tb^zt=t^{};$$
(1)   
where $`b`$ is an arbitrary positive constant. This implies, that a part of the surface smaller than the correlation length $`\xi (t)t^{1/z}`$ can be regarded as self-affine with Hurst exponent $`\alpha `$. A popular method of measuring $`\alpha `$ uses heigth-heigth correlation functions of (theoretically) arbitrary order $`q`$:

$$G(q,\stackrel{}{l},t):=\left|f(\stackrel{}{x},t)f(\stackrel{}{x}+\stackrel{}{l},t)\right|^q_\stackrel{}{x}l^{q\alpha }g(l/\xi (t)),$$
(2)   
where $`g(x)\text{const.}`$ for $`x0`$ and $`g(x)\text{const.}x^{q\alpha }`$ for $`x\mathrm{}`$. In practice, $`q=2`$ is the most common choice.

In principle, there are two different ways to measure $`\alpha `$: The local approach determines $`\alpha `$ from the initial slope of $`\mathrm{ln}(G(q,\stackrel{}{l},t))`$ versus $`\mathrm{ln}(l)`$ for small $`l`$. The global approach analyzes the dependence of the surface width $`w=\sqrt{(f(\stackrel{}{x},t)f(\stackrel{}{x},t))^2_\stackrel{}{x}}`$ in the saturation regime on the system size $`N`$: $`w_{sat}(N)N^{\alpha _g}`$. Before saturation, the surface width increases like $`wt^\beta `$, where $`\beta =\alpha /z`$. An alternative which avoids the simulation of different system sizes uses the complete functional dependence of equation 2: $`\alpha _g`$ and $`z`$ are chosen such that the curves of $`G(2,\stackrel{}{l},t)/l^{2\alpha _g}`$ versus $`l/t^{1/z}`$ collapse on a unique function $`g`$ within a large range of $`t`$ and $`l`$.

However, a careful analysis of simulation data  has shown, that several models of epitaxial growth show significant deviations from this simple picture. First, one obtains different values of $`\alpha `$ from the local than from the global approach, a phenomenon which is called anomalous scaling. Second, one often finds multiscaling: height-height correlation functions of different order yield a hirarchy of $`q`$-dependent exponents $`\alpha (q)`$, when determined from the initial power-law behaviour of $`G(q,\stackrel{}{l},t)`$.

These observations can be interpreted within the mathematical framework of multifractality: The Hölder exponent  $`h(\stackrel{}{x_0})`$ of a function $`f`$ at $`\stackrel{}{x}_0`$ is defined as the largest exponent such that there exists a polynomial of order $`n<h(\stackrel{}{x}_0)`$ and a constant $`C`$ which yield $`|f(\stackrel{}{x})P_n(\stackrel{}{x}\stackrel{}{x}_0)|C|\stackrel{}{x}\stackrel{}{x}_0|^{h(\stackrel{}{x}_0)}`$ in the neighbourhood of $`\stackrel{}{x}_0`$. The Hölder exponent is a local counterpart of the Hurst exponent: a self-affine function with Hurst exponent $`\alpha `$ has $`h(\stackrel{}{x})=\alpha `$ everywhere. However, in the case of a multiaffine function different points $`\stackrel{}{x}`$ might be characterized by different Hölder exponents. This general case is characterized by the singularity spectrum $`D(h)`$, which denotes the Hausdorff dimension of the set of points, where $`h`$ is the Hölder exponent of $`f`$.

## 3 The wavelet approach to multifractality

There is a deep analogy between multifractality and thermodynamics , where the scaling exponents play the role of energy, the singularity spectrum corresponds to entropy, and $`q`$ plays the role of inverse temperature. So, theoretically $`D(h)`$ might be calculated via a Legendre transform of $`\alpha (q)`$: $`D(h)=\text{min}_q(qhq\alpha (q)+2)`$ , a method which has been called structure function (SF) approach. However, its practical application raises fundamental difficulties: First, to obtain the complete singularity spectrum, one needs $`\alpha (q)`$ for positive and negative $`q`$. But as $`|f(\stackrel{}{x},t)f(\stackrel{}{x}+\stackrel{}{l},t)|`$ might become zero, $`G(q,\stackrel{}{x},t)`$ is in principle undefined for $`q<0`$. Therefore, only the left, ascending part of $`D(h)`$ is accessible to this method. Additionally, the results of the SF method can easily be corrupted by polynomial trends in $`f(\stackrel{}{x})`$ . It might be due to these difficulties, that - to our knowledge - no attempt to determine the singularity spectrum of growing surfaces from $`\alpha (q)`$ has ever been made. Although it has been argued that the $`\alpha (q)`$ collapse onto a single $`\alpha `$ in the limit $`t\mathrm{}`$, which characterizes the asymptotic universality class of the model , we are convinced that deeper insight into fractal growth on experimentally relevant finite  timescales can be gained from a detailed knowledge of the $`D(h)`$ spectrum.

To this end, we follow the strategy suggested by Arnéodo et. al , which circumvents the problems of the SF approach and permits a reliable measurement of the complete $`D(h)`$. Mathematically, the wavelet transform of a function $`f(\stackrel{}{x})`$ of two variables is defined as its convolution with the complex conjugate of the wavelet $`\psi `$, which is dilated with the scale $`a`$ and rotated by an angle $`\theta `$ :

$$T_\psi [f](\stackrel{}{b},\theta ,a)=C_\psi ^{1/2}a^2d^2x\psi ^{}(a^1𝐑_\theta (\stackrel{}{x}\stackrel{}{b}))f(\stackrel{}{x}).$$
(3)   
Here $`𝐑_\theta `$ is the usual 2-dimensional rotation matrix, and $`C_\psi =(2\pi )^2d^2k|\stackrel{}{k}|^2|\widehat{\psi }(\stackrel{}{k})|^2`$ is a normalization constant, whose existence requires square integrability of the wavelet $`\psi (\stackrel{}{x})`$ in fourier space. Apart from this constraint, the wavelet can (in principle) be an arbitrary complex-valued function. Introducing the wavelet $`\psi _\delta (\stackrel{}{x})=\delta (\stackrel{}{x})\delta (\stackrel{}{x}+\stackrel{}{n})`$, where $`\stackrel{}{n}`$ is an arbitrary unit vector, one obtains easily

$$T_{\psi _\delta }[f](\stackrel{}{b},\theta ,a)=C_{\psi _\delta }^{1/2}\left[f(\stackrel{}{b})f(\stackrel{}{b}+a𝐑_\theta \stackrel{}{n})\right]d^2b\left|T_{\psi _\delta }[f](\stackrel{}{b},\theta ,a)\right|^qG(q,a𝐑_\theta \stackrel{}{n}).$$
(4)   
Consequently, a calculation of the moments of the wavelet transform of the surface yields the SF approach as a special case. To avoid its weaknesses, two major improvements are necessary:

First, we use a class of wavelets with a greater number of vanishing moments $`n_\stackrel{}{\mathrm{\Psi }}`$ than $`\psi _\delta (\stackrel{}{x})`$. This increases the range of accessible Hölder exponents and improves the insensitivity to polynomial trends in $`f(\stackrel{}{x})`$. We introduce a two-component version of the wavelet transform

$$\stackrel{}{T}_\stackrel{}{\mathrm{\Psi }}[f](\stackrel{}{b},a)=\frac{1}{a^2}d^2x\left(\begin{array}{c}\mathrm{\Psi }_1(a^1(\stackrel{}{x}\stackrel{}{b}))\\ \mathrm{\Psi }_2(a^1(\stackrel{}{x}\stackrel{}{b}))\end{array}\right)f(\stackrel{}{x}),$$
(5)   
where the analyzing wavelets $`\mathrm{\Psi }_1`$, $`\mathrm{\Psi }_2`$ are defined as partial derivatives of a radially symmetrical convolution function $`\mathrm{\Phi }(\stackrel{}{x})`$: $`\mathrm{\Psi }_1(\stackrel{}{x})=\mathrm{\Phi }/x`$, $`\mathrm{\Psi }_2(\stackrel{}{x})=\mathrm{\Phi }/y`$. Then $`\stackrel{}{T}_\stackrel{}{\mathrm{\Psi }}[f](\stackrel{}{b},a)`$ can be written as the gradient of $`f(\stackrel{}{x})`$, smoothed with a filter $`\mathrm{\Phi }`$ w.r.t. $`\stackrel{}{b}`$. This definition becomes a special case of equation 3, when multiplied with $`\stackrel{}{n}_\theta =(\mathrm{cos}(\theta ),\mathrm{sin}(\theta ))^{}`$, yet allows for an easier numerical computation<sup>1</sup><sup>1</sup>1For simplicity, the irrelevant constant $`C_\mathrm{\Psi }`$ has been omitted. For example, $`\mathrm{\Phi }`$ can be a gaussian, where $`n_\stackrel{}{\mathrm{\Psi }}=1`$, or $`\mathrm{\Phi }_1(\stackrel{}{x})=(2\stackrel{}{x}^2)\mathrm{exp}(\stackrel{}{x}^2/2)`$, which has two vanishing moments.

Second, the integration over $`\stackrel{}{b}`$ in equation 4 is undefined for $`q<0`$, since the wavelet coefficients might become zero. The basic idea is to replace it with a discrete summation over an appropriate partition of the wavelet transform which obtains nonzero values only, but preserves the relevant information on the Hölder regularity of $`f(\stackrel{}{x})`$. In the following, we will give a brief outline of the rather involved algorithm and refer the reader to  for more details and a mathematical proof. The wavelet transform modulus maxima (WTMM) are defined as local maxima of the modulus $`M_\stackrel{}{\mathrm{\Psi }}[f](\stackrel{}{b},a):=|\stackrel{}{T}_\stackrel{}{\mathrm{\Psi }}[f](\stackrel{}{b},a)|`$ in the direction of $`\stackrel{}{T}_\stackrel{}{\mathrm{\Psi }}[f](\stackrel{}{b},a)`$ for fixed $`a`$. These WTMM lie on connected curves, which trace structures of size $`a`$ on the surface. The strength of each is characterized by the maximal value of $`M_\stackrel{}{\mathrm{\Psi }}[f](\stackrel{}{b},a)`$ along the line, the so-called wavelet transform modulus maxima maximum (WTMMM) . While proceeding from large to small $`a`$, successively smaller structures are resolved. Connecting the WTMMM at different scales yields the set $``$ of maxima lines $`l`$, which lead to the locations of the singularities of $`f(\stackrel{}{x})`$ in the limit $`a0`$. The partition functions

$$Z(q,a)=\underset{l(a)}{}\left(\underset{(\stackrel{}{b},a^{})l,a^{}a}{sup}M_\stackrel{}{\mathrm{\Psi }}[f](\stackrel{}{b},a^{})\right)^qa^{\tau (q)}\text{for}a0$$
(6)   
are defined on the subset $`(a)`$ of lines which cross the scale $`a`$. From the analogy between the multifractal formalism and thermodynamics, $`D(h)`$ is calculated via a legendre transform of the exponents $`\tau (q)`$, which characterize the scaling behaviour of $`Z(q,a)`$ on small scales $`a`$: $`D(h)=\text{min}_q(qh\tau (q))`$. Additionally, $`\tau (q)`$ itself has a physical meaning for some $`q`$: $`\tau (0)`$ is the fractal dimension of the set of points where $`h(\stackrel{}{x})<\mathrm{}`$, while the fractal dimension of the surface $`f(\stackrel{}{x})`$ itself equals $`\text{max}(2,1\tau (1))`$.

## 4 Results

In our simulations, we choose the parameters $`\nu _0=10^{12}/s`$, $`E_b=0.9\text{eV}`$ and $`E_n=0.25\text{eV}`$, and a temperature $`T=450\text{K}`$. To study the influence of desorption, we consider three models with different activation energies $`E_d`$: in model A desorption is forbidden, i.e. $`E_d=\mathrm{}`$. Models B and C have $`E_d=1.1\text{eV}`$ and $`E_d=1.0\text{eV}`$. We simulate the deposition of $`210^4`$ monolayers at a growth rate of one monolayer per second on a lattice of $`N\times N`$ unit cells using periodic boundary conditions, our standard value being $`N=512`$. To check for finite size effects, we have also simulated $`N=256`$. In all presented results averages over 6 independent simulation runs have been performed. Although we have used an optimized algorithm, these simulations consumed several weeks of CPU time on our workstation cluster.

First, we have checked our results for artifacts resulting from properties of the analyzing wavelet rather than from the analyzed surface by using different convolution functions $`\mathrm{\Phi }_n`$: $`\mathrm{\Phi }_0`$ is the gaussian function, $`\mathrm{\Phi }_n,n1`$ are products of gaussians and polynomials, which have been chosen in a way that the first $`n`$ moments vanish. Then, the analyzing wavelets have $`n_{\stackrel{}{\mathrm{\Psi }}_n}=n+1`$ vanishing moments. We find (figure 1 a), that the $`\tau (q)`$-curve obtained with $`\mathrm{\Phi }_0`$ deviates significantly from those obtained with $`\mathrm{\Phi }_1`$, $`\mathrm{\Phi }_2`$ and $`\mathrm{\Phi }_3`$. The latter agree apart from small differences which are mainly due to the discrete sampling of the wavelet in the numerical implementation of the algorithm. This is explained by the theoretical result  that $`d\tau (q)/dq=n_\stackrel{}{\mathrm{\Psi }}`$ for $`q<q_{crit.}<0`$ if the number of vanishing moments of the analyzing wavelet is too small. Consequently, the agreement of the other curves proves their physical relevance.

Figure 1 b shows averages of $`\tau (q)`$ curves obtained with the convolution functions $`\mathrm{\Phi }_1`$, $`\mathrm{\Phi }_2`$, $`\mathrm{\Phi }_3`$ from surfaces after $`210^4s`$ of growth on an initially flat substrate. For all our models, their nonlinear behaviour reflects the multiaffine surface morphology. From the fact that these curves are reproduced within statistical errors in simulations with $`N=256`$, we conclude that finite size effects can be neglected. Clearly, desorption reduces the slope of $`\tau (q)`$, although only a small fraction of the incoming particles is desorbed: $`0.18\%`$ in model B and $`2.57\%`$ in model C with slightly higher values at earlier times.

The corresponding singularity spectra are shown in figure 2 a. They have a typical shape whose descending part seems to be symmetrical to the ascending part and which changes at most slightly, while the whole spectra are shifted towards smaller Hölder exponents as desorption becomes more important. We emphasize that we find no evidence for a time dependence of the singularity spectra within the range $`9700st210^4s`$, so that our results do not support the idea of an asymptotic regime characterized by a single exponent $`\alpha `$. However, the accessible time range of computer simulations is limited, so that we cannot finally disprove the existence of such a regime.

The multifractal formalism has replaced the unique scaling exponent $`\alpha `$ of spatial extension in the simple picture of dynamic scaling (equation 1) with a wide spectrum of Hölder exponents. By analogy, one might find it necessary to replace the scaling exponent $`\beta `$ with a distribution of temporal counterparts of $`h`$. To answer this question, we investigate the propability distribution function (PDF) $`P(ff,t)`$ of surface heights. Dynamical scale invariance with a single $`\beta `$ demands that

$$P(ff,t)=\stackrel{~}{P}\left(\frac{ff}{t^\beta }\right)\frac{1}{t^\beta },$$
(7)   
i.e. the rescaled PDFs $`Pt^\beta `$ should collapse onto a single function $`\stackrel{~}{P}`$ when plotted as a function of $`(ff)/t^\beta `$ within a large time range.

We measure $`\beta `$ from the increase of the surface width with time, which follows a power law for $`t150s`$ in models A and B ($`\beta _A=0.19\pm 0.01`$, $`\beta _B=0.17\pm 0.01`$) respectively $`150s<t<7500s`$ in model C ($`\beta _C=0.11\pm 0.01`$), which then starts to approach the final saturation regime. The high quality of the data collapse of the PDFs shown in figure 2 b proves that the scaling form 7 holds, showing that a single exponent describes the scaling behaviour of $`P(ff,t)`$. This parallels the finding of Krug in  for the one-dimensional Das Sarma-Tamborenea model.

Finally, the WTMM method, which is a precise tool to investigate local scaling properties of surfaces might help to get some insight into the phenomenon of anomalous scaling. The conventional picture  notes the difference between the global $`\alpha _g`$ and a “local $`\alpha `$” which is determined from the power-law behaviour of $`G(2,\stackrel{}{l},t)`$ for small $`l`$, and, within the multifractal formalism, simply corresponds to a Hölder exponent on the ascending part of the singularity spectrum. We have determined the global scaling exponents $`\alpha _g`$ and $`z`$ from the data collapse of the scaled height-height correlation function $`G(2,\stackrel{}{l},t)`$ and find agreement within statistical errors between $`\alpha _g`$ and that value of the Hölder exponent $`h_m`$ which maxmizes $`D(h)`$ (table 1). This empirical result can be explained with a saddle-point argument: We calculate the surface width

$$w^2=\frac{1}{N^2}d^2x(f(\stackrel{}{x})f)^2=\frac{1}{N^2}𝑑\stackrel{~}{h}\underset{I(\stackrel{~}{h})}{\underset{}{d^2x\delta (\stackrel{~}{h}h(\stackrel{}{x}))(f(\stackrel{}{x})f)^2}}.$$
(8)   
Since $`I(\stackrel{~}{h})`$ grows like $`N^{D(\stackrel{~}{h})}`$ with the system size, in large systems the integral over $`\stackrel{~}{h}`$ will be dominated by $`I(h_m)`$. That means, that $`w`$ and therefore the global scaling properties of the surface are governed by the subset of points, which has the greatest fractal dimension. Consequently, the surface will behave like a self-affine surface with Hurst exponent $`h_m`$ on lenghtscales comparable to the system size.

## 5 Conclusions

Table 1 summarizes our results. Model A without desorption reviews the results in , which have been obtained with slightly different activation energies on smaller systems and shorter timescales. Models B and C show, that desorption is an important process, which, although it affects only a small fraction of the adsorbed particles, must not be neglected, since it alters the scaling properties of the surfaces by reducing $`\beta `$ and by shifting the singularity spectrum towards smaller Hölder exponents. Since the scaling behaviour depends strongly on the height of the energy barrier of desorption, and the singularity spectra have no measurable tendency to narrow with time, our results can not be used to make any decision on the aymptotic universality class of the investigated model. However, they show, that the paradigm of a few universality classes characterized by a small number of exponents, which are independent on details of the model, is not adequate to catch the features of kinetic roughening on experimentally relevant timescales of a few hours of growth.

We are convinced, that the application of new mathematical tools like the wavelet analysis will help to find a better description of fractal growth phenomena in the future.

We thank A. Arnéodo and J. M. López for providing us recent preprints before publication and A. Freking for a critical reading of the manuscript.