www.nature.com/scientificreports

OPEN
Received: 23 August 2016 Accepted: 20 March 2017 Published: xx xx xxxx

Distinguishing physical mechanisms using GISAXS experiments and linear theory: the importance of high wavenumbers
Scott A. Norris   1, Joy C. Perkinson2, Mahsa Mokhtarzadeh3, Eitan Anzenberg3, Michael J. Aziz   2 & Karl F. Ludwig3,4
In this work we analyze GISAXS measurements of the structure factor of Si surfaces evolving during 1 keV Ar+ ion bombardment. Using newly-developed methods sensitive to the full range of experimentally-available wavenumbers q, we extract the linear amplification rate R(q) governing surface stability over a range of wavenumbers 4–5 times larger than has previously been obtained. Comparing with theoretical models also retaining full wavenumber-dependence, we find an excellent fit of the experimental data over the full range of irradiation angles and wavenumbers. Moreover, the fitted parameter values represent experimental evaluation of the magnitudes of most physical mechanisms currently believed to be important to the pattern-formation process. In all cases, the extracted values agree well with direct observations or atomistic simulations of the same quantities, suggesting that GISAXS analysis may allow more powerful comparison between experiment and theory than had previously been thought.

Energetic ions and atoms play an important role in a wide variety of surface and thin film processes including plasma etching, sputter deposition, and ion beam assisted deposition. At low energy irradiation, where the effect of the ions is confined primarily to the surface, a variety of surface modification phenomena are observed, including ultra-smoothening on the one hand1, or alternatively, the spontaneous formation of arrays of nanostructures, with lengths as small as 10 nm2. Because ion bombardment is already ubiquitous in industrial settings, and is relatively inexpensive compared to other surface processing techniques, self-organized patterning by ion bombardment could enable a simple, economical means of inducing well-defined nanoscale structures in a variety of settings. Understanding the fundamental behavior of surfaces during ion bombardment is therefore a vital goal; however, despite many years of research on the topic, a predictive understanding of such structure formation has proved elusive.
The goal of most theoretical studies into ion-induced pattern formation is to understand the evolution of the surface height h(r, t) = h(x, y, t) – where x corresponds to the in-plane projection of the ion direction, y is the in-plane direction perpendicular to x. A conceptually powerful approach is to focus on very early times during pattern development, when the amplitude of undulations in the surface height are very small, and the evolution is well-approximated by a linear first-order differential equation of the form:

∂h(q, ∂t

t)

=

R(q)h(q,

t)

+

β(q,

t)

(1)

where h(q, t) is the Fourier transform of h(r, t), the linear coefficient R(q) is called the amplification factor or dispersion relation, and the constant β(q, t) is the Fourier transform of a stochastic noise. A positive R(q) at a given

1Department of Mathematics, Southern Methodist University, Dallas Texas, 75275, USA. 2Harvard School of Engineering and Applied Sciences, Cambridge Massachusetts, 02138, USA. 3Department of Physics, Boston University, Boston Massachusetts, 02215, USA. 4Division of Materials Science and Engineering, Boston University, Boston Massachusetts, 02215, USA. Correspondence and requests for materials should be addressed to S.A.N. (email: snorris@smu.edu) or K.F.L. (email: ludwig@bu.edu)

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

1

www.nature.com/scientificreports/

bombardment angle drives exponential amplification of modes of wavevector q leading to topographic instability. Conversely a negative R(q) dampens fluctuations and stabilizes modes of wavevector q, leading to smoothing3.
The linear framework represented by Eq. (1) becomes especially powerful when coupled to GISAXS measurements of surface scattering, which correlate to the surface structure factor S(q, t) when the surface structure amplitudes are small. As described in prior work4, 5, the early-time assumption that surface dynamics are governed by Eq. (1) leads to a solution governing the evolution of the structure factor of the form

S(q, t)

=

S(q,

0)

+

α 2R(q)



exp[2R(q)t]

−

α 2R(q)

(2)

where α is the amplitude of the noise. By fitting the experimentally-measured timeseries S(q, t) for each wave-

number q to Eq. (2), the dispersion relation R(q) can be reconstructed, and then compared directly to theoretical

models of early surface morphology dynamics, in principle providing a direct connection between theory and

experiment. This approach was used previously to compare the relative contributions of sputtered vs redistributed

atoms during ion bombardment4, 5.

However, despite their contributions, these initial explorations were constrained by a focus on small wave-

numbers |q|. Experimentally, ion-induced nanostructures often have characteristic wavelengths λ larger than the

stthhueigcglkiemnsetisstesadhc0rhaoanrfgaacettohefriniQstaimmc edoairmnpshetnohsuaisot ntfihllemesssttrnhuuactmtduberevereflQaoc p=tos rha0tSqo(pqth,atan)t

irradiated crystalline semiconductor. This naturally is limited in its magnitude. In GISAXS experiments, has the greatest intensity and variation at small wav-

enumbers, and refs 4, 5 focused on a range of q ∈ [0.1, 0.4] nm−1 considerably smaller than the full range

q ∈ [0.1, 1.5] nm−1 recorded by the detector. In theoretical models, the limited range of Q leads to a common

simplification called the longwave approximation, in which the amplification factor R(q) is expanded as a Taylor

series about q = 0. For pure materials this leads to approximations of the form

R(q) ≈ −Sx(θ)qx2 − Sy(θ)qy2 − B(θ) q 4 ,

(3)

dwihreecretioSnx,ys(,θa)nadreBi(rθr)aidsiathtieocnoaenffgicleie-dnetpoefnadneunntccoonedffiitciioennatsllfyo-rsttahbeilcizuirnvgatsuurref-adceepeennedrgeyntresularxfaatcieonremspeocnhsaeniinsmth(eutswuo-

ally presumed to be Although useful

iisnotthroepdice)v,etlhoaptmisernetqoufiprehdystoicarel ginutluaritiizoent,htehseylsotnemgwianvtehaepcparsoextihmaat tSiXo,nY

(θ) < 0. suffers from

two

very

important drawbacks in the context of experimental hypothesis testing with GISAXS. First, the approximation (3)

is quantitatively poor as Q increases. For example, in a typical experiment from recent studies on Si with 1 keV Ar+

bombardment, we observed at θ = 65° the formation of ripples with a wavenumber of about q ≈ 0.28 nm−1 4, on

a thin film with a thickness of approximately 6 nm6. The resulting value of Q ≈ 1.7 is greater than unity, and

therefore calls into question the validity of the truncated expansion (3). Second, and more perniciously, a focus

on small wavenumbers hinders the distinguishing of different physical effects. Many physical processes occur simul-

taneously during ion bombardment, including (a) the sputtering of some surface atoms7–9, (b) the relocation of

many others1, 10, 11, (c) surface diffusion9, 12, (d) the creation of stress within the film13–17, and (e) the relaxation of

strain and surface energies via viscous flow18–20. Critically, all of these processes lead to longwave linear growth

rates with the same form exhibited in Eq. (3), making an experimental estimation of the relative magnitudes of

different effects very difficult. Only at higher values of Q will different physical mechanisms be readily

distinguished.

We here present a study dedicated to the comprehensive consideration of the high-Q regime, asking (a) how

best to obtain and analyze the full range of experimentally-available data to evaluate R(Q) at high values of Q, (b)

whether models retaining full Q-dependence will provide a significantly better fit to the resulting values, and (c)

whether these new fits enable us to distinguish different contributions to the fundamental physics of ion-induced

pattern formation. Our principal findings are fourfold:

	 1.	 With careful attention to noise reduction in both the experimental and analytical stages, it is possible to extract meaningful measurements of the dispersion relation that extend to high values of Q.
	 2.	 It is only possible to fit the resulting amplification rates using theoretical models retaining full Q-dependence; therefore such models should be preferred for all tasks involving quantitative comparison with experiment.
	 3.	 Fits of better models to broader data allows evaluation of more model parameters than has previously been possible; indeed, we obtain reasonable values of the magnitudes of most mechanisms of interest.
	 4.	 There is significant evidence that ion-induced stress may play an equal or greater role in ion-induced pattern formation than competing mechanisms such as erosion and impact-induced redistribution.

These findings, which represent the first direct experimental comparison between the relative magnitude of various physical mechanisms operating during ion-bombardment of pure materials, suggest that GISAXS analysis may represent a considerably more powerful tool for the comparison of theory and experiment than had previously been thought. We therefore include as many details of our approach as possible to facilitate future efforts in this direction. Moreover, we release the software we developed to perform our analyses as an open-source library.

Results
Using the methods described below in the Methods section, we obtained our primary results in three stages. First, we performed GISAXS measurements of the structure factor S(q, t) associated with early times during the

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

2

www.nature.com/scientificreports/

Figure 1.  Results of analyzing GISAXS data over an extended range of q (squares), together with fits of that data

tiaronrgmaldeoisadθteiol∈snr[ea0tna,gi2nle0isn, og3f5fθu, l4∈l 5w,[a06v,5e2,n5u8,0m4]b5oe,vre5rd5te,hp7ee0nr,ad8ne0gn,ec8eq5(]∈couv[re−vre1ts.h)5.e,Sr1ha.on5wg]ennqmar−∈e1,([aa−n) 1de.v(2ab5lu), ae1tv.i2aol5nu]santomifoR−n1s(.qoxf,

0) for irradiation R(0, qy) for

irradiation process. Second, using a newly-developed statistical approach, we fit the parameters in Eq. (2) to the

resulting data of q than has

to extract values previously been

oofbRta(iqnxe, d0).

TanhdirRd,(0w,eqyu)sfeodr

a variety of angles of incidence, over a much wider range the extracted R(q) to fit the parameters in a composite

model containing the three most widely-studied physical mechanisms known to operate during ion-induced

pattern formation on amorphous, monotomic targets:

R(q)

=

−

γη−1(θ) 2h0

(

1

Q[sinh(2Q) − 2Q] + 2Q2 + cosh(2Q)

)

−6fACXst,rYess(θ) [1

+

2Q2

Q2 + cosh(2Q)]

−CXdi,sYpl.(θ)(1

−

exp( −

1 2

(Dq)2 )) .

(4)

In this composite model, the first term on the right-hand side is Orchard’s result for viscous flow due to surface

energy relaxation18, ond term is a result

containing the surface energy describing viscous flow due to

γ, the fluidity η−1, and amorphous film an ion-induced stress21, containing the

tfhluicxkfn, easms he0a.sTuhree

secA of

the stress generated per ion impact, and a third term is a simplification of Bradley’s

shape result

ffuonrctthioeneCffeXstc,rYetsso(fθ)iotnha-itncdauncdeidffsepr uinttethrienxg-22a,ncdonyt-adiinreinctgioannso.tThheer

schuasspeedfubenlcotwio)n, bCuXdti,sYipsl.o(θth) earnwdisaecahnareascsetenrtiisatlilcy

lengthscale D. This complete aggregate

composite model omits of known mechanisms.

surface

diffusion

(dis-

The extracted R(q), and the fitted composite models, are shown in Fig. (1), where good agreement between

experiment and model is observed. Of course, the total number of parameters in the model means that the sig-

nificance of the fit quality depends on whether or not the fitted values of the parameters match up with external

estimates where such are available, or at least with general expectations where external estimates are unavailable.

Therefore, we now turn to a more detailed consideration of the angle-dependence of the fitted parameter values.

For consistency with previously published results studying long-wavelength approximations, and comparisons

between mechanisms with the same long-wavelength form and units, we will discuss, in addition to the physical

parameters in the problem, the coefficients of the long-wave approximations of each of the three mechanisms just

described, defined via

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

3

www.nature.com/scientificreports/ Figure 2.  Comparison of (a) fitted values of the film depth, to (b) values obtained from simulation using TRI3DST. Essentially exact agreement is observed across all angles.

Figure 3. (a) Fitted values of the fluidity η−1(θ). We note the decrease toward zero as θ → 90° (as expected of a flux-driven mechanism), and the reasonable agreement between coefficient values in the x and y directions (as expected of an isotropic mechanism). (b) Comparison to simulated values of the model in Eq. (8) (scaled to show shape rather than absolute magnitude).

BOrch.(θ)

=

lim
q→0

1 q4

⋅

γ Q[sinh(2Q) − 2Q] 2ηh0 [1 + 2Q2 + cosh(2Q)]

=

γh03(θ) 3η(θ)

SXst,rYess(θ)

=

lim
q→0

1 q2

⋅

6fACXst,rYess(θ)Q2 [1 + 2Q2 + cosh(2Q)]

= 3fAh02(θ)CXst,rYess(θ)

SXdi,sYpl.(θ)

=

lim
q→0

1 q2

⋅

CXdi,sYpl.(θ)1

−

exp−

1 2

(Dq)2



=

1 2

D2CXdi,sYpl.(θ).

(5) (6) (7)

Film thickness.  In the ion-enhanced viscous flow conceptualization proposed by Umbach19 and adopted

by more recent models of stress21, 23, viscous flow is confined to a thin surface layer that is provided an enhanced

fluidity due to the energy deposited by ions in collision cascades. Sigmund’s prolate Gaussian ellipsoid approxi-

mation of energy deposition suggests that the thickness would be a decreasing function of angle24, which suggests

dtheecrfeoarsminfgotrrehn0(dθ,)fraopmpeaarroinugndin6E.5q n. m(15to).

In Fig. 2a, we plot the fitted around 3 nm as θ increases.

values

of

h0(θ),

which

do

indeed

exhibit

a

Although these values seem reasonable, they can readily be compared to simulations of ion impact using

atomistic simulations. Using TRI3DST25 via the PyCraters library26, we performed 1000 simulations of 1000 eV

Ar+ → Si at each of 9 angles. As a measure of the film depth, we take the mean plus two standard deviations of the

ion implantation distribution (see e.g., ref. 24; this approach produces very good agreement with limited exper-

imental data available in ref. 27). The results are plotted in Fig. 2b. Comparing to final results for the fitted film

depth in Fig. 2a, we find very good agreement between fit and simulation.

Fluidity.  We now turn to the effective fluidity parameter η−1. Within the ion-enhanced viscous flow conceptualization, this fluidity is assumed to be proportional to the deposited power through the ion flux. Because the flux through a plane normal to the surface decreases like cos (θ), we might expect the fitted value of η−1 to decrease with increasing θ. Moreover, because the viscous flow mechanism is isotropic, we should expect the fitted value to

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

4

www.nature.com/scientificreports/

FfAigCuXsrt,reYes4s(. θ)F.iWtteednvoatleutehseoef x(cae)ltlhenetloagnrgewemaveensttroefstshceoleaftftiecriewnittShXstt,rhYeses(oθr)et=ica3llfyA-hp0r2eCdXsitc,rYetses,dafnodrm(bs){tchoes(r2aθte),ccooesf2fi(cθi)e}n. t

be the same regardless of the direction from which it is measured. Fitted values of η−1 are shown in Fig. 3a, where

both expectations are confirmed. (We note that both properties can also be observed directly from the R(q) curves

in Fig. 1, whose slopes in the limit q → ∞ of Eq. (4) are seen to be proportional to γ/η). These findings are

intriguing, because although viscous flow is increasingly assumed to be the dominant relaxation mechanism19,

few estimates of this parameter exist in the literature.

The magnitude of our fitted values, represented by η−1(0), can be compared to the only known estimate of

fluidity, found in the supplemental material of ref. 4. There, molecular dynamics simulations in ref. 20 are extrap-

olated to the present irradiation environment, obtaining an estimate of η ∼ 2.5 × 1012 Pa s. Here, because the

constants thin film:

BOrch.

and

h0

are

independent

fit

parameters,

they

can

be

used

directly

to

evaluate

the

viscosity

in

the

η(0)

=

γh03(0) 3BOrch.(0)

=

1.5

×

1011Pa sec

where we have used a surface energy of γ = 1.36 J/m2 20, 28. This is within about an order of magnitude of the scaled simulation result obtained in ref. 4; because the latter is, itself, an estimate with some uncertainty, we consider this a reasonable level of agreement.
The angle-dependence of the fluidity, represented by η−1(θ)/η−1(0), has not received much attention in the literature. We therefore briefly propose a potential model that could be informed by BCA simulations. If one assumes that the viscous relaxation occurs continuously and globally at all points in the film, then the power could simply be divided by the film thickness to obtain the average power density:

η−1(θ)

∝

cos(θ)E(θ) h0(θ)

,

(8)

winhgetrheeEe(nθ)erisgitehseodfesppoustitteerdeednaetrogmy spearnidmrpeaflcetc.tLeidkeiohn0safbroovme,tEh(eθ)inccaonmbienegstiiomnaetnederfgryo.mSiTmRuIl3aDteSdTv, ablyuesusbotfrathctisfunction are shown in Fig. 3b, showing reasonable agreement with the fitted values in Fig. 3a. This suggests interesting questions for future efforts to understand the atomistic mechanisms of ion-enhanced fluidity.

Stress Coefficient.  We next turn to the angle-dependent coefficient of stress. For comparison with prior

fsseutaxumnpdeceriteiidmosanuetansalitnifanoglreatmvhdaseliupflfoarenetridegoni-ncwttoefafodvtrehimneap,rpeppaflrr.ooa2txm1tii:mneCtgeaXstrt,rfiYeAgossrnC(oθ,Xus)wt,rpYe=essf(fAiθ{rc)s≈otinsr(3e2Fpθ×iog)r., 1t4c0tboh−.se42T(svhθ−ae)1l}.ur.WeeTsheouefnlrtSoieXsntfte,rogYertssche(u,aθwtr)vteihenmesFaaaigryger.ue v4eseaemr.tyTehnhcetoemwnn,stitoiwhsmtetehpnaerktoeewsraeyitndihstirtntehhoceett

hapapsAaarsdeaninftfeeirnxetntehtrenaanlolgntuegslawtroadfvetehpleiemsnedifteit–ntecBdeevtchaalauunseesth,twheeeplamornaagmywuesateevrethgceromoeufftpiocifpeArneCtdXsSitX,crsYett,rsYest(shsθec).ostnrteasisnisnhth0,ewahmicohrpithsoeulfsdfielpme,nwdhs iocnh

θ, it can

be compared to a direct experimental observation using wafer curvature measurements. On the one hand, ref. 21

gives the steady thin film stress in terms of fA and η, which appear a translation from fitted coefficients to implied steady stress:

in

the

forms

for

CsLtrWess

and

COLrWchard.

This

allows

T0theory

= 6fAη =

SXstr,Yess(0) BOrch.(0)

⋅

2γh0 3

≈ 0.3 GPa.

(9)

On the other hand, the stress in the film can be inferred directly from wafer curvature measurements. Under normal incidence irradiation at 1 keV, using the same ion source and flux that was used to obtain the amplification factors in Section 6, we have measured the development of a steady-state compressive stress of approximately 0.2–0.3 GPa, in very good agreement with the values inferred from the fits. We note that this value is slightly above the 0.2 GPa attributed in ref. 29 to the amorphization of the thin film.

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

5

www.nature.com/scientificreports/

Figure 5.  Comparison of coefficients of displacement effects. (a) As fitted to Eq. (14) from the values of R(q). v(inbot)leuArmvsaeclsofominrpStuhilteiecxdo-nfdroiorfmeΩctci≈roant0ea.rr0fe2uqanntumocimt3tie,oalnanrsdgimea,ufblluuaxttisooynfmsf ma≈te1t28ri×icr.rT1ahd01eia2rtecimfoioo2nnrsseeac,ntaogslfeuassceuildsiitanattgeBPvNyisCLu.ra(alNtceoortsme26pt,hauarstiisntohgnea,unpnaacrtteosrmtoaificnty some error bars are cut off).

Model Orchard only Orchard + Bradley Orchard + Norris Orchard + Norris + Bradley

χr2ed. a [nm] b [nm] d [nm] 225 5.2 9.4 —
99.6 5.9 0.0 11
71.6 0.0 9.3 —
33.2 3.3 3.5 4.4

Table 1.  Summary of fit results for different model combinations. The parameters {a, b, d} are defined below in the Methods section. The first three entries are strictly for comparative purposes; only the final model is discussed in the Results section.

Prompt displayed

Effects.  in Fig. 5a.

We As an

conclude by external test

examining the of these values,

fwitetetdurvnaltuoetshoef“tChreatdeirspFluanccetmioennFt rtearmme,wSoXdri,sYkp,l”.(dθe).vTelhoepseedatroe

predict the angle-dependent form of these coefficients directly from the results of atomistic simulation of individ-

ual ion impacts (see general results in ref. 30, and subsequent particular results in refs 31 and 32). Using the

“PyCraters” Python library26, These are plotted in Fig. 5b.

we

obtained

values

of

SXdi,sYpl.(θ)

for

the

experimental

conditions

described

above.

Comparing these panels, we see that the fitted values show reasonable agreement with the simulated values

egmxivaaegcnnt,itwthuehdiuelen(ctwheeretnfaiiotntteteidetsvhaaalntudgersseooaufterScrXdeidsspeol.v(fθina)t,oifioosnersiwnfhrhoiecmrhenmsitmeinmuolaaudtreedoavvnearilamullepasoparprteraonaatccschiom–mptphliaefnicaiagetrdieoebnmy, alearnretgoienfrtSehYdreirsorpilr.g(bθh)atriossr)nd.eeTarhroliysf

finding lends support to a hypothesis, offered in Methods below, that redistribution may have a similar form in

the wavenumber q as does erosion. that are purely negative, and larger

cInonptarribtiuctuiloanr,siftriosmknroewdinsttrhibaut StiYdoinsplt.(hθa)tiasrceopmupreolsyepdoosfiteirvoes3i1v. eBcoothntcroibnutrtiiobnus-

tions are required to produce the curves obtained from simulation in Fig. 5b, and the essentially exact agreement

between suggests

the the

psirmesuelnacteedofSYbdiosptlh.(θer)oasniodnthaendfitrteeddisvtarliubeustiionnFiing.t h5ae

(which latter.

are

less

uncertain

than

those

for

SXdispl.(θ))

Discussion
For each of the model parameters described above, we found that fits of the extracted dispersion relations to the composite model (4) produced estimates that were in reasonable to excellent agreement with available external simulations, experiments, or theoretical predictions. This finding has several implications. First, it suggests that (4) is a reasonably complete description of the dynamics of the irradiated surface. Indeed, when any one mechanism is omitted from the model, the quality of both the overall fits, and the per-parameter agreement between fits and external estimates, decreases significantly (see Table 1). This is precisely what one would expect given a model that is “as simple as possible, but not simpler.”
Second, the favorable agreement between the fitted parameter values and available external estimates (e.g. the film thickness) builds confidence in the validity of those fitted values for which external estimates are unavailable. Indeed, the fitting process we have employed here can be viewed as the experimental inference of previously-unavailable parameter values. In particular, we have obtained values of the fluidity η−1 and stress parameter A, for which no direct experimental measurement strategies currently exist. Ultimately, this approach allows us to overcome a long-standing inability, introduced above, to distinguish physical mechanisms with the same limiting form in the long-wave limit. Here, with access to amplification rates for a much wider range of wavenumbers, such distinctions have become possible.

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

6

www.nature.com/scientificreports/

A straightforward demonstration of this last point can be made in the case of the surface energy relaxation mechanism. Early theoretical models of surface structure formation invoked atomic surface diffusion as the ref. 12, which produces linearized contributions to R(q) of the form

RMullins(q) = −µγ q 4 ,

(10)

where γ is again the surface energy, and μ is a temperature-dependent mobility. The quartic behavior of this mechanism directly informed the quartic term in the longwave approximation (3); indeed, because the competing Orchard mechanism (11) is also quartic in the longwave limit q → 0, it has often been viewed as a “drop-in” replacement for surface diffusion at low operating temperatures. Within the present context, however, we are focused on the opposing limit q → ∞. As described above, in this limit the Orchard mechanism grows linearly in q, which is drastically different than the quartic response exhibited by surface diffusion in the same limit. Therefore, the observation in Fig. 1 of nearly linear growth in R(q) speaks strongly and directly in favor of the Orchard mechanism rather than that of surface diffusion (and is the reason we did not include surface diffusion in the composite model (4)).
A more recent, and still unresolved question on mechanisms concerns the magnitude of the stress effect relative to the prompt collisional effects of erosion and redistribution21. In ref. 31, the combination of simulated values of displacement coefficients and an estimate of the strength of surface-confined viscous flow was found to produce good agreement with experimental data on pattern wavelengths for 250 eV Ar+ → Si. In ref. 21, the combination of stress and viscous flow was also found to produce good agreement with the data, this time without the need to estimate any parameters. Unfortunately, because these two mechanisms look identical in the limit q → 0, it was not previously possible to estimate the relative magnitude of these two effects. Here, by comparing Figs 4a and 5a, we see that the magnitudes of the stress coefficient significantly exceed those of the displacement coefficient at angles below θ = 45°. At higher angles, however, erosion seems to play an important role in suppressing the instability in the x-direction associated with stress, and driving its own instability in the y-direction. This result suggests that a solid understanding of both stress and displacement effects will be required for predictive models of ion irradiation, especially near grazing incidence.
We conclude by summarizing four important findings of this work:

	 1.	 With the use of hierarchical fitting methods exploiting the assumed uncorrelated nature of system noise, it is possible to extract much more information from GISAXS analysis of surface roughening than previously reported. In particular, compared to prior studies limited to the range q ∈ [0.1, 0.4] nm−1, we have obtained values of the linear dispersion relation over the full range of wavenumbers q ∈ [0.1, 1.5] nm−1 accessible to the x-ray detector in a typical beam line.
	 2.	 Although longwave approximations of the dispersion relation associated with ion bombardment are convenient for intuitive theoretical discussions, they are inappropriate for comparison with experimental results containing data at higher values of q. Instead, dispersion relations retaining full wavenumber dependence are required to obtain good agreement between experiment and theory. As a demonstration of this point, we found that the high-q data are entirely inconsistent with the use of surface diffusion as a surface energy relaxation mechanism, but highly consistent with the use of thin-film viscous flow.
	 3.	 The comparison of experimental data across the full range of available wavenumbers, to linear models retaining full wavenumber dependence, enables the extraction of significantly more model parameters than has previously been possible. In particular, we are able to extract reasonable values of (a) the angle-dependent thickness of the amorphous thin film, (b) the ion-enhanced fluidity within that film, (c) the magnitude of the stress being generated by the ion beam, and (d) the strength of prompt atomic displacement mechanisms. In each case the fitted values of these coefficients agreed well with external theoretical, experimental, and simulation results.
	 4.	 Most importantly, the approach we have employed enables the evaluation of absolute and relative strength of hypothesized physical mechanisms during ion bombardment. In a direct comparison of the relative magnitudes of the stress mechanism to the combined effects of erosion and redistribution, we find that the former is significantly larger at angles below about θ = 60°. Although questions remain regarding the accuracy of BCA for estimating ion-induced displacement effects33, the agreement of fits with the BCA in this case strongly supports the possibility that stress may be more important than impact-induced redistribution in the formation of patterns during low-energy ion bombardment of Si, and highlights the value of efforts such as34, 35 to understand it in more detail.

However, we believe the most important result we have obtained lies in the methodology itself. The overall quality of the composite model fits both to the underlying experimental data, and to external estimates of individual model parameters, both supports the validity of the overall approach, and also suggests that the models employed here encompass most of the relevant physical mechanisms in operation during ion bombardment. The ability to extract most model parameters of interest directly from experimental data represents a significant advance toward a comprehensive understanding of the ion-induced pattern formation process, and an increase in the value of GISAXS studies more generally.

Methods Experimental Methods.  Two sets of data were used in these studies. For both, samples were p-doped
Si(001) with resistivity 1–10 Ω ⋅ cm, affixed clipless to a molybdenum sample platen with extra care to minimize secondary collisions that can lead to sputtering of impurities onto the surface. The sample platen was mounted in

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

7

www.nature.com/scientificreports/

a custom vacuum chamber with base pressure 1 × 10−8 Torr. Samples were bombarded at room temperature with 1 keV Ar+ using a Physical Electronics Inc. PHI ion source with a beam diameter of approximately 1.5 cm.

Existing Data.  For x-ray experiments with projected incident x-ray beam direction parallel to the projected ion

beam direction on the sample (i.e., revealing corrugations perpendicular to the ion beam, which is typically called

the y-direction), we re-used measurements taken in refs 4 and 5. These experiments were performed using a

dedicated facility for the study of surface and thin film processes on beamline X21 at the National Synchrotron

Light Source (NSLS), using a photon energy of 10 keV (0.124 nm wavelength) that was selected by a Si(111) mon-

ochromator and a flux of approximately 1012 photons/s. The x-ray incidence angle on the sample surface was 0.82

from the surface tangent, and a 384 pixel linear detector measured the scattering pattern with an exit angle set to

the critical angle for silicon, 0.18. Samples were approximately 1 × 1 cm2 in area and were affixed to the sample

platen using silver paste, and irradiated with an ion flux of approximately perpendicular to the ion beam. We collected data at irradiation angles θ =

{20,×251,01425c,mio52n5sse,c

reckoned in a plane 70, 80}. More details

of the experiment can be found in refs 4 and 5.

New Data.  For x-ray measurements with projected incident x-ray beam direction perpendicular to the ion

beam direction (i.e., revealing corrugations parallel to the ion beam, which is typically called the x-direction),

new experiments were performed at the Cornell High-Energy Synchrotron Source (CHESS) on beamline G3 in a

facility that focuses on surface and thin film processes, using a photon energy of 11.2 keV (0.111 nm wavelength)

selected by a Si(111) monochromator with a 2 mm × 1 mm. The x-ray incidence angle on the

flux of sample

aspurpfraocxeiwmaaste0l.y893°.×Sc1a0tt1e3rmepmdho2pt⋅ohnsseoctoannsdwbeeraemmedaismuerendsiboynas

497 × 195 pixel PILATUS area detector. Samples were approximately 1.25 × 1.25 in2 in order to cover the entire

elevated area of the sample platen to further minimize incorporation of metallic impurities, and were affixed

using molten indium. The ion flux was approximately ion beam. We collected data at ion incidence angles θ

1.2 =

× {0,

12001,33cm0i2,on⋅5ss0ec,

reckoned 65}.

in

a

plane

perpendicular

to

the

Discussion.  Experimental work at the CHESS synchrotron was informed by the earlier experiments at the NSLS.

Analysis of NSLS data suggested that the nonlinear regime of nanopattern amplification may start at lower flu-

ences than previously expected, so experiments at CHESS focused on lower fluence nanopatterning. Before data

ocoxlildeec,tiaomno, rspamhipzleetshweesruerbfaocme blaayredre,danudndmeirnvimaciuzueminiwtiiathl tarafnluseienncteboefhaatvlieoars.tS8a.m6 p×les10t1h5aictomnh2satdo

remove native been sitting in

atmosphere had a critical GISAXS angle slightly different than that of samples having undergone significant bom-

bardment at vacuum, which we attribute to the removal of the native oxide layer. We ensured that each sample

was bombarded in vacuum for at least long enough to cause this shift in critical angle before using the sample for

data collection. To ensure amorphization, the chosen initial fluence was further selected to be at least one order of

magnitude greater than the amorphization fluences reported for 1 keV Ar+ irradiation of Si (111) and Ge (100),

which are in the range of 0.9 − 3.0 with [smoothing/roughening], then

×oxi1d0e14wicoamns2s

36. If data were to be removed at an angle

collected for an irradiation angle associated associated with [roughening/smoothing].

However, at CHESS an irregular beam intensity profile, as well as variations in beam position by as much as

200 m over the course of data collection for a single sample, resulted in a variable incident x-ray intensity that did

not correspond to any flux values measured upstream. This required scaling scattering patterns by observed cen-

tral peak intensity rather than measured incident beam intensity, as described in Section A of the Supplemental

Material. In addition, the larger sample size and wider beam profile at the CHESS facility resulted in a significant

fraction of the scattering signal originating outside the region of maximum ion beam intensity. Hence, although

the ion flux at CHESS was nominally higher than the ion flux at NSLS, our recorded intensity profiles strongly

suggest that the actual effective flux was lower at CHESS by a factor of about 3.3. To compensate for this, the time

scale of all CHESS experiments was decreased by a common factor of 3.3.

Analytical Methods.  The analytical methods used in this manuscript have been made available as an open-source library, called PyGLIDRE – Python GISAXS LInear Dispersion Relation Extraction. This library is described in some detail in Section B of the Supplemental Material; here we summarize the main background and capabilities.

Prior Work.  As described above, at early times in the pattern formation process, one may assume that the ampli-

tude of surface structures is small, in which case the measured GISAXS intensity corresponds to the surface

structure factor S(q). By performing Fourier analysis on the linear governing Eq. (1), one obtains a solution to

the evolution of the structure factor of the form Eq. (2), which is then fit to the measured timeseries S(q, t) for

each wavenumber q to reconstruct the dispersion relation R(q). If the x-ray beam is aligned in the same direction

as the ion beam, one obtained if the x-ray

obtains beam is

aalniggnlee-ddeppeerpnednednitcvuallauretsootfhRe (i0o,nqby)e,awmh.ereas

angle-dependent

values

of

R(qx,

0)

are

However, the analyses performed in refs 4, 5 were limited to a typical range of about q ∈ [0.1, 0.4] nm−1 for

most incidence angles, significantly narrower than the range q ∈ [0, 1.5] nm−1 available on the detector. The

reason is that above a certain value of q, the intensity is quite small initially, and almost immediately reaches the

me(aalicsnhoimosfmizthaelerl)stheserteekaesdpvyearrvayamlsumeetaeIlr∞lsr{=eId0,u−Rct2,αiRoα.n}Fsiosirnuantnhdeeerso-sbdejenettceitraimvlleyinfcueondnc,statinaondntr(dessaeutealtbssleuinte,

the problem of independently fitting very large variations in R and α as the squares in Fig. 6a,b, and compare the

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

8

www.nature.com/scientificreports/

Figure 6.  An illustration and evaluation of the hierarchical fitting approach, for irradiation of Si by 1000 eV

Ar+ ions, at an angle of incidence of 20 degrees. (a) Fitted values of R(q) for local vs. hierarchical fits. (b) Fitted

values of α for local hierarchically-fitted hierarchically-fitted

vvvsaa.lluuheeiessrooaffrctthhhieecssaaalmmfitees..[((tdoc))fCaCcooilmmitapptaaerreiissvooannluooafftiddoaantta,awsslleiiccheeassvIIe((qqp0,l,ott0t))t,effodorresvseeellreeycctfteiefddthvvavalaluuleuesesoooffftqq0,,0,taotnod

displaced successive plots downwards by 0.5 units]. These results suggest that the hierarchical fit can return

much more information than a local fit, while remaining consistent with the underlying data. We also note that

the hierarchical fits of {R, α} are consistent with the local fits over the range q ∈ [0.1, 0.3] where the local fits

produce meaningful results.

erratic behavior in (5)b for q > 0.4 with the essentially flat timeseries in (6)c for q > 0.4). Consequently, prior studies have included values of R only for those values of q where both R and α could be independently obtained. Hierarchical Fitting.  To overcome this limitation, we here exploit the expected and observed properties of the noise term. If the noise β(q, t) produced by the ion gun is uncorrelated in space and time, as long assumed37, then the term α should be simply a constant, and indeed, over the range of q where the intensity varies enough to indepeaxnepdnldiRceinttoltylbyaesfisftiut{mtIe0de, dRlot,chaαalt}ly,αtfhoiseraefiactcotehndsvtvaaalnultue,eaonofdfq.αfiTttdheoidseiasssaaipcnpcgoelema,rgptllooisbbhaeeldvaawcluoitenhosatfahαnitea.rcBarroacsshseidaclaollvnfaittltuhienissgoosbftrqsae,trwevghayitlideoenasl,clowrwiebiehndagviIne0 the Supplemental Material. An example of the result is shown by the red squares in Fig. 6a,b, where we see the improvement in the range of extracted values of R(q) and α attained by the hierarchical fits, noting that those fits are consistent with the strictly local fits when the latter provide meaningful results. Then, in Fig. 6c,d we compare cross-sections of the intensity function I(q, t) to hierarchical fits of the same, finding good agreement both for selected constant values of q, and also for selected constant values of t, despite the selection of just a single value of α for all values of q. Modeling.  We now consider the comparison of the recovered values of R(q) to various theoretical models. Because our analytical approach yields values over a much wider range of q than was previously available, we can see immediately from Fig. 1 that the longwave approximation in Eq. (3) is unlikely to be a good description of R(Q) over this expanded range. Indeed, the most striking feature of our evaluations of R(q) is that they grow at most linearly as Q → ∞. It is therefore clear that Eq. (3) must be replaced with non-truncated models to attain consistency with the observed behavior of R(Q) as Q → ∞. (Indeed, attempts to fit Eq. (3) to the data return negative values of B, a non-physical result that contradicts the explicit assumptions of models generating such terms). Thin-Film viscous flow.  We begin by considering the mechanism of surface energy relaxation, which is usually assumed to dominate surface dynamics in the limit q → ∞. The linear behavior of R(q) in this limit observed in Fig. 1 strongly suggests the mechanism of ion-enhanced viscous flow19, in which the ion beam enhances the fluidity of the amorphous layer, which then enables viscous flow in the Stokes limit of high viscosity. Such flow was shown by Orchard to exhibit the dispersion relation18

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

9

www.nature.com/scientificreports/

ROrchard(q)

=

−

γ 2ηh0

(1

Q[sinh(2Q) − 2Q] + 2Q2 + cosh(2Q)

)

(11)

where η−1 is the fluidity, h0 is the average amorphous film exhibiting quartic behavior at small wavenumbers (R ∝

thickness, and Q q4 as q → 0), this

=mehc0haqn.iAsmlthboeuhgahvewsildineleyakrlnyoawt nhifgohr

wavenumbers (R ∝ q as q → ∞). We note that although the Orchard mechanism has recently been favored over

e.g. surface-diffusion due to its consistency with studies on temperature dependence of the irradiation process,

and the results of molecular dynamics studies20, Fig. 1 represents, to our knowledge, the first experimental obser-

vation of definitively Orchard-like behavior under irradiation conditions.

Beam-Induced Stress.  Because the Orchard model is unconditionally stable, whereas our data are consistent with a well-known instability at high angles of incidence, we now need to add a mechanism capable of generating an instability, for which a full-spectrum model is available. We first consider two recent models describing the effect of ion beam-induced stress (which can be in the GPa range27) on the viscous flow in the film. The first is due to Castro and Cuerno23, and has the form

RCastro/Cuerno(q)

=

fE

CXst,rYess(θ)

[sinh(2Q) Q[1 + 2Q2 +

− 2Q] cosh(2Q)]

.

(12)

fHuenrcet,ifoEndoefscarnibglees.

the characteristic strength of an angle-dependent “effective body The second is due to Norris21, and has a dispersion relation with

tfohrecseo”m(EeBwFh)a, tadnidffCerXset,rYnests(fθo)rims a

RNorris(q)

=

6fACXst,rYess(θ) [1

+

2Q2

Q2 + cosh(2Q)]

(13)

where angle.

f is As

the ion flux and discussed in ref.

A is 21,

the the

magnitude of models differ

ianntihmepwaarytesdtrsetsrsesiss-ifnretreosdtruaciend, ,alnedadaignagintoCtXsht,rYeessd(θif)feisreanftufnocrtmiosnfoorf

the first terms in Eqs (12) and (13) (the second term in each is again the Orchard relaxation rate). In principle, the

models could be distinguished experimentally through the polynomial versus exponential decay rates asQ → ∞, which arise due to different modeling assumptions. Unfortunately this distinction is more subtle than that exist-

ing between surface diffusion and surface-confined viscous flow, and because the models have the same behavior

for small values of Q, both yield about the same fit of the existing data. For the reasons discussed in ref. 21, we will use Eq. (13) to perform the fits.

Prompt appears

Atomic Displacements.  The fitted value of to asymptote to a non-zero constant as q →

R∞(q. xT,h0i)s

for θ = 80° behavior is

is unique among the data sets in not compatible with either of the

that it previ-

ous mechanisms. Suggestively, this angular and directional regime is precisely that in which erosive effects have

been presumed to be most important21, 31. Therefore, we now consider the addition of the Sigmund model of

atomic sputtering7, 8, for which the exact dispersion relation was recently obtained by Bradley22:

RBradley(qx , qy) =

ΛJεa4 2π (αβ

B1

)3

exp−

a4cos2θ 2α2β 2B1



×cos2(θ)

−

exp−

a2 2B1

qx2

−

β2 2

qy2

+

a3 sin(θ) α2B1

iqx



×cos2(θ)

−

β

2

sin(θ) a

iqx



( ) ( )wlatheeWrraelhΛsetnJrεaugissgilanegsp,rtahonpisdorrBets1iou=nltatloiαathye2clsopinnf2sit(taθon)utr+, eaxitsβratahc2etcemodsev2a(aθnlu)ieiossnaopfgeReno(qemt)r,eawttrieiocwnalidlcliosmntaasntkacenett,w.αoainmdpβoratraentthaesslounmgpittuiodnins.aFl iarnsdt,
for simplicity, we shall abbreviate this result to the following form

RBradley (q)

≈

CXdi,sYpl.(θ)1

−

exp−

1 2

(Dq)2



(14)

twiqshxth-eidecwrilreaeevncceogtnintouhsntmsa,cnobatrleseriafo(nαafdGt≠hmaeuoβcsso,stitlaohlinfesittaohrbnaebnacrsnaeitsgviclioeaantddifoeer.npoTemihns idi0snetaenobxc1abec)ratewr,vebhiauailtbteisonaonelrlvobieswerdetinhxiagnecltaeotlslistnahenetxghhflueeibnqdicyte-tspidoteihnrneeCdcseXtdanii,msoYcpnlee.(tθ(goq)e,xbn a=eenrr d0ae)lpDdirfeeisαspeae =nncthd eβdae.rnIbacnycettteohhrneefit pSaercaomnedt,earlCthXdoi,sYuplg.(hθ)r.ef. 22 considered only the Sigmund model of erosion, we here hypothesize that the accompanying collisional mechanism of mass redistribution1, 10, though exhibiting a markedly different dependence on irradiation angle4, 11, 31, nevertheless exhibits a similar dependence on the wavenumber q. This is done in part by necessity – existing models for redistribution are inherently longwave in nature, as no true analog of Sigmund’s

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

10

www.nature.com/scientificreports/

model for erosion exists for redistribution (see discussion in ref. 30). However, erosion and redistribution are both, ultimately, caused by the same collision cascade, are occasionally both modeled using rotating Gaussian ellipsoids (see refs 11 and 38), and the expansion of the full collision cascade in moment form reveals that each erosive term has a redistributive counterpart30, 31. Therefore the assumption of similar dependence on q has a reasonable physical and theoretical basis, in addition to the statistical support observed in Figure 5.

Fit Constraints.  The combined model exhibits effective fluidity γη−1, a measure of beam stress and a lengthscale D associated with the prompt

ffeiAvffSeeXscft,rrtYee.sseH(θpo)aw,raaemvmeeert,aewsruegrerfooouufnppsdrpothemraptatnadtgtilesemp: ltpahcteienmfiglemtnottfheitifcfeekacnctsehsSsoXdfhi,sYtp0h,l.(teθhs)ee,

dpmeainrzaoemrmteoitneparutsorirsnuodefestphmeenalldelaedndeitcnlyrgefacosoreesefafiincchitehanent rgoeufslitadhruedaOla–tracwhsaeitrthdwotuaestrmcphraocnobugleipnmlegastthtihce.erIvenlaapltuaiversetiomcfuahlga0nra,inttuhddeγeaηpo−pf1,tehaaelrlaoenwffceienctog–fthhbe0yimmn iotnhvie--

ing these and D.

parameters

in

opposing

directions.

A

similar

coupling

behavior

is

seen

between

the

coefficients

SXdi,sYpl.(θ)

To limit large fluctuations in the fitted values of these parameters, we again used a hierarchical fitting scheme,

allowing γη−1, of the form

fA,

and

C

to

be

fitted

independently

for

each

angle,

but

requiring

h0

and

D

to

be

a

global

functions

h0(θ) = a + b cos(θ).

(15)

D(θ) = d

(16)

These forms represent the leading terms of a Fourier cosine series, chosen because of assumed symmetry of these parameters in the irradiation angle. Applying this overall fitting strategy, we obtain the fitted versions of Eq. (13) displayed in Fig. 1. We see that the fitted models match quite closely the extracted values of R(q). For completeness, we have also fitted the data to various subsets of the mechanisms just described. These results are summarized in Table 1, showing that all three mechanisms contribute significantly to the minimization of modified χ2 error.

References

	 1.	 Moseler, M., Gumbsch, P., Casiraghi, C., Ferrari, A. C. & Robertson, J. The ultrasmoothness of diamond-like carbon surfaces. Science 309, 1545–1548, doi:10.1126/science.1114577 (2005).
	 2.	 Chan, W. L. & Chason, E. Making waves: kinetic processes controlling surface evolution during low energy ion sputtering. J. Appl. Phys. 101, 121301, doi:10.1063/1.2749198 (2007).
	 3.	 Cross, M. & Greenside H. Pattern Formation and Dynamics in Nonequilibrium Systems. ISBN: 0521770505 (Cambridge University Press, 2009).

	 4.	 Madi, C. S., Anzenberg, E., Ludwig, K. F. Jr. & Aziz, M. J. Mass redistribution causes the structural richness of ion-irradiated surfaces.

Phys. Rev. Lett. 106, 066101, doi:10.1103/PhysRevLett.106.066101 (2011). 	 5.	 Anzenberg, E., Madi, C. S., Aziz, M. J. & Ludwig, K. F. Jr. Time-resolved measurements of nanoscale surface pattern formation

kinetics in two dimensions on ion-irradiated si. Physical Review B 84, 214108, doi:10.1103/PhysRevB.84.214108 (2011). 	 6.	 Mutzke, A., Schneider, R., Eckstein, W. & Dohmen, R. SDTrimSP version 5.0. IPP Report 12/8, Garching (2011).

	 7.	 Sigmund, P. Theory of sputtering. I. Sputtering yield of amorphous and polycrystalline targets. Phys. Rev. 184, 383–416, doi:10.1103/ PhysRev.184.383 (1969).

	 8.	 Sigmund, P. A mechanism of surface micro-roughening by ion bombardment. J. Mater. Sci. 8, 1545–1553, doi:10.1007/BF00754888 (1973).

	 9.	 Bradley, R. M. & Harper, J. M. E. Theory of ripple topography induced by ion bombardment. J. Vac. Sci. Technol. 6, 2390–2395, doi:10.1116/1.575561 (1988).

	10.	 Carter, G. & Vishnyakov, V. Roughening and ripple instabilities on ion-bombarded si. Phys. Rev. B 54, 17647–17653, doi:10.1103/ PhysRevB.54.17647 (1996).

	11.	 Davidovitch, B. P., Aziz, M. J. & Brenner, M. P. On the stabilization of ion sputtered surfaces. Phys. Rev. B 76, 205420, doi:10.1103/ PhysRevB.76.205420 (2007).

	12.	 Mullins, W. W. Flattening of a nearly plane solid surface due to capillarity. J. Appl. Phys. 30, 77–83, doi:10.1063/1.1734979 (1959). 	13.	 Volkert, C. A. Stress and plastic flow in silicon during amorphization by ion bombardment. J. Appl. Phys. 70, 3521–3527,
doi:10.1063/1.349247 (1991).

	14.	 Snoeks, E., Weber, T., Cacciato, A. & Polman, A. MeV ion irradiation-induced creation and relaxation of mechanical stress in silica.

J. Appl. Phys. 78, 4723–4732, doi:10.1063/1.359820 (1995).

	15.	

Brongersma, M. L., Applied Physics 88,

Snoeks, E., van Dillen, T. & Polman, A. 59–64, doi:10.1063/1.373624 (2000).

Origin

of

MeV

ion

irradiation-induced

stress

changes

in

SiO2.

Journal

of

	16.	 Mayr, S. G. & Averback, R. S. Ion-irradiation-induced stresses and swelling in amorphous ge thin lms. Phys. Rev. B 71, 134102,

doi:10.1103/PhysRevB.71.134102 (2005).

	17.	 Chan, WaiLun & Chason, Eric Stress evolution and defect di usion in cu during low energy ion irradiation: Experiments and

modeling. J. Vac. Sci. Technol. A 26, 44–51, doi:10.1116/1.2812432 (2008).

	18.	 Orchard, S. E. On surface levelling in viscous liquids and gels. Appl. Sci. Res. 11A, 451 (1962).

	19.	

Umbach, C. C., Headrick, R. enhanced viscous flow. Phys.

L. & Rev.

Chang, K.-C. Spontaneous nanoscale corrugation of ion-eroded Lett. 87, 246104, doi:10.1103/PhysRevLett.87.246104 (2001).

SiO2:

The

role

of

ion-irradiation-

	20.	 Vauth, S. & Mayr, S. G. Relevance of surface viscous flow, surface di usion, and ballistic effects in kev ion smoothing of amorphous

surfaces. Phys. Rev. B 75, 224107, doi:10.1103/PhysRevB.75.224107 (2007). 	21.	 Norris, S. A. Stress-induced patterns in ion-irradiated silicon: model based on anisotropic plastic flow. Phys. Rev. B 86, 235405,
arxiv:1207.5754 (2012).

	22.	 Bradley, R. M. Exact linear dispersion relation for the sigmund model of ion sputtering. Physical Review B 84, 075413, doi:10.1103/ PhysRevB.84.075413 (2011).

	23.	 Mario, Castro & Rodolfo, Cuerno Hydrodynamic approach to surface pattern formation by ion beams. Applied Surface Science 258, 4171–4178, doi:10.1016/j.apsusc.2011.09.008 (2012).

	24.	 Hofsäss, H., Zhang, K., Pape, A., Bobes, O. & Brötzmann, M. The role of phase separation for self-organized surface pattern

formation by ion beam erosion and metal atom co-deposition. Appl. Phys. A 111, 653–664, doi:10.1007/s00339-012-7285-8 (2012).

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

11

www.nature.com/scientificreports/
	25.	 Nietiadi, M. L., Sandoval, L., Urbassek, H. M. & Möller, W. Sputtering of Si nanospheres. Physical Review B 90, 045417, doi:10.1103/ PhysRevB.90.045417 (2014).
	26.	 Scott, A. Norris. Pycraters: A python framework for the collection of crater function statistics. arXiv:1410.8489. 	27.	 Madi, C. S.. Linear Stability and Instability Patterns in Ion Bombarded Silicon Surfaces. PhD thesis (Harvard University, 2011). 	28.	 Eaglesham, D. J., White, A. E., Feldman, L. C., Moriya, N. & Jacobson, D. C. Equilibrium shape of Si. Phys. Rev. Lett. 70, 1643–1646,
doi:10.1103/PhysRevLett.70.1643 (1993). 	29.	 Ishii, Y., Madi, C., Aziz, M. J. & Chason, E. Stress evolution in si during low-energy ion bombardment. Journal of Materials Research
(2014). 	30.	 Norris, S. A., Brenner, M. P. & Aziz, M. J. From crater functions to partial differential equations: A new approach to ion
bombardment induced nonequilibrium pattern formation. J. Phys. Cond. Matt. 21, 224017, doi:10.1088/0953-8984/21/22/224017 (2009). 	31.	 Norris, S. A. et al. Molecular dynamics of single-particle impacts predicts phase diagrams for large scale pattern formation. Nature Communications 2, 276, doi:10.1038/ncomms1280 (2011). 	32.	 Harrison Matt, P. & Bradley, R. Mark Crater function approach to ion-induced nanoscale pattern formation: Craters for flat surfaces are insuffcient. Physical Review B 89, 245401, doi:10.1103/PhysRevB.89.245401 (2014). 	33.	 Bukonte, L. et al. Comparison of molecular dynamics and binary collision approximation simulations for atom displacement analysis. Nuclear Instruments and Methods in Physics Research B 297, 23–28, doi:10.1016/j.nimb.2012.12.014 (2013). 	34.	 Castro, M., Gago, R., Vázquez, L., Muñoz-García, J. & Cuerno, R. Stress-induced solid flow drives surface nanopatterning of silicon by ion-beam irradiation. Physical Review B 86, 214107, doi:10.1103/PhysRevB.86.214107 (2012). 	35.	 Moreno-Barrado, A. et al. Nonuniversality due to inhomogeneous stress in semiconductor surface nanopatterning by low-energy ion-beam irradiation. Physical Review B (2015). 	36.	 Bock, W., Gnaser, G. & Oechsner, H. Modification of crystalline semiconductor surfaces by low-energy Ar+ bombardment: Si(111) and Ge(100). Surface Science 282, 333–341, doi:10.1016/0039-6028(93)90938-G (1993). 	37.	 Cuerno, R. & Barabási, A.-L. Dynamic scaling of ion-sputtered surfaces. Phys. Rev. Lett. 74, 4746–4749, doi:10.1103/ PhysRevLett.74.4746 (1995). 	38.	 Liedke, B. Ion beam processing of surfaces and interfaces: Modeling and atomistic simulations. PhD thesis (Helmholtz Zentrum Dresden Rossendorf, 2011).
Acknowledgements
S.A.N. was supported by NSF DMR-1408642. J.P. and M.J.A. were supported NSF-DMR-1409700. M.M., E.A., and K.F.L. were supported by NSF DMR-1307979. We are grateful to John Snyder, Emily Ghosh, and Gözde Erdem Rainville for performing wafer curvature measurements. Use of the National Synchrotron Light Source, Brookhaven National Laboratory, was supported by the US Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-98CH10886. This work is based partially upon research conducted at the Cornell High Energy Synchrotron Source (CHESS) which is supported by the National Science Foundation and the National Institutes of Health/National Institute of General Medical Sciences under NSF award DMR-1332208.
Author Contributions
S.A.N. developed the PyGLIDRE library, analyzed the experimental data, performed BCA simulations, and wrote the final version of the manuscript. J.C.P. performed experiments at BNL and CHESS, wrote the experimental methods section, and contributed to strategies for smoothing GISAXS data. M.M. performed experiments at CHESS and oversaw wafer curvature measurements. E.A. performed experiments at BNL and wrote the initial draft of the manuscript. M.J.A. supervised J.C.P. and co-initiated the project. K.F.L. co-initiated the project, supervised M.M. and E.A., and provided overall guidance.
Additional Information
Supplementary information accompanies this paper at doi:10.1038/s41598-017-01059-x Competing Interests: The authors declare that they have no competing interests. Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
© The Author(s) 2017

Scientific Reports | 7: 2016 | DOI:10.1038/s41598-017-01059-x

12