Quantum Electron Dynamics in Helium Ion Injection onto Tungsten Surfaces by Time-Dependent Density Functional Theory

Abstract

The neutralization process of an ion particle on a surface is key issue of the plasma-wall interaction (PWI). We investigated helium (He) ion injection onto a tungsten surface using time-dependent density functional theory (TDDFT) simulation. We developed the TDDFT code QUMASUN for this study, simulating the process of electron transfer from the surface to the He nucleus by solving the time evolution of electron wave function and the classical motion of nuclei simultaneously. Our results show that the probabilities of injected $\text{He}^{2+}$ changing into $\text{He}^{1+}$ and $\text{He}^{0}$ on the surface are approximately 40 percent and 13 percent, respectively. The captured electrons by $\text{He}^{1+}$ and $\text{He}^{0}$ predominantly occupy the 2s and 2p orbitals, corresponding to excited states. In addition, we stated some challenges for applying TDDFT to plasma-wall interactions.

Keywords: helium plasma, neutralization process, plasma-wall interaction, TDDFT, surface recombination

1 Introduction

The neutralization process of an ion particle at a solid surface is a fundamental process in plasma-wall interaction. For fusion science, it is called the surface recombination process, and is an important factor for not only edge plasma physics but also plasma confinement. Recently, the generation of heat due to surface recombination processes is relatively non-negligible in the heat load of the divertor [1].

Many experiments on the neutralization process of an ion particle on a surface have been performed [2, 3, 4]. However, there are few experiments at incident energy below 1 keV[5], which is important energy range in recent PWI. Furthermore, in the case of a helium (He) ion which is the subject of the present work, there are very few experiments of $\text{He}^{2+}$injection[6].

Theoretical explanations have claimed mainly neutralization by the Auger processes[7] and by the resonance processes[8]. However, these mixed processes have also been proposed[9], and it has not been clearly understood yet. Although it is often assumed that the neutralization occurs only in very near surface region, theory advanced to consider the change of potential for electrons due to the motion of the ion particle and detailed electronic state of surface[10].

In terms of recent simulation research, a multiscale approach that focuses on atomistic dynamics from various aspects with binary collision approximation (BCA), molecular dynamics (MD), and kinetic Monte-Carlo (KMC) and so on has been developed in PWI research field[11, 12, 13]. In BCA and MD, all particles are neutral atoms because of classical mechanics. In other words, an incident ion particle is also substituted by a neutral atom. Although density functional theory (DFT)[14] is currently often used to consider more detailed energy and potential, MD with DFT cannot treat an ion injection also. The MD with DFT on Born-Oppenheimer approximation represents the atomic dynamics in which electronic state is always a ground state. However, the system which has an incident ion particle and a solid material in one simulation box initially, even at larger distances, corresponds to an electronic excited state.

Time-dependent DFT (TDDFT)[15] can treat an ion injection onto a surface. In fact, Krasheninnikov and Miyamoto, et. al.[16, 17] have simulated a hydrogen ion injection and an argon ion injection by simultaneously solving the motions of nuclei and the dynamics of the electron wave function. In those simulations, the target surface was a monolayer of graphene which has a few electrons, probably due to the computational resources available at the time.

In the present work, we have challenged the TDDFT simulation of a single He ion injection onto the tungsten surface, in which the system has many number of electrons than that of the graphene, using recent computer resources. We deal with low-energy He ion injection because in the future we would like to commit issues of the He ash exhaust problem and the fuzz formation problem where the previous simulation suggest that sputtering by He ion is effective[18, 19]. The main purpose is to directly watch the neutralization process of the incident ion particle on a surface. For this purpose, we started with the development of the TDDFT code QUMASUN. The simulation methods are described in section 2. We discuss the results and interpretations of the simulations in section 3. In addition, as this is our first effort and not yet directly comparable to experiments, we aim to clarify what TDDFT can do and what challenges exist in dealing with this problem.

2 Simulation Method

2.1 TDDFT code QUMASUN

To treat the quantum dynamics of electrons in the PWI, we have developed the DFT and TDDFT code QUMASUN. In the present simulation with the QUMASUN, the time-dependent Kohn-Sham equation for electron and the classical equation of motion for nuclei are solved simultaneously[20].

The atomic unit system is used in the following statements. According to the time-dependent Kohn-Sham equation, the time evolution of electron wave function is governed by

$$i\frac{\partial\psi_{j}(x,t)}{\partial t}=\hat{H}_{\text{KS}}\psi_{j}(x,t),$$ (1)

$$\hat{H}_{\text{KS}}=-\frac{1}{2}\Delta+V_{\text{H}}[\rho]+V_{\text{xc}}[\rho]+V_{\text{l}}+V_{\text{nl}},$$ (2)

where $\psi_{j}(x,t)$ is the $j$ th electron orbital function, and the variable $x=(\vec{r},\sigma)$ is the set of the spatial coordinate $\vec{r}$ and the spin variable $\sigma$. The electron density $\rho$ is given by $\rho=\sum_{j}|\psi_{j}(x,t)|^{2}$. In the present simulation, the Hartree potential $V_{\text{H}}[\rho]$ is calculated by solving the Poisson equation, $-4\pi\rho=\Delta V_{\text{H}}$, under the periodic boundary condition. For the exchange-correlation potential $V_{\text{xc}}[\rho]$, the local spin density approximation (LSDA)[21] was used simply. In actual, the effect of core electrons are approximated by the pseudo-potentials of the local term $V_{\text{l}}$ and non-local term $V_{\text{nl}}$, and then only the valence electrons are solved by following the time-dependent Kohn-Sham equation (1). For the pseudo-potential $V_{\text{l}}$ and $V_{\text{nl}}$, we adopted the Morrison, Bylander, and Kleinman (MBK) pseudo-potential[22] which is a norm-conserving model based on the Vanderbilt’s ultra-soft pseudopotential[23]. Actually, the MBK pseudo-potential data set[24] were generated using the ADPACK tool for another DFT code OpenMX[25, 26]. In addition, for the present work, some parameters were modified from the original data set.

To apply the MBK pseudo-potential, which was generated in a spherical coordinate system for independent atom, into the present simulation, the local potential is converted under the periodic boundary condition as follows: First, the local term for nucleus $\alpha$, $\tilde{V}_{\text{l},\alpha}(\vec{r})$, is loaded from the pseudo-potential data-set, and then the core charge density of the nucleus $\alpha$, $\tilde{n}_{\alpha}(\vec{r})$, is calculated as $\tilde{n}_{\alpha}=\Delta\tilde{V}_{\text{l},\alpha}/4\pi$ in an independent system centered the position of nucleus $\alpha$. Second, the total core charge in simulation system with periodic boundary condition $n(\vec{r})$ is calculated as the conversion $n(\vec{r})=\sum_{\alpha}\tilde{n}_{\alpha}(\vec{r}-\vec{R}_{\alpha})$ where $\vec{R}_{\alpha}$ is the position of nucleus. This conversion correctly takes into account the folding at the boundary of the simulation box. Here, with the cutoff length $r_{\text{c}}$, $\tilde{n}_{\alpha}(\vec{r})=0$ and $\tilde{V}_{\text{l},\alpha}(\vec{r})=-{(Z_{\alpha}-Q_{\alpha})}/{|\vec{r}|}$ if $|\vec{r}|>r_{\text{c}}$, where $Q_{\alpha}$ is a valence charge. Third, the actual local potential $V_{\text{l}}(\vec{r})$ is obtained by solving the Poisson equation, $4\pi n(\vec{r})=\Delta V_{\text{l}}(\vec{r})$, in the actual simulation box under periodic boundary condition.

The equations of motions for nuclear positions $\vec{R}_{\alpha}$ and momenta $\vec{P}_{\alpha}$ are given by

$$\frac{d\vec{R}_{\alpha}}{dt} =\frac{P_{\alpha}}{m_{\alpha}},$$ (3)

$$\frac{d\vec{P}_{\alpha}}{dt} =-\frac{\partial}{\partial\vec{R}_{\alpha}}\left(E_{\text{KS}}[\rho]+E_{\text{nn}}\right),$$ (4)

where the $E_{\text{KS}}[\rho]$ is the Kohn-Sham total energy for electrons. The $E_{\text{nn}}$ is the nuclear-nuclear interaction energy. Since the present simulation deals with collisions of nuclei at higher energy than the energy range in general material physics and chemistry, the definition of $E_{\text{nn}}$ was carefully designed as follows:

$$E_{\text{nn}} =-\int V_{\text{l}}(\vec{r})n(\vec{r})d\vec{r}-E_{\text{self}}+\frac{1}{2}\sum_{\alpha}\sum_{\beta\neq\alpha}\left[-\epsilon_{\alpha\beta}(R_{\alpha\beta})+\epsilon^{\text{TF}}_{\alpha\beta}(R_{\alpha\beta})\right]\theta(2r_{c}-R_{\alpha\beta}),$$ (5)

$$\epsilon_{\alpha\beta}(R_{\alpha\beta}) =-\int\tilde{V}_{\text{l},\alpha}(\vec{r}-\vec{R}_{\alpha})\tilde{n}_{\beta}(\vec{r}-\vec{R}_{\beta})d\vec{r},$$ (6)

$$\epsilon^{\text{TF}}_{\alpha\beta}(R_{\alpha\beta}) =\frac{(Q_{\alpha}+(Z_{\alpha}-Q_{\alpha})e^{-R_{\alpha\beta}/C_{\alpha}})(Q_{\beta}+(Z_{\beta}-Q_{\beta})e^{-R_{\alpha\beta}/C_{\beta}})}{R_{\alpha\beta}}.$$ (7)

Here, the first term is calculated as the local potential $V_{l}(\vec{r})$ and the core charge density $n(\vec{r})$ introduced in the above manner. However, the first term includes self interaction of nuclei, and then the second term $E_{\text{self}}$ is the correction of the self interaction and is a constant. The third term is the correction for two-body core-core repulsion composed of $\epsilon_{\alpha\beta}(R_{\alpha\beta})$ and $\epsilon^{\text{TF}}_{\alpha\beta}(R_{\alpha\beta})$ where the distance between nuclei $R_{\alpha\beta}=|\vec{R}_{\alpha}-\vec{R}_{\beta}|$. By the step function $\theta(2r_{c}-R_{\alpha\beta})$, the third term is effective if $R_{\alpha\beta}<2r_{c}$. The $\epsilon_{\alpha\beta}(R_{\alpha\beta})$ is the term to cancel core-core interaction with smoothed core charge density generated by the local term of the pseudo-potential, and it was calculated on ellipsoidal coordinate in the simulation. The two-body repulsive potential $\epsilon^{\text{TF}}_{\alpha\beta}(R_{\alpha\beta})$ is modeled to asymptotically approach the Coulomb potential between nuclei, $Z_{\alpha}Z_{\beta}/R_{\alpha\beta}$, when $R_{\alpha\beta}\rightarrow 0$. The artificial parameters $C_{\alpha}=r_{c}Q_{\alpha}/2Z_{\alpha}$ and $C_{\beta}=r_{c}Q_{\beta}/2Z_{\beta}$. This correction for two-body core-core repulsion is important to represent the high energy collision between nuclei such as the present simulation because in general DFT the local term potential is smaller than the Coulomb potential. If the correction is not used, the projectile with high energy would pass though the target nuclei. However, there is room for improvement on the function form of $\epsilon^{\text{TF}}_{\alpha\beta}$.

In a numerical scheme, the time-dependent Kohn-Sham equation (1) and the equations of motions (3) and (4) are calculated by using the leap-frog method. The electron orbital functions $\psi_{j}(x,t)$, and densities $\rho(x,t)$ and $n(x,t)$, and potentials are represented as real-space grid data. In addition, to make better accuracy, the resolution of grid for core charge density $n(x,t)$ and potential $V_{\text{l}}$ can be set $M_{\text{HR}}$ times of the resolution of other grid data. In the present simulation, $M_{\text{HR}}=2$. Note that the TDDFT simulation is performed without k-point sampling, which means to treat only a gamma point, because the probabilistic interpretation in quantum many-body systems, discussed below, becomes very difficult.

2.2 Preparation of initial state

To perform TDDFT simulation, the initial states of electron orbital functions are prepared by DFT with common Hamiltonian operator (2), generally. However, there are problems particular to ion injection. Although the solution of DFT is a grand state, the initial state that there in an ion far from the target material is not a grand state. Therefore, we cannot simply once do the DFT calculation. In addition, there is a wide vacuum region in the simulation system for PWI, and then the cost of DFT calculation becomes higher. This is not negligible because the cost of the DFT calculation is $O(N^{3})$, which is more expensive than that of a TDDFT calculation of $O(N^{2})$. To solve these problems, in this work, we propose the following approach.

The initial states of the target material and of the incident ion were prepared separately and were combined as shown in Figure 1. Let us now assume that the target material is composed of $N_{x}\times N_{y}\times N_{z}$ supercell and that its surface is parallel to the x-y plane and perpendicular to the z-axis. The lattice vectors of the unit cell are given by $\vec{a}$, $\vec{b}$, and $\vec{c}$. Note that the unit cell here is the smallest unit of the periodic structure, which may be different from the exact unit cell in crystallography. With the width of vacuum region $W_{z}$, the size of the original simulation box is $N_{x}\vec{a}\times N_{y}\vec{b}\times(N_{z}+W_{z})\vec{c}$. In such a setup, the target material has $N_{x}$-fold periodicity and $N_{y}$-fold periodicity in the $x$- and $y$-directions even after adding the vacuum region.

As the procedure 1 in Fig. 1, using this feature, we perform the DFT calculation for a smaller $1\times 1\times N_{z}$ supercell in the small system with the size of $\vec{a}\times\vec{b}\times(N_{z}+W_{z})\vec{c}$. Instead, it is important to set the k-point sampling ranged for $N_{x}\times N_{y}\times 1$. Let put the number of electrons $n_{e}$ in the unit cell, $N_{z}n_{e}$ solutions of electron orbital functions in the small system is obtained by the DFT calculation for each k-point. Then, total number of solutions in the small system for all k-points is $N_{x}N_{y}N_{z}n_{e}$, which is same as the number of solutions required for the original system at a gamma point. By the way, the cost DFT calculation are reduced from $O(N_{x}^{3}N_{y}^{3}N_{z}^{3})$ to $O(N_{x}N_{y}N_{z}^{3})$. In addition, the structural relaxation for the tungsten atom arrangement is also performed in the small system.

The electron orbital functions $\tilde{\psi}_{\vec{k}}(\vec{r})$ for each k-point in the small system can be expanded to the electron orbital function $\psi(\vec{r})$ in the original system according to the Bloch’s theorem[27] as follows:

$$\psi(\vec{r})=e^{i\vec{G}\cdot\vec{r}}\tilde{\psi}_{\vec{k}}(\vec{r}),$$ (8)

where the reciprocal lattice vector $\vec{G}$ is defined by

$$\vec{G}=\frac{2\pi k_{x}}{V_{\text{u}}N_{x}}\vec{b}\times\vec{c}+\frac{2\pi k_{y}}{V_{\text{u}}N_{y}}\vec{c}\times\vec{a},$$ (9)

and $V_{\text{u}}=\vec{a}\cdot(\vec{b}\times\vec{c})$. Note that although $\vec{r}$ for the original system can be out of range for the small system, the electron orbital functions in the small system has periodicity, that is, $\tilde{\psi}_{\vec{k}}(\vec{r})=\tilde{\psi}_{\vec{k}}(\vec{r}+\vec{a})=\tilde{\psi}_{\vec{k}}(\vec{r}+\vec{b})$. Thus, the initial state of the target material for TDDFT simulation, which treat only a gamma point, was prepared.

As the procedure 2 in Fig 1, we prepare the initial state of the incident particle. If the incident particle is just a nucleus such as $\text{He}^{2+}$, it is simply added into the simulation box as a point particle with incident velocity $\vec{v}_{0}$. On the other hand, if the incident particle is an ion accompanied by some electrons such as $\text{He}^{1+}$, we need to prepare the electron orbital functions around the incident nucleus $\psi_{0}(\vec{r})$ by using DFT in the simulation box whose size is in common with the original system. However, the obtained electron orbital function $\psi_{0}(\vec{r})$ rests around the central nucleus. Therefore, $\psi_{0}(\vec{r})$ is converted to have the initial velocity $\vec{v}_{0}$ as follows

$$\psi_{v_{0}}(\vec{r})=e^{i\vec{v}_{0}\cdot\vec{r}}\psi_{0}(\vec{r}),$$ (10)

where the Planck constant $\hbar=1$ in the atomic unit system. The exponential factor $e^{i\vec{v}_{0}\cdot\vec{r}}$ depends on the initial velocity and does not necessarily satisfy the periodic boundary condition in the original simulation box. Therefore, it is necessarily for $\psi_{0}(\vec{r})$ to be sufficiently damped within the size of the original simulation box.

As the procedure 3, the electron orbital functions prepared in procedures 1 and 2 are combined, and the TDDFT simulations were performed. This method can be generally used in TDDFT simulations even for any kinds of incident particles and target materials.

2.3 Simulation condition

The present simulation condition is described here. The target tungsten material is composed of 48 atoms, which was prepared as $3\times 2\times 2$ supercell and was relaxed in atomic positions according to the above method. There are 576 electron orbital functions for the target material side where were the grand state at temperature $T=0$. The half of them are up-spin states and the other half are down-spin states. The lattice constant of the tungsten material is 3.16 Å = 5.97 Bohr. The size of the simulation box is $17.91\times 16.89\times 50.76$ Bohr³. The tungsten material is located at center of the simulation box, and the width of vacuum region is about 33 Bohr in $z$ direction. The tungsten material has the (110) surface. The simulation box follows periodic boundary conditions for all directions. The numbers of real-space grids in the simulation box is $60\times 60\times 180$ grids, of which the resolution corresponds to the cutoff energy of 110 Ry in the case of the DFT with plane wave basis.

The incident particle is $\text{He}^{2+}$ or $\text{He}^{1+}$. In the case of $\text{He}^{1+}$ injection, one electron orbital function is added by the above method where it is the ”1s” state with up-spin. The incident energy is 100 eV, and the incident angle is perpendicular to the $x$-$y$ surface. From the incident energy and angle, the initial velocity $\vec{v}_{0}$ was determined. We tried the cases in terms of the incident target sites (see Fig. 2). The injection onto the site C expects to bring about ”channeling”. Here, the $\text{He}^{2+}$ injections onto the site A, B and C are named the cases (a), (b), and (c), while the $\text{He}^{1+}$ injection onto the site A, B and C are named the cases (d), (e), and (f), respectively.

The time step used for the time-dependent Kohn-Sham equation (1) and the equations of motions (3) and (4) is $\Delta t=4.84\times 10^{-19}$ s. The simulations were terminated when the He nucleus is sufficiently far away from the tungsten material after reflected or penetrating. Actually, simulations required about 60000 to 80000 steps. Note that the simulations was finished about 18 to 24 hours with 320 CPU cores on Intel Xeon Gold system.

3 Results and Discussion

In the present work, the processes of He ion injections for six cases (a)-(f) were simulated. As predicted, the He ion was reflected for the cases (a), (b), (d) and (e), while the He ion penetrated in the cases (c) and (f). Since the He collided head-on with the tungsten atom for the former cases, the velocity of the He after reflected were also parallel to the z-axis. In the latter cases, the velocity of He after penetrating in was also parallel to the z-axis, which indicates the channeling. Understandably, the apparent dynamics of reflection and penetration were in common with $\text{He}^{2+}$ and $\text{He}^{1+}$ injections.

The simulations were terminated when the reflected or penetrated He was sufficiently far away from the surface or back surface of the tungsten material. Actually, at the final time, the distances between the He and the surface in the case of reflection was greater than 15.1 Bohr, and the distances between the He and the back surface in the case of penetration was greater than 9.1 Bohr.

Next, we discuss the dynamics of electrons with respect to the motion of these nuclei. Figure 3 is simulation snapshots of $\text{He}^{2+}$ injection of the case (a). There are no electron density around the He nucleus at the initial time. As the He nucleus reaches the nearest tungsten atom, the electron density are seen moving from the surface side to the near side of the He nucleus. After collision, the He moves away from the surface. The He is still accompanied by the electron density Even after it has been sufficiently far from the tungsten surface.

Similarly, in the $\text{He}^{2+}$ injection for cases (b) and (c), the reflected or penetrating He was accompanied by a surrounding electron density even after they were sufficiently far from the surface or back surface of the tungsten material.

On the other hand, in the cases of $\text{He}^{1+}$ injections for the cases (d)-(f), the change in electron density around He after reflection or penetration was hardly noticeable through visualization. To estimate electronic charge around the He nucleus $Q_{\text{He}}$, the helium peripheral region $V$ is defined as the region centered on the helium atom and within the radius of 6.0 Bohr. The electronic charge $Q_{\text{He}}$ is estimated by the spatial integral of the electron density in that region $V$ at the final time, as shown in Table 1. In fact, $Q_{\text{He}}$ is kept about 1.0 to 1.2 in the case of $\text{He}^{1+}$ injection, while it changed from initial 0 to 1.1 or greater in the case of $\text{He}^{2+}$ injection.

It is not possible to judge from the electron density whether the He ions have been neutralized or not. This is because the electron density around He is actually contributed by multiple electrons. Therefore, according to many body quantum system, we consider the probability of detecting electrons around the He nucleus. If no electron is detected around the He nucleus, the He should correspond to be observed as a $\text{He}^{2+}$. Similarly, if only one electron is detected, the He corresponds to be observed as a $\text{He}^{1+}$, and if two electrons are detected, the He corresponds to be observed as a neutral $\text{He}^{0}$. In addition, the spin states of electrons must be taken into account.

Note that, since the simulation follows the time evolution without real observation, the final state of the simulation is obtained as a superposition of those $\text{He}^{2+}$, $\text{He}^{1+}$, and $\text{He}^{0}$ states.

All-electron wave function is given as the Slater determinant of electron orbitals $\psi_{j}(x,t)$. In this simulation, the electron orbitals for $j=1,\cdots,N_{a}$ are the states corresponding to up-spin state, and those for $j=N_{a}+1,\cdots,N$ are the states corresponding to down-spin state. The probabilities of detecting the electrons of up-spin state and down-spin state are estimated independently. In our other work[28], using the helium peripheral region $V$ again, the probability of detecting $m$ electrons of the up-spin state in $V$ is derived by

$$P(m;S_{\uparrow}) =\eta(m;S_{\uparrow})-(m+1)\eta({m+1};S_{\uparrow})$$

$$\quad+\frac{(m+2)(m+1)}{2}\eta({m+2};S_{\uparrow})+O(\braket{\psi_{k_{j}}}{\psi_{\tau(k_{j})}}_{V}^{3}),$$ (11)

$$\eta(m;S_{\uparrow}) =\sum^{S_{\uparrow}}_{\{j_{1},\cdots,j_{m}\}}\sum_{\tau}\text{sgn}(\tau)\prod^{m}_{k=1}\braket{\psi_{j_{k}}}{\psi_{\tau(j_{k})}}_{V},$$ (12)

$$\braket{\psi_{k_{j}}}{\psi_{\tau(k_{j})}}_{V} =\int_{V}\psi^{*}_{k_{j}}(x,t)\psi_{\tau(k_{j})}(x,t)dx,$$ (13)

where the summation $\sum^{S_{\uparrow}}_{j_{1},\cdots,j_{n}}$ means the sum of the all subsets of $n$-combination of the set $S_{\uparrow}=\{1,\cdots,N_{a}\}$, the summation $\sum_{\tau}$ means the all permutations from the subset $\{j_{1},\cdots,j_{n}\}$ to $\{\tau(j_{1}),\cdots,\tau(k_{n})\}$, and the $\text{sgn}(\tau)$ is positive or negative sign when $\tau$ is even or odd permutations. Similarly, the probability of $m$ electrons of down-spin state in $V$ is given by $P(m;S_{\downarrow})$, where the set $S_{\downarrow}=\{N_{a}+1,\cdots,N\}$.

Furthermore, the observation probability of $\text{He}^{2+}$, $P(0)$, the observation probability of $\text{He}^{1+}$, $P(1)$, and the observation probability of $\text{He}^{0}$, $P(2)$, after the incidence or penetration are evaluated by

$$P(0) =P(0;S_{\uparrow})P(0;S_{\downarrow}),$$ (14)

$$P(1) =P(1;S_{\uparrow})P(0;S_{\downarrow})+P(0;S_{\uparrow})P(1;S_{\downarrow}),$$ (15)

$$P(2) =P(2;S_{\uparrow})P(0;S_{\downarrow})+P(1;S_{\uparrow})P(1;S_{\downarrow})+P(0;S_{\uparrow})P(2;S_{\downarrow}).$$ (16)

The evaluated observation probabilities $P(m)$ are shown in Table 1. In the case of the $\text{He}^{2+}$ injection, the probabilities of observing $\text{He}^{1+}$ after reflection and penetration was 0.405, 0.385, and 0.376 for the cases (a), (b), and (c), respectively. The probabilities of observing neutral $\text{He}^{0}$ were 0.135, 0.129, and 0.125. The dependence on the incident position was small.

On the other hand, in the cases (d)-(f) of the $\text{He}^{1+}$ injection, the probabilities that $\text{He}^{1+}$ is maintained after reflection and penetration are 0.93 to 0.98. The reason why this $\text{He}^{1+}$ hardly changes to $\text{He}^{0}$ can be explained by the spin-state limitation of the simulation method as follows: in this simulation, the spin-orbit interaction is neglected in the Hamiltonian operator in Eq. (2). That is, not only total energy and momentum but also total spin magnetic moment are conserved in the time evolution. Here, the initial tungsten material has an even number of electrons and its spin magnetic momentum is initially zero (spin-neutral). Adding $\text{He}^{1+}$ with one up-spin electron in the simulation box, the total spin magnetic moment of the system becomes +1. If the injected $\text{He}^{1+}$ changes to $\text{He}^{0}$, it captures a down-spin electron from the tungsten material. Since the total spin magnetic moment is conserved, the tungsten material after collision should be a spin magnetic moment of $-1$. However, the spin-polarized state is unstable in energy for such a small tungsten material. Therefore, we can conclude that the change from $\text{He}^{1+}$ to $\text{He}^{0}$ hardly occurred in the present simulation system.

If the size of the target material is sufficiently large, the effect of spin polarization becomes relatively small. Therefore, the problem caused by total spin magnetic moment is a simulation-specific size effect issue. It is also not an issue in experiments because the material is sufficiently large and there are inflow and outflow of electrons from the ground connection. In addition, almost the previous theoretical works[7, 8, 9, 10] did not consider the change of electronic state and energy in the material side.

From Fig. 3(g), the electron density around the reflected the He nucleus has a butterfly-like shape rather than a sphere. This implies a 2p state of atomic orbital. Then, we examine whether the electronic state of He after reflection or penetration is a ground state or an excited state. Note that we focus on the case of $\text{He}^{2+}$ injection, because significant electronic transition is not recognized in the $\text{He}^{1+}$ injection.

We consider analyzing the electron orbital functions around the He nucleus and to divide into atomic orbitals. However, the analysis of taking the inner product with the atomic orbital of an independent He atom is not adequate because the function form of the atomic orbital depends on the number of electrons around the nucleus. For example, the 2p orbital of $\text{He}^{0}$ is different from the 2p orbital of $\text{He}^{1+}$. Therefore, we now take the inner product only with the spherical harmonic function as follows.

To begin with, a spherical coordinate system $(r,\theta,\phi)$ centered at the position of He nucleus $\vec{R}$ at the final time is introduced. That is, $\vec{r}-\vec{R}$ $=$ $(r\sin\theta\cos\phi$, $r\sin\theta\sin\phi$, $r\cos\theta)$. In the spherical coordinate system, the inner product between the electron orbital function and the spherical harmonic function were estimated, where only the s and p orbitals were considered for simplicity. The spherical harmonic function of the s orbital is given by $Y_{\text{s}}(\theta,\phi)=Y_{00}(\theta,\phi)=1/\sqrt{4\pi}$, and the spherical harmonic functions of the $\text{p}_{x}$, $\text{p}_{y}$, and $\text{p}_{z}$ are given by $Y_{\text{px}}(\theta,\phi)=(1/\sqrt{2})[-Y_{11}(\theta,\phi)+Y_{1,-1}(\theta,\phi)]=\sqrt{3/4\pi}\sin\theta\cos\phi$, $Y_{\text{py}}(\theta,\phi)=(i/\sqrt{2})[Y_{11}(\theta,\phi)+Y_{1,-1}(\theta,\phi)]=\sqrt{3/4\pi}\sin\theta\sin\phi$, and $Y_{\text{pz}}(\theta,\phi)=Y_{10}(\theta,\phi)=\sqrt{3/4\pi}\cos\theta$, respectively. By taking the inner product of the $n$-th electron orbital function $\psi_{n}(x,y,z)$ and the spherical harmonic function $Y_{X}(\theta,\phi)$, the radial function $\rho_{X}(r)$ is defined as follows:

$$\braket{Y_{X}}{\psi_{n}} =\int Y_{X}(\theta,\phi)\psi(\vec{r})\sin\theta d\theta d\phi,$$ (17)

$$\rho_{X}(r) =\sum_{n}\braket{Y_{X}}{\psi_{n}}^{2},$$ (18)

where $X$ is s, $\text{p}_{x}$, $\text{p}_{y}$, or $\text{p}_{z}$. Then, integrating $\rho_{X}(r)$ in the radial direction, the charge $Q_{X}$ associated with each s, $\text{p}_{x}$, $\text{p}_{y}$, or $\text{p}_{z}$ orbital is obtained.

$$Q_{X}=\int^{\infty}_{0}r^{2}\rho_{X}(r)dr.$$ (19)

The charges $Q_{X}$ for cases (a), (b) and (c) are shown in Fig. 4. Note that the sum of $Q_{X}$ is not 1 because $Q_{X}$ has a dimension of the charge and is not a normalized probability. In the $\text{He}^{2+}$ injection of the case (a) corresponding to Fig. 3(g), the charge $Q_{\text{p}_{y}}$ is 0.68, while the charge $Q_{\text{s}}$ is 0.26. Therefore, the $\text{p}_{y}$ orbital is dominant in this case. On the other hand, in the injections for case (b) and (c), the s orbitals account for about half of the total charge. Thus, the ratio of s- and p-orbitals in the electronic state around the He nucleus differs depending on incident positions.

Furthermore, the p orbitals are clearly excited states for He. For the s orbital, we can distinguish between the ground state and excited states from the shape of the radial function $\rho_{X}(r)$. To be precise, $r^{2}\rho_{X}(r)$ was plotted in Fig. 4. If the orbital is the 1s orbital, as the radial coordinate $r$ increases, $r^{2}\rho_{\text{s}}(r)$ increases from 0 to a peak and then decreases monotonically. On the other hand, if the orbital is the 2s orbital, $r^{2}\rho_{\text{s}}(r)$ has two peaks and the one node which is a point satisfying $r^{2}\rho_{\text{s}}=0$ without the point that $r=0$. In general, $r^{2}\rho(r)$ of $n$-th s-orbital should have $n$ peaks and $n-1$ nodes. From the figure for all cases, $r^{2}\rho_{\text{s}}(r)$ has two peaks and then the 2s orbital is dominant. The reason why $r^{2}\rho_{\text{s}}(r)$ is not completely zero at the nodes, especially in the case (a), is because the component of 1s orbital is a little mixed. Consequently, the 2s and 2p orbitals are dominant in the electron orbital functions around reflected or penetrating He nucleus, and they are excited states.

One of the issues in historical researches is whether the electron transition mechanism from the surface to He nucleus is the Auger process or the resonance process. For He injection, the previous researches expected as follows: In the Auger process, the electron is emitted through a transition to the grand state of the 1s orbital. On the other hand, in the resonance process, the transition to the excited 2s and 2p orbital occurs. The fact that transitions to the excited 2s and 2p states were dominant in the present simulation suggests that the electron transition mechanism of He injection was a resonance process of the previous classifications. Certainly, the electron emission predicted by the Auger process was not occurred in the simulation. However, it may not be necessary to divide Auger and resonance processes[9]. We can now treat the continuous wave function and the motion of nuclei and study the position dependence of the surface at atomic resolution. Accumulation of the TDDFT simulations is expected to lead to further discussions.

In addition, we expect that by using the present TDDFT simulation, the electron stopping power will be also investigated. For example, in the BCA simulation, the kinetic energy loss of the projectile is composed of the kinetic energy transfer to target atoms and the energy to excite electrons. The latter loss is the electron stopping power. These, along with the threshold energy to eject atoms from the material, are the three major components of BCA. The kinetic energy transfer is determined by potential model between atoms, and it can be evaluated from quantum mechanics. Actually, we have derived the ReGZ potential[29] analytically, which is high order model of the ZBL potential[30]. The threshold energy to eject atoms can be evaluated by using DFT calculation on recent computer. However, the modeling of energy stopping power was difficult. Here, the TDDFT simulation can treat the excitation of electrons during the nuclear motion, and it is helpful for the modeling the electron stopping power. For instance, We estimated the kinetic energy of the He nuclei at final time (see in Table 1). In particular, both the cases (a) and (b) brought about the head-on collision, and then the kinetic energy transfer were commonly 8.3 eV, which is determined by masses of He and W atoms. Therefore, the difference of the kinetic energy between the cases (a) and (b) are affected by the electron stopping power. Since the kinetic energy in the case (b) is smaller than that in the case (a), the loss due to the electron stopping power in the case (b) was larger than that in the case (a). This is because the He nucleus reflected at the second layer in the case (b) is more deeply submerged in the free electrons in the metal than the case (a).

Finally, we discuss some issues to apply the TDDFT into PWI researches that we noticed in the present work. First, we should remember that the TDDFT is a mean-field approximation. In short, the trajectories of nuclei may be different between the neutralized He and the ionized He, while in this simulation, they move on the same trajectory. Second, the representativity of the excited state depends on the pseudo-potential, in addition to the issue of DFT applicability limit to grand state. Third, the numerical accuracy is a challenge for the PWI problems typically because the incident energy of 100 eV is much higher than general chemistry problems. We honestly mention that the numerical error of order $\pm 1$ eV in the current version of the QUMASUN code in the present simulation. However, there is room for further improvement, which will be reported in the future.

4 Conclusion

The He ion injection onto the tungsten surface was investigated using TDDFT simulation. We developed the TDDFT code QUMASUN for this work and simulated the process of electron transfer from the surface to the He nucleus. In addition, we proposed the method to prepare the initial state of TDDFT simulation for the target system in PWI.

Interestingly, the simulation results are a superposition of the states of $\text{He}^{2+}$, $\text{He}^{1+}$, and $\text{He}^{0}$. This indicates that the neutralization of incident ions cannot be determined only by electron charge. Instead, the probability of detecting electrons should be considered using the electron wave function. As a result, the probabilities that the injected $\text{He}^{2+}$ transforming into $\text{He}^{1+}$ and $\text{He}^{0}$ on the surface are approximately 40 percent and 13 percent, respectively. The captured electrons by $\text{He}^{1+}$ and $\text{He}^{0}$ predominantly occupy the 2s and 2p orbitals, corresponding to excited states.

Moreover, we opened the issues to apply TDDFT simulation into PWI researches. For example, the problem in terms of the spin-polarization of material was demonstrated in the case of the $\text{He}^{1+}$ injection as a small size effect. Additionally, more simulations are needed to compare with experimental results, considering variables such as position, incident angle, and energy.

Despite these challenges, TDDFT simulations offer unique insights into the neutralization process, providing a level of detail that methods like MD and BCA cannot achieve. This is particularly valuable for fusion science and the study of plasma-wall interactions (PWI), where understanding the quantum mechanical aspects of these processes is crucial. Future work will focus on refining these simulations and exploring a wider range of parameters to further elucidate the neutralization process.

Case	Projectile	Target site	$Q_{\text{He}}$	$P(0)$	$P(1)$	$P(2)$	$K_{\text{He}}$ [eV]
(a)	He2+	A	1.11	0.441	0.405	0.135	89.7
(b)	He2+	B	1.16	0.464	0.385	0.129	80.1
(c)	He2+	C	1.29	0.478	0.376	0.125	81.8
(d)	He1+	A	1.03	0.017	0.980	0.004	92.7
(e)	He1+	B	1.16	0.064	0.930	0.004	87.0
(f)	He1+	C	1.18	0.032	0.967	0.001	91.6

Table 1:

Figure 1:

Figure 2:

Figure 3:

Figure 4: