Theory of nanoscale surface ripple formation during oblique-incidence thin film deposition

In a series of fascinating experiments, Voronov and co-workers deposited multilayer Mo/Si thin films by oblique-incidence sputter deposition.¹ They discovered that the surface of the growing thin film spontaneously developed a rippled topography if the angle between the beam and the substrate normal $θ$ exceeded critical value $θc$⁠. The ripples were oriented so that, on average, their crests were perpendicular to the incident atomic flux. Similar results have recently been obtained by another group.² These findings are not just of academic interest: multilayer coatings consisting of alternating layers of high and low Z-materials have high x-ray reflectivities and are important in extreme ultraviolet lithography, x-ray imaging, high resolution x-ray spectroscopy, and x-ray astronomy. Moreover, in the fabrication of many advanced optical devices, multilayers must be deposited on strongly curved or nanostructured substrates, and this inevitably leads to oblique-incidence deposition. When multilayer blazed gratings for extreme ultraviolet radiation or soft x rays are produced, for example, multilayers are deposited on substrates with sawtooth topographies.^3,4

In the experiments of Voronov et al., the energies of the atoms arriving at the surface of the growing multilayer coating were on the order of 10 eV. This is below the sputter yield threshold for both Mo and Si, and so essentially no sputtering of atoms that had already been incorporated into the multilayer thin film occurred. Although the impinging atoms produced no sputtering, momentum transfer from the incident atoms to atoms near the solid surface would have led to inelastic displacement of atoms already resident in the solid, a phenomenon that is known as mass redistribution (MR).^5–7 Molecular dynamics simulations of the deposition of diamond-like carbon thin films, for example, show that arriving carbon atoms with energies as low as 30 eV produce craters.⁶ MR has a smoothing effect on a solid surface for $θ$ less than a critical value, but it has a destabilizing influence for higher angles of incidence.^5–7 Voronov et al. argued convincingly that MR was responsible for the ripple formation they observed in their experiments.¹ 

In this paper, we will develop a theory of surface ripples that emerge when an amorphous thin film is deposited by oblique-incidence sputter deposition. For simplicity, we will consider the deposition of a thin film composed of a single atomic species, rather than studying multilayer deposition as Voronov et al. did. We will also simplify the problem by analyzing the case in which two diametrically opposed, broad atomic beams are simultaneously incident on the substrate and focus our attention on $θ$ just above the threshold angle for ripple formation $θc$⁠. This makes the problem more tractable by ensuring that the surface slope remains small as time passes. As a result, no overhangs develop as the film grows and no part of the surface is shadowed from the incident atomic fluxes. In contrast, overhangs formed and shadowing occurred in the single-beam experiments conducted by Voronov et al. Shadowing is a nonlocal phenomenon that is particularly difficult to deal with theoretically.

Interestingly, we find that the ripples that form on the surface of the growing film develop a sawtooth form as time passes. The ripple crests and troughs are both rounded, but, in general, the extent to which the crests and the troughs are rounded differs. Thus, the surface ripples display vertical asymmetry, i.e., the appearance of the ripples changes if they are reflected about the substrate plane. At early times, the ripples are rounded and roughly sinusoidal in form, but as time passes, their wavelength and amplitude grow and they become more faceted. At long times, the ripples become highly ordered in the transverse direction and few dislocations are present. The ripple wavelength and amplitude vary in the longitudinal direction, however, and so the ripples are not perfectly ordered.

This paper is organized as follows. In Sec. II, we begin by considering the oblique-incidence sputter deposition of an amorphous thin film by a single atomic beam before moving on to the case in which two diametrically opposed beams are incident on the surface in Sec. III. An in-depth analysis of the dual-beam equation of motion (EOM) is carried out in Sec. IV. We discuss our results and place them in context in Sec. V and conclude in Sec. VI.

Consider the oblique-incidence sputter deposition of an amorphous thin film. We take the substrate surface to be nominally planar and place the origin on it. We orient the $z$ axis so that it is perpendicular to the substrate surface and so that it points away from the substrate. A single broad beam of atoms is now directed onto the substrate, and a thin film begins to accumulate on it, as shown in Fig. 1(a). For simplicity, we assume that the incident atoms all have the same energy and momentum. The angle that the beam direction makes with the $z$ axis will be called the angle of incidence and will be denoted by $θ$⁠. We orient the $x$ axis so that it is parallel to the projection of the beam direction onto the $x−y$ plane. The incident atomic flux is then $J=−Je^$⁠, where $e^≡−x^sinθ+z^cosθ$⁠.

We assume that the impinging atoms are all incorporated into the growing film, which is taken to have uniform density. The incident atoms will be taken to have an energy low enough that essentially no sputtering takes place. They are assumed to be sufficiently energetic that appreciable MR occurs, however. In contrast, in molecular beam epitaxy, the incident atoms or molecules have low enough energies that MR is completely negligible.

We will employ a continuum description of the surface dynamics in which the position of an arbitrary point on the film’s surface is given by $r=xx^+yy^+h(x,y,t)z^$⁠, where $h(x,y,t)$ is the height of the point above the $x−y$ plane at time $t$⁠. The surface height $h$ is obtained by coarse-graining the detailed microscopic surface configuration and is assumed to be a smoothly varying function of its arguments $x$⁠, $y$⁠, and $t$⁠. We assume that no overhangs are initially present or develop as the film grows, so that $h(x,y,t)$ is a single-valued function of $x$ and $y$ for all times $t$⁠. We will also assume that the slope of the surface remains small enough that shadowing of the incident atomic flux does not occur.

If there were no currents on the surface of the solid, the normal velocity of a point on the surface would simply be $vn=−ΩJ⋅n^$⁠, where $Ω$ is the atomic volume and

is the outward-pointing surface normal. Therefore, the EOM is in this case

where the subscripts on $h$ denote partial derivatives. Equation (2) is of course easy to solve for arbitrary initial surface heights $h(x,y,0)$⁠.

We now move on to consider the case in which there are surface currents. Let $j=jxx^+jyy^+jzz^$ be the atomic current on the surface of the solid. $j$ is everywhere tangent to the solid surface. The equation of motion for the surface is

The surface current $j$ depends on the angle of incidence $θ$⁠. Its value at time $t$ also depends in principle on the form of the entire surface at that time. In other words,

for $i=x$ and $y$⁠. We can simplify the EOM obtained by combining Eqs. (3) and (4) a little by setting $h(x,y,t)=ΩJtcos⁡θ+u(x,y,t)$⁠, i.e., by transforming to a frame of reference that moves upward with the spatially averaged surface velocity. We obtain

Equation (5) completely specifies the dynamics of the surface but, in general, it is an exceedingly complicated nonlinear partial differential equation. In addition, how $jx$ and $jy$ depend on all of the derivatives $ux$⁠, $uy$⁠, $uxx$⁠, $uxy$⁠, $uyy,…$ is not known for any atomic species, energy or angle of incidence. Fortunately, the surface dynamics become much simpler just above the threshold for pattern formation, i.e., when $ϵ≡(θ−θc)1/2$ is small and positive. We seek solutions to Eq. (5) of the form

where

$X$⁠, $Y$⁠, and $T$ are the so-called “slow” variables and $x$⁠, $y$⁠, and $t$ are the corresponding “fast” variables. Heuristically speaking, Eqs. (6) and (7) say that for $θ$ just above the critical angle $θc$⁠, the height of the surface disturbance varies slowly in space and time. In addition, the spatial variation in the $y$ direction is more gradual than in the $x$ direction because ripples with their wave vector pointing in the $x$ direction develop for $θ>θc$⁠. An a posteriori justification for adopting the scaling ansatz given by Eqs. (6) and (7) will be obtained once we have arrived at an EOM that is well-behaved in the $ϵ→0$ limit for the case in which diametrically opposed beams are incident on the substrate.

Using the scaling given by Eqs. (6) and (7) in Eq. (5), we obtain

Let $jx,n$ denote the derivative of $jx(θ;ux,uy,uxx,uxy,uyy,…)$ with respect to its $n$th argument after the semicolon evaluated for $ux=uy=uxx=uxy=uyy=⋯=0$⁠. We define $jy,n$ in an analogous fashion. The way in which quantities like $jx,n,m$ and $jy,n,m$ are defined should now be apparent. We expand Eq. (8) in powers of $ϵ$ and retain terms up to fourth order in $ϵ$⁠. The EOM must be invariant under transformation $Y→−Y$ and so the coefficients of terms that do not have this invariance must have coefficients equal to zero. We obtain

The constants $jx,1$⁠, $jx,1,1$⁠, $jx,1,1,1$⁠, $jy,2$⁠, $jx,3$⁠, $jx,6$⁠, and $jx,1,3$ that appear in Eq. (9) all depend on $θ$⁠.

The term $−ϵJ(sin⁡θ)UX$ that appears in Eq. (9) can be eliminated by transforming to a moving frame of reference. The resulting equation can then be recast in terms of dimensionless quantities, which yields an EOM that depends on three dimensionless parameters. This EOM does not have a well defined $ϵ→0$ limit, however. Physically, the reason for this is that the scaling ansatz given by Eqs. (6) and (7) implies that $ux=ϵUX$ is of order $ϵ$ and so is small, but if a single beam is incident on the film, surface slopes of order one develop.¹ This means that the scaling ansatz is invalid for the single beam case.

We now turn our attention to the problem in which there are two diametrically opposed atomic beams, each with flux $J$⁠. These beams both have the polar angle $θ$ but one has azimuthal angle $0°$ and the other has an azimuthal angle of $180°$⁠, as shown in Fig. 1(b). If only the beam that is incident from the left were present, the EOM would be Eq. (9). Conversely, if only the beam that is incident from the right were present, the EOM would be Eq. (9) with $X$ replaced by $−X$⁠, i.e.,

When both beams are turned on, $−Ω−1ϵ4UT$ is obtained by summing the right-hand sides of Eqs. (9) and (10). Hence

Note that in addition to being invariant under transformation $Y→−Y$⁠, the dual-beam EOM (11) is invariant under $X→−X$⁠, as it must be.

Equation (11) holds for $θ$ just above the critical angle $θc$⁠, i.e., for small $ϵ=(θ−θc)1/2$⁠. Experiments performed with a single incident atomic beam show that ripples with their wave vector in the $x$ direction form for $θ>θc$⁠, while the surface remains flat for $θ<θc$⁠.^1,2 We assume that this would also be the case if the experiments had been carried out with diametrically opposed beams. This assumption implies that $jx,1$ is negative for $θ<θc$⁠, zero for $θ=θc$⁠, and positive for $θ>θc$⁠, at least provided that $θ$ is not too large. Accordingly, for $θ$ close to the critical angle, $jx,1≅A1(θ−θc)=A1ϵ2$⁠, where $A1$ is a positive constant that does not depend on $θ$⁠. On the other hand, $jy,2$ is negative for $θ=θc$ because it has been assumed that ripples with their wave vector pointing in the $y$ direction do not develop for $θ$ just below $θc$⁠. Equation (11) may now be written

Notice that $ϵ$ does not appear in Eq. (12). Thus, with the scaling we posited in Eqs. (6) and (7), the dual-beam EOM has a well defined $ϵ→0$ limit.

We must have $jx,6>0$⁠, since otherwise arbitrarily short wavelengths are linearly unstable and the continuum description breaks down. In addition, $jx,1,1,1$ cannot be positive since if it were, the slope of the surface would grow without bound. Equation (12) becomes

where $A≡2ΩA1$⁠, $B≡2Ωjx,6$⁠, $D≡2Ω|jy,2|$⁠, $α≡Ω|jx,1,1,1|/3$⁠, and $β≡−Ωjx,1,3$⁠. The coefficients $A$⁠, $B$⁠, and $D$ in Eq. (13) are positive while $α$ is nonnegative. In addition, by replacing $U$ by $−U$ if necessary, we can arrange for $β$ to be nonnegative.

The term $−AUXX$ in Eq. (13) describes the destabilizing effect of MR above the critical angle. Thermally activated surface diffusion could produce the term $−BUXXXX$⁠,⁸ but we suspect that MR could also contribute to this term. MR is smoothing for disturbances in the transverse direction for all angles of incidence.⁷ This is the physical origin of the term $DUYY$⁠. The term $α∂XUX3$ in Eq. (13) is familiar from the theory of mounding instability that can occur during molecular beam epitaxy. There the term results from the Ehrlich–Schwoebel effect and can lead to the formation of a faceted surface.⁹ In our problem, the slope dependence of the surface current produced by MR leads to the presence of this term. The term $β∂X2uX2$⁠, on the other hand, is the so-called conserved Kuramoto–Sivashinsky (CKS) nonlinearity. Although this term was first encountered in molecular beam epitaxy,^9,10 it is also believed to play a role in the deposition of amorphous thin films.¹¹ The CKS nonlinearity tends to produce coarsening of the surface patterns, i.e., the characteristic lateral and vertical length scales increase with time. It also breaks the $U→−U$ symmetry that would be present if $β$ were zero. Since there is vacuum above the surface and solid below, there is no reason that such a symmetry should exist.

Our EOM (13) conserves mass, i.e., $∫−∞∞dX∫−∞∞dYU(X,Y,T)$ is a constant of the motion. This is a result of transforming to a frame of reference that moves upward with the spatially averaged surface velocity of the growing thin film.

How slowly must the height of the surface disturbance $u$ vary with $x$ for Eq. (13) to be valid? Let $a0$ be the characteristic size of the collision cascade, or, equivalently, the characteristic lateral length scale of the crater that the impact of a single atom produces on the surface of the growing thin film. The characteristic lateral length scale (or length scales) of the surface ripple must be large compared to $a0$ for Eq. (13) to be a good approximation. If the surface disturbance is sinusoidal, it has a single characteristic lateral length scale—its wavelength. On the other hand, if the surface disturbance is periodic but the width of the crests $w+$ differs from the width of the troughs $w−$⁠, then the characteristic lateral length scales $w+$ and $w−$ must both be large compared to $a0$ for Eq. (13) to apply.

We introduce the dimensionless quantities

Making these substitutions in Eq. (13) and dropping the tildes results in an equation of motion with a single dimensionless parameter $ψ$⁠,

where $ψ∈[0,π/2]$ is defined by the relation $tan⁡ψ=β/αB$⁠. This parameter is a measure of the relative strength of the quadratic and cubic nonlinearities.

Because the substrate surface is nominally flat, Eq. (15) can be linearized at early times. The amplitude of a sinusoidal ripple with wave vector $k=kxx^+kyy^$ increases exponentially in time with the rate

[The ripple amplitude decays exponentially if $σ(k)$ is negative.] It follows that ripples with wave number $1/2$ and with their wave vector along the $x$ direction emerge shortly after the growth of the thin film begins.

If the surface height $u$ does not depend on the transverse coordinate $y$⁠, Eq. (15) reduces to

which we will refer to as the one-dimensional (1D) EOM. For the case $ψ=0$⁠, we set $ϕ≡ux$⁠. Differentiating Eq. (17) with respect to $x$ yields

which is the 1D Cahn-Hilliard (CH) equation.¹² This leads us to the conclusion that for long times the solution to Eq. (17) with $ψ=0$ tends to a state in which most of the surface has a slope nearly equal to one of the two selected values $±1$⁠. The slope changes rapidly in interfacial regions between adjacent intervals where the slope is nearly constant. As some regions of nearly constant slope contract and then disappear and others expand, the pattern coarsens. For $ψ=π/2$⁠, on the other hand, Eq. (17) is the CKS equation, which has been studied as a model of amorphous thin film growth¹¹ and the step meandering instability on a crystal surface.¹³ A family of periodic steady-state solutions to the CKS equation exists;¹¹ these consist of concave, nearly parabolic segments that meet at “kinks.” These kinks are not discontinuities in $ux$ but are instead relatively narrow regions where $uxx$ is negative. For generic nominally flat initial conditions, a nearly periodic pattern with the linearly selected wave number emerges at early times, but at longer times coarsening occurs: kinks merge and the average size of the parabolic segments grows in time.¹³ 

Let us consider the behavior of solutions to the two-dimensional (2D) EOM (15) with $ψ=0$ on the domain in which $0≤x≤L$ and $0≤y≤L$ and apply periodic boundary conditions. We introduce the effective free energy

where

will be referred to as the effective potential. Equation (15) can be written

where $δF/δu$ denotes the variational derivative of $F$ with respect to surface height $u$⁠. Equation (21) implies that $dF/dt≤0$⁠, i.e., the effective free energy can never increase. The dynamics therefore tends to minimize the value of $F$⁠. The effective potential $f$ has minima at $(ux,uy)=(±1,0)$⁠. Therefore, the surface will tend toward a state in which most of the surface has a gradient $∇u$ nearly equal to $±x^$⁠, i.e., the surface will facet. A flat facet with one of the two selected gradients $±x^$ has a free energy equal to zero. Adjacent facets are separated by “edges” that have a positive free energy per unit length. $∇u$ changes rapidly but not discontinuously as an edge is traversed.

The EOM (15) cannot be written in the variational form (21) for $0<ψ≤π/2$⁠. Nevertheless, we will see that the surface tends to facet for $0<ψ<π/2$⁠. The selected gradients $∇u$ will be shown to be equal $±(secψ)x^$ for $0≤ψ<π/2$⁠, and so their common magnitude increases with $ψ$⁠. We will also see that there are no selected slopes for $ψ=π/2$⁠. Our first step toward these conclusions will be to study steady-state solutions to Eq. (15) that depend only on the longitudinal coordinate $x$⁠. This will be done in Sec. IV A.

We have an additional motivation for studying these steady-state solutions: our numerical integrations of the EOM (15) show that as time passes and the patterns coarsen, they grow close to steady-state solutions of increasingly long wavelength. Steady-state solutions to Eq. (15) of infinitely long wavelength are therefore of particular interest.

We will seek steady-state solutions to Eq. (17) in this subsection. Setting $ut=0$ in Eq. (17), recalling that $ϕ≡ux$⁠, and integrating the resulting ordinary differential equation (ODE) with respect to $x$⁠, we obtain

where $C$ is a constant of integration. We will look for steady-state solutions to Eq. (17) that are even functions of $x$⁠, which means that $ϕ$ must be an odd function of $x$⁠. Applying this condition to Eq. (22), it follows that $C=0$⁠. Then, making the replacements $x→t$ and $ϕ→x$ in Eq. (22), we obtain

where the dots denote time derivatives. For $ψ=0$⁠, Eq. (23) is the EOM for an undamped, undriven Duffing oscillator.

For $ψ=0$⁠, Eq. (23) can be written

where

Equation (24) is formally identical to Newton’s second law for a particle with unit mass moving in potential $V(x)$⁠. The potential has global maxima at $x=±1$ and a local minimum at $x=0$⁠. Equation (24) has the first integral

where the constant of integration $E$ can be thought of as the energy of the particle. For $0<E<1/4$⁠, the particle oscillates about the origin between the classical fixed points at $x=±x0$ that are given by the equation $V(x0)=E$⁠. For $E=1/4$⁠, Eq. (24) has the solutions

where $t0$ is an arbitrary real constant. In terms of original variables, Eq. (27) is

where $x0$ is an arbitrary constant. Equation (28) is the well-known solution to the 1D CH equation that describes the interface between regions in two different phases. Here, it gives the form of the surface between two regions in which the slope is constant.

For $ψ=π/2$⁠, we set $ut=0$ in Eq. (17), which yields

Integrating this equation twice gives $−u−uxx+ux2=k1x+k2$⁠, where $k1$ and $k2$ are constants of integration. We once again seek steady-state solutions $u=u(x)$ to Eq. (17) that are even functions of $x$⁠, and so $k1=0$⁠. Moreover, by setting $u¯=u+k2$ and then dropping the bar, we obtain

Making the replacements $u→y$ and $x→t$ gives

This can be recast in a form that is formally identical to Newton’s second law using a change of dependent variable due to Raible et al.¹¹ Setting $z=e−y$⁠, Eq. (31) becomes $z¨=−zlnz$ or

where the effective potential $V$ is given by

$V(z)$ has a global minimum at $z=1$⁠, a local maximum at $z=0$⁠, and tends to infinity for $z→∞$⁠. Equation (32) has the first integral

where the constant $E$ can be thought of as the energy of a particle of unit mass moving in potential $V$⁠.

For $E=0$⁠, Eq. (32) has the solution

where $t0$ is an arbitrary constant. In terms of the original variables, Eq. (35) is

which is a parabolic surface profile. Naturally, if we translate this surface profile up or down, we still have a steady-state solution to Eq. (17).

For $−1/4<E<0$⁠, the particle oscillates between the classical fixed points $z0$ and $z1$ that are given by the equation $E=V(zi)$ for $i=0$ and 1. We choose $z0$ to be the left classical fixed point, so that $0<z0<1<z1$ and $z0→0$ as $E→0−$⁠. The period of the oscillations increases as $E$ tends to zero from below and becomes infinite in the limit that $E=0$⁠. For $z$ close to $z0$⁠, we have $z¨≅−V′(z0)$⁠. Therefore, choosing the zero of time so that $z(0)=z0$⁠, we obtain

In terms of the original variables, a temporally periodic oscillation in $z$ is a spatially periodic steady-state solution to Eq. (17). Equation (37) becomes

where $c$ is an arbitrary constant. Equation (38) gives the form of the surface close to one of the wave crests. The radius of curvature at the crest is $|ln⁡z0|−1$⁠, which tends to zero in the limit that $z0→0$⁠, i.e., the limit in which the wavelength of the surface ripple tends to infinity. Thus, the ripple crests become increasingly sharp as the wavelength of the periodic steady-state solution grows large. An analogous analysis shows that the radius of curvature at the bottom of a trough tends to 2 as the ripple wavelength tends to infinity. In the limit that the wavelength becomes infinite, the surface takes on the parabolic form given by Eq. (36).

Equation (23) cannot be written in the form of Newton’s second law for $0<ψ<π/2$⁠. To make progress for any $ψ∈[0,π/2]$⁠, we set $x1=x$ and $x2=x˙$⁠. This gives us a system of two first-order ODEs,

and

A temporally periodic solution to this system is equivalent to a spatially periodic, steady-state solution to Eq. (17).

The system of equations given by Eqs. (39) and (40) defines a flow in the phase space with coordinates $x1$ and $x2$⁠, as illustrated in Fig. 2 for three values of $ψ$⁠. The flow has three fixed points for $0≤ψ<π/2$ which appear at $(x1,x2)$ = $(0,0)$ and $(±sec⁡ψ,0)$⁠. For $ψ=π/2$⁠, there is a single fixed point at $x1=x2=0$⁠. The phase plane has a region $R$ in which all of the trajectories are closed; outside $R$ the trajectories are unbounded. The boundary of $R$ is the separatrix.

Linearizing the systems (39) and (40) about $x1=x2=0$⁠, we find that the trajectories about this fixed point are approximately circular: $x12+x22≅A2$⁠, where $A$ is a small, positive constant. The corresponding solutions to Eq. (17) have the form $u=Acos(x−x0)$⁠, where the amplitude $A$ of the sinusoidal surface disturbance is small and $x0$ is an arbitrary constant.

Let us confine our attention to $ψ<π/2$ for the moment. The fixed points with $(x1,x2)=(±sec⁡ψ,0)$ are saddle points and are on the separatrix. Equations (39) and (40) have the solutions given by

and

where

and $t0$ is an arbitrary constant. These are heteroclinic trajectories that start on one of the fixed points at $(x1,x2)=(±sec⁡ψ,0)$ and end on the other. Eliminating $t$ from Eqs. (41) and (42), we obtain

the equation describing the separatrix. Equation (44) shows that the separatrix consists of two parabolic arcs, as shown in Figs. 2(a) and 2(b).

Inside the separatrix, there are nested closed orbits. The period of these orbits is finite, but if we look at orbits that are progressively closer to the separatrix, the period diverges. It takes infinitely long to traverse either of the heteroclinic trajectories that travel from one of the fixed points to the other.

In terms of the original variables, Eq. (41) is

where $x0$ is an arbitrary constant. The slope $ux$ tends to the constant values with magnitude $sec⁡ψ$ for $x→±∞$ and it switches sign at $x=x0$⁠. Note that for $ψ=0$⁠, Eq. (45) reduces to the result we obtained earlier, Eq. (28). Integrating Eq. (45) with respect to $x$⁠, we obtain

Equation (44) is

when expressed in terms of original variables.

Equation (46) describes the curved surface that separates two regions in which the slope $ux$ differs and takes on the constant values $±sec⁡ψ$⁠. If the surface is convex (concave), the width of the curved region is $l+≡1/k+$ (⁠$l−≡1/k−$⁠). The widths $l+$ and $l−$ are both equal to $2$ for $ψ=0$⁠. If $ψ>0$⁠, then $l−>l+$ and hence the concave curved region is wider than the convex curved region, as seen in Fig. 3. This lack of symmetry under transformation $u→−u$ stems from the presence of CKS nonlinearity in Eq. (17), which breaks the $u→−u$ symmetry. $l+$ is a decreasing function of $ψ$ and tends to zero for $ψ→π/2$⁠, whereas $l−$ is increasing and it tends to infinity in this limit.

For $ψ→π/2$⁠, one of the parabolic arcs that make up the separatrix reduces to $x2=1/2$ and the other goes off to infinity. Thus, for $ψ=π/2$⁠, the separatrix is the line $x2=1/2$⁠. Equation (39) shows that along this line $x˙1=1/2$ and so $x1=(t−t0)/2$⁠. In terms of the original variables, this is $ux=(x−x0)/2$⁠. Integrating this, we obtain

where $k3$ is an arbitrary constant. This is precisely the parabolic surface profile we found earlier for the case $ψ=π/2$⁠.

We integrated the 1D EOM (17) numerically for a selection of $ψ$ values using the fourth-order Runge–Kutta exponential time-differencing (ETD4RK) method described by Cox and Matthews.¹⁴ Simulations were carried out on the interval $0≤x≤L$⁠, and periodic boundary conditions were employed. The initial condition was spatial white noise with an amplitude of 0.001. The linear terms in the EOM were calculated in Fourier space and the nonlinear terms were computed in real space using central finite differencing accurate to fourth order in the grid spacing.

As time passes, the pattern coarsens and the radii of curvature of the peaks decrease. In our numerical work, it was important to ensure that the chosen grid spacing was always relatively small compared to the width of the peaks. The peak widths are on the order of half the linearly selected wavelength at early times. The average peak width decreases as time passes, however, and tends to $l+$ at long times. We therefore chose the spacing for the spatial grid $Δx$ to be approximately equal to $l+/10$⁠. Because $l+$ vanishes for $sin⁡ψ=1$⁠, we were not able to carry out reliable simulations for $sin⁡ψ$ values larger than about 0.8.

Figures 4(a)–4(c) show the results of a simulation of Eq. (17) with $ψ=0$⁠. At early times, the ripple wavelength is close to the linearly selected wavelength $22π$⁠, and there is little evidence of a selected slope. The ripple pattern becomes progressively more faceted as time passes. The average width of the facets grows and the number of interfacial regions that separate facets decreases. This leads to a reduction in the effective free energy for the surface,

as shown in Fig. 5. As expected, the slope distribution has two pronounced peaks at $ux=±1$ at long times, as shown in Fig. 6(a) for $t=1500$⁠. The wavelength and amplitude of the ripples vary with position—see Fig. 4(c) for example. This disorder in the pattern is not surprising, since the band of linearly unstable wave numbers is not narrow: Eq. (16) shows that it extends from $kx=0$ to $kx=1$⁠.

Even though Eq. (17) cannot be written in the variational form (21) for $0<sin⁡ψ≤1$⁠, the behavior we found for $0<sin⁡ψ<1$ is in some respects similar to what we observed for $sin⁡ψ=0$⁠. Figures 4(d)–4(i) show the time evolution of the surface for $sin⁡ψ=0.4$ and $sin⁡ψ=0.7$⁠. These figures show that the pattern coarsens with time, just as it does for $sin⁡ψ=0$⁠. In addition, the slope distribution once again has pronounced peaks at long times [see Figs. 6(b) and 6(c)], showing that the surface facets.

For $ψ<π/2$⁠, the solution to Eq. (17) grows increasingly close to the steady-state solutions with infinite wavelength as the time grows long. This is illustrated in Fig. 7. In this figure, $uxx$ is plotted vs $ux$ at time $t=1500$ for three different values of $ψ$⁠. Also included in the plots are the separatrices given by Eq. (47). The agreement is excellent.

The similarity between the form of the solutions to Eq. (17) at long times and the steady-state solutions is further reflected in the agreement between the predicted and measured selected slopes. In our simulations, the selected slope $ux∗$ was taken to be the positive value of $ux$ where the peak in the slope distribution appears. This is compared to the predicted value $sec⁡ψ$ in Fig. 8. Once again, the agreement is impressive.

In Sec. IV A, we found that the peaks on the surface are sharper than the troughs for $ψ>0$⁠. This is seen in Figs. 4(f) and 4(i), for example. In contrast, for $ψ=0$⁠, the radii of curvature of the peaks and troughs are the same, as illustrated in Fig. 4(c).

During the time evolution of the surface, there are long intervals in which little change occurs. These are punctuated by brief periods in which the surface changes rapidly and the number of wavelengths in the spatial domain $0≤x≤L$ decreases. For $ψ=0$⁠, this was already evident in Fig. 5. A free energy cannot be defined for $0<ψ≤π/2$⁠, but we can plot $I1D(t)≡∫0Lut2dx$ as a function of time to provide an illustration of the punctuated time evolution, as we have done in Fig. 9 for three representative values of $ψ$⁠. Each spike in these plot coincides with a decrease in the number of wavelengths in the domain by one.

To compute the mean wavelength $Λ$ in the simulations, we found the average inter-peak and inter-trough separation distances, averaged them over 100 runs, and plotted the results as a function of time for various values of $ψ$⁠. We also used the root-mean-square deviation of the surface height from its mean value $w$ as a measure of the interface width. The results are shown in Fig. 10.

The 1D CH equation produces coarsening, but at long times the coarsening is very slow—it is logarithmic.¹⁵ For the CKS equation, in contrast, $Λ∼t1/2$ and $w∼t$ at long times.¹³ Accordingly, we would expect that increasing $ψ$ would increase the rate that coarsening occurs. Our simulations show that the wavelength and the interface width do indeed grow more rapidly when $ψ$ is increased, as seen in Fig. 10. Our simulations do not give compelling evidence for a regime in which $Λ$ and $w$ display power-law scaling. In fact, the coarsening may be logarithmic for $0<ψ<π/2$⁠, just as it is for $ψ=0$⁠. This type of logarithmic behavior would be very difficult to discern in simulations.

We carried out simulations of the 2D EOM (15) for a selection of $ψ$ values, again using the ETD4RK method. As in the 1D case, we adopted periodic boundary conditions and used spatial white noise with an amplitude of 0.001 as the initial condition. The grid spacing in the $x$ direction was chosen as before. Because there is no instability in the $y$ direction, a coarser grid spacing could be used along that direction with a negligible loss of accuracy.

The behavior in 2D is in most respects similar to that in 1D. Simulations for three different values of $ψ$ are shown in Fig. 11. For $0≤ψ<π/2$⁠, the surface is faceted at sufficiently long times and the slope distributions have sharp peaks, as shown in Fig. 12. The selected slopes are once again close to the predicted value $sec⁡ψ$—see Fig. 8.

Cuts of the surface along the $x$ direction are shown in Figs. 12(d)–12(f) for $t=1500$⁠. As before, the amplitude and wavelength of the surface ripples vary with $x$⁠, and, for $ψ>0$⁠, the peaks are narrower than the troughs.

At long times, the surface is again close to the steady-state solutions with infinite wavelength. Evidence for this appears in Fig. 13, where $uxx$ is plotted vs $ux$ for three values of $ψ$ at time $t=1500$ and the separatrices given by Eq. (47) are also shown. Additional evidence is given in Fig. 14, where $I2D(t)≡∫0L∫0Lut2dxdy$ is plotted vs $t$ for three simulations with different values of $ψ$⁠. In each case, there is little change in the surface past time $t=1000$⁠.

In 2D, a phenomenon that is impossible in 1D occurs—annihilation of dislocations of opposite sign. An example of such an event appears in Fig. 15. This leads to faster coarsening in 2D than in 1D, as Fig. 10 shows. After an initial transient, $I2D(t)$ does not display the same type of high, sharp spikes as $I1D(t)$ does, because wavelengths in the ripple pattern can be removed by dislocation propagation and annihilation in 2D. At long times, no dislocations remain, and there is very little variation of the surface height in the transverse direction (see Fig. 11). However, as already noted, the amplitude and wavelength vary with $x$⁠, and so the patterns are not perfectly ordered.

In terms of the original, unscaled variables, Eq. (13) is

where $A≡Aϵ2=A(θ−θc)$⁠, and Eq. (28) is

Equation (51) shows that the unscaled selected slope $S$ is $(A/α)1/2$⁠. Interestingly, $S$ is independent of the coefficient of the CKS term $β$⁠. It also tends to zero as $(θ−θc)1/2$ as $θ$ approaches $θc$ from above.

The widths of the curved interfacial regions between adjacent facets are $(B/A)1/2l±≡w±$ in the original units. These widths diverge as $(θ−θc)−1/2$ as $θ$ approaches $θc$ from above. Because $l±$ depends on $ψ$ and $ψ$ in turn depends on $β$⁠, the widths $w±$ do depend on $β$⁠. In particular, $w+$ and $w−$ coincide for $β=0$⁠. The width of a convex curved region $w+$ is a decreasing function of $β$ and tends to zero for $β→∞$⁠. In contrast, the width of a concave interfacial region $w−$ increases with $β$ and tends to infinity for $β→∞$⁠.

As we have seen, the case in which diametrically opposed atomic beams are simultaneously incident on the growing thin film is much simpler than the single beam case because in the former case, the EOM must be invariant under transformation $x→−x$⁠. In an experiment, it would actually not be necessary for two atomic beams to be used. Instead, the sample could be rotated periodically through $180°$ increments about the $z$-axis while a single atomic beam was obliquely incident upon it. If the time between rotations were made sufficiently small, the effect would be essentially the same as if diametrically opposed beams were concurrently incident on the film.

In addition to carrying out experimental work on thin film deposition by a single obliquely incident atomic beam, Voronov and co-workers introduced a theory of ripple formation that they observed.¹ In their theory, the Carter–Vishnyakov model of MR⁵ was adopted. The Carter–Vishnyakov model, however, is unlikely to be correct for arbitrary local angles of incidence. In addition, Voronov et al. neglected the effect of the CKS term and considered only surface morphologies that do not depend on transverse coordinate $y$⁠.

Surface ripples have also been observed in oblique-incidence epitaxial growth of the (100) surfaces of metals.^16–20 In contrast to the problem studied in this paper, the energies of the incident atoms are too low to induce appreciable mass redistribution in oblique-incidence epitaxy, and shadowing effects appear to be responsible for ripple formation.^19,20

Equations that are closely related to the EOM we derived for the surface of the growing thin film with two incident beams (50) have been encountered in other physical contexts. For example, EOMs that describe the epitaxial growth of single crystal thin films and that include the Ehrlich–Schwoebel and CKS nonlinearities have been studied by Golubović and co-workers.⁹ Our EOM (50) is a special case of the EOM that they used to model the epitaxial growth of a (110) surface,^21–24 although they did not analyze the case in which $D>0$ and there is no instability in the transverse direction.

When a (100) GaAs surface is maintained at a temperature in excess of its recrystallization temperature and is bombarded with a normally incident argon ion beam, highly ordered ripples form.^25–27 Ou et al.²⁵ modeled the time evolution of the GaAs surface using a variant of Eq. (50) in which the anisotropic terms $−Buxxxx$ and $β∂x2ux2$ are replaced by their isotropic counterparts $−B∇2∇2u$ and $β∇2(∇u)2$⁠. The term $Duyy$ in our EOM (50) suppresses variations in the surface height in the transverse direction, and so the differences in the time evolution produced by Eq. (50) and the equation studied by Ou et al. are expected to be modest.