Multiscale simulation of the focused electron beam induced deposition process

Focused electron beam induced deposition (FEBID) is a powerful technique for 3D-printing of complex nanodevices. However, for resolutions below 10 nm, it struggles to control size, morphology and composition of the structures, due to a lack of molecular-level understanding of the underlying irradiation-driven chemistry (IDC). Computational modeling is a tool to comprehend and further optimize FEBID-related technologies. Here we utilize a novel multiscale methodology which couples Monte Carlo simulations for radiation transport with irradiation-driven molecular dynamics for simulating IDC with atomistic resolution. Through an in depth analysis of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {W(CO)}_6$$\end{document}W(CO)6 deposition on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {SiO}_2$$\end{document}SiO2 and its subsequent irradiation with electrons, we provide a comprehensive description of the FEBID process and its intrinsic operation. Our analysis reveals that simulations deliver unprecedented results in modeling the FEBID process, demonstrating an excellent agreement with available experimental data of the simulated nanomaterial composition, microstructure and growth rate as a function of the primary beam parameters. The generality of the methodology provides a powerful tool to study versatile problems where IDC and multiscale phenomena play an essential role.


S1. MOLECULAR ADSORPTION ENERGY AND SIMULATION OF SURFACE DIFFUSION
The diffusion coefficient of molecules on a surface strongly depends on temperature T and the activation energy for diffusion E diff [1]. The latter is related to the adsorption energy E ads for the molecule on the substrate, leading to the following relation: where D 0 is the diffusion coefficient at very large temperatures and k is the Boltzmann's constant. The adsorption energy E ads depends on the strength of the adhesion forces between the molecule and the substrate. The molecule-substrate forces (arising in our simulations from van der Waals interactions between pairs of atoms) are modelled in molecular dynamics (MD) by means of the Lennard-Jones potential, which needs two parameters (interatomic distance and depth of the potential well). For W, C, Si, O and H atoms, different sets of van der Waals parameters can be found in the literature [2][3][4][5]. In the present study, these parameters have been combined in different ways, producing several sets which allow screening a wide region of potential adsorption energies, until reaching a diffusion coefficient in reasonable agreement with the only available (to the best of our knowledge) empirical information (estimated from the adjustment of model calculations to experimental data on nanostructure growth) [6]. Note that the choice could have been also performed relying on ab initio calculations of the adsorption energy, although we preferred to rely here on empirical information.
The diffusion coefficients and adsorption energies of W(CO) 6 molecules on hydroxylated SiO 2 have been obtained from MD simulations of 2 ns duration at T = 300 K for each set of parameters, as explained in [2]. The results are depicted in Fig. S1 by symbols; the solid line represents the best fit by means of Eq. (S1). The experimentally determined value [6] is shown by a dashed line (the adsorption energy was not empirically obtained). An adsorption energy around 2.5 eV yields a surface diffusion coefficient of 7.71 µm 2 /s, fairly close to the experimentally determined value of 6.64 µm 2 /s [6]. Parameters for C, H and O atoms come from [2], those for Si from [3], while for W they come from [5].

S2. MONTE CARLO CODE SEED FOR ELECTRON TRANSPORT AND INTER-ACTION CROSS SECTIONS FOR SILICA
The details of the SEED (Secondary Electron Energy Deposition) Monte Carlo (MC) code for energetic electron transport in solids are explained in detail in [7]. It is based on the calculation of (i) the differential inelastic scattering cross sections obtained by using the dielectric formalism [8,9], (ii) the electron-phonon quasi-elastic scattering cross-sections computed by the use of the Fröhlich theory [10] and (iii) the differential elastic scattering cross sections performed by the relativistic partial wave expansion method (RPWEM) [11] including the Ganachaud and Mokrani empirical correction for low electron energies (≤ 20-30 eV) [12]. It should be noted that coherent scattering effects due to the long De Broglie wavelength for very low energy electrons are not yet considered in any MC code. Nonetheless, the good comparison of the results obtained by means of SEED with experimental and reference data (see below or [7] for further examples) demonstrates that its accuracy is good enough for the purposes of the present investigation. From experiment [6] FIG. S1. Determination of molecular adsorption energy and diffusion coefficient. W(CO) 6 surface diffusion coefficient on hydroxylated SiO 2 as a function of the adsorption energy for different sets of van der Waals parameters. An empirically estimated value (the adsorption energy was not determined) is shown by a dashed line [6].
The dielectric formalism, very successful for the calculation of electronic interaction cross sections in condensed matter, requires the knowledge of the energy-loss function (ELF) of the target material, Im[−1/ (k, ω)] (where (k, ω) is the complex dielectric function), accounting for its electronic excitation spectrum for excitations of given energy and momentum, ω and k, respectively. An estimate for the mean binding energy of the outer-shell electrons of the target is needed in order to disentangle the processes of ionisation and electronic excitation [8]. The optical ELF ( k = 0) of SiO 2 (mass density 2.19 g/cm 3 ) has been taken from the compilation of [13] and extended to finite momentum transfers by means of the Mermin Energy Loss Function-Generalised Oscillator Strengths (MELF-GOS) method [14], tested for many condensed-phase materials. A mean binding energy of 12.2 eV has been estimated from data taken from [15]. The effect of surface hydroxylation in the electronic properties of SiO 2 is not deemed to affect much the electron propagation through the substrate, which mostly occurs within the bulk of the material. The calculated total inelastic mean free path is shown by a solid line in Fig. S2(a) and compared to the available experimental data [16][17][18][19], finding a rather good agreement in a wide energy range.
The elastic mean free path for electrons in SiO 2 calculated by meas of the RPWEM is shown in Fig. S2(b) (dotted line). This method is known to provide unreasonable small values (shorter than interatomic distances) for T < 10 eV. To amend this tendency, the Ganachaud and Mokrani empirical correction is used [12], which requires the setting of the Ganachaud-Mokrani parameter α. Figure S2(b) shows the effect of different choices of α in the resulting elastic mean free path. Furthermore, the Fröhlich theory [10] for electronphonon scattering requires the knowledge of the characteristic phonon energy, W ph , whose typical values are in the range 0.01 − 0.1 eV.
All elastic, inelastic and electron-phonon interactions have an impact in the trajectories of (especially low energy) electrons, affecting their ability to escape from the substrate and hence determining the secondary electron (SE) yield. Since inelastic and elastic (high energy) cross sections are fully determined (the former from the dielectric theory and the latter from the RPWEM calculations), only the setting of the free parameters for the Ganachaud-Mokrani and Fröhlich theories remains. These parameters have been set in order to correctly simulate the experimentally known SE yield for SiO 2 [20]. Figure S3   Values of α = 5 and W ph = 0.15 eV (solid lines) provide a very reasonable agreement with the experimental data (symbols) [20]. The reliability of the MC simulations has been tested by calculating the energy spectra as well as energy loss spectra of the reflected electrons (i.e., the reflection electron energy loss spectra, REELS) and comparing them to reference and experimental data. Figure S4(a) compares the calculated (solid histogram) energy spectrum of ejected SE and backscattered electrons (BSE) for a 100 eV primary electron (PE) beam impinging in SiO 2 with reference simulations (dashed histogram) [21]. Figures S4(b) and (c) compare the simulated REELS spectra (red lines) with experimental data (black symbols) [22] for PE beams of 470 eV and 2000 eV, respectively. Simulations qualitatively agree with the energy spectrum from [21] and are in an excellent agreement with the experimental REELS spectra [22].

S3. W(CO) 6 MOLECULE FRAGMENTATION CROSS SECTIONS
For the W(CO) 6 molecule, experimental measurements exist for the relative cross sections for electron-impact dissociative ionisation (DI) [23] and for lower energy channels [24], although no information about the absolute cross sections is available.  6 , scaled from experiment Fe(CO) 5 , estimated from Ref. [27] Dissociation cross section (Å 2 ) Electron energy (eV) FIG. S5. Electron-impact fragmentation cross section of W(CO) 6 . Comparison of the scaled low energy dissociation cross section for W(CO) 6 (this work) with that of Fe(CO) 5 (estimated from experiments and calculations from [27]).
The experimental data on DI relative cross sections reveals that almost every ionising collision (≥ 97%) leads to W(CO) 6 fragmentation to some extent. Thus, it is possible to approximate the total DI cross section, σ DI (T ), by estimating the total ionisation cross section employing the dielectric formalism [8,9] (see Supplementary Information S2). The optical ELF of W(CO) 6 has been estimated from a parametric approach for organic materials [8,25] and a binding energy for the outer-shell electrons of 8.47 eV has been estimated from [23]. The calculated DI cross section is shown in Fig. 1(c) of the main text.
The absolute value of the lower energy dissociation channels cross section, σ low−T (T ), can be fixed by calculating the experimentally known decomposition cross section for W(CO) 6 molecules on SiO 2 by a PE beam of energy T 0 : where dN (T )/dT is the energy spectrum of PE and the generated SE and BSE crossing the surface near the beam area (further normalised to the number N PE of PE) and σ frag (T ) = σ DI (T ) + σ low−T (T ). The decomposition cross section has been experimentally determined to be 1.2Å 2 for 30 keV PE [26]. The energy spectrum is obtained from MC simulations (see Fig. 1(b) of the main text). The experimental relative low energy dissociation cross section σ low−T (T ) [24] has been scaled so the integral of Eq. (S2) gives a value of 1.2Å 2 . The resulting low energy cross section is shown by a solid line in Fig. S5, and is compared with that for the similar metal carbonyl Fe(CO) 5 . For the latter molecule, measurements of the relative cross sections also exist, as well as calculations of its absolute values [27], which allow estimating the absolute cross section (dashed line in Fig. S5). The scaled cross section for W(CO) 6 is comparable with the absolute cross section for Fe(CO) 5 , presenting the main characteristic peak (corresponding to the loss of a single CO ligand [24]), of a similar height, at T ∼ 0 eV. Other dissociation channels observed for W(CO) 6 at 3.3 and 4.7 eV (corresponding to the loss of two and three CO ligands, respectively [24]) seem to be much weaker in Fe(CO) 5 . The good comparison validates the scaling procedure employed.

S4. IRRADIATION DRIVEN MOLECULAR DYNAMICS SIMULATIONS OF THE FEBID PROCESS WITH MBN EXPLORER
The fundamentals of irradiation driven molecular dynamics (IDMD) are summarised in "Methods" and full details can be found in [2]. Here, the main aspects influencing simulations within the current investigation are explained.

A. Scaling of primary electron beam currents in simulations
The molecular fragmentation rate is influenced (see "Methods") by the PE flux, i.e., the number N PE of delivered PE per unit time t and unit area S at the irradiated spot: where I 0 = e dN PE /dt is the PE beam current, S 0 = πR 2 is the PE beam area (with R being the beam radius) and e the elementary charge. Typical FEBID experiments use a beam radius of several nanometres. Here, it was set to R = 5 nm. For simplicity, uniform PE fluxes J 0 were simulated by the MC code SEED (Supplementary Information S2) within the beam area S 0 . In all MBN Explorer simulations, the 20×20 nm 2 substrate is covered by 1-2 layers of W(CO) 6 molecules (density 5.4 molecules/nm 2 ), which guarantees the full coverage of the substrate while keeping a layer thin enough to not significantly affect PE trajectories. Currents of the order I 0 ∼ pA-nA for irradiation (dwell) times τ d ∼ µs or longer are commonly used in experiments. However, currently it would be computationally very costly to perform MD simulations for such long times. Instead, in this work a number of PE per unit area and per dwell time similar to experiments was sought. For this purpose, the simulated currents (and hence the fluxes) are scaled in the following manner [2]: where λ = τ exp d /τ sim d ; the super-indexes "sim" refer to simulation beam parameters while "exp" to experimental beam parameters. In such a way, simulated beam currents I sim 0 ∼ µA give the same amount of PE per simulated dwell time τ sim d ∼ ns (which is feasible for MD simulations) as experimental beam currents I exp 0 ∼ nA for experimental dwell times τ exp d ∼ µs. In Table S1, the PE beam parameters corresponding to each simulation (this work) and each experiment from [28] (to which we compare) are summarised.

B. Energy deposited to fragmenting bonds
Bond dissociation events are simulated by the deposition of an average energy per electron-molecule collision, E dep , in the atoms forming the bond, so that the atomic velocities increase fulfilling the requirements of momentum conservation [29]. For the sake of TABLE S1. Beam parameters corresponding to each simulation (present work) and to each experiment reported in [28]. Simulations  simplicity, in this work it is assumed that every fragmentation event leads to the cleavage of a single W-C bond, while the much stronger C-O bonds will stay intact [29]. Further collisions of the fragments with the environment may lead to subsequent cleavage of additional W-C bonds [29]. The choice of E dep influences the kinetics of the electron-driven chemical reactions, since low values favour metal-ligand recombination after dissociation, while larger values effectively put the metal atom and the ligand apart. A value of E dep = 325 kcal/mol (∼ 14 eV) has been chosen in simulations, since it is consistent with average values obtained from mass spectrometry experiments [24,30] and dedicated W(CO) 6 fragmentation simulations [29].
Unfortunately, the kinetics of W(CO) 6 molecules fragmentation has not been measured on SiO 2 substrates, although it is known for gold substrates [31]. In [31], 1-2 layers of W(CO) 6 molecules were deposited on Au, as in present simulations in SiO 2 , and the effective decomposition cross section σ decomp (T 0 ) (see Supplementary Information S3) was measured for a T 0 = 500 eV broad PE beam. On the one hand, this dense packing of molecules on the substrate prevents surface diffusion to some extent. On the other hand, the broad macroscopic beam (of ∼ 1 cm 2 area [31]) characterised by a decomposition cross section σ decomp (T 0 ) produces a uniform fragmentation probability P = J 0 σ decomp (T 0 ) over a nanometric area. As a consequence, it might be possible to, at least qualitatively, compare simulations on SiO 2 to the reported experiments on a Au substrate [31], despite their potential differences in terms of surface diffusion and electron emission, provided that the empirical decomposition cross section is employed for these specific simulations.
The measured fraction of CO ligands remaining on the Au substrate as a function of the PE dose is shown by symbols in Fig. S6. IDMD simulations have been performed for an spatially uniform beam, employing the experimentally determined decomposition cross section σ decomp (T 0 ) for W(CO) 6 molecules on gold [31]. Note that, in this particular case, there is no need to rely on MC simulation results, since the fragmentation rate P = J 0 σ decomp (T 0 ) is constant all over the simulated surface. The value E dep has been scanned in the range 70-400 kcal/mol, consistent with gas-phase findings [24,30]. As can be seen from the simulated curves, a value of 300-325 kcal/mol (∼13-14 eV) reproduces fairly well the main features of the experimental data, this value reasonably agreeing with the average deposited energies estimated from gas-phase mass spectrometry [24,30]. This comparison further supports our choice for E dep .
C. Models for the chemistry of reactive sites within the growing nanostructures Once molecules are fragmented and chemically reactive sites start to appear in the system, atoms with dangling bonds will start to form new bonds as they encounter other atoms with unsaturated valencies, following the rules of the reactive force field [32]. In a previous work [2], two models for the formation of new chemical bonds were introduced: model A only accounts for the chemistry between newly added precursors and the growing deposit, without restructuring of the dangling bonds within the latter (i.e., the search for reactive neighbours is done only among the atoms located beyond the molecular structure to which a chosen atom belongs). Model B allows for the formation of new bonds within the deposit itself (i.e., the search is performed over all atoms with dangling bonds in the simulation box, including the molecular structure to which a chosen atom belongs). As can be seen in the main text,