Abstract

High-temperature superconducting material is a promising candidate to fabricate superconducting magnet for magnetic confinement fusion reactors. The DPA number of the 1 µm thick superconducting layer in a high temperature superconducting tape under neutron irradiation needs to be calculated to predict the property changes. The DPA cross sections, which ignore the spatial distribution of vacancies caused by PKAs, are commonly used to obtain the results of the damage energy and DPA. However, for geometric models with the thickness as small as 1 µm, the energy and angular distribution of PKAs reveal that a significant number of PKAs with relatively high energy tend to scatter forward and cross the boundary of model, so the thickness of model has the potential to affect the number of displaced atoms. In this paper, we developed a method based on Geant4 and SRIM to evaluate the deviation of the traditional analytic method caused by the thickness. Geant4 is used to obtain the location, direction, and energy of PKAs, while SRIM is used to track every PKA and obtain damage energy and the number of displaced atoms. The radiation damage calculation of simple thin plate models with different thicknesses and the tape model are conducted with the neutron energies from 1 to 14 MeV. The results show that PKAs need to be tracked continuously for models with thickness less than 10 µm and the deviation of the analytic formulas increases rapidly with the decrease of thickness. For the superconducting layer composed of four different elements in the tape, the deviation also depends on the proportion of each atomic species and the neutron-atom interaction cross sections under different incident neutron energy.

1. Introduction

Interactions of incident particles with materials may induce heavy recoils to move away from their original lattice positions resulting from the kinetic energy transferred towards the target atoms. The recoils, directly generated in the collision between neutrons and target atoms, are known as primary knock-on atoms (PKAs). PKAs with the sufficient initial energy above the threshold displacement energy (Ed) have the potential to cause the displacement of lattice atoms and create considerable amounts of vacancies and interstitials, also called Frenkel Pairs (FPs). The accumulation of the FPs and other consequences of irradiation damage can significantly degrade the performance and operating lifetime of functional materials. The displacement per atom (DPA), which is the average number of times an atom is displaced from the original lattice, is widely used as an exposure parameter to evaluate the atomic-level structural damage in irradiated materials.

In the irradiation effect study on superconducting materials, DPA number has been used to predict the properties change of material [1, 2]. There is evidence that DPA can be used to correlate the radiation effect, such as the change of the critical current and the critical temperature of superconductors, under different irradiation particles [3]. In the study of ion irradiation effects on the Yttrium Barium Copper Oxide (YBCO) high temperature superconducting films or tapes, SRIM [4] is commonly used to calculate DPA [58]. For uncharged particles, only DPA calculating method of YBCO superconductors under gamma irradiation has been reported [9]. In the recent years, the research on neutron irradiation effects of YBCO superconducting tapes has been of interest [10, 11]. The study on the DPA number calculation of YBCO superconducting tapes under neutron irradiation has practical significance.

The Norgett–Robinson–Torrens (NRT) [12] DPA has been the standard evaluation parameter for several decades. DPA calculations in this paper only involve the NRT-DPA, which measures the primary radiation damage and does not account for the annealing of displacement cascades.

For neutron radiation damage calculation, total integrated damage energy is often obtained by folding the energy-dependent irradiation spectrum with energy-dependent neutron displacement cross sections σd. The most common method to obtain σd is to use the nuclear data processing code NJOY [13], which applies the Lindhard, Scharff, and Schiøtt (LSS) partition theory [14] to evaluate the damage energy directly. Then, the number of atomic displacements and DPA can be calculated according to the value of damage energy by NRT formula. SPECTRA-PKA provides more useful information about the contribution of each different reaction channel to the final DPA number obtained by convoluting the energy-dependent cross-sectional data of recoils and neutron flux spectrum [15]. However, the above DPA calculating methods ignore the spatial distribution and scattering angle distribution of PKAs and assume that recoils are stopped in the material of interest. Under this assumption, the detailed spatial distribution of damage energy deposited by every individual recoil is not taken into consideration and only the total value of damage energy is obtained. However, for models with small thickness, spatial distribution of damage energy deposited by PKAs is of importance and may affect the final DPA rate results seriously.

Besides the high temperature superconducting (HTS) tape model, the multilayer nanocomposites model in reactors [16] and the silicon drift detector model in deep space [17] both contain the geometric structures with thickness from micrometer to nanometer scale. In fact, it is necessary to study on the applicability of radiation damage method for models with small thickness under neutron irradiation.

The applicability of the analytic formulas assuming that all recoils are stopped in the material of interest for primary radiation damage calculation of thin plate models under neutron irradiation is first analyzed. Geant4 [18], a commonly applied open-source code package for the simulation of interactions between the particle and the matter, is used to analyze the energy and angular distribution of PKAs. In addition, Stopping and Range of Ions in Materials (SRIM) is used to predict the FPs number and damage energy induced by primary recoils. In present work, the influence of thickness on the primary radiation damage calculation of thin plate model under neutron irradiation with energy from 1 MeV to 14 MeV is analyzed by using Geant4 and SRIM. The deviation of analytic method for HTS tape model with different neutron energies is also calculated.

2. Methodology

2.1. PKAs Counts and Analysis

Geant4, an open-source toolkit written in C++, is developed for the simulation of the transportation of particles through matter. Geant4 tracks every individual particle from creation to disappearance by random sampling methods. The types of interactions and the parameters of recoils depend on the corresponding probability distributions based on nuclear databases. Geant4 only allows using the evaluated nuclear data libraries in the G4NDL format. G4NDL4.6 based on JEFF-3 [19] is developed by Mendoza et al. [20] using their self-produced program. In this research, Geant4 4.10.6 and G4NDL4.6 neutron cross sections were applied for obtaining the angular, energy, and spatial distributions of recoil ions.

The predefined physics list QGSP_BIC_HP was used without any modification. QGSP_BIC_HP is suitable for neutron processes with the energy below 20 MeV. G4NeutronHP package is included in QGSP_BIC_HP. The package is used in simulations involving neutrons with the energy from 0.025 eV to 20 MeV. Three different processes of neutrons are included: elastic, capture, and inelastic processes.

In the process of elastic scattering, the relationship between the energy and angular direction of recoils can be expressed in the following theoretical equation:where and denote the energy of recoil and the energy of incident neutron, respectively. is the atomic mass of recoil and is the angle between the direction of recoil and the direction of the incident neutron.

We rewrote the UserSteppingAction function to record the location, direction, and energy of every recoil in a ROOT file [21].

In our work, thin plate models with different materials are calculated. The neutron source is monoenergetic and the incidence direction is parallel to the depth direction of the model. The spatial distribution of neutron source is uniformed. The material outside the target box is vacuum. To illustrate the effect of thickness on displacements calculation, materials with different atomic number and atomic mass are chosen.

The information recorded with Geant4 is the location, direction, and energy of each PKA. When the thickness of target decreases to the micrometer scale, it will cost a lot of computing time to collect enough PKAs. For example, only about 103 PKAs can be recorded in a 1 µm thick Si model after 108 neutrons with the energy of 14 MeV passing through the material. If the thickness is in nanometer scale, computing time will be unacceptable. In fact, when the thickness of model is much smaller than the mean free path of neutrons, all neutrons collide almost randomly and tend to collide just once. So, for models with the thickness less than dozens of micrometers, we assume that the spatial distribution of PKAs is uniform along the depth of the model and all PKAs are generated from the first collision of neutrons based on the properties of the Monte Carlo method.

In summary, we recorded the location, direction, and energy of the PKAs generated in a plate model with the thickness of dozens of micrometers and then removed the part that was created by secondary collisions. The last step was to scale down the coordinate values of PKAs along the depth direction of the model to a required thickness. By this method, it is feasible to get a sufficient number of PKAs even for a target with the thickness of several nanometers.

2.2. Displacement Damage Calculation
2.2.1. SRIM Simulation

SRIM-2008 (Stopping and Range of Ions in Matter) is a broadly used Monte Carlo code for computing the range of ions and the radiation damage exposure parameters such as the number of displaced atoms (Nd) and DPA. There are two primary computing methods in SRIM: Ion Distribution and Quick Calculation of Damage method and Detailed Calculation with Full Damage Cascades (F-C) method. The former just tracks every PKA and assumes that there is a linear relationship between Nd and the damage energy deposited by PKAs, which is based on the Kinchin-Pease (K-P) method [22]. Meanwhile the latter method tracks every PKA and all particles obtain energy by collisions until the energy falls below the threshold energy and the number of Nd can be obtained by summing the numbers of vacancies, interstitials, and replacements.

F-C method is recommended in SRIM manual. However, Stoller et al. recommend using K-P method and setting the lattice binding energy to zero, then extracting the damage energy from a specific output file, and finally calculating the corresponding Nd according to NRT equation [23]. F-C method was not recommended based on the fact that Nd calculated by F-C is about twice as large as that by K-P according to vacancy.txt output file, while their damage energy values are close. In fact, when the K-P option is chosen and the lattice binding energy is set to zero, Nd calculated by NRT equation is very close to the corresponding values in vacancy.txt. Before SRIM calculation, we extracted the location, direction, and energy values of PKAs from .ROOT file and transferred them to a specific format by a self-developed data processing ROOT function according to the input template file TRIM.DAT as the ion source, as seen in Figure 1.

2.2.2. Lindhard–Robinson Analytic Method

The damage energy induced by neutrons can be also directly assessed according to the formulas described by Robinson based on Lindhard theory [12]:where is the damage energy, which represents the energy available to generate atomic displacements by the elastic collisions between atoms in the lattices. is a parameter presenting the constants of the Thomas-Fermi description of atomic interactions. Both k and are the functions of the mass and atomic number of the incident atom and target atom. T denotes the initial energy of PKA.

It should be noted that damage energy and corresponding Nd are determined from the analytic formula assuming that recoiling atoms are stopped in the material of interest.

The biggest difference between the Monte Carlo method and the analytic method is that the analytic method ignores the detailed spatial distribution of damage energy. However, the analytic method is widely used in traditional software for DPA calculation, like the damage cross section module in NJOY or DPA calculation module in SPECTRA-PKA.

3. Results and Discussion

3.1. Thin Plate Model with the Single Material
3.1.1. PKA Analysis

In order to further analyze PKAs, we simulated and extracted the angle and energy of PKAs generated by neutrons with different energy using Geant4 software packages. The number of incident particles that we simulated is 108. The efficiency of the calculation depends on the number of threads and the speed of processors. The common silicon is chosen as the material of the model with a thickness of 80 µm to illustrate the properties of PKAs, as shown in Figure 2. Hydrogen and helium recoils are generated in the inelastic neutron scattering. They can only produce few displaced atoms because of small Tdam but can have a great influence on the shape of energy distribution of PKAs. For model with small thickness, they have little effect on the result of displacements. So, the hydrogen and helium recoils are removed manually in this figure. PKAs are separated into two parts, elastic PKAs and inelastic PKAs, according to the processes that generate them. Elastic PKAs are the recoils that generated in elastic collision process. Other PKAs are denoted as inelastic PKAs. It is shown that the proportion of elastic PKAs varies with the incident neutron energy in Table 1.

The results in Table 1 show that the proportion of the elastic PKAs tends to decrease as the neutron energy increases for silicon, while the tendency for inelastic PKAs is opposite. In fact, the probabilities of different reaction types depend on the neutron cross sections of silicon. As the energy of incident neutron is relatively low, the elastic collision process is predominant and the number of elastic PKAs is higher than that of inelastic PKAs. When the neutron energy decreases to 1 MeV, only the elastic scattering happens. When the energy of neutron increases to 10 MeV, the number of inelastic PKAs is more than that of elastic PKAs. It indicates that the number of inelastic PKAs is considerable when the target atoms are collided by neutrons with relatively high energy.

In Figure 2, the PKA counts are 64860, 72997, and 73074 for neutrons with energies of 14 MeV, 10 MeV, and 6 MeV, respectively. The scattering angle in Figure 2 is the cosine of the angle between the direction of PKA and the direction of incident neutrons. The energy-angle distribution of elastic PKAs shows that the energy and angle of PKAs change along a specific curve. As shown in Figures 2(a), 2(c), and 2(e), the vast majority of elastic PKAs have relatively low energy and the scattering angles are between 78 and 90°, which means the directions of PKAs are nearly perpendicular to the direction of incident neutron. The numerical relationship between the energy and direction of elastic PKAs agrees well with equation (1). Moreover, the maximum value of energy and the uniformity of energy-angle distribution will decrease if the neutron energy increases.

However, most of inelastic PKAs have relatively higher energy and their directions are close to the incident direction of neutrons, as shown in Figure 2. When the energy of incident neutron is 14 MeV, the cosine of the angle between the direction of the incident neutron and the corresponding PKA is even larger than 0.8 and the energy is about 1 MeV for most circumstances, as shown in Figure 2(b).

It can be predicted that the inelastic PKAs have the potential to move across the boundaries if the thickness of a model is small enough. If a PKA leaves the model, the true deposited damage energy will be less than Tdam calculated by equation (2). Therefore, it is necessary to evaluate the applicability of traditional Tdam calculating method for models with small thickness.

3.1.2. Displacement Damage Calculation

PKAs are generated randomly in the collisions between neutrons and target atoms, so every PKA has different position, angular direction, and energy. It is hard to evaluate the valid part of Tdam by an analytic method directly. PKAs cannot be tracked directly with Geant4 unless the corresponding screened Coulomb scattering code is developed and well verified, because all the energy of the PKA is considered to be deposited in the position at which the particle is generated. In our work, SRIM is applied, since it can simulate the spatial distribution of deposited damage energy caused by every PKA. Since the parallel processing is not supported by SRIM, the efficiency of tracking PKAs only depends on the speed of processors. The influence of the thickness on displacement damage calculation in thin plate models with the material of silicon and niobium under different neutron energies is analyzed in Figure 3. The values of Ed of silicon and niobium are 37 eV and 78 eV, respectively [24].

Ns denotes Nd calculated by SRIM in a plate model with a limited thickness, while NI denotes Nd produced under the assumption that PKAs stay in the model. NI is consistent with the result of analytic formula. The value of Ns/NI represents the deviation between the results considering the effect of thickness and the result of analytic formula. So, higher Ns/NI represents better suitability of analytic method for displacement damage calculation. The corresponding error bars are shown in the figure. The errors are mainly from the sampling of the PKAs and the statistical error of SRIM. The results indicate that Ns/NI decreases rapidly when the thickness decreases to 1 µm and the deviation increases with increase of the neutron energy in most instances for silicon and niobium plate model.

When the thickness decreases, the number of vacancies caused by PKAs is about half of NI under the radiation of 14 MeV neutrons, as shown in Figure 3(a). The errors are smaller than 0.5%, but there is no apparent change tendency for the silicon model. The neutron cross section is a function of energy, so the proportion of inelastic PKAs varies with incident neutrons energy. As mentioned before, inelastic PKAs have relatively higher energy and larger cosine of scattering angle. If the inelastic PKAs become predominant, more PKAs will tend to leave the model, which make the assumption of theoretical formulas invalid. So, if the neutron cross sections of inelastic processes increase with the increase of neutron energy for the material, Ns/NI will decrease. The results show that the thickness of model has obvious influence on the value of Ns/NI. When the thickness of silicon model is 50 nm, Ns/NI is less than 0.5 even if all PKAs are produced in elastic scattering process.

The results of the niobium plate model are shown in Figure 3(b), which is similar to Figure 3(a) in shape. But the results of Ns/NI in Figure 3(b) are larger than those in Figure 3(a). In fact, the angular and energy distributions of PKAs depend on the nature of the incident particle and the properties of the material. The ranges of PKAs also have a great influence on the final results. For example, when the target atoms are collided by the 14 MeV incident neutrons and the elastic scattering occurs, the maximum energies of silicon and niobium recoils are 1.86 MeV and 0.59 MeV, respectively. In addition, the range of silicon recoil is larger than the range of niobium recoil when they have the same energy. The range of 1 MeV silicon recoil in silicon target is 1.13 µm and the range of 1 MeV niobium recoil in niobium target is only 0.22 µm. For the heavier recoils, it is harder to leave the plate model, so the thickness has less effect on the results of displacements calculation. Errors are smaller than 0.8% in the niobium plate model. The error bars show that the errors for 1 MeV are relatively larger than the errors in other energy points. It is due to small NI for this case. NI is less than 50 vacancies/ion for niobium, while NI for silicon is more than 200 vacancies/ion under 1 MeV neutron irradiation. NI will be smaller if the neutron energy is lower or Ed is larger. In fact, the same deviation of Ns may make the error of Ns/NI larger when NI becomes smaller. So, the value of NI has effect on the error of the result.

In summary, Ns/NI varies with incident neutron and thickness for models with small thickness; the applicability of analytic methods needs to be evaluated according to specific material. When the thickness is less than 1 µm, the detailed distribution of damage energy should be taken into consideration.

3.2. HTS Tape Model

Multilayer thin plate models are common in neutron irradiation environments. HTS materials have the potential to be applied in fusion reactor to produce the required magnetic field. The model of the 2nd generation composite HTS tape is shown in Figure 4 [25]. The YBCO tape consists of upper Cu layer of 20 µm thickness, Ag layer of 2 µm thickness, YBCO layer of 1 µm, buffer layers of about 70 nm, Hastelloy C276 layer of 50 µm, lower Ag layer of 2 µm, and lower Cu layer of 20 µm. The stoichiometric coefficients for yttrium, barium, copper, and oxygen are 1, 2, 3, and 7, respectively, in YBCO layer. The damage condition of HTS layer, which directly relates to the lifetime of superconducting facilities, has been an issue of great interest. Copper oxide superconductors have been found to be very sensitive to damage caused by irradiation. It is critical to calculate the DPA accurately as exposure parameters in superconducting layer to predict property changes of the functional material.

The density of YBCO layer is 6.54 g/cm3. In the calculation, the incident neutrons are monoenergetic and the direction is shown in Figure 4. Since the total thickness is much smaller than the mean free path of neutrons, each neutron tends to collide just once in the whole volume of the sample. Although the collision happens randomly, the total number of collisions in each layer depends on the neutron cross sections of the material. For the multilayer model, we evaluate Ns/NI in HTS layer under the irradiation of neutrons with different energies. PKAs are generated in the reactions between neutrons and target atoms. PKAs, which come from yttrium, barium, copper, and oxygen lattice atoms, are recorded separately by Geant4. The PKAs generated by neutrons are recorded at 14 different energy points from 1 MeV to 14 MeV. The counts of PKAs generated by neutrons with different energies are shown in Figure 5. The x-axis represents the energy of the incident neutron. The y-axis represents the number of PKAs.

The results show that the contribution of PKAs with different atomic species to the total number of PKAs varies obviously with the energy of incident neutrons. It depends on the neutron-atom interaction cross sections of JEFF-3 in different energy point.

For the compound material, each type of atoms may have different Ed and nk. Ed and nk of different atomic species in YBa2Cu3O7 are shown in Table 2 and nk is the relative fraction of the k-atom in its crystalline sublattice [9].

To evaluate the contribution of PKAs with different atomic species to the total value of Nd, Nd is calculated separately and the results of Ns/NI are compared in Figure 6.

The x-axis represents the energy of incident neutron. The y-axis represents the results of Ns/NI. For this compound material, YBa2Cu3O7, Ns/NI caused by yttrium, barium, and copper recoils tends to decrease when the neutron energy increases, as the variation trend of Si and Nb mentioned before. However, the results of oxygen recoils show that Ns/NI increases obviously when the neutron energy increases to 4 MeV or 12 MeV. In fact, it depends on the ratio of the elastic scattering cross sections to the total neutron cross sections of oxygen at different energy points. The results of final Ns/NI and the corresponding statistical errors in YBCO layer are shown in Figure 7.

Final Ns/NI and the errors in the compound material are affected by the number, the value of Ns/NI, and the corresponding error for each type of PKAs. The results of oxygen recoils have the biggest effect on the final result because Ns/NI of oxygen is much lower than that of other types of PKAs, and relatively more oxygen recoils are generated, as shown in Figure 5. The errors are less than 0.4%. The results show that Ns/NI in 1 MeV and 6 MeV is much lower than that at other energy points, which is different from the trend of results of oxygen recoils, because the number of oxygen recoils is relatively larger, as shown in Figure 5. So, the actual damage energy deposited in superconducting layer is obviously lower than the damage energy calculated with the assumption that PKAs stay in the material of interest and the value of Ns/NI is 0.84 to 0.93 with the neutron energy from 1 MeV to 14 MeV for the thin multilayer model.

In fact, the widely used analytic formulas for neutron radiation damage could overestimate the value of damage energy and DPA for models with small size, because the assumption of analytic formulas (all recoils are stopped in the material of interest) is ignored in most instances. The effect of the size of model on the damage energy calculation is complicated, because every PKA has different energy, direction, and location. In compounds, the neutron cross sections of different nuclide also affect the final results. So, the effect of the size and incident energy on the calculation of DPA number needs to be evaluated according to the specific radiation environment. The number of PKAs that can be tracked is less than 105, which is limited by the internal setting of SRIM.

4. Conclusion

If the analytic formulas based on Lindhard–Robinson theory are directly used to calculate the DPA of the thin plate model in neutron radiation environment, the validity of the assumption that the energies of all PKAs are deposited in the material of interest will be affected by the size of model and the energy of incident neutrons.

The spatial and angular distributions of PKAs need to be taken into account and tracked continuously to evaluate the applicability of traditional analytic formulas for models with thickness in micrometer scale. When the thickness is less than 10 µm, the deviation of analytic formulas increases rapidly for the smaller thickness. For different incident neutron energies, the deviation of analytic formulas depends on the neutron-atom interaction cross sections, which determine the ratio of the number of elastic PKAs to the total number. For the HTS tape model with the 1 µm superconducting layer, the final results are strongly affected by the neutron cross sections of oxygen atoms, which needs to be evaluated according to the specific radiation environment.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this study.

Acknowledgments

This study was supported by Program of Natural Science Foundation for Distinguished Young Scholars of Anhui Province (no. 2008085J28), the National Natural Science Foundation of China (General Program) (no. 52077211), and Outstanding Member of the Youth Innovation Promotion Association of the Chinese Academy of Sciences (YIPA, CAS) (no. Y201979).