Height transitions, shape evolution, and coarsening of equilibrating quantum nanoislands
Mikhail Khenner





Department of Mathematics and Applied Physics Institute, Western Kentucky University, Bowling Green, KY 42101








Morphology evolution and coarsening of metal nanoislands 

is computed 
within the framework of a surface diffusion-type model that includes the effects of the electron energy confinement within the film, the charge spillage at the film/substrate interface, 
the energy anisotropy of the film surface and the surface stress. The conditions that result in large
islands with flat tops, steep edges, and strongly preferred heights are determined. A strong influence of the film height on the coarsening dynamics and final morphologies is found; the conditions leading to interrupted coarsening are highlighted. The dependence of the geometric parameters of the equilibrium island on the film height and on the island initial volume is computed.





Keywords: Heteroepitaxial ultra-thin metal films, quantum size effect, surface diffusion, coarsening, solid-on-solid wetting



Height transitions, shape evolution, and coarsening of equilibrating quantum nanoislands
    Mikhail Khenner
    December 30, 2023
========================================================================================








§ INTRODUCTION



Formation of a nanoscale islands on the surfaces of an ultrathin metal films was studied extensively, especially since it was discovered that a low temperature (<140K)
deposition of a nonwetting Ag film on GaAs(110), followed by prolonged annealing at moderate temperature (300-400K), results in the atomically flat film above a certain critical film height (thickness) <cit.>.

Soon after this discovery it was shown <cit.> that the quantum size effect (QSE), which is rooted in the quantization of the electronic energy in the direction across the 
metal film, is mainly responsible for this new “electronic" or “quantum" growth mode. 

Subsequent exploration of the quantum growth in heteroepitaxial Ag/Si(111) and Pb/Si(111) systems “revealed the presence of large
islands with flat tops, steep edges, and strongly preferred heights". (Quoted from Ref. <cit.>; 

also see the review <cit.> and references therein.) Resolved Scanning Tunneling Microscopy (STM) images of a 3D island morphologies can be seen, for instance, 
in Fig. 2 of Ref. <cit.> and in Figures 2 and 17 of Ref. <cit.>. A detailed understanding of the kinetic pathways to form such “magic"  island morphologies remains incomplete <cit.>.


Different models were proposed <cit.> to help explain the experimental findings in some metal/substrate systems (Pb/Si, Ag/AlPdMn, and Ag/NiAl, respectively). 
Ref. <cit.> presents the rate equation-type model. The islands are approximated as circular mesas with radii R_i and heights h_i, and the evolution equation for the number of atoms in the i-th island is constructed and solved numerically, assuming the diffusion of adatoms on the substrate. The QSE is incorporated through a term in the island chemical potential, in a fashion very similar to this paper.
That model, despite its apparent success at predicting the time-evolution of quantities such as the area fraction of the islands of a particular height, is not capable of tracking
neither the surface morphology as a whole, nor the morphology of the individual islands. Besides, since it is limited by its starting assumption of the circular islands, it is unclear how the deviations
from this shape would affect the model results. In Ref. <cit.> a step dynamics model is developed of the growth of a 
single multilayer island by incorporation of atoms deposited within its capture zone. The island is assumed to have a “wedding cake" morphology, i.e. i circular layers with 
the radii r_i(t). The conditions for creation of new top layers are prescribed, and the evolution of the radii are determined by the net attachment fluxes of diffusing adatoms to
the edges of the circular layers. The fluxes, in turn, are obtained by solving the appropriate deposition-diffusion equations for adatom density on each terrace. The QSE enters
rather indirectly through prescribing the enhanced adsorption energy on the top of the second layer, with the result that the growth of the 3rd layer is enhanced and the growth of the
4th layer is inhibited. The model thus emphasizes the growth of three-layer circular islands, with such height selectivity indeed dominating in the Ag/AlPdMn system. In Ref. 
<cit.> the Kinetic Monte Carlo simulations are performed of the atomistic lattice-gas model of film growth. The Density Functional Theory (DFT) calculations provide adatom
adsorption energies (which reflect QSE promoting a bilayer island formation), interaction energies and key diffusion barriers. Lastly, Monte Carlo simulations of the post-deposition evolution of an
ensemble of Pb islands are presented in Ref. <cit.>. Here, the (i,j) island is assumed in the form of a rectangle of height h_i and width j. The growth is driven by 
the competition of probabilities π_1 and π_2, where the former probability is one of an atom jump from a wetting layer onto island's sidewall (and thus the corresponding
flux increases the island's width), and the
Arrhenius probability π_2 is one of an atom jump from a sidewall onto the top of the island (thus the corresponding
flux increases the island's height). 
The cited models of the islands ensembles <cit.> are most closely related to the study in this paper. However, these models do not allow studies of the dynamics of the morphology evolution and the islands coarsening, 
since the individual island shape is pre-determined and is a constant of the model (i.e., a rectangle). Nonetheless, these models (see also Refs. <cit.>) are quite successful in predicting the 
equilibrium distributions of the islands heights and widths as functions of the parameters, such as various diffusion and attachment barriers on the island's 
sidewalls and the top surface.
Other models
are even more restrictive by the assumption of a single island of a pre-determined shape <cit.>. 

In order to enable the studies of the morphology evolution of a metal film surface without the above-mentioned constraints, 
the author recently introduced a new model <cit.> that may be seen as partially complementary to Ref. <cit.> The new model is a specification to material systems with strong QSE of the classical surface diffusion model 

that was very successfully used to study the morphology evolution in a 
semiconductor-on-semiconductor systems <cit.>. 
Such specification is primarily achieved by using the form of the film surface energy that very roughly mimics its oscillation 
with the film height due to QSE in a real metal/substrate system <cit.>. Precisely, the oscillation of the surface energy is taken as simply periodic 
and with a constant period. This assumption does not hold for many ultrathin metal films (Pb is the most extensively studied case), since in such films the surface energy often displays a complicated multi-layer periodicity
(sometimes with a superimposed beating pattern) <cit.>.
Thus the model at this stage is generic and phenomenological. Still, the numerical studies 
of the (nonlinear) islands dynamics that we present in this paper reveal many minutiae details of the islands formation, equilibration, and coarsening that 
are seen in experiment <cit.>. 

It is important that
the new model predicts the final quantized island heights and associated island morphologies based (at large) on the form of the dominant QSE 
surface energy that directly follows from the quantum theories <cit.>, with the physical parameters in this 
expression (Eq. (<ref>) below) that have clear physical meaning, and without any fitting or otherwise completely or partially ad-hoc treatment. Some other models that we are aware of <cit.> assume, tacitly or explicitly, that the diffusion along the island sidewalls is negligible in comparison with the peripheral diffusion along the edge of the island. 
This assumption helps in obtaining (or preserving) certain “magic" island heights that are detected in the experiment. The new model does not employ this assumption or similar assumptions, thus 
the “magic" island heights emerge spontaneously through the standard surface diffusion-driven dynamics - given only the initial surface morphology 
(which often does not manifest islands, rather it is composed of bumps and pits). We caution however, that the numerical values of these heights that we calculate
and compute in this paper cannot be taken literally as the approximations to their counterparts in a real system (due to the above-mentioned approximate nature of the QSE oscillation employed by the model and values of most parameters that are at best typical for most metal films). 
In fact, the “magic" island heights measured in experiment are the integer-multiples of a monolayer (ML) height, while such heights in this paper 
are estimated as the non-integer numbers of monolayers. Nonetheless, the equilibration and coarsening pathways that we compute and that lead to these heights are generally valid
and persist in any quantum system, irrespective of its particular height-selection protocol imposed by a special QSE oscillation and the material parameters.

Another remark is that the annealing temperature of <400K is sufficiently low to 
preserve QSE as the major effect <cit.>, and also sufficiently high to provide the necessary adatom mobility for the morphology evolution on a very large time scale of 24 hours or longer. In fact, in Pb films on Ni(111) the pronounced QSE was observed at even higher temperatures (up to 474K) <cit.>. This justifies the use of a continuum model.
A caveat is that it is universally accepted that continuum models are useful for thick films (i.e., with heights h>∼ 50 ML); it is less known and appreciated that such models are very successful in the 
description of the dynamics of very thin (liquid or solid) films (the dynamics in either case is governed by a fourth-order PDE). For example, in Ref. <cit.>
the dewetting of a sub-4nm thick liquid film is captured by a thin-film model (which, it is worth mentioning, is designed using the lubrication approximation; this is the analog of
thin-film scaling employed in this paper), and in Ref. <cit.> the extended Mullins-type model is used for solid films with heights as small as 1 ML to model cluster growth
during dewetting. More generally, in the mature field of solid-state dewetting the continuum models based on the Mullins-type approach deliver results that are in the semi-quantitative agreement with
experiment <cit.>; these models typically consider a retracting edge of the film, with one of the boundary conditions h=0 there. 
Other examples include modeling of films rupture <cit.>. Thus even the infinitesimally small film 
thickness is not, per se, a factor that invalidates the continuum model approach.

To frame our study in the context of current and future applications, we notice that self-organized metal nanoisland arrays are important for emerging technologies based on nonlinear optics, 
plasmonics, photovoltaics, and photocatalysis. Such arrays hold large promise in the design of spectrally selective absorbers <cit.>,
enhanced fluorescence and infrared spectroscopic devices <cit.>, polarizers and spectral filters <cit.>, and molecular detectors and biosensors <cit.>.
They are also of interest for manufacture of the next-generation solar cells <cit.>. Our modeling study is a step toward design of the film growth protocols that would result in regular arrays that are in demand by these applications.





§ THE MODEL FOR SURFACE MORPHOLOGY EVOLUTION



In this section we briefly describe the model that enables the computation of the post-deposition dynamics of metal nanoislands. Also we provide the summary of key results of the linear stability analysis (LSA) for
the trivial equilibrium (the planar film surface); for more details, see Ref. <cit.>.

The model is based on the Mullins' surface diffusion equation: 

    h_t=√(1+h_x^2+h_y^2) 𝒟∇_s^2 μ,

where h(x,y) is the height of the surface of  the thin film above the substrate, x, y the coordinates in the 
substrate plane, z be the coordinate normal to the substrate plane, ∇_s^2 
the surface Laplacian, 𝒟=Ω^2 D ν/k T the diffusion constant (Ω is the adatom volume, D the adatom diffusivity, 
ν the surface density of the adatoms, k T the Boltzmann's factor),
and μ=δ/δ h∫γ  dS the chemical potential;
here δ/δ h is the functional derivative, the integral is over the film surface, and

    γ = γ_0 + 
    γ^(QSE) + γ^(S) +γ^(Anis)

is the surface energy. In Eq. (<ref>) 
γ_0 is the nominal (constant) surface energy, and the contributions labeled “QSE", “S" and “Anis" are due to quantum size effect, the stress and the anisotropy, respectively. 






Due to crystalline anisotropy, the surface energy depends on the surface orientation n, and QSE and the stress introduce a dependence on h, i.e. γ = γ( n,h).
Then the calculation of the functional derivative leads to a four-term expression <cit.> (see also Refs. <cit.>)

    μ = μ_κ+μ_wet+μ_anis+μ_h.o.t.,

with

    μ_κ   =   γκ,  μ_wet = n_3 ∂_h γ, 
    μ_anis   =    -2n_3[(h_x h_xx+h_yh_xy)∂_h_xγ+(h_y h_yy+h_xh_xy)∂_h_yγ]-∂_x∂_h_xγ+
    ∂_y∂_h_yγ/n_3, 
    μ_h.o.t.   =   ∂_xx∂_h_xxγ/n_3+∂_xx∂_h_yyγ/n_3+
    ∂_yy∂_h_xxγ/n_3+∂_yy∂_h_yyγ/n_3,

where κ is the curvature, n_3 the z-component of n, and “wet", “anis", and “h.o.t." are the short-hand notations for wetting, anisotropy and higher order terms, respectively.

The non-interacting electron gas model <cit.> and the DFT calculations <cit.> suggest <cit.> the following generic form of γ^(QSE) in Eq. (<ref>)
for many metal films, including Ag, Pb,  Al,  Na, and Be <cit.>:


    γ^(QSE)(h) = γ_0^(QSE) + g_0 s^2/(h+s)^2cosη h-g_1 s/h+s,   g_0, g_1 > 0,

where the first (constant) term and the second term are due to quantum confinement of the electrons in a thin film, and the third term is due to the
electrons spilling out to the film/substrate interface (the charge spilling effect). 
Eq. (<ref>) results in the limit of thick film of a quantum model, meaning that N≫ 1, where N is the number of energy bands occupied by the electrons 
in the z-direction. γ_0^(QSE) = k_F^2 E_F/(80 π) ∼ 400-1000 erg/cm^2 and η = 2k_F=4π/λ_F∼ 20 nm^-1 for Ag, Au, Mg, Al and Pb. 
(Here E_F is the energy at the bulk Fermi level.)
The energies g_0 and g_1 also are proportional to k_F^2 E_F (see Ref. <cit.>). 
A small wetting length s = 0.4 ∼ 1.2 ML  (where the estimate of the ML value is based on the atomic diameter of 0.35 nm for Pb) 
prevents divergence as h→ 0 <cit.>. 

The continuous form (<ref>) is the result of the interpolation of the set of data points produced by the quantum models <cit.>, whose z-coordinates are the 
integer multiples of the monolayer height.

Thus the characteristic amplitudes decays as ∼ 1/h^2 and ∼ 1/h also are the features of the cited quantum models; see Eq. (25) in section 3.3.5, 
and section 3 in Appendix D of Ref. <cit.>, and also section IIC and Fig. 5 in Ref. <cit.>. The decay of the dimensionless second derivative of γ^(QSE)(h) 
with increasing h can be seen in Fig. <ref>(d), and the surface energy itself decays with the same rate. We remark here that QSE oscillation is persistent only in relatively 
thin metal films, as seen in Fig. <ref>(d); Ref. <cit.> cites the experimental evidence of QSE oscillation in films as thick as 30 ML.

With the choice of s and γ_0 for the height and energy scales, the dimensionless form of γ^(QSE)(h) reads:

    γ^(QSE)(H) = R_γ_0+G_0/(1+H)^2cosρ H-G_1/1+H,

where R_γ_0=γ_0^(QSE)/γ_0,  G_0=g_0/γ_0,  G_1=g_1/γ_0, H=h/s and ρ=sη. Notice that 0 < R_γ_0, G_0, G_1 < 1, and ρ∼ 10
for cited values of s and η.

We assume that the lattice mismatch between the film and the substrate (and more generally, 
the interaction energy of a film and substrate atoms)  has a negligible effect on the
stress at the surface; in other words,
a thin metal film on a substrate is considered free-standing for the purpose of calculating the stress at the surface. Thus the only source of stress in the model is the 
intrinsic surface stress. Clearly, this is justified only for such nearly perfectly lattice-matched quantum growth systems as 
Ag/NiAl(110) <cit.>, Pb/Ni(111) <cit.> and Ag/Fe(100) <cit.> (also see the discussion in Ref. <cit.>). But using large values of the strain parameter (up to the limit of the numerical method) allows the computations to mimic strongly stressed films  due to the lattice mismatch, i.e. Ag/GaAs(110) <cit.>, Ag/Si(111) <cit.>, Pb/Si(111) <cit.> and Pb/Cu(111) <cit.>.
Following the seminal paper by Hamilton & Wolfer <cit.>, in Ref. <cit.> we derived the closed-form analytical expression for the surface stress energy.
The dimensionless form of this expression is

    γ^(S)(H)=-ΣΓ[Γ + M (H-2)]/[2 Γ + M (H-2)]^2,

where Σ=4 Γ P ϵ_*^2/(γ_0 ℓ^2). 

Also Γ (M) is the sum of the surface (bulk) elastic moduli, P the total number of atoms in the film, ϵ_* 
the strain, and ℓ the (large) lateral
dimension of the film (the same along the x and y axes).


The anisotropic part of the surface energy (<ref>) is derived naturally using a thin film scaling <cit.>. 
Thus the small parameter α=s/ℓ is introduced, and then the
dimensionless anisotropic surface energy reads <cit.>

    γ^(Anis)(H_X, H_Y, ...)    =   ϵ_γ[a_00 +  α^2(a_20 H_X^2 + α a_30  H_X^3 + α^2 a_40 H_X^4 + a_02 H_Y^2 + 
          α a_12 H_X H_Y^2 .. 
       +   .. α^2 a_22 H_X^2 H_Y^2 + 
          α^2 a_04 H_Y^4)],

where ϵ_γ < 0 is the anisotropy strength, a_ij's are the anisotropy coefficients, and X=x/ℓ, Y=y/ℓ are the dimensionless in-plane coordinates.
In this paper we limit the discussion of the anisotropy to [111] surface; for that surface  a_00=1/3, a_20=4/3, a_30=-2√(2)/3, a_40=-5/2, a_02=4/3, a_04=-5/2, a_12=2√(2), a_22=-5. <cit.>

Finally, the evolution PDE for the time-dependent film surface morphology is derived by first substituting Eqs. (<ref>), (<ref>) and (<ref>) 
into Eq. (<ref>).
Then the total energy γ is substituted into the system (<ref>), (<ref>) - (<ref>), where the adimensionalization is completed using 
the time scaling t=(ℓ^2/D)t̂ and the height and in-plane coordinates scalings as shown above. Lastly, a thin film scaling is introduced and the equation is expanded in the powers of α, with the dominant 
contributions retained <cit.>. This results in the final PDE that is shown in the Appendix.

Primarily for reference, in Fig. <ref>(a-c) the perturbation growth rate ω from the one-dimensional (1D) linear stability analysis of the PDE is shown for the case of
the isotropic surface energy and zero stress (ϵ_γ=Σ=0); in the Figure, k is the perturbation wavenumber and H_0
the height of the planar film immediately after the low-temperature deposition (the base height). (The expressions for ω in 1D and 2D cases are listed in Ref. <cit.>) From  Fig. <ref>(a,b)  it can be seen that as H_0 increases, the planar film surface z=H_0 is alternately unstable or stable. 
Such multi-height stability is unusual in a continuum, PDE-based model of a solid film, and it is noted here for the first time.
Larger unstable heights are less unstable
than smaller unstable heights, since the interval of the unstable wavenumbers decreases as H_0 increases.

Manifestly, the stable (unstable) H_0 values belong to H-intervals where  
Δγ^(QSE)(H)≡ d^2 γ^(QSE)/d H^2 > 0 (d^2 γ^(QSE)/d H^2 < 0), see Fig. <ref>(d) (the thermodynamic stability/instability criterion <cit.>).
Thus we emphasize here the agreement between the linear stability analysis and the thermodynamic analysis.



















§ NONLINEAR DYNAMICS OF NANOISLANDS



In this section we numerically study the 1D island dynamics and their transient and equilibrium morphologies. 
We compare the dynamics for the situations when only the QSE is operative; the QSE and anisotropy
are operative; and all three effects, i.e. the QSE, the anisotropy, and the stress are operative.
This way, the understanding of how each physical mechanism contributes to the morphology evolution is obtained. In our model, switching the effect on or off is accomplished simply by taking a non-zero or zero value of the associated governing parameter, 
i.e. ϵ_γ or Σ.
The 1D setup allows to perform a large number of simulations, and as the result we
present, for the first time in a modeling paper, the complete rationalization of the height and morphology evolution in a small, generic, and idealized quantum system. 
This rationalization is expected to fully carry over to a more realistic 2D system. The parameters set that is chosen for computations is typical. Moderate changes of major parameters
would yield correspondingly moderate, quantitative-only changes of the results, without affecting the overall conclusions and rationalizations. For instance,
slightly changing the amplitudes G_0 and G_1 of the confinement and spilling terms in the QSE surface energy, Eq. (<ref>), would result in the correspondingly moderate 
shifts of the base heights H_0, as seen in Fig. <ref>, but the computed island morphologies, heights, volumes, and the coarsening scenarios will not be greatly affected. 



 §.§ The effects of the variation of the base height H_0



To elucidate the differences in the morphology evolution at varying H_0, here we employ the simple and unique initial condition
in all simulations: two Gaussian-shaped islands on the base surface z=H_0, that have the same constant width at half-maximum and the different (but also constant) heights.
Specifically, the initial height (measured from the base surface z=H_0) and the width-at-half-maximum 
of the left island are 5 and 8, respectively; these parameters for the right island are 7 and 8.
Thus it can be noticed that the initial island heights 
extend across several stable and unstable height intervals as plotted in Figures <ref>(a,b). 

The horizontal extent of the base surface (along the X-axis) is also kept constant at 20λ_max, where λ_max=2π/k_max is
the wavelength of the most unstable perturbation from the LSA; here k_max is the wavenumber at which the growth rate ω  
attains a maximum value. We used value of k_max from Fig. <ref>(c). The base film height H_0 is varied in a wide range (from 1.5 to 6). 
For most simulations we chose the linearly stable H_0 values listed in the caption to Fig. <ref>.
The boundary conditions for H(x,t) are periodic at X=0 and X=20λ_max.







In Fig. <ref> the typical coarsening process can be seen, until the smaller island disappears and the larger one evolves into the equilibrium shape with the flat
top surface. The height of this surface is H=6.62, which is one of the stable heights according to Fig. <ref>; the corresponding equilibrium island height
is thus H_e=6.62-5.3=1.32, which translates into roughly 1.6 ML (for Pb). Notice that this computation was done without the anisotropy and stress. 
Only for comparison, dewetting of islands is shown in panel (c) at ρ=0, i.e. in the case without QSE oscillation.  




Notice that the total initial volume (the one of the bulk film plus the islands) is conserved in Fig. <ref> as well as in other simulations.

When the anisotropy and stress are not present, the typical features of the evolution toward a single equilibrium island seen in Fig. <ref> were also observed in the computations with 
all stable H_0 values larger than 2.8. These computations, in particular, show that “if the lateral size of an unstable island ... is large enough, it will not decay, but rather evolve
to a stable height" <cit.>; indeed, in Fig. <ref>(a) the larger island develops a larger base when it evolves. However, when H_0=2.1 was chosen, the flat top surface disappeared from the island equilibrium shape, i.e. the top became rounded,
and when H_0=1.5 was chosen the “reverse" coarsening was observed - that is, the small island grew at the expense of the large island, and the final equilibrium shape 
was also rounded. Due to the absence of a wetting layer for the chosen parameters (see Fig. <ref>(e)), the atypical reverse coarsening at smaller film heights can be attributed only to a complex nonlinear interaction of the oscillatory QSE terms in the governing PDE (Eq. (<ref>) in Appendix). We confirmed that, as expected, the dynamics without QSE oscillation 
either in the absence of a wetting layer (Fig. <ref>(c)), or when it is 
present <cit.> does not show a tendency to reverse coarsening.



When the anisotropy and/or stress effects are turned on, the likelihood of an interrupted coarsening or a reverse  coarsening increases <cit.>. In particular, when the moderate anisotropy is activated, 
the reverse coarsening was found not only at H_0=1.5, but also at H_0=2.1; with the stress activated, the reverse coarsening is seen at H_0=2.8, and the interrupted coarsening
at H_0=2.1 and 1.5; and with both anisotropy and stress activated, the interrupted coarsening is seen at H_0=1.5, 2.1 and 2.8, and the reverse coarsening is seen at 
H_0=3.1 and 3.4. Fig. <ref> shows the reverse coarsening sequence for QSE+anisotropy+stress case at H_0=3.1 (the unstable base surface). 
Movie 1 (see the Supplemental Information)
shows the sequence of shape transitions leading to interrupted coarsening for QSE+anisotropy+stress case at H_0=2.8. In the movie, the small island first reaches the stable height
H=5.3 (measured from the substrate) and the flat top surface emerges, but because the large island is not stable, the small island increases its height to H=8 at the expense of the lateral shrinking; next, the large island decreases its height
until it reaches the stable value H=4.1 (and simultaneously the height of the small island increases); the large island then starts to grow again while the small island 
starts getting smaller, etc. These cycles of the slight height growth or incomplete decay seem to repeat, and the equilibrium state of a single island is not reached in a reasonable computation time.




For the cases where the single island equilibrium was achieved, Fig. <ref> shows the equilibrium island height, the width of the top facet, and the width at the 
base as a function of H_0.
In the QSE case the height of the island rapidly decreases and reaches the constant value 1.3 (∼1.56 ML for Pb) at H_0=4.1; simultaneously, W_top and W_base
also reach constant values 125 and 200, respectively. Notice that the flat top does not emerge until H_0=2.8, and then W_top first grows then decays. 
From H_0=2.8 to H_0=6 the dynamics of W_base is “in phase" with W_top.
In the QSE+anisotropy+stress case the equilibrium island does not emerge until
H_0=3.1 (for this H_0 the island is shown in Fig. <ref>), and its height monotonically reaches the same constant value as in the QSE case; the flat top emerges simultaneously with the equilibrium island and its width
tends to the same constant value as in the QSE case, but slower and in the oscillatory “tail" fashion. A similar trend can be conjectured for W_base.
Information on the geometric parameters of the islands is not presented in the published experimental papers, but the just discussed predictions of the model may be useful for
future experiments.




Notice that in Fig. <ref> the base surface recedes from its initial unstable height at H_0=3.1 to a nearest lower stable height H=2.8; 
this also happens without the anisotropy and the stress. (Similar decomposition of an initial unstable blanket coverage into a stable blanket coverage and stable islands on top
is observed in the experiment <cit.>.) Obviously, in the computation such surface retraction does not prevent the islands from evolving and coarsening.
This is attributed to either of the two factors.  First, if a small deformation at the base of the island emerges on a horizontal length scale which is 
much smaller than λ_max, then such short-wavelength perturbation quickly decays since ω<0 for small perturbation wavelengths, see Fig. <ref>(c).    
Second, if such deformation is of an unstable wavelength, then it stops growing (ω becomes zero) as soon as its amplitude reaches 
the stable film height immediately below H_0. In this case the film surface will fairly quickly recede to this new stable height, as seen in Fig. <ref>.
Since deformations of a second kind are characteristic of a nonlinear island dynamics even when H_0 is stable, the surface retraction was noted for several such 
H_0 values (H_0 = 1.5, 2.1, 2.8 and 3.4). Movie 2 (see the Supplemental Information) shows the details of the retraction process for H_0=3.4; in a process similar 
to the coarsening sequence shown in Fig. <ref>, the smaller left island
in the movie later is absorbed into the right island. 
At larger H_0 the retraction is arrested, which seems to indicate that the 

QSE again plays a role. The small satellite islands that are seen on the receded surface (Fig. <ref>) on a path to equilibrium are another hallmark of the two-islands dynamics on initial unstable base surfaces. 
In fact, in the beginning of the coarsening the size of these islands may be only slightly less than the size of either left or right island. For comparison, such satellite islands
do not emerge on the receding stable base surface, as seen in the Movie 2.







 §.§ The effects of the variation of the island volume


In this section we study the height transitions of the equilibrium island as a function of the increasing island volume, at fixed H_0 (and only for QSE case). 
Real islands have various volumes at different stages on a path to equilibrium; to our knowledge, the island height transitions due to volume variations were not systematically reported in the experimental publications.
We chose H_0=5.3 to eliminate the base surface retraction; for the example see Fig. <ref>.
The initial condition is again two Gaussian-shaped islands, which mimics the real-life situation - since in the experiment there is always multiple islands or none.

The initial height (measured from the base surface z=5.3) and width-at-half-maximum 
of the left island are 1.3 and 8, respectively, and the height of the right island is 3.2; its width at half maximum is increased from 8 to produce the
increasing volume V.

The smaller left island disappears fast, and the equilibrium height of the right island is plotted vs. 

V.



From Fig. <ref> it can be concluded that the new equilibrium island height, if measured from the substrate, is selected from the set of stable heights between 6 and 7.8 (see the caption to Fig. <ref>) every time V reaches a new critical value; notice that the transitions are sharp. It takes 83% volume growth to increase H_e from 6 to 6.6 (or from 0.7 to 1.3, as in the Figure); 71% volume growth to increase H_e from 6.6 to 7.2 (from 1.3 to 1.9); and 51% volume growth to increase H_e from 7.2 to 7.8 (from 1.9 to 2.5).

The 

QSE again exerts a strong influence through value of the film height H_0. Decreasing H_0 from 5.3 to 1.5 while keeping all other parameters constant 
(including the initial heights of both islands) resulted in the 
data points arrangement that is visually similar to Fig. <ref>, however, it takes only 60% volume growth to increase H_e from 2.8 to 3.4 (if measured from the substrate; or from 1.3 to 1.9, if measured from the base surface); 26% volume growth to increase H_e from 3.4 to 4.1 (from 1.9 to 2.6); and 36% volume growth to increase H_e from 4.1 to 4.7 (from 2.6 to 3.2). 

The slope of the linear fit to the data points for this case is 0.012, and thus the conclusion is that at smaller H_0 the height of the equilibrium island increases faster with the increase of the island volume. 








 §.§ Coarsening of a pitted surface into the “magic" islands


In this section we discuss the results of the computations of the surface coarsening, assuming that the starting base surface is “pockmarked" by irregularly shaped pits.
The depth of all pits is chosen randomly, such that in a given run some pits initially may be shallow and other pits may be deep. 

Also, the pits number density along the base surface is a constant for all runs, but their positioning and shapes are random.
This is a rather common initial condition in the experiments <cit.>. Our primary goal is to understand qualitatively the major factors that govern the outcomes of the coarsening process. We will show that there is two such factors: 
first, the stability or instability of the base surface, and second, the presence of deep pits.

In Fig. <ref> we collected the results of four representative computations. In the panels (a)-(c) the base surface is stable at H_0=5.3, 
and a few deepest pits are of increasing depth from (a) to (c). 
In the panel (a), where the pits are shallow, the coarsening was computed until two large flat-top islands were formed. They “sit" on the stable surface H=4.1 and
their top surface is at the stable level H=5.4. Thus the islands extend vertically across one stable domain with the center at H=4.7.  When the depth of the deepest
pit increased slightly from d=2.3 (panel (a)) to d=3.3 in panel (b), the coarsening rate sharply increased - notice that the single flat-top island formed by t=6× 10^7. 
The height of this island also reacted to the increase of the pit depth by stabilizing its base at a lower stable height H=2.8, making the island extend across three stable domains with 
the centers at H=3.4, 4.1, 4.7.
When the pits are even deeper (d=4.3) 

as in panel (c), the outcome of the coarsening changes dramatically. The islands now form on a stable surface
whose position roughly coincides with the tip of the deepest pit, and most important, the flat-top “magic" islands do not emerge. Instead a fairly narrow and very high islands with a curved top and sidewalls emerge and coarsen.

In panel (d) the base surface is at the unstable height H_0=3.1, and the pits are shallow. 

The flat-top islands are formed
similar to the scenarios seen in panels (a) and (b). As in panel (a), two flat-top islands 
extend vertically across one stable domain. Characteristically, the island's top surface is at the nearest (to H_0=3.1) larger stable height H=3.5 (despite that another, smaller 
stable height H=2.8 is also available; thus it is easier to increase the island height rather than to decrease it and increase the width). 

Lastly, the anisotropy and the stress accelerate coarsening; indeed, the single large island shown by the dashed line in panel (a) has been formed
by t_f=4.2× 10^8, vs. t_f=4.9× 10^9 for the QSE case. However, these factors do not affect the height of the island and the vertical positioning of its base and top surfaces.

In Fig. <ref> the quantitative measures of coarsening in Figs. <ref>(a-c) are provided. It can be seen that in these three cases 
the root-mean-square roughness first decreases fast (linearly in time) for a very short period, resulting in the formation of pre-“magic" islands with the curved sidewalls and
the top surface, as shown by the green line in Figs. <ref>(a,b). For the cases of initially shallow and moderately deep pits
(Figs. <ref>(a,b)) this is followed by the power laws dynamics until roughness reaches a constant value, which coincides with the 
emergence of four flat-top “magic" islands (Fig. <ref>(c)). Further coarsening of these islands into the single equilibrium island 
shown in Figs. <ref>(a,b) proceeds with islands not changing their heights, thus the roughness is unchanged. 


This scenario is contrasted to the case of initially deep pits in Fig. <ref>(c). In that case, after the initial decrease, the roughness grows linearly in time (Fig. <ref>(b)).
This growth is capped only by the termination of the coarsening, i.e. when a single, large “non-magic" island emerges, as seen in Fig. <ref>(c).







The results in this section point out again that the pits in the film surface serve as the initiation sites for the morphological transition, strongly 
mediated by the QSE, toward the flat-top islands <cit.>. Also we remark that these results and the results in Sec. <ref> re-emphasize and cast some new light on the importance of the QSE strength (through the film height) for the formation and dynamics 
of the metal nanoislands <cit.>.












§ CONCLUSIONS


A model that incorporates the QSE, anisotropy, and surfaces stress has been used to compute evolution to equilibrium of 
the pre-existing nanoscale metal islands and the coarsening of a severely rough (pitted) surface. 

Some of the most common qualitative 
features of the morphological evolution of metal nanoislands 
reported in the experiments on various metal/substrate systems are reproduced. The kinetic pathways in a formation and coarsening of the quantum islands are clarified.

We showed that the nonlinear dynamics of the morphology evolution is complicated and that it depends on factors such as 
the initial (base) film height H_0 after the deposition, the levels of anisotropy and stress and, for a pitted surface, whether the pits 

are deep.
The dependence on the base height is not only on a height value itself, which sets the strength of the QSE,

but also whether H_0 belongs to one 
of the stable or unstable height intervals imposed by the QSE. In most situations the surface evolves into the “magic" islands that have the shape of a truncated pyramid 
with a wide, flat top surface and a well-defined height, which is a difference of two stable base heights (again as imposed by the QSE). Interrupted 
coarsening
is often found at a small base height, and/or active anisotropy and strong stress. 
The dependence of the geometric parameters of the equilibrium island on the base height and on the island initial volume was computed. The equilibrium island height,
width of the base, and width of the flat top surface reach constant values as the base height increases. The gradual increase of the island initial volume at a constant base height leads to the abrupt increases of the equilibrium island height. 



It should be again noted that we employed approximations or assumptions in order to obtain a concise and numerically tractable model of the complicated multi-physics phenomena.
They are: a small aspect ratio of the film (height/lateral dimension), small lattice mismatch, and 1D modeling. 

The first assumption holds with overwhelming accuracy in the thin film systems such as discussed in the paper, since the film height and the lateral dimension typically are of the order of 10 ML and 10-100 microns, respectively. Our mathematical treatment merely formalizes this vast difference of the spatial scales by expanding the extended Mullins equation (<ref>) (after the adimensionalization) in powers of a small parameter, which is the above ratio, and retaining the dominant contributions. The result is the PDE that contains all physics of the original equation, but unlike the original equation, it allows not only the direct computation,
but also the stability analysis of the equilibria. As we show in the paper, that analysis, when it goes hand-in-hand with computation, is instrumental and indispensable 
for understanding the formation of “magic" islands (see the discussions of Figures <ref>, <ref> and <ref>). <cit.>

In regard to strain/stress effects, we chose in this paper to limit the consideration to the intrinsic surface stress. This is because this effect is important and was cited as the primary
source of stress in many systems of interest <cit.> and also, because a theoretical treatment so far has been insufficient.
Our model of the surface stress paves way to continuum models of surface stress-dominated morphology evolution of thin films (either metal or semiconductor).  
In a comprehensive model that is presently under development, this treatment is augmented by the traditional model of the lattice mismatch stress <cit.>, thus both sources of 
stress are considered and their combined effects on nanoscale metal islands can be computed. Since QSE strongly dominates and defines (at large) how such islands form and evolve,
the augmentation by the lattice mismatch stress is expected to only provide a numerical correction and not affect the qualitative features of the dynamics of islands' morphology.

Lastly, 1D modeling is useful to distill the key features of the 3D model, extract the dynamics that is at least semi-quantitatively correct, and perform a parametric
analysis in order to explore a range of experimental conditions, such as the varying levels of stress and anisotropy (as is done in this paper). Fine features of the islands' shape
and the corrections to coarsening rates will be in future computed with the help of a 3D model, which is unavoidably computationally intensive. Thus only a single set of the material parameters for a particular metal/substrate system can be employed by such model.

By accounting for the lattice mismatch stress, growth, and the
complicated, non-trivial QSE oscillation (which is unique to a chosen metal/substrate combination), the comprehensive model will allow computation of the morphology evolution
and coarsening of several tens of 3D metal nanoislands. 
We expect that this will further advance our understanding of these technologically important systems and help design the film growth protocols that would result in better height, 
shape and position selectivity of the nanoislands.







§ APPENDIX


The 1D version of the governing PDE reads

    H_t    =    B[G_1/1+H(H_XXXX-2/(1+H)^2H_XX-2/1+HH_XH_XXX -((1+H) H_XX-2H_X^2)H_XX/(1+H)^2+6 /(1+H)^3H_X^2). 
       -   (1+R_γ_0)H_XXXX+G_0/(1+H)^2(6-ρ ^2(1+H)^2/(1+H)^2H_XX-H_XXXX+4 /1+HH_XH_XXX.
       -   .(6 H_X^2
    -2(1+H)H_XX -(1+H)^2H_X^2 ρ ^2 )H_XX/(1+H)^2+( 6(1+H)^2ρ ^2-24 )H_X^2/(1+H)^3)
    cosρ  H
       +   G_0ρ/(1+H)^2(4/1+HH_XX+2 H_XH_XXX-(4H_X^2 -(1+H) H_XX)H_XX/1+H+((1+H)^2ρ
    ^2-18 )H_X^2/(1+H)^2) sinρ  H
       -   ϵ _γ(a_00+2a_20)H_XXXX-4ϵ _γa_20(2H_XXX^3+6H_XH_XXH_XXX+H_X^2H_XXXX)
       +   ΣΓ  (Γ  +(H-2) M)/(2 Γ  +(H-2) M)^2( H_XXXX -2M/2 Γ  +(H-2) MH_XX^2-4M /2Γ  +(H-2) M H_XH_XXX.
       -   .6M^2 /(2 Γ  +(H-2) M)^2H_X^2H_XX)+3 ΣΓ  (H-2)M^3/(2 Γ  +(H-2) M)^4(4M/2 Γ  +(H-2) MH_X^2- H_XX)
       +   ΣΓ M/(2 Γ  +(H-2) M)^2(H_XX^2+2H_XH_XXX-6 M^2/(2 Γ  +(H-2) M)^2H_X^2 +M/2
    Γ  +(H-2) MH_XX.
       -   ..4M/2Γ  +(H-2) MH_X^2H_XX)].

Here B=Ω^2 νγ_0/ℓ^2 α kT is Mullins number; 

in the computations
B is set equal to one. Also, t stands for the dimensionless time (that is, the hat over t̂ is dropped).
The terms proportional to cosρ H and sinρ H originate in the oscillatory contribution to the surface energy due to QSE (Eq. <ref>). Also we remark that
in the typical case Γ > M there is no singularity of the stress terms, as 2Γ + (H-2)M>0 at all H≥ 0.



200

Smith Smith A R, Chao K J, Qiu N and Shih C K 1996 Formation of atomically flat silver films on GaAs with a "silver mean" quasi periodicity
Science 273 226

Yu Yu H, Jiang C S, Ebert Ph, Wang X D, White J M, Niu Q, Zhang Z and Shih C K 2002 Quantitative determination of the metastability of flat Ag overlayers on GaAs(110)
Phys. Rev. Lett. 88 016102

Zhang Zhang Z, Niu Q and Shih C K 1998 “Electronic growth" of metallic overlayers on semiconductor substrates
Phys. Rev. Lett. 80 5381

Ozer Ozer M M, Wang C Z, Zhang Z and Weitering H H 2009 
Quantum size effects in the growth, coarsening, and properties of ultra-thin metal films and related nanostructures
J. Low Temp. Phys. 157 221-251

Han1 Han Y, Unal B, Jing D, Thiel P A, Evans J W and Liu D J 2010 
Nanoscale “quantum" islands on metal substrates: Microscopy studies and electronic structure analyzes
Materials 3 3965

UFTE Unal B, Fournee V, Thiel P A and Evans J W 2009 Structure and Growth of Height-Selected Ag Islands on Fivefold i-AlPdMn
Quasicrystalline Surfaces: STM Analysis and Step Dynamics Modeling Phys. Rev. Lett 102 196103

HUQJJLTE Han Y, Unal B, Qin F, Jing D, Jenks C J, Liu D-J, Thiel P A and Evans J W 2008 Kinetics of Facile Bilayer Island Formation at Low Temperature: Ag/NiAl(110)
Phys. Rev. Lett. 100 116105

YJEWWNZS Yu H, Jiang C S, Ebert Ph, Wang X D, White J M, Niu Q, Zhang Z and Shih C K 2002 Quantitative Determination of the Metastability of Flat Ag Overlayers on GaAs(110) 
Phys. Rev. Lett. 88 016102 

Li Li M, Evans J W, Wang C Z, Hupalo M, Tringides M C, Chan T L and Ho K M 2007 Strongly-driven coarsening of height-selected Pb islands on Si(111)
Surface Science 601 L140-L144

Kuntova Kuntova Z, Chvoj Z, Tringides M C and Yakes M (2008) Monte Carlo simulations of growth modes of Pb nanoislands on Si(111) surface Eur. Phys. J. B
64 61–66

Kuntova1 Kuntova Z, Hupalo M, Chvoj Z and Tringides M C (2007) Bilayer-ring second-layer nucleation morphology in Pb/Si(111)-7x7
Phys. Rev. B 75 205436

Kuntova2 Kuntova Z, Tringides M C and Chvoj Z (2008) Height-dependent barriers and nucleation in quantum size effect growth Phys. Rev. B 78 155431

JEM_Khenner Khenner M 2017
Interplay of quantum size effect, anisotropy and surface stress shapes the instability of thin metal films
J. Eng. Math. 104 77-92






SDV Spencer B J, Davis S H and Voorhees P W 1993 Morphological instability in epitaxially strained dislocation-free solid films: Nonlinear evolution Phys. Rev. B
47 9760-77 





Chiu Chiu C H 2004 Stable and uniform arrays of self-assembled nanocrystalline islands
Phys. Rev. B  69 165413





GolovinPRB2004 Golovin A A, Levine M S, Savina T V and Davis S H 2004
Faceting instability in the presence of wetting interactions: A mechanism for the formation of quantum dots
Phys. Rev. B 70 235342





Korzec Korzec M D, Münch A and Wagner B 2012
Anisotropic surface energy formulations and their effect on stability of a growing thin film
Interfaces and Free Boundaries 14 545-567

Korzec1 Korzec M D and  Evans P L 2010
From bell shapes to pyramids: A reduced continuum model for self-assembled quantum dot growth
Physica D 239 465-474

KTL Khenner M, Tekalign W T and Levine M 2011 
Stability of a strongly anisotropic thin epitaxial film in a wetting interaction with elastic substrate
Eur. Phys. Lett. 93 26001

LLRV Li B, Lowengrub J, Ratz A and  Voigt A 2009
Geometric evolution laws for thin crystalline films: modeling and numerics
Commun. Comput. Phys. 6  433-482

Aqua Aqua J-N and Frisch T 2010 Influence of surface energy anisotropy on the dynamics of quantum dot growth
Phys. Rev. B 82 085322

Hirayama Hirayama H 2009 Growth of atomically flat ultra-thin Ag films on Si surfaces
Surf. Sci. 603 1492-1497

Han Han Y and Liu D J 2009 
Quantum size effects in metal nanofilms: Comparison of an electron-gas model and density functional theory calculations
Phys. Rev. B 80 155404

Wu Wu B and Zhang Z 2008 Stability of metallic thin films studied with a free electron model
Phys. Rev. B 77 035410

OPL Pierre-Louis O 2016 Solid-state wetting at the nanoscale Progress in Crystal Growth and Characterization of Materials 62 177-202

CHPPM Calleja F, Hinarejos J J, Passeggi M C G Jr, Vazquez de Parga A L and Miranda R 2007 Thermal stability of atomically flat metal nanofilms on metallic
substrates Appl. Surf. Sci. 254 12-15

Golubovic Constantinescu A, Golubovic L and Levandovsky A 2013 Beyond the Young-Laplace model for cluster growth during dewetting 
of thin films: Effective coarsening exponents and the role of long range dewetting interactions Phys. Rev. E 88 032113

BGZP Bollmann T R J, van Gastel R, Zandvliet H J W and Poelsema B 2011 Anomalous decay of electronically stabilized lead mesas on Ni(111)
Phys. Rev. Lett. 107 136103

Becker Becker J, Grun G, Seemann R, Mantz H, Jacobs K, Mecke K R and Blossey R 2003 Complex dewetting scenarios captured by thin-film models
Nature Materials 2 59-63

Cheynis Cheynis F, Leroy F and Muller P  2013 Dynamics and instability of solid-state dewetting  C. R. Physique 14 578-589

Srolovitz 
Wang Y, Jiang W, Bao W Z and Srolovitz D J (2015) Sharp interface model for solid-state dewetting problems with weakly anisotropic surface energies
Phys. Rev. B 91 045303

Vaynblat Vaynblat D, Lister J R and Witelski T P 2001 Rupture of thin viscous films by van der Waals forces: Evolution and self-similarity 
Phys. Fluids 13 1130-1140 

Held Held M, Stenzel O, Wilbrandt S, Kaiser N and Tunnermann A 2012 Manufacture and characterization of optical coatings with incorporated copper island films
Appl. Optics 51 4436-4447

Aslan Aslan K, Malyn S N and Geddes C D 2008 Angular-dependent metal-enhanced fluorescence from silver island films
Chem. Phys. Lett. 453 222-228

Jensen Jensen T R, van Duyne R P, Johnson S A and Maroni V A 2000 Surface-enhanced infrared spectroscopy: A comparison of metal island films with discrete and nondiscrete surface plasmons Appl. Spectroscopy 54 371-377

Heger Heger P, Stenzel O and Kaiser N 2004 Metal island films for optics Proc. SPIE 5250 21

Anker Anker J N, Hall W P, Lyandres O, Shah N C, Zhao J and Van Duyne R P 2008 Biosensing with plasmonic nanosensors 
Nature Materials 7 442–453

Liao Liao H , Nehl C L, Hafner J H 2006 Biomedical applications of plasmon resonant metal nanoparticles Nanomedicine 1 201-208

Santbergen  Santbergen R, Temple T L, Liang R, Smets A H M, van Swaaij R A C M M and Zeman M 2012 Application of plasmonic silver island films in thin-film silicon solar cells J. Optics 14 024010



MQA Man K L, Qiu Z Q and Altman M S 2010 Quantum size effect driven thermal decomposition of Ag films on Fe(100) in the presence
of pinhole-growth morphological defects Phys. Rev. B 81 045426

HW Hamilton J C and Wolfer W G 2009 Theories of surface elasticity for nanoscale objects
Surf. Sci. 603 1284-1291

WettingL A wetting layer can be introduced by increasing G_0 while keeping G_1=0.1; say, G_0=1/15 results in a wetting layer of the dimensionless height one.
This follows from the expression for the critical film height <cit.>,  H_0c=3G_0/G_1-1, above which all heights are linearly unstable and below which all heights are stable.
At the chosen value G_0=0.01 for Fig. <ref> and other figures in the paper, H_0c<0 and thus all film heights are unstable - as shown in Fig. <ref>(e).




FTN This approach, which bears
the names longwave expansion, thin film scaling, or lubrication approximation has been fruitfully used in the theoretical/modeling thin film community for decades; see Refs. <cit.> for some of the first applications to solid and liquid films, respectively.

WD Williams M B and Davis S H 1982 Nonlinear theory of film rupture
J. Colloid and Interface Sci. 90 220-228





