Study on the Applicability of Neutron Radiation Damage Method Used for High-Temperature Superconducting Tape Based on Geant4 and SRIM

High-temperature superconducting material is a promising candidate to fabricate superconducting magnet for magnetic confinement fusion reactors. -e 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. -e 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. -e radiation damage calculation of simple thin plate models with different thicknesses and the tape model are conducted with the neutron energies from 1 to 14MeV. -e 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.


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. e 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 (E d ) have the potential to cause the displacement of lattice atoms and create considerable amounts of vacancies and interstitials, also called Frenkel Pairs (FPs). e accumulation of the FPs and other consequences of irradiation damage can significantly degrade the performance and operating lifetime of functional materials. e 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]. ere 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 [5][6][7][8]. 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]. e study on the DPA number calculation of YBCO superconducting tapes under neutron irradiation has practical significance.
e 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 energydependent neutron displacement cross sections σ d . e 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. en, 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. e 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. e deviation of analytic method for HTS tape model with different neutron energies is also calculated.

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. e 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. e 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. e package is used in simulations involving neutrons with the energy from 0.025 eV to 20 MeV. ree 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 E r and E n denote the energy of recoil and the energy of incident neutron, respectively. A is the atomic mass of recoil and u 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. e neutron source is monoenergetic and the incidence direction is parallel to the depth direction of the model. e spatial distribution of neutron source is uniformed. e 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. e 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 10 3 PKAs can be recorded in a 1 µm thick Si model after 10 8 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. e 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. method. e former just tracks every PKA and assumes that there is a linear relationship between N d 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 N d can be obtained by summing the numbers of vacancies, interstitials, and replacements.

Displacement
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 N d according to NRT equation [23]. F-C method was not recommended based on the fact that N d 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, N d 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.

Lindhard-Robinson Analytic Method.
e damage energy induced by neutrons can be also directly assessed according to the formulas described by Robinson based on Lindhard theory [12]: where T dam is the damage energy, which represents the energy available to generate atomic displacements by the elastic collisions between atoms in the lattices. k is a parameter presenting the constants of the omas-Fermi description of atomic interactions. Both k andg(ε) 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 T dam and corresponding N d are determined from the analytic formula assuming that recoiling atoms are stopped in the material of interest.
e 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.

in 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. e number of incident particles that we simulated is 10 8 . e efficiency of the calculation depends on the number of threads and the speed of processors. e 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. ey can only produce few displaced atoms because of small T dam 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.
e 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   respectively. e scattering angle in Figure 2 is the cosine of the angle between the direction of PKA and the direction of incident neutrons. e 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. e 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 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 T dam calculated by equation (2). erefore, it is necessary to evaluate the applicability of traditional T dam calculating method for models with small thickness.

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 T dam 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. e 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. e values of E d of silicon and niobium are 37 eV and 78 eV, respectively [24]. N s denotes N d calculated by SRIM in a plate model with a limited thickness, while N I denotes N d produced under the assumption that PKAs stay in the model. N I is consistent with the result of analytic formula. e value of N s /N I represents the deviation between the results considering the effect of thickness and the result of analytic formula. So, higher N s /N I represents better suitability of analytic method for displacement damage calculation. e corresponding error bars are shown in the figure. e errors are mainly from the sampling of the PKAs and the statistical error of SRIM. e results indicate that N s /N I 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 N I under the radiation of 14 MeV neutrons, as shown in Figure 3(a). e errors are smaller than 0.5%, but there is no apparent change tendency for the silicon model. e 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, N s /N I will decrease.
e results show that the thickness of model has obvious influence on the value of N s /N I . When the thickness of silicon model is 50 nm, N s /N I is less than 0.5 even if all PKAs are produced in elastic scattering process. e 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 N s /N I 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. e 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. e 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. e error bars show that the errors for 1 MeV are relatively larger than the errors in other energy points. It is due to small N I for this case. N I is less than 50 vacancies/ion for niobium, while N I for silicon is more than 200 vacancies/ion under 1 MeV neutron irradiation. N I will be smaller if the neutron energy is lower or E d is larger. In fact, the same deviation of N s may make the error of N s /N I larger when N I becomes smaller. So, the value of N I has effect on the error of the result.
In summary, N s /N I 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 Science and Technology of Nuclear Installations 5 detailed distribution of damage energy should be taken into consideration.

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. e model of the 2nd generation composite HTS tape is shown in Figure 4 [25]. e 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. e stoichiometric coefficients for yttrium, barium, copper, and oxygen are 1, 2, 3, and 7, respectively, in YBCO layer. e 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. e density of YBCO layer is 6.54 g/cm 3 . 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 N s /N I 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. e PKAs generated by neutrons are recorded at 14 different energy points from 1 MeV to 14 MeV. e counts of PKAs generated by neutrons with different energies are shown in Figure 5. e x-axis represents the energy of the incident neutron. e yaxis represents the number of PKAs.
e 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 E d and n k . E d and n k of different atomic species in YBa 2 Cu 3 O 7 are shown in Table 2 and n k 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 N d , N d is calculated separately and the results of N s /N I are compared in Figure 6. e x-axis represents the energy of incident neutron. e y-axis represents the results of N s /N I . For this compound material, YBa 2 Cu 3 O 7 , N s /N I 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 N s / N I increases obviously when the neutron energy increases to    Figure 7. Final N s /N I and the errors in the compound material are affected by the number, the value of N s /N I , and the corresponding error for each type of PKAs. e results of oxygen recoils have the biggest effect on the final result because N s / N I of oxygen is much lower than that of other types of PKAs, and relatively more oxygen recoils are generated, as shown in Figure 5. e errors are less than 0.4%. e results show that N s /N I 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 N s /N I 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. e 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. e number of PKAs that can be tracked is less than 10 5 , which is limited by the internal setting of SRIM.

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    Science and Technology of Nuclear Installations material of interest will be affected by the size of model and the energy of incident neutrons. e 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
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this study.