Uniwersytet Jagiello«ski w Krakowie Wydział Fizyki, Astronomii i Informatyki Stosowanej Dawid Maci¡»ek Nr albumu: 1064310 Badanie efektw morfologicznych podczasbombardowaniapowierzchni atomami oraz klastrami o energiach keV-owych Praca doktorska na kierunku Fizyka Praca wykonanapod kierunkiem Prof. dr. hab. ZbigniewaPostawy Krakw 2020 O±wiadczenie autora pracy Jani»ejpodpisanyDawid Maci¡»ek(nr indeksu: 1064310) doktorantWydziału Fizyki, AstronomiiiInformatyki Stosowanej Uniwersytetu Jagiello«skiego o±wiadczam, »e przedło»ona przeze mnie rozprawa doktorska pt. „Badanie efektw morfologicznych podczas bombardowania powierzchni atomami oraz klastrami o energiach keV-owych” jest oryginalnai przedstawia wyniki bada« wykonanych przeze mnie osobi±cie,pod kierunkiem prof. dr. hab. ZbigniewaPostawy. Prac¦ napisałem samodzielnie. O±wiadczam, »e moja rozprawa doktorska została opracowana zgodnie z Ustaw¡ o prawie autorskimiprawachpokrewnych z dnia4lutego 1994r. (Dziennik Ustaw1994nr 24poz.83 wrazzpniejszymi zmianami). Jestem ±wiadom, »e niezgodno±¢ niniejszego o±wiadczenia zprawd¡ ujawnionawdowolnym czasie, niezale»nieod skutkw prawnych wynikaj¡cych z ww. ustawy, mo»espowodowa¢ uniewa»nienie stopnianabytegonapodstawietej rozprawy. ................................ .......................................... Krakw, dnia Podpis autora pracy O±wiadczenie kieruj¡cego prac¡ Potwierdzam, »e niniejsza praca została przygotowanapodmoim kierunkiemikwalifkuje si¦do przedstawieniajejwpost¦powaniuo nadanietytułu zawodowego. ................................ ................. ........................... Krakw, dnia Podpiskieruj¡cego prac¡ Podzi¦kowania Przede wszystkim chciałbympodzi¦kowa¢ mojemu promotorowi prof. dr hab. Zbigniewowi Postawie za zaanga»owanie oraz wsparcie we wszystkich moich naukowych d¡»eniach. Szczegnie doceniam udost¦pnione mi mo»liwo±ci rozwoju takie jak sta»e, wyjazdy oraz wspprace mi¦dzynarodowe. Chciałbymrnie»podzi¦kowa¢ kolegomz zespołu, Michałowi oraz Mikołajowi, za wspnie sp¦dzony czas na doktoracie a w szczegno±ci za inspiruj¡ce dyskusje naukowe. Pragn¦rnie»podzi¦kowa¢rodzinie za wsparcie na wszystkich szczeblach mojejedukacji, kte doprowadziło mnie do tego momentu. Spis tre±ci 1 Wst¦p 2 1.1Forma pracydoktorskiej ............................ 2 1.2 Zagadnieniebadawcze ............................. 3 1.3Kontekstnaukowy ............................... 4 2 Opis bada« składaj¡cych si¦ na prac¦ doktorsk¡ 11 2.1 Wpływ morfologiipowierzchni na procesbombardowania . . . . . . . . . . 11 2.1.1 Modelowanie profli gł¦boko±ciowych . . . . . . . . . . . . . . . . . 11 2.1.2 Modelowanie profli gł¦boko±ciowych rozszerzone o reakcje chemiczne 16 2.1.3 Efektychemicznepodczasbombardowania wywołane morfologi¡powierzchni ................................ 19 2.2 Wpływbombardowania na morfologi¦powierzchni . . . . . . . . . . . . . . 21 2.2.1 Wpływ implantowanych atomw gazw szlachetnychna morfologi¦ powierzchniwformalizmie funkcji krateru . . . . . . . . . . . . . . 25 2.2.2 Porwnanie metodydynamiki molekularnej oraz metodyMonte Carlo w ramachformalizmu krateru ..................... 28 2.2.3 Prosty model szorstko±ci dla metody Monte Carlo . . . . . . . . . . 32 2.2.4 Intuicyjny model opisuj¡cy zjawisko indukowania morfologii powierzchniwi¡zk¡jonow¡ ........................ 34 2.3Podsumowaniepracy .............................. 40 3 Przedruki artykułw 48 Rozdział1 Wst¦p 1.1 Forma pracy doktorskiej Niniejsza rozprawa została napisana jako podsumowanie wcze±niejszych działa« naukowych opublikowanych w nast¦puj¡cych artykułach [1–7]: 1. Dawid Maci¡»ek,RobertJ.Paruch,ZbigniewPostawa, BarbaraJ. Garrison„Microand Macroscopic Modeling of Sputter Depth Profling“ J. Phys. Chem. C, 2016, 120, 25473-25480 2. Sloan J. Lindsey, Gerhard Hobler, Dawid Maci¡»ek, Zbigniew Postawa „Simple model of surface roughness for binary collision sputtering simulations“ Nucl. In-strum. Methods Phys. Res.,2017, 393, 17-21 3. Hua Tian, Dawid Maci¡»ek, Zbigniew Postawa, Barbara J. Garrison, Nicholas Winograd „C-O Bond Dissociation and Induced Chemical Ionization Using High Energy (CO2 )n+ Gas Cluster Ion Beam“ J. Am. Soc. Mass Spectrom., 2018, 30, 476-481 4. Gerhard Hobler, Dawid Maci¡»ek, ZbigniewPostawa „Crater function moments: Role of implanted noble gas atoms“ Phys. Rev. B, 2018, 97, 155307 5. Gerhard Hobler, Dawid Maci¡»ek, ZbigniewPostawa „Ion bombardment induced atom redistribution in amorphous targets: MD versus BCA“ Nucl. Instrum. Methods Phys. Res., 2019, 447, 30-33 6. NunzioTuccitto, Dawid Maci¡»ek,ZbigniewPostawa,Antonino Licciardello„MD-Based Transport and Reaction Model for the Simulation of SIMS Depth Profles of MolecularTargets“. Phys. Chem. C, 2019, 123, 20188-20194 7. Dawid Maci¡»ek, Michał Ka«ski, Zbigniew Postawa „Intuitive Model of Surface Modifcation Induced by Cluster Ion Beams“ Anal. Chem., 2020, 92, 10, 7349–7353 Pracapodzielona jest na trzy cz¦±ci.W pierwszej cz¦±ci przedstawiampodejmowane przeze mnie zagadnienie badawcze wraz zkontekstem naukowym.Wdrugiej cz¦±cipokrce opisuj¦ artykuły składaj¡cesi¦nat¦ prac¦.W ostatniej cz¦±ci zamie±ciłem przedruki opisywanychpublikacji, kte umo»liwiaj¡ dogł¦bne zapoznaniesi¦z tre±ci¡ opisywanychprac. 1.2 Zagadnienie badawcze Celem niniejszej pracy doktorskiej jest zbadanie efektw morfologicznych (kształt powierzchni prki)pojawiaj¡cychsi¦podczasbombardowaniapowierzchnipociskami atomowymi oraz klastrowymi. Poniewa» przyczynowo±¢ tych efektw mo»e działa¢ w dwie strony, opis moichdziała« naukowychpodzieliłem na dwie sekcje tematyczne, kte zgł¦biaj¡ nast¦puj¡ce zagadnienia: • Jak morfologia prki wpływa na zachodz¡cy na jejpowierzchniproces bombardowania i emisji cz¡stek. • Jak bombardowanie prki wpływa na wytwarzan¡ na jejpowierzchni morfologi¦. Zrozumienie wy»ej przedstawionychzagadnie« jest bardzo istotne,poniewa»bombardowanie wi¡zk¡jonow¡jest fundamentaln¡ cz¦±ci¡wielu technik badawczych oraz procesw technologicznych[8]. Przykładem takiej techniki badawczej jest spektrometria masowa jonw wtnych (SIMS), kta jest najczulsz¡ metod¡ badania składuchemicznego powierzchni [9]. W tym przypadku wi¡zka jonw słu»y do emisji atomw lub molekuł prkiz wybranego obszaru najejpowierzchni.Wyemitowane na skutek uderzeniapocisku zjonizowane cz¡stki trafaj¡ nast¦pniedospektrometru masowego, gdzie mierzone jest widmo masowe. Schematyczne przedstawienie tej techniki znajduje si¦ na rysunku 1.1. Wstandardowym układzie mierzymytylko te cz¡stki, kte uległy jonizacji podczas uderzenia. Dost¦pna jestrwnie»odmiana tej techniki, w ktej wyemitowane cz¡stki oboj¦tne s¡ jonizowane zapomoc¡ wi¡zki laserowej [9].Kolejnym przykładem zastosowania wi¡zki jonw jest usuwanie materiału z badanejprki, aby odsłoni¢ jej wn¦trze. Je±lipoł¡czymytenprocesz technikami takimi jak spektroskopia fotoelektronw (XPS) [8], spektroskopia elektronw Augera (AES) [8] lub wspomniany SIMS, to uzyskamymo»liwo±¢ zbierania in-Rysunek 1.1: Schemat działania aparatury formacji o składzie chemicznym z całej ob-SIMS z detektorem czasu przelotu. j¦to±ci prki. Podej±cie takie nazywa si¦ proflowaniemgł¦boko±ciowym.Najlepszym przykłademwykorzystania wi¡zkijonowejw przemy±le jest proces domieszkowaniapprzewodnikw, kty jest jednymz krokwpodczas wytwarzania układw scalonych[8], kte s¡podstaw¡dzisiejszej cywilizacji. 1.3 Kontekst naukowy Wprocesiebombardowaniapowierzchnipociskami atomowymi lub klastrowymi dochodzi do szeregu zjawisk zachodz¡cych na powierzchni bombardowanej prki pod wpływem uderzeniapocisku.Do najwa»niejszych z nich nale»¡ rozpylanie, jonizacja oraz redystrybucja masy. Rozpylanie opisuje zjawisko erozji prki z jednoczesn¡emisj¡ cz¡stek stano-wi¡cychjej integraln¡ cz¦±¢ przed uderzeniem. Przyczyn¡ tego zjawiska jest deponowanie energii kinetycznejpociskuw rozpylanym materiale. Energiata mo»e by¢zdeponowana bezpo±rednio do układu j¡drowego na skutek elastycznych zderze«pomi¦dzy atomami pocisku oraz prki, jaki do struktury elektronowej układu [10]. Je±li ilo±¢ elastycznych zderze« wywołanych bombardowaniem jest na tyle mała, »e mo»emy zało»y¢ i» zderzenia nast¦puj¡ wył¡czniepomi¦dzy atomemb¦d¡cymwspoczynkua atomemobdarzonym energi¡ kinetyczn¡,tomwimywtedyo liniowejkaskadzie zderze«[10].Kiedyto zało»enie przestaje obowi¡zywa¢, a do zderze« dochodzi rwnie» mi¦dzy atomami w ruchu, to mwimywtedyo efektachnieliniowych[10].Wtrakcie propagacjikaskady zderze«, na skutek wzbudze« elektronowychcz¦±¢ atomw lub molekuł ulega jonizacji [10]. Je±li taka cz¡stka zostanie wyemitowania nim utraci uzyskany ładunek, to nazywamy j¡ jonem wtnym. Propaguj¡casi¦ wprcekaskada zderze«prowadzirwnie»dotrwałegoprzemieszczenia si¦ materiału w jej obr¦bie w procesie nazywanym redystrybucj¡ masy. Wa»nym czynnikiem dyktuj¡cym przebieg procesubombardowania jesttyp stosowanegopocisku.Pociskami mog¡by¢pojedyncze atomy,molekuły,klastry atomowe oraz klastry molekularne.Wybpocisku dla danego eksperymentu ma wpływ na wspczynnik rozpylania, prawdopodobie«stwo jonizacji oraz fragmentacj¦ molekuł bombardowanego podło»a [11].Wa»nym momentemw historii rozwoju dział jonowychgeneruj¡cychwspomnianepociskibyło wprowadzenie dział klastrowych[11]. Redukcja fragmentacji molekuł podło»a oraz mechanizm „samoczyszczenia“ (self-cleaning) wporwnaniudopociskw mono-atomowych miała bardzo du»y wpływ na jako±¢ uzyskiwanych profli gł¦boko±ciowych [11]. Morfologia bombardowanej powierzchni ma znacz¡cy wpływ na zachodz¡ce na niej procesy. Wielko±¢ amplitudy, jakii kształt struktur napowierzchni wpływaj¡ na wspczynnik rozpylania [2], rozkładk¡towy emitowanych cz¡steczek [12], prawdopodobie«stwo dysocjacji molekułpocisku [3] oraz jako±¢ uzyskiwanychprofli gł¦boko±ciowych[12]. Prka mo»eposiada¢ okre±lon¡ morfologi¦ przed rozpocz¦ciem procesubombardowania, jak rwnie» morfologia mo»eby¢ na jejpowierzchni indukowana na skutek tego»bombardowania. Bardzo ciekawym przykładem takiego zjawiska jestpowstawanie samoorganizuj¡cych si¦ struktur wykazuj¡cych uporz¡dkowanie dalekiego zasi¦gu [13]. Fascynuj¡ce w tymprocesiejestto,»ezperspektywypowierzchnimiejsce uderzeniapociskujestcałkowicie losowe, natomiastpowstaj¡ce na skutektychpojedynczychuderze« struktury cechuj¡ si¦uporz¡dkowaniemwskali znacznie wi¦kszej ni»skalapojedynczego uderzenia. Rne przykłady takiego zjawiska przedstawione s¡ na rysunku 1.2. Strukturami, kte mog¡ powsta¢ s¡: nano-zmarszczki [14–16] (nanoripples) nazywane rwnie» riplami, kte wy-gl¡dem przypominaj¡ makroskopowe falkipowstaj¡ce na piaskupodwpływem działania wiatru, nano-dziury (nanohole) [17], nano-kropki (nanodots) [14] oraz nano-pilary [18]. Poniewa»procesbombardowaniajest bardzo zło»ony,jegoopiszapomoc¡modelianalitycznychniedajepełnego obrazu sytuacji.Dopełnego zrozumieniazjawiskzachodz¡cych podk¡temθ=0◦ Ek=3keV [18] napowierzchni prki niezb¦dne jest zastosowanie metod numerycznych, takich jak symulacjekomputerowe. Przykładem metody, kta wykorzystywana jest do modelowania zjawiskabombardowania jest dynamika molekularna.W dynamice molekularnej modelowane atomys¡ reprezentowane jako punktymaterialne o okre±lonej masie, ktych ruch dyktowany jest drug¡ zasad¡ dynamiki Newtona [19]. Sprowadza si¦ to do tego, »e aby uzyska¢ informacj¦opoło»eniachtychatomww czasie, co nazywamytrajektori¡ układu, nale»y rozwi¡za¢ sprz¦»ony układklasycznych rwna« ruchu Newtona: d2r~i X ~ mi = Fij, i,j =1, ..., N, (1.1) dt2 i0“ wprowadza ograniczenie na taki przypadek.Trzecim sprawdzanym sposobemkonstrukcji trajektoriipo zderzeniu jest TRIM („noτ“),dla ktegopozycje atomw pozderzeniu nies¡korygowane napodstawie całki czasowejτ (rysunek 2.12b). Rysunek 2.12: Schematyczne przedstawienie trajektorii zderzaj¡cych si¦ atomwpocisku A (kolor zielony) oraz tarczy B (kolor czerwony) obliczonymi na podstawie całki czasowej. Czarn¡ strzałk¡ oznaczono drog¦ swobodnego przelotu. A0 oznacza pozycj¦ atomu pocisku przed zderzeniem, natomiast pozycja atomu tarczy oraz atomu pocisku po zderzeniu,a przedkontynuowaniem dalszej drogi toa) A0 1 , B1 0 po uwzgl¦dnieniupoprawki wynikaj¡cej z przeci¦cia asymptot wchodz¡cej oraz wychodz¡cej, b) A1, B1 bez uwzgl¦dnienia tejpoprawki. Metodologia Aby jak najlepiej odwzorowa¢ rzeczywist¡ struktur¦powstał¡ wewn¡trzbombardowanej prki krzemu, w symulacjachmetod¡ dynamiki molekularnej wykorzystałem prk¦ krzemu uzyskana na skutek bombardowania wi¡zk¡ 2 keV Kr θ = 60◦ w poprzedniej publikacji. W badaniach interesował mnie wpływ pierwotnej energii kinetycznej nadawanej atomowi prki na pierwszy moment funkcji krateru. Badane energie kinetyczne znajdowały si¦ w zakresie 5 eV do 100 eV. Dla ka»dej badanej energii kinetycznej przeprowadziłem 300 symulacji, w ktych atom, ktemu przypisywałem pierwotn¡ energi¦ kinetyczn¡ losowany był z regionu 3-6 nm pod powierzchni¡ prki. Ka»da symulacja wykonywanabyła nanowejkopii badanej prki, dzi¦ki czemu wszystkie przeprowadzone symulacje były wykonane na identycznym układzie. U±rednione wyniki symulacji uzyskane metod¡ dynamiki molekularnej dla rnych energii pierwotnych przedstawione s¡ na rysunku 2.13a, wraz z wynikami metody Monte Carlo dla trzech badanych metod poprawek pozycji po zderzeniu. W symulacjach Monte Carlo, parametr Ed został tak dobrany, abyuzyska¢ najlepsze dopasowanie do wynikw symulacji metod¡ dynamiki molekularnej. Jak wida¢, symulacje wykorzystuj¡ce metod¦„˙p>0“ oraz TRIM dokorekcji pozycji atomwpo zderzeniu, uzyskały bardzo dobre dopasowaniedo wynikw uzyskanych metod¡ dynamiki molekularnej. Rysunek 2.13: a)Porwnanie wynikw trzechmetodwprowadzeniapoprawkidopozycji atomwpo zderzeniu, domy±lnyspos modelu IMSIL, „˙p>0“ oraz TRIM.Warto±ci Ed zostały tak dobrane,abyuzyska¢ najlepsze dopasowaniedo wynikw dynamiki molekularnej.b)Zale»no±¢ pierwszego momentu M(1) odpocz¡tkowejgł¦boko±ci pierwotnego atomu o energii 100 eV dla metody „˙p>0“ dla trzech zwrotw wektora pr¦dko±cipocz¡tkowej wzgl¦dem kierunku rwnoległego dopowierzchni. Nast¦pnie zbadałem wpływ gł¦boko±ci na jakiejpoło»onyjest atom zpocz¡tkow¡ energi¡ nawarto±¢ pierwszego momentu.Wyniki dla energii 100eV oraz trzech rnych orientacjiwektorapr¦dko±ci przedstawia rysunek 2.9b. Je±li o± x skierowana jest rwnolegle dopowierzchni, to wybranymi kierunkamibyły: rwnoległy do x,podk¡tem 60◦ do x w kierunkupowierzchni orazpodk¡tem60◦ odpowierzchniw gł¡b prki.Poniewa» pierwszy moment obliczany jest w kierunku x, to aby mo»nabyłoporwna¢ otrzymane wyniki, na rysunku 2.9b wyniki dla kierunkww stron¦powierzchni orazodpowierzchni pomno»one s¡ razy dwa(cos(60◦)=0.5). Ka»dy punkt na rysunku 2.9b został obliczony napodstawie 300 symulacji. Jak wida¢ na wykresie, dla gł¦boko±cipowy»ej4 nmwarto±ci pierwszego momentu nie zale»¡ od kierunku. Natomiast im bli»ej powierzchni znajduje si¦ atom, ktemu nadajemypocz¡tkow¡pr¦dko±¢,tym wi¦ksza jest rnicaw otrzymanych momentach w zale»no±ci od kierunku pr¦dko±ci. Spowodowane jest to tym, »e przy powierzchnikolektywnyruch atomw mo»e propagowa¢si¦wjej kierunku zmieniaj¡cjej kształt (efektynieliniowe). Dlakierunkuw stron¦powierzchni, dlapomiaru najbli»ejpowierzchni(0.5nm)warto±¢ pierwszego momentunaglespada,cospowodowanejesttym,»e du»a cz¦±¢ atomwzostajepo prostu rozpylonainie ma wkładudo procesu redystrybucji. Dodatkowo, dlaporwnania, na wykres zostały naniesione wyniki metody Monte Carlo o parametrach „˙p>0“, Ed=5eV. Dla tej metody, dla gł¦boko±cipowy»ej2.5 nm, obserwujemy stał¡warto±¢ pierwszego momentu, natomiast dla gł¦boko±ci bli»ejpowierzchni nast¦puje spadekwarto±ci pierwszego momentu.Ponowniejesttospowodowanetym,»e metoda Monte Carlo nie jest w stanie symulowa¢ efektw nieliniowych (odkształcanie powierzchni),a spadek tejwarto±cipodyktowanyjesttym, »e cz¦±¢ atomw zostaje rozpylonai nie ma wkładudo momentu. Dlagł¦boko±cipowy»ej4 nmwarto±ci otrzymane obiema metodami zbiegaj¡ si¦. Wnioski Otrzymane wynikipokazały,»eodpowiednio sparametryzowana metoda Monte Carlo, wykorzystuj¡ca model IMSIL mo»e uzyska¢ wyniki zgodne z symulacjami metod¡ dynamiki molekularnej. Dla badanego zakresu niskichenergii(Ek <100 eV), aby model IMSIL generowałpoprawne wyniki,koniecznebyło wprowadzenie modyfkacji do metody obliczania pozycji atomwpo zderzeniu. 2.2.3 Prosty model szorstko±ci dla metody Monte Carlo Sekcja oparta na publikacji: „Simple modelof surfaceroughness for binarycollision sputtering simulations“ Motywacja Jak przedstawiłem to we wcze±niejszych publikacjach, metoda Monte Carlo jest wa»nym narz¦dziemwykorzystywanymdo bada«procesubombardowaniaprkipociskamiatomowymi. Jednak dla du»ychk¡tw padania(θ>80◦)morfologiabombardowanejpowierzchni zaczyna mie¢ dominuj¡cy wpływ na dokładno±¢ uzyskiwanych t¡ metod¡ wynikw [49]. Bombardowaniapowierzchnipod takimik¡tami jest bardzo istotne w wytarzaniu nanostruktur zapomoc¡ techniki FIB [50]. Celem opisywanej tutajpublikacjibyło zbadanie czterech modeli szorstko±ci zaimplementowanych napotrzeby symulacji metod¡ Monte Carlo.Trzyz nich, nazywane modelami geometrycznymi, wprowadzaj¡ szorstko±¢bezpo±redniopoprzez zmian¦ kształtupowierzchni modelowanej prki: sinus,podwjnysinus oraz kształt trjk¡tny. Ka»dy z tych modeli zdefniowany był poprzez dwa parametry, amplitud¦ oraz długo±¢ fali. Zapomoc¡ dynamiki molekularnej wyznaczyłem rzeczywist¡ długo±¢ fali struktur, ktepowstaj¡ napowierzchni krzemupodczasbombardowaniapod k¡tami ±lizgowymi. Czwartymodel szorstko±ci modelowałchropowato±¢tylkopoprzez szeroko±¢warstwy przej±ciowejpomi¦dzy prni¡a prk¡.Wtejwarstwie g¦sto±¢ materiału zmieniałasi¦odzera, napowierzchni,do g¦sto±cimateriału wła±ciwegow gł¦biprbki. Metodologia Rysunek 2.14: a) Powierzchnia krzemu symulowana metod¡ dynamiki molekularnej, utworzona na skutek 2100 uderze«pociskiem5keV Gapodk¡tem θ=89◦.Kolor atomu reprezentuje wysoko±¢ atomu Si w kierunku normalnym do powierzchni. b) Porwnanie wynikw eksperymentalnych z metod¡ Monte Carlo wykorzystuj¡c¡ rne modele szorstko±ci powierzchni: i) wspczynnik rozpylania w zakresie 0◦-90◦, ii) wspczynnik rozpylania w zakresie 80◦-90◦, iii) wspczynnikodbiciapocisku. Dla tak du»ychk¡tw padania, krytycznym parametrem, jaki nale»y zapewni¢ w symulacjach metod¡ dynamiki molekularnej, jestodpowiednipotencjałoddziaływaniapomi¦dzypociskiema atomamipodło»a. Jakopotencjałoddziaływania galuz krzemem wybrałempotencjał Lennarda-Jonesa[25] ogł¦boko±ci studnipotencjału =0.9 eV.Warto±¢ ta została wybrana napodstawie szeregu niezale»nych symulacji wielokrotnegobombardowania krzemu galemo energii5 eVpodk¡tem89◦.Wsymulacjachzmieniałem warto±¢ ,a»do momentu,wktym obliczony napodstawiesymulacjiwspczynnikrozpylania, zgadzał si¦zwarto±ci¡ eksperymentaln¡(Yexp=1.7 [2]).Tak znalezion¡ warto±¢ parametru  zweryfkowałem symulacjami wykonanymidlak¡ta padania 85◦. Otrzymana z tych symulacji warto±¢ wspczynnika rozpylania Ysim=6.49 pokrywała si¦ niemal idealnie z warto±ci¡ eksperymentaln¡ Yexp=6.66 [2]. Potwierdza to, »e warto±¢  jest wybranapoprawnie. Napodstawie stanupowierzchnipobombardowaniupodk¡tem85◦ (Rysunek 2.14a), obliczyłem długo±¢ faliodpowiadaj¡c¡ strukturom napowierzchni, kta wyniosłaokoło λ=3nm.Warto±¢ta zostaławykorzystanawmodelachgeometrycznychw metodzieMonte Carlo. Optymalna amplituda chropowato±ci dla modeli geometrycznych oraz w modelu gradientowym została wyznaczona napodstawie wielu symulacji z rnymi warto±ciami tej amplitudy. Rysunek 2.14b przedstawia uzyskane wyniki wspczynnika rozpylenia dla badanych modeli z zadan¡ optymaln¡ amplitud¡ chropowato±ci. Dla du»ychk¡twpadania(rysunek 2.14b.ii)symulacje Monte Carlo wykorzystuj¡ce modele opisuj¡ce chropowato±¢ prki uzyskuj¡ znacznie dokładniejsze warto±ci wspczynnika rozpylania ni» symulacje Monte Carlo bazuj¡ce napodstawowym modelupłaskiej powierzchni. Jest to spowodowane tym, »e w modelach implementuj¡cych chropowato±¢ powierzchni, atompocisku ma mniejsz¡ szans¦ naodbiciesi¦odprkibezspowodowania rozpylenia(rysunek 2.14b.iii). Wnioski Zastosowanie modelichropowato±cipowierzchniw symulacjachmetod¡ Monte Carlo przeło»yło si¦ na dokładniejsze obliczenia wspczynnika rozpylania. Zwi¦kszenie dokładno±ci otrzymywanych wynikwbyło szczegnie widoczne dla du»ychk¡tw padania(θ>80◦), dla ktychobliczeniabezchropowato±cidawały bardzo niedokładne wyniki. Opracowany w tej publikacji model gradientowy dzi¦ki swojejprostocie (jeden parametr) oraz wydaj-no±ci obliczeniowej mo»eby¢ efektywniewykorzystywany wsymulacjachMonte Carlo. 2.2.4 Intuicyjny model opisuj¡cy zjawisko indukowania morfologiipowierzchni wi¡zk¡ jonow¡ Sekcja opartanapublikacji: „An IntuitiveModelof SurfaceModifcation Inducedby Cluster Ion Beams“ Motywacja Istnieje wiele modeli opisuj¡cychzjawiskopowstawania riplipodczasbombardowaniapowierzchni wi¡zk¡ jonw [13]. Modele te jednak maj¡ bardzo zło»on¡ defnicj¦ matematyczn¡, co nie jestkorzystn¡ cech¡,poniewa» im wi¦kszypoziom abstrakcji matematycznej, tym ci¦»ej jest zrozumie¢ wyst¦puj¡ce zjawiska. Dobrze demonstruje to wykorzystywany w wy»ej opisywanych publikacjach formalizmu krateru, gdziepowstawanie ripli tłumaczone jest napodstawiewarto±ci parametru C, natomiastparametr ten jestw zło»ony spos obliczany na podstawie momentw funkcji krateru. Je±li czytelnik nie ma wystarczaj¡co obszernej wiedzy na tematy matematyki wykorzystywanej przez model, to uzyskiwane napodstawie tegopodej±cia wnioski mog¡ wydawa¢ si¦ nie doko«ca uzasadnione. Taki stan rzeczy zmotywował mnie do zaproponowania modelu, kty bazowałby nabezpo±redniej fzycznej interpretacji, daj¡c łatwy do zrozumienia obraz zjawiskapowstawania ripli. Dodatkowym wymogiem dla tworzonego przeze mnie modelu było to, bypozwalał przewidywa¢, przy jakichparametrachwi¡zki (energiapocisku,k¡t padania) mo»emyspodziewa¢ si¦powstania ripli. Metodologia My±l¡ przewodni¡ proponowanego przeze mnie modelu jest zwrcenie uwagi na to jak lokalnyk¡t padania wi¡zki,a co zatym idzie lokalna redystrybucja masy wywołana uderzeniempocisku, przekładasi¦wskali całejpowierzchniprki na zwi¦kszanie lub zmniejszanie morfologii indukowanejbombardowaniem. Przy czym opisywany tutajlokalnyk¡t padania to taki, kty wst¦puje w miejscu padania pocisku, a jego warto±¢ zale»y od globalnegok¡ta padaniapocisku oraz od lokalnego nachyleniapowierzchni.Zkolei globalnymk¡tem padania jestk¡t eksperymentalny, czyli taki jaki padaj¡ca wi¡zkatworzy zpowierzchni¡ prki w skali makroskopowej. Rysunek 2.15 przedstawia wizualizacj¦ przykładowych globalnych k¡tw padania ΩA, ΩB (przerywana linia na wykresie M(θ) oraz towarzysz¡cych im lokalnychk¡tw padania(θi, θj )oraz(θk, θp)). Przedmiotem analizyjest wpływ lokalnej redystrybucji masy na morfologi¦powierzchni, dladwch przypadkw globalnegok¡tapadania,korzystaj¡czuproszczonego schematu przedstawionego na rysunku 2.15. W pierwszym przypadku globalny k¡t padania ΩA znajduje si¦ w regionie, w ktym warto±¢ funkcji redystrybucji masy M(θ) ro±nie wraz z lokalnym k¡tem padania. Oznacza to, »e redystrybucja masy dla k¡ta padania θi jest Rysunek 2.15: Wizualizacja wpływu kształtu k¡towej zale»no±ci redystrybucji masy (M(θ))wywołanej uderzeniempocisku na ewolucj¦ morfologiipowierzchni.W przedstawionych rozwa»aniachpowierzchni¦ traktuj¦ jako seri¦dyskretnych w¦złw, mi¦dzyktrymidochodzido wymianymasy. Wielko±¢ transferu masy uzale»niona jestodlokalnego k¡ta padania.Na ilustracjiwarto±¢ta zobrazowanajest wielko±ci¡ ¢wiartkikoła. Czerwony kolor reprezentuje wypływ masy z w¦zła, a zielonyprzypływ masy do w¦zła. mniejsza ni» dlak¡ta padania θj, je»eli θi >θj. ‚ledz¡c lokalne przemieszczanie si¦ masy, obserwujemyerozj¦ masy na wzniesieniach oraz akumulacj¦ masy w depresjach. Przeprowadzaj¡c takie samo rozumowanie dla przypadku, gdy globalnyk¡t padania ΩB znajduje si¦ w regionie, gdzie M(θ) maleje wraz ze wzrostem θ, otrzymamy odwrotnyefekt. Oznacza to, »e b¦dziemy obserwowa¢ akumulacj¦ masy na szczytach, natomiast erozj¦ masy w depresjach. Ł¡cz¡c oba efekty, mo»emy wnioskowa¢, »e wygładzanie lub zwi¦kszanie szorstko±cipowierzchni zale»yodk¡ta krytycznego θC , dla ktegofunkcja redystrybucji masy mawarto±¢ maksymaln¡. Je±lik¡t padania wi¡zki w eksperymencie znajdujesi¦ pomi¦dzy 0◦ a θC , to bombardowanie b¦dzie miało wygładzaj¡cy efekt. Natomiast je±li b¦dziepowy»ejθC , tob¦d¡powstawały riple. Abysprawdzi¢ przewidywania mojego modelu wybrałem dwa eksperymentyjako punkty odniesienia.Wtychdo±wiadczeniachkrzem [15] oraz złoto [16]byłybombardowane wi¡zk¡ 30keVAr3000.S¡ toprzykłady idealne,poniewa» autorzy przeprowadziliprocesbombardowania dla rnychk¡tw,apo zako«czonym eksperymencie, dlaka»dego przypadku, zmierzyli obrazy AFM, na ktychpodstawie mo»na wnioskowa¢ czypojawiły si¦ riple czy te» nie. Zapomoc¡ dynamiki molekularnej technik¡ wielokrotnych uderze« obliczyłem funkcj¦ redystrybucji masy dla rnychk¡tw padania(rysunek 2.16). Dla obu przypadkw maksimum znajdowało si¦ dokładnie w miejscu, dla ktego w układacheksperymentalnych zachodziło przej±cie z gładkiej powierzchni do pojawienia si¦ ripli, co jednoznaczniepotwierdzało skuteczno±¢ mojego modelu. Dodatkowo sprawdziłem, jak funkcja redystrybucji masy zale»e¢ b¦dzie od energii kinetycznejpocisku(rysunek 2.17).W tym przypadku badan¡ prk¡był krzembombardowanypociskiem Ar3000 o energiach10keV, 15keV, 20keV, 25keVi30keV. Ka»da energia to osobna seria symulacji wielokrotnego bombardowania dla rnych k¡tw padania.Poło»enie maksimum dla najmniejszej energiipocisku Ek=5 keV znajdowało si¦ w θc = 23◦, w miar¦ zwi¦kszania tej energiipoło»enie maksimum przesuwało si¦ w kierunku wy»szychk¡tw θc, osi¡gaj¡c najwi¦ksz¡warto±¢ θc = 45◦ dla najwy»szej badanej energii Ek=30keV. Oznacza to, »e wraz ze wzrostem energii kinetycznejpocisku,poło»enie maksimum przesuwasi¦w kierunku wi¦kszychk¡tw.Taki trendbył obserwowany wpracy[51].W tym eksperymencie prk¦organiczn¡bombardowano wi¡zk¡klastrow¡ Ar5000 o rnych energiach (5 keV, 10 keV, 20 keV). Dla energii 5 keV zaobserwowano, »ejako±¢ proflu (rozdzielczo±¢gł¦boko±ciowa mierzonejwarstwy delta)pogarszasi¦ wraz z gł¦boko±ci¡, a degradacja ta była spowodowana powstaniem morfologii powierzchni (riple), natomiast dla energii 10keV oraz 20keVjako±¢ proflupozostaje stała (efekt wy-gładzaj¡cy).Poniewa»k¡t padania dla wszystkich energiibył taki sam θ = 45◦, oznacza to, »e dla energii wi¡zki pomi¦dzy 5 keV a 10 keV powinna nast¡pi¢ zmiana poło»enia maksimum na funkcji redystrybucji z θc < 45◦ do θc > 45◦ . Wnioski W opisywanej pracypokazałem, »ebombardowaniepowoduj¡ce redystrybucj¦ masy zawsze powoduje wzmacnianie lub wygaszanie chropowato±ci powierzchni. To, w ktym z tych dwch re»imw znajduje si¦ bombardowana prka, zale»y od dwch głwnych czynnikw:pod jakim globalnymk¡tem wi¡zka pada napowierzchni¦i jaka jestk¡towa zale»no±¢ redystrybucji masy dla badanego układu. Redystrybucja masy musi mie¢ warto±¢ zero dla skrajnych k¡tw padania θ =0◦ (całkowita symetria uderzenia, czyli tyle samo masy jest przemieszcznew prawoi w lewo) oraz θ = 90◦ (pociskporusza si¦ prostopadle dopowierzchni, wi¦c nie mo»epowodowa¢ »adnej redystrybucji masy). Maksimum tej funkcjimusi wi¦c le»e¢ gdzie±pomi¦dzy tymiwarto±ciami.W prezentowanych badaniach warto±cik¡ta krytycznego zawierały si¦w przedziale 23◦ − 50◦.Fakt ten ma du»e znaczenie dla technik, w ktrych bombardowanie wykorzystywane jest do warstwowego usuwania prki. W tych przypadkach zale»y nam na tym, aby by¢ w regionie wygładzaj¡cym(θ<θC ).Na przykładw układachSIMS, wi¦kszo±¢komercyjnej aparaturyma na stałe zamontowane działo jonowe pod k¡tem θ = 45◦ do normalnej do powierzchni badanych układw. Co oznacza, »e dla cz¦±ci prek b¦dziemy indukowa¢ niepo»¡dan¡ morfologi¦powierzchni. Dzi¦ki zrozumieniu zjawiska dostarczanemu przezniniejszy model widzimy, »e zmianak¡ta padania wi¡zki naweto kilka stopniw kierunku mniejszych warto±ci skutkowałobyprzej±ciemdo re»imu wygładzaj¡cego.A cozatym idzie, znaczn¡ popraw¡ otrzymywanych profli gł¦boko±ciowych. Dodatkowo, zwi¦kszenie energii padaj¡cej wi¡zki przesuwak¡tkrytyczny w kierunku mniejszych warto±ci. Oznacza to, »e dla przypadkw kiedy mamyproblemzchropowato±ci¡,a nie mo»emyzmieni¢k¡ta padania, to jednymz mo»liwych rozwi¡za« jest zwi¦kszenie energii kinetycznej padaj¡cej wi¡zki. 2.3 Podsumowanie pracy Celem niniejszej pracy było zbadanie efektw morfologicznych zachodz¡cych podczas bombardowaniapowierzchnipociskami atomowymi oraz klastrowymi.W pierwszej cz¦±ci przedstawiłem prace, w ktych badania koncentrowały si¦ na wpływie morfologii powierzchni naprocesbombardowaniaiemisji cz¡stek.Wpublikacji„Micro-and Macroscopic Modeling of Sputter Depth Profling“ wprowadziłem model słu»¡cy do symulowania profli gł¦boko±ciowychuzyskiwanych w technice SIMS napodstawie danycheksperymentalnych oraz informacji uzyskanych z symulacji komputerowych. Stworzone za pomoc¡ modelu profle gł¦boko±ciowe dobrze odtwarzały wyniki eksperymentalne. Analizuj¡c wyniki generowane przez model, mo»na zaobserwowa¢, »e na kształt otrzymywanych profli dominuj¡cy wpływ ma chropowato±¢ prki, natomiast wpływ mieszania jonowego ma wpływ drugorz¦dny.Wkolejnej pracy „MD-Based Transport and Reaction Model for the Simulation of SIMS Depth Profles of MolecularTargets“ , wspomniany model rozwijamo efekty chemiczne zachodz¡cepodczasbombardowania. Badanym układembyłpolistyren bombardowanyklastrami argonowymi orazC60. Układ ten przechodzi proces sieciowania, kty mo»na redukowa¢, doprowadzaj¡c cz¡steczki NO nabombardowan¡powierzchni¦. Model dobrzeodwzorowałjako±ciowywpływ ilo±cipodawanychcz¡steczekNOna kształt profli gł¦boko±ciowych, jednak najwa»niejsz¡ wnioskiem z tej pracy było uzyskanie ilo±ciowej informacji, »ewydajno±¢ cz¡steczekNOw reakcjizwolnymirodnikamipowstaj¡cymi napowierzchnipodczasbombardowaniapociskiemC60 wynosi 10%.W ostatniej publikacjiporuszanej w tej cz¦±ci pracy „C-O Bond Dissociation and Induced Chemical Ionization Using High Energy (CO2 )n+ Gas Cluster Ion Beam“ za pomoc¡ symulacji metod¡ dynamiki molekularnejpokazałem, jak morfologiapowierzchni wpływa na proces fragmentacji molekuł CO2 podczasbombardowania klastrami (CO2)n. Druga cz¦±¢ niniejszej pracykoncentrowałasi¦ na wpływiebombardowania napowstaj¡c¡ morfologi¦powierzchni.Wtym celuw cz¦±ci pracwykorzystywałem formalizm funkcji krateru, kty przewiduje efekt wygładzania lub zwi¦kszenia chropowato±ci w oparciu o momentyfunkcji krateru obliczane napodstawie symulacji metod¡dynamiki molekularnej lub metod¡ Monte Carlo. W pierwsze pracy „Crater function moments: Role of implanted noble gas atoms“ pokazałem, »e implantowane atomy pocisku maj¡ istotny wkład do pierwszego momentu funkcji krateru. Jest to istotny wniosek,poniewa» wszystkie dotychczasowe badaniaignorowałytenwpływ.Dodatkowoporwnaniewynikwzsymulacji metod¡dynamiki molekularnejzmetod¡MonteCarlo,pokazało,»etadruganie doszacowuje wkładu do pierwszego momentu odredystrybucji atomw prki. Problem ten został zbadanyw publikacji „Ionbombardment inducedatomredistributionin amorphous targets: MD versus BCA“, w ktej badałem wpływ parametrw metody Monte Carlo na dokładno±¢ otrzymywanych za jejpomoc¡ pierwszych momentw funkcji krateru.Warto±ciami odniesienia dlatychsymulacjibyływarto±ci uzyskane zapomoc¡ dynamiki molekularnej. Wynikipokazały, »e mo»liwe jest takie dobranie parametrw metody Monte Carlo, aby rozwi¡za¢ problem niedoszacowywania wkładuodredystrybucji atomw prki.Wkolejnej opisywanej publikacji „Simple model of surface roughness for binary collision sputtering simulations“ badałem efekt wprowadzenia rnych modeli chropowato±ci w metodzie Monte Carlona dokładno±¢wynikwuzyskiwanychzapomoc¡tej metody, ze szczegnym naciskiem dla ±lizgowychk¡tw padania wi¡zki. Zastosowanie modelichropowato±ci, kte były parametryzowane na podstawie symulacji metod¡ dynamiki molekularnej, skutkowało znacznym zwi¦kszeniem dokładno±ci otrzymywanych wynikw, w szczegno±ci dla k¡tw ±lizgowych.W ostatniej opisywanej publikacji „An Intuitive Model of Surface Modifcation Induced by Cluster Ion Beams“ stworzyłem nowy model, kty w przyst¦pny spos tłumaczy zjawisko powstawania lub wygaszania morfologii powierzchni podczas bombardowania wi¡zk¡klastrow¡. Przewidywania modelu opieraj¡ si¦ napoło»eniuk¡ta krytycznegowk¡towej zale»no±ci redystrybucji masy wywołanej uderzeniem.Wpublikacjipokazałem, »eaby bombarduj¡ca wi¡zka miaławygładzaj¡cy efekt napowierzchni¦, k¡t padania musi by¢ mniejszy ni» k¡t krytyczny. Fakt ten ma bardzo du»e znaczenie w technikachproflowania gł¦boko±ciowego.Dodatkowopokazałem, »e zwi¦kszenie energii kinetycznejpocisku przesuwapoło»eniek¡takrytycznegow kierunku mniejszychwarto±ci. Bibliografa [1] Dawid Maciazek, Robert J. Paruch, Zbigniew Postawa, and Barbara J. Garrison. Micro-and macroscopic modeling of sputter depth profling. The Journal ofPhysical ChemistryC, 120(44):25473–25480, October 2016. [2] Sloan J. Lindsey, Gerhard Hobler, Dawid Maci¡»ek, and Zbigniew Postawa. Sim-ple model of surface roughness for binary collision sputtering simulations. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials andAtoms, 393:17–21,February 2017. [3] Hua Tian, Dawid Maci¡»ek, Zbigniew Postawa, Barbara J. Garrison, and Nicholas Winograd. C-obond dissociation and inducedchemical ionization using high energy (co2)n gas cluster ionbeam. Journal of TheAmericanSociety for MassSpectrometry, 30(3):476–481, November 2018. [4] Gerhard Hobler, Dawid Maci¡»ek, andZbigniewPostawa. Crater function moments: Role of implanted noble gas atoms. PhysicalReviewB, 97(15),April 2018. [5] G. Hobler, D. Maci¡»ek, and Z.Postawa. Ionbombardmentinduced atom redistribution in amorphous targets: MD versus BCA. Nuclear Instruments and Methods in PhysicsResearchSectionB:Beam Interactions with MaterialsandAtoms,447:30–33, May2019. [6] Nunzio Tuccitto, Dawid Maciazek, Zbigniew Postawa, and Antonino Licciardello. MD-based transport and reaction model for the simulation of SIMS depth profles of molecular targets. The Journal of Physical ChemistryC, 123(33):20188–20194, July 2019. [7] Dawid Maci¡»ek, Micha Ka«ski, and Zbigniew Postawa. Intuitive model of surface modifcation inducedbycluster ionbeams. Analytical Chemistry, 92(10):7349–7353, April 2020. [8]FaridF. Umarov andAbdiravufA. Dzhurakhalov. Ionbombardment-induced surface e˙ectsin materials.InWaldemarA.Monteiro, editor, Radiation E˙ects in Materials, chapter 14. IntechOpen, Rijeka, 2016. [9] J. C. Vickerman. ToF-SIMS: surface analysisby massspectrometry. IM, Chichester, 2001. [10]Peter Sigmund. Sputteringby ionbombardment theoretical concepts. In Topics in Applied Physics, pages 9–71. Springer Berlin Heidelberg, 1981. [11] Nicholas Winograd. Gas cluster ionbeams for secondary ion mass spectrometry. AnnualReview ofAnalytical Chemistry, 11(1):29–48, June 2018. [12] Dawid Maciazek, Michal Kanski, Lukasz Gaza, Barbara J. Garrison, and Zbigniew Postawa. Computer modeling of angular emission from ag(100) and mo(100) surfaces due to arn cluster bombardment. Journal of Vacuum Science & Technology B, Nanotechnology and Microelectronics: Materials, Processing, Measurement, and Phenomena, 34(3):03H114, May2016. [13] Scott A. Norris and Michael J. Aziz. Ion-induced nanopatterning of silicon:Toward a predictive model. Applied PhysicsReviews, 6(1):011311, March2019. [14] B. Ziberi, F.Frost, and B. Rauschenbach.Formation of large-area nanostructures on si and ge surfaces during low energy ionbeam erosion. Journal of Vacuum Science &Technology A:Vacuum, Surfaces, and Films, 24(4):1344–1348, July 2006. [15] Omar Lozano,Q.Y. Chen,B.P. Tilakaratne,H.W. Seo,X.M.Wang,P.V.Wadekar, P. V. Chinta, L. W. Tu, N. J. Ho, D. Wijesundera, and W. K. Chu. Evolution of nanoripples on silicon by gas cluster-ion irradiation. AIP Advances, 3(6):062107, June 2013. [16] BuddhiTilakaratne,QuarkChen,andWei-KanChu. Self-assembledgold nano-ripple formationbygas cluster ionbeambombardment. Materials, 10(9):1056, September 2017. [17] MonikaFritzsche, Arndt Muecklich, andStefanFacsko. Nanohole pattern formation on germanium inducedbyfocusedionbeamand broadbeamGa+ irradiation. Applied PhysicsLetters, 100(22):223108, May2012. [18] Elzbieta Trynkiewicz, Benedykt R. Jany, Dominik Wrana, and Franciszek Krok. Thermally controlled growth of surface nanostructures on ion-modifed AIII-BV semiconductor crystals. Applied Surface Science, 427:349–356, January 2018. [19] M. P. Allen. Computer simulation of liquids. Clarendon Press Oxford University Press, Oxford England NewYork, 1987. [20] Steve Plimpton.Fast parallel algorithms for short-range molecular dynamics. Journal of Computational Physics, 117(1):1–19, March 1995. [21] William C. Swope, Hans C. Andersen, Peter H. Berens, and Kent R. Wilson. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. The Journal of Chemical Physics, 76(1):637–649, January 1982. [22] Barbara J. Garrison and Zbigniew Postawa. Computational view of surface based organic mass spectrometry. Mass SpectrometryReviews, 27(4):289–315, July 2008. [23] Michael F. Russo, Zbigniew Postawa, and Barbara J. Garrison. A computational investigation of c60 depth profling of ag: Molecular dynamics of multiple impact events. The Journal of Physical ChemistryC, 113(8):3270–3276, January 2009. [24] James F. Ziegler and JochenP. Biersack. The stopping and range of ions in matter. In Treatise on Heavy-Ion Science, pages 93–129. Springer US, 1985. [25] On the determination of molecular felds. —II. from the equation of state of a gas. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 106(738):463–477, October 1924. [26] Murray S. Daw and M. I. Baskes. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals. Physical Review B, 29(12):6443–6453, June 1984. [27] J.Terso˙. New empirical approachfor the structure and energy of covalent systems. PhysicalReviewB, 37(12):6991–7000, April 1988. [28]AdriC.T.van Duin, Siddharth Dasgupta,Francois Lorant,and WilliamA.Goddard. ReaxFF: a reactive force feld forhydrocarbons. The Journal of Physical Chemistry A, 105(41):9396–9409, October 2001. [29] Zbigniew Postawa, Bartlomiej Czerwinski, Marek Szewczyk, Edward J. Smiley, Nicholas Winograd, and Barbara J. Garrison. Microscopic insights into the sputtering of ag{111} inducedbyc60and gabombardment. The Journal of Physical Chemistry B, 108(23):7831–7838, June 2004. [30] S. Sun, C. Szakal, T. Roll, P. Mazarov, A. Wucher, and N. Winograd. Use of c60 cluster projectiles for sputter depth profling ofpolycrystalline metals. Surface and InterfaceAnalysis, 36(10):1367–1372, 2004. [31] NunzioTuccitto, Gabriella Zappalà, Stefania Vitale, AlbertoTorrisi, and Antonino Licciardello. A transport and reaction model for simulating cluster secondary ion massspectrometry depth proflesof organic solids. The Journal of Physical Chemistry C, 120(17):9263–9269, April 2016. [32] Robert J. Paruch, Zbigniew Postawa, Andreas Wucher, and Barbara J. Garrison. Steady-state statistical sputtering model for extracting depth profles from molecular dynamics simulations of dynamic SIMS. The Journal of Physical Chemistry C, 116(1):1042–1051, December 2011. [33] Robert J. Paruch, Barbara J. Garrison, and Zbigniew Postawa. Partnering analytic models and dynamic secondary ion mass spectrometry simulations to interpret depth profles due to kiloelectronvolt cluster bombardment. Analytical Chemistry, 84(6):3010–3016, March2012. [34] Robert J.Paruch, BarbaraJ. Garrison, and ZbigniewPostawa. Computed molecular depth profle for c60 bombardment of a molecular solid. Analytical Chemistry, 85(23):11628–11633, November 2013. [35] S. Sun, A. Wucher, C. Szakal, and N. Winograd. Depth profling of polycrystalline multilayers using aBuckminsterfullerene projectile. Applied Physics Letters, 84(25):5177–5179, June 2004. [36] Greg Gillen, MarlonWalker, Phillip Thompson, andJoe Bennett. Use of anSF[sub 5][sup +] polyatomic primary ion beam for ultrashallow depth profling on an ion microscope secondary ion massspectroscopyinstrument. Journal ofVacuum Science &Technology B: Microelectronics andNanometer Structures, 18(1):503, 2000. [37] S. Klaumunzer, Q.Q. Zhu, W. Schnabel, and G. Schumacher. Ion-beam-induced crosslinkingofpolystyrene stillan unsolved puzzle. Nuclear Instruments and Methods in PhysicsResearchSectionB:Beam Interactionswith Materials andAtoms, 116(14):154–158, August 1996. [38] D. Rading,R.Moellers, H.-G. Cramer,andE. Niehuis. Dualbeam depth proflingof polymer materials: comparisonof c60and ar cluster ionbeams for sputtering.Surface and InterfaceAnalysis, 45(1):171–174, July 2012. [39] R. Havelund, A. Licciardello, J. Bailey, N.Tuccitto, D. Sapuppo, I. S. Gilmore, J. S. Sharp, J. L. S. Lee, T. Mouhib, and A. Delcorte. Improving secondary ion mass spectrometry c60n+ sputter depth proflingofchallengingpolymers with nitricoxide gas dosing. Analytical Chemistry, 85(10):5064–5070, May2013. [40] Hua Tian, Dawid Maci¡»ek, Zbigniew Postawa, Barbara J. Garrison, and Nicholas Winograd. CO2 cluster ionbeam, an alternative projectile for secondary ion mass spectrometry. Journalof TheAmericanSociety for MassSpectrometry, 27(9):1476– 1482, June 2016. [41] Scott A Norris, Michael P Brenner, and Michael J Aziz. From crater functions to partial di˙erential equations:a new approach toionbombardment induced nonequilibrium pattern formation. Journal of Physics: Condensed Matter, 21(22):224017, May2009. [42] Matt P. Harrison and R. Mark Bradley. Crater function approach to ion-induced nanoscale pattern formation: Craters for fat surfaces are insuÿcient. PhysicalReview B, 89(24), June 2014. [43] Zhangcan Yang, Michael A. Lively, and Jean Paul Allain. Kinetic monte carlo simulationof self-organized pattern formationinducedby ionbeam sputtering using crater functions. PhysicalReviewB, 91(7),February 2015. [44] G. Hobler. Monte carlo simulationoftwo-dimensional implanted dopantdistributions at mask edges. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials andAtoms, 96(1-2):155–162, March1995. [45] J. F. Ziegler. Ion implantation : science and technology. Academic Press, Orlando, 1984. [46] G. Hobler and G. Betz. On the useful range of application of molecular dynamics simulations in the recoil interaction approximation. Nuclear Instruments and Methods in PhysicsResearchSectionB:Beam Interactionswith Materials andAtoms, 180(14):203–208, June 2001. [47] D.W. Oh, S.K. Oh, H.J. Kang, H.I. Lee, and D.W. Moon. In-depth concentration distribution of ar in si surface after low-energy ar+ ion sputtering. Nuclear Instruments and Methods in PhysicsResearch Section B: Beam Interactions with Materials andAtoms, 190(1-4):598–601, May2002. [48] Mark Robinson and IanTorrens. Computer simulationof atomic-displacement cascades in solids in the binary-collision approximation. Physical Review B, 9(12):5008– 5024, June 1974. [49] S. Lindsey and G. Hobler. Sputtering of silicon at glancing incidence. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials andAtoms, 303:142–147, May2013. [50] Wico C L Hopman, Feridun Ay, Wenbin Hu, Vishwas J Gadgil, Laurens Kuipers, MarkusPollnau, and RenéMde Ridder.Focused ionbeam scan routine,dwell time and dose optimizations for submicrometreperiodplanar photonic crystal components and stamps in silicon. Nanotechnology, 18(19):195305, apr 2007. [51] E. Niehuis, R. Mlers, D. Rading, H.-G. Cramer, and R. Kersting. Analysis of organic multilayers and 3d structures using ar cluster ions. Surface and Interface Analysis, 45(1):158–162, July 2012. Rozdział3 Przedruki artykułw pubs.acs.org/JPCC
Micro-and Macroscopic Modeling of Sputter Depth Profiling Dawid Maciazek,†
Robert J. Paruch,‡
Zbigniew Postawa,†
and Barbara J. Garrison*,‡
†Smoluchowski Institute of Physics, Jagiellonian University, ulica Lojasiewicza 11, 30-348 Krakow, Poland ‡Department of Chemistry, Penn State University, 104 Chemistry Building, University Park, Pennsylvania 16802, United States *Supporting
Information
S
ABSTRACT: A model for predicting depth profiles due to energetic particle bombardment based on the RMS roughness of the system and the sputtering yield is proposed. The model is an extension of the macroscopic transport model proposed previously [Tuccitto, N.; Zappala, G.; Vitale, S.; Torrisi, A.; Licciardello, A. J. Phys. Chem. C 2016, 120, 9263−9269]. The model is used to reconstruct the experimental depth profiles of a NiCr heterostructure due to bombardment by C60,SF5,O2, and Ga. 1. INTRODUCTION Efforts have been undertaken toward understanding the factors involved in depth profiling of atomic and molecular solids due to bombardment by energetic projectiles. Typical model systems include embedded delta layers and multilayer heterostructures. The first models introduced the three fundamental factors for depth profiling, that is information depth of sputtered material, ion-beam mixing, and surface roughness.1,2
These models propose analytic functions for fitting to experimental data and are useful in many applications. The procedure, however, of how to incorporate microscopic information on the atomic or molecular motion, such as comes from atomistic molecular dynamics (MD) simulations, is not clear. A model has been developed recently to take information from MD simulations and predict the depth profile of a delta layer.3−6
This model is valuable in providing insight as to how the physics from the MD simulations fits with the concepts of roughness, information depth, and mixing. The computational complexity of this model and the underlying MD simulations, however, make it prohibitive to consider large depths and long times. As a result, modeling depth profiling of systems with 100s of nm depth and inclusion of long time scale processes, such as thermal diffusion, is not tractable. It is also not clear how to connect experimental quantities with the calculated depth profile. Recently, a continuum transport and reaction model of depth profiling due to energetic ion-beam bombardment or sputtering has been proposed by Tuccitto et al.7
The physics associated with sputtering, specifically the sputtering yield and the ionbeam damage or displacements, are associated with the macroscopic processes of advection and diffusion, respectively. A reaction term has also been incorporated to include bombardment-induced chemistry. This transport model is capable of correctly reproducing qualitative features of a depth profile. It is not clear, however, how the parameters of this model map to the microscopic quantities of MD simulations or experimental data. On the other hand, incorporating quantities such as thermal diffusion and modeling large sample depths is straightforward and easy to compute. We have long been proponent of analytic models that incorporate input from MD simulations, especially for the ionbeam bombardment process of secondary ion mass spectrometry (SIMS). These include the mesoscale energy deposition footprint (MEDF) model8,9
and the steady-state statistical sputtering model (SS-SSM).4,5
In the MEDF model, the energy deposition profile from short-time MD simulations of cluster particle bombardment are used to provide input into an analytic model to predict sputtering yields as a function of incident energy. The SS-SSM uses input from repetitive bombardment MD simulations10
to predict the corresponding depth profile of a delta layer. For C60 and Au3 bombardment of a Ag surface, it was found that the best depth profiles resulted when the sputtering yield was high relative to the amount of ion-beam mixing or damage.5
In this study, we propose how to use the SS-SSM and underlying MD simulations to determine how to interpret and choose input into the transport model. Specifically we will use the examples of C60 and Au3 bombardment of Ag(111)4,5
and C60 bombardment of a molecular solid of octane.6
On the basis of this development, we propose simple relationships of the input quantities to the transport model based on only the experimental quantities of sputtering yield and RMS roughness. The predicted depth profiles from these parametrizations will be compared to the depth profiles predicted by the SS-SSM. Finally, we predict depth profiles for C60+,Ga+,SF5+, and O2+ depth profiling of a NiCr heterostructure and explain the beam type dependence of the experimental depth profiles.11−13
The Received: September 12, 2016 Revised: October 12, 2016 Published: October 19, 2016 © 2016 American Chemical Society 25473 DOI: 10.1021/acs.jpcc.6b09228
J. Phys. Chem. C 2016, 120, 25473−25480 The Journal of Physical Chemistry C ultimate target to be addressed in future studies is explaining the temperature dependence of depth profiling of Irganox delta layers.14−16
The ability to incorporate the experimental data directly into the transport model makes this model useful in direct interpretation of experimental results. 2. MICRO-AND MACROSCOPIC DESCRIPTIONS OF SPUTTER DEPTH PROFILING The first step is to establish a procedure for using the results of the repetitive bombardment MD simulations and the SS-SSM analysis to understand and quantify the input quantities of the transport model in terms of the basic physics and motion of the bombardment process. On the basis of this analysis and comparison, we formulate an empirical approach based on only experimental quantities for determining the input quantities of the transport model. 2.1. Transport and Reaction Model. This model is described in detail in the original reference7
and is based on the advection−diffusion−reaction equation for mass transfer above the average surface level, peaks in the roughened surface, as there are relatively few particles here. Second, the distributions from the average surface level (50% occupancy) and into the substrate are normalized by the occupancy of each layer. Thus, the distributions now appear as they are from a flat surface. The roughness information remains, however, in the widths of the distributions. In addition, the displacement distributions in the SS-SSM were for movement up and down. These have been converted to movement by one layer, by two layers, and by three layers, where each of these quantities is known as a function of layer position. Finally, the number of atoms per layer that do not move vertically are determined. This step is essential for determination of the value of the ionbeam mixing term. Of note is that for each layer the sum of the number of atoms that do not move, that move one layer, that move two layers, etc., is approximately proportional to the surface area of the master sample used in the MD simulation. Later, we will use the explicit surface area of the master sample for the calculation of time between impacts, and it is important to realize that this information will cancel out. ⎛⎜⎝ Dx () ⎞⎟⎠ + Rx () 2.2.1. Velocity. The velocity v that the mass is traveling with is the rate at which the surface recedes and is associated with ∂∂∂∂c c c =−+ v ∂∂x ∂∂ (1) t x x the sputtering yield, Y,in nm3 and experimental fluence, F,in ions/(nm2 s) as The independent variables are t = time [s] and x = depth [nm]. The dependent variable is c = c(x;t) = concentration profile of v=FY (2) mass (or number of atoms/molecules). The parameters are v = dvelocity the mass is traveling with [nm/s], D(x) = depth dependent diffusivity [nm2/s], and R(x) = depth dependent 2.2.2. Sampling Deepth profile is determ pth.
In
the
original
formulation,7
the
ined
only
by
the
concentration
of
the
dreaction term [1/s], a quantity not utilized in this study. The elta-layer component as it approaches the surface. In reality, hinitial condition is c0 = c(x;t=0). A second equation of the same owever, also particle are being s located below the surface eform, with its own set of the parameters, can be added to model jected by a single proj ectile impact. The depth distribution of a nthe depth profiling of the second component as is needed for umber of atoms ejected by a single impact is known as a sthe NiCr heterostructure. The initial concentration profile of a ampling or informati on depth. This quantity can be easily ddelta layer is modeled as a double sigmoid. It corresponds to a etermined from the s puttering distributions given by the MD silayer of mass embedded in the sample at given depth. As the mulations. We norm alize this function to a unit area as the total time passes, the mass travels toward the surface of the sample, sputtering yield is incorporated in the velocity. The nx = 0. During this process the profile changes its shape due to ormalized sampling as a depth, S(x), is subsequently used wthe effect of D(x). In the original model, the depth profile is eighting factor for the concentration profile, c(x), to calculate edetermined only by the concentration of the delta layer jected mass in the p redicted depth profile. In addition to cscomponent
at
the
surface
as
time
is
progressed
in
the
numerical
integration.7
Below
we
will
propose
a
method
to
include
ontaining information puttered particles, it i about the depth of origin of the mplicitly contains information about the RMS roughness. information depth in the transport model based on MD simulation results. 2.2.3. Ion-Beam M ixing and Diffusion. The ion-beam m2.2. Information Available in MD Simulations and the ixing and diffusion h as the form of SS-SSM.
The
first
step
is
to
assess
the
information
available
in
the
MD
simulations
and
the
SS-SSM4−6
analysis
to
make
a
=+D Dx D () total diff (3) wconnection to the transport model. The critical quantities that here D(x) is the de pth-dependent ion-beam mixing term. Tcarise
in
a
straightforward
analysis
from
the
repetitive
bombard
ment
MD
simulations10
are
the
yield
expressed
in
volume
per
impact
and
the
RMS
roughness.
The
SS-SSM4−6
provides
a
rue thermal diffusion Conceptually, there omponent that can b is represented by Ddiff and is a constant is a challenge to tie ion-beam mixing to e temperature dependent. dformalism for analyzing the results of the MD simulation and iffusion. Diffusion is generally considered to be a long time palso predicts depth profiles of delta layers. The SS-SSM divides rocess whereas the io n-beam mixing occurs on a time scale up tthe MD system into layers and calculates quantities on a per o tens of picosecond s. Even by examining MD simulations, player basis. The surface is roughened thus all quantities are roviding a precise tim e scale is challenging. We start, however, brelative to the average surface level. The three quantities of y implementing the f ormula used for calculating a diffusion crelevance are the sputtering yield per layer, displacements of onstant in MD simu an infinite onelations. Diffusivity for datoms or molecules from one layer to another, and the average irectional system is d escribed as occupation of each layer. The sputtering, displacement, and occupancy distributions as =⟨1 disp( ) ⟩x 2 determined in the SS-SSM from the MD simulations are for a →∞Dx () 2 lim t t (4) roughened surface, whereas the transport model assumes a flat Tsurface. The first step in making these two approaches he sputtering process and associated ion-beam mixing do not ecompatible is to ignore the information from the volume xactly fit with the ab ove description, but in order to find a 25474 DOI:
10.1021/acs.jpcc.6b09228
J. Phys. Chem. C 2016, 120, 25473−25480 The Journal of Physical Chemistry C connection between the macroscopic model for depth profiling and the microscopic motions we assume a similar relation as 1 ⟨x ⟩ disp( )2 Dx = () 2 Δt (5) where Δt is a time interval and thus we proceed from here. To calculate the quantity ⟨disp(x)2⟩, we use the information contained in the displacement distributions from the SS-SSM analysis as described above modified to account for the flat surface. For each depth we calculate the number of atoms that are displaced by one layer, by two layers, three layers, etc., and the number of atoms that are not displaced. The displacement term is thus a sum of the number of atoms that are displaced one layer times the layer spacing squared plus the number of atoms displaced two layers times twice the layer spacing squared, etc., divided by the total number of atoms. As noted above, the total number of atoms per layer is approximately proportional to the area of the master sample used in the MD simulations. The quantity Δt is assumed to be the time between impacts on the master sample or the reciprocal of the area of the master sample times the fluence F. The final situation is that we have used the SS-SSM displacement distributions to calculate a D(x) value. The time factor has an area of the master sample which effectively cancels out the total number of atoms in a layer in the MD simulation, a factor that is not physically relevant. Overall, then we can write D(x)as Dx Dx (6) () =F′() where D′(x) has units of nm4 per impact. This quantity D′(x) does not depend on the fluence, the time between impacts, or the size of the master sample. It is a quantity related to an average over single impacts much like the sputtering yield is an average over individual impact events. Now the ion-beam mixing term in the macroscopic transport equation has a microscopic quantity associated with it. At this point we are not prepared to make a precise physical interpretation of the D′(x) term. What is clear, however, is that D′(x) only depends upon properties associated with the individual impact events and not any long time diffusion-like behavior. This quantity should depend upon the projectile beam type, the incident kinetic energy, KE, of the projectile beam, and the material properties including displacement energies. Following the convention that sputtering yields17
are proportional to the ratio of the KE to the cohesive energy, U0, we assume the integral of D′(x) over depth to be proportional to KE/U0 when we discuss below the application of the transport model to the experimental NiCr heterostructures. To summarize, we used the displacement distributions from the MD simulations and the SS-SSM analysis to define a microscopic ion-beam mixing diffusion term, D(x). This relation results in the ion-beam mixing term to be a product of the fluence times a term that corresponds a quantity, D′(x), related to single impacts. For convenience, we assume the fluence to be 1 ion/(nm2 s) which corresponds to an ion current density of 6.4 nA per (200 μm)2. These conditions are in the experimental regime for C60 cluster sources. In fact, if there is no time-dependent quantity such as diffusion in the transport equation, the fluence does not matter in this model or in the experiment if the depth profile is shown as a function of depth rather than time. For comparisons with experimental data shown below, the apparent fluences are taken from the experimental data. 2.2.4. Verification of the Transport Model Depth Profiles vs the SS-SSM Depth Profiles. The first step is to define the S(x) and D′(x) quantities from the MD simulations and the SSSSM quantities, incorporate them into the transport model, and compare the depth profiles. Repetitive bombardment MD simulations have been performed for three systems and the results analyzed by the SS-SSM with the results presented in previous publications. • 20 keV C60/Ag(111).4−6,18
The sputtering yield is 6.4 nm3 (or 373 atoms) per impact,6
and the RMS roughness of the system is 2.35 nm.5
• 20 keV Au3/Ag(111).5,18
The sputtering yield is 3.69 nm3 (or 216 atoms) per impact, and the RMS roughness of the system is 2.69 nm. • 10 keV C60/octane.6
The total sputtering yield is 147 nm3 (or the equivalent of 703 molecules),6
and the RMS roughness is 3.3 nm. Only the yields and RMS values are given explicitly here as these are the quantities most directly comparable to information that can be extracted in experiments. As the differential equation solver that we use requires an analytical form of the input distributions sigmoidal functions are fit to the sampling depth distributions, S(x), as well as the ion-beam mixing terms, D′(x), obtained from the SS-SSM model. The results are shown in Figure
1
with the sigmoid function equation given in the Supporting
Information
and the parameters given in Table
1. The sampling distributions S(x) shown in Figure
1a are normalized to unit area since they are weighting functions used to obtain sputtered signal from concentration profiles. The distributions for the three systems are similar to basically the Figure 1. Distributions for use in the transport model for the three systems. (a) Sampling depth, S(x), from SS-SSM values and the RMS model. (b) Ion-beam mixing, D′(x), from SS-SSM values and RMS model; right-hand side scale is for octane, and left-hand side scale is for Ag. 25475 DOI: 10.1021/acs.jpcc.6b09228
J. Phys. Chem. C 2016, 120, 25473−25480 The Journal of Physical Chemistry C Table 1. Parameters for the Transport Model Based on the MD Results and SS-SSM Analysis sampling depth, S(x) ion-beam mixing, D′(x) velocity (nm/s) amplitude slope (1/nm) inflection (nm) amplitude (nm4/impact) slope (1/nm) inflection (nm) 20 keV C60/Ag(111) 20 keV Au3/Ag(111) 10 keV C60/octane 6.4 3.7 147 0.94 0.74 0.42 1.02 1.05 1.03 0.66 1.09 2.28 13.5 12.2 105.8 0.88 0.64 0.65 3.48 4.50 3.90 entire curve at distances less than the twice the RMS roughness. The ion-beam mixing, D′(x), functions are shown in Figure
1b. The ion-beam mixing curves are broader than the sampling distributions, consistent with the information from the MD simulations and SS-SSM analysis which shows that the ionbeam mixing occurs at a deeper depth than the sputtering process. The ion-beam mixing distributions are similar in magnitude for C60 and Au3 bombardment of Ag and a factor of 7 larger for C60 bombardment of octane. Explaining the large difference is beyond the scope of this study as we need to develop a better physical or computational picture of D′(x). Even though the RMS roughness is not explicitly included in the model, information about the RMS roughness does implicitly appear in the distributions in Figure
1. The depth profiles calculated with the SS-SSM, which were interpreted previously,5,6
are shown in Figure
2a. In the Figure 2. Depth profiles for a delta layer, represented by the blue vertical bar, of thickness 1 nm at a depth of 13.2 nm in the solid: (a) SS-SSM; (b) transport model with input from MD results and from the RMS model. transport model, the delta layer is represented by a double sigmoid as described in the Supporting
Information. The depth profiles calculated by the transport model using input from the MD results with SS-SSM analysis are shown in Figure
2b. The SS-SSM distributions do not have an associated time scale; thus, the transport model distributions are shown as a function of depth with the velocity being the conversion factor. In our opinion, there is qualitative agreement between the depth profiles predicted by the transport model and those from the SS-SSM. The worst agreement is for the octane system in which the leading edge of the depth profile from the transport equation is too narrow. In examining the SS-SSM sputtering distribution for this system,6
it is apparent that a significant amount of material is ejected from above the average surface level, and the assumption stated earlier that we can ignore the peaks of the roughened surface is not valid. 2.2.5. Development of an Empirical RMS Model. We feel that we have shown the correspondence between the terms in the transport equation and the microscopic quantities of the MD simulations. Our approach of using information from repetitive bombardment MD simulations provides the quantitative link and the conceptual definition of the impact related ion-beam mixing term, D′(x). This approach, however, is not practical for use with experimental data. Consequently, using insights from the previous calibration studies, we developed a simple model that uses the RMS roughness value and the sputtering yield, both quantities that are available from MD simulations or experimental data. This model, called the RMS model for input to the transport equation, is based on following assumptions. • The inflection point in the D′(x) distribution is proportional to the RMS value. • The inflection point in the S(x) distribution is smaller than D′ inflection point by a constant value, d, since the analyses of the MD results by the SS-SSM show this to be the case for the model systems. • The area under the curve of D′(x) is given, for now, by the area under the curve from the MD results. This assumption is the main area needed for improvement of the RMS model and, in fact, understanding the basic physics underlying D′(x). • For each the S(x) and D′(x) sigmoid functions, ignoring the amplitude, we assume that 90% of the amplitude change of the function is between distances of zero and twice the inflection point, a distance denoted “cut”. Since the sigmoid function is antisymmetric around the Table 2. Parameters for the Transport Model Based on the RMS Model sampling depth, S(x) ion-beam mixing, D′(x) velocity (nm/s) amplitude slope (1/nm) inflection (nm) amplitude (nm4/impact) slope (1/nm) inflection (nm) 20 keV C60/Ag(111) 6.4 0.72 0.89 1.01 13.46 0.89 3.31 20 keV Au3/Ag(111) 3.7 0.54 0.77 1.50 12.18 0.77 3.81 10 keV C60/Octane 147 0.39 0.64 2.28 105.8 0.64 4.58 25476 DOI: 10.1021/acs.jpcc.6b09228
J. Phys. Chem. C 2016, 120, 25473−25480 The Journal of Physical Chemistry C Table 3. Experimental Conditions sampling depth D′(x) RMS (nm) velocity (nm/s) Y (nm3) amplitude slope (1/nm) inflection (nm) amplitude (nm4/impact) slope (1/nm) inflection (nm) 15 keV C60+/NiCr 2.5 0.0518 2.95 0.628 0.835 1.23 6.78 0.835 3.53 15 keV Ga+/NiCr 100 0.0122 0.174 0.007 07 0.0209 138.7 0.2 0.0209 141 3 keV SF5+/NiCr 10 0.201 0.174 0.082 0.209 11.8 0.181 0.209 14.1 3 keV O2+/NiCr 62 0.441 0.087 0.011 0.0337 85.1 0.0644 0.0337 87.4 inflection point, this means that at a distance of twice the inflection point the value of the function should be 0.05 and the value at zero should be 0.95. This allows us to determine the slope to be ln(1/cut −1) slope =inflection point (7) We have found that a good agreement between the depth profiles predicted by this simplified RMS model and the depth profiles using MD input can be obtained, if the inflection point for D′(x) = 1.4·RMS and d = 3.2 nm as shown in Figure
2b. The calculated sigmoid function parameters of the RMS model for the three systems are given in Table
2, and the S(x) and D′(x) distributions are shown in Figure
1. Regardless of its simplicity, the RMS model reproduces the depth profiles of the SS-SSM method as well as the functions based on the MD results. We have, of course, used information from the MD results to give the area under the D′(x) curve for the RMS model. The effect of thermal diffusion on the Ag depth profiles can be easily tested in the transport model by setting the Ddiff parameter. For Ag, self-diffusion has been measured over the temperature ranges of 454−777 K19
and 903−1208 K,20
and temperature-dependent diffusion constants have been determined. If these relations are extrapolated to 300 K, the diffusion constant is on the order of 10−18 nm2/s. This value can be added to the calculation, but this value is 18 orders of magnitude smaller than the ion-beam mixing values; thus, it is obvious that no significant thermal diffusion occurs in Ag. The big point, however, is that both ion-beam mixing and thermal diffusion can be accommodated within the transport model. 3. INTERPRETATION OF EXPERIMENTAL DATA We now take the RMS model one step further and try to reproduce experimental depth profiles of NiCr heterostructures. It will be an important test of model versatility as these depth profiles were obtained with four very different beam conditions. The NiCr heterostructure sample has nine alternating layers of Cr and Ni with the five Cr layers 53 nm thick and the four Ni layers 66 nm thick for a total depth of 529 nm. Two groups have performed depth profiling experiments on this standardized sample. First, Gillen et al. have used 3 keV ++ SF5 and O2 beams and measured the Ni+ and Cr+ signals.13
Second, Sun et al. have performed experiments using 15 keV C60+ and Ga+ beams and measured the Ni and Cr signals with multiphoton resonance ionization.11,12
The objective of this analysis is to use the available experimental data to predict the depth profiles using the transport model using the RMS model to generate input functions. The three main quantities in the model are the sputtering yield, the displacement quantity D′(x), and the fluence. The quantities that appear in the experiment are the time to depth profile through the 529 nm of material, the RMS roughness, and the sputtering yields for the C60+ and Ga+ beams. The experimental time to depth profile through the sample was used to calculate the average velocity in the model. The obtained values are given in Table
3. In the original + publications,11,12
for C60 and Ga+ bombardment, the yields were given in atoms per impact, and we converted these to nm3 by assuming the single crystal atomic density. For Ni and Cr both the atomic densities and yield are within 9% of each other so we used the averaged values for these quantities. The effective fluence is the velocity divided by the yield. This value accounts for the fluence of the beam and also the effect of rastering the beam to produce a depth profile. The yields for SF5+ and O2+ bombardment were determined by simulation using the experimental values for C60+ and Ga+ as reference points. The calculation of absolute sputtering yields is challenging because none of the interaction potentials are completely accurate and because some physics is missing such as energy loss to electronic effects. Thus, our strategy is to use comparable simulations to estimate relative yields and then determine an estimate of the experimental yields for SF5+ and O2+ bombardment based on the experimental C60+ and Ga+ yields. Previously we determined from MD simulations the yields for C60 and Ga bombardment of Ag at 15 keV to be 331 and 21 atoms per impact.21
No potentials exist for bombardment of SF5 and O2 on Ag, and we do not want to include too much potential variation so we tried to construct hydrocarbon analogues for these two projectiles. We chose projectiles of neopentane C(CH3)4 and C2 (with mass of O rather than C). One hundred impacts were calculated for each projectile on the flat surface at an energy of 3 keV and the experimental impact angle of 52°. The calculated yields are 21.6 ± 1.5 and 10.2 ± 1.0 atoms/impact for C(CH3)4 and C2, respectively. The neopentane calculated yield is almost the same as the Ga calculated yield, and the C2 yield is about half the Ga calculated yield on Ag; thus, we scale the experimental yields for Ga+ on NiCr for SF5+ and O2+ as given in Table
3. The effective fluence is calculated from the experimental average velocity. The RMS model is used with the experimental RMS values to determine the S(x) and D′(x) functions as shown in Figure
3. For the area under the D′(x) curve we feel that the appropriate calculation should be for the metal substrate, and since D′(x) for both the C60 and Au3 projectiles have similar areas, we chose the 20 keV Au3/Ag system, which has a value of 55.8 nm5. This value was scaled for each set of beam conditions by KE/U0, as it is commonly used to scale sputtering yields.22
We assume that it is also a logical scaling factor for displacements. The value of U0 for Ag is 2.93 eV,23
and the average value for the NiCr system is 4.27 eV.23
The four experimental conditions break into two groups depending on the RMS value. The RMS values for the C60+ and SF5+ bombardment are 2.5 and 10 nm, respectively. For these two systems the axes at the left and bottom of Figure
3
should be 25477 DOI: 10.1021/acs.jpcc.6b09228
J. Phys. Chem. C 2016, 120, 25473−25480 The Journal of Physical Chemistry C used. The S(x) and D′(x) functions are relatively short-ranged and intense. In contrast, the RMS values for O2+ and Ga+ have larger RMS values of 62 and 100 nm, respectively. The right and top axes should be used for their S(x) and D′(x) functions. These functions are longer-ranged and less intense than the ones for C60+ and SF5+. Of note is that the NiCr layers are 53 and 66 nm wide, distances larger than the C60+ and SF5+ functions, smaller than the Ga+ functions, and comparable to the O2+ functions. The predicted depth profiles are shown in Figure
4
overlaid on the experimental data that was digitized from the original literature. The double sigmoid used to describe the NiCr heterostructure is described in the Supporting
Information.We have maintained the plot representation, linear vs log scale, of the experimental publications. The intensities of the individual Ni and Cr signals have a multitude of experimental factors so we used two normalization factors per plot. Overall, we are pleased with the agreement between the transport model depth profiles using the RMS model for predicting the S(x) and D′(x) functions with the experimental data. The RMS roughness values for C60+ and SF5+ bombardment are less than the heterostructure layers, and consequently the depth profiles show discrete layers. For the O2+ and Ga+ bombardment there is an induction time to develop the large RMS roughness values, a factor not included in the transport model. Even with this omission, we feel that the predicted depth profiles qualitatively reproduce the experimental distributions. We did try not using KE/U0 to scale the D′(x) area. The resulting depth profiles were slightly broader but not qualitatively different. There is one additional piece of information available in the + experimental study on the C60 bombardment of NiCr;11
Figure 4. Dependence of signals of ions (a, b) and neutral atoms (c, d) on the sputter time for a nine-layer Ni:Cr multilayer stack bombarded by (a) 3 keV SF5+, (b) 3 keV O2+, (c) 15 keV C60+, and (d) 15 keV Ga+. Black curves depict experimental signal of Cr (solid line) and Ni (dash-dotted line) atoms, while red curves depict analogous data obtained from the RMS model. The dashed black line depicts experimental intensity of the substrate Si signal which was not modeled. 25478 DOI: 10.1021/acs.jpcc.6b09228
J. Phys. Chem. C 2016, 120, 25473−25480 The Journal of Physical Chemistry C namely, the intensity of the NiCr dimer is measured as a function of time. The emission of these dimers should only occur where both Ni and Cr are present in the sample, that is, at the interfaces between the layers. The experimental data clearly show that the widths of the NiCr peaks increase as the material is depth profiled. This observation is also slightly visible in the experimental depth profile shown in Figure
4c. The transport model does not show the same broadening, and in fact, it cannot since the S(x) and D′(x) functions become vanishingly small at about one-half the layer width. This discrepancy implies to us that some experimental condition may not be fully controlled or understood. One possibility could be, for instance, nonuniform material removal, which would result in a formation of a crater with a bottom that becomes less parallel to the sample surface with time. 4. CONCLUSIONS Using results from molecular dynamics simulations and the steady-state statistical sputtering model for predicting depth profiles due to energetic particle bombardment, we have defined the quantities in a macroscopic model for predicting depth profiles based on the transport equation. Specifically, we developed a microscopic definition for the ion-beam mixing or diffusion term, a quantity that is associated with the amount of displacements in the individual impact event. In addition, we proposed a protocol to incorporate sampling or information depth calculated from individual impacts, a term that was not present in the original model. The transport model for depth profiling was calibrated for three test cases for which there are associated repetitive bombardment MD simulations and predicted depth profiles from the SS-SSM. Using the functions developed for these calibration studies, we developed a simple RMS model that predicts the functions for the transport model from just the RMS roughness value and the sputtering yield. We show that we can reproduce the depth profiles of a NiCr heterostructure for four beam conditions that range from RMS values of 2.5 to 100 nm. In addition, we interpret the observed broadening of the depth profile for C60 bombardment, the system with the smallest RMS value of 2.5 nm, to be due to some unknown experimental condition. The main open issue is how to define theoretically the impact related ion-beam mixing term and develop a protocol for calculating it. *Supporting Information ■ ASSOCIATED CONTENT S The Supporting Information is available free of charge on the ACS
Publications
website
at DOI: 10.1021/acs.jpcc.6b09228. Procedure used to solve eq
1;definition of sigmoid functions used in our study (PDF) ■ AUTHOR INFORMATION Corresponding Author *E-mail bjg@psu.edu. Notes The authors declare no competing financial interest. ■ ACKNOWLEDGMENTS The authors gratefully acknowledge the financial support from the National Science Foundation, Grant CHE-1212645, and the Polish National Science Center, Program 2013/09/B/ST4/ 00094. ■ REFERENCES (1) Hofmann, S. Atomic Mixing, Surface Roughness and Information Depth in High-Resolution AES Depth Profiling of a GaAs/AlAs Superlattice Structure. Surf. Interface Anal. 1994, 21, 673−8. (2) Dowsett, M. G.; Rowlands, G.; Allen, P. N.; Barlow, R. D. Analytic Form for the SIMS Response Function Measured from Ultra-Thin Impurity Layers. Surf. Interface Anal. 1994, 21, 310−315. (3) Krantzman, K. D.; Wucher, A. Fluence Effects in C60 Cluster Bombardment of Silicon. J. Phys. Chem. C 2010, 114, 5480−5490. (4) Paruch, R. J.; Postawa, Z.; Wucher, A.; Garrison, B. J. Steady-State Statistical Sputtering Model for Extracting Depth Profiles from Molecular Dynamics Simulations of Dynamic SIMS. J. Phys. Chem. C 2012, 116, 1042−1051. (5) Paruch, R. J.; Garrison, B. J.; Postawa, Z. Partnering Analytic Models and Dynamic Secondary Ion Mass Spectrometry Simulations to Interpret Depth Profiles Due to Kiloelectronvolt Cluster Bombardment. Anal. Chem. 2012, 84, 3010−6. (6) Paruch, R. J.; Garrison, B. J.; Postawa, Z. Computed Molecular Depth Profile for C60 Bombardment of a Molecular Solid. Anal. Chem. 2013, 85, 11628−11633. (7) Tuccitto, N.; Zappala, G.; Vitale, S.; Torrisi, A.; Licciardello, A. A Transport and Reaction Model for Simulating Cluster Secondary Ion Mass Spectrometry Depth Profiles of Organic Solids. J. Phys. Chem. C 2016, 120, 9263−9269. (8) Russo, M. F.; Szakal, C.; Kozole, J.; Winograd, N.; Garrison, B. J. Sputtering Yields for C60 and Au3 Bombardment of Water Ice as a Function of Incident Kinetic Energy. Anal. Chem. 2007, 79, 4493− 4498. (9) Russo, M. F.; Garrison, B. J. Mesoscale Energy Deposition Footprint Model for Kiloelectronvolt Cluster Bombardment of Solids. Anal. Chem. 2006, 78, 7206−7210. (10) Russo, M. F., Jr.; Postawa, Z.; Garrison, B. J. A Computational Investigation of C60 Depth Profiling of Ag: Molecular Dynamics of Multiple Impact Events. J. Phys. Chem. C 2009, 113, 3270−3276. (11) Sun, S.; Wucher, A.; Szakal, C.; Winograd, N. Depth Profiling of Polycrystalline Multilayers Using a Buckminsterfullerene Projectile. Appl. Phys. Lett. 2004, 84, 5177−9. (12) Sun, S.; Szakal, C.; Roll, T.; Mazarov, P.; Wucher, A.; Winograd, N. Use of C60 Cluster Projectiles for Sputter Depth Profiling of Polycrystalline Metals. Surf. Interface Anal. 2004, 36, 1367−72. (13) Gillen, G.; Walker, M.; Thompson, P.; Bennett, J. Use of an SF5+ Polyatomic Primary Ion Beam for Ultrashallow Depth Profiling on an Ion Microscope Secondary Ion Mass Spectroscopy Instrument. J. Vac. Sci. Technol., B: Microelectron. Process. Phenom. 2000, 18, 503− 508. (14) Mao, D.; Wucher, A.; Brenes, D. A.; Lu, C. Y.; Winograd, N. Cluster Secondary Ion Mass Spectrometry and the Temperature Dependence of Molecular Depth Profiles. Anal. Chem. 2012, 84, 3981−3989. (15) Mao, D.; Lu, C.; Winograd, N.; Wucher, A. Molecular Depth Profiling by Wedged Crater Beveling. Anal. Chem. 2011, 83, 6410− 6417. (16) Shard, A. G.; et al. Argon Cluster Ion Beams for Organic Depth Profiling: Results from a VAMAS Interlaboratory Study. Anal. Chem. 2012, 84, 7865−7873. (17) Nastasi, M.; Mayer, J.; Hirvonen, J. K. Ion-Solid Interactions: Fundamentals and Applications; Cambridge University Press: Cambridge, 1996; p 145. (18) Paruch, R.; Rzeznik, L.; Russo, M. F.; Garrison, B. J.; Postawa, Z. Molecular Dynamics Study of the Effect of Surface Topography on Sputtering Induced by 20 keV Au3 and C60 Clusters. J. Phys. Chem. C 2010, 114, 5532−5539. (19) Lam, N. Q.; Rothman, S. J.; Mehrer, H.; Nowicki, L. J. Self-Diffusion in Silver at Low-Temperatures. Phys. Status Solidi B 1973, 57, 225−236. (20) Tomizuka, C. T.; Sonder, E. Self-Diffusion in Silver. Phys. Rev. 1956, 103, 1182−1184. (21) Postawa, Z.; Czerwinski, B.; Szewczyk, M.; Smiley, E. J.; Winograd, N.; Garrison, B. J. Enhancement of Sputtering Yields Due 25479 DOI: 10.1021/acs.jpcc.6b09228
J. Phys. Chem. C 2016, 120, 25473−25480 The Journal of Physical Chemistry C to C60 Versus Ga Bombardment of Ag{111} as Explored by Molecular Dynamics Simulations. Anal. Chem. 2003, 75, 4402−4407. (22) Paruch, R. J.; Postawa, Z.; Garrison, B. J. Seduction of Finding Universality in Sputtering Yields Due to Cluster Bombardment of Solids. Acc. Chem. Res. 2015, 48, 2529−36. (23) Kittel, C. Introduction to Solid State Physics, 8th ed.; John Wiley & Sons, Inc.: Hoboken, NJ, 2005; p 50. 25480 DOI: 10.1021/acs.jpcc.6b09228
J. Phys. Chem. C 2016, 120, 25473−25480 Nuclear
Instruments
and
Methods
in
Physics
Research
B
393
(2017)
17–21
Contents lists available at ScienceDirect
Nuclear Instruments and Methods in Physics Research B journal
homepage:
www.elsevier.com/locate/nimb
Simple model of surface roughness for binary collision sputtering simulations ˛ _ Sloan J. Lindsey a, Gerhard Hobler a,⇑, Dawid Maciazek b, Zbigniew Postawa b
a Institute of Solid-State Electronics, TU Wien, Floragasse 7, A-1040 Wien, Austria b Institute of Physics, Jagiellonian University, ul. Lojasiewicza 11, 30348 Krak, Poland article info abstract
Article history: Received 19 July 2016 Received in revised form 16 September 2016 Accepted 29 September 2016 Available online 7 October 2016 Keywords: Sputtering Surface roughness Binary collision simulations Molecular dynamics It has been shown that surface roughness can strongly influence the sputtering yield – especially at glancing incidence angles where the inclusion of surface roughness leads to an increase in sputtering yields. In this work, we propose a simple one-parameter model (the ‘‘density gradient model”) which imitates sur-face roughness effects. In the model, the target’s atomic density is assumed to vary linearly between the actual material density and zero. The layer width is the sole model parameter. The model has been implemented in the binary collision simulator IMSIL and has been evaluated against various geometric surface models for 5 keV Ga ions impinging an amorphous Si target. To aid the construction of a realistic rough surface topography, we have performed MD simulations of sequential 5 keV Ga impacts on an initially crystalline Si target. We show that our new model effectively reproduces the sputtering yield, with only minor variations in the energy and angular distributions of sputtered particles. The success of the density gradient model is attributed to a reduction of the reflection coefficient – leading to increased sputtering yields, similar in effect to surface roughness.  2016 Elsevier B.V. All rights reserved. 1. Introduction It is well known that ion bombardment roughens the target sur-face [1], which in turn may influence the sputtering yield [2].Ina recent study [3]
we demonstrated that the simulation of sputtering yields for grazingly incident ions requires the consideration of sur-face roughness. Grazing incidence conditions are typically found during transmission electron microscopy (TEM) sample preparation, one of the most important applications of focused ion beams (FIB) [4,5]. Glancing incidence angles may also occur during FIB milling of holes [6]
or during irradiation of nanowires, for instance, when the nanowires bend towards the beam [7]. Recently, several groups have developed Monte Carlo binary collision (BC) codes that are capable of simulating ion bombardment of 2D and 3D micro-and nanostructures [8–14]. They all lack models of surface roughness, prohibiting meaningful simulation of glancing angle effects. BC simulation studies of the effect of surface roughness on sputtering have mostly used geometrically modified flat surfaces, employing square waves [2], sinusoidal waves [3], or by applying fractal surface models [15,16]. All these models would be difficult ⇑ Corresponding author. E-mail address: gerhard.hobler@tuwien.ac.at
(G. Hobler). to implement in 2D and 3D simulations, since the simulation needs to store a rough surface geometry interspersed with any mesoscale topographies of interest. Yamamura
et
al.
[2]
have shown that reducing the target den-sity in a surface layer has a similar effect as a geometrically defined rough surface. However, they used a monoatomic surface layer only, which limited the degree of roughness that could be introduced. In the present work, we generalize Yamamura’s idea by using a surface layer with a density that decreases linearly towards the surface, which we henceforth call density gradient model. The sole parameter of this model is the thickness of the layer, which may be larger than monoatomic. This removes ambiguity from the model as compared to, e.g., the sinusoidal model where fits to the experimental data are available across multiple pairings of wavelengths and amplitudes [3]. Equally important, the density gradient model is easily implemented in 2D and 3D BC simulations. To
validate
the
model,
we
compare
sputtering
yields
and
energy
and angular distributions of sputtered atoms with the predictions of three geometrically defined models with wave vectors parallel to the projection of the ions’ incidence direction to the surface. To aid the construction of a realistic rough surface, we have also performed molecular dynamics (MD) simulations of sequential Ga impacts on a Si surface. This
study
is
carried
out
for
5
keV
Ga
ions
impinging
on
Si.
The
simulations are compared to experimental data obtained at FEI http://dx.doi.org/10.1016/j.nimb.2016.09.028
0168-583X/ 2016 Elsevier B.V. All rights reserved. S.J. Lindsey et al. / Nuclear Instruments and Methods in Physics Research B 393 (2017) 17–21 Company [17]. A lower ion energy is used than in our previous study [3]
mainly to facilitate the MD simulations. 2.
Simulation specifics 2.1.
MD modeling MD modeling was carried out at Jagiellonian University using a modified version of LAMMPS [18]. The simulation cell was 120 Å 120 Å 80 Å initially filled with single-crystalline (100)Si with a (2 1) reconstructed surface. The number of Si atoms in this cell was 57112. Periodic boundary conditions were used in the lateral directions. Stochastic and rigid layers, 7 Å and 3 Å thick, respectively, were used at the bottom to simulate the thermal bath that kept the sample at the required temperature and to keep the shape of the sample. The simulations were run at 0 K temperature. The target was sequentially bombarded with 5 keV Ga ions at polar angles of 89 and 85 and an azimuthal angle of 35 with respect to the cell edge which is a [010] direction. The latter was chosen as to minimize possible artifacts of the periodic boundary conditions due to the passage of the ions over regions they have previously interacted with while slowly glancing off the surface. Each impact was simulated for 2 ps. The resulting structure was used as initial condition for the subsequent impact after removal of all sputtered atoms and any excess kinetic energy from the system. The latter was achieved by an energy quenching procedure that involved application of gentle viscous damping forces to the entire sample for 0.2 ps. A Tersoff-3 potential [19]
was used for Si-Si interactions, and the ZBL potential [20]
splined with a Lennard-Jones (LJ) potential for Ga-Si interactions. 2.2.
BC modeling Monte Carlo simulations using the BC approximation were carried out at TU Wien using the simulator IMSIL [21,22]. As in our earlier work [3]
we use the ZBL interatomic potential [20], the Oen– Robinson model for electronic stopping [23]
with a cutoff energy of 10 eV, and a planar surface potential with a surface binding energy between Si and Si of 4.7 eV. Since the refraction of the incident ions at the surface potential is significant under the conditions studied, the choice of the surface binding energy between Ga and Si is also critical. We use a value of 2.82 eV [24]. All of the BC simulations were carried out using the static mode of IMSIL, wherein the target starts in an amorphous state as pure silicon, and its modification by the implanted gallium is not taken into account over the course of the simulation. 25,000 impacts were carried out for sputter yield calculations and 5 million impacts for determining the energy and angular distributions of sputtered atoms. 2D geometries are specified in IMSIL by polygons which are converted to a signed distance function defined on a Cartesian grid covering the simulation domain [13]. IMSIL was adapted for this research through the addition of a periodic geometry mode, which allowed to set the lateral size of the simulation domain equal to one wavelength. The vertical size was chosen 300 Å plus the roughness amplitude, which led to a negligible forward sputtering yield of 2 104, thus indicating sufficient thickness to simulate an infinitely thick target. In order to exclude any significant discretization errors, we used 2000 segments for the polygons and 1 million cells for the internal grid. Three geometric surface models were used: Cosine, triangular, and double cosine (Fig.1). The double cosine function is defined by the superposition of two cosine functions with wavelengths differing by a factor of three: Fig. 1. Geometric roughness models used in this study, shown for a wavelength of 30 Å: Cosine, triangular, and double cosine. All surface models are shown with bestfit parameters. A 2p 6p zðxÞ¼ cos x þcos x; ð1Þ 2 kk where A is the amplitude and k is the period. The choice of the wavelength k is not very critical and will be estimated from the results of the MD simulations. The amplitude A will be determined by fitting to experimental sputtering yield data. The density gradient model is implemented by creating collision partners at the end of the free flight paths with probabilities 8>1 for d P w < p ¼ d=w for 0 < d < w ð2Þ > : 0 for d 6 0; where d denotes the signed distance from the surface (negative values outside the target) and w is the width of the density gradient layer. Effectively this means that the density of target atoms decreases linearly towards the surface within the layer 0 < d < w, as illustrated in Fig.
2a. The parameter w will be determined by fitting to the experimental sputtering yield data. Note that in the limit of zero wavelengths the geometric surface models correspond to reduced-density layers. These densities are shown in Fig.
2b and are compared with the density gradient model. Fig. 2. (a) Schematic drawing of the target density in the density gradient model. Filled circles represent atoms. (b) Comparison of the density with equivalent density profiles of the geometric roughness models. The height axis in (b) has been shifted with respect to (a) so its origin is at the center of each roughness layer. All surface models are shown with best-fit parameters. S.J. Lindsey et al. / Nuclear Instruments and Methods in Physics Research B 393 (2017) 17–21 3. Results MD simulations were performed for an incidence angle of 89 with different binding energies of the LJ potential. The best fit to the experimental sputter yield of 1.7 was obtained with a binding energy of 0.9 eV. In these simulations, the sputtering yield was averaged between the 1100th impact, when the sputtering yield had stabilized, and the 2191st (last) impact. With the same binding energy, impacts at an incidence angle of 85 were simulated resulting in a sputtering yield of 6.49, which compares favorably with the experimental value of 6.66, thus giving confidence in the simulations. The RMS roughness amplitudes of the surface after the last impact are 2.8 Å and 3.1 Å for incidence angles of 89 and 85, respectively. Inspection of the surface at the end of the simulation reveals intertrough/intercrest distances of approximately 20– 50 Å (Fig.
3). We therefore chose a wavelength of k ¼ 30Å for our geometric roughness models. BC simulations were then performed with all models, and least squares fitting was used to determine the best fit of the amplitudes A to the experimental sputtering yields at incidence angles of 86.6, 88.1, and 89, resulting in A ¼ 3:42Å, 3:12Å, and 3:04 Å for the cosine, triangle, and double cosine model, respectively. In the same way the best-fit layer width w ¼ 7:49 Å of the density gradient model was obtained. The surface models with these optimized parameters are represented in Figs.
1
and 2. Sputtering yields obtained with these models are shown in Fig.
4a and b, and are compared with the experimental data. Good fits are observed for all roughness models. Interestingly, the density gradient model features the best fit (see Fig.
4b). In contrast, large deviations between simulation results and experimental data are observed when a flat surface is assumed in BC simulations of incidence angles larger than 82. Comparing this to our earlier study of 30 keV Ga ions impinging on Si [3], where flat surface simulations have been found to deviate from the experimental data at angles larger than 86, it may be concluded that the range of conditions at which sputtering yields Fig.3.Topviewofthesiliconsampleafter2100sequentialhitsby5keVGaionssimulatedbyMDforapolarangleof89andanazimuthalangleof35.Atomsarecoloredaccordingtotheirzcoordinate;red:thepositionofthetopmostatoms,white:3Åbelowthetoplevel,andblue:6Åbelowthetoplevel.Theedgelengthofthecellshownis120Å.(Forinterpretationofthereferencestocolorinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.) Fig. 4. (a,b) Sputtering yields as a function of incidence angle: Experimental data [17]
(filled circles), simulation results for a flat surface (long-dashed line) and for the various roughness models with the amplitude fitted to the last three experimental data points, (b) is a magnification of the range between 80 and 90, (c) shows the reflection coefficients obtained with all models. are affected by surface roughness, increases with decreasing ion energy. The
reflection
coefficients
(the
fraction
of
the
incident
ions
that
are reflected or backscattered) in these simulations are shown as a function of incidence angle in Fig.
4c. The results indicate that all models that include surface roughness increase the sputtering yield by means of reducing the reflectivity of the target. The relation between sputtering yield and reflection coefficient is easily understood, as reflected ions interact weaker with the target and therefore produce fewer recoils than ions that are scattered into the material. While
it
is
reassuring
that
the
density
gradient
model
fits
the
experimental sputtering yield data over the whole range of incidence angles with a single value of the layer width, higher order yield statistics such as the energy and angular distribution of the sputtered atoms are also important, especially when the simulation data is used as an input to topography simulations [17,26,27]. Modeling of surface roughness clearly plays a role in the resulting energy and angular distributions as shown in Fig.
5a and b, respectively. For this comparison, energy and angular distributions have been simulated for an incidence angle of 87.7, the approximate crossover point of the sputtering yield curves obtained with the four roughness models (see Fig.
4b). Nearperfect agreement of the energy distributions is observed between all roughness models, while the energy distribution of atoms sputtered from the flat surface (black long-dashed line) is significantly different. For the angular distribution the agreement between the roughness models is less marked, but all four roughness models yield a more pronounced first-knock-on peak than the flat surface. This peak in the angular distribution is attributed to high-energy ejected target atoms that undergo relatively few collisions before leaving the target and thus retain more energy and a greater fingerprint of the incoming ion trajectory. The larger role of primary recoils in sputtering from a rough surface can also be induced from S.J. Lindsey et al. / Nuclear Instruments and Methods in Physics Research B 393 (2017) 17–21 Fig. 5. (a) Energy distributions and (b) angular distributions of sputtered particles obtained by BC simulations at an incidence angle of 87.7: Flat surface (long-dashed line) compared to the four roughness models. The gray dashed circle in (b) indicates a cosine distribution, which is predicted by theory [25]
for well developed cascades. the energy distribution which is weaker than E 2, since a dependence as E 2 is expected when the collision cascades are well developed [25]. Despite the favorable results presented so far, there is one discrepancy: The RMS amplitudes obtained from the MD simulations are larger than those of the geometric roughness models fitted to the experimental data. To shed light on a possible explanation, we have fitted the amplitudes of our three geometric roughness models as a function of the wavelength. The results for the ampli- Fig. 6. Best fits of roughness profile amplitudes to the experimental data as a function of assumed wavelength. While the upper three, thicker lines show the amplitude A, the lower three, thinner lines represent the respective RMS values. Half of the best-fit layer thickness w of the density gradient model is shown for comparison (horizontal line labelled ‘‘gradient”). tudes are presented in Fig.
6
(upper curves) together with the corresponding RMS values (thinner broken lines in the lower part of the figure). It can be seen that the amplitudes rise towards small wavelengths. This can tentatively be explained by the increasing role of redeposition: With decreasing wavelength the aspect ratio of the topography increases, which increases the amount of redeposition. Redeposited atoms are not sputtered, so to fit the model to the experimental sputtering yield, the roughness amplitude has to be increased. When using a roughness model with a wavelength k > 10 Å in the BC simulations, short-wavelength components that are present in the MD simulations are suppressed, therefore redeposition is underestimated, and smaller amplitudes are sufficient to fit the experimental data. It should be noted, however, that other factors may play a role such as the planar surface potential model used in the BC simulations that becomes questionable on a rough surface. We note that all wavelengths except for extremely small ones ( 1 Å) give sputtering yields that are in good agreement with the experimental data over the whole range of incidence energies (not shown). 4. Conclusions We have demonstrated that the newly proposed ‘‘density gradient” model can be used to fit BC simulations to experimental data on sputtering yields as a function of incidence angle. We have also shown that a variety of roughness models is capable of describing the data reasonably with the density gradient model providing the best fit in the particular case studied. That one and the same roughness model can describe sputtering conditions for all incidence angles, is not self-evident – different surface roughness could develop as a result of different incidence angles. Our MD simulations show that for 5 keV Ga bombardment of Si the RMS roughness amplitude is only slightly different for incidence angles of 85 and 89. We have not investigated surface roughness for more-perpendicular ion impact. However, moderate beaminduced surface roughness as investigated in this work, affects sputtering yields only weakly at most incidence angles (<82 for 5 keV Ga incident on Si). So even a discrepancy of the roughness model with the actual target topography at these less-oblique incidence angles would not inhibit successful simulation of sputtering effects. We have also shown that all roughness models investigated lead to almost identical energy distributions of the sputtered atoms and to similar angular distributions and reflection coefficients. It may be concluded that the main effect of a surface roughness model for sputtering simulations must be a reduction in the reflection coefficient, while the actual shape of the assumed sur-face is of less importance as long as the correct sputtering yield is fitted. This explains why a physically questionable model – densities approaching zero as in our density gradient model are not possible in condensed matter – may be so successful. The new density gradient model is computationally efficient and easy to implement in BC codes. This is true even for 2D and 3D topographies, provided a signed distance function field is available such as in our IMSIL simulator. It is hoped that the new roughness model enables a greater degree of realism in BC sputtering simulations with very little computational overhead. Acknowledgments MK
and
ZP
would
like
to
gratefully
acknowledge
financial
sup
port from the Polish National Science Center, Program No. 2013/09/B/ST4/00094. S.J. Lindsey et al. / Nuclear Instruments and Methods in Physics Research B 393 (2017) 17–21 References [1]
E.A.
Eklund,
R.
Bruinsma,
J.
Rudnick,
R.S.
Williams,
Phys.
Rev.
Lett.
67
(1991)
1759. [2]
Y.
Yamamura,
C.
Mossner,
H.
Oechsner,
Radiat.
Eff.
103
(1987)
25. [3]
S.
Lindsey,
G.
Hobler,
Nucl.
Instrum.
Methods
B
303
(2013)
142. [4]
L.
Giannuzzi,
F.
Stevie,
Micron
30
(1999)
197. [5]
J.
Mayer,
L.A.
Giannuzzi,
T.
Kamino,
J.
Michael,
MRS
Bull.
32
(2007)
400. [6]
W.C.L.
Hopman,
F.
Ay,
W.
Hu,
V.J.
Gadgil,
L.
Kuipers,
M.
Pollnau,
R.M.D.
Ridder,
Nanotechnology
18
(2007)
195305. [7]
M.
Bettge,
S.
MacLaren,
S.
Burdin,
R.T.
Haasch,
D.
Abraham,
I.
Petrov,
M.-F.
Yu,
E.
Sammann,
Nanotechnology
23
(2012)
175302. [8]
I.
Bizyukov,
A.
Mutzke,
R.
Schneider,
A.M.
Gigler,
K.
Krieger,
Nucl.
Instrum.
Methods
B
266
(2008)
1979. [9]
C.
Borschel,
C.
Ronning,
Nucl.
Instrum.
Methods
B
269
(2011)
2133. [10] W.
Möller,
Nucl.
Instrum.
Methods
B
322
(2014)
23. [11] R.
Timilsina,
S.
Tan,
R.
Livengood,
P.D.
Rack,
Nanotechnology
25
(2014)
485704. [12] G.
Hobler,
D.
Kovac,
Nucl.
Instrum.
Methods
B
269
(2011)
1609. [13] G.
Hobler,
Nucl.
Instrum.
Methods
B
352
(2015)
22. [14] Y.G.
Li,
Y.
Yang,
M.P.
Short,
Z.J.
Ding,
Z.
Zeng,
J.
Li,
Sci.
Rep.
5
(2015)
18130. [15] T.
Kenmotsu,
Y.
Yamamura,
T.
Muramoto,
N.
Hirotani,
Nucl.
Instrum.
Methods
B
228
(2005)
369. [16] A.
Hu,
A.
Hassanein,
Nucl.
Instrum.
Methods
B
281
(2012)
15. [17] S.J. Lindsey, PhD Thesis, TU Wien, Vienna, Austria, 2015. [18] S.
Plimpton,
J.
Comput.
Phys.
117
(1995)
1. [19] J.
Tersoff,
Phys.
Rev.
B
39
(1989)
5566. [20] J.F.
Ziegler,
J.P.
Biersack,
U.
Littmark,
The
Stopping
and
Range
of
Ions
in
Solids,
Pergamon
Press,
New
York,
1985. [21] G.
Hobler,
Nucl.
Instrum.
Methods
B
96
(1995)
155. [22] C.
Ebm,
G.
Hobler,
Nucl.
Instrum.
Methods
B
267
(2009)
2987. [23] O.S.
Oen,
M.T.
Robinson,
Nucl.
Instrum.
Methods
132
(1976)
647. [24] Y.
Kudriavtsev,
A.
Villegas,
A.
Godines,
R.
Asomoza,
Appl.
Surf.
Sci.
239
(2005)
273. [25] P.
Sigmund,
Phys.
Rev.
184
(1969)
383. [26] H.-B.
Kim,
G.
Hobler,
A.
Lugstein,
E.
Bertagnolli,
J.
Micromech.
Microeng.
17
(2007)
1178. [27] H.-B.
Kim,
G.
Hobler,
A.
Steiger,
A.
Lugstein,
E.
Bertagnolli,
Nanotechnology
18
(2007)
245303. BAmerican Societyfor MassSpectrometry, 2018 J.Am.Soc. MassSpectrom. (2018) DOI:10.1007/s13361-018-2102-z RESEARCH ARTICLE C-OBondDissociation andInducedChemicalIonization UsingHighEnergy(CO2)n+ Gas Cluster Ion Beam HuaTian,1 DawidMaciążek,2 Zbigniew Postawa,2 Barbara J.Garrison,1 Nicholas Winograd1 1ChemistryDepartment, ThePennsylvaniaState University, 215ChemistryBuilding, UniversityPark, PA 16802, USA 2SmoluchowskiInstitute of Physics, Jagiellonian University, ulicaLojasiewicza 11, 30-348, Krakow,Poland Abstract. A gas cluster ion beam (GCIB) source,consisting of CO2clustersandoperating with kineticenergiesofup to60keV, has been developed for the high resolution and highsensitivityimagingofintactbiomolecules. TheCO2moleculeis an excellent moleculeto employ in a GCIB source due to its relative stability and improved focusing capabilities, especiallywhen comparedto the conventionallyemployedArclustersource.Here we re port on experimentsaimedto examinethebehaviorofCO2clusters astheyimpactasurfaceunderavarietyof conditions. Clusters of(CO2)n+(n =2000~10,000) with varying sizes andkineticenergies were employedto interrogateboth anorganic andinorganicsurface.TheresultsshowthatC-Obonddissociationdid notoccur whenthe energypermoleculeislessthan5 eV/n,butthat oxygenadductswere seeninincreasingintensity as theenergyisabove5eV/n,particularly,drastic enhancement upto100timesofoxygenadductswasobserved onAusurface.ForIrganox1010, an organicsurface,oxygencontaining adducts wereobservedwith moderate signalenhancement.MoleculardynamicscomputersimulationswereemployedtotestthehypothesisthattheCObondisbrokenathighvaluesof eV/n.ThesecalculationsshowthatC-Obonddissociation occursat eV/n valueslessthantheC-Obond energy(8.3eV)byinteractionwith surfacetopologicalfeatures.Ingeneral,the experiments suggest that the projectiles containing oxygen can enhance the ionization efficiency of surface molecules via chemically induced processes, and that CO2 can be an effective cluster ion source for SIMS experiments. Keywords: Gas cluster ion beam, C-O bond dissociation, Carbon dioxide cluster, Irganox 1010, Au film, Moleculardynamics computersimulations Received: 13 September2018/Revised:31October2018/Accepted:1November2018 Introduction G as cluster ion beam (GCIB)sources have gained popularity in conjunction with secondaryion mass spectrometry (SIMS) experiments.These beams provide an effective means of eroding an organic material without damage accumulation, Electronic supplementary material Theonline versionofthisarticle(https:// doi.org/10.1007/s13361-018-2102-z)contains supplementary material, which is availableto authorizedusers. Correspondenceto: HuaTian; e-mail: hut3@psu.edu resulting in improved molecular depth profiling behavior. Moreover,whenemployedasaSIMSsource,theGCIByields mass spectra that exhibit improved detection of larger intact molecular ions and reduced formation of molecular fragment ions relative to atomic ion and C60 clusters conventionally employed[1–4]. The most commonlyemployed GCIB consists of clusters composed of Ar, generallyin the size range of n= 100 to 10,000 atoms. At a typical acceleration energy of 20 keV, each Ar atom carries between2 and 200 eV/n. This eV/n value spans the range of energies associated with chemical bonds andis the primary reason that less fragmentation is observedduringdesorption,particularlyforthelarger clusters. H.Tian etal.:C-OBondDissociationInducedbyHighEnergyCO2-GCIB AnadvantageoftheGCIBdesignisthatthe compositionof the cluster can be varied at will. For example, it has been possible to tune the chemistryof theGCIBto enhance ionization probability, and hence, expand the range of applications amenable to study. We have employed the conceptof a mixed cluster whereby a hydrogen-containing small molecule is doped into the Ar cluster to create an environment suitable for protonation of desorbed molecules[5]. This approach has been particularly successful when using HCl as a dopant, especially when H2O is present either in the vacuum or on the sample surface. The waterpresumablyenhances theformationofH3O+ andtheformationof[M+H]+,whereMisthe molecule of interest. In a related effort, large water clusters have been generateddirectly, andhave demonstrated ion yield enhancements of greater thana factor of 10[6–8]. Moreover, the HCl-dopedAr cluster beam can overcome salt suppression by promotingthe protonated molecule ion, mitigating the pronounced matrix effect in biological samples[7]. In an effort to search for the ideal gas candidate for the cluster, we explored thepossibilityof employingclustersof CO2 as theGCIB.In earlierwork, we demonstrated that this cluster is more stable than Ar clusters and hence is produced usinglower pressuresintheexpansionregionofthe GCIB source [9].Moreover,thisstability resultsinbetter focusing properties than Ar clusters. For example, using a 20-keV beam,Arclusterscouldbe focusedtoa20-μmspot, whileasimilarCO2beam couldbefocusedtoa7-μmspot. This resultis obviouslycriticalformolecule-specificimaging experiments where it has been a goal to reach submicron spatial resolution, a value important for imaging of biologicalcells. Andfinally, bothmolecular dynamicscomputer simulations and experiments suggested that the behaviorofCO2wasvery closetothatofArwhen energyper molecule<5 eV/n, where theCO2 molecule actsasa pseudo-atom of mass 44 amu[9]. To improve the spatial resolutionof GCIB sources forbioimaging applications,an ion gun has been constructed with an acceleration voltage of 60 keV, allowing submicron spot sizes to be achieved. However, this expanded energy range opens the possibility that CO2 molecules will be dissociated during ion impact. Akizukietal.havereportedthatSiO2film formsonsilicon surface upon the CO2 cluster beam irradiation with minimumenergy per molecule of40 eV/n,resulting from the high densityenergydepositionwithinlocalized surfacearea [10]. Here we show from experiment and computer simulations that for energy per molecule above 5 eV/n dissociation is observed as evidenced by ejected molecular ions thathaveincorporated oxygen atoms.Wepresentasystematic study over the kinetic energy range of 20 to 60 keV, and for cluster sizes of 2000 to 10,000 CO2 molecules to examine the magnitude of the effect. As a consequence of this newchemistry,theintensityofmolecularions containing oxygen can be enhanced by more than a factor of 10 with organic molecules and more drastically on metal surfaces. In general, we show that this ion source offers new opportunities forGCIB SIMS. Experimental Sample Preparation The Irganox 1010 thin film was purchasedfrom the National Physical Laboratory (NPL, Huntington, UK). The film was depositedontofinelypolished silicon wafers (1.0×1.0× 0.5 mm3). The thickness of the film was monitored during thedepositionprocess atNPLbyaquartz crystalmicrobalance calibratedbyspectroscopicellipsometry.Specifically,thesample measured in this workhad a thickness of 49.5 nm. Agold film was purchasedfromElectron MicroscopySciences(Hatfield,USA)withathickness of 18 μm. Thefilmwas pre-sputtered using a 20-keV Ar1000 + cluster beam with a current of 100pA for 30 min to remove the organic contaminants on the surface. SIMS Characterization Theinitial motivationforthis workstemsfrom anobservation of intense gold oxide ions from a Au surface under bombardmentby(CO2)n+.Duetothelimited supplyof oxygeninthe vacuum chamber where the sample is analyzed, it was speculatedthatOisproducedbydissociationofCO2molecules.This study is aimed toward elucidating the conditions where C-O bonds break andif the oxidization occurs on organic surfaces as well. The gas cluster ion beam utilized here (GCIB SM from Ionoptika,UK)isshownschematicallyinFigure 1.Thesystem is designed to operate at kinetic energies of up to ~ 60 keV, Lens 1 Lens 2 Expansion chamber Figure 1. Schematic of GCIM SM. The GCIB SM is a gas cluster ion beam source developed by Ionoptika (Chandler’s Ford,UK),for molecularimagingSIMS.Itoperatesat a maximumaccelerationpotentialof70kV,usingatwolenssystemto form a spot size of 1 μm. Running CO2 gas, the system is capable of operating with cluster sizes of <30,000. Cluster generation occursashighpressureinputgaspassesthrough a de Laval nozzle, condensing into clusters at the exit due to rapidexpansionandcoolingofthegas.Theneutralgasbeamis then ionizedbyelectronbombardment, before beingfocused throughthecolumn.AWienfilterallowsformassfiltering,anda narrowercluster sizedistribution. Timeofflightmeasurements enableprecise tuningoftheclusterssize H. Tian etal.: C-OBond Dissociation Inducedby High Energy CO2-GCIB withaclustersizeofupto30,000 molecules.TheGCIBSM was installed on a buncher-TOF SIMS instrument, J105 3D ChemicalImager(Ionoptika,UK)[11]andisincidentuponthe sample with an angle of 45°. The spectra from fresh Irganox and Au surfaces were acquired using (CO2)n+ with varying cluster sizes and energies. The kinetic energy of the clusters has a range of 2.5~30 eV/n, determined by the availability of the cluster from the GCIB SM. Specifically, (CO2)n+ clusters consisting of 2000, 4000, 6000, 8000, and 10,000 molecules with an energy range of 20~60 keV were used. The cluster size was selected using a Wienfilter.Aspreadof ~±5 μs(FWHM)ToF wasmeasured, which,forinstance, correspondstoaspread±370CO2 moleculesfor50keV(CO2)2000 + and±830for50keV(CO2)10000 + cluster projectiles. On the Irganox film, the spectra were acquired over area of 50×50 μm2 using10pAof eachclusterionbeamwith32×32 pixels and100shots per pixel. The ion dose for each spectrum is 2.5×1013 ions·cm−2. Five parallel measurements on a fresh Irganox 1010 film surface were conducted. On the gold sur-face, 10 parallel measurements were conducted on the same spot,butonlythelastfivespectrawereusedfordataprocessing to avoidsurface contaminationbyorganicadventitious sources of organic molecules. Control experiments were performed using 50 keV (Ar)n+ withaclustersizeof2000,4000,6000,8000, and10,000with the same condition as corresponding(CO2)n+ clusters on both the Irganox 1010 andtheAu film. DataProcessing The software Ionoptika Image Analyser (Version: 1.0.8.14) was usedto readtheintensityof selectedions,[M+2O]− at m/z1207.77(Δm/z 0.6) for Irganox 1010 and m/z 228.96(Δm/z 0.2)forAu.The averageintensityandthestandarddeviationof selected ions were plotted against energy per molecule of the cluster ions beams. Computer Simulation A detailed description of the molecular dynamics computer simulations used to modelcluster bombardment can be found elsewhere[12]. Briefly, the motion of the particles is determined by integrating Hamilton’s equations of motion. The forces among the particles are described by a blend of pairwise additive and many-bodypotential energy functions. The interaction amongcluster projectile molecules consisting ofC and O atoms is described by the ReaxFF-lg [13] potential splined with a ZBL potential[14]toproperlydescribehighenergy collisions. For interaction amonggold atoms, the EAM potential[15]is usedwhiletheinteractionbetweengold atoms and cluster projectile atoms is described by a ZBL potential. ThecalculationsareperformedwithaLAMMPScode[16]that was modified for a more efficient modeling of sputtering phenomena. The experiments were performed on an organic substrate (Irganox) and a metal (gold) with the measured quantities being attachment of O atoms to the ejected molecule or atom.The more straightforwardof thesamplestomodelis AubecausetheonlysourceofOintheejectedspeciesisthe incident cluster beam. Thus, the dissociation of the CO2 molecules that occurs early in the collision cascade can be used as a signal that Au-O adducts can be formed. Two types of Au samples were studied, a flat Au(100) surface and an Au sample with already developed morphology whichwas createdinaprevious study[17].In the caseof flat Au(100), a hemisphere with a diameter of ~ 40 nm, consistingof ~1 millionatoms, was cut outof perfectgold crystal. For the roughened surface, a cylinder with a diameter of ~ 60 nm and height of 25 nm, consisting of ~ 2.4 millionatoms,waschosen.In both cases,rigidandstochasticregions withathicknessof0.7and2.0nm, respectively, were used around the sample. These boundary conditions simulate the thermal bath that keeps the sample at the required temperature and helps inhibit the pressure wave reflectionfromthe system boundaries[18]. Thesesamples werebombardedby50keV(CO)n clusterprojectiles consisting of 4000, 6000, 8000, and 10,000 molecules at 45°impact angle to reproduce experimental conditions. For the roughened surface, two types of impacts were chosen, one that aimed at a hill or protrusion into the vacuum and one that aimed at a depression or hole position. The simulations are calculated for a 0-K target temperature. Each impact eventis evaluatedto8ps, whichissufficientlylong to determinethe numberofCO2 moleculesthatdissociate. Results andDiscussions Investigation of Oxygen Adducts from Organic andInorganic Surfaces UnderBombardment of(CO2)n+ Beams Duringbombardmentbythe(CO2)n+clusterionbeam, oxygen adducts [M+nO]− were observed to be emitted from the sur-face of Irganox 1010 andAu, as shown inFigure S1.Several oxide ions are seen from Irganox 1010, described as [MIrganox H+nO]− n =1~4. Correspondingly, m/z values at m/z = 1191.77, 1207.77, 1223.76, and 1239.76 are observed. The intensity of the oxygen adducts decreases slightly with the additionofOatoms.OntheAu surface,[MAu+nO]− n=1~3 are the most intense ions, with m/z values of 212.96, 228.96, and 244.95. Meanwhile, oxygen adducts of gold clusters are observed with composition of [2MAu+nO]− and [3MAu+ nO]− n= 1~3, with intensities that are lower when compared to [MAu +nO]−. Only the bi-oxygen adducts, [MIrganox-H+ 2O]− and[MAu +2O]− are discussed as they are representative of the other oxide products. The StudyontheFormationofOxygenAdducts atVaryingKinetic Energyof(CO2)n+ Beams To understand the formation of oxygen adducts, a systematic study was performed using (CO2)n+ beams with H.Tian etal.:C-OBondDissociationInducedbyHighEnergyCO2-GCIB varying cluster size and energy. The ion intensities of [MIrganox-H+2O]− and[MAu +2O]− are plotted versus the kinetic energy per molecule, eV/n, in Figure 2a, c, respectively. It is clear that oxygen adducts occur beyond a threshold value ~5 eV/n.Above thethreshold, theintensity of the oxygen adducts increases along with the increase of the energy per molecule. In general, for a given energy per molecule value, the beam with the higher kinetic energy yields higher oxygen adduct intensities than the beamwithlowerkineticenergyatthe same energyper molecule. For example, at energyper molecule of 10 eV/n, theintensity of [MIrganox-H+2O]− from Irganox 1010 is 30% higher using 60 keV (CO2)6000 + than using 40 keV (CO2)4000 +, and 63% higher compared with 20 keV (CO2)2000 +as shownin Figure 2.Asimilartrendof[MAu+ 2O]− from the Au surface is observed in term of the intensity changes of oxygen adducts with varying eV/n of the analyzing beams, except that there is more dramatic enhancement of oxygen adducts at higher kinetic energy and eV/n. These results certainly suggest that oxygen originates from the (CO2)n+ cluster itself. Above 3 eV/n according to our observations, the C-O bonds dissociate during impact on the sample surface. The energy carried by each CO2 molecule allows the free O atoms to interact with the sample surface instead of escaping, yielding the various oxygen adducts. The ratio of C-O bond Intensity of m/z 228.96 (Counts) Intensity of m/z 1207.77 (Counts) 5 2x10 0 7 2x10 0 dissociation increases along with the increasing kinetic energy and energy per molecule of the beams, resulting in the higher intensity of oxygen adducts. The very low intensity of oxygen adducts on the Au surface using 20keV(CO2)n+ beams could result fromthe lowersputter rate of metal atoms. At higher kinetic energy, 30, 40, and 50 keV, the sputtering yield increases and more free O atoms promote the formation of Au oxide. This concerted actionleads to the enhancementof the Au oxygen adducts at higher kinetic energy. To rule out the possibility of other sources of free O atoms, the control measurements were carried out using (Ar)n+ beams at the same cluster size. Here, only the measurements using 50 keV (Ar)n+ beams are plotted as in Figure 2b, d to simplify the presentation. It is clear that only16~30% of oxygen adducts, [MIrganox-H+2O]−,are observed on the Irganox 1010 surface using 50 keV Arn+ compared with corresponding (CO2)n+ beams when energy per molecule is over 8 eV/n as shown in Figure 2b. The source of the O atoms for these oxide products could be fragmented from the Irganox 1010 molecule during bombardment. While on the Au surface, the oxygen adducts, + [MAu +2O]−, are barely detectable using 50 keV Arn beamsasshown in Figure 2d. The oxygen adducts of Au are drastically increased by up to 850-fold using 50 keV (CO2)n+ beams. It is evident that (CO2)n+ beams with Figure 2. Ionintensitiesasafunctionofkineticenergypermolecule,eV/nfordifferenttotalenergiesofCO2orArcluster.(a)Signal of m/z 1207.77, assigned as [MIrganox-H+2O]− from Irganox 1010 using 20–60 keV (CO2)n+ cluster at energy per molecule of 2.5~30 eV/n.(b)Comparison ofsignalofm/z 1207.77fromIrganox1010using50keV(CO2)n+ andArn+atenergypermoleculeof5– 25 eV/n.(c)Signalofm/z 228.96, assigned as [MAu +2O]− fromAufilm using20–50keV(CO2)n+clusterat energyper moleculeof 2.5~25 eV/n.(d)Comparisonofsignalofm/z228.96fromAufilmusing50keV(CO2)n+ andArn+atenergypermoleculeof5–25eV/n H. Tian etal.: C-OBond Dissociation Inducedby High Energy CO2-GCIB 50 40 30 20 10 0 E (eV/n) Figure3. Percentage of dissociatedCO2moleculesfor50 keV (CO2)n+impact ondifferentconfigurationsofAu surface higher energy per molecule favor the oxygen adducts, leading to the enhanced chemical ionization by providing an oxygen atom source through the dissociation of C-O bonds in the cluster beam itself. ComputerSimulationofC-OBondsDissociation atEnergyper Molecule below5eV/N The percentages of CO2 molecules that dissociate in the simulations as a function of energy per molecule are shown in Figure 3 for impacts on the flat Au surface, the hill position, and the hole position. Examining first the flat surface results (black line andpoints), there is no dissociation below eV/n~ 8eV, consistentwiththeCObond energyinCO2of ~8.3eV. The resultsforthe roughened surface(red andbluelines and points)indicatethattheCO2 moleculedissociatesatincident energies as low as 5 eV/n. The side views of the 50 keV (CO2)10000 + cluster-surface impact event at the flat and roughened surface are shown in Figuer 4 for three times at the beginning of the impact event for the flat surface and a hill configuration.For these impact events, the initialkinetic energyper moleculeis5eV/n.Smallballsdepictintact molecules, while large balls represent molecules that have undergone dissociation. The impact at the bottom of the valley is not shown but leads to the same conclusions. The color of the Percentage of dissociationof C-O bond (%) species is based on the current center of mass kinetic energy wherewhiteisthe original valueof5eV,blue moleculeshave less energy, and red molecules have more kinetic energy. The initial sequence of events is similar for the flat and roughened surfaces. The cluster projectile flattens during the impact, whichleads to a buildupof pressure inside the projectile. The pressure is subsequently relieved by material ejection in a lateral direction. As indicated by the red color, many molecules during the decompression phase are accelerated to a kinetic energy higher than the initial value. The fate of accelerated molecules is, however, different on the flat and roughened surfaces. On the flat surface, most of these molecules slide alongthe surface and are reflectedintothe vacuum intact.Forthe roughened surface,however, someof accelerated (red) molecules collide with the nearby irregular surface. Now the kinetic energy per molecule can be as high as 10 eV, sufficient to cause dissociation on the second collision. We wouldexpectthesameprocesstooccuronaroughenedorganic substrate such as Irganox. The effect of the buildup in pressure and subsequent increase in the kinetic energy of a molecule during the impact eventshouldbegreaterthehigherthetotalincident energyfora given initial energy per molecule value.This effect is seen in theexperimentaldatashowninFigure2c.There areno[MAu+ 2O]− ions observed for a total energy of 20 keVup to energy per molecule of 10 eV/n. On the other hand, for a total energy of50keV, the[MAu+2O]− ions are observed at less than energy permolecule of5eV/n. Conclusion andOutlook Asystematic studywas performedto investigatethe dissociation of C-Obonds in a CO2 GCIB upon impact with a sample surface, consistingof an Irganox 1010 thin film and a Au film. Theclusterswith varyingeV/natdifferentkinetic energywere employedtointerrogatethesample surface.Thedissociationof C-O bonds and the observation of oxygen adducts occurs above ~5 eV/n for both samples. The moderate increase ~2 times of oxidation was observed on organic Irganox 1010 film along with the energy per molecule and kinetic energy of the Figure4. Timelapseimagesofdistributionofenergypermoleculein50keV(CO2)10,000 clusterduringtheimpact ontheflat(top panels)and roughened(bottompanels)Au(100)surface.ThesmallballsrepresentintactCO2 molecules.Thelargerballsrepresent dissociated molecules.The colorscaleisthe centerof masskineticenergyof each particle.Theinitialkineticenergyis energyper molecule of5eV/n,denotedbythe white color H.Tian etal.:C-OBondDissociationInducedbyHighEnergyCO2-GCIB beam, and the enhancement tends to reach a plateau when energy per molecule reaches 25 eV/n. The same trend was observed on Au films but with more drastic enhancement up to14-fold comparedthebeamsat5 and20eV/nwithkinetic energy of 60 keV; moreover, the trend points to further increases if the beam energy couldbe increasedbeyond 60 keV. Theunexpectedly observeddissociationof(CO2)n+ below 8eV/n is shown via computer simulation to arise from surface topography.The resultsdemonstratethat ona rough surface, CO2 molecules may be compressed during interaction with surface features, yielding higher kinetic energies during the decompression phase. Our findings provide new insight into the cluster-sample interaction, revealing a new chemical ionization pathway that couldlead to enhanced sensitivity. Moreover, operating in a high kinetic energy regime opens the possibilityof findingadditional gas candidates that may more effectively exploit the chemical ionization mechanism. Acknowledgments This work was supportedby NIHgrants 5R01GM113746-21. DM andZP gratefully acknowledge financialsupportfrom the Polish National Science Centre, Grant No. 2015/19/B/ST4/ 01892.MD simulationswere performed at thePLGrid Infrastructure. References 1. Tian,H.,Sparvero,L.J.,Amoscato,A.A.,Bloom,A.,Bayr,H.,Kagan, V.E., et al.: Gas cluster ion beam time-of-flight secondary ion mass spectrometry high-resolution imaging of cardiolipin speciation in the brain: identification of molecular losses after traumatic injury. Anal. Chem. 89,4611–4619(2017) 2. Rabbani,S.,Barber,A.M.,Fletcher,J.S.,Lockyer,N.P.,Vickerman,J.C.: TOF-SIMS withargon gas clusterionbeams: acomparisonwith C60+. Anal.Chem.83, 3793–3800(2011) 3. Ninomiya, S., Nakata, Y., Ichiki, K., Seki, T., Aoki, T., Matsuo, J.: Measurements of secondary ions emitted from organic compounds bombarded with large gas cluster ions. Nucl. Instrum. Methods Phys. Res. B. 256,493–496(2007) 4. Ninomiya, S., Ichiki, K., Yamada, H., Nakata, Y., Seki, T., Aoki, T., etal.:Preciseandfastsecondaryion massspectrometrydepthprofilingof polymermaterialswithlargeAr clusterionbeams. RapidCommun. Mass Spectrom. 23, 1601–1606(2009) 5. Wucher, A., Tian, H., Winograd, N.: A mixed cluster ion beam to enhancetheionizationefficiencyin molecularsecondaryionmassspectrometry.RapidCommun. MassSpectrom. 28,396–400 (2014) 6. Tian, H., Wucher, A., Winograd, N.: Reducing the matrix effect in organic cluster SIMS using dynamic reactive ionization. J. Am. Soc. Mass Spectrom.27, 2014–2024(2016) 7. Tian,H.,Wucher,A.,Winograd,N.:Reducethematrixeffectinbiologicaltissue imaging usingdynamicreactiveionizationandgascluster ion beams. Biointerphases.11, 02A320(2016) 8. Tian,H.,Wucher,A.,Winograd,N.:Dynamic reactiveionizationwith clustersecondaryionmassspectrometry.J.Am. Soc. MassSpectrom.27, 285–292(2016) 9. Tian,H.,Maciazek,D.,Postawa, Z.,Garrison,B.J., Winograd,N.: CO2 clusterionbeam, analternativeprojectileforsecondaryionmassspectrometry.J.Am. Soc.Mass Spectrom. 27,1476–1482(2016) 10. Akizuki,M.,Matsuo,J.,Yamada,I.,Harada,M.,Ogasawara,S.,Doi,A.: SiO2filmformationat roomtemperaturebygasclusterionbeam oxida-tion.Nucl.Instr.Meth. Phys. Res. Section B.112,83–85 (1996) 11. Fletcher,J.S.,Rabbani, S.,Henderson,A.,Blenkinsopp, P., Thompson, S.P.,Lockyer,N.P., etal.:A newdynamicinmassspectralimaging of singlebiological cells.Anal.Chem. 80, 9058–9064(2008) 12. Garrison,B.J.,Postawa,Z.:Computationalviewofsurfacebasedorganic mass spectrometry.MassSpectrom.Rev. 27,289–315 (2008) 13. Liu, L.C., Liu, Y., Zybin, S.V., Sun, H., Goddard, W.A.: ReaxFF-lg: correctionoftheReaxFFreactiveforcefieldforLondondispersion, with applications to the equations of state for energetic materials. J. Phys. Chem.A. 115, 11016–11022(2011) 14. Biersack,J.P.:Theeffectofhighchargestatesonthestoppingandranges ofionsin solids.Nucl.Instrum.Meth.B. 80-1,12–15(1993) 15. Williams,P.L.,Mishin,Y.,Hamilton,J.C.:Anembedded-atompotential for the Cu-Ag system. Model. Simul. Mater. Sci. Eng. 14, 817–833 (2006) 16. Plimpton, S.: Fast parallel algorithms for short-range molecular-dynamics. J.Comput. Phys.117,1–19 (1995) 17. Maciazek, D., Paruch, R.J., Postawa, Z., Garrison, B.J.: Micro-and macroscopic modeling of sputter depth profiling. J. Phys. Chem. C. 120, 25473–25480(2016) 18. Postawa,Z.,Czerwinski,B.,Szewczyk,M.,Smiley,E.J.,Winograd,N., Garrison, B.J.: Enhancement of sputteringyields duetoC60versusGa bombardmentofAg{111}as exploredby moleculardynamicssimulations.Anal. Chem. 75,4402–4407(2003) PHYSICAL REVIEW B 97, 155307 (2018) Crater function moments: Role of implanted noble gas atoms Gerhard Hobler,1,*
Dawid Maciążek,2
and Zbigniew Postawa2
1Institute of Solid-State Electronics, TU Wien, Floragasse 7, A-1040 Wien, Austria 2Institute of Physics, Jagiellonian University, ul. Lojasiewicza 11, 30348 Krak, Poland (Received 12 December 2017; revised manuscript received 5 April 2018; published 24 April 2018) Spontaneous pattern formation by energetic ion beams is usually explained in terms of surface-curvature dependent sputtering and atom redistribution in the target. Recently, the effect of ion implantation on surface stability has been studied for nonvolatile ion species, but for the case of noble gas ion beams it has always been assumed that the implanted atoms can be neglected. In this work, we show by molecular dynamics (MD) and Monte Carlo (MC) simulations that this assumption is not valid in a wide range of implant conditions. Sequential-impact MD simulations are performed for 1-keV Ar, 2-keV Kr, and 2-keV Xe bombardments of Si, starting with a pure single-crystalline Si target and running impacts until sputtering equilibrium has been reached. The simulations demonstrate the importance of the implanted ions for crater-function estimates. The atomic volumes of Ar, Kr, and Xe in Si are found to be a factor of two larger than in the solid state. To extend the study to a wider range of energies, MC simulations are performed. We find that the role of the implanted ions increases with the ion energy although the increase is attenuated for the heavier ions. The analysis uses the crater function formalism specialized to the case of sputtering equilibrium. DOI: 10.1103/PhysRevB.97.155307
I. INTRODUCTION The
crater
function
formalism
[1–3] has recently established itself as an important technique in the theory of ion-beam induced spontaneous pattern formation. Earlier approaches used Sigmund’s theory of sputtering
[4]
together with a surface smoothing mechanism to derive equations of motion for
the
surface
height
function
[5–8], some including simplified models of atom redistribution due to momentum transfer to target
atoms
[9–13]. Other work explains pattern formation in
terms
of
mechanical
stress
[14–18]. In contrast, in the crater function formalism, the coefficients of the equation of motion are calculated in terms of the moments of the so-called crater function. The crater function is defined as the average response
of
the
surface
to
a
single
ion
impact
[19]
and
may
be
obtained
by
molecular
dynamics
(MD)
[19–22] or binary collision
Monte
Carlo
(MC)
simulations
[23–26]. In this way, the crater function formalism allows to include the results of numerical simulations in the theory of spontaneous pattern formation. In most previous work, using the crater function formalism, only sputtering and atom redistribution have been taken into account. The effect of the implanted ions on pattern formation has
been
investigated
only
recently
[27,28]. It has been found that ion implantation has a destabilizing effect on the surface along the projected beam direction, if the incidence angle exceeds a critical value, while it always has a stabilizing effect in transverse direction. The impact of the implant contributions is expected to be particularly significant when erosion and redistribution are moderate such as for energetic light ions, and
when
the
surface
binding
energy
is
high
[28].
*gerhard.hobler@tuwien.ac.at It is commonly assumed that the effect of the implanted ions on pattern formation is negligible, if the implanted atom species is a noble gas (NG). This is motivated by the fact that NG
atoms
have
a
tendency
to
leave
the
target
[29],
so
the
retained fluence is lower than for nonvolatile ions. Radiation induced transport of NG atoms towards the surface is well established
experimentally
[30]
and
has
been
named
“rapid
relocation”
[31,32] in the absence of a clear understanding of the mechanism on an atomistic level. MC simulations cannot explain
rapid
relocation
[33],
thus
it
seems
to
be
a
nonballistic
effect. The rationale for neglecting implanted NG atoms in pat-tern formation, however, is flawed. First, the effect strongly depends on the system studied. For instance, while after 1-keV Ar bombardment of SiO2 virtually no Ar is present in
the
oxide
[32],
for
5-keV
Ne
bombardment
of
ta-C
[28] the experimentally found areal density of Ne agrees with that predicted by MC simulations. Second, the NG ions do not disappear instantaneously when they come to rest, but require additional ion bombardment to be transported (“rapidly relocated”) towards the surface. It will be shown that the separation of implantation and removal of the NG atoms leads to a contribution to the crater function. The purpose of this work is to demonstrate that the implanted NG ions indeed have an important effect on the crater function moments. As a target we choose Si, which is a commonly studied material. Selected low-energy ion bombardments are investigated by MD simulations. MC simulations are used to study the effect in a wider range of impact energies. In the absence of a model for rapid relocation, MC simulations cannot predict the contribution of NG atom redistribution within the target. However, we will show that the net effect of the volumes added by the implanted ions and removed by their erosion is well predicted by MC simulations. 2469-9950/2018/97(15)/155307(13) 155307-1 ©2018 American Physical Society HOBLER, MACIĄŻEK, AND POSTAWA PHYSICAL REVIEW B 97, 155307 (2018) Our analysis will be based on the assumption of sputtering equilibrium. We therefore present the crater function formalism specialized to sputtering equilibrium in Sec. II.
Details
of the MD and MC methods employed are given in Sec. III.
In Sec. IV,
the
MD
and
MC
results
are
presented.
Finally,
in
Sec. V
we discuss the limitations of our approach and argue that they do not affect our qualitative conclusions. II. CRATER FUNCTION FORMALISM We start with a summary of the crater function formalism in order to define the nomenclature (Secs. II
A
and II
B).
Assumptions that are either normally made or are introduced by us, will be discussed in Sec. II
C.
Their
consequences
for
the
calculation of the crater function moments will be presented in Secs. II
D
and II
E.
A. The crater function In the crater function formalism, the equation of motion of the surface is written in terms of the moments of the crater function F(x,y), where x is parallel to the projection of the beam direction to the originally flat surface, y is the surface coordinate perpendicular to x, and x =y =0 corresponds to the impact point. The crater function is the average change in height due to the impact of a single ion. Different signs have been used for the crater function in previous work. When only erosion is considered, it is natural to define F(x,y)asthe negative
change
in
surface
height
[2,28]so F(x,y) is always positive. Here, we consider implantation and erosion on an equal footing. Therefore we prefer to adhere to the original convention introduced by Kalyanasundram et al. [19],
which
has also been used by Norris et al. [1,3], where erosion leads to negative and implantation to positive contributions to the crater function. B. Crater function moments and surface stability The ith moment with respect to x is defined by  ∞ (i) i M=xF(x,y) dxdy. (1) −∞ The crater function is composed of contributions from different atom species X and different mechanisms (erosion, redistribution, and implantation), so we can write  F(x,y) = FX(x,y) (2) X with FX(x,y) =FX,eros(x,y) +FX,redist(x,y) +FX,impl(x,y). (3) The same partitioning applies to the moments:  (i)(i) M=M(4) X X with (i)(i)(i)(i) M=M+M(5) X X,eros X,redist +MX,impl. Some of the contributions vanish by definition; for the zeroth moment, the contribution of atom redistribution is zero if one assumes that the atomic volume does not change during relocation: (0) M=0. (6) X,redist Moreover, the target atom species X =NG, of course, have no implant contributions: (i) M=0. (7) X=NG,impl In the work, introducing the extended crater function formalism
[2],
the
moments
defined
in
Eq.
(1) were written M(i) to distinguish them from moments with respect to y or curvature. However, in the final result, neither the moments with respect to y nor the moments with respect to curvature of order higher than zero appear. Zeroth-order moments are physically the same, independent of which quantity is considered variable. We have therefore dropped the index “x” in order to avoid confusion with X, which we use to denote atom species. According to Harrison and Bradley, surface stability is determined by the curvature coefficients C11 and C22 given by
[2]
C11 =∂ ∂θ (M(1) cos θ)    +∂ ∂K11 (M(0) cos θ)    (8) S11 T11 and C22 =cot θ (M(1) cos θ)    +∂ ∂K22 (M(0) cos θ) , (9) S22   T22  where θ denotes the incidence angle with respect to the surface normal and K11, K22 the surface curvatures in x and y direction, respectively. The derivatives with respect to the curvatures are evaluated
at
zero
curvature.
Compared
to
Ref.
[2],
all
terms
appear with the opposite sign here due to the different sign conventions for the crater function. However, the meaning of C11 and C22 is
the
same
as
in
Ref.
[2];
instabilities
occur
when C11 and/or C22 is negative. Apart from the signs, the magnitudes of C11 and C22 are also of interest since the growth rate and velocity of the patterns in the linear regime depend on
them
[5].
As
will
be
shown,
M(0) is largely independent of the implanted ions in the cases studied, so the effect of the implanted ions on the stability of the surface is via S11 and S22 only
[compare
Eqs.
(8)
and
(9)]. C. Assumptions We will use the following assumptions as appropriate. (1) Each atom species X has a fixed atomic volume X that is known or can be determined and does not change during ion bombardment. (2) Sputtering equilibrium has been established, which means that the average rates of addition and removal of the ion species are equal. (3) The implanted ions do not affect the spatial characteristics of the collision cascades. (4) The spatial distribution of NG ejection with respect to the ion’s impact point is unaffected by rapid relocation. Assumption 1 has first been made by Norris et al. [1]to
arrive at a simplified algorithm for determining the moments from MD simulations (see Sec. IID).
This
algorithm
is
also
the only known way to determine the moments from MC 155307-2 CRATER FUNCTION MOMENTS: ROLE OF IMPLANTED … PHYSICAL REVIEW B 97, 155307 (2018) simulations
[23].
It
is,
however,
not
self-evident,
in
particular
for NG atoms: when NG atoms are implanted, they are likely to initially end up in an interstitial position, while they usually
form
bubbles
upon
further
bombardment
[34].
The
NG atoms do not necessarily have the same volume in these configurations, and the atomic volume is also not known a priori. Section IVB
will address this in some detail. Assumption 2 seems safe: sputtering equilibrium is reached for
sputter
depths
on
the
order
of
the
ions’
projected
range
[35],
while pattern formation usually requires higher fluences. Assumption 3 can be valid only approximately: the different mass and volume of the implanted ions changes the stopping power of the target and therefore the spatial characteristics of the collision cascades. The effect is weaker the smaller the mass difference and the ion concentration. For the cases studied in this work (Ar, Kr, Xe in Si) we find from our MD simulations that assumption 3 is well fulfilled. Beyond these cases, we expect the effect of the implanted ions to increase with increasing implant energy, at lower implant angles, and when the rapid relocation effect is weaker, as all these conditions increase the implanted ion concentrations. The rationale for assumption 4 is the following: the implanted NG atoms do not disappear immediately; otherwise, no NG concentration would be measured in the target. Instead, a thin near-surface region depleted from NG atoms is observed
[31,32]. It is therefore reasonable to assume that the ejected NG atoms originate near this layer. Since the layer is thin, the properties of the collision cascade at these positions are similar as at the surface, and the spatial distribution of the ejected atoms is similar as in conventional sputtering. D. Calculation and interpretation of the moments With assumption 1, the crater functions may be approximatedbysumsof δ functions
[1,3]  FI nn XX   F (r) =Xδ(r −rF) −δ(r −rI) , (10) XiXiXi=1 i=1 where X is the atomic volume of atom species X (NG or Si, in FFF F our case), r=(xXi,y ) denote the final positions of the n XiXiX II atoms of species X after the impact, rI =(x ,y ) the initial Xi XiXi positions of the n I atoms of species X before the impact, and X r =(x,y) the position where the crater function is evaluated. The
atoms
considered
in
Eq.
(10)
include
only
atoms
in
the
target. This means, n F includes the implanted ion if it ends X up in the target, and the redistributed atoms. n I includes the X eroded and redistributed atoms. The angle brackets denote the average over a sufficient number of impacts. Inserting
Eq.
(10)
into
Eq.
(1), one obtains  M(0) = X n F X −n I X (11) X   n F X n I X M(1) = X x F Xi − x I Xi (12) X i=1 i=1 The crater function as well as the moments may be decomposed into the contributions by the different atom species and physical mechanisms in an obvious way. The contributions of implantation and erosion to the mo- F ments have simple physical meaning: n impl is one minus NG the reflection coefficient rNG, thus (0) MNG,impl =NG(1 −rNG), (13) (0) I while M=0, see Eq. (7).
n is the partial X=NG,impl X eros sputtering yield YX, thus (0) M=−XYX. (14) X,eros Next, we observe that the contribution of ion implantation to the crater function is always positive, i.e., FNG,impl(x,y)  0 for all (0) x and y. FNG,impl(x,y)/MNG,impl may therefore be interpreted as the probability density that the end point of an ion trajectory has coordinates (x,y). It follows that (1) (0) xNG,impl =MNG,impl/MNG,impl (15) is the mean x coordinate of the ion trajectory end points. Similarly, FX,eros(x,y)  0 for all x and y, and FX,eros(x,y)/M(0) , X,eros which is always positive, defines the probability density that the origin of a sputtered atom of type X has coordinates (x,y). Therefore (1) (0) xX,eros =M/M(16) X,erosX,eros is the mean x coordinate of the origins of sputtered X atoms. Henceforth, we will call xNG,impl and xX,eros the mean projected implant and erosion distance, respectively, where distance is meant to be defined with respect to the impact point. Note that xNG,impl and xX,eros are determined solely by the geometry (0) (1) of the collision cascades, while Mand Mare NG,eros NG,eros proportional to the near-surface concentrations of the NG atoms. E. The moments in sputtering equilibrium Assumption 2 means that the average number of implanted atoms per incident ion, which is one minus the reflection coefficient, equals the partial sputtering yield of the ion species: 1 −rNG =YNG. (17) From
Eqs.
(13)
and
(14) follows (0) (0) M=−M(18) NG,eros NG,impl, and
because
of
Eqs.
(5)
and
(6) the contribution of the implant species (NG) to the zeroth moment vanishes: (0) MNG =0. (19) The total zeroth moment is therefore given by  (0) (0) (0) M=M=M, (20) X X,erosX=NG X=NG where
in
the
last
step
Eqs.
(6)
and
(7) have been used. Thus M(0) is simply minus the sum of the partial sputtering yields of the target atoms times the respective atomic volumes. (1) (1) Introducing Mfrom
Eqs.
(15)
and
(16), NG,impl and MNG,eros respectively,
in
Eq.
(5)
and
using
Eq.
(18), the contribution of the implanted ions to the first moment can be written: (1) (0) (1) MNG =(xNG,impl −xNG,eros)MNG,impl +MNG,redist. (21) 155307-3 HOBLER, MACIĄŻEK, AND POSTAWA PHYSICAL REVIEW B 97, 155307 (2018) For the redistributed atoms, the number of initial and final positions within the target are equal (nX,redist), and the contribution of redistribution of atom species X in
Eq.
(12)
can
be
written:
 nX,redist  (1) FI M=(x −x ) . (22) X,redist X Xi Xi i=1 The sum is over all x components of the displacement vectors of the redistributed atoms. Because of momentum conservation, the average displacement must be in the direction of the (1) incident ion. Therefore the sum is positive, and so is M X,redist. Specializing to X =NG, (1) M(23) NG,redist > 0 follows,
and
from
Eq.
(21),
(1) (0) MNG > (xNG,impl −xNG,eros)MNG,impl. (24) Equation
(24)
specifies
a
lower
limit
to
the
contribution
of
the implanted NG atoms to the first moment. The importance of the assumption of sputtering equilibrium lies in the fact that all quantities on the
right-hand
side
of
Eq.
(24)are
independent of the implanted ion concentration, if this is the case for the geometry of the collision cascade (assumption 3) and the spatial distribution of NG ejection (assumption 4). This means that we now may use static MC simulations with a target containing a small constant concentration of NG atoms to determine at least a lower limit to the contribution of the implanted NG atoms to the first moment, even in the absence of a rapid relocation model. The only missing term in
the
NG
contribution
to
the
first
moment,
Eq.
(21),
then
(1) is MNG,redist, which is proportional to the NG content. If, in addition, the contributions of the target atoms to the moments are independent of the NG content, then we may also estimate the other moment contributions using static MC simulations. Note that an imbalance of the implant and erosion related contributions to the first moment, and thus a nonvanishing right-hand
side
of
Eq.
(24),
implies
that
implantation
and
ero
sion are separate effects that have different spatial distributions. Our MD results presented in Sec. IV
C
validate this picture. III. SIMULATION METHODS A. Molecular dynamics simulations MD simulations were carried out at Jagiellonian University using LAMMPS [36].
The
simulation
cell
had
a
surface
area
of
28a ×20a(≈152A×109 =A is the lattice A), where a 5.431 constant of Si, and the sample thickness was adjusted as to accommodate most of the collision cascades. The cell was initially filled with single-crystalline (100)-Si with a (2 ×1) reconstructed surface. Periodic boundary conditions were used in the lateral directions. Stochastic and rigid layers, 7 A and 3 A thick, respectively, were used at the bottom to simulate the thermal bath that kept the sample at the required temperature and to keep the shape of the sample. The simulations were run at 0 K temperature. The target was sequentially bombarded with 2-keV Kr ions at polar angles of 40◦,50◦,60◦,70◦,80◦, and 85◦and an azimuthal angle of 0◦with respect to the cell edge, which is a [010] direction. In addition, one simulation was carried out each for 1-keV Ar and 2-keV Xe impacts at an incidence angle of 60◦. Each impact was simulated for 2 ps. The resulting structure was used as initial condition for the subsequent impact after removal of all sputtered atoms and any excess kinetic energy from the system. The latter was achieved by an energy quenching procedure that involved application of gentle viscotic forces to the entire sample for 5 ps. A Tersoff-3 potential
[37]
was
used
for
Si-Si
interactions,
and
the
ZBL
potential
[38]
for
all
other
interactions
(Kr-Kr,
Ar-Ar,
Xe-Xe,
Kr-Si, Ar-Si, Xe-Si). Since the ZBL potential is known to overestimate the interaction at large interatomic separations, we also performed a simulation with the Süle potential for the Kr-Si
and
Kr-Kr
interaction
[39,40], which has been fitted to ab initio calculations. The Kr-Kr Süle potential is close to the well
established
HFD-B2
potential
[41]
in
the
eV
energy
range.
The sequential impact simulations take several weeks to run on a few dozens of CPUs for each incidence angle. It is quite obvious that such simulations currently cannot be performed for much higher impact energies, e.g., 200 keV, at reasonable expense.
According
to
SRIM
[42],
the
projected
range
Rp of Kr ions increases by a factor of 23 when going from 2 to 200 keV. To estimate the simulation time required for the MD simulation of a 200-keV Kr impact we assume that the cell size has to be scaled proportional to Rp in each dimension. This means that the number of atoms in the simulation and thus, roughly, the simulation time per impact increases as R3, p i.e., by a factor of 233 ≈1.2 ×104. In addition, the number of impacts required to reach sputtering equilibrium increases by about the same factor: to reach a certain fluence, the number of impacts scales with R2 corresponding to the change in surface p area. To reach sputtering equilibrium, a depth of ∼Rp has to be
sputtered
[35].
Since
the
sputtering
yield
is
only
weakly
dependent
on
the
energy
in
the
keV
energy
range
[43],
the
fluence required to reach sputter equilibrium scales as Rp, and the number of impacts as R3. This motivates the use of MC p simulations to investigate a wider range of impact energies. B. Monte Carlo simulations The sequential MD impact simulations correspond to a “dynamic mode” in MC parlance in that the changes to the target induced by one impact are taken into account in the next impact. Dynamic simulations are also performed with MC for the comparison of the retained fluence (Sec. IV
A).
For the MC simulations we use the IMSIL code
[44,45]. The dynamic mode has been added to IMSIL some time ago based on the algorithm implemented in TRIDYN [46].
In
this
approach,
the substrate is subdivided into slabs, whose thicknesses are adjusted periodically as to relax the atom densities to their equilibrium values. In our simulations, we used an initial slab thickness of 1 A and perform target relaxation after every 100 −2 impacts, corresponding to a fluence increment of 1012 cm. For reasons explained in Sec. IIE,
the
MC
simulations
for
the
moments calculations are performed in static mode with a Si target containing a constant NG concentration of 2%. In the MC simulations, atoms interacted via the Ziegler-Biersack-Littmark (ZBL)
potential
[38].
Electronic stop-ping was calculated using a mixed Lindhard/Oen-Robinson model
[47,48] with Si parameters given in
Ref.
[49]
and
similar parameters for the NG. At low energies, electronic stopping plays only a minor role, since the Lindhard (nonlocal) 155307-4 CRATER FUNCTION MOMENTS: ROLE OF IMPLANTED … PHYSICAL REVIEW B 97, 155307 (2018) part
is
small
at
low
energies
[49],
and
the
apsis
of
collision
used in the Oen-Robinson model is always large. To obtain realistic
sputtering
yields,
a
planar
potential
barrier
[50,51] corresponding to a surface binding energy of Es =4.7eVis assumed for Si atoms and of Es =0.25 eV for the NG atoms. The former is chosen independent of the NG concentration contrary to the default model of IMSIL, as in reality we expect the
near-surface
region
to
deplete
from
the
NG
[32].
The
surface binding energy for the NG atoms was chosen so that reduction towards smaller values did not change sputtering yields and Es  Ed (see below) was fulfilled. Unnecessarily low values of Es should be avoided, since for the sake of a consistent model, trajectories are simulated down to a cutoff energy Ef equal to the minimum surface binding energy of all atoms, which leads to large simulation times if Es is small. For the displacement energy Ed [38,52], a value of 0.25 eV was used in order to fit the first redistributive moments of MC to those of the MD simulations as much as possible without having to accept excessive simulation times. The use of a very low value of Ed is in line with earlier MC simulations within the
crater
function
formalism
[23].
C. Calculation of moments and derivatives Since the crater functions are defined as the average response to a single ion, the moments have to be averaged over a number of impacts. In the MD simulations, we use 750 impacts for averaging, corresponding to a fluence increment of about 5 ×1014 cm−2. The 750 impacts are either chosen after the first 375 impacts (“low fluence”) or as the last 750 impacts of a simulation (“steady state”). In the MC simulations, averaging is done over all impacts (usually 100 000) since there is no transient in the simulation. As stated in Sec. II
A,
the
origin
of
the
coordinate
system
is defined by the impact point. Determination of the impact point in the MD simulations is complicated by the fact that surface roughness develops during bombardment (the typical RMS amplitude in our MD simulations is ∼2–3 A). In MD, we define the surface by moving a probe atom with radius 2.1 Aover the sample and taking the coordinates of the sample atom that is touched by the probe as the surface position. The impact point is taken as the intersection of the incoming ion direction with the average surface position. In the MC simulations, we define the surface position at the plane that divides the half spaces where target atoms are randomly generated or not. We note that the definition of the impact point is relevant to the values of the first moment contributions by implantation (1) (1) and erosion (Mand M, respectively), but it is NG,impl X,eros(1) irrelevant to the sum of the two NG contributions (M NG,impl + (1) M)
in
sputtering
equilibrium,
see
Eq.
(24).
When
com- NG,eros paring the individual contributions, however, e.g., between MD and MC, it is important to use consistent definitions. The derivative with respect to the incidence angle θ occurring
in
the
first
term
of
Eq.
(8)
is
calculated
by
fitting
parabolas
through three adjacent θ values, which is second-order accurate in θ [70].
Derivatives
of
the
sputtering
yield
with
respect
to
curvature are calculated by simulating sputtering from cylinders of radius R =5a, where a is the mean depth of energy deposition at normal incidence, and taking the derivative equal FIG. 1. Retained NG fluence as a function of implanted fluence for 1-keV Ar, 2-keV Kr, and 2-keV Xe bombardment of Si at an incidence angle of 60◦as obtained with MD (solid lines) and MC simulations (dashed lines). Note the significantly higher steady-state retained fluence in the MC simulations than in MD. to −R times the difference in sputtering yield between cylinder and
flat
target
[53].
It is difficult to do this with MD, since the dependence of the sputtering yield on the curvature becomes nonlinear for
relatively
little
changes
in
the
yield
[71].
So,
in
order
to approximate the derivative with respect to curvature by a difference quotient, small changes have to be evaluated, which require statistics that is expensive to obtain with MD. We will therefore use MC results for the terms T11 and T22 when reporting results for C11 and C22. IV. RESULTS A. Retained fluence Figure 1
shows the retained NG fluence as a function of implanted fluence for 1 keV Ar, 2 keV Kr, and 2 keV Xe bombardment of Si at an incidence angle of 60◦. While all curves start at the origin with a slope close to unity (dotted line), indicating that little reflection occurs at this angle, they saturate at different levels: the MC steady-state values are more than a factor of three larger than the MD values. This can be assigned
to
the
rapid
relocation
effect
discussed
in
Ref.
[31],
which is not included in the MC simulations. Figure 2
shows the Ar concentration depth profiles after 1 keV Ar bombardment of Si as predicted by MC and MD, compared to experimental data obtained by medium energy ion
scattering
(MEIS)
[54].
Obviously,
the
MD
simulations
describe a large fraction of the rapid relocation effect. To check whether the remaining difference is due to insufficient simulation times, we annealed one sample at 600 K for 10 ns, and found a reduction in the retained NG fluence of only ∼10% (this was done for 2-keV Kr at an incidence angle of 60◦). This modification seems low enough to be neglected in our study. It also might be a real effect due to the elevated temperature. 155307-5 HOBLER, MACIĄŻEK, AND POSTAWA PHYSICAL REVIEW B 97, 155307 (2018) FIG. 2. Steady-state Ar concentration profiles after 1-keV Ar bombardment of Si as obtained with MD (solid line), MC (dashed line),
and
MEIS
[54]
(dashed
line
with
symbols).
B. Atomic volume of the implanted species According
to
Eqs.
(10)to(12) the effect of the implanted ions on the crater function and its moments is proportional to their atomic volume. As a first guess, it might be plausible to use
the
densities
in
the
liquid
or
solid
phase
[55].
However,
bubbles are known to form during NG implantation, and the atomic density in the bubbles is uncertain. We have therefore determined the atomic density of the NG atoms by assuming 3 a fixed value of Si =20.4 Afor the atomic volume of amorphous
Si
[56]
and
fitting
the
NG
atomic
volume
NG to the densities observed in the MD simulations. This is done by counting the NG and Si atoms (nNG and nSi, respectively) in a layer of thickness 4a ≈21.7 A at a depth of 1–4 nm below the surface. The known volume V of the slab must equal V =nNGNG +nSiSi. (25) Plotting nSi versus nNG after each impact provides a point cloud to
which
Eq.
(25)
can
be
fitted
by
adjusting
the
only
unknown
NG. Since we are only interested in the steady state values, only the last 750 impacts are used for the fit. In Fig. 3(a),the
nSi(nNG) data are shown for 1-keV Ar, 2-keV Kr, and 2-keV Xe bombardment of Si at 60◦. Each data set starts with nNG =0, nSi =17 920, the situation before the first impact, and evolves towards the right. The first few impacts lead to a reduction in target density (strong decrease in nSi with only moderate increase in nNG), which is subsequently reversed. Closer inspection of the data shows that relaxation of the initial dilution occurs after 50–100 impacts (50 for Ar and −2 100 for Xe), corresponding to fluences around 5 ×1013 cm. This fluence is on the order of the amorphization threshold. The initial density reduction can be explained by the fact that target atoms are ejected from the near-surface region, where we measure the density, to either the vacuum or deeper parts of the target
[57].
The
subsequent
increase
of
the
target
density
is
due
to relaxation during amorphization. The last 750 data points in each set lead to the mean values and standard deviations of the atomic volumes given in the inset. The standard deviations might be slight underestimates caused by the limited impact FIG. 3. Number of Si atoms vs number of NG atoms in a constant volume V (a) after each impact for 1-keV Ar, 2-keV Kr, and 2-keV Xe bombardment of Si at an incidence angle of 60◦, and (b) after each of the last 750 impacts of 2-keV Kr at various incidence angles. The dotted lines are fits of the NG atomic volume to the last 750 impacts at 60◦assuming the atomic volume of Si to equal the experimental value of Si =20.4 A3 [56].
intervals used for averaging. Choosing different 750-impact intervals after the first ∼5 ×1014 cm−2 (not shown), the atomic volumes are stable to within ∼10%. Equation
(25)
with
Kr fitted to the 60◦data is shown in Fig. 3(b)
together with the MD data for incidence angles between 40◦and 85◦(last 750 impacts). For incidence angles up to 70◦the fit is excellent. For larger angles (80◦and 85◦), the MD data lie above the fit, which means that the material is denser than described by the fit. A possible explanation is that at the lower incidence angles bubbles form which have a somewhat lower density than small NG clusters and interstitial NG atoms. Bubbles are less likely to emerge at large incidence angles where the sputtering yield is high and the reflection coefficient is considerable. The results for the atomic volumes NG are plotted in Fig. 4
together with published experimental data of solid state atomic volumes
[55].
Our
data
exceed
the
solid
state
data
by
about
a
factor of two. To exclude an artifact of the ZBL interatomic potential, we have repeated the Kr simulation with the Süle potentials. These potentials are weaker at large interatomic separations, so one would expect to obtain a smaller atomic volume. This is not the case within the statistical error, see the red symbol in Fig. 4.
We
conclude
that
the
error
introduced
by
the ZBL potential, if any, is quite moderate. 155307-6 CRATER FUNCTION MOMENTS: ROLE OF IMPLANTED … PHYSICAL REVIEW B 97, 155307 (2018) C. Effect of the implanted ions on the crater function moments: MD results MD results for the contributions of implantation and erosion to the zeroth moment M(0) are shown in Fig. 5(a)
as a function of incidence angle for 2-keV Kr bombardment of Si. As mentioned in Sec. II,
atom
redistribution
does
not
contribute
to the zeroth moment. Both “low-fluence” (dashed lines) and “steady-state” results (solid lines) as defined in Sec. III
C
are shown. Notably, the contributions of Si erosion and Kr implantation have hardly any fluence dependence, while the contribution of Kr erosion only gradually builds up as the Kr ions are implanted. In the high-fluence case, the erosive contribution of Kr is approximately the negative of its implant contribution, indicating that steady state has indeed been reached to a good degree. Given these results, it is not surprising that the total zeroth moment in steady state [blue line labeled “total” in Fig. 5(b)]
agrees well with the values of Si erosion alone (green curve labeled “Si only”). The remaining small difference is probably due to the fact that sputtering equilibrium has not completely been reached. The added Si and implant contribution (red curve labeled “Si + Kr implant”) is not a good approximation of the total zeroth moment in steady state. As for the zeroth moment M(0), the contributions of Si erosion and Kr implantation to the first moment M(1) hardly depend on the fluence, see Fig. 6(a).
The
same
holds
true
for
the
contribution by Si redistribution. Kr redistribution plays only a minor role, while the effect of Kr erosion on M(1) depends on fluence as expected: the influence of Kr erosion develops only gradually as Kr ions are implanted. However, in contrast to its contribution to M(0), the effect of Kr erosion is not completely compensated for in steady state by the effect of implantation. FIG. 5. MD results for the zeroth moment M(0) of the crater function of Si bombarded with 2 keV Kr vs incidence angle. (a) Contributions by implantation and erosion comparing low fluence and steady state; (b) sum of contributions in steady state. The difference between the blue solid line and the green dashed line in panel (b) should vanish in sputtering equilibrium. According
to
Eq.
(21),
this
would
be
the
case
if
the
mean
projected implant distance xKr,impl equaled the mean projected (1) (1) erosion distance xKr,eros. Rather, MKr,impl > |MKr,eros|can be read from Fig. 6(a),
and
therefore
xKr,impl >xKr,eros validating, at least qualitatively, assumption 4. As a consequence, the first moment in steady state does not agree well with the contributions of Si alone [solid blue and dashed green line, respectively, in Fig. 6(b)].
This
is
an
important conclusion: the contribution of the implanted NG ions to the first moment is significant, on the order of 50% of the Si contribution in the present case. Adding the Kr implant contribution to the Si contribution overestimates the total first moment M(1) [dotted red line in Fig. 6(b)].
Further
adding
the
negative
contribution
of
Kr
erosion, thus neglecting only the redistributive contribution of the NG atoms, gives a lower limit to the total first moment M(1) (magenta
line
with
triangles),
see
Eq.
(24).
For
large
angles
it
is a good approximation of M(1), while the deviation increases with decreasing incidence angle. In any case, it is a better approximation to M(1) than neglecting the Kr contributions altogether (dashed green line). In Fig. 6(c),
the curvature coefficient S11 is plotted as calculated from the data given in Fig. 6(b).Asfor
M(1),the contribution of Kr to S11 is significant. Considering only the implant and erosive contributions (magenta line with triangles) is a reasonable approximation to the total values, although the deviation increases with decreasing incidence angle. 155307-7 HOBLER, MACIĄŻEK, AND POSTAWA PHYSICAL REVIEW B 97, 155307 (2018) FIG. 6. MD results for the first moment M(1) of the crater function [(a) and (b)] and for the curvature coefficient S11 as defined in Sec. II
(c) for Si bombarded with 2 keV Kr as a function of incidence angle. (a) Contributions by implantation, erosion, and redistribution; [(b) and (c)] sum of all contributions. For the meaning of the line styles see Fig. 5.
In
addition,
results
neglecting
only
the
Kr
erosive
contribution
[Eq.
(24)]
are
shown
by
the
magenta
lines
with
triangles.
The differences between the solid blue line and the dashed green line in (b) and (c) indicate the role of the implanted Kr ions. Similar results for both the zeroth and first moment are also obtained for Ar and Xe ions, see Fig. 7.
In
all
cases,
the
contribution of the implanted ions is significant, and the lower limit
according
to
Eq.
(24)
is
the
best
approximation
to
the
full
calculation. Figure 8
shows three versions of the curvature coefficients C11 and C22 for the Kr case: according to the original crater function formalism considering only flat targets and
Si
atoms
[1,21], according to the extended crater function formalism considering in addition the effect of surface curvature
[2],
and
considering
in
addition
the
contributions
of
the Kr atoms as proposed in this work. As is obvious from FIG. 7. Analogous to Figs. 5(b)
and 6(b),
but
as
a
function
of
ion
species for an incidence angle of 60◦. For the meaning of the line styles see Fig. 6.
the plot, all contributions are significant. For C11, the effect of surface curvature (the difference between the red dashed line and the green dotted line) partly compensates for the effect of the Kr atoms (difference between blue solid line and red dashed line), while for C22 both contributions work in the same direction. D. MC versus MD results As discussed in Sec. IIIA,
the
computational
expense
of
MD simulations increases dramatically with increasing impact energy. To investigate the energy dependence of the role of the implanted ions, we have therefore performed MC simulations as described in Sec. III
B.
We
recall
that
it
is
not
possible
to
reliably predict the NG concentration in the target using MC simulations. As a consequence, there is uncertainty for the redistributive and erosive NG contributions to the first moment. The redistributive contribution has turned out to be small in the cases studied [Figs. 6(a)
and 7(b)]
and
will
therefore
be neglected. The erosive contribution can be approximated assuming sputtering equilibrium and that the mean erosion distance does not depend on the details of the NG profile in the target. We use a Si target with a constant Kr concentration of 2% for our calculations. In Fig. 9,
MD
and
MC
results
of
the
contributions
to
the
first
moment
according
to
Eq.
(24)
are
compared
(magenta
triangles
connected by solid and dashed lines, respectively), showing good agreement. Towards lower incidence angles this approximation increasingly underestimates the total Kr contribution [blue circles, Fig. 9(a)]
as
the
Kr
redistributive
contribution
155307-8 CRATER FUNCTION MOMENTS: ROLE OF IMPLANTED … PHYSICAL REVIEW B 97, 155307 (2018) increases [compare also Fig. 6(b)].
For
comparison,
the
Si
redistributive contributions are also shown. The MD and MC results are similar, but again somewhat diverge towards smaller incidence angles. Figure 9(b)
shows similar agreement for the lighter Ar and the heavier Xe ions. The ion mass dependence of the Si redistributive contribution is somewhat underestimated by
MC,
which
may
be
explained
by
a
spike
effect
[58]
not
included in the MC simulations. The
two
critical
quantities
in
Eq.
(24)
are
the
mean
projected
distances of the implant positions xNG,impl and the erosion points xNG,eros from the impact point. The distributions of these distances for 60◦Kr impacts at 2 keV are plotted in Fig. 10(a).
Good agreement between the MD and MC results is observed. Very clearly, the average implant distance is larger than the average erosion distance, i.e., implantation takes place considerably further away from the impact point than sputtering. This is plausible as nuclear energy is preferably deposited in a narrow region close to the impact point (see Fig. 15 of
Ref.
[59]).
According
to
Eq.
(21), xNG,impl >xNG,eros > 0 means that the added contributions of NG erosion and implantation are positive but smaller than the contribution by implantation alone. The agreement between the MD results, which include the rapid relocation effect, and the MC results, which do not, is remarkable also because it confirms the picture proposed in as- (b) FIG. 9. Contributions of the NG (blue circles) and Si (green squares) to the first moment M(1) as obtained by MD (solid lines) and MC simulations (dashed lines). In addition, the implant plus erosive contributions of Kr are shown by the magenta triangles. (a) 2-keV Kr as a function of incidence angle; (b) for an incidence angle of 60◦as a function of ion species (1-keV Ar, 2-keV Kr, and 2-keV Xe). sumption 4 of Sec. II
C:
while
the
NG
atoms
are
driven
towards
the surface by rapid relocation, their eventual sputtering has similar spatial properties as in the absence of rapid relocation. In Fig. 10(b),
these
quantities
are
compared
as
a
function
of incidence angle. Excellent agreement is found for angles up to 60◦, while MC and MD results moderately diverge at larger angles. However, all that matters to the added effect of implantation and erosion is the difference between xNG,impl and xNG,eros, which agrees well between MC and MD also for the larger angles. This also explains why the added contributions of Kr implantation and Kr erosion to the first moment (magenta lines in Fig. 9)
agree
well
between
MD
and
MC.
The
MD
results in Fig. 10(b)
also show that the mean projected implant and erosion distances hardly depend on the fluence, compare the solid and dotted lines. This confirms that xNG,impl and xNG,eros are robust quantities that may be calculated without exact knowledge of the NG concentration in the target. 155307-9 HOBLER, MACIĄŻEK, AND POSTAWA PHYSICAL REVIEW B 97, 155307 (2018) (a) (b) FIG. 10. (a) Distribution of the projected distance from the impact point (=x coordinate) for 2-keV Kr bombardment of Si at an incidence angle of 60◦; (b) mean projected distance from the impact point as a function of incidence angle for 2 keV Kr bombardment of Si. Red lines: implanted ions; green lines: sputtered Kr atoms. Solid lines: MD in steady state; dashed lines: MC; dotted lines: MD at low fluence. In (a), the MD results have been averaged over the last 1500 ion impacts. E. MC results for a wider range of impact conditions The ion species and energy dependence of xNG,impl and xNG,eros as calculated by our MC simulations for an incidence angle of 60◦is shown in Fig. 11.
Note
that
always
xNG,impl > xNG,eros. According
to
Eq.
(24)
this
means
that
the
added
contributions of NG implantation and erosion on the first moment are always positive. Both quantities increase with energy, and this is also true of their ratio xNG,impl/xNG,eros. While xNG,impl/xNG,eros is only 1.33, 1.57, and 1.88 at 200 eV for Ar, Kr, and Xe ions, respectively, it exceeds a factor of three at 200 keV. The importance of the implanted ions to the crater function must of course be related to the contributions of the target atoms. In Fig. 12,
we
therefore
compare
the
curvature
coefficients C11 and C22 with and without consideration of FIG. 11. MC results for the mean projected distances from the impact point for the implanted ions (xNG,impl, red lines and symbols) and for the sputtered atoms (xNG,eros, green lines and symbols) as a function of impact energy. Data for the three ion species Ar, Kr, and Xe at an incidence angle of 60◦are shown. the NG atoms (solid and dashed lines, respectively), where the contribution of the NG atoms has been estimated using Eq.
(24).
Data
are
given
as
a
function
of
incidence
angle
for
three ion species (Ar, Kr, and Xe) and four energies (0.2, 2, 20, and 200 keV). The effect of the implanted NG ions is significant in all cases, increases with impact energy although the increase is attenuated for the heavier ions. In all cases the ions have a destabilizing effect on the surface with respect to parallel mode ripples (contribution to C11 < 0) for incidence angles larger than 40◦–50◦up to at least 85◦. They do not change the lower critical angle for parallel mode ripple formation [zero of C11(θ)] significantly, but move the upper critical angle to higher values. On the other hand, the ions always have a stabilizing effect on the surface with respect to perpendicular mode ripples (contribution to C22 > 0). These results are qualitatively the same
as
for
self-implantation
[27,28]. V. DISCUSSION AND CONCLUSION From the results presented in the previous section, we conclude that the implanted NG ions in general play an important role in determining crater function moments. The rationale may be summarized as follows: (i) the average projected distance of the implanted ions from the impact point, measured along the surface, is larger than the corresponding average distance of the sputtered atoms. This leads to incomplete compensation of the effects of ion implantation and erosion on the first moment. (ii) Redistribution of the implanted ions due to subsequent impacts is always away from the impact point and therefore can only add to (rather than compensate for) the effect of ion implantation. (iii) The average atomic volumes of the NG atoms have been found in our MD simulations to be approximately twice those of the NG atoms in the solid state. This means that their effect on the crater function moments is twice that one would obtain with the solid state volumes. There are some uncertainties in the simulations; however, they do not compromise our conclusion. First, MD may only 155307-10 CRATER FUNCTION MOMENTS: ROLE OF IMPLANTED … PHYSICAL REVIEW B 97, 155307 (2018) FIG. 12. MC results for the curvature coefficients (a) C11 and (b) C22, comparing simulations with and without consideration of the effect of the implanted NG ions (solid and dashed lines, respectively). Data are shown for Si bombarded with Ar, Kr, and Xe at energies of 0.2, 2, 20, and 200 keV as a function of incidence angle. The results have been scaled with the square root of the impact energy to better represent them in the plots. describe effects that occur on the timescale of the collision cascades. Thermally activated processes that mainly operate between the collision cascades are not included. While a systematic understanding of the effect of temperature on pattern formation
is
still
lacking
[60],
there
are
indications
that
tem
perature does have an effect, even not too far from room temperature
[61–65]. Also, thermally activated NG diffusion could explain the differences between the Ar concentration profiles obtained by MD and experimentally, see Fig. 2,
although
our
annealing simulation did not provide enough evidence for a thermal effect. Anyway, thermal effects do not influence the range of the implanted ions and thus the implant contribution to the first crater function moment. Moreover, NG redistribution would be reduced, if the NG concentrations in the target were decreased by thermal diffusion, but our argument is not based on the redistributive NG contribution, which is small in our MD simulations and has been neglected in the MC simulations. Second, MC simulations have their inaccuracies due to the binary collision approximation. As a result, the ion mass dependence of redistribution is underestimated by MC, see the Si redistribution results in Fig. 9(b).
This
is
reminiscent
of
the ion mass dependence of damage formation in crystalline Si, which is explained by the increasing relevance of thermal spikes
with
increasing
ion
mass
[58].
However,
the
effect
is
moderate and is not expected to increase with ion energy: higher impact energies add high-energy portions to the ion trajectories, which are well described in the binary collision approximation. Another consequence of the binary collision approximation is its failure to describe the rapid relocation effect that leads to an overestimation of the retained NG fluence, see Fig. 2.The
retained
fluence
is
expected
to
increase
with
ion
energy
[31].
Therefore one may expect an increasing contribution of NG redistribution to the first crater function moment with increasing energy. This would further enhance the importance of the implanted NG atoms beyond the analysis given in Sec. IV
E.
The argument is not so straightforward, however, since the other contributions to the first moment also increase with energy. The issue could be clarified by using a rapid relocation model
[33]
in
the
MC
simulations
that
is
fitted
to
experimental
data on the retained fluence where such data is available. Again, NG redistribution would only add to the effects of implantation and erosion, so it would only increase the importance of the implanted ions. 155307-11 HOBLER, MACIĄŻEK, AND POSTAWA PHYSICAL REVIEW B 97, 155307 (2018) As mentioned in Sec. II
C,
our
MC
simulations
are
based
on the assumption that it is not necessary to know the exact concentration of the implanted NG atom (or that the actual concentration is close to the assumed one). This assumption is not perfectly valid. For instance, the sputtering yield increases by 25% from a pure Si target to steady state for 140-keV or 270-keV Xe bombardments
[66].
This
is
explained by a reduction in ion and recoil range by the implanted heavy Xe atoms. This also means that the mean projected implant and erosion distances (xNG,impl and xNG,eros) both shrink, which would not invalidate our analysis. It should also be taken into account
that
the
ion
species
and
energy
in
this
study
[66]
represent extreme cases in view of the present study, and that this experiment was done at normal incidence, while we are more interested in oblique incidence where NG gas retention is reduced. Before finishing, a few words on how our calculated curvature coefficients compare to experiments seem appropriate. From a qualitative point of view, all levels of sophistication of the crater function formalism [considering planar surfaces only
[1],
including
the
effect
of
curvature
[2], and including the effect of the implanted NG atoms (this work)] agree with experiments
performed
under
impurity-free
conditions
[60],
in
that they predict parallel mode ripple formation (C11 < 0) for incidence angles above a critical angle around 45◦up to another critical angle at grazing incidence. On closer examination, several discrepancies manifest themselves. For instance, the very careful experiments on 2-keV Kr bombardment of Si by Engler et al.
[67]
yielded
a
lower
critical
angle
of
58◦, significantly larger than 45◦. For Ar ions, it has been found that
no
ripples
occur
for
ion
energies
of
3
to
10
keV
[26]in
contrast to the predictions of the crater function formalisms. Finally, perpendicular mode ripples have been measured for 10-keV Xe bombardment at an incidence angle of 80◦[68],
while our MC results do not predict perpendicular mode ripple formation (C22 < 0) at any of the investigated conditions. Unfortunately, we cannot report resolution of these discrepancies due to consideration of the effect of the implanted ions. On the other hand, the crater function formalism is intuitive, and MD and MC simulations are well-established techniques. We therefore believe that the cause of the discrepancy lies outside the crater function formalism, possibly in stress effects operating in the amorphous layer on a scale not covered
by
MD
simulations
[15–18]. Finally, we wish to comment on a recent study of pattern formation
by
1-MeV
Au
in
Ti
and
its
alloy
TiAlV
[69].
Applying our MC analysis to these conditions, we arrive at exactly
the
same
conclusions
as
the
authors
of
Ref.
[69].
The
relative importance of the implanted ions is close to negligible (results not shown). Au is about 50% heavier than Xe, which might be the cause why the implanted ions can be neglected in spite of the high impact energy. This also indicates that our conclusion about the importance of the implanted NG ions might not apply to Rn, the only NG heavier than Xe, which is even heavier than Au. Au, of course, is not a NG, and Au implanted at this high energy is substantially incorporated in the target as demonstrated
in
Ref.
[69]
by
RBS
measurements.
Such
conditions
should be analyzed by dynamic MC simulations, which would predict the implanted ion depth profile. In these simulations the redistributive contribution of the implanted ions to the first moment would be significant. This contribution could be determined with confidence, since depth profiles of nonvolatile ions are well predicted by MC simulations; the “rapid relocation” effect is a phenomenon associated with NG atoms only. Our main finding in this work is that due to the difference between the mean projected implant and erosion distances (xNG,impl −xNG,eros) the implanted NG ions usually play a role comparable to those of sputtering and atom redistribution as evaluated by crater functions. This is the cases even though NG atoms are rapidly removed from the target so that their concentration is low at any time. While there are still open questions in the modeling of ion-target interaction that must be solved to make quantitative predictions, the fact that xNG,impl −xNG,eros agrees well between MD and MC simulations (Fig. 10)
is strong support for this conclusion to be true not only at the conditions investigated by MD but also at higher energies. ACKNOWLEDGMENTS We gratefully acknowledge R. M. Bradley for helpful discussions. D.M. and Z.P. would like to gratefully acknowledge the financial support from the Polish National Science Center, Program 2015/19/B/ST4/01892. [1] S. A. Norris, M. P. Brenner, and M. J. Aziz, J. Phys.: Condens. Matter
21,
224017
(2009).
[2] M. P. Harrison and R. M. Bradley, Phys.
Rev.
B
89,
245401
(2014).
[3] S. A. Norris, J. Samela, M. Vestberg, K. Nordlund, and M. J. Aziz, Nucl.
Instrum.
Methods
B
318,
245
(2014).
[4] P. Sigmund, J.
Mater.
Sci.
8,
1545
(1973).
[5] R.M.Bradley andJ.M.E.Harper, J.
Vac.
Sci.
Technol.
A
6,
2390
(1988).
[6] R. Cuerno and A.-L. Barabási, Phys.
Rev.
Lett.
74,
4746
(1995).
[7] M. A. Makeev, R. Cuerno, and A.-L. Barabasi, Nucl. Instrum. Methods
Phys.
Res.
B
197,
185
(2002).
[8] J. Muñoz-García, R. Cuerno, and M. Castro, Phys.
Rev.
B
78,
205408
(2008).
[9] G. Carter and V. Vishnyakov, Phys.
Rev.
B
54,
17647
(1996).
[10] T. Aste and U. Valbusa, New
J.
Phys.
7,
122
(2005).
[11] B. Davidovitch, M. J. Aziz, and M. P. Brenner, Phys.Rev.B
76,
205420
(2007).
[12] O. Bobes, K. Zhang, and H. Hofsäss, Phys.
Rev.
B
86,
235414
(2012).
[13] S. Numazawa and R. Smith, J.
Phys.:
Condens.
Matter
25,
095003
(2013).
[14] N. V. Medhekar, W. L. Chan, V. B. Shenoy, and E. Chason, J.
Phys.:
Condens.
Matter
21,
224021
(2009).
[15] M. Castro, R. Gago, L. Vázquez, J. Muñoz-García, and R. Cuerno, Phys.
Rev.
B
86,
214107
(2012).
[16] M. Castro and R. Cuerno, Appl.
Surf.
Sci.
258,
4171
(2012).
[17] S. A. Norris, Phys.
Rev.
B
86,
235405
(2012).
155307-12 CRATER FUNCTION MOMENTS: ROLE OF IMPLANTED … PHYSICAL REVIEW B 97, 155307 (2018) [18] S. A. Norris, J. C. Perkinson, M. Mokhtarzadeh, E. Anzenberg, M. J. Aziz, and K. F. Ludwig, Sci.
Rep.
7,
2016
(2017).
[19] N. Kalyanasundaram, M. Ghazisaeidi, J. B. Freund, and H. T. Johnson, Appl.
Phys.
Lett.
92,
131909
(2008).
[20] N. Kalyanasundaram, J. B. Freund, and H. T. Johnson, J. Phys.: Condens.
Matter
21,
224018
(2009).
[21] S. A. Norris, J. Samela, L. Bukonte, M. Backman, F. Djurabekova, K. Nordlund, C. S. Madi, M. P. Brenner, and M. J. Aziz, Nat.
Commun.
2,
276
(2011).
[22] M. Z. Hossain, K. Das, J. B. Freund, and H. T. Johnson, Appl. Phys.
Lett.
99,
151913
(2011).
[23] H. Hofsäss, Appl.
Phys.
A
114,
401
(2014).
[24] H. Hofsäss, Appl.
Phys.
A
119,
687
(2015).
[25] O. El-Atwani, S. A. Norris, K. Ludwig, S. Gonderman, and J. P. Allain, Sci.
Rep.
5,
18207
(2015).
[26] H. Hofsäss, O. Bobes, and K. Zhang, J.
Appl.
Phys.
119,
035302
(2016).
[27] R. M. Bradley and H. Hofsäss, J.
Appl.
Phys.
120,
074302
(2016).
[28] H. Hofsäss, K. Zhang, and O. Bobes, J.
Appl.
Phys.
120,
135308
(2016).
[29] P. Blank, K. Wittmaack, and F. Schulz, Nucl.
Instrum.
Methods
132,
387
(1976).
[30] N. Menzel and K. Wittmaack, Nucl.
Instrum.
Methods
B
7,
366
(1985).
[31] K. Wittmaack, Nucl.
Instrum.
Methods
B
267,
2846
(2009).
[32] K. Wittmaack, A. Giordani, R. Umbel, and J. L. Hunter, Jr., J.
Vac.
Sci.
Technol.
A
34,
051404
(2016).
[33] A. Mutzke and W. Eckstein, Nucl.
Instrum.
Methods
B
266,
872
(2008).
[34] K. Wittmaack and H. Oppolzer, Nucl.
Instrum.
Methods
B
269,
380
(2011).
[35] M. Nastasi and J. W. Mayer, Ion Implantation and Synthesis of Materials (Springer, Berlin, Heidelberg, 2006). [36] S. Plimpton, J.
Comput.
Phys.
117,
1
(1995).
[37] J. Tersoff, Phys.Rev.B
39,
5566
(1989).
[38] J. F. Ziegler, J. P. Biersack, and U. Littmark, The Stopping and Range of Ions in Solids (Pergamon Press, New York, 1985). [39] P. Süle and K.-H. Heinig, J.
Chem.
Phys.
131,
204704
(2009).
[40] P. Süle, arXiv:1401.1612
[cond-mat.mtrl-sci]. [41] A. K. Dham, W. J. Meath, A. R. Allnatt, R. A. Aziz, and M. J. Slaman, Chem.
Phys.
142,
173
(1990).
[42] J. F. Ziegler, SRIM-2013 software package, available online at http://www.srim.org.
[43] W. Eckstein, in Sputtering by Particle Bombardment,editedby R. Behrisch and W. Eckstein, Topics in Applied Physics, Vol. 110 (Springer, Berlin, Heidelberg, 2007), pp. 33–189. [44] G. Hobler, Nucl.
Instrum.
Methods
B
96,
155
(1995).
[45] C. Ebm and G. Hobler, Nucl.
Instrum.
Methods
B
267,
2987
(2009).
[46] W. Möller, W. Eckstein, and J. Biersack, Comput. Phys. Com-mun.
51,
355
(1988).
[47] J. Lindhard and M. Scharff, Phys.
Rev.
124,
128
(1961).
[48] O. S. Oen and M. T. Robinson, Nucl.
Instrum.
Methods
132,
647
(1976).
[49] G. Hobler and C. Murthy, in 2000 International Conference on Ion Implantation Technology Proceedings (IEEE, Piscataway, 2000), pp. 209–212. [50] P. Sigmund, Phys.
Rev.
184,
383
(1969).
[51] W. Eckstein, Computer Simulation of Ion-Solid Interactions (Springer, Berlin, 1991). [52] G. H. Kinchin and R. S. Pease, Rep.
Prog.
Phys.
18,
1
(1955).
[53] H. M. Urbassek, R. M. Bradley, M. L. Nietiadi, and W. Möller, Phys.
Rev.
B
91,
165418
(2015).
[54] D. W. Oh, S. K. Oh, H. J. Kang, H. I. Lee, and D. W. Moon, Nucl.
Instrum.
Methods
B
190,
598
(2002).
[55] WebElements, http://www.webelements.com,
accessed
20
June,
2017. [56] J. S. Custer, M. O. Thompson, D. C. Jacobson, J. M. Poate, S. Roorda, W. C. Sinke, and F. Spaepen, Appl.
Phys.
Lett.
64,
437
(1994).
[57] G. Hobler, in Process Physics and Modeling in Semiconductor Technology, edited by G. R. Srinivasan, C. S. Murthy, and S. T. Dunham (The Electrochem. Soc., Pennington, 1996), pp. 509– 521. [58] M.-J. Caturla, T. Díaz de la Rubia, L. A. Marqués, and G. H. Gilmer, Phys.Rev.B
54,
16683
(1996).
[59] G. Hobler, R. M. Bradley, and H. M. Urbassek, Phys.Rev.B
93,
205443
(2016).
[60] J. Muñoz García, L. Vázquez, M. Castro, R. Gago, A. Redondo-Cubero, A. Moreno-Barrado, and R. Cuerno, Mater.
Sci.
Eng.
R
86,
1
(2014).
[61] G. Carter, V. Vishnyakov, Y. V. Martynenko, and M. J. Nobes, J. Appl.
Phys.
78,
3559
(1995).
[62] G. Carter, V. Vishnyakov, and M. J. Nobes, Nucl. Instrum. Methods
B
115,
440
(1996).
[63] C. C. Umbach, R. L. Headrick, and K.-C. Chang, Phys. Rev. Lett.
87,
246104
(2001).
[64] R. Banerjee, S. Hazra, and M. K. Sanyal, J.
Phys.
D:
Appl.
Phys.
41,
055306
(2008).
[65] D. Chowdhury, B. Satpati, and D. Ghose, Mater.
Res.
Express
3,
125003
(2016).
[66] P. Blank and K. Wittmaack, J.
Appl.
Phys.
50,
1519
(1979).
[67] M. Engler, S. Macko, F. Frost, and T. Michely, Phys.Rev.B
89,
245412
(2014).
[68] K. Zhang, F. Rotter, M. Uhrmacher, C. Ronning, H. Hofsäss, and J. Krauser, Surf.
Coat.
Technol.
201,
8299
(2007).
[69] M. A. Garcia, J. Rickards, R. Cuerno, R. Trejo-Luna, J. Cañetas-Ortega, L. R. de la Vega, and L. Rodríguez-Fernández, Phys. Rev. Appl.
8,
064027
(2017).
[70] This algorithm is also used in the gradient function of NUMPY, see http://www.numpy.org.
[71] For example, for 20-keV Ar in Si, nonlinearity with respect to curvature becomes significant when the sputtering yield has changed by ∼10%,
see
Fig.
4(a)
of
Ref.
[53].
155307-13 Nuclear Inst. and Methods in Physics Research B 447 (2019) 30–33
ContentslistsavailableatScienceDirect NuclearInst.andMethodsinPhysicsResearchB journal homepage: www.elsevier.com/locate/nimb Ionbombardmentinducedatomredistributioninamorphoustargets:MD versusBCA G.Hoblera,⁎,D.Maciążekb,Z.Postawab aInstituteofSolidStateElectronics,TUWien,Gußhausstraße25-25a,A-1040Wien,Austria bInstituteofPhysics,JagiellonianUniversity,ul.Lojasiewicza11,30348Kraków,Poland ARTICLEINFO Keywords: SpontaneouspatternformationAtomredistribution BinarycollisionapproximationMoleculardynamics ABSTRACT Atomredistributionisanimportantmechanismofionbombardmentinducedspontaneouspatternformationinamorphousoramorphizedtargets.Itmaybecharacterizedbyeithermoleculardynamics(MD)simulationsorbyMonte Carlo (MC) simulations based on the binary collision approximation (BCA). In this work we analyzeproblemsoftheBCAapproachinpredictingatomredistributionbycomparingMCandMDsimulationsofrecoileventsinamorphizedSi.WefindthattheMCresultscriticallydependonthedisplacementenergyused,onthechoice of the free flight paths and maximum impact parameters, and on the construction of the trajectories. Moreover,thenetatomredistributionas determinedbyMDissignificantlylargerforrecoilsstartingnear the surfacethanforbulkrecoils.TheeffectisnotreproducedbytheMCsimulations. 1. Introduction Spontaneouspatternformationbyenergeticionbeamshasreceivedsignificantimpetusinrecentyearsfromprogressintheunderstandingofitsmechanisms[1–12].Animportantcontributiontothefieldistheso-calledcraterfunctionformalism [2,7,12],whichderivesthecoefficientsoftheequationofsurfacemotionfromthemomentsof acrater function. These moments can be calculated by molecular dynamics(MD) simulations or by Monte Carlo (MC) simulations based on thebinarycollisionapproximation(BCA). The crater function moments have contributions from sputtering[13],atomredistributionwithinthetarget [14],andionimplantation [10,12].Often,theeffectofatomredistributiondominatesoratleastisan important contribution. It enters the formalism through its contributiontothefirstcraterfunctionmoment (1) Here, a Si target has been assumed and the contribution of relocatedimplantedionshasbeenneglected. denotestheatomicvolumeofSi, thenumberofSiatomsredistributedwithinthetarget,and their displacement component in the direction of the projection of the ionbeamtothesurface( axis).Theredistributivecontributiontothefirstmoment thus describes the net displacement of the target atoms indirection. The purpose of this paper is to analyze problems of the BCA approach in determining the sum of the displacements appearing in Eq. (1). Shortcomings of the BCA compared to MD with respect to atomredistributionhavebeenpointedoutbefore[15,16].Ontheotherhand,MC simulations using the BCA have successfully been used in combinationwiththecraterfunctionformalism [17,18].Thediscrepancyislikely due to different parameters and implementations of the MCcodes.Inthe present work wefocus ourattentionontherolesof dis-placement energy, free flight path selection, and trajectory construction. 2. Methodology We assess the qualtity of BCA implementations by comparison ofMC and MD simulations of atom redistribution resulting from recoilsstartingwithgivenenergy ,directionwithrespecttothesurface,and depth belowthesurface(Fig.1).ForMDthechoiceoftheinitialatomconfigurationiscritical.Atfluencessufficienttoproducepatterns,Siisamorphized. A not well relaxed sample may lead to unrealistic displacements[4,16].Oneapproachthereforeistoexerciseparticularcarewhen relaxing the initial atom configuration [16]. We have alternativelyusedsequentialionimpactsonaninitiallyvirgintargetuntilasteady state has been reached [12]. We use the target from Ref. [12]whichhasbeenproducedby2keVKrimpactswithanincidenceangleof 60° on a Si cell with dimensions approximately equal to ⁎ Correspondingauthor.E-mailaddress:gerhard.hobler@tuwien.ac.at(G.Hobler). https://doi.org/10.1016/j.nimb.2019.03.028 Received15August2018;Accepted15March2019 Available online 27 March 2019 0168-583X/ © 2019 Elsevier B.V. All rights reserved. G.Hobler,etal. Nuclear Inst. and Methods in Physics Research B 447 (2019) 30–33 Fig.1.Visualrepresentationofthemastersampleusedinthesimulationsoftherecoil events. Si atoms are represented as a yellow transparent isosurface, Kratoms as red spheres. A chosen primary recoil is represented by the blacksphere. ,compare Fig.1.Thetargetcontains afewpercent of Kr atoms with a profile decaying towards the surface and thebulk.WeassumetheirinfluenceontheredistributionofSiatomstobe negligible.Thetargetalsoexhibitssurfaceroughnessdevelopedduringthesequentialimpacts. Simulationsoftheprimaryrecoileventswithenergiesbetween5eVand100eV areperformedforthreeorientationsofthevelocityvector (Fig.1),paralleltothe axis,at60°pointingtowardsthesurface,andat60°pointingintothesample.Foranygivenenergyandorientation,300 simulations are carried out where atoms with initial velocity are chosen randomly from a region 60 Å thick starting from the surface. Eachsimulationstartswitha perfectcopyofthemastersampleat0K temperatureandis runfor1ps.Theappliedboundaryconditionsand interatomicpotentials arethesame asinRef.[12]. FortheMCsimulations weusetheIMSILcode [19].Unlessotherwise noted we use the same models and parameters as in Ref. [12]. Unlikein[12],weneglectelectronicstoppinginordertoallowamoreconsistent comparison with the MD results. In addition, we introduce variationsintotheBCAmodel asfollows:First, in [12]themaximum impactparameterischosensuchthatnocollisionswithenergytransferexceeding the displacement energy are missed. This requires very large impact parameters on the order of 3–4Å at low recoil and displacement energies. The value of is calculated from the minimumenergytransferandionenergyasproposedin [20].Thefree flightpath ischosenaccordingto (2) inordertousethecorrectatomicdensity [21].Incontrast,oftensimply a constant free flight path of is used together with Eq. (2) [16,22] . For Si this leads to a fixed maximum impact parameter of,which weconsiderasanoptioninoursimulations. Second,ionandrecoiltrajectoriesmaybeconstructedindifferentways.Themostaccuratetreatmentof abinarycollisioncalculatestheintersectionpointoftheasymptotesusingthe“timeintegral” [23,24]. InFig.2 thispointisdenoted1’,anditsdistancefromthefootoftheimpactparameterontheincomingasymptoteisdenoted .Puttingtheion at this point after a collision is the standard treatment in IMSIL,which we therefore label as “IMSIL default”. Note that the initiallychosenfreeflightpath isthedistancebetweentheion(1)andthefootof the impact parameter on the incoming asymptote. The actual free flightpathisreducedby toffp asshowninFig.2a. Atlowenergies, maybequitelarge[24],andatthesametimethe free flight paths are small due to the energy transfer criterion.Thereforetheactualfreeflightpathsffpmaybecomenegative,meaningthat the ion is moved in the opposite direction of its momentum(cf.Fig.2b).Toavoidthisunphysicalsituation,thefreeflightpathsffp werelimitedtonon-negativevaluesinRef.[12].Wedenotethismodel as“ ”.It comesattheexpenseofunphysicallyshiftingtheoutgoingasymptoteintheforwarddirection. Fig.2.Asymptotesoftheion(1)andrecoil(2)trajectories.Theionisassumedtoinitiallymovetotheright.Inthestandardmodel,afterthecollision,ionandrecoil are placed into the intersection points of incoming and outgoingasymptotes. This normally leads to positive actual free flight paths ffp for a givenchosenfreeflightpath (panela),butcanalsoleadtomovementinthedirectionoppositetotheinitialionmomentumatlowenergies(panelb). As a third option, the turning point of the ion trajectory may beassumedatthefootoftheimpactparameterontheincomingasymptote( ,whichwewillrefertoas“no ”or“TRIM”)asisdoneinmostMCcodesforamorphoustargets.Thisavoidsthediscrepancybetweenactual and chosen free flight path (ffp vs. ), again at the expense of constructingawrongoutgoingasymptote.Thelattermay seemasthe lesser of the two evils. However, consider the case when the ion trajectoryendsandtheimpactparameter andthedistance arelarge.Thenthemodelwhichputstheionattheintersectionoftheasymptotesterminatesthetrajectoryconsiderablybeforethetargetatom(2),whichisqualitativelycorrect.Incontrast,the“no ”modelputstheionclose toatom2,whichisunphysical. 3. Resultsanddiscussion 3.1. Bulksimulations In our recoil simulations we first only consider recoils starting atdepthsbetween30Å and60Å inordertoexcludesurfaceeffects.Afirst setofresultsforrecoilsstartingparalleltothesurfaceissummarizedinFig.3.InallpanelsthefirstmomentaccordingtoEq.(1)isshownasa function of the recoil energy. Results obtained with different BCAmodels(dashedlines)arecomparedtotheMDresults(blackdotswitherrorbarsconnectedbyblacklines).Differentdisplacementenergies areindicatedbydifferentcolors.Inthetoprowtheenergydependentmodel for the maximum impact parameter is used which guaranteesthat no energy transfers above are missed. With the default IMSIL model(leftcolumn)theMCsimulationscannotbefittedtotheMDdataby varying the displacement energy; the best result is obtained with eV. Oddly, the first moment can become negative for smalldisplacementenergies( eVand eV),i.e.,thesumoftheatom displacements is opposite to the direction of the initial recoilmomentum, which is clearly unphysical. This does not happen whentheactualfreeflightpathsare restrictedtonon-negativevalues asdescribed in Section 2 (middle column). Then a displacement energy of5eVyieldsanexcellentfittotheMDdata.WiththeTRIMmodel(rightcolumn)theproblemofnegativemomentsdoesnotoccureither,butanoptimum fit is now obtained with a displacement energy near 15eV. WiththeTRIMmodelthelongestfreeflightpathsarechosenamongthethreemodels,whichiscompensatedforbyasmallernumberofrecoilsresulting from a larger displacement energy. In the bottom row ofFig. 3, the corresponding MC results using a fixed maximum impactparameterof1.53Å areshown.NowtheMDdata arewellfittedwith displacementenergies near15eV. Fig.4showsthedisplacement oftheprimaryrecoilinadditiontothesumofalldisplacementsforthethreemodelswiththebestfitof 31 G.Hobler,etal. Nuclear Inst. and Methods in Physics Research B 447 (2019) 30–33 Fig.3.FirstmomentasafunctionoftherecoilenergyusingvariousBCAmodels.Toprow:usingtheenergydependentmaximumimpactparameter;bottomrow:usingafixedmaximumimpactparameterof1.53Å.Columns1–3:differenttrajectoryconstructionmodels,seeSection2.TheMDresultsrepresentedbytheblack symbolswitherrorbars,connectedbylines,arethesameinallpanels.Recoilsarestartedparalleltothesurface. the displacementenergy.Whileboththe modeland theTRIM model are in excellent agreement with the MD data taking all recoilsinto account, the TRIM model overestimates the displacement of theprimaryrecoil.Sincethedisplacementoftheprimaryrecoilisoneterminthesumdefiningthefirstmoment,Eq.(1),thisiscompensatedforbyfewer terms in the sum, which is accomplished by the higher dis-placement energy. Thus, the TRIM model partitions the first momentunphysically into contributions of the primary recoils and those ofhigher-order recoils. The model therefore appears to be preferable. 3.2. Theeffectofthesurface First moments resolved with respect to the starting depth of theprimary recoils are shown in Fig. 5 for a recoil energy of 100eV. Focussing first on recoils startingparallelto the surface (green lines),one observes a significant increase in the first moment towards thesurfaceincaseoftheMDresults,whilethefirstmomentcalculatedbyMC decreases towards the surface. The decrease in the MC results is easily explained: Some recoils escape into the vacuum and cannotproducehigher-orderrecoilsthere,thusreducingthenumberoftermsinEq.(1).TheincreaseintheMDresultsindicatesareducedresistance ofthesurfaceatomsagainstrelocation.Astheeffectisquitesubstantial,thisreducedresistancelikelyextendstoconsiderabledepth,indicatingcollectivemotionofatomsinanear-surfacelayerofthetarget. 32 G.Hobler,etal. Nuclear Inst. and Methods in Physics Research B 447 (2019) 30–33 Fig.5alsocontainssimulationresultsforrecoilsstartingat60°withrespecttothe axiseithertowardsthesurface(redlines) ortowards thebulk(bluelines).Atlargeenoughdepths,underisotropicconditionsthe calculated first moments should be of the results for theinitialdirectionparalleltothesurface.Theresultsfortheinclinedinitial directions have therefore been multiplied with a factor of two. Thegoodagreementoftheresultsforthethreeinitialconditionsdeeperinthesampleindicatestheabsenceof(artificial)anisotropies. The first moments calculated by MD for primary recoils directedtowardsthesurface(redlines)show astrongincreaseasafunctionofdecreasing initial depth, until close to the surface where they arestrongly lowered. A pronounced decrease towards the surface is alsoobservedintheMCresults.ThelargerdecreaseinboththeMDandMCresults, compared to when the primary recoils start parallel to thesurface, is due to the fact that recoils have a larger probability of escapingintothevacuumwhentheprimaryrecoilisdirectedtowardsthesurface. The strong increase in the MD result at intermediate depthsindicatesthatthecollectivemotionisfacilitatedwhen alargerpartof therecoilcascadeisnearthesurface. Whentheprimaryrecoilsaredirectedtowardsthebulk(bluelines),there is virtually noinfluence oftheprimaryrecoil depth on the first momentintheMCresults,indicatingthatveryfewrecoilsescapeintothevacuumevenwhentheprimaryrecoilstartsatthesurface.Thisisconfirmed,e.g.,byaMCsputteringyieldof forprimaryrecoils startingatadepthof2Å withacomponentintothebulk,comparedto for the same primary recoils starting parallel to the surface.TheMDresultsshowasubstantialincreasetowardsthesurfacewitha maximumatthesurfacethatexceedseventhemaximumfortherecoils directedtowardsthesurface. 4. Conclusions We have shown that atom redistribution as calculated by MC simulationsstronglydependsontheimplementationoftheBCA,i.e.,onthe displacement energy, the choice of the free flight path and themaximum impact parameter, and the construction of the trajectories.Totheknowledgeoftheauthors,inallpreviousworktargetedatpredicting spontaneous pattern formation using MC simulations, thesedetails have not completely been specified. This means that none ofthesepublicationsisreproducible. Forrecoilsstartingdeepinthebulk,wehaveobtainedbestresultsbylimitingthefreeflightpathstonon-negativevaluesandchoosingadisplacement energy of eV. This value is considerably higher thanthevalueusedin ourimpactsimulations[12].Fortheimpactsimulations, the choice of eV results in considerably underestimatedcontributionstothefirstmoment.The reasonisthesurface effectdescribedinSection3.2:Thetargetatomsaremuchmoremobilenearthesurface,probablyduetosomecollectivemotion.Modelingthissurfaceeffectisanimportantchallenge,ifrealisticMCsimulationsareaspired. Acknowledgments DMandZPgratefullyacknowledgefinancialsupportfromthePolishNational Science Centre, Grant No. 2015/19/B/ST4/01892. MD simulationswereperformedatthePLGridcomputerfacilities. References [1] N.Kalyanasundaram,M.Ghazisaeidi,J.B.Freund,H.T.Johnson,Singleimpactcraterfunctionsforionbombardmentofsilicon,Appl.Phys.Lett.92(13)(2008)131909 . [2] S.A.Norris,M.P.Brenner,M.J.Aziz,Fromcraterfunctionstopartialdifferentialequations:a newapproachtoionbombardmentinducednonequilibriumpatternformation,J.Phys:Condens.Matter21(2009)224017 . [3] M.Z.Hossain,K.Das,J.B.Freund,H.T.Johnson,Ionimpactcraterasymmetrydeterminessurfacerippleorientation,Appl.Phys.Lett.99(15)(2011)151913 . [4] S.A.Norris,J.Samela,L.Bukonte,M.Backman,F.Djurabekova,K.Nordlund, C.S.Madi,M.P.Brenner,M.J.Aziz,Moleculardynamicsofsingle-particleimpactspredictsphasediagramsforlargescalepatternformation,NatureCommun.2(2011)276. [5] P.D.Shipman,R.M.Bradley,Theoryofnanoscalepatternformationinducedbynormal-incidenceionbombardmentofbinarycompounds,Phys.Rev.B84(8)(2011)085420, https://doi.org/10.1103/PhysRevB.84.085420. [6] M.Castro,R.Gago,L.Vázquez,J.MuñozGarcía,R.Cuerno,Stress-inducedsolidflowdrivessurfacenanopatterningofsiliconbyion-beamirradiation,Phys.Rev.B86(21)(2012)214107 . [7] M.P.Harrison,R.M.Bradley,Craterfunctionapproachtoion-inducednanoscalepatternformation:cratersforflatsurfaces areinsufficient,Phys.Rev.B89(24) (2014)245401 . [8] J.MuñozGarcía,L.Vázquez,M.Castro,R.Gago,A.Redondo-Cubero,A.MorenoBarrado,R.Cuerno,Self-organizednanopatterningofsiliconsurfacesbyionbeamsputtering,Mater.Sci.Eng.R:Reports86(2014)1–44. [9] D.A.Pearson,R.M.Bradley,Theoryofterracedtopographiesproducedbyobliqueincidenceionbombardmentofsolidsurfaces,J.Phys.:Condens.Matter271(2015)015010, https://doi.org/10.1088/0953-8984/27/1/015010. [10] R.M.Bradley,H.Hofsäss,Nanoscalepatternsproducedbyself-sputteringofsolidsurfaces:Theeffectofionimplantation,J.Appl.Phys.120(7)(2016)074302 . [11] S.A.Norris,J.C.Perkinson,M.Mokhtarzadeh,E.Anzenberg,M.J.Aziz,K.F.Ludwig,DistinguishingphysicalmechanismsusingGISAXSexperimentsandlineartheory:theimportanceofhighwavenumbers,Sci.Rep.7(1)(2017)2016. [12] G.Hobler,D.Maciążek,Z.Postawa,Craterfunctionmoments:roleofimplantednoblegasatoms,Phys.Rev.B9715(2018)155307, https://doi.org/10.1103/ PhysRevB.97.155307. [13] R.M.Bradley,J.M.E.Harper,Theoryofrippletopographyinducedbyionbombardment,J.Vac.Sci.Technol.A6(4)(1988)2390–2395. [14] G.Carter,V.Vishnyakov,Rougheningandrippleinstabilitiesonion-bombardedSi,Phys.Rev.B54(24)(1996)17647–17653. [15] L.Bukonte,F.Djurabekova,J.Samela,K.Nordlund,S.Norris,M.Aziz,Comparisonofmoleculardynamicsandbinarycollisionapproximationsimulationsforatomdisplacementanalysis,Nucl.Instrum.Meth.Phys.Res.B297(2013)23–28,https://doi.org/10.1016/j.nimb.2012.12.014. [16] A.Lopez-Cazalilla,A.Ilinov,L.Bukonte,K.Nordlund,F.Djurabekova,S.Norris, J.C.Perkinson,Simulationofredistributiveanderosiveeffectsina-SiunderAr+irradiation,Nucl.Instrum.Meth.Phys.Res.B414(2018)133–140,https://doi.org/ 10.1016/j.nimb.2017.11.019. [17] H.Hofsäss,Surfaceinstabilityandpatternformationbyion-inducederosionand massredistribution,Appl.Phys.A114(2)(2014)401–422. [18] H.Hofsäss,O.Bobes,K.Zhang,Argonionbeaminducedsurfacepatternformation onSi,J.Appl.Phys.119(3)(2016)035302 . [19] G.Hobler,MonteCarlosimulationoftwo-dimensionalimplanteddopantdistributionsatmaskedges,Nucl.Instrum.Meth.Phys.Res.B96(1995)155. [20] G.Hobler,A.Simionescu,Accelerationofbinarycollisionsimulationsincrystallinetargetsusingcriticalanglesforionchanneling,Nucl.Instrum.Meth.Phys.Res.B102(1995)24–28,https://doi.org/10.1016/0168-583X(94)00792-6. [21] J.F.Ziegler,J.P.Biersack,U.Littmark,TheStoppingandRangeofIonsinSolids,PergamonPress,NewYork,1985. [22] W.Möller,TRI3dyn–Collisionalcomputersimulationofthedynamicevolutionof3-dimensionalnanostructuresunderionirradiation,Nucl.Instrum.Meth.Phys.Res.B322(2014)23–33, https://doi.org/10.1016/j.nimb.2013.12.027. [23] M.T.Robinson,I.M.Torrens,Computersimulationofatomic-displacementcascadesinsolidsinthebinary-collisionapproximation,Phys.Rev.B9(12)(1974)5008–5024. [24] W.Eckstein,ComputerSimulationofIon-solidInteractions,Springer,Berlin,1991. 33 Cite This: J.
Phys.
Chem.
C
2019,
123,
20188−20194
pubs.acs.org/JPCC
MD-Based Transport and Reaction Model for the Simulation of SIMS Depth Profiles of Molecular Targets Nunzio Tuccitto,†
Dawid Maciazek,‡
Zbigniew Postawa,‡
and Antonino Licciardello*,†
†Dipartimento di Scienze Chimiche, Universita di Catania, viale A. Doria, 6, 95125 Catania, Italy ‡Smoluchowski Institute of Physics, Jagiellonian University, 30-348 Krakow, Poland *Supporting
Information
S ABSTRACT: We present a novel theoretical approach allowing us to model erosion and chemical alteration of organic samples during depth profiling analysis by secondary-ion mass spectrometry with cluster projectile ion beams. This approach is able to take into account all of the cumulative phenomena occurring during such analysis, including ion-beam-induced reactions and atomic/molecular mixing by means of the numerical solution of an advection−diffusion−reaction (ADR) differential equation. The results from the singleimpact molecular dynamics computer simulations are used as a source for the input parameters into the ADR model. Such an approach is fast and allows turning the phenomenological model into a more quantitative tool capable of calculating molecular secondary-ion mass spectrometry depth profiles. The model is used to describe phenomena taking place during depth profiling of polystyrene samples by 20 keV C60,Ar872, and Ar1000 projectiles. It is shown that theoretical findings are in good agreement with the experimental results. The model is also used to determine the overall efficiency of nitrogen monoxide molecules in eliminating the radicals responsible for polystyrene cross-linking induced by analyzing ion beams. ■ INTRODUCTION In recent years, secondary-ion mass spectrometry (SIMS), a well-established technique for the characterization of solid surfaces and thin films of both inorganic and organic materials,1,2
has undergone a drastic evolution, related to the introduction of polyatomic primary ion beams,3
which boosted new perspectives in the analysis of organic materials and polymers. A big effort has been made to model the impact of primary cluster ions on the organic target by means of molecular dynamics (MD) simulations.4−9
These simulations are used to model the phenomena involved in an individual cluster-target impact that typically occurs over a time scale of the order of a few tens of picoseconds.10
In some aspects, MD simulations can be considered as the theoretical counterpart of static-SIMS experiments. In contrast, the theoretical modeling of the physicochemical phenomena occurring during dynamic-SIMS (D-SIMS) depth profiling of organic and polymer targets under cluster ion beams is less developed. In modeling D-SIMS experiments, time scales much larger than those involved in MD simulations must be considered to take into account slower and cumulative phenomena. Such phenomena, including ion-beam-induced reactions and atomic/molecular mixing, can substantially modify the nature of the target. Wucher et al.11,12
were the first to develop a model (statistical sputtering model, SSM) aimed to correlate the information obtained from molecular dynamics simulations to depth profiles. Additional models based on the SSM concept appeared later.4,5
However, all of these approaches were computationally intensive. Recently, we simulated D-SIMS experiments by means of the numerical solution of an advection−diffusion−reaction differential equation.13
Such a © 2019 American Chemical Society 20188 transport−reaction (TR) model takes into account the effects of ion-beam-induced mixing and reactivity, which occur during the depth profiling experiments. This model is able to simulate a complete D-SIMS experiment without requiring large computational capabilities. The values for the model input parameters can be estimated either from experimental data, as in our previous work, or from the output of computational simulations. Recently, it has been shown that the results of MD simulations can be used successfully for determining these parameters in the case of nonreactive inorganic samples.14
In this paper, we present the integration of MD simulations results into the transport−reaction model for depth profiling of polymer films, also in the case when ion-beam-induced chemical damage cannot be neglected. The MD simulation outputs are used as input parameters into the model. This allowed us to turn the previously reported phenomenological transport−reaction (TR) model into a quantitative model for the prediction of molecular D-SIMS depth profiles by the MD results. Henceforth, we call this model MD-TR. ■ EXPERIMENTAL METHODS The molecular dynamics (MD) computer simulations are used to model cluster bombardment of polystyrene, silicon, and diamond surfaces by 20 keV C60 at a 45° impact angle. Briefly, the motion of the particles is determined by integrating Hamilton’s equations of motion. The forces among carbon atoms in the system are described by the ReaxFF-lg force Received: February 20, 2019 Revised: June 19, 2019 Published: July 24, 2019 DOI: 10.1021/acs.jpcc.9b01653
J. Phys. Chem. C 2019, 123, 20188−20194 The Journal of Physical Chemistry C field,15
which allows for the creation and breaking of covalent bonds. This potential is splined at short distances with a ZBL potential16
to describe high-energy collisions accurately. The forces among Si−Si are described by the Tersoff-3 potential.17
The interactions between Ar atoms in the projectile and between the projectile atoms and all other atoms in the systems are represented by the Lennard-Jones potential splined at short distances with the KrC potential,18
which is better than ZBL to describe energetic collisions of noble gas atoms. Energetics of bond breaking, essential for our studies (C−Hor C−C bond cleavage), was calculated with the ReaxFF potential. The values were compared with the reference data to verify whether this potential gives correct predictions. The results are presented in the Supporting
Information.We observe some deviations between the bond energies predicted by ReaxFF and the reference data. However, the differences always remained below 10% of the reference data. A detailed description of the MD method can be found elsewhere.19
The shape and size of the samples are chosen based on visual observations of energy transfer pathways stimulated by impacts of C60 projectiles. As a result, hemispherical samples with diameters of 40, 32, and 26 nm for polystyrene, silicon, and diamond, respectively, are used. These samples contain 1 621 182, 432 449, and 823 719 atoms. The calculated densities of these samples are 1.097, 2.334, 3.54 g/cm3 and agree well with the experimental values of 0.96− 1.04, 2.33, and 3.5−3.53 g/cm3 for polystyrene, silicon, and diamond, respectively. Rigid and stochastic regions near the external boundaries of the sample with thicknesses of 0.7 and 2.0 nm, respectively, were used to simulate the thermal bath that keeps the sample at the required temperature to prevent reflection of pressure waves from the boundaries of the system and to maintain the shape of the sample.20,21
The simulations are run at the 0 K target temperature in an NVE ensemble and extend up to 50, 20, and 20 ps for polystyrene, silicon, and diamond, respectively, which is long enough to achieve saturation in the ejection yield vs time dependence. Impacts of 20 keV C60,Ar872, and Ar1000 were modeled at a 45° incidence angle relative to the surface normal to comply with the experimental conditions used in ref 22. The two different Ar clusters, differing about 10% in size, were chosen to get an idea on the possible effects of non-monodisperse size distribution of the clusters since in real experiments clusters exhibit a distribution around the nominal size. Simulations are performed with the Large-scale Atomic/Molecular Massively Parallel Simulator code,21
which was modified to describe sputtering conditions better. Numerical solution of the equation involved in the proposed model was performed by means of a Python-based script, developed on purpose. We simulated the D-SIMS depth profiles of 100 nm-thick PS thin films deposited onto a Si substrate. The primary ion beam current was 1 nA, rastered over 500 × 500 μm2. Experimental results have been replicated as in ref 22. ■ RESULTS AND DISCUSSION In the MD-TR model, the position of the bombarded surface is kept invariant so that the erosion process is represented as a “travel” of the underlying material toward the surface. The sputtering rate is represented by the travel velocity (v) of the inner target layers toward the surface. When the material moving toward the surface enters the layer altered by the ion beam, it is redistributed by ion-beam mixing and, also, ionbeam-induced chemistry is triggered. The partial differential equation describing the evolution of the concentration profile, C(z,t), of a certain species during a sputter-profile experiment can be written as ∂C ( DC ) −∇·( vC =∇·∇) +R ∂t (1) where C is the volumetric atomic concentration of the species composing the target material, expressed in atoms/nm3. If the sputtering area is kept constant during the simulations, only the traveling direction z (normal to the surface plane) must be considered; thus, the unit of C is atoms/nm. For the reconstruction of the depth profile, the model is using the so-called sampling depth, which takes into account the fact that atoms ejected by a single projectile impact can be initially located also below the surface. SIMS intensities (I(t)) are calculated by the following equation ∞ −z /( ) λ I() t ∝∫C(, zt )e ·d z (2) 0 where λ is the extent of the sampling depth distribution, which is known to be of the order of interatomic distances. The quantity (v) is equivalent to the erosion rate expressed in nm/s. If we express the flux (φ) of projectiles in ions nm−2 s−1 and the volumetric sputtering yield (Y)innm3/projectile, we can write ÄÉ Ä 3 ÉÄ É nm nm projectile v =Y ·φ 2 Çs ÖÇprojectile ÖÅ ÇÖ(3) Å Ñ Å Ñ nm s Ñ The dominant quality of the transport and reaction model is the capability to simulate depth profiles of reactive target layers, such as those constituted by organics or polymers. For example, during depth profiling of the organic target by monatomic beams or during the bombardment of polystyrenelike polymers by means of C60 primary beams, the intensity of fragments characteristic of the original polymer is readily lost.23
Even in the case of sputtering of thick samples by means of large argon clusters, damaging is not negligible so far.24
Focusing on the specific case of the depth profiling of a polystyrene layer deposited onto a silicon substrate, based on the simplified assumption that the instantaneous erosion rate is equal to the weighted mean between the erosion rates of the intact and fully damaged material, we can write v =vC +vC +vC PS PS c damage Si Si (4) where vPS, vc, and vSi are the erosion velocities of the pristine polystyrene, damaged material, and silicon substrate, respectively. CPS, Cdamage,and CSi are the normalized relative concentrations of target atoms that, inside the altered layer, pertain to undamaged polystyrene repeating units, to the damaged material, and to the silicon substrate, respectively. In the present application of the MD-TR model, the volumetric sputtering yields of YPS, Yc, and YSi are obtained by MD simulations performed on polystyrene, diamond (chosen as a representative of the fully damaged, carbonlike material formed upon beam irradiation), and silicon, bombarded with C60, Ar872,or Ar1000 projectiles. It should be noted that the choice of diamond as a representative of a fully damaged polymer is probably too extreme, while an amorphous carbon or a:C−H cross-linked network would be more realistic. However, the generation of a truly amorphous system for MD simulation is very difficult and, on the other hand, the good agreement of 20189 DOI: 10.1021/acs.jpcc.9b01653
J. Phys. Chem. C 2019, 123, 20188−20194 The Journal of Physical Chemistry C our MD-TR simulations with experimental results (see infra) indicates that the choice of diamond is acceptable. The volumetric sputtering yield values, calculated by MD simulations, are reported in Table
1. Table 1. Volumetric Sputtering Yield (nm3/ion) from MD Simulations C60 Ar872 Ar1000 PS 84 99 93 Si 6 2 0.8 diamond 1.9 0.011 0.011 The quantity D in eq
1
is a function that takes into account the ion-beam-induced mixing in the altered layer region. In the present simulation, surface roughening is neglected. The mixing term D in eq
1
is related to the ion-beam-induced motion of particles in the collision cascade, which can be described by analogy with the self-diffusivity process. Strictly speaking, we should also consider the lateral component to the diffusion. In simulating a depth profile of a laterally homogeneous system, however, we can neglect the lateral component and, consequently, the model is reduced to a unidimensional one along the traveling direction z (normal to the surface plane). By analogy with the unidimensional random walk model for diffusion, the units of D is nm2/s. To extract the per impact ion-beam-induced diffusion (D′) from the MD simulations, we divide each simulated sample into the equally spaced slices along the z-axis. For each slice, we calculated a sum of the total square displacements along the z-axis of all atoms located initially in a given slice multiplied by the volume occupied by the displaced atom. For a slice i of a width dz at depth zi, this can be expressed by the following equation N 2 Dz d zV /d z ′[]=∑i i jj j =1 (5) where Ni is the total number of atoms in the given slice, dzj is the total displacement of the jth atom along the z-axis, and Vj is the atomistic volume occupied by this atom calculated from a modeled sample density. Final diffusion can be obtained with a simple multiplication by the flux DDφ(6) =′ Figure
1
shows such distribution, obtained from simulations, after the impact of each primary ion (Ar1000,Ar872, and C60)on the three target materials (PS, Si, and diamond, respectively). Note that, in the case of Ar clusters, a mass spread of ±10% does not have a significant impact on the D′ distribution. The function R in eq
1
incorporates beam-induced reactions occurring in the region involved in the interaction with the beam. Since the aim of the model is that of simulating the molecular depth profile, we are interested in estimating the evolution, with projectile fluence, of the surface concentration of the undamaged material. In view of this, R can be interpreted as a numeric function representing the sink of the reacted portion (Creact) of the original material in the simulated volume at each simulation step Δt. It is worth noting that at each iteration step of the simulation (involving the solution of the transport and reaction differential equation) the concentration of the unreacted material changes due to the formation of damaged material so that R will change with ion fluence. We assume that the ion-beam-induced reactions occur in the same region where the beam-induced mixing is active. Thus, we write Dz () ΔC react Rz =−Dz z d Δt () ∫0 ∞()(7) () where Dz is a factor representing the normalized numeric ∫∞Dz z ()d 0 function D, obtained from MD simulations as described above. The negative sign accounts for the fact that R is a sink term, which decreases the amount of pristine material. To estimate the atomic fraction of material that underwent some modification, we count, at the end of the MD simulation, (i) the number of C−H bonds that were turned into C−C bonds (indicating the formation of cross-links during the time scale of the simulation) and (ii) the number of atoms that exhibit a number of bonds lower than those they had in the original polymer. The latter quantity provides an indication about the concentration of reactive species (such as radicals) that, in turn, can evolve in damaged material on a time scale longer than that considered in the MD simulation. Table
2
reports the Table 2. Atomic Portion of the Reacted Material per Single Ion Beam Impact C60 Ar872 Ar1000 freeHatoms 24 3 0 reactive carbon atomsa
675 180 104 C−H bonds converted to C−C200 Creact 701 183 104 aC atoms that at the end of MD simulations display less bonds than the pristine material. atomic portion of the reacted material per single ion beam impact, in the case of C60,Ar872, and Ar1000 projectiles. Results indicate, in agreement with experimental results,25−27
that C60 induces much greater damage into polystyrene compared to Ar clusters. We observe that the amount of atoms participating to newly formed C−C bonds or that remain reactive (carbon Figure 1. Distribution along the depth of ion-beam-induced diffusion stimulated by the impact of each primary ion (C60,Ar872, and Ar1000) in (a) silicon, (b) polystyrene, and (c) diamond targets. Note that, according to MD results, the considered Ar clusters do not modify the diamond target. 20190 DOI: 10.1021/acs.jpcc.9b01653
J. Phys. Chem. C 2019, 123, 20188−20194 The Journal of Physical Chemistry C Article Figure 2. Simulated depth profiles of polystyrene on silicon obtained by Ar872,Ar1000, and C60 primary ions as a function of ion fluence (a) and depth (b). atoms with “missing” bonds or free H atoms) represents a small fraction (0.512% for C60, 0.057% for Ar872, and 0.036% for Ar1000) of the atoms that experienced the primary ion impact (please refer to Figure
1) and remained in the target at the end of the MD simulation (i.e., that were not sputtered away). However, as we will see in the following (Figure
2), the 1 order of magnitude difference between C60 and argon clusters can account for the different behavior of the two kinds of a projectile in terms of damage since the reactive species (most of them, presumably, radicals) can produce cross-linking in polystyrene. This is consistent also with previous findings, obtained in a different experimental context, showing that a few ion-beam-induced cross-links can induce an irreversible sol−gel transition in polystyrene.28−30
As it is evident from the simulated profiles reported in Figure
2, the model is able to reproduce the essential characteristics of the experimental profiles,25
namely, the initial drop of molecular signal intensity, which is very pronounced in the case of C60 and much lighter in the case of the largest argon cluster, as well as the strong differences in the primary ion dose needed for reaching the interface. Table
3
reports the average volumetric sputtering yields for a 100 nm-thick PS film, obtained by means of the MD-TR Table 3. Volumetric Sputtering Yields (Average) and Depth Resolution Calculated by the MD-TR Model for a 100 nm- Thick PS Film C60 Ar872 Ar1000 Y (nm3/ion) 2 77 87 depth
resolution
(nm)a
118 16 19 aDepth resolution was calculated by the usual31
16−84% intensity method. model using the rise of the substrate signal. Considering that the single-impact sputtering yields (as obtained from MD simulations, see Table
1) are rather similar for C60 and the two considered Ar clusters, it is clear that the much smaller average sputtering yield obtained in the case of C60 is due to the damage produced in the target by C60 ions. The value of 2 nm3/ion obtained by the MD-TR model for polystyrene sputtering with C60 primary ions is in reasonable agreement with the experimental finding of 0.5 nm3/ion for PS25
and <0.2 nm3/ion for a similar system.26
In other words, the MD-TR model is able to fill the gap between the experimentally determined yields and those calculated by MD simulations, thanks to the introduction of the reactivity term. Interestingly, also in the case of Ar clusters, the MD-TR model gives rise to lower average sputtering yields (ca. 30 and 7% lower than those reported in Table
1
for Ar872 and Ar1000, respectively), but in reasonable agreement with the figures (107 nm3/ion for Ar872 and 94 nm3/ion for Ar1000) obtained by extrapolation, according to the equation proposed by Seah,32
of experimental data obtained by Rading et al. in slightly different conditions.25
Also, the predicted trend of reduction of sputtering yield (larger for the smaller cluster) is in agreement with the expected increase of damage at higher energy per component atom of the cluster. Moreover, we observe that a deterioration of the depth resolution at the interface accompanies the case (C60 projectile) where large damage accumulation occurs, as observed in systems with PS-like behavior.26
According to Table
2, the main contribution to the ionbeam-induced damage is the production of reactive species such as H atoms or C species with “dangling bonds” that may undergo successive reactions. In the case of polystyrene, these reactions produce an increasingly cross-linked material. Many of these reactive species can be regarded as radicals. Recently, some of us demonstrated that nitrogen monoxide (NO), a well-known radical scavenger, is able to reduce strongly the damage of PS-like polymers in C60-SIMS depth profiling experiments.22,26
Since the MD-TR model can include any chemical reaction occurring inside the altered layer, we included NO radical scavenging. We assumed that the reactive radical species R• produced during a cluster impact, whose amount is estimated from the MD simulation (see Table
2), are quenched by the reaction with NO molecules with the formation of stable, unreactive species R +NO →R −NO accordingly to a 1:1 stoichiometry. Thus, eq
5
becomes Dz () Δ( C −C ) react NO Rz =− () ∫∞Dz z Δt ()d 0 (8) 20191 DOI: 10.1021/acs.jpcc.9b01653
J. Phys. Chem. C 2019, 123, 20188−20194 The Journal of Physical Chemistry C Figure 3. Simulated C60-SIMS depth profiles of a 100 nm-thick film of polystyrene on silicon. Polystyrene (a) and silicon (b) intensities were calculated at different levels of CNO/Creact ratios. Figure 4. (a) PS sputtering yield (20 keV C60 projectile) as a function of CNO/Creact. The top X-axis shows the NO partial pressure computed as described in the main text. (b) Linear regression of the theoretical to experimental NO partial pressure. where CNO represents the number of nitric oxide molecules reaching the surface during the simulation step Δt. Figure
3
reports the simulated depth profiles of silicon and polystyrene in the case of C60 projectiles at various CNO/Creact ratios. For the sake of simplicity, we assumed that (i) the probability that a NO molecule reacts with a beam-produced radical is unitary and (ii) such probability is independent of the location of the radical inside the (remaining) modified material. Also, we assume (iii) that the quenched macromolecular radicals will continue to behave, under ion impacts, in the same way as the original polymer.13
In agreement with the data obtained in NO-assisted C60 depth profiling of PS and PS-like polymers,22,26
simulation shows that the higher the NO amount, the faster the erosion process and the better the shape of the molecular depth profile of PS. To correlate theoretical findings with experimental results, we estimated the number of NO molecules hitting the surface per time unit by means of impingement rate at several NO pressures molecules NP A = 2 cms 2 π MRT (9) where NA is Avogadro’s number, P is the NO partial pressure in the analysis chamber, M is the NO molar mass, R is the ideal gas constant, and T is the temperature. 20192 Figure
4a shows the trend of PS sputtering yield as a function of the NO/radical ratio as obtained from the simulation. In the same figure, the top X-axis shows the corresponding NO partial pressures, computed according to eq
9.In Figure
4b, we report the correlation between theoretical and experimental pressures needed for obtaining the same sputtering yield volume. In the range of NO pressures considered, correlation is linear, with a slope of about 0.1. This means that to obtain a certain effect on the sputtering yield the experimental NO pressure must be increased about 10 times more than expected from the simulation. We believe that this discrepancy is related to the approximation that NO molecules can access quantitatively all of the radicals produced in the altered layer so that the predicted overall efficiency of NO in “killing” the radicals is overestimated (unity instead of ca. 0.1). Also, we must note that the linear correlation between calculated and experimental NO pressures is verified when significant amounts of NO (CNO/Creact > 0.8) are considered. Indeed, the nonzero value of intercept in Figure
4b has no physical meaning and it indicates that the linear behavior cannot be extrapolated at very low NO pressures. Clearly, although the qualitative agreement between the previsions and the experiment is promising, additional elaboration and refinements are needed on this point. In conclusion, the integration of the results of molecular dynamics simulations in a transport/reaction model was used DOI: 10.1021/acs.jpcc.9b01653
J. Phys. Chem. C 2019, 123, 20188−20194 The Journal of Physical Chemistry C for the prediction of the SIMS depth profile. In particular, the model allows the simulation of the SIMS depth profiles of organic samples under cluster primary beam irradiation and is able to take into account ion-beam-induced reactions and the effect of reactive gas dosing. *Supporting Information ■ ASSOCIATED CONTENT S The Supporting Information is available free of charge on the ACS
Publications
website
at DOI: 10.1021/acs.jpcc.9b01653. List and energetics of the selected chemical reactions essential for our studies leading to C−Hor C−C bond cleavage calculated with the ReaxFF potential (PDF) ■ AUTHOR INFORMATION Corresponding Author *E-mail: alicciardello@unict.it. ORCID Nunzio Tuccitto: 0000-0003-4129-0406
Zbigniew Postawa: 0000-0002-7643-5911
Antonino Licciardello: 0000-0001-5146-8971
Notes The authors declare no competing financial interest. ■ ACKNOWLEDGMENTS N.T. and A.L. acknowledge “Piano della Ricerca di Ateneo” from the University of Catania for financial support. Z.P. and D.M. would like to gratefully acknowledge the financial support from the Polish National Science Centre, Programme 2015/19/B/ST4/01892. Molecular dynamics computer simulations were performed at the PLGrid Infrastructure. ■ REFERENCES (1) Benninghoven, A. Chemical Analysis of Inorganic and Organic Surfaces and Thin Films by Static Time-of-Flight Secondary Ion Mass Spectrometry (TOF-SIMS). Angew. Chem., Int. Ed. 1994, 33, 1023− 1043. (2) Giamblanco, N.; Zhavnerko, G.; Tuccitto, N.; Licciardello, A.; Marletta, G. Coadsorption-dependent Orientation of Fibronectin Epitopes at Hydrophilic Gold Surfaces. Soft Matter 2012, 8, 8370− 8378. (3) Gillen, G.; Roberson, S. Preliminary Evaluation Of An SF5 + Polyatomic Primary Ion Beam For Analysis of Organic Thin Films by Secondary Ion Mass Spectrometry. Rapid Commun. Mass Spectrom. 1998, 12, 1303−1312. (4) Paruch, R. J.; Garrison, B. J.; Postawa, Z. Partnering Analytic Models And Dynamic Secondary Ion Mass Spectrometry Simulations to Interpret Depth Profiles due to Kiloelectronvolt Cluster Bombardment. Anal. Chem. 2012, 84, 3010−3016. (5) Paruch, R. J.; Garrison, B. J.; Postawa, Z. Mixed MD Simulation -Analytical Model Analysis of Ag(111), C60 Repetitive Bombardment in The Context of Depth Profiling for Dynamic SIMS. Surf. Interface Anal. 2013, 45, 154−157. (6) Garrison, B. J.; Schiffer, Z. J.; Kennedy, P. E.; Postawa, Z. Modeling Dynamic Cluster SIMS Experiments. Surf. Interface Anal. 2013, 45,14−17. (7) Czerwinski, B.; Delcorte, A. Molecular Dynamics Study of Fullerite Cross-Linking under keV C60 and Arn Cluster Bombardment. J. Phys. Chem. C 2013, 117, 3595−3604. (8) Czerwinski, B.; Postawa, Z.; Garrison, B. J.; Delcorte, A. Molecular Dynamics Study of Polystyrene Bond-Breaking and Crosslinking under C60 and Arn Cluster Bombardment. Nucl. Instrum. Methods Phys. Res., Sect. B 2013, 303,23−27. (9) Czerwinski, B.; Delcorte, A. Chemistry and Sputtering Induced by Fullerene and Argon Clusters in Carbon-Based Materials. Surf. Interface Anal. 2014, 46,11−14. (10) Delcorte, A.; Vanden Eynde, X.; Bertrand, P.; Vickerman, J. C.; Garrison, B. J. Kiloelectronvolt Particle-Induced Emission And Fragmentation Of Polystyrene Molecules Adsorbed On Silver: Insights From Molecular Dynamics. J. Phys. Chem. B 2000, 104, 2673−2691. (11) Wucher, A.; Krantzman, K. D. A Statistical Approach to Delta Layer Depth Profiling. Surf. Interface Anal. 2012, 44, 1243−1248. (12) Wucher, A.; Krantzman, K. D.; Lu, C.; Winograd, N. A Statistical Interpretation of Molecular Delta Layer Depth Profiles. Surf. Interface Anal. 2013, 45,39−41. (13) Tuccitto, N.; Zappala, G.; Vitale, S.; Torrisi, A.; Licciardello, A. A Transport And Reaction Model For Simulating Cluster Secondary Ion Mass Spectrometry Depth Profiles Of Organic Solids. J. Phys. Chem. C 2016, 120, 9263−9269. (14) Maciazek, D.; Paruch, R.; Postawa, Z.; Garrison, B. J. Micro-and Macroscopic Modeling of Sputter Depth Profiling. J. Phys. Chem. C 2016, 120, 25473−25480. (15) Liu, L. C.; Liu, Y.; Zybin, S. V.; Sun, H.; Goddard, W. A. ReaxFF-Lg: Correction of the ReaxFF Reactive Force Field for London Dispersion, with Applications to the Equations of State for Energetic Materials. J. Phys. Chem. A 2011, 115, 11016−11022. (16) Ziegler, J. F.; Biersack, J. P.; Littmark, U. The Stopping and Range of Ions in Matter; Pergamon: New York, 1985. (17) Tersoff, J. Modeling Solid-State Chemistry -Interatomic Potentials for Multicomponent Systems. Phys. Rev. B 1989, 39, 5566− 5568. (18) Wilson, W. D.; Haggmark, L. G.; Biersack, J. P. Calculations of Nuclear Stopping, Ranges, and Straggling in Low-Energy Region. Phys. Rev. B 1977, 15, 2458−2468. (19) Garrison, B. J.; Postawa, Z. Computational View of Surface Based Organic Mass Spectrometry. Mass Spectrom. Rev. 2008, 27, 289−315. (20) Postawa, Z.; Czerwinski, B.; Szewczyk, M.; Smiley, E. J.; Winograd, N.; Garrison, B. J. Enhancement of Sputtering Yields Due to C60 Versus Ga Bombardment of Ag{111} as Explored by Molecular Dynamics Simulations. Anal. Chem. 2003, 75, 4402−4407. (21) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular-Dynamics. J. Comput. Phys. 1995, 117,1−19. (22) Havelund, R.; Licciardello, A.; Bailey, J.; Tuccitto, N.; Sapuppo, D.; Gilmore, I. S.; Sharp, J. S.; Lee, J. L. S.; Mouhib, T.; Delcorte, A. Improving Secondary Ion Mass Spectrometry C60 N+ Sputter Depth Profiling Of Challenging Polymers With Nitric Oxide Gas Dosing. Anal. Chem. 2013, 85, 5064−5070. (23) Möllers, R.; Tuccitto, N.; Torrisi, V.; Niehuis, E.; Licciardello, A. Chemical Effects in C60 Irradiation of Polymers. Appl. Surf. Sci. 2006, 252, 6509−6512. (24) Yamamoto, Y.; Ichiki, K.; Seki, T.; Aoki, T.; Matsuo, J. Ion-Induced Damage Evaluation with Ar Cluster Ion Beams. Surf. Interface Anal. 2013, 45, 167−170. (25) Rading, D.; Moellers, R.; Cramer, H.-G.; Niehuis, E. Dual beam depth profiling of polymer materials: comparison of C60 and Ar cluster ion beams for sputtering. Surf. Interface Anal. 2013, 45, 171− 174. (26) Zappala, G.; Motta, V.; Tuccitto, N.; Vitale, S.; Torrisi, A.; Licciardello, A. Nitric Oxide Assisted C60 Secondary Ion Mass Spectrometry For Molecular Depth Profiling of Polyelectrolyte Multilayers. Rapid Commun. Mass Spectrom. 2015, 29, 2204−2210. (27) Mahoney, C. M. Cluster secondary ion mass spectrometry of polymers and related materials. Mass Spectrom. Rev. 2010, 29, 247− 293. (28) Puglisi, O.; Licciardello, A.; Calcagno, L.; Foti, G. Molecular Weight Distribution and Solubility Changes in Ion-Bombarded Polystyrene. Nucl. Instrum. Methods Phys. Res., Sect. B 1987, 19−20, 865−871. 20193 DOI: 10.1021/acs.jpcc.9b01653
J. Phys. Chem. C 2019, 123, 20188−20194 The Journal of Physical Chemistry C (29) Licciardello, A.; Puglisi, O.; Calcagno, L.; Foti, G. Crosslinking Yield in Particle Beam Irradiated Polystyrene. Nucl. Instrum. Methods Phys. Res., Sect. B 1990, 46, 338−341. (30) Licciardello, A.; Puglisi, O.; Calcagno, L.; Foti, G. UV Absorption and Sol-Gel Transition in Ion-Bombarded Polystyrene. Nucl. Instrum. Methods Phys. Res., Sect. B 1989, 39, 769−772. (31) ASTM E-42, Standard Terminology Relating to Surface Analysis E 673-91c; ASTM: Philadelphia, 1992. (32) Seah, M. P. Universal Equation for Argon Gas Cluster Sputtering Yields. J. Phys. Chem. C 2013, 117, 12622−12632. 20194 DOI: 10.1021/acs.jpcc.9b01653
J. Phys. Chem. C 2019, 123, 20188−20194 pubs.acs.org/ac
Intuitive Model of Surface Modification Induced by Cluster Ion Beams Dawid
Macia̧zek,
Micha
Kanski,
and
Zbigniew
Postawa*
Cite This: Anal. Chem. 2020, 92, 7349−7353 Read
Online
Metrics
&
More
Article
Recommendations
ACCESS ABSTRACT: Topography development is one of the main factors limiting the quality of depth profiles during depth profiling experiments. One possible source of topography development is the formation of self-organized patterns due to cluster ion beam irradiation. In this work, we propose a simple model that can intuitively explain this phenomenon in terms of impact-induced mass transfer. By coupling our model with molecular dynamics simulations, we can predict the critical incidence angle, which separates the smoothening and roughening regimes. The results are in quantitative agreement with experiments. It is observed that the problems arising from topography development during depth profiling with cluster projectiles can be mitigated by reducing the beam incidence angle with respect to the surface normal or increasing its kinetic energy. S urface-sensitive techniques such as Secondary Ion Mass Spectrometry (SIMS), X-ray Photoelectron Spectroscopy (XPS), and Auger Electron Spectroscopy (AES) combined with ion-induced material removal can be used to create spatial maps of the chemical composition of materials in a process called depth profiling. This approach has been successfully applied to many systems both inorganic and organic.1−3
The main problem of depth profile acquisition is the degradation of resolution with a time of analysis.4
There are two main factors in play here: topography development5,6
and ion-beaminduced material alteration.4,7,8
The latter issue is virtually solved by the introduction of large gas cluster ion beams (GCIB).9−11
The topography development problem can be mitigated by sample rotation.5
However, this approach limits analysis to depth dimension only because lateral spatial information is averaged out. One possible source of topography development during ion irradiation is so-called spontaneous pattern formation where self-organized nanoscale ripples appear on the surface of the sample during the bombardment.12,13
Over the years, there has been a substantial theoretical effort to develop models capable of qualitative and quantitative prediction of the surface evolution.12−21
However, all available models have two main drawbacks concerning depth profiling with cluster projectiles. One is that none of these models have been applied to this type of projectiles or organic samples. All theoretical and almost all experimental research has been done with monatomic projectiles and inorganic samples. Only few experiments with cluster projectiles have been performed so far.22−25
Another problem is that proposed theoretical models have a complicated mathematical formulation, which makes it difficult to develop adequate physical intuition regarding the phenomenon in question. Without proper understanding, countering the emergence of undesirable surface roughness during depth profiling is difficult. In this manuscript, a simple model, based on a concept of mass transfer, is proposed to predict the conditions favorable for the formation of ripples on the solid surfaces bombarded by cluster projectiles. The unique beauty of this model is that, despite its simplicity, it can predict these conditions with good accuracy. Furthermore, it makes it easy to comprehend why surface roughness, which we equate to ripple formation, will increase or decrease under given experimental conditions. ■ SIMULATION DETAILS The molecular dynamics computer simulations are used to model the effect of argon cluster bombardment of gold and silicon samples. General information about MD simulations can be found elsewhere.26
Briefly, the motion of particles is determined by integrating Hamilton’s equations of motion. Forces among particles are described by following potentials: the Lennard-Jones potential splined with KrC to accurately Received: March 20, 2020 Accepted: April 21, 2020 Published: April 21, 2020 © 2020 American Chemical Society https://dx.doi.org/10.1021/acs.analchem.0c01219
7349 Anal. Chem. 2020, 92, 7349−7353 Analytical Chemistry pubs.acs.org/ac
Figure 1. (a) Idealized effect of the incidence angle θ on the mass transfer function M(θ) for a bombardment by cluster projectiles and schematic visualizations of the effect of mass transfer on the sample topography in regions where M(θ) increases (b) or decreases (c) with θ. Green and red quarters represent the global mass inflow and outflow for a given node, respectively. Black tilted arrows represent directions of the impacting projectiles, and the symbols, θi, θj, θk, θp, describe local incidence angles. The global incidence angles are θv and θw. describe high energy collisions is used to describe Ar−Ar interactions.27
The Tersoff potential splined at a short distance with ZBL potential28
is used to describe forces between Si atoms,29
while EAM force field is used for Au atoms.30,31
Finally, the ZBL potential is used to represent collisions between Ar−Au and Ar−Si atoms.28
The simulations are performed with the LAMMPS code.32
The values of the mass transfer function M(theta) for different projectile parameters (kinetic energy, incidence angle) are extracted from MD trajectories. The divide and conquer approach developed for modeling depth profiling is used. This approach has been described in detail previously.33,34
Briefly, a master sample is prepared as a cuboid with periodic boundary conditions imposed in x and y directions. The size of the master sample for silicon and gold equals to 40 × 40 × 30 nm and 42 × 42 × 32 nm, respectively. These samples contain 2.45 and 3.4 million atoms, respectively. Then, an impact point on the surface is selected randomly. A hemispherical region with a radius of 19 nm centered at this impact point is subsequently cut out from the master sample and used to simulate an impact of an Ar3000 cluster with a given kinetic energy and incident angle for 25 ps. The sample is subsequently quenched for 10 ps in order to maintain the desired temperature; it was equal to 0 K in the present study. After each simulation, all sputtered atoms are removed, the mass transfer function is calculated, and the remaining (nonsputtered) atoms are reinserted into the master sample. Subsequently, a new point of impact is randomly selected, and the cycle is repeated. For each combination of the kinetic energy and incident angle, a series of 50 simulations are performed, which corresponds to a fluence of approximately 3 × 1012 impacts/cm2, and the final value of the mass transfer function is calculated as an average. Four test studies corresponding to the incidence angles of 30°,35°,45°, and 50° are performed on a silicon sample with a 10 times larger number of impacts to probe the effect of the projectile fluence on the mass transfer function dependence on the angle of incidence. While the amplitude of this dependence varies with the number of impacts, its shape (position of maximum) remained the same. Repetitive bombardment setup has been chosen to avoid strain introduced by a single impact. The effect of projectile parameters (kinetic energy, incidence angle) on mass transfer can be obtained from molecular dynamics simulations. The mass transfer function is defined as M =∑Vdx i i (1) where V is the volume occupied by an individual atom, dxi is the displacement of the ith atom in the x-direction caused by the projectile impact, where the x-direction is the azimuthal direction of the incoming projectile, as proposed in ref 19. The ion beam density at the bombarded surface decreases with the angle of incidence θ due to the increase of the ion beam spot size on the irradiated surface. To account for this phenomenon, the results from the MD simulations should be normalized by the following formula: M () =M ()θ cos() θθ (2) MD where M(θ) is the angle-dependent mass transfer function, and MMD(θ) is the mass transfer calculated from the molecular dynamics simulations. ■ RESULTS AND DISCUSSION The process of sputtering is believed to have a minor effect on the ripple formation during cluster projectile bombard 35,36
ment.It is postulated that the mass transfer stimulated by impact of a projectile is of critical importance in this 35−37
case. In order to understand its influence, we begin by posing a question: what effect will a mass transfer near a surface have on its topography? To answer this question, we introduce a simple two-dimensional model, where the surface is represented by a set of discrete nodes, as shown in Figure
1. The coordinate system in this model is selected in such a way that the azimuthal direction of the ion beam is the x-axis, as depicted in Figure
1. Incoming projectiles are directed toward the macroscopic surface at the same global angle of incidence. This angle is equivalent to the experimental incidence angle. However, the global incidence angle usually does not represent the actual angle of incidence at a given node, because the real surface is never flat. In real samples, the local incidence angle should be defined relative to the surface normal at the point of projectile impact. It depends on the global angle of incidence 7350 https://dx.doi.org/10.1021/acs.analchem.0c01219
Anal. Chem. 2020, 92, 7349−7353 Analytical Chemistry pubs.acs.org/ac
Figure 2. Angle-dependent mass transfer functions for (a) silicon and (b) gold samples sputtered by 30 keV Ar3000. The AFM images show experimentally measured topography of the sample23,24
at respective angles. The dashed green line depicts a location where the smoothing/ roughening transition should be observed under reasoning presented in the text. The AFM images are reproduced with permission from refs 23
and 24
and are licensed under a Creative Commons Attribution (CC BY) license. and the local surface inclination. The magnitude of these effects is proportional to the inclination of the M(θ) function for a given global incidence angle, so, for the cases where a plateau is reached, the mass redistribution will not influence surface topography. After the projectile impact, the matter is transferred from one node to another. Our model is based on global mass transfer along the x-axis. For a cluster projectile impact at a flat surface along the surface normal, the same amount of material is transferred, on average, to the left and right neighbors of the bombarded node.35,36
Therefore, from the point of view of a global mass transfer along the x-axis, the resulting mass transfer is zero. For the off-normal angle of incidence, the mass is transferred along an azimuth of the incoming beam.36
The global mass transfer is no longer zero along the x-axis. The amount of mass transferred from a bombarded node increases initially with an increase of global incidence angle. However, the mass transfer must drop eventually, because no mass is transferred for the 90° incidence angle since the projectile never hits the sample. Based on these considerations, the value of the mass transfer function at a given node M(θ) should depend on the local angle of incidence θ, as shown in Figure
1a. There is a specific angle θc, which will be called a critical angle, where M(θ) has a maximum. This angle separates two regions, where M(θ) increases and decreases with the incidence angle. Three possible surface morphologies should be considered to investigate the temporal evolution of the bombarded surface. The surface can be flat, convex (hill), or concave (hole), near a point of projectile impact. An example of a surface exhibiting all these cases is shown in Figure1b,c. The total amount of material transferred into and outside a given node is depicted as radii of the quarter-circles. Green and red quarters represent the mass inflow and outflow for a given node, respectively. For the off-normal incidence, the mass is transferred in the x-direction. In our model, this means that for a given node, the mass inflow occurs only from a node on the left, while the material is moved to the node on the right. A flat section of a surface represented by node A in Figure
1b,c is the simplest situation to discuss. In this case, the local angle of incidence at nodes A-1 and A is identical. As a result, the amount of material incoming and outcoming from node A is the same, as represented by quarters of the equal radii. Consequently, the height of node A is not affected relative to its neighbors, and the morphology does not change at this point. The situation is different for convex and concave surfaces represented by nodes B and C, respectively. In the case of a convex surface, the local incidence angle θj or θp at node B is larger than the local incidence angle θi or θk at node B-1. The opposite situation occurs for node C. The final mass balance at nodes B and C depends on the shape of the M(θ) function. In the case where, for a given global incidence angle, M(θ) increases with θ, as depicted in Figure
1b, the mass inflow from the node B-1 is smaller than the outflow from node B, because θj > θi, therefore, M(θj)> M(θi). The mass is removed, and the height of node B decreases, as indicated by downward-pointing vertical arrows. For a concave surface (node C), the situation is the opposite. More mass is transferred into node C than removed from it. As a result, the depth of this depression decreases. Considering both these effects, the global surface roughness will decrease for the impacts presented in Figure
1b. Similar reasoning can be applied for the impact conditions where the M(θ) function decreases with θ. This case is depicted in Figure
1c. However, now the conclusions will be opposite to the situation shown in Figure
1b. While the relationship between θp and θk remains the same, that is, θp > θk, now M(θp)< M(θk). Therefore, the global transfer to node Bin Figure
1c is positive, while the mass is removed from node C. As a result, the elevation of the hill increases while the depression becomes more profound, which leads to an increase of the surface roughness. It is evident that the incidence angle corresponding to the maximum of the M(θ) function separates the surface smoothing and roughening regimes. Experimental results describing the effect of the incidence angle on the ripple formation during the Ar3000 bombardment of silicon and gold surfaces23,24
are used to verify the predictions of the proposed model. The shapes of the mass transfer functions calculated from molecular dynamics (MD) computer simulations are shown in Figure
2. The AFM images presenting experimentally measured topography of the sample at respective global incidence angles are shown below the graphs of individual M(θ) function. The experimental data show that the ripples begin to form above 40° and 30° at silicon and gold surfaces, respectively. The critical angles obtained from the calculated M(θ) functions for the same systems are 45° and 40°, respectively. It is evident that the transition from a smooth to a rough surface correlates very well with the value of the critical angle. This agreement proves that our model, regardless of its simplicity, is working correctly and 7351 https://dx.doi.org/10.1021/acs.analchem.0c01219
Anal. Chem. 2020, 92, 7349−7353 Analytical Chemistry pubs.acs.org/ac
corroborates the hypothesis that the mass transfer is predominantly responsible for an angle-dependent roughening/smoothing phenomenon during cluster projectile bombardment. So far, only the effect of the incidence angle on the surface roughness was discussed. However, the kinetic energy of cluster projectiles may also influence the surface roughness. In fact, it has been shown that, at very low energies per atom, ripples appear in organic samples bombarded by gas cluster projectiles.25
This process leads to fast degradation of the depth resolution with depth.6,25
The problem was eliminated by increasing the projectile kinetic energy (energy per projectile atom). Unfortunately, the RMS roughness has not been directly measured in that paper.25
Nevertheless, is has been established that the buildup of the surface roughness is the main factor limiting the achievable depth resolution.6,11
Therefore, a possibility to achieve a better depth resolution indicates that the surface roughness must decrease with the increase of the kinetic energy per atom. Angle-dependent mass transfer functions calculated for several kinetic energies of Ar3000 projectiles, bombarding the silicon surface, are shown in Figure
3
to probe the effect of the Figure 3. Angle-dependent mass transfer functions for silicon, bombarded by Ar3000 with different kinetic energies. For each incidence energy, the location of the maximum is highlighted by an arrow. primary kinetic energy. Silicon is very different from the systems investigated in ref 25
where organic multilayers were analyzed. Therefore, no quantitative agreement between theoretical and experimental data can be expected. Nevertheless, similar trends should be observed in both these studies. It is evident that the value of the critical angle, which separates smoothing from roughening conditions, indeed shifts toward larger incidence angles with the increase of the projectile kinetic energy. It is expected, therefore, that at low projectile kinetic energy, the critical angle is below the incidence angle used in the experiment (45°), and the surface will be roughened. However, when kinetic energy is increased, the value of the critical angle θc shifts toward larger values, leading to a decrease of roughness due to the elimination of ripple formation. ■ CONCLUSION We have shown that ripple formation during cluster ion beam irradiation, which we equate to surface roughening, can be explained by impact-induced mass transfer with a simple and intuitive model. The main conclusion from our model is that the transition between smoothing and roughening regimes can be explained by the location of a maximum, called the critical angle, of an angle-dependent mass transfer function. If the incidence angle is lower than the critical angle, the mass transfer will have a smoothening effect, and if it is higher, the roughness will increase. Due to the nature of this function, the location of the maximum is expected to be around 45°. This finding has substantial practical importance because the incidence angle of a sputtering beam is 45° in a majority of commercially available apparatuses. We suggest that shifting the incidence angle closer to the surface normal by even a few degrees will have a beneficial impact on the quality of acquired depth profiles. Furthermore, we have shown that the increase of the cluster impact energy shifts the critical angle toward larger values. This observation suggests that the roughening problem, for some cases, can be resolved by increasing the kinetic energy of the beam. Finally, it should be mentioned that the proposed model works both for inorganic and organic samples. In this paper, we have focused on inorganic materials because results discussing the effect of the angle of incidence on the ripple formation are only available in the literature for such samples. We need such results to validate the predictions of our model. However, it should also be emphasized that our model can only be used to predict the onset of ripple formation during surface bombardment by cluster projectiles. Both the experimental results and computer simulations show that for these projectiles process of mass transfer prevails over sputtering.23,24,34,36
As a result, sputtering can be ignored as it was done in our model. A similar approach cannot be applied to atomic missiles for which mass transport is much smaller due to the significantly smaller projectile size and momentum. Both the existing theories and experimental data indicate that the effect of sputtering cannot be ignored for these projectiles.12−21
■ AUTHOR INFORMATION Corresponding Author Zbigniew Postawa − Smoluchowski Institute of Physics, Jagiellonian University, 30-348 Krakow, Poland; orcid.org/
0000-0002-7643-5911; Email: zbigniew.postawa@uj.edu.pl
Authors Dawid Macia̧ zek − Smoluchowski Institute of Physics, Jagiellonian University, 30-348 Krakow, Poland Micha Kanski − Smoluchowski Institute of Physics, Jagiellonian University, 30-348 Krakow, Poland Complete contact information is available at: https://pubs.acs.org/10.1021/acs.analchem.0c01219
Notes The authors declare no competing financial interest. ■ ACKNOWLEDGMENTS This work was supported by the Polish National Science Center Program No. 2019/33/B/ST4/01778. MD simulations were performed at the PLGrid Infrastructure. ■ REFERENCES (1) Fletcher, J. S. Biointerphases
2015,
10,
018902.
(2) McGettrick, J. D.; Speller, E.; Li, Z.; Tsoi, W. C.; Durrant, J. R.; Watson, T. Org.
Electron.
2017,
49,85−93.
(3) Hourani, W.; Gorbenko, V.; Barnes, J. P.; Guedj, C.; Cipro, R.; Moeyaert, J.; David, S.; Bassani, F.; Baron, T.; Martinez, E. J.
Electron
Spectrosc.
Relat.
Phenom.
2016,
213,1−10.
7352 https://dx.doi.org/10.1021/acs.analchem.0c01219
Anal. Chem. 2020, 92, 7349−7353 Analytical Chemistry pubs.acs.org/ac
(4) Wucher, A.; Cheng, J.; Winograd, N. J.
Phys.
Chem.
C
2008,
112,
16550−16555.
(5) Sjovall, P.; Rading, D.; Ray, S.; Yang, L.; Shard, A. G. J.
Phys.
Chem.
B
2010,
114,
769−774.
(6) Maciazek, D.; Paruch, R. J.; Postawa, Z.; Garrison, B. J. J.
Phys.
Chem.
C
2016,
120,
25473−25480.
(7) Cheng, J.; Wucher, A.; Winograd, N. J.
Phys.
Chem.
B
2006,
110,
8329−8336.
(8) Shard, A. G.; Green, F. M.; Brewer, P. J.; Seah, M. P.; Gilmore, I. S. J.
Phys.
Chem.
B
2008,
112,
2596−2605.
(9) Ninomiya, S.; Ichiki, K.; Yamada, H.; Nakata, Y.; Seki, T.; Aoki, T.; Matsuo, J. Rapid
Commun.
Mass
Spectrom.
2009,
23,
1601−1606.
(10) Mahoney, C. M. Cluster Secondary Ion Mass Spectrometry: Principles and Applications; Wiley: Hoboken, NJ, 2013. (11) Winograd, N. Annu.
Rev.
Anal.
Chem.
2018,
11,29−48.
and
references
therein.
(12) Chan, W. L.; Chason, E. J.
Appl.
Phys.
2007,
101,
121301.
and
references
therein.
(13) Norris, S. A.; Aziz, M. J. Appl.
Phys.
Rev.
2019,
6,
011311.
and
references
therein.
(14) Bradley, R. M.; Harper, J. M. E. J.
Vac.
Sci.
Technol.,
A
1988,
6,
2390−2395.
(15) Makeev, M. A.; Cuerno, R.; Barabasi, A. L. Nucl.
Instrum.
Methods
Phys.
Res.,
Sect.
B
2002,
197,
185−227.
(16) Aste, T.; Valbusa, U. New
J.
Phys.
2005,
7,
122.
(17) Munoz-Garcia, J.; Cuerno, R.; Castro, M. Phys.
Rev.
B:
Condens.
Matter
Mater.
Phys.
2008,
78,
205408.
(18) Chini, T. K.; Datta, D. P.; Bhattacharyya, S. R. J.
Phys.:
Condens.
Matter
2009,
21,
224004.
(19) Norris, S. A.; Brenner, M. P.; Aziz, M. J. J.
Phys.:
Condens.
Matter
2009,
21,
224017.
(20) Norris,S.A.; Samela,J.; Bukonte, L.;Backman,M.; Djurabekova, F.; Nordlund, K.; Madi, C. S.; Brenner, M. P.; Aziz, M. J. Nat.
Commun.
2011,
2,
276.
(21) Norris, S. A.; Samela, J.; Vestberg, M.; Nordlund, K.; Aziz, M. J. Nucl.
Instrum.
Methods
Phys.
Res.,
Sect.
B
2014,
318,
245−252.
(22) Toyoda,
N.;
Mashita,
T.;
Yamada,
I.
Nucl.
Instrum.
Methods
Phys.
Res.,
Sect.
B
2005,
232,
212−216.
(23) Tilakaratne, B. P.; Chen, Q. Y.; Chu, W. K. Materials
2017,
10,
1056.
(24) Lozano, O.; Chen, Q. Y.; Tilakaratne, B. P.; Seo, H. W.; Wang, X. M.; Wadekar, P. V.; Chinta, P. V.; Tu, L. W.; Ho, N. J.; Wijesundera, D.; Chu, W. K. AIP
Adv.
2013,
3,
062107.
(25) Niehuis, E.; Mollers, R.; Rading, D.; Cramer, H.-G.; Kersting, R. Surf.
Interface
Anal.
2013,
45,
158−162.
(26) Garrison,
B.
J.;
Postawa,
Z.
Mass
Spectrom.
Rev.
2008,
27,
289−
315.
and
references
therein.
(27) Aziz, R. A.; Slaman, M. J. Mol.
Phys.
1986,
58,
679−697.
(28) Ziegler, J. F.; Biersack, J. P.; Littmark, U. The Stopping and Range of Ions in Matter; Pergamon: New York, 1985. (29) Tersoff, J. Phys.
Rev.
B:
Condens.
Matter
Mater.
Phys.
1989,
39,
5566−5568.
(30) Daw, M. S.; Baskes, M. I. Phys.
Rev.
B:
Condens.
Matter
Mater.
Phys.
1984,
29,
6443−6453.
(31) Foiles, S. M.; Baskes, M. I.; Daw, M. S. Phys.
Rev.
B:
Condens.
Matter
Mater.
Phys.
1986,
33,
7983−7991.
(32) Plimpton, S. J.
Comput.
Phys.
1995,
117,1−19.
(33) Russo, M. F.; Postawa, Z.; Garrison, B. J. J.
Phys.
Chem.
C
2009,
113,
3270−3276.
(34) Paruch, R.; Rzeznik, L.; Russo, M. F.; Garrison, B. J.; Postawa, Z. J.
Phys.
Chem.
C
2010,
114,
5532−5539.
(35) Postawa, Z.; Czerwinski, B.; Szewczyk, M.; Smiley, E. J.; Winograd, N.; Garrison, B. J. J.
Phys.
Chem.
B
2004,
108,
7831−7838.
(36) Paruch, R. J.; Postawa, Z.; Wucher, A.; Garrison, B. J. J.
Phys.
Chem.
C
2012,
116,
1042−1051.
(37) Carter, G.; Vishnyakov, V. Phys.
Rev.
B:
Condens.
Matter
Mater.
Phys.
1996,
54,
17647−17653.
7353 https://dx.doi.org/10.1021/acs.analchem.0c01219
Anal. Chem. 2020, 92, 7349−7353