text:
# A variational approach to nonlinear dynamics of nanoscale surface modulations

## 1 Introduction

In the recent years, nanoscale material structures have been fabricated on solid surfaces using various patterning techniques which, in some cases, exploit the phenomenon of stress-driven self-assembly. Potential application of these structures as functional elements in nano- and micro-electronics has led to considerable interest in understanding their stability and the processes that play a role in determining their morphological evolution. Due to their small size, the time scales involved in the processing and fabrication of these structures are small. Their small size also implies that the surface energy is a significant part of the overall energy of individual structures and, hence, is important in determining their morphological evolution. For example, transport of mass via surface diffusion can lead to smaller surface areas, thereby reducing surface energy for orientation independent surface energy density. Surface diffusion can also result in transformation of surface orientations with large surface energies to certain preferred low-energy orientations. Therefore, a fundamental understanding of the kinetics of surface transport, including the role of orientation-dependent surface energy, is important in the fabrication of nanostructures.

In the case of patterned features on the surfaces of amorphous solids or on crystalline surfaces above the thermodynamic roughening temperature, the surface energy is independent of the surface orientation. In either situation, surface diffusion is driven solely by the tendency to reduce the surface area and is described within the framework established through the classic work of Herring  and Mullins . According to their theory, a small sinusoidal perturbation on a flat surface retains its shape while the amplitude decays at a rate that is determined by the wavelength of the perturbation. In the case of a crystal below its thermodynamic roughening temperature, on the other hand, the situation is more complicated. Because a surface with an orientation that is close to a high symmetry orientation is made up of steps, the dependence of surface energy density on orientation acquires a cusp at such a high symmetry orientation.

Experiments have shown that sinusoidal perturbations about these high symmetry surface orientations do not retain the sinusoidal shape. Instead, the peaks in the profile become flattened, with the extent of the flat tops increasing as the amplitude of the profile decays \[3-10\]. The strong nonlinearity associated with the cusp in the surface energy causes difficulties in extending the theory of Mullins to faceted crystals. In this case, the decay of periodic surface profiles has been theoretically studied using various continuum models \[3-5,11-16\], by explicitly accounting for the dynamics of individual steps  and by using Monte-Carlo simulations \[19-25\]. While useful insights have been gained, each method has its shortcomings and the problem of decay remains open \[26-28\]. In order to motivate the present work, which is based on a variational framework for handling the nonlinear and singular aspects of surface evolution of crystals below the roughening temperature, we briefly summarize here the key results obtained in these studies.

The cusp in the surface energy leads to partial differential equations for surface evolution that involve Dirac-delta functions or other nonanalytic functions of the surface slope. In the case of both unidirectional and bidirectional sinusoidal perturbations, these equations have been solved by smoothing out the singularities \[3-5,11\]. Even though calculations performed using this procedure show flattening of profile shapes due to facet formation, in agreement experimental observations, the approach is ad hoc and it is not clear that different smoothing schemes would lead to identical results. For unidirectional modulations, an alternate approach was taken by Villain, Zangwill and coworkers \[12-14\], who argued that the surface chemical potential for surface diffusion should not involve the singular term that arises from the cusp. Retaining only the nonsingular contributions arising from step interactions, they showed that the sinusoidal profile does not flatten as it decays but, to the contrary, the top of the profile becomes sharper with time. According to their model, the curvature of the profile diverges at its extremum points. In distinct contrast to the above approaches, Spohn developed an ingenious method to deal with the singular nature of the evolution equations for unidirectional modulations by transforming them to a free-boundary problem . In his approach, if the location of the facet is known in advance (which is true for highly symmetric profiles), an ordinary differential equation for the length of the facet can be coupled to nonsingular partial differential equations for the evolution of stepped regions to obtain a coupled problem that is well-posed. Numerical solution of these equations in the case of unidirectional profiles shows pronounced flattening of the tops and bottoms of the periodic profiles.

In order to resolve the inconsistent predictions of the various continuum theories, more atomistic based approaches, particularly Monte Carlo (MC) simulations \[19-25\] and models that incorporate the dynamics of straight steps , have been pursued by several investigators. While some of the MC simulations show greater tendency to formation of facets than do others, none of the simulations indicate a sharpening of the profile as predicted in  and . Since kinetics of surface diffusion within the MC simulations is slow, most of these simulations have been confined to small sample sizes which are sensitive to step fluctuations. Furthermore, most of the MC simulations do not account for the long-range interactions of surface defects and, therefore, cannot capture the inverse-square elastic or electrostatic interactions between steps. This problem does not occur in step dynamics models where step interactions can be included explicitly. However, this approach also has limitations, since the decay kinetics and profile shapes depend strongly on the details of the short-range step interactions which are poorly understood. For example, if an attractive interaction is present between unlike steps at the top of the sinusoidal profile, the profile tends to be flatter and decays more rapidly . Also, in the case of bidirectional profiles, step dynamics models become very complicated and tractable solutions have been achieved only for cylindrically symmetric surface perturbations that involve circular steps ; only approximate analytic results for the case of surface evolution due to evaporation kinetics have been obtained for bidirectional sinusoidal perturbations . The more complicated and experimentally relevant case of surface diffusion limited decay of sinusoidal profiles has not been treated in a quantitative way using discrete step-dynamics models.

In this paper, we propose a variational approach to solve the nonlinear evolution equations for surface evolution and address the problem of decay of both one- and two-dimensional sinusoidal profiles on crystals below the roughening temperature. In the past, variational approaches have been used to study surface evolution involving complex geometries and shapes, but attention was limited to functional forms of surface energy that are well-behaved. For example, a variational model for evolution of a surface with an isotropic surface energy was introduced by Needleman and Rice  to study the problem of grain boundary transport and cavitation. In recent years, Suo  has applied the variational formulation to consider surface evolution during motion of defects such as cracks, voids and inclusions. The formulation of the variational approach adopted here is similar, and attention is focused on addressing the difficulties posed by nonanalytic forms of the surface energy. The basic idea of the formulation is to expand the surface shape in terms of a complete set of basis functions and to use a variational principle derived from the continuum evolution equations to obtain coupled nonlinear ordinary differential equations for the expansion coefficients. The advantage of this method is that even though the evolution equations are singular, the ordinary differential equations that determine the expansion coefficients are well-behaved and can be integrated by standard numerical techniques. The only inputs required to solve them are the orientation-dependent surface energies and the surface diffusion constants; there are no other ad hoc parameters in the model. The conclusions drawn through application of this approach are briefly summarized here and the details are described in the sections to follow.

In general, the main strength of the variational approach is that it provides a basis for extracting approximate solutions of a set of equations representing a physical phenomenon without the need for a priori choices on which features of those equations to retain or discard. The approximation made in this approach is usually the introduction of a lower limit on resolution of features of the solution, so the approximate solution found is the ”best” available within this restricted range of possible solutions. The lower limit on resolution can be reduced indefinitely, in principle, and convergence is expected. However, convergence can be proved in only the simplest circumstances, and we are normally content with establishing essential features of behavior, with the physics of the system providing guidance on the level of detail needed. As will be demonstrated below, extensive numerical testing with ever finer resolution is a good indicator of the level of detail required in any particular problem.

In the case of decay of unidirectional sinusoidal profiles, when the step creation energy is finite, our results agree with the predictions of Spohn . Here, we find that the decay rate of the amplitude increases with the magnitude of step formation energy. On setting the step creation energy to zero (or if the singular term in the chemical potential is ignored), we find perfect agreement with the results of Villain and Zangwill . This implies that the term arising from the cusp in the surface energy represents a singular perturbation; the profile decay in the presence of the this term, however small, is qualitatively different from the decay when this term is absent. While our model is in agreement with the free-boundary approach for symmetric profiles, it has several advantages from both computational and conceptual viewpoints. For example, the free boundary approach requires the knowledge of the location of facets at the beginning of a calculation, but this is not known for the case of asymmetric profiles. Within the variational formulation, this information is not needed a priori and the approach can therefore be applied to study the evolution of arbitrary profile shapes. The numerical procedure for solving the free-boundary problem involves both partial and ordinary differential equations, and these solutions must be matched through certain continuity conditions at the edges of the facet. The present variational formalism involves only solution of coupled ordinary differential equations and is therefore more straightforward in its implementation. Perhaps the most significant advantage of the variational approach is in handling the evolution of two-dimensional profiles, where the formulation of the free-boundary problem becomes mathematically very involved. We are not aware of any calculations of two-dimensional decay of sinusoidal profiles, other than by regularizing the nonanalytic features in the surface energy. Once again, the present approach does not require ad hoc parameters and is not fundamentally different from the one-dimensional case.

The ability to handle two-dimensional profiles allows us to address an important and outstanding issue concerning the limiting behavior of highly elongated two dimensional profiles, that is, profiles with a wavelength of modulation in one of the coordinate directions that is much larger than the wavelength in the other direction. The question we would like to answer is whether the one-dimensional decay behavior is recovered as one of the modulating wavelengths of the two-dimensional profiles diverges, or whether one-dimensional decay has qualitatively distinct features that cannot be described as a limiting case of the two-dimensional behavior. It has been suggested by Rettori and Villian  that the later case is true, who point out the facets form in the case of two dimensional profiles because of the line tension contribution of steps in the surface chemical potential. Since there is no contribution of the step formation energy in their chemical potential for the 1D profiles, they conclude that there is no facet formation in the 1D case. Within the present approach, where the line tension contribution of the steps is explicitly accounted for, we are able to analytically establish that the limiting behavior of the 2D profiles approach the 1D case. In particular, we show that the evolution equations of the variational parameters that express the shape of the 2D profiles approach the corresponding equations for the 1D profiles as the aspect ratio of the 2D profiles (the ratio of wavelengths in the coordinate directions) approaches infinity. Our analytical results are confirmed by numerical calculations which show that the 2D profile decay is already close to the 1D behavior when the aspect ratio is between 100 and 1000.

In most experimental situations, the study of decay of the 1D sinusoidal perturbations is influenced by the presence of widely spaced steps that run in a direction perpendicular to the direction of the corrugation. These steps arise as a result of a very small miscut from the singular orientation and are difficult to eliminate111STM images that show the presence of these steps on patterned Si and Au surfaces are given in .. In the present work we will show that, as the miscut angle tends to zero, the decay profiles approach the shape of the sinusoidal 1D profile on the singular surface. In contrast to earlier expectations , our studies indicate that the tendency to form facets persists even as the miscut slope goes to zero.

This paper is organized in the following away. Sec. 2 is devoted to the development of the variational formulation for both one and two-dimensional profiles. Here, we will derive evolution equations in the case where mass transport is controlled by terrace-diffusion. Analytical calculations establishing the limiting one-dimensional behavior of highly elongated two dimensional profiles will also be presented in this section. In Sec. 3, we will present numerical results for the decay of one-dimensional and two-dimensional profiles and profiles on surfaces that are miscut with respect to a singular surface. Particular emphasis is given to the analysis of the decay rate of two-dimensional profiles as a function of their aspect ratios and the behavior of sinusoidal profiles on miscut surfaces as the angle of miscut vanishes. In the case of one-dimensional profiles, we will also consider the decay of asymmetric profile shapes and present numerical results on the convergence properties of the variational formulation. Some the predictions of our calculations are compared with recent experiments . A summary of key results is given in Sec. 5.

## 2 Variational Formulation

We start by considering the free energy $`F`$ of the surface $`S`$, written in terms of the orientation-dependent surface energy $`\gamma `$ as

$$F=_S\gamma \left(_Sh\right)𝑑S,$$
(1)
where $`h`$ denotes the height of the surface measured relative to a high symmetry reference plane and $`_S`$ is the interior gradient operator on the surface. Using $`v_n`$ to denote the normal velocity of the surface relative to the material points instantaneously on that surface, the rate of change of the total free energy can be expressed in terms of the surface chemical potential $`\mu `$ as

$$\dot{F}=_S\mu v_n𝑑S.$$
(2)
The above equation provides a definition for the surface chemical potential, which is a surface field representing the driving force for change of shape of the crystal by diffusive mass transport over its surface. The chemical potential can be viewed as a continuum representation of the sensitivity of the free energy of the system to alterations in surface shape. The mass flux on the surface, here denoted by the surface vector field $`𝐣`$, is assumed to be related to the surface chemical potential through a kinetic relationship of the form

$$𝐣=C_S\mu ,$$
(3)
where $`C`$ is a positive mobility parameter that, among other things, can depend on the local surface slope, temperature and the concentration of the diffusing species.222Here, for convenience, we assume that surface diffusion is isotropic, so that the kinetic relationship is written using only one mobility parameter. The variational formulation that will be developed in this section can be generalized to include diffusion anisotropy in a straightforward manner. Conservation of mass at each point on the surface requires that the net rate of accumulation of mass at that point on the surface, determined by the normal speed $`v_n`$, is the negative of the divergence of the surface mass flux at that point, or

$$v_n=_S𝐣.$$
(4)
When the surface energy is sufficiently smooth, eqns (1-4) can be solved by standard numerical techniques to obtain the evolution of the surface shape. Alternatively, the weak form of the equations can be solved by employing a finite-element formulation or other variational approaches. However, if the dependence of a surface energy on surface shape is not analytic, a direct solution of the evolution equations becomes difficult due to the singularities in the chemical potential. In such cases, variational approaches provide a powerful tool for determining the evolution of the surface.

The development of the variational formulation involves two key ingredients. First, suppose that $`𝐪`$ is any possible mass flux field, sharing properties of periodicity and smoothness with $`𝐣`$, but is otherwise arbitrary. If we take the inner product of each side of (3) with $`𝐪`$ and integrate each side over $`S`$, after the application of a vector identity, the result can be written in the form

$$_S_S(\mu 𝐪)𝑑S_S\mu _S𝐪𝑑S+_S\frac{1}{C}𝐣𝐪𝑑S=0.$$
(5)
If the fields and the surface shape are spatially periodic and if $`S`$ includes exactly one period of the system, then the first term in the above equation vanishes. This follows from the fact that $`\mu `$ would have the same value on opposite sides of a segment of $`S`$ due to periodicity, and the mass flux out of $`S`$ would exactly cancel mass flux out on the opposite side for the same reason. This result provides the basis for an extremum principle which can be expressed as a maximum principle or a minimum principle; the two forms are equivalent and we adopt the latter form here.

Next, define a functional over the range of admissible trial fields $`𝐣^{}`$ as

$$\mathrm{\Phi }[𝐣^{}]=_S\left[\mu _S𝐣^{}+\frac{1}{2C}𝐣^{}𝐣^{}\right]𝑑S.$$
(6)
Then (5) implies that the functional is stationary under arbitrary variations of the surface flux $`𝐣^{}`$ when the flux field is the actual flux field, that is, when $`𝐣^{}=𝐣`$. This follows immediately from (5) if $`𝐪`$ is interpreted as $`\delta 𝐣^{}`$ in the usual variational notation. The arbitrariness of $`\delta 𝐣^{}`$ and the fundamental theorem of the calculus of variations require that $`𝐣`$ must satisfy (3) pointwise on $`S`$. The fact that $`\mathrm{\Phi }[𝐣^{}]`$ is not only stationary but also an absolute minimum under variations about $`𝐣`$ can be seen by considering $`\mathrm{\Phi }[𝐣^{}+\delta 𝐣^{}]\mathrm{\Phi }[𝐣^{}]`$ to second order in $`\delta 𝐣^{}`$. Since the mobility parameter $`C`$ is positive, it follows from (6) that

$$\mathrm{\Phi }[𝐣^{}+\delta 𝐣^{}]\mathrm{\Phi }[𝐣^{}]=_S\frac{1}{2C}\delta 𝐣^{}\delta 𝐣^{}𝑑S>0,$$
(7)
thus ensuring that $`𝐣^{}=𝐣`$ is indeed a minimum of $`\mathrm{\Phi }[𝐣^{}]`$.

The main advantage of the variational formulation is that instead of solving the conventional evolution equations, we can expand the mass flux in terms of a complete set of basis functions and obtain the time dependence of the expansion coefficients by minimizing the functional $`\mathrm{\Phi }`$. There is no need to deal explicitly with the singular chemical potential while carrying out this minimization process. This can be seen by considering the first term in (6) which can be seen to be equal to $`\dot{F}`$ if we replace the normal velocity $`v_n`$ in favor of $`_S𝐣^{}`$ in (4); $`\dot{F}`$ can be directly evaluated in terms of $`𝐣^{}`$ without using the chemical potential. The variational functional can then be conveniently expressed as

$$\mathrm{\Phi }[𝐣^{}]=\dot{F}[𝐣^{}]+_S\frac{1}{2C}𝐣^{}𝐣^{}𝑑S.$$
(8)
Below, we provide details of the implementation of the variational formulation to study the decay of periodic one and two dimensional profiles; the one-dimensional case is considered first.

### 2.1 One-dimensional profiles

The orientation-dependent surface energy that applies to one-dimensional modulations can be written as

$$\gamma (h_x)=\beta _1\left|h_x\right|+\frac{\beta _3}{3}\left|h_x\right|^3,$$
(9)
where the subscript $`x`$ denotes the derivative with respect to the spatial coordinate and the terms proportional to $`\beta _1`$ and $`\beta _3`$ represent step formation and step interaction energies, respectively. This form of the surface energy expression presumes that $`|h_x|1`$. If we focus attention on periodic profiles with wavelength $`\lambda `$, then free energy in one period and the position dependent surface chemical potential can be written as

$$F=_0^\lambda \gamma (h_x)𝑑x,$$
(10)
and

$$\mu (x)=\frac{d}{dx}\left(\frac{d\gamma }{dh_x}\right)=\beta _1\delta \left[h_x\right]h_{xx}2\beta _3\left|h_x\right|h_{xx},$$
(11)
respectively, where $`\delta []`$ denotes the Dirac delta function. The singular term in the chemical potential can be ignored if the curvature vanishes at the points where the slope goes to zero. However, since the step interaction term is cubic in the surface slope, it can be shown that the curvature diverges like $`(xx_0)^{1/2}`$ around an extremum located at $`x=x_0`$ . This implies that the first term plays a crucial role in determining the decay shape.

If the surface shape is symmetric about the origin, it can be expanded in terms of a Fourier cosine series as

$$h(x,t)=\underset{n=1}{\overset{N}{}}a_n(t)\mathrm{cos}(nkx),\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}x\lambda ,$$
(12)
where $`k=2\pi /\lambda `$ and the number of Fourier modes $`N`$ is finite but otherwise unrestricted.333Numerical studies on the sensitivity of the results to the number of Fourier coefficients and the convergence properties of the present formulation will be presented in Sec. 3. If attention is focused on small perturbations that satisfy the criterion $`|h_x|1`$, the normal velocity of the surface can be written in terms of surface shape as $`v_n=h_t`$, where subscript $`t`$ denotes differentiation with respect to time. We can now invoke conservation of mass (4) to express the surface mass flux as

$$j(x,t)=\underset{n=1}{\overset{N}{}}\dot{a}_n(t)\frac{\mathrm{sin}(nkx)}{nk}.$$
(13)
The variational formulation provides first-order coupled ordinary differential equations for the Fourier coefficients $`a_n(t)`$; these equations will be numerically integrated to obtain the evolving surface shape for $`t>0`$ using information on the surface shape at $`t=0`$.

In order to write the functional $`\mathrm{\Phi }`$ (refer to (8)) in terms of the variational parameters, which in the present case are the $`\dot{a}_n`$s, we first have to specify the dependence of the surface mobility parameter on the surface slope. This dependence can be obtained by considering the mechanisms involved in determining mass transport on stepped surfaces, namely, the attachment/detachment of atoms from step-edges and the diffusion of adatoms on the terraces separating the steps. In this paper we will restrict attention to the case where the mass transport is limited by terrace diffusion; the mobility parameter in this situation does not depend the surface slope.444In case of attachment-detachment kinetics, a slope dependent mobility parameter  can be used to derive the evolution equations. These equations will be reported elsewhere.

Using the expression for mass flux in (13), we find the variational functional for the terrace-diffusion limited case to be

$$\mathrm{\Phi }[\dot{a}_1,\mathrm{},\dot{a}_N]=\underset{n=1}{\overset{N}{}}H_n[a_1,\mathrm{},a_N]\dot{a}_n(t)+\frac{\pi }{2k^3C}\underset{n=1}{\overset{N}{}}\frac{\dot{a}_n^2}{n^2},$$
(14)
where

$$H_n[a_1,\mathrm{},a_N]=nk\beta _1_0^\lambda \mathrm{Sgn}\left[h_x\right]\mathrm{sin}(nkx)𝑑x+nk\beta _3_0^\lambda \mathrm{Sgn}\left[h_x\right]\left|h_x\right|^2\mathrm{sin}(nkx)𝑑x,$$
(15)
in which Sgn\[$`p`$\] $`=d|p|/dp`$. Minimizing the functional by setting $`\mathrm{\Phi }/\dot{a}_n=0`$ for $`n=1,\mathrm{},N`$, we obtain the evolution equations

$$\dot{a}_n(t)=\frac{n^2k^3C}{\pi }H_n[a_1,\mathrm{},a_N]$$
(16)
for the Fourier coefficients in (12). If the initial shape is known, approximate solutions to these coupled first-order ordinary differential equations can be found using standard numerical integration techniques.

If the problem under consideration is the decay of a single Fourier mode of amplitude $`a_0`$, the initial conditions satisfy

$$a_n(0)=\{\begin{array}{cc}a_0\hfill & n=1\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}$$
(17)
An important point to note is that the evolution equations are well-behaved and do not directly inherit any singularities that are present in the chemical potential. However, the influence of the physical features that give rise to the singularities are incorporated.

Insights into the nature of the solutions can be gained by expressing (16) in terms of dimensionless variables. Introducing the rescaled variables

$$\stackrel{~}{h}=kh,\stackrel{~}{a}_n=ka_n,\stackrel{~}{x}=kx,\stackrel{~}{t}=\frac{C\beta _3k^4}{\pi }t$$
(18)
we find the Fourier coefficients satisfy

$$a_n(t)=\frac{1}{k}\stackrel{~}{a}_n(\stackrel{~}{\beta }_1,\stackrel{~}{t})$$
(19)
where $`\stackrel{~}{\beta }_1=\beta _1/\beta _3`$ and the functional form of the solution $`\stackrel{~}{a}_n`$ must be determined by solving (16). We also point out that when the surface shape is not symmetric about the origin, a complete Fourier expansion that involves both cosine and sine terms should be adopted; the evolution equations for the expansion coefficients can be derived by following the steps outlined in this section.

### 2.2 Two-dimensional profiles

The derivation of evolution equations for two-dimensional surface modulations closely follows the procedure used in the one-dimensional case. In the case of symmetric profiles, the surface shape can be expressed in terms of the Fourier series

$$h(𝐱,t)=\underset{m=1}{\overset{N}{}}\underset{n=1}{\overset{N}{}}A_{mn}(t)\mathrm{cos}(mk_1x_1)\mathrm{cos}(nk_2x_2),$$
(20)
where $`k_{1(2)}=2\pi /\lambda _{1(2)}`$ is the wavenumber of the modulation in the $`x_1(x_2)`$-direction. If the amplitude of the modulation is small, we can invoke conservation of mass, expressed in terms of the surface shape as

$$\frac{h(𝐱,t)}{t}+𝐣(𝐱,t)=0,$$
(21)
to obtain the components of the surface mass flux

$$(j_1,j_2)=\underset{m}{}\underset{n}{}(\dot{a}_{mn}(t)\frac{\mathrm{sin}(mk_1x_1)\mathrm{cos}(nk_2x_2)}{mk_1},\dot{b}_{mn}(t)\frac{\mathrm{cos}(mk_1x_1)\mathrm{sin}(nk_2x_2)}{nk_2}),$$
(22)
where

$$\dot{a}_{mn}(t)+\dot{b}_{mn}(t)=\dot{A}_{mn}(t).$$
(23)
For crystal surfaces that are below the roughening temperature, the surface energy can be written as

$$\gamma (h_{x_1},h_{x_2})=\beta _1(h_{x_1}^2+h_{x_2}^2)^{1/2}+\frac{\beta _3}{3}(h_{x_1}^2+h_{x_2}^2)^{3/2},$$
(24)
where $`\beta _1`$ and $`\beta _3`$ have the same meaning as in (9). If we restrict attention to the terrace-diffusion limited surface kinetics, the mobility parameter in (3) becomes independent of surface slope. In this case, using (22) and (24), the variational functional can be written as

$$\mathrm{\Phi }[\dot{a}_{11},\dot{a}_{12},\mathrm{},\dot{b}_{11},\mathrm{}]=\underset{n=1}{\overset{N}{}}\underset{m=1}{\overset{N}{}}(\dot{a}_{mn}+\dot{b}_{mn})H_{mn}+\underset{n=1}{\overset{N}{}}\underset{m=1}{\overset{N}{}}\frac{\pi ^2}{2k_1k_2C}\left[\frac{\dot{a}_{mn}^2}{m^2k_1^2}+\frac{\dot{b}_{mn}^2}{n^2k_2^2}\right],$$
(25)
where

$$H_{mn}=_0^{\lambda _1}_0^{\lambda _2}\left[mk_1\frac{\gamma }{h_{x_1}}\mathrm{sin}(mk_1x_1)\mathrm{cos}(nk_2x_2)+nk_2\frac{\gamma }{h_{x_2}}\mathrm{cos}(mk_1x_1)\mathrm{sin}(nk_2x_2)\right]𝑑x_1𝑑x_2.$$
(26)
Minimizing the functional with respect to $`\dot{a}_{mn}`$ and $`\dot{b}_{mn}`$, we obtain the first-order ordinary differential equations

$$\dot{a}_{mn}(t)=\frac{Cm^2k_1^3k_2}{\pi ^2}H_{mn},\dot{b}_{mn}(t)=\frac{Cn^2k_1k_2^3}{\pi ^2}H_{mn}.$$
(27)
Adding these two contributions according to (23), we find that the evolution equation for the Fourier amplitude $`A_{mn}`$ in (20) is

$$\dot{A}_{mn}(t)=\frac{C(m^2k_1^3k_2+n^2k_1k_2^3)}{\pi ^2}H_{mn}.$$
(28)
If the amplitudes of the Fourier modes at $`t=0`$ are known, (28) can be integrated to obtain the surface shape for $`t>0`$.

The dependence of the decay rate on the aspect ratio of the profile $`\alpha \lambda _1/\lambda _2`$ can be understood by using the following dimensionless variables:

$$\stackrel{~}{x}_{1(2)}=k_{1(2)}x_{1(2)},\stackrel{~}{h}=k_1h,\stackrel{~}{A}_{mn}=k_1A_{mn},\stackrel{~}{t}=\frac{C_1\beta _3k_1^4}{\pi ^2}t,\stackrel{~}{H}_{mn}=\frac{k_2}{\beta _3}H_{mn}.$$
(29)
The evolution equations for the scaled Fourier amplitudes in terms of these rescaled variables are

$$\frac{d\stackrel{~}{A}_{mn}}{d\stackrel{~}{t}}=(m^2+\alpha ^2n^2)\stackrel{~}{H}_{mn},$$
(30)
where

$$\stackrel{~}{H}_{mn}=_0^{2\pi }_0^{2\pi }\frac{(\stackrel{~}{\beta }_1+\stackrel{~}{h}_{\stackrel{~}{x}_1}^2+\alpha ^2\stackrel{~}{h}_{\stackrel{~}{x}_2}^2)}{\sqrt{\stackrel{~}{h}_{\stackrel{~}{x}_1}^2+\alpha ^2\stackrel{~}{h}_{\stackrel{~}{x}_2}^2}}\left[m\stackrel{~}{h}_{\stackrel{~}{x}_1}\mathrm{sin}(m\stackrel{~}{x}_1)\mathrm{cos}(n\stackrel{~}{x}_2)+n\alpha ^2\stackrel{~}{h}_{\stackrel{~}{x}_2}\mathrm{cos}(m\stackrel{~}{x}_1)\mathrm{sin}(n\stackrel{~}{x}_2)\right]𝑑\stackrel{~}{x}_1𝑑\stackrel{~}{x}_2.$$
(31)
As noted earlier, the profile decay for the case of $`\beta _1=0`$ is qualitatively different from the case when this parameter is non-zero but arbitrarily small. It is also of interest to examine the limiting behavior of the decay of the two-dimensional profiles as $`\alpha 0`$, which corresponds to the situation of one-dimensional behavior determined as a limit of two-dimensional behavior. To this end, we proceed to show that the two-dimensional evolution equations precisely approach the one-dimensional equations as $`\alpha 0`$.

### 2.3 Decay of elongated 2D profiles: approach of the 1D limit

The approach to the one-dimensional limit can be studied by looking at the decay behavior of $`h(x_1,0,t)`$, given by555Since the variational approach is strictly exact only when all the Fourier modes are retained, the upper limit in the sum over the Fourier modes will be taken to be infinity. This allows us to use closed form results for certain infinite series; these results are required for establishing the limiting 1D behavior.

$$h(x_1,0,t)=\underset{m}{}\left(\underset{n=1}{\overset{\mathrm{}}{}}A_{mn}(t)\right)\mathrm{cos}(mx_1)$$
(32)
when $`\alpha 0`$; the goal is to establish that the quantity $`A_n^{\mathrm{}}(t)_{m=1}^{\mathrm{}}A_{nm}(t)`$ obeys the same evolution equations as the one-dimensional expansion coefficients, $`a_n(t)`$ in (12). It follows from (30) and (31) that when $`\alpha 0`$, the scaled Fourier amplitudes satisfy

$$\frac{d\stackrel{~}{A}_m^{\mathrm{}}}{d\stackrel{~}{t}}=m^2\stackrel{~}{H}_m^{\mathrm{}},$$
(33)
where

$$\stackrel{~}{H}_m^{\mathrm{}}=m\underset{n=1}{\overset{\mathrm{}}{}}_0^{2\pi }_0^{2\pi }(\stackrel{~}{\beta }_1+\stackrel{~}{h}_{\stackrel{~}{x}_1}^2)\mathrm{Sgn}[\stackrel{~}{h}_{\stackrel{~}{x}_1}]\mathrm{sin}(m\stackrel{~}{x}_1)\mathrm{cos}(n\stackrel{~}{x}_2)𝑑\stackrel{~}{x}_1𝑑\stackrel{~}{x}_2.$$
(34)
The integral in (34) can be explicitly evaluated by noting that

$$\mathrm{Sgn}[\stackrel{~}{h}_{\stackrel{~}{x}_1}]=\mathrm{Sgn}[\stackrel{~}{h}_{\stackrel{~}{x}_1}^{\mathrm{}}]\mathrm{Sgn}[\mathrm{cos}(\stackrel{~}{x}_2)],$$
(35)
where $`\stackrel{~}{h}_{\stackrel{~}{x}_1}^{\mathrm{}}=_pp\stackrel{~}{A}_p^{\mathrm{}}\mathrm{sin}(p\stackrel{~}{x}_1)`$. In what follows we will separately calculate the contribution to $`\stackrel{~}{H}_m^{\mathrm{}}`$ due to terms that arise from step formation and interactions.

Using the results

$$_0^{2\pi }\mathrm{Sgn}[\mathrm{cos}(\stackrel{~}{x}_2)]\mathrm{cos}(n\stackrel{~}{x}_2)𝑑\stackrel{~}{x}_2=\{\begin{array}{cc}\frac{4}{n}\hfill & n=1,5,9,\mathrm{}\hfill \\ & \\ \frac{4}{n}\hfill & n=3,7,11,\mathrm{}\hfill \\ & \\ 0\hfill & n=2,4,6,\mathrm{}\hfill \end{array}$$
(36)
and

$$\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n}{(2n+1)}=\frac{\pi }{4},$$
(37)
the contribution to (34) due to the step formation term is shown to be

$$\stackrel{~}{H}_m^{\mathrm{}}|_f=\pi m\stackrel{~}{\beta }_1_0^{2\pi }\mathrm{Sgn}[\stackrel{~}{h}_{\stackrel{~}{x}_1}^{\mathrm{}}]\mathrm{sin}(m\stackrel{~}{x}_1)𝑑\stackrel{~}{x}_1.$$
(38)
The contributions due to step interactions can be evaluated by employing the expansion

$$\stackrel{~}{h}_{\stackrel{~}{x}_1}^2=\underset{p,q,r,s}{}pq\stackrel{~}{A}_{pr}\stackrel{~}{A}_{qs}\mathrm{sin}(p\stackrel{~}{x}_1)\mathrm{sin}(q\stackrel{~}{x}_1)\mathrm{cos}(r\stackrel{~}{x}_2)\mathrm{cos}(s\stackrel{~}{x}_2)$$
(39)
and the result

$$\underset{l=1}{\overset{\mathrm{}}{}}_0^{2\pi }\mathrm{Sgn}[\mathrm{cos}(\stackrel{~}{x}_2)]\mathrm{cos}(l\stackrel{~}{x}_2)\mathrm{cos}(m\stackrel{~}{x}_2)\mathrm{cos}(n\stackrel{~}{x}_2)𝑑\stackrel{~}{x}_2=\pi ,$$
(40)
which can be proved for any integers $`m`$ and $`n`$ by invoking certain trigonometric identities. Using the fact that the sum in (40) is independent of both $`m`$ and $`n`$, it can be shown that the contribution to (34) due to the step interactions is

$$\stackrel{~}{H}_m^{\mathrm{}}|_i=\pi m_0^{2\pi }(\stackrel{~}{h}_{\stackrel{~}{x}_1}^{\mathrm{}})^2\mathrm{Sgn}[\stackrel{~}{h}_{\stackrel{~}{x}_1}^{\mathrm{}}]\mathrm{sin}(m\stackrel{~}{x}_1)𝑑\stackrel{~}{x}_1.$$
(41)
Combining (38) and (41) according to (34) , we find that

$$\frac{dA_n^{\mathrm{}}}{dt}=\frac{C_1\beta _3n^3k_1^4}{\pi }_0^{\lambda _1}\left(\stackrel{~}{\beta }_1+(h_{x_1}^{\mathrm{}})^2\right)\mathrm{Sgn}[h_{x_1}^{\mathrm{}}]\mathrm{sin}(nx_1)𝑑x_1,$$
(42)
which shows that the coefficients $`A_n^{\mathrm{}}`$ satisfy the same evolution equations as the 1-dimensional expansion coefficients (refer to (16)) if we identify $`C_1`$, $`k_1`$ and $`\lambda _1`$ with $`C`$, $`k`$ and $`\lambda `$ respectively. Thus, in the limit of vanishing wavelength aspect ratio, the two-dimensional case reduces to the independently established one-dimensional case. Numerical results for profile decay as a function of the aspect ratio are discussed in the following section.

### 2.4 Sinusoidal profiles on vicinal surfaces

As noted earlier, experiments on 1D sinusoidal modulations are always affected by the presence of steps due to small miscuts in directions in a direction perpendicular to the direction of the surface modulations. In this case of a modulated vicinal surface, the surface height can be expressed as

$$h(x_1,x_2)=\underset{n=1}{\overset{N}{}}a_n(t)\mathrm{cos}(nkx_1)s_0x_2,$$
(43)
where $`\mathrm{tan}^1(s_0)`$ is the miscut angle. At time $`t=0`$, we assume that all the Fourier coefficients except the longest wavelength mode ($`n=1`$)have vanishing amplitudes. A schematic of this initial surface shape is given in Fig. 1. We will now present the evolution equations for the Fourier coefficients and discuss their behavior as $`s_00`$.

Using the procedure used to derive the evolution equations for the 2D profiles, we find the Fourier coefficients satisfy the evolution equations

$$\dot{a}_n(t)=\frac{n^2k^3C}{\pi }H_n[a_1,\mathrm{},a_N],$$
(44)
where

$$H_n[a_1,\mathrm{},a_N]=nk_0^\lambda \frac{h_{x_1}}{\sqrt{h_{x_1}^2+s_0^2}}(\beta _1+\beta _3h_{x_1}^2)\mathrm{sin}(nkx_1)𝑑x_1.$$
(45)
It can be seen that these equations have the same form as the evolution equations of the Fourier coefficients for the 1D profiles given by (16), except for the fact that the $`\mathrm{Sgn}`$ function of the surface slope in the expression for $`H_n`$ is replaced by an analytic function in (45). Furthermore, since

$$\underset{s_00}{lim}\frac{h_{x_1}}{\sqrt{h_{x_1}^2+s_0^2}}=\mathrm{Sgn}[h_{x_1}],$$
(46)
the evolution equations for the Fourier coefficients in (43) approach the corresponding equations for the for the 1D modulations derived earlier. We therefore expect the profile shapes on the miscut surfaces to coincide with the 1D decay profiles of the as the angle of miscut vanishes. The details of the numerical calculations that allow us to look at this limiting behavior are given in the following section.

## 3 Numerical Results

In this section, we present results on the decay toward morphological equilibrium of surface modulations obtained by numerically integrating the evolution equations derived in Sec. 2. We will first consider the decay of one dimensional sinusoidal perturbations for different values of step-formation and step-interaction energies. Numerical studies on the sensitivity of the results to the number of Fourier components that are used to express the surface shape will then be presented. This is followed by an example of the decay of an anisotropic profile where the location of the facets are not easily determined at the beginning of the calculation. In the particular case that we consider, some of the facets that appear during the early stages of profile evolution, eventually disappear as other neighboring facets grow at faster rates. Next, we consider the decay rates of two-dimensional profiles with different aspect ratios and diffusion constants. Finally, we focus on the evolution of sinusoidal modulations on a vicinal surface that is slightly misoriented with respect to the flat singular surface.

### 3.1 1D sinusoidal profiles

The evolution equations given by (16) were integrated using a Runge-Kutta integrator with adaptive step-size control discussed in detail in Ref. . The time dependence of the amplitude decay, for different values of rescaled step-formation energies ($`\stackrel{~}{\beta }_1`$), is shown in Fig. 2. The profile shapes at different stages of morphological equilibration are shown, for certain representative values of $`\stackrel{~}{\beta }_1`$, in Fig. 3. When the step formation energy vanishes, we find that the profile sharpens at its extrema and the decay rate shows the $`1/\stackrel{~}{t}`$ scaling behavior predicted by Ozdemir and Zangwill  (refer to inset (a) in Fig. 2). On the other hand, for any non-zero value of $`\stackrel{~}{\beta }_1`$, we find that the profile decays with the formation of facets which, after an initial transient period of growth, remain nearly constant in size as seen in Fig. 3. The amplitude during this self-similar decay stage, decreases as a linear function of time, (refer to Fig. 2), in agreement with the results of Hager and Spohn  which were obtained using the free-boundary approach. With increasing value of $`\beta _1`$, the lengths of the facets in the self-similar shape increase in magnitude and the initial transient periods decrease in their duration. These results indicate that the case $`\stackrel{~}{\beta }_10`$ is a singular perturbation of the special case with $`\stackrel{~}{\beta }_1=0`$; any non-zero value of $`\stackrel{~}{\beta }_1`$ leads to facet formation during the morphological equilibration.

The sensitivity of the profile shapes and decay rates on the number of Fourier modes used in the expansion of the surface shape (refer to (12)) is considered in Fig. 4, where the decay rates obtained by retaining 15, 20 and 30 modes are given. The insets in the figure show the evolution of the amplitudes of the individual Fourier modes for the three cases considered. We find that the amplitude decay remains virtually unchanged as the number of Fourier modes increase from 15 to 30. Fig. 4 shows that the dominant Fourier components in the surface shape are the long wavelength modes and that the magnitudes of the different Fourier components are also seen to be independent of the total number of modes that are used in the analysis. The insets also show that all the Fourier components decay linearly in time after an initial transient period; in this regime, the surface shape remains self-similar as noted earlier. We have investigated the time dependence of the decay shapes for several other values of the step-formation energies and initial amplitudes of the sinusoidal profiles and found that accurate results are obtained by using between 15 and 20 Fourier modes. As is the case with any type of spectral method, if more basis functions are retained, the time steps used in numerical integration have to be decreased, which leads to longer simulation times. The optimum choice for a particular application is determined by a trade-off between accuracy and simulation time.

### 3.2 Asymmetric profiles

Next, we consider the evolution of an asymmetric profile shown in Fig. 5(a), where the surface shape can be expressed by using both cosine and sine functions in the Fourier expansion as

$$h(x,t)=\underset{m}{}a_m(t)\mathrm{cos}(nkx)+\underset{n}{}b_n(t)\mathrm{sin}(nkx).$$
(47)
The amplitudes of the non-vanishing Fourier components in the initial shape were chosen as follows: $`ka_1=0.2,ka_3=0.05,kb_2=0.1`$ and $`kb_4=0.025`$. We find that during the early stages of morphological evolution, facets labeled 1, 2 and 3 appear as shown in Fig. 5(b). While these facets are situated close to the extrema in the initial profile, their actual locations are not determined very easily. As the shape evolves in time, the facet labeled 3 shrinks in size, while the other two facets grow is size; in Fig. 5(c) this facet has almost disappeared. In the late stages of the evolution, all the Fourier sine components with short wavelengths vanish, so that the shape closely resembles the decay profile of the symmetric cosine mode with the largest wavelength.

The above example illustrates that the present approach to surface evolution can be used to study arbitrary profile shapes in a straightforward manner. While the free-boundary approach, in principle, can be used in this case, the actual implementation can be quite tedious. First, the locations of the facets are difficult to guess at the outset and secondly, application of boundary conditions at the points at which the facets disappear require special consideration. On the other hand, there is no fundamental difference in the analysis of symmetric and asymmetric profiles within the variational approach.

### 3.3 2D sinusoidal profiles

In the case of bidirectional modulations, we start by looking at the decay behavior of a bidirectional profile whose wavelengths of modulation in the coordinate directions are equal. The amplitude evolution and profile shapes for different values of $`\stackrel{~}{\beta }_1`$ are given in Fig. 6 and Fig. 7, respectively. Just as in the 1D case, when $`\stackrel{~}{\beta }_1=0`$, we find that the profile decays without formation of facets. Furthermore, the amplitude of the modulation decays with the $`1/\stackrel{~}{t}`$ scaling behavior, which is also observed in the decay of 1D profiles. It is also evident from Fig. 7 that the surface profile close to extrema becomes shaper than the parabolic shape, which once again resembles the behavior in 1D. For any non-zero $`\stackrel{~}{\beta }_1`$, the profile decays with the formation of facets, with the decay rate increasing with $`\stackrel{~}{\beta }_1`$. The decay rate, however, does not follow the $`1/\stackrel{~}{t}`$ scaling behavior; instead, it shows a linear behavior after an initial transient as indicated in Fig. 7; this trend is also similar to the decay of a 1D profile. These results indicate that there is no fundamental difference in the profile decay behavior in one- and two- dimensions; the decay behavior is governed primarily by the ratio of the step formation and interaction energies.

In Sec. 2, we presented analytical arguments that establish that the decay of elongated the 2D profiles approaches the 1D limit as the aspect ratio of the profile diverges. We now discuss numerical calculations that allow us to understand the way the profile shapes evolve with increasing aspect ratios. In Fig. 9, the decay rates for profiles with different values of $`\alpha `$ are given along with the curve for the case $`\alpha =0`$, calculated with the 1D evolution equations (refer to (16)). We find that as $`\alpha `$ falls below 1, the decay rates initially decrease continuously until $`\alpha `$ reaches about $`0.001`$. As $`\alpha `$ decreases further, there is no significant change in the profile decay, so that in the limit of very small $`\alpha `$ all the decay curves converge to a single curve. It can be seem from Fig. 9 that this curve coincides with the decay profile of the 1D case, in agreement with the analytical results of Sec. 2. We have also looked at the convergence properties for several other values of $`\stackrel{~}{\beta }_1`$ and initial amplitudes and have found that the 1D limit is achieved when $`\alpha `$ lies between 0.01 and 0.001. The evolution of the surface shape during the decay of an elongated profile with $`\alpha =0.1`$ is illustrated in Fig. 8. It can be seen that well-defined straight features that are aligned along the direction in with the profile is elongated, appear in the decay shape. The cross-section of the profile, particularly at the extrema of the elongated shape, are found to closely resemble the decay profile of a 1D modulation.

### 3.4 Sinusoidal profiles on vicinal surfaces

As a final example of the application of the variational approach, we consider the decay of sinusoidal profiles on vicinal surfaces that are slightly misoriented from the flat singular surface. Our goal is to understand the limiting behavior of profile decay as the miscut slope of the surface approaches zero. To this end, we consider the decay of a sinusoidal modulation with scaled amplitude, $`ka_1=0.1`$ in Fig. 10. The step formation energy was taken to be $`\stackrel{~}{\beta }_1=1.0`$. The figure shows decay curves for several values of the misorientation slopes, along with the decay profile for the 1D modulation calculated using (16). We find that the profile decay progresses at a faster rate as the misorientation slope decreases from 0.1 to 0.01. As the slope is decreased further, all the decay curves approach a limiting curve which coincides with the 1D decay profile. In the particular example that we consider, the decay behavior is already very close to the 1D case when the misorientation slope is about $`0.001`$. We have also numerically investigated the decay behavior for several other values of step formation energies and initial amplitudes and found that the decay curves approach the 1D case as the misorientation slope approaches zero.

Our results differ from the earlier results of Lancon and Villain , who point out that facets form on surfaces only if the misorientation slope is larger than some critical value. For the range of parameters considered in Fig. 10, they estimate the critical value to be $`s_0=0.001`$. In contrast to their predictions, our numerical studies clearly show that facets form for arbitrarily small values of the miscut slope. This observation is consistent with the fact that the evolution equations for the Fourier coefficients that determine the shape modulations on miscut surfaces approach the corresponding equations for 1D coefficients as the miscut slope vanishes.

### 3.5 Comparison to experiments on Si(001)

We now turn our attention to recent experimental work of Erlebacher and coworkers , who investigated the decay of sputter ripples on Si(001) surfaces. They found that amplitude of the ripples, with $`\alpha =\lambda _1/\lambda _2=0.1`$, exhibited a decay behavior consistent with the $`1/t`$ scaling behavior for 1D modulations predicted by Zangwill and Villain and coworkers . In this section, using available experimental and theoretical data on step formation and interaction energies for Si(001), we will numerically investigate the decay of these elongated ripples.

The parameters $`\beta _1`$ and $`\beta _3`$, for the Si(001) surfaces with alternating SA and SB steps relevant to the experiments, can be estimated from experimentally determined step formation energies reported by Zandvliet  and the atomistic calculations of Poon and coworkers . Using the step formation energies of SA (8 meV/Å) and SB (16 meV/Å) steps, we find that $`\beta _1`$ 12.5 meV/Å2. From the atomistic calculations of step interactions on Si(001) , $`\beta _3`$ is estimated to be $``$1200 meV/Å2. The profile decay for the elongated ripples with $`\stackrel{~}{\beta }_1=0.01`$ is shown in Fig. 11. As noted earlier, when $`\beta _1=0`$, both 1D and 2D profiles show a $`1/t`$ decay behavior. The figure shows that in distinct contrast to the $`\beta _1=0`$ case, even for a relatively small non-zero value of $`\beta _1`$, the amplitude decay does not exhibit the $`1/t`$ scaling behavior. The observed behavior is therefore not consistent with diffusion limited decay of 2D ripples. The results can perhaps be explained by including anisotropies in surface diffusion and the effects of force-monopoles at the steps due to the 2$`\times `$1 and 1$`\times `$2 reconstruction of the terraces between the steps. It has also been suggested that the surface transport at the temperatures of interest is limited by attachment-detachment kinetics . We leave a detailed consideration of these issues for future study.

## 4 Concluding Observations

In conclusion, we have developed a numerical method to study the singular equations that govern the dynamics of crystal surfaces below the roughening temperature. The method was applied to study the morphological equilibration of both unidirectional and bidirectional modulations. In all these cases, the surface evolution was determined by solving well-behaved coupled nonlinear ordinary differential equations for the modal expansion coefficients of the surface shape. Analysis of highly elongated 2D sinusoidal profiles and sinusoidal profiles on vicinal surfaces with vanishing miscut angles shows that the limiting behavior in both of these cases agrees with the 1D decay behavior accompanied by facet formation. In the following paragraph, we point to directions of further work that will allow us to make closer connection to experiments.

In this paper we restricted attention to surface diffusion kinetics limited by terrace diffusion, while many experimental systems are in a regime where the mass transport is limited by attachment and detachment at the steps . For the latter case, Nozieres has derived a kinetic relation between the gradient of the surface chemical potential and the mass flux . In the future we plan to derive a variational formulation appropriate to this kinetic law to analyze morphological evolution in the attachment-detachment limited regime. Also, in some experimental situations, anisotropies associated with surface diffusion as well as surface energies become can be significant. As noted earlier, these effects can be included in the present formulation and will be the subject of forthcoming publications.

## Acknowledgments

The research support of the National Science Foundation through grants CMS-0093714 and CMS-0210095 and the Brown University MRSEC Program is gratefully acknowledged. VBS thanks the Graduate School at Brown University for providing computational support through the Salomon Research Award.

## References

1. C. Herring, The use of classical macroscopic concepts in surface energy problems, Structure and Properties of Solid Surfaces, edited by R. Gomer and C. S. Smith, University of Chicago Press, (1953) 5.
2. W. W. Mullins, J. Appl. Phys. 28 (1957) 333.
3. H. P. Bonzel, E. Preuss and B. Steffen, Appl. Phys. A 35 (1984) 1.
4. H. P. Bonzel, E. Preuss and B. Steffen, Surf. Sci. 145 (1984) 20.
5. H. P. Bonzel and E. Preuss, Surf. Sci. 336 (1995) 209.
6. C. C. Umbach, M. E. Keeffe and J. M. Blakely, J. Vac. Sci. Technol. A 9 (1991) 1014.
7. M. E. Keeffe, C. C. Umbach and J. M. Blakely, J. Phys.Chem. Solids 45 (1994) 965.
8. S. Tanaka, C. C. Umbach, J. M. Blakely, R. M. Tromp and M. Manos, in Structure and Evolution of Surfaces, Edited by R. C. Cammarata, E. H. Chason, T. L. Einstein and E. D. Williams, MRS Symposium Proceedings, Vol.440 (Materials Research Society, Pittsburg, 1997) p. 25.
9. S. Tanaka, C. C. Umbach, J. M. Blakely, R. M. Tromp and M. Manos, J. Vac. Sci. Technol. A 15 (1997) 1345.
10. J. Blakely, C. Umbach and S. Tanaka, Dynamics of Crystal Surfaces and Interfaces, Edited by P. M. Duxbury and T. J. Pence (Plenum, New York, 1997) p. 23.
11. M. V. R. Murthy, Phys. Rev. B 62 (2000) 17004.
12. A. Rettori and J. Villain, J. Phys. France 49 (1988) 257.
13. F. Lancon and J. Villain, Phys. Rev. Lett. 64 (1990) 293.
14. M. Ozdemir and A. Zangwill, Phys. Rev. B 42 (1990) 5013.
15. H. Spohn, J. Phys. I France 3 (1993) 69.
16. J. Hager and H. Spohn, Surf. Sci. 324 (1995) 365.
17. N. Israeli and D. Kandel, Phys. Rev. Lett. 80 (1998) 3300.
18. N. Israeli and D. Kandel, Phys. Rev. B 62 (2000) 13707.
19. Z. Jiang and C. Ebner, Phys. Rev. B 40 (1989) 316.
20. M. A. Dubson and G. Jeffers, Phys. Rev. B 49 (1994) 8347.
21. W. Selke and P. M. Duxbury, Phys. Rev. B 52 (1995) 17468.
22. P. C. Searson, R. Li and K. Sieradzki, Phys. Rev. Lett. 74 (1995) 1395.
23. Z. Jiang and C. Ebner, Phys. Rev. B 53 (1996) 11146.
24. M. V. R. Murthy and B. H. Cooper, Phys. Rev. B 54 (1996) 10377.
25. J. D. Erlebacher and M. J. Aziz, Surf. Sci. 374 (1997) 427.
26. A. Chame, F. Lancon, P. Politi, G. Renaud, I. Vilfan and J. Villain, Int. Jnl. Mod. Phys. B 11 (1997) 3657.
27. H. P. Bonzel, Interface Science, 9 (2001) 21.
28. A. Pimpinelli and J. Villain, Physics of crystal growth, (Cambridge University Press, Cambridge, 1998) Chap. 8.
29. A. Needleman and J. R. Rice, Acta Met. 28 (1980) 1315.
30. Z. Suo, Adv. Appl. Mech. 33 (1997) 193.
31. J. Erlebacher, M. J. Aziz, E. Chason, M. B. Sinclair and J. A. Floro, Phys. Rev. Lett. 84 (2000) 5800.
32. P. Nozieres, J. Phys. France 48 (1987) 1605.
33. W. H. Press, S. A. Teukolsky and W. T. Vetterling and B. P. Flannery, ”Numerical recipes: The art of scientific computing”, (Cambridge University Press, Cambridge, 1992) pages 708-715.
34. H. J. W. Zandvliet, Rev. Mod. Phys. 72 (2000) 593.
35. T. W. Poon, S. Yip, P. S. Ho and F. F. Abraham, Phys. Rev. B 45 (1992) 3521.
36. N. Israeli and D. Kandel, Phys. Rev. Lett. 88 (2002) art. no. 169601.


format: text

text:
# Euroattractor: a brief introduction to Iterated Function Systems

## I Introduction

Any function $`f`$ which maps a phase space $`\mathrm{\Omega }`$ into itself may be considered as a classical dynamical system. An arbitrary point $`x\mathrm{\Omega }`$ determines the trajectory uniquely, $`x_n=f^n(x_0)`$, where the natural index $`n`$ plays the role of time. Dynamical system defined by $`f`$ is deterministic. There exist different ways to introduce randomness into the system. One possibility relays on taking into account an additive noise, $`x_{n+1}=f(x_n)+\xi _n`$, where $`\xi _n`$ are uncorrelated random variables drawn according to a prescribed probability distribution $`P(\xi )`$. The deterministic case is obtained if $`P(\xi )`$ tends to the Dirac delta distribution, $`P(\xi )=\delta (\xi )`$.

In this work we analyze another class of stochastic systems, in which the choice of a system is random at each step. Such systems are called Iterated Function System (IFS), and may be considered as a generalization of the notion of a classical dynamical system.

An IFS is defined Ba  by a set of $`k`$ functions $`f_i:\mathrm{\Omega }\mathrm{\Omega }`$, which represent dynamical systems in the classical phase space $`\mathrm{\Omega }`$. The functions $`f_i`$ act randomly with given probabilities $`p_i`$. They characterize the likelihood of choosing a particular map at each step of the time evolution of the system. In general, the probabilities may be place dependent.

Analyzing a given iterated function system, one may ask, how an initial point $`x_0\mathrm{\Omega }`$ is transformed by the random process. One may also consider a more general question, how does a probability measure $`\mu `$ on $`\mathrm{\Omega }`$ change with respect to the Frobenius-Perron operator $`M`$ associated with the IFS. If the functions $`f_i`$ are strongly contracting, then there exists unique invariant measure $`\mu _{}`$ of $`M`$ – see for instance K81 ; E87 ; Ba ; BDEG88  and references therein.

Interestingly, for a large class of IFS’s the invariant measure $`\mu _{}`$ is characterized by a fractal structure. If the invariant measure is attracting, such IFS’s may be used to generate fractal sets in the space $`\mathrm{\Omega }`$. In particular, there exist iterated function systems leading to popular fractal sets: Cantor set, Sierpinski carpet or Sierpinski gasket Ba . In fact by considering a slightly generalized class of IFS’s it is possible to design a random system, the invariant measure of which possesses any desired properties. In this work we shall visualize the statement, by presenting an Euroattractor - a system, the attractor of which has the shape resembling the contours of Europe.

This paper is organized as follows. In the next section we recall the definition and basic properties of the standard IFS. Section III is devoted to a weakly contracting IFS, which exhibit nontrivial invariant measures, although they do not fulfill the standard contractivity assumptions. In section IV we define a generalized class of the systems and present the definition of an Euroattractor.

## II Iterated Function Systems

An iterated function system (IFS) is defined by a set of $`k`$ functions (dynamical systems) $`f_i`$, $`i=1,\mathrm{},k`$, defined on a compact space $`\mathrm{\Omega }`$ and the set of $`k`$ positive probabilities $`p_i`$, which sum to unity. In the IFS of the first kind the probabilities are constant, while in the more general, IFS of the second kind BDEG88 ; SKZ00 , they are functions of the position, $`p_i=p_i(x)`$. For any $`x\mathrm{\Omega }`$ the condition $`_{i=1}^kp_i(x)=1`$ has to be fulfilled. In general, IFS $``$ is defined by a set $`=\{\mathrm{\Omega },f_i:\mathrm{\Omega }\mathrm{\Omega },p_i:\mathrm{\Omega }[0,1]\}`$ with $`i=1,\mathrm{},k`$.

Any IFS generates a Frobenius–Perron operator $`M`$, which describes the time evolution of the probability densities on $`\mathrm{\Omega }`$. In the simplest case, in which the space $`\mathrm{\Omega }`$ is an interval in $``$ and the maps $`f_i`$ are invertible, any initial probability density $`\rho (x)`$ is mapped into $`\rho ^{}(x)`$ by the Frobenius–Perron operator LM94

$$\rho ^{}(x)=M[\rho (x)]=\underset{i=1}{\overset{k}{}}p_i\left(f_i^1(x)\right)\rho \left(f_i^1(x)\right)\left|\frac{df_i^1(x)}{dx}\right|.$$
(1)
If $`\rho (x)`$ is a non-singular probability density, so is $`\rho ^{}(x)`$. However, in the limit $`n\mathrm{}`$ the image $`M^n[\rho (x)]`$ needs not to be a smooth probability density, so it is advantageous to work in the larger space of probability measures on $`\mathrm{\Omega }`$. This space includes also also singular measures like the Dirac delta distributions and their convex combinations.

If an IFS is sufficiently contracting, i.e.

a) the functions $`f_i`$ satisfy the Lipschitz condition, $`d(f_i(x),f_i(y))L_id(x,y)`$) for any $`x,y\mathrm{\Omega }`$. Here $`d(x,y)`$ denotes the distance between both points and the Lipschitz constants $`L_i`$ are assumed to be smaller than unity,

b) the probabilities $`p_i(x)`$ are continuous functions, strictly positive for any $`x\mathrm{\Omega }`$, the IFS is called hyperbolic and there exists a unique invariant probability measure $`\mu _{}`$ satisfying the equation $`M\mu _{}=\mu _{}`$ LM94 .

Interestingly, the conditions a) and b) are not necessary to assure the existence of a unique invariant probability measure. Some other less restrictive sufficient assumptions which regard also IFS of the second kind, with place dependent probabilities, were analyzed in K81 ; E87 ; BE88 ; BDEG88 ; LY94 ; JL95 ; E96 ; Sz00 .

The invariant measures of an IFS display often fractal properties. For instance, consider an IFS of the first kind defined on an interval and consisting of $`k=2`$ affine transformations,

$$_1=\{\mathrm{\Omega }=[0,1],f_1(x)=x/3,f_2(x)=x/3+2/3;p_1=p_2=\frac{1}{2}\}.$$
(2)
Since both functions are continuous contractions with Lipschitz constants $`L_1=L_2=1/3<1`$, and the probabilities $`p_i`$ are constant and nonzero this IFS is hyperbolic. Thus there exists a unique, attracting invariant measure $`\mu _{}`$. It is not difficult to show Ba  that $`\mu _{}`$ is concentrated uniformly on the Cantor set of the fractal dimension $`D=\mathrm{ln}2/\mathrm{ln}3`$.

Although $`\mu _{}`$ is singular, it is possible to give a general prescription, how to compute integrals with respect to this invariant measure. This is due to the fact that the measure $`\mu _{}`$ is attractive, i.e., $`M^n\mu `$ converges weakly to $`\mu _{}`$ if $`n\mathrm{}`$. Therefore, in order to obtain the exact value of an integral $`_\mathrm{\Omega }h(x)𝑑\mu _{}`$ for any continuous function $`h(x)`$ it is sufficient to find the limit of the sequence $`h_n:=_\mathrm{\Omega }h(x)d(M^n\mu )`$ for an arbitrary initial measure $`\mu `$. This method of computing integrals over the invariant measure $`\mu _{}`$ is purely deterministic Ba ; E96 ; Sl97 ; SKZ00 . Alternatively, one may employ a random iterated algorithm by generating by the IFS a random sequence $`x_j\mathrm{\Omega }`$, $`j=0,\mathrm{},n1`$ which originates from an arbitrary initial point $`x_0`$. Due to the the ergodic theorem for IFS’s K81 ; E87 ; IG90  the mean value $`\frac{1}{n}[_{j=0}^{n1}h(x_j)]`$ converges in the limit $`n\mathrm{}`$ to the desired integral, $`_\mathrm{\Omega }h(x)𝑑\mu _{}`$, with probability one.

To show an example of an IFS of the second kind let us change probabilities in $`_1`$ and define

$$_2=\{\mathrm{\Omega }=[0,1],f_1(x)=x/3,f_2(x)=x/3+2/3;p_1(x)=x,p_2(x)=1x\}.$$
(3)
The probability $`p_2`$ vanish at the point $`x=1`$, so this IFS is not hyperbolic. In spite of this fact there exists the unique invariant measure $`\mu _{}`$ concentrated in a non-uniform way on the Cantor set SKZ00 . In this case $`\mu _{}`$ displays multifractal properties, since the generalized dimension of the measure depends on the Rényi parameter BS93 .

## III Weakly contracting IFS

Contraction conditions a) and b) imply existence of a unique invariant measure for the IFS. However, these conditions are too strong, and there exist IFSs with an attracting invariant measure which do not obey them. If the space $`\mathrm{\Omega }^m`$ and the functions $`f_i`$ are linear, the Lipschitz condition a) is equivalent to the statement that the moduli of all elements of the Jacobi matrices $`J_i:=f_i(x_1,\mathrm{},x_m)/(x_1,\mathrm{},x_m)`$ are strictly smaller then one.

This is not the case for the following IFS defined on the square

$$_3=\left\{\begin{array}{ccc}\hfill \mathrm{\Omega }=[0,1]\times [0,1],k& =& 2;p_1=p_2=1/2\hfill \\ \hfill f_1(x,y)& =& (xy,y)\hfill \\ \hfill f_2(x,y)& =& ((x1)y+1,y+\delta )\hfill \end{array}\right\},$$
(4)
since the Jacobi matrices

$$J_1=\left(\begin{array}{cc}y& x\\ 0& 1\end{array}\right)\mathrm{and}J_2=\left(\begin{array}{cc}y& x1\\ 0& 1\end{array}\right)$$
(5)
contain entries equal to unity. This reflects the fact that both functions $`f_1`$ and $`f_2`$ are not contracting in the $`y`$–direction. Numerical experiments suggest that in spite of this fact for irrational values of $`\delta 1`$, the associated Frobenius–Perron operator $`M`$ possess a unique invariant measure $`\mu _{}`$. Such an IFS, not fulfilling the condition a) or b) will be called weakly contracting.

Each random trajectory of the IFS $`_3`$ exhibits a slow diffusion in the $`y`$–direction, the speed of which is governed by the rotation parameter $`\delta `$. Due to periodic boundary conditions, $`(x,y)=(x|_{\mathrm{mod1}},y|_{\mathrm{mod1}})`$, every trajectory is confined to the square $`\mathrm{\Omega }`$. If the parameter $`\delta `$ is chosen to be irrational, there are no periodic trajectories. If $`\delta `$ is sufficiently small than the contraction in the $`x`$–direction will dominate the diffusion in $`y`$, so the trajectory will be sticked to the support of the invariant measure $`\mu _{}`$. Figure 1 shows a typical trajectory consisting of $`6000`$ points generated by this IFS. To present the asymptotic properties of the invariant measure the first $`100`$ points were omitted.

Note that if the parameter $`\delta `$ is set to zero, the variable $`y`$ would be fixed to the initial value $`y_0`$. Such an IFS would have an invariant measure concentrated on the interval $`y=y_0`$. Analyzing the contraction constant we see that for $`y_0=1/3`$ its support is equal to the Cantor set generated by the Cantor IFS (2). In general, for any $`y_0<1/2`$ its dimension would be

$$D(y_0)=\frac{\mathrm{ln}2}{\mathrm{ln}y_0}.$$
(6)
We conjecture that this formula gives also the fractal dimension of the horizontal cross-section at $`y=y_0`$ of the invariant set of the IFS (4) in the limit $`\delta 0`$.

Another example of an invariant measure of a nonlinear, weakly contracting IFS is presented in Fig.2. The system is defined by

$$_4=\left\{\begin{array}{ccc}\hfill \mathrm{\Omega }=[0,1]\times [0,1],k& =& 2;p_1=p_2=1/2\hfill \\ \hfill f_1(x,y)& =& (\frac{1}{2}x\mathrm{sin}^2(\pi y),y)\hfill \\ \hfill f_2(x,y)& =& (\frac{1}{2}(x1)\mathrm{sin}^2(\pi y)+1,y+\delta )\hfill \end{array}\right\}.$$
(7)
As in the previous example each random realization of the IFS generates a trajectory confined in the torus, and for small values of the shift parameter $`\delta `$ the contraction in the horizontal direction dominates the vertical diffusion. For $`\delta =0`$ the IFS reduces to one-dimensional, and the fractal dimension of its invariant set reads

$$D(y_0)=\frac{\mathrm{ln}2}{\mathrm{ln}22\mathrm{ln}\mathrm{sin}(\pi /y_0)}.$$
(8)
This number is smaller than one for all initial values apart of $`y_0=1/2`$. Numerical results performed for $`\delta 1`$ suggest that the same formula describes the dimension of the cross-section at $`y=y_0`$ of the invariant set of the IFS (7) in the limit $`\delta 0`$. In other words, the support of the invariant measure $`\mu _{}`$ of the IFS $`_4`$ may be considered as a collection of Cantor sets, with a continuously varying dimension $`D(y_0)`$, placed horizontally, one parallel to another.

## IV Repellers and generalized IFS

In the original definition of the IFS, presented in section II, all functions $`f_i`$ map the space $`\mathrm{\Omega }`$ into itself. One may, however, relax this condition, and consider generalized IFS constructed of functions which map some fragments of $`\mathrm{\Omega }`$ out of this space. Such functions are called repellers, since they do not conserve the probability, and during the time evolution the probability flux leaks out of the system.

A simple $`1`$-D repeller is obtained from the standard logistic map, $`f(x)=rx(1x)`$, if the control parameter $`r`$ exceeds the critical value, $`r_c=4`$ Ra89 . In such systems a randomly selected initial value leads to a trajectory, which escapes the system with probability one. The complementary set of initial points which generate infinite trajectories is of measure zero and often forms a fractal set. The same property is characteristic to a large class of chaotic 1-D maps with a gap - a subset of $`\mathrm{\Omega }`$ mapped out of $`\mathrm{\Omega }`$. An analysis of the fractal dimension of the invariant set for the tent map with a gap is presented in ZB99 .

In this chapter we will consider an IFS defined on the square, $`\mathrm{\Omega }=[0,1]\times [0,1]`$, consisting of $`k`$ repelling functions. Assume that each function $`f_i`$ restricted to a given set $`X_i\mathrm{\Omega }`$ acts as an affine transformation defined by real matrices $`J_i`$ of size $`2`$ and two components translation vectors $`\stackrel{}{\kappa }_i`$, but the complementary subset $`\mathrm{\Omega }X_i`$ is mapped out of $`\mathrm{\Omega }`$,

$$f_i\left(\begin{array}{c}x\\ y\end{array}\right)=\{\begin{array}{c}J_i\left(\begin{array}{c}x\\ y\end{array}\right)+\stackrel{}{\kappa }_i\mathrm{for}xX_i\\ \mathrm{out}\mathrm{of}\mathrm{\Omega }\mathrm{for}xX_i\end{array}$$
(9)
Such functions $`f_i`$ are repelling, so instead of analyzing single trajectories, generated by a random initial point, it is convenient to work with initially smooth probability measures. Let us take for the initial measure $`\mu _0`$ the uniform (Lebesgue) measure on $`\mathrm{\Omega }`$. To iterate the repelling IFS we find the images of the measures covering each of the subsets $`X_i`$ by the maps $`f_i`$, $`i=1,\mathrm{},k`$ and add all contributions with prescribed weights – probabilities $`p_i`$. More formally, we use the the dimensional analogue of the Frobenius-Perron operator (1)

$$\rho ^{}(x,y)=M[\rho (x,y)]=\underset{i=1}{\overset{k}{}}p_i\left(f_i^1(x,y)\right)\rho \left(f_i^1(x,y)\right)\left|\frac{f_i(x,y)}{(x,y)}\right|^1,$$
(10)
where the contributions from the function $`f_i`$ have to be taken into account only if the analyzed point $`(x,y)`$ belongs to $`f_i(X_i)`$. Then the preimage $`f_i^1(x,y)`$ exists, so the $`i`$-th term in the expression (10) is meaningful.

As an example of such a generalized, repelling system we consider the IFS $`_5`$ consisting of $`k=13`$ functions with equal weights, $`p_i=1/13`$. The detailed form of the affine functions $`f_i`$ as well as the definitions of the sets $`X_i`$ are provided in the Appendix. Figure 4 shows time evolution of the initially uniform measure $`\mu _0`$. The right hand side of the figure shows the measures $`\mu _n=M^n[\mu _0]`$, while the left hand side shows the sets, where the measure $`\mu _n`$ is concentrated. Initially uniform measure $`\mu _0`$ converges fast to the asymptotic measure $`\mu _{}`$, and already after $`n=6`$ iterations can hardly be distinguished from it. The shape of the support of $`\mu _{}`$ explains why the name Euroattractor has been concocted for this system.

## V Closing Remarks

The purpose of this work was to demonstrate usefulness of the concept of Iterated Function System. This generalization of the concept of a dynamical system incorporates deterministic and stochastic behavior and may be applied in various fields of physics, and in general, of nonlinear science. A possible extension of the concept of IFS for quantum mechanics was recently worked out in LZS02 .

The IFS’s belong to a larger class of random systems studied in YOC91 ; PSV95 . From a mathematical point of view several questions concerning the rigorous theory of IFS remain open. For instance, the sufficient conditions, under which an IFS possesses a unique, attracting invariant measure are known, but the formulation of the necessary conditions remains still an exciting mathematical challenge LY94 ; G92 ; GB95 ; E96 . In this work we have introduced a weakly contracting IFS’s $`_3`$ and $`_4`$, which do not satisfy the standard contractivity conditions, but do posses attracting invariant measures.

For any given set $`S`$ in the plane one may try to find an IFS, such that its invariant measure $`\mu _{}`$ has support on $`S`$. This versatility of IFS’s allows one to apply them for image synthesis, image processing and image encoding Ba ; BE88 . For instance, the attractor of a generalized IFS $`_5`$, described in this work, has the shape resembling the contours of Europe111Analyzing political changes which took place in Europe during the last decade, one may ask, whether European Union has enough attractive power to become a global attractor for all countries in this continent. The theory of nonlinear systems prevents us from predicting the time evolution of an unstable dynamical system in a long run. Should we than just observe the events waiting passively? Or perhaps, taking the lesson from the Lorenz’s butterfly, could it be sufficient to flap the wings gently in the very right moment?.

It is a pleasure to thank Radosława Bach and Wojciech Słomczyński for constructive remarks and fruitful interaction. This paper was presented at the European Interdisciplinary School on Nonlinear Dynamics Euroattractor 2002 organized in Warsaw in June 2002. We are grateful to Włodzimierz Klonowski for an inspiring name of the conference as well as for the opportunity to present this work during that event.

## Appendix A An Euroattractor

In this appendix we provide an explicit definition of an Euroattractor – a generalized, repelling IFS $`_5`$, the invariant measure of which is shown in Fig.4. The space $`\mathrm{\Omega }`$ is equal to an rectangle. To address each pixel in a simple way its size is taken to be $`500\times 600`$.

The IFS consists of $`k=13`$ affine functions $`f_i`$ defined in (9), and all the weights are constant, $`p_i(x,y)=1/13`$. The parameters defining the corners of the quadrangles $`X_i`$, and the elements of the transformation matrices $`J_i`$ and the translation vectors $`\stackrel{}{\kappa }_i`$ are collected in Table 1. Observe that some entries of the Jacobi matrix $`J_6`$ have modulus larger than one, so this IFS does not satisfy the contractivity condition.

Table 1. Parameters defining an Euroattractor – the IFS $`_5`$.

| $`i`$ | $`X_i`$ | $`J_i`$ | $`\stackrel{}{\kappa }_i`$ |
| --- | --- | --- | --- |
| 1 | $`\left(\begin{array}{c}90\\ 370\end{array}\right)\left(\begin{array}{c}170\\ 520\end{array}\right)\left(\begin{array}{c}220\\ 350\end{array}\right)\left(\begin{array}{c}300\\ 500\end{array}\right)`$ | $`\left(\begin{array}{cc}0.37& 0.66\\ 0.74& 0.66\end{array}\right)`$ | $`\left(\begin{array}{c}276\\ 371\end{array}\right)`$ |
| 2 | $`\left(\begin{array}{c}520\\ 430\end{array}\right)\left(\begin{array}{c}340\\ 320\end{array}\right)\left(\begin{array}{c}420\\ 530\end{array}\right)\left(\begin{array}{c}240\\ 420\end{array}\right)`$ | $`\left(\begin{array}{cc}0.33& 0.73\\ 0.80& 0.59\end{array}\right)`$ | $`\left(\begin{array}{c}626\\ 355\end{array}\right)`$ |
| 3 | $`\left(\begin{array}{c}170\\ 350\end{array}\right)\left(\begin{array}{c}265\\ 480\end{array}\right)\left(\begin{array}{c}230\\ 320\end{array}\right)\left(\begin{array}{c}325\\ 450\end{array}\right)`$ | $`\left(\begin{array}{cc}0.52& 0.38\\ 0.01& 0.68\end{array}\right)`$ | $`\left(\begin{array}{c}16\\ 218\end{array}\right)`$ |
| 4 | $`\left(\begin{array}{c}340\\ 300\end{array}\right)\left(\begin{array}{c}240\\ 420\end{array}\right)\left(\begin{array}{c}500\\ 400\end{array}\right)\left(\begin{array}{c}400\\ 520\end{array}\right)`$ | $`\left(\begin{array}{cc}0.34& 0.54\\ 0.61& 0.51\end{array}\right)`$ | $`\left(\begin{array}{c}192\\ 814\end{array}\right)`$ |
| 5 | $`\left(\begin{array}{c}380\\ 500\end{array}\right)\left(\begin{array}{c}500\\ 300\end{array}\right)\left(\begin{array}{c}250\\ 450\end{array}\right)\left(\begin{array}{c}370\\ 250\end{array}\right)`$ | $`\left(\begin{array}{cc}0.17& 0.95\\ 0.96& 0.08\end{array}\right)`$ | $`\left(\begin{array}{c}812\\ 31\end{array}\right)`$ |
| 6 | $`\left(\begin{array}{c}240\\ 460\end{array}\right)\left(\begin{array}{c}340\\ 550\end{array}\right)\left(\begin{array}{c}260\\ 440\end{array}\right)\left(\begin{array}{c}360\\ 530\end{array}\right)`$ | $`\left(\begin{array}{cc}1.44& 1.05\\ 0.23& 1.26\end{array}\right)`$ | $`\left(\begin{array}{c}497\\ 24\end{array}\right)`$ |
| 7 | $`\left(\begin{array}{c}200\\ 490\end{array}\right)\left(\begin{array}{c}50\\ 590\end{array}\right)\left(\begin{array}{c}150\\ 270\end{array}\right)\left(\begin{array}{c}0\\ 370\end{array}\right)`$ | $`\left(\begin{array}{cc}0.31& 0.47\\ 0.98& 0.22\end{array}\right)`$ | $`\left(\begin{array}{c}625\\ 97\end{array}\right)`$ |
| 8 | $`\left(\begin{array}{c}510\\ 430\end{array}\right)\left(\begin{array}{c}460\\ 540\end{array}\right)\left(\begin{array}{c}305\\ 330\end{array}\right)\left(\begin{array}{c}245\\ 440\end{array}\right)`$ | $`\left(\begin{array}{cc}0.47& 0.032\\ 0.39& 0.99\end{array}\right)`$ | $`\left(\begin{array}{c}485\\ 438\end{array}\right)`$ |
| 9 | $`\left(\begin{array}{c}230\\ 280\end{array}\right)\left(\begin{array}{c}230\\ 450\end{array}\right)\left(\begin{array}{c}460\\ 280\end{array}\right)\left(\begin{array}{c}460\\ 450\end{array}\right)`$ | $`\left(\begin{array}{cc}0.26& 1\\ 0.60& 0.61\end{array}\right)`$ | $`\left(\begin{array}{c}20\\ 267\end{array}\right)`$ |
| 10 | $`\left(\begin{array}{c}440\\ 280\end{array}\right)\left(\begin{array}{c}440\\ 450\end{array}\right)\left(\begin{array}{c}245\\ 280\end{array}\right)\left(\begin{array}{c}245\\ 450\end{array}\right)`$ | $`\left(\begin{array}{cc}0.51& 0.70\\ 0.87& 0.41\end{array}\right)`$ | $`\left(\begin{array}{c}408\\ 339\end{array}\right)`$ |
| 11 | $`\left(\begin{array}{c}250\\ 280\end{array}\right)\left(\begin{array}{c}310\\ 190\end{array}\right)\left(\begin{array}{c}220\\ 190\end{array}\right)\left(\begin{array}{c}270\\ 100\end{array}\right)`$ | $`\left(\begin{array}{cc}0.22& 0.37\\ 0.33& 0.11\end{array}\right)`$ | $`\left(\begin{array}{c}439\\ 456\end{array}\right)`$ |
| 12 | $`\left(\begin{array}{c}330\\ 0\end{array}\right)\left(\begin{array}{c}200\\ 190\end{array}\right)\left(\begin{array}{c}410\\ 100\end{array}\right)\left(\begin{array}{c}280\\ 290\end{array}\right)`$ | $`\left(\begin{array}{cc}0.60& 0.51\\ 0.47& 0.38\end{array}\right)`$ | $`\left(\begin{array}{c}359\\ 187\end{array}\right)`$ |
| 13 | $`\left(\begin{array}{c}420\\ 10\end{array}\right)\left(\begin{array}{c}420\\ 90\end{array}\right)\left(\begin{array}{c}330\\ 10\end{array}\right)\left(\begin{array}{c}330\\ 90\end{array}\right)`$ | $`\left(\begin{array}{cc}0& 0.31\\ 0.66& 0\end{array}\right)`$ | $`\left(\begin{array}{c}243\\ 30\end{array}\right)`$ |


text:
# Secure communication using coherent states

## Acknowledgments

This work was supported by DARPA/AFRL under grant # F30602-01-2-0528.


text:
# Eta physics at threshold

## 1 Manifestation of the $`\eta `$-nucleon-nucleon interaction

In the last decade large experimental as well as theoretical efforts were concentrated on the study of the creation of $`\eta `$ and $`\eta ^{}`$ mesons via the hadronic interactions. Measurements of the production of these mesons in the elementary nucleon-nucleon collision have been performed in the vicinity of the kinematical threshold where only one partial wave in both initial and final state is expected to contribute to the production process. For example in case of the proton-proton collision the dominance of the $`{}_{}{}^{3}P_{0}^{}Ss`$ transition is expected up to an excess energy of about 40 MeV and 100 MeV for $`\eta `$ and $`\eta ^{}`$ meson, respectively review . This simplifies significantly the interpretation of the data, yet still appears to be challenging due to the three particle final state system with a complex hadronic potential. The determined energy dependences of the total cross section for $`\eta ^{}`$ etap\_data ; hibou  and $`\eta `$ hibou ; eta\_data  mesons in proton-proton collisions are presented in figures 1 and 2. Comparing the data to the arbitrarily normalized phase-space integral (dashed lines) reveals that the proton-proton FSI enhanced the total cross section by more than an order of magnitude for low excess energies. One recognizes also that in the case of the $`\eta ^{}`$ meson the calculation –assuming that the on-shell proton-proton amplitude exclusively determines the phase-space population– describes the data very well (solid line). This indicates that the proton-$`\eta ^{}`$ interaction is too small to be observed within the present accuracy swave . In the case of the $`\eta `$ meson the increase of the total cross section for very low and very high energies is much larger than expected from the final state interaction between protons. The excess at higher energies can be assigned to the onset of higher partial waves, and the enhancement at threshold can be plausibly explained by the influence of the attractive interaction between the $`\eta `$ meson and the proton.
Though the simple phenomenological treatment –based on factorization of the transition amplitude into the constant primary production and the on-shell incoherent pairwaise interaction among particles– works well for the energy dependence of the total cross section, it fails completely as far as the description of the differential cross section is concerned (see the thin solid line in figure 3).

This discrepancy is rather too large to be utterly caused by the underestimation of the S-wave proton-$`\eta `$ interaction. An explanation could be a contribution from P-wave proton-proton interaction hanhart  or a significant influence of the off-shell effects of the interaction between outgoing particles. Indeed, a much better description of the differantial cross section shown in figure 3 is achieved when the enhancement due to the proton-proton FSI is modeled by the inverse of the squared Jost function instead of the on-shell amplitude of the elastic proton-proton scattering (see thick solid line in fig. 3). The Jost function accounts for the off-shell effects of the interacting particles, however it gives wrong energy dependence of the total cross section, and unfortunately it depends rather strongly on the potential model baru ; review  and applied formalism kle01 . Therefore in the framework of this phenomenology it is not possible –with a satisfactory accuracy– to separate the contribution of the proton-$`\eta `$ interaction from other effects in a model independent way. Thus, more sophisticated theoretical calculations are required. An estimation of the nucleon-$`\eta `$ scattering parameters has successfuly been performed by comparing the close-to-threshold data of the $`\eta `$ meson photoproduction on a deuteron hejny  to the calculations performed to the first order rescattering in the two-body subsystems elster . The recent investigation of Fix and Arenhövel fix  shows, however, that the three-body dynamics contribute substantially to the total and differential cross sections of the $`\eta `$ photoproduction on the deuteron, and thus it cannot be neglected in the quantitative analysis.

Data of the $`\eta `$-NN system created via the hadronic interaction of nucleons still awaits a strict theoretical explanation. The hitherto performed analysis of the energy dependence of the total cross section reveals the sensitivity of the data to the $`\eta `$-nucleon interaction. Clearly, further comparative studies of the $`\eta `$-nucleon-nucleon system created in photo- and hadro-production should lead to a better determination of the $`\eta `$-nucleon potential and specifically should be helpful in distinguishing between effects due to the production dynamics and the final state interaction among the produced particles.

## 2 Production dynamics

It is rather well established that close to threshold the energy dependence of the total cross section is due to the interaction among the outgoing particles and that the entire production dynamics manifests itself only in a single constant which determines the magnitude of the total cross section review . Therefore in spite of the precise measurements of the total cross section for the creation of $`\eta `$ and $`\eta ^{}`$ mesons in proton-proton– as well as of $`\eta `$ meson in proton-neutron collisions there is still a lot of ambiguities in the description of the mechanism underlying the production process. It is generally anticipated eta\_theo ; eta\_rho ; eta\_all  that the $`\eta `$ meson is produced predominantly via the excitation of the S11 baryonic resonance $`N^{}(1535)`$, whose creation is induced through the exchange of the virtual $`\pi `$, $`\eta `$, $`\rho `$, $`\sigma `$ and $`\omega `$ mesons, however at present it is still not established what are the relative contributions originating from a particular meson. Measurements of the total cross section in different isospin channels put some more limitations to the models, yet still the $`\eta `$ meson production in the $`pppp\eta `$ and $`pnpn\eta `$ can be equally well described by eg. assuming the $`\rho `$ meson exchange dominance eta\_rho  or by taking contributions from the pseudoscalar and vector meson exchanges eta\_all . Therefore for a full understanding of the production dynamics the determination of polarisation observables is mandatory. Already the precise measurement of the beam analyzing power could exclude one of the above mentioned possibilities. The first measurement of that quantity have been recently performed poleta , but for a conclusive inference a better accuracy of the data is required.

Until now it was also not possible to determine unambiguously the mechanism of the production of the $`\eta ^{}`$ meson review . Model uncertainties and in practice only one number –namely the magnitude of the total cross section of the $`pppp\eta ^{}`$ reaction– serving as input for theory are by far not sufficient to distinguish between the mesonic, nucleonic or resonant production currents etap\_theo . The understanding of that mechanism on the hadron and quark-gluon level will shed a light on the structure of the $`\eta ^{}`$ meson. A possibly large glue content of the $`\eta ^{}`$ wave function and the dominant flavour–singlet combination of its quarks component may cause that the dynamics of its production process in nucleon–nucleon collisions is significantly different from that responsible for the production of other mesons. A creation of that meson via a fusion of gluons excited in the interaction region should lead to the same production yield in proton-proton and proton-neutron collisions because gluons do not distinguish between flavours bass . Therefore it is interesting to determine experimentally the ratio $`R_\eta ^{}=\sigma (pnpn\eta ^{})/\sigma (pppp\eta ^{})`$ and to compare it to the already known value of $`R_\eta 6.5`$ calen2667 . In the extreme case when the $`\eta ^{}`$ meson is created only from glue excited in the interaction region the ratio $`R_\eta ^{}`$ should be close to unity after the correction for the final and initial state interactions between the interacting nucleons bass . On the contrary a production mechanism dominated by the isovector meson exchange would yield a significantly larger ratio as already observed in the case of the $`\eta `$ meson. The appropriate experiments are in preparation at the COSY-11 facility pnetapc11  and the first feasibility study has already been accomplished sucessfully.

## Acknowledgement

The work has been partly supported by the European Community - Access to Research Infrastructure action of the Improving Human Potential Programme