E Prakash

Bitzek

J Guénolé

A Prakash

E Bitzek

Atomistic simulations of focused ion beam machining of strained silicon

Atomistic simulations of focused ion beam machining of strained silicon

Introduction

Over the last decades, focused ion beam (FIB) techniques have become widely used in the material science community [START_REF] Tseng | Recent developments in micromilling using focused ion beam technology[END_REF]. Typical applications include imaging [START_REF] Ishitani | Origins of material contrast in scanning ion microscope images[END_REF], sample preparation for electron microscopy [START_REF] Kang | Ultra-thin TEM sample preparation with advanced backside FIB milling method[END_REF] or atom probe tomography [START_REF] Larson | Field-ion specimen preparation using focused ion-beam milling[END_REF], the fabrication of micro or nano scale specimens for mechanical testing [START_REF] Kiener | A further step towards an understanding of size-dependent crystal plasticity: In situ tension experiments of miniaturized single-crystal copper samples[END_REF], and its combined application with digital image correlation (DIC) to determine intrinsic stresses [START_REF] Korsunsky | Focused ion beam ring drilling for residual stress evaluation[END_REF][START_REF] Ahmed | Assessment of stress relaxation experiments on diamond coatings analyzed by digital image correlation and micro-raman spectroscopy[END_REF]. FIB milling has, however, been known to induce irradiation damage in the machined material; the influence of irradiation-induced defects on the material properties is a major drawback of the method [START_REF] Volkert | Focused ion beam microscopy and micromachining[END_REF]. The damage induced by FIB irradiation has been extensively studied, both by experiments [START_REF] Lugstein | Method to characterize the three-dimensional distribution of focused ion beam induced damage in silicon after 50 keV Ga + irradiation[END_REF][START_REF] Wang | Surface damage induced by focused-ion-beam milling in a Si/Si pn junction cross-sectional specimen[END_REF][START_REF] Kiener | FIB damage of Cu and possible consequences for miniaturized mechanical tests[END_REF][START_REF] Cullen | Assessment of surface damage and sidewall implantation in AlGaN-based high electron mobility transistor devices caused during focused-ion-beam milling[END_REF] and simulations [START_REF] Pastewka | Surface amorphization, sputter rate, and intrinsic stresses of silicon during low energy Ga + focused-ion beam milling[END_REF][START_REF] Tong | Molecular dynamic simulation of low-energy FIB irradiation induced damage in diamond[END_REF][START_REF] Xiao | Annealing recovery of nanoscale silicon surface damage caused by Ga focused ion beam[END_REF], as well as using combined approaches [START_REF] Lane | Validating molecular dynamics with direct imaging of radiation damage debris[END_REF]. In these studies, the authors usually focus on the influence of the ion beam energy, incidence angle, ion type etc. on the sputter process and irradiation in different target materials. In particular for materials like silicon, which upon irradiation exhibit structural amorphization, this damage induces stress and strain, as shown by atomistic simulation [START_REF] Pastewka | Surface amorphization, sputter rate, and intrinsic stresses of silicon during low energy Ga + focused-ion beam milling[END_REF] and experiments [START_REF] Mayr | Ion-irradiation-induced stresses and swelling in amorphous ge thin films[END_REF]. Materials undergoing irradiation can furthermore exhibit intrinsic stresses up to 1.5 GPa [START_REF] Treml | High resolution determination of local residual stress gradients in single-and multilayer thin film systems[END_REF], as those used in strain engineering [START_REF] Cheng | Electron mobility enhancement in strained-si n-mosfets fabricated on sige-on-insulator (sgoi) substrates[END_REF][START_REF] Jacobsen | Strained silicon as a new electro-optic material[END_REF] or in the case of multi-layer structures [START_REF] Treml | High resolution determination of local residual stress gradients in single-and multilayer thin film systems[END_REF][START_REF] Lu | In situ control of atomic-scale si layer with huge strain in the nanoheterostructure nisi/si/nisi through point contact reaction[END_REF]. Only a few atomistic studies have looked into the relation between the intrinsic stress state and low-energy irradiation processes in silicon [START_REF] Pastewka | Surface amorphization, sputter rate, and intrinsic stresses of silicon during low energy Ga + focused-ion beam milling[END_REF][START_REF] Beardmore | Direct simulation of ion-beam-induced stressing and amorphization of silicon[END_REF][START_REF] Korsunsky | Quantifying eigenstrain distributions induced by focused ion beam damage in silicon[END_REF]. A thorough understanding of the complex interplay between the stress or strain state and irradiation induced damage and sputtering process is, however, still lacking. This fundamental question is of particular concern for so-called FIB-DIC techniques which combine FIB milling with DIC in order to measure the intrinsic stress from the release of surface strains. In this context, residual stresses up to 6 GPa have been determined previoulsy in the literature [START_REF] Korsunsky | Focused ion beam ring drilling for residual stress evaluation[END_REF].

Here, we present a detailed atomistic study of FIB milling on strained silicon using gallium ions. The manuscript is organized as follows: We first present the particular methods we use for the classical atomistic modeling of ionic irradiation. The results for samples with and without strain are then described and discussed in the light of the calculation of the threshold displacement energy (TDE) and the surface binding energy (SBE) of strained silicon.

Simulation methods

Sample geometry

Most of the previous work on FIB induced damage in silicon used slab-like samples (corresponding to an infinitely extended surface) to study the influence of the beam incidence angle and the ion dose [START_REF] Pastewka | Surface amorphization, sputter rate, and intrinsic stresses of silicon during low energy Ga + focused-ion beam milling[END_REF][START_REF] Xiao | Annealing recovery of nanoscale silicon surface damage caused by Ga focused ion beam[END_REF][START_REF] Russo | Trench formation and lateral damage induced by gallium milling of silicon[END_REF][START_REF] Krasheninnikov | Ion and electron irradiation-induced effects in nanostructured materials[END_REF][START_REF] Das | Mechanisms of material removal and mass transport in focused ion beam nanopore formation[END_REF] on the sputter process and the resulting damage. Only few studies considered other sample geometries, like lamella [START_REF] Holmström | Atomic-scale effects behind structural instabilities in si lamellae during ion beam thinning[END_REF].

In the present study, we choose two sample geometries. The first is a slab-like geometry (Fig. 1a) similar to the ones used in the aforementioned studies. It consists of one (0 0 1) free surface and an ion beam with a normal incidence angle θ [001] = 1 • , with θ being the angle between the ion beam and the normal to the targeted (0 0 1) surface. This small deviation of 1 • is chosen to reduce ion channeling. The second geometry consists of two orthogonal {0 0 1} free surfaces and an ion beam centered at the common corner of these two surfaces (Fig. 1b). The (0 0 1) surface is hence irradiated under normal incidence with θ [001] = 1 • , whereas the (1 0 0) surface is irradiated under an extreme grazing incidence with θ [100] = 89 • . Note that the {0 0 1} family of surfaces chosen in this study is typical for experiments on silicon wafers [START_REF] Wang | Surface damage induced by focused-ion-beam milling in a Si/Si pn junction cross-sectional specimen[END_REF][START_REF] Ishii | Stress evolution in si during low-energy ion bombardment[END_REF][START_REF] Shariq | Investigations of fieldevaporated end forms in voltage-and laser-pulsed atom probe tomography[END_REF][START_REF] Giri | Studies on the surface swelling of ion-irradiated silicon: Role of defects[END_REF]. The sample dimensions are chosen to have a balance between the large volume required for correct stress-strain relaxation and reasonable computational times: 60×20×74 a 0 (33×11×43 nm 3 ) in the x, y and z directions, respectively, with a 0 = 0.5431 nm being the lattice parameter at 0 K of the chosen Si potential [START_REF] Pizzagalli | A new parametrization of the stillinger-weber potential for an improved description of defects and plasticity of silicon[END_REF]. This results in initial samples with approximately 0.7 million atoms. Periodic boundary conditions (PBC) are only applied in the y-[0 1 0] direction (Fig. 1a,b). The boundary conditions in the z-[0 0 1] and x-[1 0 0] are either free or 2D-tethered, as shown in Fig. 1(a,b). Tethered boundary conditions (TBC) do not prohibit but rather restrain the displacement of atoms. The displacement restriction is controlled by a spring which links a particle to its original position (here we used a stiffness of k = 10 9 N • m -1 ). Note that in the case of 2D-TBC, the tethered atoms can still move freely in the plane of the corresponding surface to allow the dilatation/contraction of the whole sample.

Interatomic potentials

For the Ga-Si interaction, the Ziegler-Biersack-Littmark (ZBL) potential is used [START_REF] Ziegler | The stopping and range of ions in matter[END_REF]. This universal repulsive potential is designed to model high energy atomic collisions, when the influence of the chemistry becomes negligible. After the collision cascade, the chemistry of the Ga + ion should, however, be considered and suitable inter-atomic models for Si-Ga and Ga-Ga interactions used. Since currently no well-tested binary potentials exist, the Ga + ions are usually kept inside the sample with an approximated interaction model ignoring atomic bonding (e.g. the ZBL potential), or they are removed from the system after the end of the collision cascade. The latter approach was chosen in this work, since it has been previously shown to not significantly influence the defect formation in comparison to retaining the Ga + ions with ZBL interactions [START_REF] Pastewka | Surface amorphization, sputter rate, and intrinsic stresses of silicon during low energy Ga + focused-ion beam milling[END_REF].

We use a recent parametrization of the Stillinger-Weber potential (SW m ) to model the Si-Si interaction [START_REF] Pizzagalli | A new parametrization of the stillinger-weber potential for an improved description of defects and plasticity of silicon[END_REF]. This parametrization improves the description of defect-related properties in silicon, in particular the threshold displacement energies which is a critical quantity for irradiation cascades simulations. The properties of the amorphous phase, often observed in irradiation simulations, are also improved [START_REF] Pizzagalli | A new parametrization of the stillinger-weber potential for an improved description of defects and plasticity of silicon[END_REF]. During a collision cascade, some Si atoms, particularly the primary knock-on atoms (PKA), will have a very high energy. To represent such Si-Si interaction at high energies, the SW m potential is smoothly merged to the ZBL potential for short inter-atomic distances [START_REF] Holmström | Atomic-scale effects behind structural instabilities in si lamellae during ion beam thinning[END_REF]. The two potentials are merged over a distance between 1.7 -2.0 Å; the details of the joining procedure are described in [START_REF] Nordlund | Molecular dynamics simulations of threshold displacement energies in fe[END_REF].

Molecular dynamics simulations

The MD simulations presented in this work are performed with the LAMMPS software (version 7 May 2015) [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF], using the Verlet algorithm with a default timestep of 1.0 fs. Simulations in the NVT and NPT ensemble use Nosé-Hoover thermostat and barostat [START_REF] Hoover | Canonical dynamics : Equilibrium phase-space distributions[END_REF][START_REF] Shinoda | Rapid estimation of elastic constants by molecular dynamics simulation under constant stress[END_REF]. Samples are cut from the bulk material and statically relaxed using the conjugate gradient algorithm. The structures are then homogeneously scaled to fit the lattice parameter at 300 K, and thermalized at 300 K by using successively the NVE and NVT thermodynamic ensembles for a total time of 20 ps. Successive cascade simulations were finally run from this initial structure.

During a collision cascade, the heat spike regime requires sufficient energy dissipation without perturbing the trajectories of the high velocity particles. The kinetic energy provided by the fired ions is dissipated using a thermostatted region of 2 nm thickness, parallel to the ion beam and placed as far as possible from the initial position of the fired ion (Fig. 1a,b). Consequently, the position of the thermostatted layer is different for each collision cascade. In this layer, the temperature is controlled by a Berendsen thermostat [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF]. Even if the thermostatted region acts as a barrier for the ion propagation, the influence on the effective damage formation is negligible and a similar scheme has also been successfully used previously [START_REF] Holmström | Atomic-scale effects behind structural instabilities in si lamellae during ion beam thinning[END_REF].

The position of Ga + ions is randomly determined at a distance of 2 nm from the (0 0 1) surface by following a normal distribution with a standard deviation of 2 nm typical for atomistic simulations [START_REF] Russo | Trench formation and lateral damage induced by gallium milling of silicon[END_REF][START_REF] Holmström | Atomic-scale effects behind structural instabilities in si lamellae during ion beam thinning[END_REF][START_REF] Das | Mechanisms of material removal and mass transport in focused ion beam nanopore formation[END_REF]. Figure 1c shows an example of such a distribution with 500 ions. A velocity corresponding to a kinetic energy of 5 kV in the direction of the beam is attributed to the ion. An MD simulation -the collision cascade -is then performed for a maximum time of 15 ps with a variable timestep to ensure that no atoms will be displaced by a distance larger than 0.01 nm during one MD step. Because of the high energy involved, sputtering can occur. An atom is considered as sputtered and then removed from the simulation if it leaves the simulation box, i.e., the distance between the atom and the closest surface of the original sample being larger than two times the cutoff of the potential. For all ions, the collision cascade is finished by the end of these 15 ps, with a final global temperature below 400 K. A final MD simulation is performed for 15 ps to slowly cool down the system from the current temperature to 300 K, and obtain the structure ready to run the next collision cascade.

As reference, silicon bulk samples with amorphous structure are obtained by the classical melt-quench procedure. A bulk crystalline Si sample is melted and quenched to 300 K at a quench rate of 10 13 K s -1 , typical for atomistic simulations [START_REF] Pizzagalli | A new parametrization of the stillinger-weber potential for an improved description of defects and plasticity of silicon[END_REF][START_REF] Guénolé | Plasticity in crystalline-amorphous core-shell si nanowires controlled by native interface defects[END_REF][START_REF] France-Lanord | Atomistic amorphous/crystalline interface modelling for superlattices and core/shell nanowires[END_REF]. This amorphous structure is maintained at 300 K for 60 ps, with a barostat ensuring a globally stress free sample within the NPT thermodynamic ensemble.

Calculation of threshold displacement energy and surface binding energy

The TDE is an intrinsic material property. It corresponds to the energy required to move an atom from its original stable position to a meta-stable position, creating a point defect [START_REF] Was | Fundamentals of Radiation Materials Science[END_REF]. The TDE is a critical input parameter for simulations of collision cascades in, e.g., the binary collision approximation (BCA) framework [START_REF] Pizzagalli | A new parametrization of the stillinger-weber potential for an improved description of defects and plasticity of silicon[END_REF][START_REF] Holmström | Threshold defect production in silicon determined by density functional theory molecular dynamics simulations[END_REF]. In crystals, this threshold is highly anisotropic. In addtion to the crystallographic direction it also depends on the thermal noise which slightly displaces the atoms from their ground state. In this context, the determination of an average TDE can be performed by using a well-constructed statistical sample over N = 1000 randomly chosen crystallographic directions [START_REF] Holmström | Threshold defect production in silicon determined by density functional theory molecular dynamics simulations[END_REF]:

TDE =   2π θ=0 π φ=0 E TDE (θ, φ)   × N -1 ,
where θ, φ are angles defining the spatial direction of E TDE and N the total number of TDE calculations.

In our work, we compute the TDE by using a scheme similar to the one suggested by Nordlund et al. [START_REF] Nordlund | Molecular dynamics simulations of threshold displacement energies in fe[END_REF]. A perfect bulk crystal of 8000 atoms was thermalized at 36 K and zero pressure. A velocity v θ,φ according to an energy of 2 eV is prescribed to a randomly chosen atom in a random direction [θ, φ], and a MD simulation is performed in the NVE ensemble for 6 ps. The velocity v θ,φ is then increased by steps of 2 eV until the formation of a defect (in our case Frenkel pairs) is observed at an energy E TDE (θ, φ). Defects are identified using the Voronoi cell method [START_REF] Rycroft | Voro++: A three-dimensional voronoi cell library in c++[END_REF][START_REF] Stukowski | Structure identification methods for atomistic simulations of crystalline materials[END_REF].

The SBE characterizes the bonding strength of surface atoms. Commonly likened to the sublimation energy [START_REF] Kudriavtsev | Calculation of the surface binding energy for ion sputtered particles[END_REF][START_REF] Behrisch | Sputtering by Particle Bombardment: Experiments and Computer Calculations from Threshold to MeV Energies[END_REF], it is possible to determine this quantity from atomistic simulations [START_REF] Behrisch | Sputtering by Particle Bombardment: Experiments and Computer Calculations from Threshold to MeV Energies[END_REF]. The SBE for the SW m potential used in our work is computed in a scheme similar to the one of Yang et al. [START_REF] Yang | Atomic scale calculations of tungsten surface binding energy and beryllium-induced tungsten sputtering[END_REF]. Considering the free surface of a perfect crystal at 0 K, a velocity parallel to the outward-pointing normal of this surface is prescribed to an atom, and a MD simulation is performed in the NVE ensemble for 2 ps. Being large enough to tear the atom from the surface, this velocity is decreased step by step till the atom remains attached to the surface. The energy corresponding to the step right before this final stage is then defined as the SBE. Note that the system being at 0 K without any thermal noise, all the non-reconstructed (0 0 1) surface atoms are equivalent, and no statistical sampling is required.

Structural analysis

The atomistic structures are visualized using Ovito [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovitothe open visualization tool[END_REF]. The analysis of irradiation induced damage is mostly performed using the diamond structure identification method [START_REF] Maras | Global transition path search for dislocation formation in ge on si(001)[END_REF]. Briefly, this method analyzes the crystallographic environment of each atom to establish a classification. Atoms in a perfect diamond cubic structures -first and second nearest neighbors located on diamond cubic lattice sites -can be distinguished from atoms in a perfect hexagonal diamond configuration and from atoms in other environments. Moreover, atoms on diamond lattice sites having first or second nearest neighbors on non-diamond lattice sites, can also be distinguished from defect atoms. A detailed description of the method is provided by Maras et al. [START_REF] Maras | Global transition path search for dislocation formation in ge on si(001)[END_REF].

Monte-Carlo simulations

Monte-Carlo simulations are performed within the binary collision approximation (BCA) framework by using the SRIM/TRIM package [START_REF] Ziegler | The stopping and range of ions in matter[END_REF][START_REF] Ziegler | SRIM the stopping and range of ions in matter[END_REF]. In this commonly used framework, each collision between the ion and an atom in the sample, or between two atoms in the sample, is considered as being independent. Probabilities based on intrinsic material properties indicate whether target atoms are displaced or sputtered: the TDE is used for atomic displacements and the SBE for sputtering. In this framework, the geometry of the targets is limited to an infinite thin film with a defined thickness. Furthermore, any anisotropy in material properties is neglected, i.e., the material is assumed to be amorphous. For the calculations, the default parameters of the SRIM/TRIM package [START_REF] Ziegler | Interactions of ions with matter[END_REF] for Si irradiated by Ga + (TDE E disp = 15 eV, lattice binding energy E latt = 2 eV, SBE E surf = 4.7 eV) as well as parameters determined in the present work using MD (TDE E disp = 38.3 eV, SBE E surf = 10.6 eV) are used.

Results

TRIM simulations

To compare TRIM and MD, collision cascade simulations using Si samples irradiated by Ga + ions at 5 kV are shown in Fig. 2(a) and Fig. 2 (b). The TRIM calculation uses the parameters determined in the present work, see sect. 2.4, and shows the path of displaced atoms. The MD simulation on the other hand displays the crystalline defects at the end of the irradiation processes. Note that the single surface geometry is used for comparison, since TRIM simulations cannot be used with more complex geometries. The corresponding damage distributions and implanted Ga + are presented on Fig. 2(c), by considering a dose of 500 Ga + at 5 kV on Si targets. The visible extension of the damaged zone, Fig. 2 (b), and the location of the recoil atoms, Fig. 2(a), as well as the location of the maximum number of implanted ions are comparable for the TRIM and MD simulation. The MD simulations, however, show a larger spread of the implanted ions to larger depths than the TRIM calculations, see Fig. 2(c). When the TRIM calculation is performed using the standard values, the Ga + implantation depth remains unchanged compared to the TRIM calculation with the MD-derived parameters. This is due to the fact that in the TRIM calculations the Ga + -Si interactions depend only on the velocity and mass of the particles. For the sake of clarity, these results are therefore not included in Fig. 2(c). Using the default parameters in the TRIM calculation, however, significantly increases the number of recoil atoms compared to the calculations with the lower, MD-determined TDE value, see Fig. 2(c).

Irradiation-induced damage

The single surface or slab geometry is commonly used in MD and allows reliable comparisons with our results. The dual surface sample is closer to the typical trench geometry and thus allows to study the influence of strain on FIB milling process during specimen preparation. Regarding the irradiation-induced damage, we observe similar results with both sample geometries. In the following, only the results obtained with the dual surface geometry will be presented.

The evolution of the irradiation-induced damage with increasing number of launched ions is shown in Fig 3 . It is worth mentioning that only stable damage configurations as obtained after thermal equilibration are considered here. During a typical collision cascade simulation performed in this work, the region where the collisions occur is in the so-called heat spike regime. The very high temperatures reached at the atomic scale during this regime lead to local sample melting, and introduce a significant amount of damage. However, the thermal energy is fully dissipated before the occurrence of the next collision cascade and the molten region quickly re-crystallizes before undergoing further collisions cascades. The remaining structural damage is then stable at room temperature, the temperature considered in this study.

Before irradiation, only atoms belonging to the reconstructed surfaces are classified as defect atoms, see Fig. 3(c) at 0 fired ions. After 5 and 10 ions, a damaged region close to the paths of the ions is formed. Note that the cavity formed at the surface is very localized 3(b). Subsequent collision cascades melt the region several times and homogenize the surface structure. After 50 ions, the damaged area has significantly enlarged and a clear amorphous zone can be observed. With further increase in the ion dosage, the damaged region is significantly enlarged; a clear amorphous zone surrounded by smaller disordered regions is observed. The rounding of the edge is characteristic for both the sputtering and the amorphization.

From the number of sputtered atoms, a sputter rate S ≈ 5.3 atoms per incident ion could be determined. For comparison, irradiation of the slab geometry with a normal incidence angle results in a sputter rate of S ≈ 1.1 atoms/ion.

To characterize the damaged zone, the radial distribution function (RDF) g(r) and the atomic volume were calculated by considering only atoms enclosed in a sphere of radius R = 50 Å located entirely in the damaged region. The RDF was calculated before and after irradiation and also for a fully amorphous structure obtained by the melt-quench method [START_REF] Pizzagalli | A new parametrization of the stillinger-weber potential for an improved description of defects and plasticity of silicon[END_REF], see Fig. 4(a). The discrete peaks obtained for the nonirradiated structure are characteristic for a perfectly crystalline structure. The continuous curves obtained after irradiation or from the melt-quench prepared sample are characteristic of an amorphous structure, exhibiting a short range order (SRO). The distribution of the atomic volume, as calculated by Voronoi tessellation, inside the same sphere used to compute the RDF is shown in Fig. 4(b). The crystalline structure (b) distribution of atomic volumes inside this sphere, both before (grey=crystalline) and after irradiation (blue). As reference, also the RDF and the distribution of atomic volume of an amorphous sample generated by the melt-quench route is included (red). The corresponding atomic densities are also provided in the plots.

is characterized by an almost perfectly mono-disperse distribution. The structure after irradiation and the melt-quenched structure show a wide spread of atomic volumes, with a slight shift of the peak to volumes larger than the atomic volume in the crystalline case. Again, this large spread is characteristic for amorphous structures. However, the distribution of atomic volume and the RDF show clear differences between the amorphous sample generated by the typical melt-quench approach and the amorphous zone created by the irradiation. Most obvious is the 1% higher density of the irradiated zone (48.7 at/nm 3 ) when compared to the melt-quenched sample (48.2 at/nm 3 ). As the crystalline sample (49.4 at/nm 3 ) has only an about 2% higher density than the melt-quenched sample, this difference is significant. To study the influence of quenching rate on density, we also quenched the sample with an ultra-high quench rate comparable to the heat dissipation rate in a heat spike. In this case, with a quench rate of 3 • 10 14 K/s, the density is about 48.6 at/nm 3 .

Effect of strain

To study the influence of strain on irradiation-induced damage and sputtering, several strained samples are irradiated and analyzed. Compressive and tensile strains ε yy of 2% and 4% are applied on non-irradiated samples along the y-[0 1 0] periodic direction by homogeneously scaling the simulation box. During the subsequent simulations of ion irradiation the box-size is kept constant. Note that these levels of strain are comparable to those measured experimentally [START_REF] Korsunsky | Focused ion beam ring drilling for residual stress evaluation[END_REF][START_REF] Treml | High resolution determination of local residual stress gradients in single-and multilayer thin film systems[END_REF][START_REF] Chu | Strain: A solution for higher carrier mobility in nanoscale mosfets[END_REF][START_REF] Lu | In situ control of atomic-scale si layer with huge strain in the nanoheterostructure nisi/si/nisi through point contact reaction[END_REF].

The number of sputtered atoms as function of impacting ions is shown for the different strained samples in Fig. 5. While the unstrained sample (ε = 0) displays the highest sputter rate, the amount of sputtered atoms decreases with increasing absolute value of the strain, with the compressed sample at ε = -4% showing the lowest sputter yield. The sputter yield S was defined as the slope of linear fits to the plots shown in Fig. 5. The values of S are provided in Tab. 1.

To characterize the extent of the irradiation damage, we count atoms which, according to the diamond structure identification method, are not classified as perfectly crystalline, see Tab. 1. Note that Defect atoms shown by the sample before irradiation (Tab. 1, first column) are surface atoms. The percentage of atoms which form the irradiation-induced damaged region (Tab. 1, Irradiation defect atoms) is the ratio of the sum of Defect atoms and Sputtered atoms, reduced by removing surface atoms, over all atoms. This ratio shows an opposite trend than the sputter rate S. Without applied strain, fewer of the remaining atoms are part of a defect structure than with applied strain. The percentage of irradiation-induced defect atoms, and consequently the extent of the damaged region, is highest for the sample at 4% tensile strain, see Tab. 1.

A more detailed mechanical characterization of the defect region is not trivial due to the superposed uniaxial strain and their different extent, surface area, etc. In particular, the elastic constants of the amorphous zones cannot be directly determined. However, by performing quasi-static tensile tests on the pre-strained irradiated and un-irradiated samples one can assess whether irradiation under strain changes the elastic material response. The results are shown in Fig. 6. The procedure is akin to the standard way of determining elastic constants from the second derivative of the strain-energy curve [START_REF] Finnis | Interatomic forces in condensed matter[END_REF], where the samples, however, have a well-defined volume and are defect-free.

Interestingly, in case the sample is irradiated (at a dose of 300 ions) without any applied strain, the relative change in energy with applied tensile strain is identical to an unirradiated sample, see Fig. 6, where the solid circles for ε = 0 overlap with the open circles of the unirradiated sample. This is, however, not the case when the sample has been irradiated while strained. For the sample under compressive pre-strain, the energy decreases when a tensile strain is applied; however, the rate of change is decreased when the sample has been irradiated. For the sample under tensile pre-strain, the energy increases with further tensile strain, but again, the rate of change is decreased when the sample has been irradiated.

Threshold displacement energy and surface binding energy

The TDE is commonly regarded as a key parameter to determine and/or model the damage produced by a collision cascade [START_REF] Gibson | Dynamics of radiation damage[END_REF][START_REF] Erginsoy | Dynamics of radiation damage in a body-centered cubic lattice[END_REF]. While the TDE of unstrained silicon has already been determined by atomistic simulations and by ab-initio methods [START_REF] Pizzagalli | A new parametrization of the stillinger-weber potential for an improved description of defects and plasticity of silicon[END_REF][START_REF] Holmström | Threshold defect production in silicon determined by density functional theory molecular dynamics simulations[END_REF], we are not aware of any such calculations for strained Si. A recent work on BCC Fe, however, already pointed out the importance of considering the strain environement when computing the TDE [START_REF] Beeler | Effect of strain and temperature on the threshold displacement energy in body-centered cubic iron[END_REF]. The average TDE for Si subjected to various hydrostatic and uniaxial [0 1 0] strains is shown in Tab. 2. The uniaxially strained samples are stress free in the two other 1 0 0 directions. The errors correspond to the standard error of the data set, and are calculated over the population N (Tab. 2).

It is important to note that the average TDE reflects neither the anisotropy due to the underlying crystal structure, nor the anisotropy with respect to the uniaxial straining direction.

Few calculations or measurements of the SBE exist for silicon [START_REF] Pizzagalli | A new parametrization of the stillinger-weber potential for an improved description of defects and plasticity of silicon[END_REF][START_REF] Holmström | Threshold defect production in silicon determined by density functional theory molecular dynamics simulations[END_REF][START_REF] Feil | Random and ordered defects on ion-bombarded si(100)-(2x1) surfaces[END_REF]. The SBE depends inherently on the crystallographic orientation and structure of the surface. Unlike for metals, covalently bonded materials like silicon undergo short range surface reconstruction, which can influence the sputter yield [START_REF] Smith | kev particle bombardment of semiconductors: A molecular-dynamics simulation[END_REF][START_REF] Rubio | Molecular dynamics simulation of amorphous silicon sputtering by ar+ ions[END_REF]. These have, however, not yet been directly taken into account in the determination of the SBE [START_REF] Behrisch | Sputtering by Particle Bombardment: Experiments and Computer Calculations from Threshold to MeV Energies[END_REF]. For the 0 0 1 surfaces considered here, the reconstruction follows the π-bonded chain model of Pandey [START_REF] Pandey | New π-bonded chain model for si(111)-(2 × 1) surface[END_REF], named p(2×1) for symmetry reasons, and is characterized by the formation of bonds between surface atoms in the 0 1 1 directions. Instead of being two-fold coordinated, the reconstructed surface atoms are three-fold coordinated and are therefore expected to be more difficult to sputter. The SBE for as-cut and reconstructed surfaces at several strains is shown in Tab. Table 3: Surface binding energy (SBE) of (0 0 1) silicon surface for different values of uniaxial strain along the [0 1 0] axis. The surface is either bulk-terminated (as-cut) or reconstructed. A reference value obtained from a recent empirical model [START_REF] Kudriavtsev | Calculation of the surface binding energy for ion sputtered particles[END_REF] is added. The discrepancy between the reference and the present value is discussed in the text.

Discussion

Comparison of BCA and MD simulations

A comparison of the extent of the damaged zone in the MD simulations with the displaced (or recoil) atoms in the TRIM calculations, Fig. 2(a)-(c), shows good agreement between the two methods, nearly independent of the used parameter values for the TRIM calculations. This is rather surprising, given the underlying vastly simplified assumptions in the BCA as well as the large differences in the SBE and TDE determined for the used potential, Tabs. 2 and 3, and the TRIM default parameters for Si [START_REF] Ziegler | Interactions of ions with matter[END_REF], see sect. 2.6. The wider spread of Ga + implantation observed in MD leads to slightly larger average penetration depth, which can be attributed to channelling effects neglected in the BCA. Nonetheless, concerning the penetration depth of Ga-ions in Si, our results show that the peak of the distribution obtained from MD can be reproduced in the BCA, even using very different values of the SBE and TDE (Fig. 2c). Regarding the calculation of the maximum in the distribution of the Ga + ions, this points to a relative insensitivity of the BCA approach on the exact values of the SBE and TDE.

The number of recoil atoms, however, clearly depends on the exact value of the TDE, see Fig. 2(c). Using the standard values, the number of generated recoil atoms in the TRIM calculation is significantly larger than with the MD-informed TDE and SBE, and similar to the number of defect atoms in the MD calculations. Compared to the MD simulations, the lower number of recoil atoms in the MD-informed TRIM calculations, however, seems more realistic, as one recoil atom leads statistically to many defect atoms. A more quantitative comparison, e.g, similar to the work of Bukonte et al. [START_REF] Bukonte | Comparison of molecular dynamics and binary collision approximation simulations for atom displacement analysis[END_REF], who used an optimization procedure to match the BCA parameters to their MD simulations with an environment-dependent interatomic potential (EDIP), is outside the scope of this paper.

Our results, however, demonstrate that for detailed Monte-Carlo studies of sputtering, crater formation, surface stability or topographic pattern formation, TDE and SBE values determined from atomistic simulations might be necessary, in particular for strained samples. Furthermore, the influence of surface reconstruction on the SBE should not be neglected, as it can lead to an increase of the SBE by more than a factor of two [START_REF] Kudriavtsev | Calculation of the surface binding energy for ion sputtered particles[END_REF][START_REF] Guénolé | Influence of intrinsic strain on irradiation induced damage: the role of threshold displacement and surface binding energies[END_REF].

Irradiation damage in unstrained samples

The sample studied here models the typical situation encountered in FIB-milling of trenches. It combines a surface irradiated under normal incidence with the irradiation of a surface under grazing angle. As ions with normal incidence lead to relatively few sputter events, the sputter rate is dominated by the impacts under grazing angle. By considering the combination of normal and grazing incidences, the sputter yield of S = 5.3 atoms/impact in our simulations is therefore comparable to S ≈ 6.5 atoms/impact determined by Pastewka et al. for a beam energy of 5 keV using MD simulations under grazing angle, however, with a different Si potential [START_REF] Pastewka | Surface amorphization, sputter rate, and intrinsic stresses of silicon during low energy Ga + focused-ion beam milling[END_REF].

The formation of an amorphous zone in Ga-irradiated silicon is expected, and was also found in other MD simulations [START_REF] Pastewka | Surface amorphization, sputter rate, and intrinsic stresses of silicon during low energy Ga + focused-ion beam milling[END_REF][START_REF] Xiao | Annealing recovery of nanoscale silicon surface damage caused by Ga focused ion beam[END_REF][START_REF] Russo | Trench formation and lateral damage induced by gallium milling of silicon[END_REF][START_REF] Satake | Molecular dynamics simulation of ga+ ion collision process[END_REF]. The comparison of the density, atomic volume distribution and RDF, however, reveals clear differences to amorphous silicon produced by the conventional melt-quench approach, see Fig. 4. As expected [START_REF] Laaziri | High resolution radial distribution function of pure amorphous silicon[END_REF], the amorphous structure generated by quenching a melt with a typical quench rate of 10 13 K/s is about 2% less dense than crystalline Si (Fig. 4). However, the irradiated structure is only about 1% less dense. Also, the RDF of the irradiation-induced amorphous structure is going faster to zero before the first peak than the RDF of the melt-quench structure. Conversely, the second peak is slightly wider. The irradiation-induced amorphous structure thus shows proportionally less close-packed nearest neighbors than the melt-quench structure, but closer second nearest neighbors. Ion irradiation thus clearly leads to a more dense amorphous zone than the overall density of typical melt-quenched samples. During an irradiation cascade, the energy dissipation -the quenching -is much faster than the quench rate usually used in the production of glass structures in MD simulations [START_REF] De La Rubia | Role of thermal spikes in energetic displacement cascades[END_REF]. This contributes to the increased density, as shown by the density of the ultra-fast quenched sample which differs only by 0.2% from the irradiated structure. Furthermore, during irradiation, the locally molten region is constrained by the surrounding crystal [START_REF] Trinkaus | Viscoelastic model for the plastic flow of amorphous solids under energetic ion bombardment[END_REF]. The melt, and subsequently the quenched-in amorphous region, cannot freely expand during the irradiation cascade, which might also lead to higher density compared to the meltquench process under zero pressure boundary conditions. This difference in the amorphous Si structures between irradiation-induced amorphization and melt-quenched samples has important consequences: typical amorphous Si samples cannot be used as a general model for irradiation-induced damage in Si.

Sputtering and irradiation damage in strained samples

Our results show a clear influence of uniaxial strain on the sputter yield (Fig. 5 and Tab. 1), the amount of induced damage (Tab. 1), and the elastic response (Fig. 6) of the irradiated sample. Considering very lowenergy irradiation (< 150 eV) of amorphous silicon, Koster and Urbassek [START_REF] Koster | Atomistic simulation of stress effects in a-si due to low-energy si impact[END_REF] already observed an increased damage production and a larger stress relaxation in a sample compressed to -3.3%, when compared to a stress-free sample. However, for crystalline silicon under tension and irradiated with more energetic particles (< 5 keV), Pastewka et al. [START_REF] Pastewka | Surface amorphization, sputter rate, and intrinsic stresses of silicon during low energy Ga + focused-ion beam milling[END_REF] reported no effect of strain on the sputter rate or on the amount of damage. In contrast to the present study, Pastewka et al. used only tensile strains up to +1.3% and focused on monitoring the thickness of the amorphous layer formed during irradiation under grazing angle. The use of significantly larger -but not unrealistic [START_REF] Chu | Strain: A solution for higher carrier mobility in nanoscale mosfets[END_REF][START_REF] Lu | In situ control of atomic-scale si layer with huge strain in the nanoheterostructure nisi/si/nisi through point contact reaction[END_REF] -tensile and compressive strains, as well as more detailed analysis methods in our work therefore seem to be key to observe the effects of strain on irradiation and sputtering, in MD simulations of Si.

In general, (uniaxial) applied strain can affect the irradiation cascade and subsequent irradiation damage in many different ways. First, strain changes the density and thereby directly affects the ion-irradiation [START_REF] Koster | Atomistic simulation of stress effects in a-si due to low-energy si impact[END_REF]. Uniaxial strain leads, furthermore, to a direction-dependence of the material properties in addition to the crystallographic anisotropy. Strain also directly affects the strength of bonds and the cohesive energy. Strained bonds should therefore be more susceptible to bond breaking by interactions with high-energy ions. Due the anharmonic nature of the Si bonds at large strains, this effect is expected to be different in tension than in compression. Furthermore, strain also influences the energy barriers experienced during the motion of knock-on atoms, as well as energetic costs to create a defect like an interstitial. For uniaxial strain, these effects will be anisotropic.

The TDE captures the aspects of bond breaking, the energy barriers for migration of knock-on atoms and the energy to form Frenkel pairs. For hydrostatic compression, the average TDE shows a clear increase compared to the unstrained sample, whereas for tensile hydrostatic strain states, the TDE is significantly reduced, see Tab. 2. The nearly linear relationship between TDE and hydrostatic strain points towards the energy to create a self-interstitial as the dominating effect for the strain-dependence of the TDE. A similar trend was recently observed in a study considering BCC iron, with a decrease of the TDE while applying hydrostatic tensile strain [START_REF] Beeler | Effect of strain and temperature on the threshold displacement energy in body-centered cubic iron[END_REF]. The trend is similar, but less pronounced for uniaxial strain. Here the averaged TDE does not include information about the direction dependence, but relaxation in the directions orthogonal to the uniaxial strain reduce its effect. That the effect of strain is stronger on the formation of the resulting damaged structure than on the bond-breaking can also be seen by the negligible effect it has on the SBE, see Tab. 3.

The strain dependence of the TDE, however, explains neither the increase in the number of defect atoms for tensile and compressive strains, nor the reduced sputter rate for (compressive and tensile) strained samples, Tab 1. For identical SBE, and by comparison to unstrained sample, the creation of damage and sputtered atoms should be higher in tension and lower in compression, where the TDE is lower and higher than in unstrained sample, respectively. This is, however, not the case. A more comprehesive study on the relation between TDE, SBE and strains has been recently published [START_REF] Guénolé | Influence of intrinsic strain on irradiation induced damage: the role of threshold displacement and surface binding energies[END_REF]. Nonetheless, the present results already clearly show that: i) the average TDE is not a good material parameter to characterize the tendency of a (strained) covalent material to form irradiation damage under multiple ion impact; ii) similarly, the widely considered notion, that the sputter yield is inversely proportional to the SBE [START_REF] Behrisch | Sputtering by Particle Bombardment: Experiments and Computer Calculations from Threshold to MeV Energies[END_REF] does not hold for strained silicon; and that iii) tensorial rather than scalar parameters are required to characterize or model irradiation damage in samples subjected to non-hydrostatic strains.

A general problem of the SBE and TDE for materials, which, like Si, show amorphization under irradiation is, that they are generally determined using crystalline samples. After the first hundred or so impacting ions, further ions strike already amorphized regions. Describing their interactions with the material, e.g. using the BCA, would in principle require values for the TDE and SBE determined in the amorphous state. As pointed out above, the correct density of the amorphous state should be considered in that case. This might, in particular, be of importance in the case of strained specimens. Here, the temporarily molten localized region can expand according to the constraints of the (strained) crystal surrounding the collision cascade. For instance, under tensile strain, the amorphous region can expand further and thus relax part of the stresses. The change of density and stress state within the irradiated region with strain could explain the change in elastic response for samples irradiated under strain, see Fig. 6).

Conclusions

We performed MD simulations of Ga + irradiation of silicon using two geometries: a slab with one (0 0 1) surface and normal incident beam and a model for FIB-milling of trenches with two intersecting {0 1 0} surfaces, where one surface is irradiated under a normal incident angle and the other at a grazing angle. Tensile and compressive uniaxial strains were applied in the latter geometry. Furthermore, the TDE and SBE for silicon were calculated for the first time as function of strain. The simulations revealed that:

• The amorphized regions in unstrained Si have a lower density than amorphous Si produced by the usual melt-quench route. This can be attributed to the high quench rates and the confining effect of the crystal surrounding the collision cascade. As a consequence, amorphous Si generated by other means than irradiation should only be used with care as model system for irradiation-induced amorphous Si.

• The extent of irradiation-induced damage and the sputter rate is influenced by (uniaxial) strains ≥ 2%. Unstrained samples show the lowest amount of damage and the highest sputter rate. Compressed samples exhibit the lowest sputter rates. Tensile strains lead to similar amount of damage as compressive strains and intermediate sputtering rates.

• Irradiation under strain affects the elastic response of the irradiated material stronger than in the unstrained case.

• The average TDE shows a linear relationship with hydrostatic strain, increasing for compression and decreasing for tension. The same trend is observed for uniaxial tension, but the effect is much smaller.

• The SBE is not influenced by tensile or compressive uniaxial strains < 4%. However, surface reconstruction can significantly increase the SBE. Empirical models therefore can underestimate the SBE of Si by upto a factor of two.

The strain dependence of the sputter rate and the extent of irradiation damage cannot be easily explained with the strain dependence of the TDE and SBE, the widely used key parameters for modeling irradiation damage, e.g., in the BCA. One shortcoming of these material parameters might be that they are usually determined for the crystalline state. Most of the sputtering and damage in Si, however, takes place in already damaged/amorphized regions. Advanced models of irradiation damage in Si might therefore benefit from calculations on amorphous systems with densities corresponding to irradiation-amorphized Si. Furthermore, tensorial rather than scalar material parameters will be required to model irradiation damage and sputtering in samples with non-hydrostatic strains.

The strain levels used in our simulations are comparable to residual strains, e.g., in multilayer thin films [START_REF] Korsunsky | Focused ion beam ring drilling for residual stress evaluation[END_REF][START_REF] Treml | High resolution determination of local residual stress gradients in single-and multilayer thin film systems[END_REF], nano-heterostructures [START_REF] Lu | In situ control of atomic-scale si layer with huge strain in the nanoheterostructure nisi/si/nisi through point contact reaction[END_REF], or strained silicon MOSFETs [START_REF] Chu | Strain: A solution for higher carrier mobility in nanoscale mosfets[END_REF]. The present findings should therefore have direct relevance to the modeling of FIB-machining of strained silicon structures, and more generally for the modeling of irradiation damages in materials showing irradiation-induced amorphization. There, improved models of irradiation damage under strain could contribute to the determination of crater functions, the prediction of surface stability and topographic pattern formation, as well as models for the strain relaxation in the FIB-DIC approach.

Figure 1 :

 1 Figure 1: Sketch and atomistic sample for (a) the 2-surfaces geometry and (b) the 1-surface geometry. (c) sample and 500 ions (green color) with a normal distribution of 4 nm emulating an ion beam focused at the edge. Atoms are colored according the diamond structure identification method implemented in Ovito (see Sec. 2.5): blue -diamond cubic, cyan -diamond cubic structure with 1st and or 2nd neighbor in non-diamond structure environment, e.g. surfaces, white -non-diamond structure.

Figure 2 :

 2 Figure 2: Collision cascades obtained by (a) BCA simulation with the TRIM package using MD-informed parameters, and (b) molecular dynamics simulation, in a Si target irradiated by Ga + at 5 kV. In (a) orange represents displaced atoms, green atoms are stopped atoms, white represents ion paths. In (b) only atoms in a non-perfect diamond cubic environment are shown, for the color code, see Fig. 1. (a) and (b) are displayed at equal scale. (c) distribution of the recoil atoms (TRIM), stable damage (MD) and implanted ions (TRIM and MD), for a dose of 500 Ga + at 5 kV on a Si target. The TRIM (MD based) results has been obtained by using the TDE and SBE values calculated in the present work.

Figure 3 :

 3 Figure 3: Induced damage in collision cascade simulations of 5 kV Ga + ions in Si at 300 K. (a) perspective view and (b) projected view of a slice with a thickness of 4 nm located at the simulation box center, i.e., with the highest concentration of incoming ions. (c) perspective view where atoms in a perfect diamond cubic structure environment (until the second neighbor) are removed to highlight the irradiation-induced damage. For the color code, see Fig.1.

Figure 4 :

 4 Figure 4: (a) radial distribution function for atoms in a spherical volume within the damaged region, both before (grey=crystalline) and after irradiation (blue).(b) distribution of atomic volumes inside this sphere, both before (grey=crystalline) and after irradiation (blue). As reference, also the RDF and the distribution of atomic volume of an amorphous sample generated by the melt-quench route is included (red). The corresponding atomic densities are also provided in the plots.

Figure 5 :

 5 Figure 5: Number of sputtered atoms as function of impacting ions for different values of uniaxial strain (compressive = negative, tensile=positive strain).

Figure 6 :

 6 Figure 6: Relative change in potential energy determined from uniaxial quasi-static tensile tests for the initial and irradiated samples (after 300 ions) at different pre-strains. The potential energy per atom is normalized to the one at the beginning of the tensile test. Note that in case of samples with ε < 0, the tensile tests are performed on compressed samples. The samples thus evolving toward the equilibrium volume, the potential energy exhibits a decreasing trend. Dose 0 300 Strain (%) --4 -2 0 +2 +4 Perfect atoms (#) 689,818 624,817 625,765 630,670 628,990 624,310 Sputtered atoms (#) 0 1,2151,360 1,546 1,452 1,408 Defect atoms (#) 25,962 89,748 88,655 83,564 85,338 90,062 Irradiation defect atoms (%) 0.0 9.1 9.0 8.3 8.5 9.2 S (atoms/impact) -4.4 4.7 5.3 5.0 5.0

Table 1 :

 1

	1,360	1,546	1,452	1,408

Number of atoms not characterized as belonging to any defect ('perfect atoms' ), sputtered atoms, atoms belonging to defects, fraction of atoms belonging to irradiation-induced damage and sputter yield S as function of uniaxial strain along the [010] direction.

Table 2 :

 2 [START_REF] Kang | Ultra-thin TEM sample preparation with advanced backside FIB milling method[END_REF]. Average threshold displacement energy TDE for silicon subjected to various hydrostatic (hydro.) and uniaxial [0 1 0] strains, calculated over a population N of randomly chosen directions. Reference values obtained from classical[START_REF] Pizzagalli | A new parametrization of the stillinger-weber potential for an improved description of defects and plasticity of silicon[END_REF] and ab-initio[START_REF] Holmström | Threshold defect production in silicon determined by density functional theory molecular dynamics simulations[END_REF] atomistic simulations are added.

	Strain	TDE (eV)	N
	0 %	38.3 ± 0.5	1000
	Ref.[31] 0 %	37.7 ± 0.4		-
	Ref.[41] 0 %	36.0 ± 2.0		-
	Hydro. -4 %	56.6 ± 1.7	216
	Hydro. -2 %	46.7 ± 0.9	457
	Hydro. +2 %	25.3 ± 0.2	1919
	Hydro. +4 %	17.7 ± 0.1	1912
	[0 1 0] -4 %	39.6 ± 0.5	1190
	[0 1 0] -2 %	39.3 ± 0.4	1453
	[0 1 0] +2 %	37.7 ± 0.5	1000
	[0 1 0] +4 %	34.3 ± 0.4	1000
	Strain	SBE (eV ±0.01)
		As-cut Reconstructed
	Ref.[44] 0 %	4.73	-
	-4 %	6.97	10.59
	-2 %	6.97	10.61
	0 %	6.97	10.63
	+2 %	6.97	10.65
	+4 %	6.93	10.63

Acknowledgments

The authors thank Christoph Breuning for preliminary calculations of threshold displacement energies. EB acknowledges the support of the Cluster of Excellence Engineering of Advanced Materials (EAM) and the Central Institute of Scientific Computing (ZISC). This work was supported by the European Union, within the large collaborative project ISTRESS, Grant Agreement No. 604646.