Modeling of aluminum erosion by neutral particles in dedicated EAST experiments using the 3D-GAPS code

To evaluate the erosion effects of charge exchange (CX) neutrals on first wall materials, aluminum (Al) coated samples installed in a cylindrical sleeve were exposed to EAST plasmas with the Material and Plasma Evaluation System (MAPES). The 3D-GAPS code was firstly used to simulate the dedicated experiments with the neutral energy spectrums measured by the Low Energy Neutral Particle Analyzer located at the MAPES. The modeled Al erosion rates with isotropic angular distributions for incident D neutrals are consistent with the measured erosion thicknesses of Al layer by Rutherford backscattering spectrometry before and after the plasma exposure. The effects of neutral reflection coefficients and the heights of cylindrical sleeve on material erosion are also studied. The measured material erosion rate by CX neutrals in EAST verified by the modeling can help to us understand the effects of neutrals on first wall erosion, which is an important erosion process in future fusion reactors.


Introduction
In tokamaks, plasma-wall interaction will cause material erosion of the plasma-facing components (PFCs), affecting the material lifetime of the PFCs, and eroded particles entering into plasma will degrade the plasma performance [1]. In addition, eroded particles may be migrated and co-deposited with fuel particles on the PFCs, which can increase fuel retention in * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. fusion devices. Due to safety problems, tritium (T) retention has to be limited in ITER and future fusion reactors. In addition to ions that move along the magnetic field lines, neutrals generated by the recombination of ions and electrons or charge exchange (CX) reactions are not constrained by magnetic field lines and directly bombard PFCs, especially at the so-called magnetic shadowed areas that are inaccessible to ions. For ITER and future fusion reactors, the neutral flux and energy will be much higher than in the current machines, which can significantly affect the material lifetime and T retention. However, the fluxes, energy, and angular distributions of neutrals on the first-wall materials are still unclear so it is hard to make accurate predictions for future reactors. Therefore, it is necessary to study material erosion induced by neutrals in present tokamaks.
In JET [2][3][4][5], several material exposure experiments with samples installed in first-wall tiles between inner wall guard limiters were carried out during 1995-2014 and strong net erosion induced mainly by CX neutrals was found. In ASDEX Upgrade [6], a material exposure experiment with four heat shield tiles coated with tungsten was carried out, which showed higher erosion rates than calculated ones indicating that CX neutrals played an important role in material erosion. After the experiments, campaign-averaged erosion rates were obtained but could not truly reflect the erosion and deposition processes induced by neutrals because the results could be influenced by many factors during the long-term experiments, such as disruptions, wall-conditioning, etc.
In EAST, a material migration experiment with an ITER first-wall panel proxy was carried out [7], which showed a net erosion in the magnetic shadowed area. This indicates material erosion contributed by neutrals. After that, dedicated experiments for studying material erosion by neutrals were carried out for the first time in EAST [8]. The spatial distributions of neutral-induced material erosion rates were obtained by comparing the Rutherford backscattering spectrometry (RBS) measurements of the material coating thickness before and after the exposure. Recently, the Low Energy Neutral Particle Analyzer (LENPA) was successfully installed at the mid-plane of sector H in EAST, which can provide the neutral energy spectrum [9,10]. Based on the experimental results and the LENPA diagnostics, it is possible to carry out simulations for the dedicated experiments. In this work, the 3D-GAPS code has been firstly used to simulate the material erosion by neutral particles in the EAST experiments.
The paper is organized as follows: section 2 gives a brief introduction on the EAST experiments and the LENPA diagnostics. Section 3 describes the 3D-GAPS code and the model setup. In section 4, the modeling results are shown in comparison with the experimental ones as well as the key impacting parameters. At last, a summary is given in section 5.

Experiments at EAST
The dedicated material erosion experiments by neutral particles in EAST are briefly summarized here and more details can be found in [8]. Due to the toxicity of Beryllium (Be) which is used as the ITER first-wall material, Aluminum (Al) was used to serve as a proxy for Be in the EAST experiments [11]. In the experiments, an Al coating with a thickness of 20-30 nm was deposited on a square silicon (Si) substrate with an area of 10 × 10 mm 2 by magnetron sputtering. The samples were placed at the center of a Molybdenum (Mo) cylindrical sleeve with an inner diameter of 50 mm and a height of 25 mm, which can protect the sample from incidence of charged particles as shown in figure 1. Two sets of samples were exposed in 2018 and 2019 campaigns, which are called the 1st and 2nd experiment in the following, respectively. With the help of Material and Plasma Evaluation System (MAPES) [7] located at the mid-plane of sector H, samples can be moved flexibly along the radial direction, which, for example, helps to avoid the impact by wall-conditioning. The samples were installed at 50 mm and 30 mm radially behind the limiter in the 1st and 2nd experiment for 849 s and 1014 s, respectively. The spatial distributions of Al erosion rates were obtained precisely at five locations in the shape of a cross by RBS measurement before and after exposures. The RBS measurement was performed on a tandem accelerator with 2 MeV 4 He [8]. The backscattered 4 He was detected by a gold-Si surface barrier detector placed at a scattering angle of 165 • . The thickness of Al coating was calculated by the yield of 4 He scattered from the Al layer. The averaged net erosion rates were 1.18 × 10 13 cm −2 s −1 and 4.64 × 10 13 cm −2 s −1 with an estimated uncertainty of about 4% in the 1st and 2nd experiment, respectively.
In the last phase of the second campaign in 2019, a timeof-flight LENPA was successfully installed at the mid-plane of sector H where the MAPES is located, which can measure neutral energy E and flux Γ p to the first wall in the energy range of 20-3000 eV [10]. The LENPA works in pulse-counting mode and the flight time of neutrals is used to calculate the neutral energy E. Therefore, the neutral energy spectrums dΓ p (E)/(dEdΩ) can be obtained and the solid angle of the LENPA detection is 6.9 × 10 −6 sr.
However, since the LENPA was not ready yet during the dedicated material erosion experiments, the neutral energy spectrums for the particular experiments could not be measured. Therefore, the LENPA data measured in the discharges with similar plasma parameters as the material erosion experiments but performed later are used for the modeling analysis. For the 1st and 2nd experiment, the averaged plasma density is around 3.8 × 10 19 m −3 and 5.0 × 10 19 m −3 and the heating power is around 1.6 MW and 3.1 MW, respectively [8]. In this work, the LENPA diagnostics in the discharge #101358 with the line-averaged density of 2.5 × 10 19 m −3 and the heating power of 2.0 MW is selected as representative for the 1st experiment. For the 2nd experiment, the discharge #103957 with the line-averaged density of 5.4 × 10 19 m −3 and the heating power of 3.3 MW is used. Figure 2 shows the neutral energy spectrums measured by the LENPA diagnostics for the two discharges. The mean energy of neutral particles is around 400 eV for #101358 and 500 eV for #103957. According to the recent analysis on the LENPA data in EAST [12], the integrated neutral flux increases with increasing heating power almost linearly, while it is less sensitive to the variation of the plasma density when the line-averaged density is higher than about 2.5 × 10 19 m −3 .

Model setup
The 3D-GAPS code is a Monte-Carlo neutral transport code, which can simulate impurity transport, material erosion and deposition processes in three-dimensional geometry including magnetic shadowed areas. Therefore, the 3D-GAPS code is also suitable for modeling of dedicated material erosion experiments by neutral particles in EAST. The surfaces of PFCs and simulation box boundaries are described by an arbitrary number of geometrical objects in the code. The code can simulate the transport of neutrals and their interaction with material surfaces including physical sputtering, chemical erosion, reflection, and deposition. It is assumed that neutrals simulated by test particles are emitted from a certain source and move in the simulation volume along straight lines without energy decay until they escape from the simulation box or reach a certain material surface (if elastic collisions with residual neutral gas are not considered). When particles reach material surface, they may be either reflected or deposited and cause surface erosion. To obtain material erosion and deposition rates, the sputtering yields are calculated based on the Bohdansky-Yamamura formula [13][14][15]. Table 1 gives the used parameters for sputtering yields calculated with the surface binding energy E s , the mass ratio of material to neutral M 2 /M 1 , the atom number of neutral Z 1 , the atom number of material Z 2 , the Thomas-Fermi energy E TF , the sputtering threshold energy E th , the fitting parameter defining the maximum of the sputtering yield curve Q and the target density n. Sputtered particles are assumed to leave the material surface with cosine angular distribution. The reflection coefficients and angles are calculated based on the SDTrimSP code data [16] or given by predefined constant values with the cosine angular distribution. The sputtering yields of Be and Al by D bombardment and self-sputtering calculated based on the Bohdansky-Yamamura formula are shown in figure 3. Pure Al is assumed for the sputtering of Al. The oxidation formation on Al surface may change the sputtering yield, which is not taken into account in the simulations. Figure 4 shows the reflection coefficients of D on Be and Al surfaces calculated based on the SDTrimSP code data. Figure 5 shows the model geometry of the Al sample located in the Mo cylindrical sleeve with the neutral injection plane above it. The neutral particles are launched from the injection plane with the energy distribution according to the measurements by LENPA diagnostic that can be seen from figure 2. The area of neutral injection plane is assumed to be large enough to take all the neutrals that could bombard the Al sample into account. Based on the neutral energy spectrum, the neutral injection flux is given by By integrating over the incident energy and the solid angle, the neutral injection fluxes for the 1st and 2nd experiment are 1.43 × 10 15 cm −2 s −1 and 5.88 × 10 15 cm −2 s −1 , respectively. It is assumed that the incident energy distribution is independent on the incident angle. The average neutral flux on the Al sample, Γ s , can be expressed as where Ω i is the integrated solid angle of neutrals to the Al sample considering the angular distribution. In the model, neutrals are assumed to be launched with isotropic or cosine angular distributions. At the center of the Al sample, the integrated solid angle of neutrals to the sample with the angular distribution of isotropy and cosine can be expressed as respectively, where b is the radius of projection area on the injection plane, r is the radial distance between micro element area dσ and the center of the projection area on the injection   plane, d is the distance between the sample and the neutral injection plane, φ and θ are the azimuth and polar angles of emitted neutrals, respectively. Meanwhile, the radius of projection area on the injection plane can be expressed as where a = 25 mm and h = 25 mm are the radius and height of the Mo cylindrical sleeve, respectively. By integrating over the radial distance and the azimuth angle, the solid angles of neutrals to the center of the Al samples with isotropic and cosine angular distributions can be calculated by equations (6) and (7), respectively, From the calculations, the integrated solid angles of neutrals to the center of the Al sample are independent of the distance between sample and neutral injection plane. For the whole Al sample, the distribution of calculated integrated solid angles on it with isotropic and cosine angular distributions are shown in figure 6. Figures 7 and 8 show the neutral fluxes on Al    respectively. The modeling results are in good agreement with the calculated neutral fluxes by equation (2).
It can be seen from figures 7 and 8 that the neutral fluxes on Al samples with cosine angular distribution are higher by more than 70% and more concentrated on the center than those with isotropic distribution, which are consistent with the integrated solid angles that can be seen from figure 6. Neutral fluxes with the two angular distributions show a gradual decrease from the center to the edge. Considering the size of Al sample and the ratio of height to radius of the Mo cylindrical sleeve, the neutrals can reach the Al sample only if the incident angle is smaller than 52 • . Neutrals with the incident angle smaller than 52 • based on cosine distribution contribute to more of the total injection fluxes than those based on isotropic distribution, which leads to more neutral bombardment on Al samples. As it can be seen from figure 6, the integrated solid angles at the center of Al samples are larger than those at the edge and the ratio of integrated solid angle between the center and the edge of Al samples with cosine distribution is larger than that with isotropic distribution, which leads to a higher gradient of neutral fluxes from the center to the edge with cosine distribution.

Modeling results and discussion
With the neutral energy spectrum and the sputtering yield, the 3D-GAPS code is then used to simulate the Al sputtering rates in comparison with the experiments. The Al erosion rates for the 1st and 2nd experiments are simulated with both isotropic and cosine angular distributions, as shown in figures 9 and 10. For the 1st experiment, the averaged simulated Al erosion rates are 1.51 × 10 13 cm −2 s −1 and 2.42 × 10 13 cm −2 s −1 with isotropic and cosine angular distributions, while they are 5.48 × 10 13 cm −2 s −1 and 8.79 × 10 13 cm −2 s −1 for the 2nd experiment, respectively. The simulation results with the isotropic angular distribution show better agreement with the measurements, which are 1.18 × 10 13 cm −2 s −1 and 4.64 × 10 13 cm −2 s −1 for the 1st and 2nd experiments, respectively. Figure 11 shows comparison of toroidal and poloidal profiles of Al sputtering rates between RBS measurement and the 3D-GAPS simulations with different assumptions of angular distribution. Obviously, the modeled profiles with the isotropic angular distribution are closer to the RBS results. The Al erosion rates at the center of the sample are lower than those at the edge from the RBS measurements while the modeling results show a gradual decrease from the center to the edge. The difference still needs to be further investigated. The averaged sputtering yields are very similar in all cases, which are from 0.03 to 0.035. The isotropic angular distribution is proven to be a more reasonable assumption. Therefore, the following discussions are based on the simulations with the assumption of isotropic angular distribution.
The Al erosion rates for the 2nd experiment are significantly higher than for the 1st experiment due to higher neutral flux on Al samples according to the LENPA measurement, that can be seen from figure 2. The averaged sputtering yields for the 2nd experiment are about 8% lower than the 1st experiment. This is probably because neutrals with lower energy below about 60 eV and higher energy above about 900 eV contribute more to the total injection flux in the 2nd experiment compared to the 1st experiment. Since the sputtering yields for the incident energy between 60 eV and 900 eV are higher than the other incident energy range, the sputtering yield for the 1nd experiment is slightly higher. This is illustrated in figure 12 that shows the contribution of particles with different incident energies to the total Al sputtering in each experiment. This contribution first increases with the incident energy and then decreases, with the peak value at around 75 eV and 65 eV for the 1st and 2nd experiments, respectively. Due to the shape of the energy spectra as shown in figure 2, the Al erosion is dominated by neutral particles with energy in the range from 20 to 500 eV. The D neutrals with energies below 500 eV contribute about 60.6% and 47.1% to the total Al erosion for the 1st and 2nd experiments. In the case of the 2nd experiments with a higher heating power, the fraction of higher energy neutrals is increased.

The effect of neutral reflection coefficients
In addition to the neutrals from the plasma, the neutrals reflected from the sides of the sleeve also can reach the   sample surface. Therefore, the neutral reflection coefficient can change the total flux of neutrals incident on the sample, thereby changing the sputtering rates. To study the effect of neutral reflection, the reflection coefficients of D on different materials are assumed with different constant values in the 3D-GAPS simulations. Figure 13 shows Al sputtering rates as a function of reflection coefficients. For the two experiments, the calculated Al sputtering rates from the SDTrimSP code data are consistent with the results based on an average constant reflection coefficient of about 0.6. The Al sputtering rates do not change significantly with increasing the reflection coefficient. When the reflection coefficient is set to 0, the Al sputtering rates caused by the first bombardment with neutrals can be calculated and these contribute to 89.4% and 88.3% of the  total erosion rates calculated based on the SDTrimSP code data for the 1st and 2nd experiments, respectively. Note that for the 2nd experiment, the Al sputtering rate increases more dramatically with higher reflection coefficients than that for the 1st experiment. This can be explained by higher mean energy of neutrals for the 2nd experiment. Although the probability of reflected neutrals bombarding on Al samples for the two experiments are similar with the same reflection coefficient, higher neutral energy can lead to a higher sputtering yield after several times of reflection. Since the Al sample is located at the bottom of the Mo cylindrical sleeve and the area of the Al sample exposed to reflected neutrals is much smaller than that of the Mo cylindrical sleeve, there is a very low probability for bombardments of reflected neutrals on the Al sample. Therefore, the effect of reflected neutrals on Al sputtering rate is weak.

The effect of height of Mo cylindrical sleeve
The integrated solid angle of neutrals to the sample depends on the ratio of height to radius of the Mo cylindrical sleeve. Therefore, changing the height of the Mo cylindrical sleeve will affect the incident neutral flux and the total reflection times of neutral particles. To identify this effect on material erosion, a scan of heights of the Mo cylindrical sleeve from 25 to 105 mm is performed by the 3D-GAPS modeling. Figure 14(a) shows the modeled Al sputtering rate as a function of sleeve height for the 1st and 2nd experiments. The Al sputtering rates decrease with increasing the height of the Mo cylindrical sleeve, which are consistent with the height dependence of integrated solid angle as shown in figure 14(b). Although neutrals are more difficult to escape from the higher sleeve, the average times of reflections per neutral do not increase significantly. However,  since the sample is located more deeply at the bottom with higher sleeve, the average time of bombardments on the Al sample by reflected neutrals decreases. The height of the Mo cylindrical sleeve not only affects the incident neutral flux on the Al sample, but also changes the angular distributions of incident particles. With a higher Mo cylindrical sleeve, the angular distribution of incident neutrals that can reach the Al sample moves towards the normal incidence, as shown in figure 15. This may reduce the sputtering yields slightly, which can be seen from figure 3(b).

Modeling of Be erosion in comparison with Al
Although Al can serve as a proxy for Be in the erosion studies [11], there are still differences in the sputtering yield and reflection coefficient between Be and Al that can be seen from figures 3 and 4. The sputtering yield of Be by D bombardment is lower with high incident energy. Here, Be erosion is calculated with the 3D-GAPS code using the same input as the previous Al modeling. The averaged simulated Be sputtering rates are 1.06 × 10 13 cm −2 s −1 and 3.64 × 10 13 cm −2 s −1 with isotropic angular distribution for the 1st and 2nd experiment, respectively. The D neutrals below 500 eV contribute 70.0% and 59.8% to the total Be erosion for the two experiments. The Be sputtering rates are slightly lower than those of Al, but are contributed more by lower energy neutrals, which is consistent with the difference between the sputtering yield of Be and Al.

Summary
In this work, the 3D-GAPS code has been used to simulate the Al erosion by neutral particles in dedicated EAST experiments. The neutral energy spectrums measured by the LENPA diagnostics during the discharges with similar plasma conditions as in the material exposure experiments are used for the modeling. The 3D geometry of the Mo cylindrical sleeve and different angular distributions of neutrals have been implemented in the 3D-GAPS code. The simulated incident fluxes with both isotropic and cosine angular distributions are consistent with the theoretic calculations. The simulated Al erosion rates with isotropic angular distributions are in good agreement with the experimental measurements. The Al erosion is dominated by neutral particles with energy in the range from 20 to 500 eV. The effects of neutral reflection coefficients and the heights of the Mo cylindrical sleeve are also discussed. The Al erosion rates decrease significantly for higher ratios of sleeve height to radius, mainly due to the reduced solid angle, from which neutral particles can reach the sample. For a hypothetical Be sample, simulated Be erosion rates under the same irradiation conditions are slightly lower than for the Al sample. In ITER, the CX neutral flux is more than one order of magnitude higher than the present EAST experiments, which leads to Be erosion rates also more than one order of magnitude higher according to the recent modeling results [17]. Be erosion rates by neutrals are in the same order of magnitude as those by ions. The good agreement between experiments and modeling verifies the LENPA measurements on neutral energy spectrums and the sputtering models by neutrals. In future work, systematic data analysis will be performed to study the effects of plasma conditions on the neutral distribution and the neutral-induced material erosion.