Coulomb explosion and steadiness of high-radioactive silicate glasses

The paper is devoted to the theoretical study of Coulomb explosion in silicate glasses with low ionization potential under internal alpha-irradiation. The phenomenon was studied in the way of computer simulation, namely by the molecular dynamics (MD) method; parameters of the so-called lavalike fuel-containing materials (LFCM), were chosen as input parameters for the model due to its practical importance. LFCM are high-radioactive glasses, which were formed during an active stage of a well known heavy nuclear accident, occurred on Chornobyl NPP in 1986. Computer simulation revealed that Coulomb explosion really may occur in the LFCMs and leads to additional radiation damages under internal alpha-irradiation. The total quantity of atomic displacements produced in the way of Coulomb explosion from each alpha-particle track is 40000 to 80000, which exceeds radiation damages from alpha-particle and heavy recoil nuclei altogether (about 3500) more than one order.


Introduction
Coulomb explosion is the phenomenon of rapid spreading of positively charged ions in condensed matter as a result of coulomb repulsion when high-density ionization of ensemble of neighbouring atoms takes place.This process has a threshold-like dependence on ionization density, and when the energy of coulomb repulsion for ensemble of ions surpasses the chemical bound energy, the process will develop in explosion-like way, that is why this phenomenon is named a Coulomb explosion.
For the first time the phenomenon was detected in ionic crystals and solid rare gases, where electrostatic interaction dominates over all other types of chemical bonds [1].
Coulomb explosion generally requires the formation of high-density cluster of positively charged ions much faster than charge relaxation time for current media.

V.Zhydkov
Typical size of such a cluster is about 1 nm and typical time for this process is 200 fs.Due to a relatively long duration necessary to transform a potential energy into a kinetic energy the phenomenon occurs almost exclusively in dielectrics and semiconductors, where ionic cluster cannot dissipate fast enough due to low conductivity.
The formation of ion cluster sufficient for the Coulomb explosion can be provided by intensive and short-time ionization energy pulses, which can be delivered to the matter by means of a strong laser pulse [2] or high-energy charged particle [3].
The Coulomb explosion effects many phenomena; the most significant of them are the surface sputtering and radiation damages.The Coulomb explosion in solids can lead to a surface sputtering, defects, and clusters of defect formation.The Coulomb explosion occurring as a result of charged particle flight leads to giant sputtering yield value [4,5] and to a large amount of radiation damages [2,6], mainly atomic displacements, which vastly increases DPA (Displacements Per Atom) value in comparison to traditional collision model.
Each Coulomb explosion act produces a significant number of displaced atoms [2].Those atomic displacements are rather close-placed and form a large cluster of radiation defects.This type of defects generally does not provide a visible effect on macroscopic parameters of materials to a great extent, but rather serves as cracked nuclei, which can significantly reduce their macroscopic mechanical steadiness if a great number of such microcracks are formed.Such a cluster may dissolve in the material volume with time, but high density of defects usually prevents this.Such a process significantly increases a DPA quantity and leads to degradation of materials undergoing self-irradiation due to internal agents.
An energy necessary for Coulomb explosion will be partially spent on ionization, so the increase of ionization potential leads to a further increase of energy density threshold.The high ionization energy is the main reason why Coulomb explosion normally does not occur in pure quartz glass, except for the case of heavy multicharged ions [3], which can provide a large amount of potential energy.In glasses having low ionization energy for electrons, however, the phenomenon may be significant.Practical importance of glassy compositions, where an electron transport dominates, stimulates the proposed design-theoretical study of Coulomb explosion process for this type of materials.
First of all, one may notice, that Coulomb explosion effect is very hard to calculate in analytical way.It needs taking into consideration a lot of factors, the majority of which are significantly non-linear (potential, charge distribution).The most convenient statistical and thermodynamic methods are generally unacceptable here because the system is far from equilibrium state (where the methods are valid) at the moment of explosion.All the analytical solutions for such a phenomenon can be made only with a number of assumptions, which do make the final result just a rough estimation.
Experimental study encounters even the worse hardships.It is impossible to observe a Coulomb explosion in "real time", because typical duration for such a phenomenon is less than 1 ps, which is a bit of challenge even for the best of time-resolve techniques; but even those methods cannot be used, since all the processes always take place somewhere beneath the material surface.One may experimentally observe an indirect evidence of Coulomb explosion only -the defects produced in the depth of material and sputtered particles emitted by the surface in both.The most successful way of observing is step-by-step slicing of investigated material with further etching and electron microscopy, which is non-productive.
Coulomb explosion study in computer simulations gives at least better results than "pure" analytical approach, because it permits to take into consideration a non-linearity and other things which are hard to calculate analytically.A conventional method for Coulomb explosion study is the Molecular Dynamics (MD) method, which is based on classical approach to atom movements in given potentials.This is quite correct approach if the potentials are established correctly.It will provide almost accurate results, if inter-atomic interactions are calculated using their wavefunctions (Car-Parinello approach).A time-dependent simulation of Molecular Dynamics permits to study the physical and chemical properties changeability with time, which makes it convenient to understand.This method is non-contradictive and the results can be verified experimentally.
The MD method was chosen as the basic for evaluations in the current work, as the most promising and giving the most detailed representation of processes taken place.

High-radioactive glasses as the object for simulation
Until now all the known manifestations of Coulomb explosion have been observed mainly on the surface of various materials under the action of external short highenergy pulses and particle bombardment.Meanwhile, there are materials containing radioactive atoms dissolved in their volume, where the Coulomb explosion may occur if their properties (charge relaxation time, ionization energy) are suitable.
LFCM (Lava-like Fuel-Containing Materials) are the materials formed as a result of well-known heavy nuclear accident occurred at Chornobyl NPP facility in 1986.As it is accepted now [7], the LFCM are glass ceramic alkaline-earth silicate compositions (devitrified glasses) containing 5-10 mass percent of irradiated uranium nuclear fuel in their volume accompanied by high-radioactive fission products.Quite up-to-date review of the main LFCM physical properties is represented in [7].Approximately 1,000 tons of this material, accompanied by both the other core debris and the building constructions of the destroyed Chornobyl NPP 4-th unit, is located in the so-called "Shelter" facility, which had been erected quickly after 1986 accident as a forced measure directed on primary prevention of radioactive substance dissemination in the environment.The "Shelter" facility is not a hermetically sealed construction and there are no doubts that a large quantity of high-radioactive compositions, such as LFCM and irradiated nuclear fuel itself do have a direct air contact with the environment.LFCM themselves look like the coloured glass ceramics and one can classify them like brown (8-10% fuel content) and black (4% fuel as usually).
Recent experimental research indicates, that LFCM are disordered semiconductors having energy gap width about 1.5 eV and Fermi level lieing 0.75 eV below the mobility edge of conductivity band.These results were obtained mainly from the measurements of LFCM electric conductivity on temperature dependence [8].An ionic conductivity in LFCM turned out to be suppressed by Ca 2+ and Mg 2+ ion presence and the electric transport in them provided mainly by electrons.
A LFCM specific radioactivity corresponds to the fuel content, but there are slight deviations in isotopic proportions: due to some distinguishing features of lava formation [9], LFCM as a whole turned out both to be depleted twice or more by volatile fission products, such as Cs-137, and correspondingly enriched with some transuranium daughter products (Pu+Am) in comparison with an irradiated fuel of a similar burnup.Thus, up-to-date LFCM specific activity is about 20 MBk/cm −3 (α-activity) and can achieve 1 GBk/cm −3 (β-emitters, mainly Sr-90).Total estimated absorbed radiation dose in LFCM volume achieves of the order of 10 MGy up to the moment, which is a significant level for dielectrics, where the energy dissipation mechanisms through free electrons are negligible.
Another important property in LFCM physical behaviour is self-sputtering of their surface.In the recent years a very interesting experimental fact was discovered [10] -both the LFCM and irradiated uranium fuel surface have one remarkable property: spontaneous dust generation capability, which means the capability of a surface, without any external impact, to generate and disseminate into the surrounding media the high-disperse phase (with the size of the particles mainly below 1µm), which in technical practice should be named as a submicronic dust.Such a phenomenon had been preliminary investigated in workbench experiments, both by providing smear tests [10] and by collecting the dust on special collectors in a high vacuum condition.Recently the spontaneous dust emission phenomenon has been independently detected in workbench experiments in Karlsruhe [11] where it was identified as the emission of submicronic dust particles from recrystallized UO 2 fuel of a high burnup.
A certain possible physical mechanisms to be responsible for the above phenomenon are unclear yet.No doubt, however, that particle emission from the surface means the displacement of surface atoms (clusters) and should originate from the processes quite similar to those leading to radiation damages.In [22] the possible physical mechanisms for the observed phenomenon were under discussion, but the giant sputtering yield, detected in experiment, was left out the satisfactory explanation.

The way for molecular dynamics simulation
The first step for any simulation is elaboration of input parameters and establishing the physical model.As it was noticed earlier, one of the most important things for simulation by the molecular dynamics method is a potential choice.The current choice was taken between L-J (Lennard-Jones) potential, and material-specific potentials, namely BKS (van Beest-Kramer-van Santen), TTAM and three-body VKRE potential.
L-J was the first potential chosen for the research [12], but later it turned out that the potential isn't good enough for simulation of condensed matter containing various sorts of atoms.So, the L-J potential was rejected and it was decided to choose more precise potentials.Silica is the base of LFCM, as it was noticed earlier, so it was decided to apply silica-specific potentials: BKS [13], Tsuneyuki, Tsukada, Aoki and Matui (TTAM) [14] and Vashista, Kalia, Rino, Ebbsjo(VKRE) [15].BKS and TTAM do have common view in exp-6 form: where r ij is the distance between interacting atoms, q i and q j are their effective charges, and i,j -their numbers.HereA ij , B ij , C ij are the specified potential constants, which depend on the sort of interacting atoms.The constants for BKS and TTAM [13,14] are listed in Appendix.Effective charges are supposed to be −1.2e for O atoms and 2.4e for Si atoms as suggested by authors.
As one may notice, the BKS and TTAM potentials do have an unphysical property of diverging to minus infinity at small inter-atomic distances.At high temperatures or some other conditions, it is a possible situation in the model simulated when an individual atom may achieve sufficient kinetic energy to overcome the potential barrier and fuse with another atom.For this reason, the potentials were replaced with a simple harmonic potential when atoms become closed within a range 0.12 nm(for Si-O interaction) and 0.17 nm (for O-O interaction).The Si-Si interaction is not considered here because for this case both potentials have coefficient at the diverging component (1/r 6 ) equal to zero.
The TTAM potential was developed for the accurate evaluation of properties in a crystalline state [14].Derived from ab-initio Hartree-Fock self-consistent-field calculations, Tsuneyuki et al. have demonstrated that four of the known crystalline polymorphs are stable under this potential.They have also used this potential to show the α to β structural phase transition at 850 K [16].Hemmati and Angell have shown that the BKS potential is superior to the TTAM potential for certain properties of amorphous silica [17].However, we include the TTAM potential in our study for comparison purposes.
The BKS potential is based on ab initio calculations and experimental data [13].The authors claimed this as an improvement over the TTAM potential and compared several structural observations between them.The BKS potential has the same functional form as the TTAM potential.The only important difference between them is that for BKS potential the Si-Si interaction has no short-range component.
The VKRE potential contains an additional three-body component in order to simulate bond stretching and bending effects.
where ϕ (2) ij is two-body potential part, which includes a Coulomb interaction, repulsion connected with finite dimensions of ions, and a charge-dipole interaction originating from the large electronic polarizability of O atoms; and ϕ (3) ijk is a threebody part that has the form where i,j,k are the atomic numbers.The B ijk constants and formula for f (r ij , r jk ) function and average bond angle θ ijk values are given in Appendix as well.By this part of potential the "bend" strength is taken into account, based on the estimated strength of bond interaction and average bond angle θ.A graphical comparison of two-body parts of potentials is presented in figure 1.In numerical experiments simulating a collision of amorphous silica structure with separate Si and O atoms having a kinetic energy of 10, 20, and 50 eV, the BKS potential has been identified as more stable and providing smaller energy divergence, compared with TTAM and VKRE potential, and therefore was chosen as the most appropriate for simulations, which include collisions.Further tests and [18] in both allow a suggestion, that BKS potential is more suitable for simulation of amorphous silica.Regarding the above properties, BKS potential seems the most appropriate one for simulations of defect formation in silica within the current research.
Coulomb part of potentials was calculated taking into account a Debye screening.The Debye-Huckel screened Coulomb potential was estimated using the expression ϕ(r) = q i q j r exp(−k 0 r).
Here k 0 is the decay parameter, which was calculated from the known relation where 0 is a dielectric constant and k B is the Boltzmann constant.As far as the free electron concentration n is less than10 11 cm −3 at a room temperature, as it was estimated for LFCM [8], value is too low (below 3•10 5 m −1 ) to provide effective screening at scales, where Coulomb explosion takes place (nanometers), but it was taken into account for the sake of precision.For the cluster of positively charged atoms, which is formed at alpha-particle flight, the screening effect was excluded, because effective Debye screening is impossible in the area where there is a lack of electrons.

Calculation for the charge density threshold
As it will be shown below, the energy linear density to be sufficient for Coulomb explosion for alpha-particle is about 300 eV/nm.Such an energy loss rate corresponds to the energy of alpha-particle of 0.5-1.2MeV [19].The energy losses for energies higher than 1.2 MeV is lower because of weak interaction of particle with media, and for energies less than 0.5 MeV the energy losses are also lower due to low alpha-particle energy.Linear energy losses function for alpha-particle was determined from the reference table for free path of alpha-particles in various media [19].Then the ionization threshold for LFCM, regarding the existing data [8], was suggested as 1 eV without any exaggeration.Cascade of secondary electrons from the particle collisions was suggested to be insufficient, because at ionization energy of 1 eV the major part of alpha-particle energy dissipates into a direct ionization process.
Finally, the charge distribution was calculated by simulation of a direct ionization process and checked for correctness under assumption, that total electrostatic potential energy increase should be equal to the energy losses.The last suggestion seems to be realistic, due to the system being adiabatic.The calculated threshold for local density of ionized atoms is about 50 nm −3 at the track center, which corresponds to the linear energy losses of 300 eV/nm.The estimated ionization track width is about 1.8 nm.

The simulation procedure
The next important point after setting the type of potential and initial conditions is choosing the correct equations of motion in order to calculate the atom positions from forces with time.Theoretically, the most correct one for general purposes of molecular dynamics is Gear's predictor-corrector algorithm [20].It turned out, however, that the corrections used in this algorithm produce additional errors when simulating collisions and other rapid changes in the system.Further, the most correct results were obtained when using the lowest order of the Gear's predictorcorrector algorithm.That is why the decision was to reject the Gear's algorithm, and use the Verlet integration algorithm with smaller timestep instead of it.
Simulation was set for a lattice of 9000 atoms with translational periodical conditions for 2 nm in "depth" along the alpha-particle track and 8 nm in "width" (perpendicular to the track).The system under simulation underwent an instant heating to 7000 K, and then was quickly cooled in order to form an amorphous silica structure similar to LFCM.Some of the calculated parameters for the amorphous structure are demonstrated in figure 2. After stabilization of the whole structure the simulation of alpha-particle flight has been performed.Then simulation was run by the established rules of molecular dynamics with a timestep of 0.05 fs.During the simulation process the system was being constantly checked for the energy conservation criteria and if the divergence turned out to be too large (> 0.01%), the timestep was decreased.The total energy conservation was also the crucial criteria to check the correctness of simulation process.Total energy deviation starting from 0 fs and until 200 fs, when all the atomic displacement had taken place, did not exceed the 10 eV quantity for the whole system, which is a satisfactory accuracy in comparison with alpha-particle linear energy loss (300 eV/nm).
The main results of simulation are represented in figures 3, 4.
In figure 3 one can see that majority of atomic displacements fall below the 0.2 nm range.Such small displacements cannot really be considered as "true" ones for this case as far as they cannot provide the stable point defects due to amorphous material structure and therefore cannot be classified as defects influencing on the material properties.So, for the correct DPA calculations one should take into account the only atoms being displaced on a distance more than size of crystalline silica cell (near 0.38 nm).Those atoms only have been accounted as real displacements.The Coulomb explosion has practically 100% chance to occur, when the particle linear (ionization) energy losses will exceed 300 eV/nm values.If energy losses will be 100÷300 eV/nm, the Coulomb explosion still has significant chances to occur, but its total power will be of one order smaller.
If to calculate DPA in the above way, it will result in about 40000÷80000 displacements per each alpha-particle track depending on variations in material structure.This quantity is much more, than radiation damage quantity produced by direct collisions of alpha-particle with atoms (300) and the ones produced by heavy recoil nuclei (3000) [21].The exact parameters for the radiation steadiness of LFCM are unknown yet, but Coulomb explosion accounting in estimated radiation damages quantity will no doubt have the significant influence on LFCM behavior prognosis with time.The total estimated radiation damage quantities in LFCM regarding the Coulomb explosion phenomenon are listed in the table 1.Additionally, each alpha-particle flight will be accompanied by 100000÷200000 broken silicate cycles (where Si bound with 3 oxygen atoms only) at 100 fs from the moment of explosion, which may turn LFCM into molecular-pinhole porous system.That is the additional factor for LFCM structure degradation, allowing oxygen from air and water to penetrate into LFCM volume.
Naturally, the Coulomb explosion phenomenon does not exclude the existence of "thermal spike" phenomena [23], because all the kinetic energy from Coulomb explosion will finally be converted into a thermal (macroscopic) energy and can form the "thermal spike" following just after Coulomb explosion and providing an additional input in radiation damages.
Figure 4 illustrates the spatial distribution of kinetic and potential energy with time.We can see multiple atoms in area of alpha-particle track gain high potential energy, which is rapidly transformed into kinetic energy (extremely high effective temperature), which rapidly diffuse and dissipate their energy in surrounding media, transforming it into the "thermal spike".

The practical outcomes
The results obtained here are the basis for further radiation damages quantity calculations in silicate glasses having low ionization potential.These results may provide the satisfactory explanation for high-radioactive dielectrics self-sputtering [10][11][12] phenomenon.In [22] the possible mechanisms of self-sputtering were under discussion.General explanation for observable giant sputtering yield was based on considerations about electron sputtering [22], where numerous electronic excitations may stimulate atomic displacements.It was not clear, however, in what namely way the low-energy electronic excitations can lead to atomic displacements, which needs much higher energy being beyond a certain threshold (16 eV for crystalline SiO 2 ).The Coulomb explosion phenomenon seems to be the namely physical mechanism for such an energy transfer.It is not difficult to predict that Coulomb explosion intensity will be enlarged greatly, if we apply an external electric field to the whole system; regarding a high electric field occurring in the near-surface layer of radioactive dielectrics, one can suggest that the near-surface Coulomb explosion should be the main mechanism responsible for self-sputtering phenomena.That case, however, is more difficult for quantitative modelling as well as to be the subject for further research.
Another very important result is the established correlation (for dielectrics) between the electron ionization energy and their potential radiation steadiness.The latter is of doubtless importance in creating a scientifically based forecast of LFCM behaviour.
As far as LFCM belongs to a class of silicate glasses containing dissolved radioactive substances, the result may be universal for similar silicate matrices initially appointed for α-radioactive materials immobilization in technologies where the vitrification of solid radioactive waste is under suggestion.

Figure 1 .
Figure 1.The potential energy identified for BKS, TTAM and VKRE potentials.Only two-body part of VKRE potential is shown.Coulomb interaction is omitted.

Figure 2 .
Figure 2. Radial distribution of distances between atoms in the material simulated.Results for O and Si atoms for the amorphous silica are shown.

Figure 3 .
Figure 3. Atomic displacements on range distribution as the result of Coulomb explosion.

Figure 4 .
Figure 4. Distribution of kinetic and potential energy in material in the chosen moments of time after alpha-particle flight.Left picture is the kinetic energy, right is the potential energy and bottom is time passed after alpha-particle flight.The alpha-particle flight direction is perpendicular to a picture surface.Energy value corresponds to levels of gray, as shown on the scale by the left.

Table 1 .
Estimated accumulation of radiation damages in LFCM with time.