CC BY
0Contents lists available at SciVerse ScienceDirect
Nuclear Instruments and Methods in Physics Research B
journal homepage: www.elsevier.com/locate/nimb
Hybrid continuum-atomistic modelling of swift heavy ion radiation damage in germanium
S.L. Daraszewicz, D.M. Duffy *
Department of Physics and Astronomy and The London Centre for Nanotechnology, University College London, Cower Street, London WC1E 6BT, United Kingdom
ARTICLE INFO
Article history:
Received 19 July 2012
Received in revised form 14 November 2012
Available online 31 December 2012
Keywords: Swift heavy ions Germanium Inelastic Thermal spike Radiation damage Molecular dynamics
ABSTRACT
The response of germanium to swift heavy ion irradiation is simulated using a hybrid continuum-atomistic approach. The continuum part of the model, which characterises the electronic excitations is an extension of the inelastic thermal spike based on an approximation to the Boltzmann transport equation; while the atomistic part is represented with molecular dynamics. This integrated method can realistically account for the non-equilibrium carrier dynamics in band-gap materials under irradiation, unlike earlier developments based on the two-temperature approach. The model is used to obtain temporal and spatial evolution of carrier density, electronic temperature and lattice temperature for germanium irradiated with carbon cluster ions. Good agreement with experimental data of amorphised latent track radii for different stopping powers is obtained by fitting a constant value for the electron-phonon coupling strength - the only parameter treated as free in the model.
© 2013 Elsevier B.V. All rights reserved.
1. Introduction
The physical mechanisms behind the ion track creation have been a subject of a debate since the late 1950s, when the first experimental evidence for the cylindrical-shape damage were produced [1,2]. It was later realised that such damage is caused by huge electronic excitations, rather than by direct projectile-target collisions. Such excitations indirectly lead to changes in the atomistic structure through the electron-phonon coupling mechanism. The energy losses in ion-matter interaction are quantified by the so-called stopping power, S = @x, a measure of the energy loss per unit distance. For the fullerene beam bombardment considered here, the energy of the projectile is lost mainly to the inelastic collisions with the electrons (electronic stopping, Se) and this energy dissipation channel dominates over the elastic collisions.
Ion tracks were originally created in insulators, however nowadays with increasing beam energies reaching the MeV regime, these can be generated in metals [3,4], alloys [5] and semiconductors [6,7]. Variations of the thermal spike model [8,9] are typically invoked when attempting to explain the mechanisms of track formation. The model assumes that two distinct temperatures can be assigned to describe the energy of ionic (Tl) and electronic (Te) subsystems. Typically, both temperatures are evolved by a heat diffusion equation and coupled via a term proportional to the temperature difference between the electrons and ions. Such
* Corresponding author. E-mail addresses: szymon.daraszewicz.09@ucl.ac.uk (S.L. Daraszewicz), d.duf-fy@ucl.ac.uk (D.M. Duffy).
0168-583X/$ - see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.nimb.2012.11.027
two-temperature (2T) formulation was also successfully applied to other ultrafast phenomena, such as materials response to laser excitation, in both continuum and more recently in continuum-atomistic frameworks. The thermal spike was primarily applied to metals and its application to band-gap materials (both insulators [10,11] and semiconductors [12]) was subject to criticism [13-15]. Following an earlier development in ultrafast dynamics in lasers [16] (and more recently in [17]), we adapted an energy transport model based on an approximation to the Boltzmann transport equation [18] to correctly account for the dynamics of excited carriers in band-gap materials [14,15]. This model is linked to molecular dynamics (MD) providing us with an improved description of thermomechanical properties and detailed atomistic evolution of the irradiated system. In this work, we attempt to describe the behaviour of germanium under swift heavy ion (SHI) irradiation using this integrated framework.
2. The hybrid model
2.1. Extending the inelastic thermal spike model
In the extended inelastic thermal spike model [15], the carrier (electron-hole) density (N) evolution is described by
— + V-J = Ge - Re (1)
with the carrier current (J) defined as
J= -D(VN+weVEg+¿S • (2)
The current pre-factor (D) is the ambipolar diffusion coefficient, following the assumption that the model is locally charge neutral. A similar development, which accounts for decoupled electrons and holes, can be found in [13]. In Eq. (1) the carrier generation (Ge) is related to the initial energy deposited by the ion track and impact ionisation, while (Re) represents Auger recombination or other sink terms [14]. The total carrier energy (U) is composed of a potential and a kinetic part U = N(3kBTe + Eg), which leads to the carrier specific heat capacity of Ceh = 3NkB + N@T-. Here Eg, kB denote the bandgap and the Boltzmann constant, respectively. In the extended model, the ambipolar energy current is given by W = (Eg + 4kBTe) - (Ke + Kh)VTe, where Ke, Kh are the electron and hole thermal conductivities, respectively. This current consists of two terms: one relating to the carrier flow and the other one dependent on the electron-hole temperature gradient. The rate equation of the carrier energy can be expressed as
@UU + V • W = A(r[v], t)-Uep• (3)
In our simulations, the source term A(r[v], t) is associated with the SHI electron-hole generation and Uep is a sink term coupled to the ions, defined as Uep = Ct (Te - Tl) with sep being the characteristic electron-phonon relaxation time constant. This represents the only parameter in the model that is treated as free, as the remaining ones are relatively well-known for semiconductors (following the assumption that these would remain the same under SHI irradiation conditions). The model parameters for the electronic subsystem are presented in Table 1. Note that we relate the thermal carrier conductivities to the experimentally measured carrier mobility values via a relationship suggested in [13] (Eqs. 19-20). Furthermore we rescale the carrier mobility for elevated electronic temperatures according to [19,20].
The source term in the Eq. (3) corresponds to a C60 (carbon cluster) irradiation of a specific ion energy of 0.07 MeV/u. The mean absorption radius (rmean) can be obtained by applying Bohr's principle of adiabatic invariance [13,21]. For non-relativistic ions it can be expressed as
rmean = hv/2Eg, (4)
where v is the velocity of the imparting ion and h is the reduced Planck constant. We assume a spatial Gaussian deposition of r = rmean and we choose A(r[v], t) = AD(r)a-at with a = 1/s, where s is a characteristic deposition time taken to be 1 fs [22,23]. Normalisation of A(r[v],f) is constrained so that the spatial and temporal integration equal the assumed deposited electronic stopping power (Se); in spherical coordinates
{•5s /*rmax
Se = / A(r[v], t)2nrdrdt, (5)
Jt=0 Jr=0
where rmax is the maximum range of carriers projected in the lateral direction and is assumed to be 5r. A more accurate method of calculating the exact energy deposited in the electronic structure can be made based on calculations found in [24], which were described
Table 1
Model parameters for the electronic subsystem of germanium.
Property Symbol Germanium
Ambipolar diffusivity [17] D 65 •(T,/300K)"1'5 (cm2/s)
Band-gap [17] Eg 0.803 - 3.9 • 10-4T (eV)
Auger recombination coeff. [17] c 2.0 • 10-31 (cm6/s)
Electron mobility at 300 K [26] le 63900 (cm2 V-1 s-1)
Hole mobility at 300 K [26] 1h 61900 (cm2 V-1 s-1)
more recently in [25]. In this model we assume that a fraction of the energy is required to create carriers (NEg) and the remaining part goes to their kinetic energy (3NkBTe).
2.2. Hybrid continuum-atomistic approach
A hybrid continuum-atomistic approach allows us to access atomistic behaviour via trajectory of the atoms, which has several advantages to continuum models. It was shown that in an SHI irradiation scenario the ion temperature can reach above the melting temperature with no track formation [27]. Also, continuum-only models typically neglect lattice straining and emission of shock waves and do not take into account the volume change induced by a phase-transition. More sophisticated continuum-only methods, which incorporate the transport of momentum and mass, exist [28], however the above effects are naturally included in MD at the level of the interatomic potential. We also note that a similar integrated model for band-gap materials was recently employed in laser irradiation studies [29,30].
The coupling to molecular dynamics is performed by replacing the equation for the lattice temperature found in [15] with modified MD equations of motions, according to [31-33]:
m = F i(t)-CiVi + Fi(t), (6)
where i runs over all atoms and Fi is the Newtonian force on an atom, while F~i and civi are additional driving and friction forces, based on an Langenvin thermostat formulation, that incorporate the electron-phonon coupling. The energy gain represented by a stochastic force Fi(t) has a random magnitude and orientation with Fi(t) = y/TAi(t), where by construction r = 2yimikBTe with yi being the friction term. The main difference from the original formulation of the combined 2T+ MD [32,33] is that the energy exchange form of Uep = Ceh (Te - Tl), forces the electron-phonon friction (yep = 1 /y) to be rescaled according to %/ep = yep(N/Na) at each iteration timestep in order to conserve the total energy, where N(Na) is the number of carriers (atoms) in a coarse-grained temperature cell.
For simplicity the energy exchange is only permitted in one way: from the carriers to the ions. Earlier calculations, allowing for a two-way exchange showed that less than 1% of the total energy transfer goes back to the carriers. This so-called cooling inhibition [12] represents a physical scenario, where the carriers had already diffused away and therefore there is minimal coupling between the hot ions and the carriers. This leads to a fundamental difference in simulations of the band-gap and metallic materials, where in the latter case the energy from the hot lattice can be exchanged back with the electronic subsystem to subsequently diffuse away.
The hybrid simulation setup is as follows. We have used a modified version of the DL_POLY (4.01) code [34] to perform all the calculations. Due to the cooling inhibition we have used stochastic boundary conditions at 300 K in the lateral direction to the ion excitation to create a mechanism for energy dissipations by pho-nons. The MD cell boundary conditions are periodic to represent a bulk condition. The simulations were performed with a constant time-step of 1 fs in the MD part of the code, having checked that it conserves energy under such strong non-equilibrium conditions. The corresponding electronic structure solver has semi-open (Robin's) boundaries in the lateral direction to represent energy dissipation into the bulk and Neumann boundary conditions applied to the sides perpendicular to the SHI impact direction. We used a cell containing 200 k Ge atoms, with a size of 28.36, 28.36, 5.67 nm and corresponding 25 x 25 x 5 coarse grained temperature cells overlay. The electronic temperature cells extended over the MD system in the lateral directions and consisted of 50 x 50 x 5 grid points. This simulation configuration is schematically presented is Fig. 1.
Swift heavy ion
Ion temperature cells
Fig. 1. Simplified simulation set-up. The energy from the imparting SHI projectile is initially given to the electronic structure grid points and is exchanged with the corresponding ion temperature cells. The energy can diffuse away since the grid points extend over the MD simulation cell and semi-open boundaries are applied to the edge cells in the lateral directions. Stochastic boundary conditions are applied in the lateral directions of the MD simulation cell.
The total simulation time was 50 ps. The atomistic configurations were pre-equilibrated in an NPT ensemble. We employ a Tersoff potential for germanium [35].
In the continuum part of the model, the spatial and temporal evolutions of carrier energy (1) and density (2) were solved using
a forward in time, centred in space explicit time-stepping (Euler) method with a time-step of ~1as, which corresponds to 1000 steps over one MD iteration. The continuum solver was disabled, once the majority of the energy was transferred to the ions, which occurred typically at ~1 ps. An alternative semi-implicit algorithm was recently proposed in [29].
We neglect athermal effects, such as the impact of the elevated electronic temperature on the interatomic potential. Furthermore, the model is locally charge neutral meaning that the formation of an electrostatic field (i.e. Coulomb explosion) at the initial stages of the irradiation event is neglected. The thermalisation time of electrons is not taken into account and it is assumed that the subsystem can exchange energy immediately, i.e. as early as during the impact of an SHI. Furthermore, owing to the complexity of the model, we neglect any variations of the band-gap, which implies that the carrier confinement (discussed in [36]) cannot occur. We also recognise the limitations of the interatomic potential employed, which overestimates the melting temperature of germanium [37].
3. Results and discussion
Using the hybrid continuum-atomistic model, we performed simulations of germanium for various electron-phonon relaxation
~(-t )
i*|(b)I
vXWX-'SvXw/'
0.05 0.075 0.1 0
O experiment
--t = 0.05ps
t = 0.10ps
--t =0.15ps JI
, , , , fP, , , , , r.I.f. 7
.v.v.v.wWw.v.v.v.
10 20 30 40 50 Stopping power [keV/nm]
Fig. 2. (a) Calculated track radii versus electronic stopping power for sep = 0.05,0.10 and 0.15 ps (dashed lines). The relaxation time of sep = 0.05ps is in a reasonable agreement with sparse experimental data (0) for carbon cluster ions [39]. The inset shows the dependence of the track radii on sep for Se = 50 keV/nm; (b) an exemplar track created at 35 keV/nm (xy-plane slice). Latent track radii are well-defined and are determined by visual inspection.
20 25 30 Time [ps]
Fig. 3. Evolution of the lattice temperature for a sample sep = 0.10 ps simulation at 50 keV/nm at different radii from the impact centre. This SHI irradiation event produced a latent track with a 2 nm radius. The inset shows the evolution of the electronic temperature and the carrier density for the first 1 ps of the simulation. The Auger recombination (predominantly) and ambipolar diffusion are responsible for the rapid carrier density decrease during the first 1 ps after the SHI impact.
times (sep = 0.05—0.15 ps) for stopping powers in the range of 10- [2 70 keV/nm. We achieved a reasonable agreement with limited [3 experimental data for sep = 0.05 ps - a relaxation time value lower than the one reported in the literature (cf. sep ~ 0.3 ps [38 ] ). An [4 exemplar latent track is presented in Fig. 2(b). Lattice temperature and the corresponding carrier density and temperature temporal evolution - at different radial distances from the track centre - [6 for a typical simulation are presented in Fig. 3.
The impact of the collapse of the band-gap in the molten region on carrier diffusivity (i.e. carrier confinement) could enhance the coupling between the two subsystems leading to formation of lar- [8 ger ion tracks. Preliminary continuum-only results for silicon indi- [9 cated that such effect could lead to 10-20% increase in track sizes. [10 We have additionally performed basic parameter sensitivity analysis, which showed that the model is quite sensitive to the choice of the ambipolar diffusivity (a factor of 0.5 decrease in D, results in an ~1.5 factor increase in ion track radii) and relatively insensitive to [12 the variations in the remaining input parameters. Furthermore, the melting temperature could decrease as a function of the electronic [13 temperature due to athermal changes in the interatomic potential [14 (as for silicon [40 ] ) and thus the melt track radii could be underestimated for such strong electronic stopping powers considered here. Results for silicon under heavy electronic excitations generated by SHI irradiation, which include the athermal effects, would be presented elsewhere.
4. Conclusion
References
[1] D.A. Young, Nature 182 (1958) 375-377.
[18 [19
Based on the premise that the two-temperature and continuum-only description of the behaviour of band-gap materials under SHI irradiation is incomplete, we built a hybrid continuum-atomistic approach based on the extended inelastic thermal spike model. We presented exemplar simulations for germanium, which showed that the model is consistent with the experimentally seen amorphous radii of latent tracks produced by SHIs using one fitting [24 parameter - the electron-phonon relaxation time. The model provided us with temporal and spatial evolutions of the ion temperature, carrier density and energy. We believe that this framework [26 can be applied to lower impact energy regimes, combining both the inelastic and elastic projectile-target scattering types and, in the future, it could employ new potentials parameterised with respect to the electronic temperature to distinguish between the thermal and athermal processes.
[28 [29 [30 [31
Acknowledgements [32
The authors acknowledge the use of the UCL Legion High Performance Computing Facility, and associated support services, in the [34 completion of this work. S. Daraszewicz acknowledges funding [35 from EPSRC under the Materials Modelling and Materials Science ^
Industrial Doctoral Training Centre and from the Culham Centre [37 for Fusion Energy (CCFE). [38
E.C.H. Silk, R.S. Barnes, Philosophical Magazine 4 (1959) 970-972.
A. Dunlop, D. Lesueur, P. Legrand, H. Dammak, J. Dural, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 90 (1994) 330-338.
H. Dammak, A. Dunlop, D. Lesueur, A. Brunelle, S. Della-Negra, Y. Beyec, Physical Review Letters 74 (1995) 1135-1138.
C. Trautmann, S. Andler, W. Bruchle, R. Spohr, M. Toulemonde, Radiation Effects and Defects in Solids 126 (1993) 207-210.
A. Dunlop, G. Jaskierowicz, S. Della-Negra, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 146 (1998) 302-308.
F. Komarov, P. Gaiduk, L. Vlasukova, A. Didyk, V. Yuvchenko, Vacuum 70 (2003) 75-79.
M. Toulemonde, C. Dufour, E. Paumier, Physical Review B 46 (1992) 1436214369.
G. Szenes, Physical Review B 51 (1995) 8026-8029.
A. Meftah, J. Costantini, N. Khalfaoui, S. Boudjadar, J. Stoquert, F. Studer, M. Toulemonde, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 237 (2005) 563-574. M. Toulemonde, W. Assmann, C. Dufour, A. Meftah, F. Studer, C. Trautmann, Matematisk-Fysiske Meddelelser 52 (2006) 263-292.
A. Chettah, H. Kucal, Z. Wang, M. Kac, A. Meftah, M. Toulemonde, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 267 (2009) 2719-2724. S. Klaumunzer, Matematisk-Fysiske Meddelelser 52 (2006) 293-328.
D. Duffy, S. Daraszewicz, J. Mulroue, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 277 (2012) 21-27.
S. Daraszewicz, D. Duffy, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 269 (2011) 1646-1649.
H. van Driel, Physical Review B 35 (1987) 8166-8176.
J. Chen, D. Tzou, J. Beraun, International Journal of Heat and Mass Transfer 48 (2005) 501-509.
J. Drabble, H. Goldsmid, Thermal Conduction in Semiconductors, Pergamon Press, Oxford, 1961.
T. Grasser, T. Tang, H. Kosina, S. Selberherr, Proceedings of the IEEE 91 (2003) 251-274.
G. Baccarani, M. Wordeman, Solid-State Electronics 28 (1985) 407-416.
A. Mozumder, The Journal of Chemical Physics 60 (1974) 1145.
B. Gervais, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 88 (1994) 355-364.
M. Toulemonde, J. Costantini, C. Dufour, A. Meftah, E. Paumier, F. Studer, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 116 (1996) 37-42. M. Waligorski, R. Hamm, R. Katz, Nuclear Tracks and Radiation Measurements 11 (1986) 309-319.
C. Dufour, V. Khomenkov, G. Rizza, M. Toulemonde, Journal of Physics D -Applied Physics 45 (2012).
C. Jacoboni, F. Nava, C. Canali, G. Ottaviani, Physical Review B 24 (1981) 10141026. Parameters for Germanium collected at: http://www.ioffe.ru/SVA/NSM/ Semicond/Ge/electric.html.
D.M. Duffy, N. Itoh, A.M. Rutherford, A.M. Stoneham, Journal of Physics: Condensed Matter 20 (2008) 082201.
M. Jakas, E. Bringa, R. Johnson, Physical Review B 65 (2002).
Y. Gan, J. Chen, Computer Physics Communications 183 (2012) 278-284.
Y. Gan, J.K. Chen, Applied Physics A 105 (2) (2011) 427-437.
A. Caro, M. Victoria, Physical Review A 40 (1989) 2287-2291.
D.M. Duffy, A.M. Rutherford, Journal of Physics: Condensed Matter 19 (2007)
016207.
A.M. Rutherford, D.M. Duffy, Journal of Physics: Condensed Matter 19 (2007) 496201.
I.T. Todorov, W. Smith, K. Trachenko, M.T. Dove, Journal of Materials Chemistry 16(2006) 1911.
J. Tersoff, Physical Review B 39 (1989) 5566-5568.
S.K. Sundaram, E. Mazur, Nature Materials 1 (2002) 217-224.
S. Cook, P. Clancy, Physical Review B 47 (1993) 7686-7699.
M. Gallant, H. van Driel, Physical Review B 26 (1982) 2133-2146.
A. Colder, O. Marty, B. Canut, M. Levalois, P. Marie, X. Portier, S. Ramos, M.
Toulemonde, Nuclear Instruments and Methods in Physics Research Section B:
Beam Interactions with Materials and Atoms 174 (2001) 491-498.
L. Shokeen, P.K. Schelling, Journal of Applied Physics 109 (2011) 073503.
