Coulomb Spike Model of Radiation Damage in Wide Band-Gap Insulators

Abstract

A novel Coulomb spike concept is applied to the radiation damage induced in LiF and SiO₂ with about the same mass density (~2.65 g cm⁻³) by ${}_{28}{}^{60}\mathrm{N}\mathrm{i}$ and ${}_{36}{}^{84}\mathrm{K}\mathrm{r}$ ions of 1.0-MeV u⁻¹ energy for about the same electronic energy loss (~10 MeV µm⁻¹). This is an alternative concept to the already known models of the Coulomb spike and inelastic thermal spike for the damage induced by swift heavy ion irradiations. The distribution of ionizations and electrostatic energy gained in the electric field by the ionized atoms is computed with the PHITS code for both targets. Further, the atomic collision cascades induced by these low-energy hot ions of about 500 eV are simulated with the SRIM2013 code. It is found that melting is reached in a small volume for SiO₂ due to the energy deposition in the subthreshold events of nuclear collisions induced by the Si and O ions. For LiF, the phonon contribution to the stopping power of the lighter Li and F ions is not sufficient to induce melting, even though the melting temperature is lower than for SiO₂. The formation of amorphous domains in SiO₂ is likely after fast quenching of the small molten pockets, whereas only point defects may be formed in LiF.

1. Introduction

The modeling of radiation damage in insulating materials by charged particles is a key issue for nuclear applications in order to understand the modifications induced by fission fragments in the nuclear fuels [1] or in the inert matrices for actinide transmutation [2]. It is also a matter of concern for space applications regarding the behavior of the SIMOX (separation by implanted oxygen) electronic devices containing buried silica layers under cosmic-ray irradiations with high-energy protons [3].

Various models and numerical simulations were developed to comprehend the evolution of wide band-gap insulators under swift ion irradiations, such as fission fragments or high-energy cosmic rays, in the electronic energy loss regime. In particular, the inelastic thermal spike (i-TS) model was applied to explain the track formation in oxides or halides [4]. The i-TS model, based on the two-temperature (TT) hypothesis, was first devised for the case of the laser irradiation of metals as a means of heat treatment, eventually leading to melting or vaporization [5]. It assumes that the energy deposited into the electronic subsystem on a short time scale (≤1 fs) by photons or charged particles is transferred to the lattice atoms by the electron–phonon coupling interaction on a longer time scale (≥1 ps) [6]. For the single excitations with photons, another approach was based on Boltzmann’s kinetic equations [7]. For the collective excitations with heavy ions, the lattice can reach the melting or vaporizing temperature above a given threshold electronic stopping power [4], and extended defects, such as ionization tracks, can be induced in the solid matrix after rapid cooling of the melt or vapor. Melting or vaporization is reached when the energy deposited in the solid is equal to the enthalpy of melting or vaporization, respectively. This approach was criticized by some authors claiming that the actual reached temperatures in metals are much lower than the computed values by the i-TS model [8]. An alternative approach was also proposed by these authors with the Monte Carlo (MC)–Molecular dynamics (MD) TREKIS model [9]. Moreover, there was also a recent development of the TT model by using the electronic excitations as an input for the MD simulations of metals such as tungsten under swift heavy ion irradiations [10,11]. In contrast, the original i-TS model is not relevant for metals and semiconductors owing to the high electronic thermal conductivity quickly spreading out the deposited energy in the electron gas.

Another type of model is based on the Coulomb spike or Coulomb explosion concept that was proposed a long time ago by Fleischer et al. [12]. This approach is generally disregarded nowadays. An analytical model of the Coulomb explosion was developed to interpret radiation damage, as well as phase changes, in metals such as titanium under swift heavy ion irradiations [13,14], which could not be readily explained by the standard i-TS model, as said above.

Among many other insulators, amorphous tracks are formed in SiO₂ quartz by swift heavy ion irradiations above a threshold electronic stopping power of ~2 MeV µm⁻¹ and 5 MeV µm⁻¹ for ion velocities of 0.2 and 5 MeV u⁻¹, respectively [15,16]. In contrast, for LiF, only saturated point defects (F centers) are induced in crystalline tracks up to the recombination volume [17]. A large swelling is actually induced in LiF above a threshold stopping power of 4.2 MeV µm⁻¹ [18]. The i-TS model was unable to explain the formation of such defective crystalline tracks in alkali halides. It is based on the strong assumptions explained above, which may not be relevant in the case of such wide band-gap ionic insulators. An approach tried to extend the self-trapped exciton (STE) model of point defect formation by the laser excitation of insulators such as alkali halides to the track formation by swift heavy ions [19]. However, this approach remained at a qualitative level.

Therefore, our scope was to revisit the Coulomb spike model of radiation damage by tracking down the ionized atoms induced by the passage of ${}_{28}{}^{60}\mathrm{N}\mathrm{i}$ and ${}_{36}{}^{84}\mathrm{K}\mathrm{r}$ ions with the 1.0-MeV u⁻¹ energy, corresponding to the fission-fragment energy per unit mass, in these two well-studied wide band-gap insulators, namely LiF and SiO₂. Either crystalline or amorphous tracks are definitely formed in both materials with the present Ni and Kr ions for stopping powers (Table 1) higher than the respective thresholds [16,17,18].

However, the effect of the electrostatic energy induced by swift heavy ions among atoms on a short time scale (<1 ps) was disregarded. Electrostatic charge effects cannot be overlooked in such high-k dielectrics. However, we did not follow the same line as for the Coulomb explosion model by Lesueur et al. [13]. Simulations with the MC RITRACKS computer code have already been used to model the jet-like sputtering of amorphous ice induced by 1.0-MeV protons [20]. In the latter case, a large part of the ionized atoms with high kinetic energies did not escape from the target surface but remained trapped inside the bulk.

We show that the Si and O ions generated in the trails of swift ions can induce the melting of SiO₂ in a very small volume by nuclear collisions, due to the subthreshold elastic events. Amorphous domains can be generated by the quenching of the molten solid after this elastic thermal spike. We have applied the same criterion of melting as for the inelastic thermal spike that melting is reached when the energy deposited in the solid is equal to the fusion enthalpy. By contrast, melting is not reached in LiF by the nuclear collisions induced by the lighter Li and F ions. Only some point defects may be formed after thermal quenching of the solid phase.

Table 1. Characteristics of the ion irradiations of LiF and SiO₂ computed with the SRIM2013 code [21]: electronic (S_e) and nuclear (S_n) stopping power, mean projected range (R_p), range straggling (∆R_p), and number of vacancies per ion.

Target	Ion	E (keV)	Se (keV µm−1)	Se (*) (keV µm−1)	Sn (keV µm−1)	Rp (nm)	∆Rp (nm)	Vacancies
LiF	Li	0.5	12.89		75.26	4.0	2.1	4.8
	F	0.5	17.19		250.1	2.0	0.9	5.5
	Ni	6 × 104	9.34 × 103	10.5 × 103	21.45	11.81 × 103	448
	Kr	8.4 × 104	1.17 × 104	1.32 × 104	36.68	11.93 × 103	399.5
SiO2	Si	0.5	24.02		281.9	2.0	0.9	7.5
	O	0.5	21.97		176.7	2.4	1.3	7.0
	Ni	6 × 104	8.52 × 103	9.62 × 103	24.31	11.86 × 103	489
	Kr	8.4 × 104	1.07 × 104	1.23 × 104	41.64	12.04 × 103	435.5

(*) Computed with the CASP 5.0 code [21].

Table 1. Characteristics of the ion irradiations of LiF and SiO₂ computed with the SRIM2013 code [21]: electronic (S_e) and nuclear (S_n) stopping power, mean projected range (R_p), range straggling (∆R_p), and number of vacancies per ion.

Target	Ion	E (keV)	Se (keV µm−1)	Se (*) (keV µm−1)	Sn (keV µm−1)	Rp (nm)	∆Rp (nm)	Vacancies
LiF	Li	0.5	12.89		75.26	4.0	2.1	4.8
	F	0.5	17.19		250.1	2.0	0.9	5.5
	Ni	6 × 104	9.34 × 103	10.5 × 103	21.45	11.81 × 103	448
	Kr	8.4 × 104	1.17 × 104	1.32 × 104	36.68	11.93 × 103	399.5
SiO2	Si	0.5	24.02		281.9	2.0	0.9	7.5
	O	0.5	21.97		176.7	2.4	1.3	7.0
	Ni	6 × 104	8.52 × 103	9.62 × 103	24.31	11.86 × 103	489
	Kr	8.4 × 104	1.07 × 104	1.23 × 104	41.64	12.04 × 103	435.5

(*) Computed with the CASP 5.0 code [21].

2. Methodology

In order to calculate the electrostatic energy of ions on the trajectories of projectile ions, a calculation procedure similar to that developed in [20] was used. The ionization of the constituent atoms in the targets was calculated by ITSART [22], one of the track-structure models of PHITS [23], for one primary particle history. Even though the calculation is on one history, the electrostatic energy calculated afterward is based on the distance between several ten thousand cations which is statistically significant. Track structure models allow for the calculation of the transport of charged particles by taking into account each atomic reaction explicitly, unlike the conventional condensed history method or the continuous slowing-down approach. When the projectiles undergo ionization reactions in the target, the spatial coordinates of the reactions are scored at each time. The secondary electrons above 1 keV and those below 1 keV were transported by using the EGS5 and ETSART [24] models, respectively. Using the recorded coordinates as the position of the cations, the electrostatic potential among them was calculated offline. This calculation was performed for LiF and SiO₂ irradiated with 60-MeV ⁶⁰Ni and 84-MeV ⁸⁴Kr ions.

Energetic ionized Li and F or Si and O atoms are produced in the trails of these swift ions that cross the target thickness of 2 µm. The parameters, such as electronic (S_e) and nuclear (S_n) stopping powers, projected range (R_p), and range straggling (∆R_p), are computed with the MC SRIM2013 code [25] for all these ions in LiF and SiO₂ (Table 1). In this range of kinetic energies, SRIM is more accurate than PHITS. The radial distributions of the secondary electrons (δ-rays) were computed for both Kr and Ni ion irradiations in LiF and SiO₂ targets (Figure 1). These electrons were knocked out by the primary ion or by energetic secondary electrons and transported by either EGS or ETSART codes until they slowed down to the cut-off energy of 10 eV up to a time scale of 150 fs (Figure 2). The steps found for the Ni ion irradiation of SiO₂ derive from the ionizations of the various electronic shells. A quadratic dependence in r⁻² is found for large radii (r), as found by the older analysis of [25]. The 3D distributions of cations computed for one history in these four cases are shown by the denser points corresponding to the ionizations induced along the ion trajectories and by discrete points corresponding to the ionization events arising from secondary electrons or photons far away from the ion path (Figure 3a–d). The primary ions are deflected by elastic Rutherford nuclear scattering and inelastic electron scattering. Almost all of the ionization events are localized on the ion path.

There is the challenging issue of a possible cancellation of the electric potential between ionized atoms by the major part of δ-rays in the track core (Figure 1). However, at the short time scale (fs) of the electro-magnetic interactions, these ejected electrons have still kinetic energies above the band-gap energy (E_G) of about 10 eV for both materials (E_G = 14 eV for LiF and E_G = 9.5 eV for α-SiO₂). They are in the conduction band in quantized plasmonic states (Figure 2), where the behaviors of electrons and nuclei can be decoupled in the adiabatic approximation. In the classical approach, the electric field is partially screened by the electron gas with the corresponding dielectric constant at the plasmon resonance frequency ε(ω_p). For SiO₂, the oxygen bulk plasmon energy is ħω_p = 23.1 eV [26], for the frequency 5.58 × 10¹⁵ Hz and period 0.179 fs, whereas for LiF, ħω_p = 25.3 eV [27]. For those insulators, the real (ε₁) and imaginary (ε₂) parts of ε(ω) tend towards the free-space permittivity (ε₀) for high frequencies [28]. The plasmon resonance occurs for ε₁ (ω_p) = 0 [27]. Since ħω_p >> E_G the electrons are quasi-free for this energy. As such, the screening by these free electrons is negligible on the fs time scale. Plasmon decay to single excitations (electron-hole pairs) takes place later on the 10 fs time scale (Figure 2).

The higher energy electrons are flying far away from the ion path (Figure 1), yet they are not likely trapped back onto the ionized atoms at the time scale of the interactions owing to the low mobility of carriers in these dielectric media. For a standard electronic diffusivity of D = 0.1 cm² s⁻¹ for trapped electrons, an upper-boundary estimation of the mean free path of the stopped electrons would be of $\mathsf{\lambda }=\sqrt{\mathrm{D}∆\mathrm{t}}$ ~10⁻⁷ cm = 1 nm for ∆t = 100 fs after the electron cascade cooling on the 150 fs time scale (Figure 2). These secondary electrons may induce a few negligible atomic displacements by elastic collisions for the most energetic ones that are far away from the track cores.

In order to estimate the damage induced by these trapped energetic ions with a very small range, SRIM2013 simulations were performed for 10⁵ ions by using the following threshold-displacement energies of E_d (Li) = E_d (F) = 25 eV for LiF (mass density of ρ₀ = 2.64 g cm⁻³) and E_d (Si) = 15 eV, E_d (O) = 28 eV for SiO₂ (ρ₀ = 2.65 g cm⁻³). The relative contributions of ionization, vacancies, and phonons of the energy loss were obtained for the primary and secondary ions and their recoils.

3. Results

The histograms of electrostatic energy distributions of ionized Li and F and Si and O atoms are shown for LiF (Figure 4 and Figure 5) and SiO₂, respectively (Figure 6 and Figure 7). This is the energy gained in the electrostatic field generated by the incoming swift ions in the targets. For LiF, similar distributions are fitted with three broad Gaussian profiles centered at 180, 450, and 610 eV for the Ni ions and 180, 480, and 660 eV for the Kr ions. The low energy events below 100 eV, corresponding to strong statistical fluctuations, are discarded. For SiO₂, a most prominent Gaussian profile centered at 620 eV and 455 eV is found for the Ni and Kr ions, respectively, with much higher intensities than for LiF and minor contributions centered at 550 and 350 eV, respectively, as deduced from fits of the distributions. The band centers and standard deviations of the Gaussian profiles are summarized in

Table 2. Fitted Gaussian profiles of the energy distributions of ionized atoms computed with the PHITS code [23].

Target	Ion	Center (eV)	Standard Deviation (eV)	Center (eV)	Standard Deviation (eV)	Center (eV)	Standard Deviation (eV)
LiF	Ni	180	80	450	100	610	75
	Kr	180	80	480	100	650	45
SiO2	Ni	550	40	620	40
	Kr	350	50	455	60

.

Maxwellian distributions of these ions are computed to approximately match the maxima of these distributions (Figure 4, Figure 5, Figure 6 and Figure 7). It is seen that the ion energy distributions are much narrower than the broad thermal equilibrium distributions for the corresponding high temperatures. This means that these energetic atoms are not in thermal equilibrium at the end of the ionization process.

The projected range and number of vacancies per ion are deduced from the SRIM2013 simulations of atomic collision cascades for a typical kinetic energy of 500 eV of the incident ions (Table 1). The relative contributions in % of ionization, vacancies, and phonons for the primary ions and recoils are also given (

Table 3. Results of the quick damage simulations of ion irradiations computed with the SRIM2013 code [29].

Target	Ion	Primaries Ionization (%)	Vacancies (%)	Phonons (%)	Recoils Ionization (%)	Vacancies (%)	Phonons (%)
LiF	Li	10.2	1.8	19	11.4	0.6	57
	F	5.1	1.9	16.6	13.1	0.8	62.5
	Ni	99.73	/	0.01	0.17	/	0.09
	Kr	99.60	/	0.01	0.26	/	0.11
SiO2	Si	7.6	1.9	12.1	12.5	1.4	64.5
	O	10.3	1.8	13.0	11.8	1.25	61.8
	Ni	99.66	/	0.01	0.20	0.01	0.12
	Kr	99.55	/	0.01	0.28	0.01	0.16

). It is clearly seen that the major part (~80%) of the ion energy is dissipated in phonons and mainly by the recoils for energies below E_d. Very few atoms (~7 per incident ion) are displaced by the nuclear collisions. A smaller part (~20%) of the total stopping power is dissipated in ionizations and electronic excitations.

The energy deposited by the primary swift Ni and Kr ions is of δE ~20 MeV for both ions inside the 2 µm target thickness (Table 1). Almost all of the energy loss is homogeneously deposited in ionizations by the primary ions and recoils across the target, with a negligible contribution of phonons by the recoils (~0.1%) and almost no atom displacement (Table 3). Calculations with the CASP 5.0 computer code [30] yield values for the electronic stopping power of S_e = 1.721 × 10⁻¹² eV cm² and S_e = 2.163 × 10⁻¹² eV cm² per molecule for the ²⁸Ni and ³⁶Kr ions in LiF, respectively, corresponding to ionizations of the 1s, 2s, and 2p shells. This gives S_e = 10.5 MeV µm⁻¹ and S_e = 13.2 MeV µm⁻¹, respectively, for A = 26 g and ρ₀ = 2.64 g cm⁻³, which are in rather good agreement with the values of S_e = 9.34 MeV µm⁻¹ and S_e = 11.7 MeV µm⁻¹ computed with the SRIM2013 code (Table I). The other values computed with CASP 5.0 for SiO₂ (for A = 60 g and ρ₀ = 2.65 g cm⁻³), corresponding to the ionizations of 1s, 2s, 2p, 3s, and 3p shells, are given in comparison (Table I). Similar deviations by about 10% between both computer codes were already found for a lower velocity ion (82-keV u⁻¹ Si ions) in SiC [30] for a similar mass density (ρ₀ = 3.21 g cm⁻³).

4. Discussion

The present combined simulations of PHITS and SRIM codes show that the energetic ionized atoms of the targets generate a great number of phonons in a very small volume due to their low kinetic energies below the threshold displacement energies. The primary Coulomb spike generates a secondary elastic collision thermal spike. This is similar to the elastic collision spike used for modeling the ion sputtering by Sigmund [31]. This process must be distinguished from the i-TS model, in which the energy is mainly imparted to the electron subsystem through ionization and electronic excitations and further transferred to the atomic subsystem by the so-called electron–phonon coupling constant [4], as explained above.

The energy deposited by the primary swift Ni and Kr ions leads to the formation of a number of electron–hole pairs ${\mathrm{N}}_{\mathrm{e}-\mathrm{h}}$ ~ ${\displaystyle \frac{\mathsf{\delta }\mathrm{E}}{{3\mathrm{E}}_{\mathrm{G}}}}$ [32], which is about 600 per ion for the band-gap energy of E_G ~10 eV. Few point defects are likely generated after the STE formation and subsequent relaxation of the atomic lattice [33]. Instead, ionization tracks are induced by these swift ions for electronic stopping powers (Table 1) that are higher than the thresholds of ~1 MeV µm⁻¹ and 4 MeV µm⁻¹ in SiO₂ and LiF, respectively. However, the ionization processes of atoms and further steps of ionized atom motion at a very short time scale (<1 ps) were disregarded in the analysis of track formation by the i-TS model.

Here, we consider the becoming of these ionized atoms with a low electronic energy loss compared to the phonon energy loss (Table 3). We calculate the energy density deposited in heat in the ion range compared to the melting enthalpy. The energy ∆E is deposited in a volume of $\mathrm{V}={\displaystyle \frac{4}{3}}\mathsf{\pi }{\mathrm{R}}_{\mathrm{p}}^3$ ~32 nm³ for R_p ~2 nm. Considering that about 80% of the total energy loss $\left(-{\displaystyle \frac{\mathrm{d}\mathrm{E}}{\mathrm{d}\mathrm{x}}}\right)={\mathrm{S}}_{\mathrm{e}}+{\mathrm{S}}_{\mathrm{n}}$ of ions is deposited in phonons according to Equation (1):

$$∆\mathrm{E}~0.8(-{\displaystyle \frac{\mathrm{d}\mathrm{E}}{\mathrm{d}\mathrm{x}}}){\mathrm{R}}_{\mathrm{p}}$$

(1)

Hence the deposited energy density in a volume of ${\displaystyle \frac{4\mathsf{\pi }{\mathrm{R}}_{\mathrm{p}}^3}{3}}$ is given in Equation (2):

$$\mathrm{D}={\displaystyle \frac{∆\mathrm{E}}{\mathrm{V}}}~0.8{\displaystyle \frac{3}{4\mathsf{\pi }{\mathrm{R}}_{\mathrm{p}}^3}}(-{\displaystyle \frac{\mathrm{d}\mathrm{E}}{\mathrm{d}\mathrm{x}}})$$

(2)

This gives ∆E~489 eV and D~14.6 eV nm⁻³ for R_p~2 nm by using $(-{\displaystyle \frac{\mathrm{d}\mathrm{E}}{\mathrm{d}\mathrm{x}}})$~306 keV µm⁻¹ = 306 eV nm⁻¹ for the Si ions in SiO₂ (Table 1). The energy per atom (E_a) is calculated from:

$${\mathrm{E}}_{\mathrm{a}}={\displaystyle \frac{∆\mathrm{E}\mathrm{v}}{\mathrm{n}\mathrm{V}}}$$

(3)

where v is the unit cell volume, and n is the number of atoms per unit cell. For hexagonal α-SiO₂ quartz, n = 9 (three Si and six O atoms), with the lattice parameters of a = 0.49 nm and c = 0.54 nm. This gives $\mathrm{v}=3{\mathrm{a}}^2\mathrm{c}{\displaystyle \frac{\sqrt{3}}{2}}$ ~0.34 nm³ and E_a ~0.6 eV at⁻¹. However, the contribution of O ions must be also added with ∆E ~318 eV and E_a ~0.4 eV at⁻¹ by using $(-{\displaystyle \frac{\mathrm{d}\mathrm{E}}{\mathrm{d}\mathrm{x}}})$ ~199 keV µm⁻¹ (Table 1). This yields a total deposited energy per atom of ~1 eV for both Si and O ions.

For the F ions in LiF, this gives ∆E~427 eV and D~13.3 eV nm⁻³ for R_p~2 nm, by using $(-{\displaystyle \frac{\mathrm{d}\mathrm{E}}{\mathrm{d}\mathrm{x}}})$ ~267 eV nm⁻¹ (Table 1). For LiF with the fcc NaCl rock-salt structure, n = 8 (4 Li and 4 F atoms) with the lattice parameter of a = 0.40 nm and v = a³ ~0.06 nm³, E_a~0.1 eV at⁻¹. The contribution of the Li ions is much smaller, for R_p ~4 nm and $(-{\displaystyle \frac{\mathrm{d}\mathrm{E}}{\mathrm{d}\mathrm{x}}})$ ~88 eV nm⁻¹, increasing E_a by ~10% only.

The energy density deposited in the phonon system ultimately yields a high temperature at the time scale (${\mathsf{\tau }}_{\mathrm{D}}={\mathsf{\nu }}_{\mathrm{D}}^{-1}$ ~1 ps) corresponding to the Debye frequency (${\mathsf{\nu }}_{\mathrm{D}}$ ~10¹² s⁻¹) (Figure 2). The local lattice ion temperature (T) (in Kelvins) may be estimated and compared to the melting temperature (T_m) from Equation (4):

$${\int }_{300}^{\mathrm{T}}\mathsf{\rho }\mathrm{V}{\mathrm{C}}_{\mathrm{p}}\left(\mathrm{T}\right)\mathrm{d}\mathrm{T}=∆\mathrm{E}$$

(4)

where m = ρV is the mass of matter corresponding to the heated volume V, and C_p is the specific heat at constant pressure which depends on temperature with a polynomial function (Shomate’s equation) in Equation (5):

$${\mathrm{C}}_{\mathrm{p}}=\mathrm{A}+\mathrm{B}\mathrm{T}+\mathrm{C}{\mathrm{T}}^2+\mathrm{D}{\mathrm{T}}^3+\mathrm{E}{\mathrm{T}}^{-2}$$

(5)

where A, B, C, D, and E are coefficients [34]. The integration of Equation (4) leads to a fourth-degree polynomial equation to solve for T by knowing the five coefficients. The same equation as Equation (4) can be used for the melting temperature by replacing T with T_m and ∆E with ∆H_m. A similar equation as Equation (5) holds for the standard enthalpy ∆H with respect to the value at 298 K after the integration of Equation (5):

$$∆\mathrm{H}=\mathrm{A}\mathrm{t}+{\displaystyle \frac{\mathrm{B}}{2}}{\mathrm{T}}^2+{\displaystyle \frac{\mathrm{C}}{3}}{\mathrm{T}}^3+{\displaystyle \frac{\mathrm{D}}{4}}{\mathrm{T}}^4-\mathrm{E}{\mathrm{T}}^{-1}+\mathrm{F}-\mathrm{G}$$

(6)

For solid LiF, A = 41.75837, B = 18.71110, C = 0.693674, D = −0.992621, E = −0.487055, F = −631.8342, and G = −616.9308, where T is expressed in mK, for 298 K ≤ T ≤ 1121 K [34]. For liquid LiF, A = 64.18298, B = 4.195423 × 10⁻¹¹, C = −2.302271 × 10⁻¹¹, D = 3.924242 × 10⁻¹², E = 1.150237 × 10⁻¹², F = −617.7885, and G = −598.6509 for 1121 K ≤ T ≤ 3000 K [34]. The values of C_p and ∆H are plotted versus T, with positive jumps corresponding to the melting point at T_m = 1117.8 K under normal pressure (Figure 8). This yields calculated values of C_p = 42 J mol⁻¹ K⁻¹ and ∆H = 0.09 kJ mol⁻¹ for T = 300 K and C_p = 61.81 J mol⁻¹ K⁻¹ and ∆H = 44.03 kJ mol⁻¹ for T = 1121 K (Figure 8). This is consistent with the scarce literature data, giving a value of C_p = 43.4 J mol⁻¹ K⁻¹ for 400 K [35]. For ∆H = E_a ~0.1 eV at⁻¹, this gives ~0.2 eV per LiF molecule ~20 kJ mol⁻¹, one finds that T ~700 K (Figure 8).

For solid SiO₂, similar calculations can be made with the two sets of parameters: A = −6.076591, B = 251.6755, C = −324.7964, D = 168.5604, E = 0.002548, F = −917.6893, and G = −27.96962 for 298 K ≤ T ≤ 847 K and A = 58.75340, B = 10.27925, C = −0.131384, D = 0.025210, E = 0.025601, F = −929.3292, and G = −910.8568 for 847 K ≤ T ≤ 1996 K (Figure 8) [34]. There is the α → β phase change at T = 847 K, with a negative jump of C_p, but a weak change in ∆H. The values of C_p and ∆H are plotted versus T (Figure 8). These plots yield C_p = 44.8 J mol⁻¹ K⁻¹ and ∆H = 0.08 kJ mol⁻¹ for T = 300 K and C_p = 78.95 J mol⁻¹ K⁻¹ and ∆H = 119 kJ mol⁻¹ for T = 1996 K (Figure 8). These calculated values are in good agreement with the drop calorimetry data for quartz yielding C_p = 44.7 J mol⁻¹ K⁻¹ and ∆H = 0.083 kJ mol⁻¹ for T = 300 K and C_p = 75.873 J mol⁻¹ K⁻¹ and ∆H = 118.226 kJ mol⁻¹ for T = 2000 K [36]. For ∆H = E_a ~1 eV at⁻¹. This gives ~3 eV per SiO₂ molecule ~300 kJ mol⁻¹; one finds that T >> T_m (Figure 8) within the approximations of our calculations. When crossing the melting point for T = T_m, a part of this energy (E_a) is spent in the enthalpy of fusion of about 10 kJ mol⁻¹ [36].

After the energy deposition in heat at the 1 ps time scale, cooling of the embedded liquid pocket occurs at a later stage by thermal conduction in the neighboring solid host. Even though the lattice thermal conductivity of SiO₂ quartz is low (κ = 1.35 W m⁻¹ K⁻¹ at T = 300 K and κ = 2.52 W m⁻¹ K⁻¹ at T = 1100 K [37]), fast cooling is expected, since the molten volume (V ~32 nm³) is very small. The formation of amorphous domains is likely after fast thermal quenching of the liquid pockets. For example, fast cooling of liquid metals for high quenching rates yields amorphous metallic glasses [38].

The lattice thermal diffusivity (D) is the important parameter involved in spreading out the heat flow:

$$\mathrm{D}={\displaystyle \frac{\mathsf{\kappa }}{\mathsf{\rho }{\mathrm{C}}_{\mathrm{p}}}}$$

(7)

where κ is the thermal conductivity coefficient. In the one-dimensional approximation, the time-dependent heat diffusion equation is written:

$${\displaystyle \frac{\partial \mathrm{T}}{\partial \mathrm{t}}}=\mathrm{D}{\displaystyle \frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{x}}^2}}$$

(8)

For solid α-SiO₂ quartz at T = 1100 K, ρ = ρ₀ (1 − α T) ~ ρ₀ = 2.65 g cm⁻³, since α = 0.75 × 10⁻⁶ K⁻¹ [39]. One finds a value of D = 7 × 10⁻³ cm² s⁻¹, for C_p = 69 J mol⁻¹ K⁻¹ = 1.35 J g⁻¹ K⁻¹ (Figure 8), and κ = 0.025 W cm⁻¹ K⁻¹ at T = 1100 K [37]. The mean diffusion time (τ_D), taking into account the characteristic diffusion length of R_p = 2 nm, is given by:

$${\mathsf{\tau }}_{\mathrm{D}}={\displaystyle \frac{{\mathrm{R}}_{\mathrm{p}}^2}{\mathrm{D}}}~5.7\times {10}^{-12} \mathrm{s}$$

(9)

This rough estimation based on the available experimental data at 1100 K means that heat is quickly flowing out of the molten zones at T_m = 2000 K to the matrix at T₀ = 300 K. This yields a very high cooling rate of ${\displaystyle \frac{∆\mathrm{T}}{{\mathsf{\tau }}_{\mathrm{D}}}}={\displaystyle \frac{{\mathrm{T}}_{\mathrm{m}}-{\mathrm{T}}_0}{{\mathsf{\tau }}_{\mathrm{D}}}}={\displaystyle \frac{1700}{5.7\times {10}^{-12}}} {\mathrm{K}\mathrm{s}}^{-1}$ ~3.0 × 10¹⁴ K s⁻¹. For such a very high quenching rate, the liquid will turn to an amorphous silica phase. Complete amorphization in the ion track can be reached after the overlap of these small amorphous domains.

For solid LiF, ρ = ρ₀ (1 − α T) ~ ρ₀ = 2.64 g cm⁻³, since α = 37 × 10⁻⁶ K⁻¹ from 273 K up to 373 K [40], and κ ~14 W m⁻¹ K⁻¹ at T = 300 K and κ ~5 W m⁻¹ K⁻¹ at T = 1000 K [41]. This yields D = 8.1 × 10⁻³ cm² s⁻¹, for C_p = 59.7 J mol⁻¹ K⁻¹ = 2.31 J g⁻¹ K⁻¹ (Figure 8) and κ ~0.05 W cm⁻¹ K⁻¹ for T = 1000 K. Assuming a value of τ_D = 4.87 × 10⁻¹² s deduced from Equation (9), for T = 1000 K, the cooling rate from the maximum reached temperature of T = 700 K is of ${\displaystyle \frac{∆\mathrm{T}}{{\mathsf{\tau }}_{\mathrm{D}}}}~{\displaystyle \frac{\mathrm{T}-{\mathrm{T}}_0}{ {\mathsf{\tau }}_{\mathrm{D}}}}~{\displaystyle \frac{400}{4.87\times {10}^{-12}}} {\mathrm{K}\mathrm{s}}^{-1}$ ~8.2 × 10¹³ K s⁻¹. For such a high quenching rate, the solid will retain at 300 K the thermal vacancies formed at high temperature. The formation energy of the Li vacancy (V_Li) is E_f = 0.74 ± 0.04 eV at T = 680 K [42]. Recent DFT calculations yield consistent values of E_f ~1 eV for V_Li and E_f ~4 eV for V_F for an expected Fermi energy of E_F ~6 eV [43]. This means that Li vacancies will be preferentially generated by this thermal process after quenching. A rough estimation of the V_Li concentration for E_f = 0.7 eV is $\left[{\mathrm{V}}_{\mathrm{L}\mathrm{i}}\right]={\mathrm{e}}^{-{\displaystyle \frac{{\mathrm{E}}_{\mathrm{f}}}{{\mathrm{k}}_{\mathrm{B}}\mathrm{T}}}}$ ~6.2 × 10⁻⁴ for T = T_m and $\left[{\mathrm{V}}_{\mathrm{L}\mathrm{i}}\right]$ ~9.2 × 10⁻⁶ for T = 700 K. This can contribute to the observed swelling [18] deriving from point-defect build up, even though LiF is not amorphized. Although this is not the classical interpretation of radiation damage by electronic excitations in wide band-gap insulators, the present analysis shows that the contribution to the radiation damage of the hot ionized atoms cannot be overlooked. A graphical model of the Coulomb spike model is displayed in Figure 9.

5. Conclusions

The radiation damage in LiF and SiO₂ by 1.0 MeV u^{−1 28}Ni and ³⁶Kr ions for about the same electronic stopping power of ~10 keV nm⁻¹ was modeled by using a novel Coulomb spike concept. Our methodology was developed to help understand the radiation damage induced by fission fragments in nuclear fuels. The ionization of constituent atoms of the target was calculated on an event-by-event basis using the ITSART computer code for the track-structure calculation model of radiation transport code PHITS. Based on the coordinates of the calculated ionizations, the electrostatic energy between ions was calculated offline. The atomic collisions’ cascades induced by these ions of about 500 eV in kinetic energy were then simulated with the SRIM2013 code. It is found that the major part of energy (80%) is deposited in phonons in a small volume due to the subthreshold events that are unable to induce atomic displacements by nuclear collisions. As a result, the elastic thermal spike induced by the Si and O ions leads to the melting of SiO₂ in a small volume, whereas melting is not reached in LiF by the lighter Li and F ions. The subsequent fast cooling of the molten pockets likely leads to small amorphous silica zones embedded in crystalline SiO₂ quartz.