Numerical Interpretation of Hydrogen Thermal Desorption Spectra for Iron with Hydrogen-Enhanced Strain-Induced Vacancies

We obtained thermal desorption spectra of hydrogen for a small-size iron specimen to which strain was applied during charging with hydrogen atoms. In the spectra, a shoulder-shaped peak in the high-temperature side was enhanced compared with the spectra of the specimen to which only strain was applied. We also observed that the peak almost disappeared by aging processes at ≥ 373 K. Then, assuming that the shoulder-shaped peak results from hydrogen atoms released by vacancies, we simulated the thermal desorption spectra using a model incorporating the behavior of vacancies and vacancy clusters. The model considered up to vacancy cluster V9\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{V_9}}$$\end{document}, which is composed of nine vacancies, and employed the parameters based on atomistic calculations, including the H trapping energy of vacancies and vacancy clusters that we estimated using the molecular static calculation. As a result, we revealed that the model could, on the whole, reproduce the experimental spectra, except two characteristic differences, and also the dependence of the spectra on the aging temperature. By examining the cause of the differences, the possibilities that the diffusion of clusters of V2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${V_2}$$\end{document} and V3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${V_3}$$\end{document} is slower than the model and that vacancy clusters are generated by applying strain and H charging concurrently were indicated.


I. INTRODUCTION
HYDROGEN embrittlement (HE) is known as one cause of delayed fracture of high-strength steels and cold cracking of weld alloys. [1][2][3][4] As is widely known, in HE, H atoms gather at the location of stress concentration resulting from tensile stress or residual stress and lead to toughness degradation in materials by promoting the production of defects and/or by weakening the strength of atomic bonds, and then the toughness degradation causes cracks or fractures. [5,6] In recent years, it has been reported that defects are formed in tempered martensitic steels to which elastic stress is applied during charging with H atoms (H charging); this is referred to as H-enhanced lattice defect formation [7] . From the measurement of the positron annihilation lifetime (PAL), it is reported that the formed defects are considered as vacancies. The report [8] also mentions the effect of H atoms on the vacancy formation in deformed steel based on PAL measurement. Similar vacancy formation is also observed in tempered martensitic steels prestressed cyclically during H charging [9] , in Inconel 625, iron [10] , and nickel alloys [11] , which are deformed and H charged concurrently. In the report [11] , PAL is also used for estimating vacancy formation. Because ductility loss is observed in a tensile test of the steels and of the alloys even after H atoms charged for the vacancy formation are removed, as shown in Figure 1, it is implied that the vacancies themselves bring about HE, and it is mentioned that the phenomenon may support the H-enhanced strain-induced vacancy (HESIV) model that is proposed as a mechanism of HE [12][13][14] . However, how many vacancies bring about HE is still unclear, and the amount and size of vacancy clusters that can be generated by vacancy agglomeration relating to HE are also unknown.
Vacancy formation can also be confirmed by the change of thermal desorption spectra measured by thermal desorption spectrometry (TDS), which is a widely employed powerful experimental method for estimating H-trapping states, in which H atoms are trapped by defects such as dislocations, vacancies, and grain boundaries [15,16] . By deforming a specimen during H charging in the TDS spectrum, a new desorption peak appears partly overlapping a peak of the defects, which are regarded mainly as dislocations. Then, by annealing the specimen at 473 K, the newly appeared peak disappears, while the peak of dislocations remains. Thus, the peak is formed by H atoms released from vacancies [13] . In the reports [9,10] , the vacancy formation is also inferred from the change of TDS spectra. Hence, by numerically simulating the thermal desorption spectra of a specimen strained during H charging considering the behavior of vacancies and vacancy clusters, the quantitative estimation of their formation becomes possible. In addition, such estimation can promote quantitative discussion about the relation of vacancies and vacancy clusters to HE.
Understanding and interpretation of experimental thermal desorption spectra using a numerical model have been reported [17][18][19][20][21][22] . Although various kinds of models are employed in those reports, the models can broadly be classified mainly into three types [23] . Among them, the model based on the McNabb-Foster theory [24] can be applied in principle to the simulation of desorption spectra without restriction because it includes all processes of diffusion, trapping, and detrapping. In addition, since these processes are clearly separated in the model that is described by reaction-diffusion equations, it is considered that the behavior of vacancies and vacancy clusters can be incorporated comparatively easily into the model. Thus, the model based on the McNabb-Foster theory is appropriate for simulating thermal desorption spectra including the peak for vacancies and vacancy clusters.
The experimental thermal desorption spectra of the tempered martensitic steels and alloys with vacancy formation have the peak for vacancies partly overlapping the peak for the other defects [9,10,13] , as described above. In more detail, the H desorption resulting from vacancies merely broadens the high-temperature side of the peak for the other defects because of using comparatively large-size specimens in which the influence of H diffusion becomes large. Therefore, it is difficult to separate and identify the peak for vacancies. However, by measuring H desorption from a small-size specimen using a low-temperature TDS (L-TDS) apparatus [25,26] that can heat a specimen from a low temperature of 77 K to 473 K, thermal desorption spectra in which the peak for vacancies is separated from the peak for other defects can be obtained [27] . From the viewpoint of modeling thermal desorption spectra including a desorption peak for vacancies and vacancy clusters, the structure in martensitic steels is so complicated that it is difficult to model the influence of the structure on the diffusion of H atoms, vacancies, and vacancy clusters because of many unknown factors. Meanwhile, it is considered that the modeling of thermal desorption spectra from an iron specimen is comparatively easy since its structure is simpler than that of steels.
Therefore, in the current study, we carried out L-TDS measurements using a small-size specimen with 0.3 mm thickness cut from an iron specimen in which strain was applied and charged with H atoms concurrently. Then, we simulated the experimental spectra using a numerical model incorporating the model for the migration of vacancies and vacancy clusters and the growth and dissociation of vacancy clusters [28] . In the simulation, we employed the H trapping energy of vacancies and vacancy clusters, which was evaluated by the molecular static calculation using an embedded atom method (EAM) potential. We also applied the model to the simulation of the variation of the amount of H atoms, vacancies, and vacancy clusters in several processes before the L-TDS measurement, such as H degassing, annealing, polishing, and H charging processes.

A. Material and Specimen
We used pure iron; Table I shows its chemical composition. It was annealed for 1800 seconds at 1173 K to reduce defects. Figure 2 schematically shows the shape of the specimens for the tensile test and L-TDS measurement.

B. Procedure
The procedure for the experiment for TDS measurement is illustrated in Figure 3.
To introduce vacancies, after pre-charging with H atoms to the specimen for 1800 seconds, we applied the tensile strain to the specimen (Figure 2(a)) at 8:3 Â 10 À6 s À1 strain rate until the total tensile strain  became 0.25 and was charged with H atoms by the cathodic electrolysis method concurrently. In the cathodic electrolysis charging, we applied 50 A m À2 in the electric current to the specimen during immersing in a 0.1 N NaOH þ NH 4 SCN 5 g L À1 solution at 303 K. After cutting out a square plate specimen from the tensile test specimen, we removed H atoms from the cut-out specimen by holding it at 303 K for 3600 s. In addition, we mechanically and chemically polished the square plate specimen at room temperature to obtain the specimen shown in Figure 2(b). It took about 5400 s for those processes. We labeled this specimen as the H + strain specimen. For comparison, we also prepared another tensile test specimen to which strain was applied without H pre-charging and H charging and then from that cut out the square specimen, and the specimen shown in Figure 2(b) was obtained by the same processes as those for the H + strain specimen. Although this specimen was not H charged, the annealing process at 303 K for 3600 s was also applied for comparing with the case of the H + strain specimen. The polishing process was the same as for the H + strain specimen. We labeled this specimen as the strain specimen. By immersing both specimens in a 0:1 N; NaOH þ NH 4 SCN 20 g L À1 solution at 303 K for 3600 s, the specimens were H charged as a tracer for L-TDS measurement. After the specimens were held in the liquid nitrogen, we measured the thermal desorption spectrum for each specimen using the L-TDS apparatus mainly with 60 K h À1 as the heating rate.

A. Model
We employed a model based on the McNabb-Foster theory, which is represented by the following reaction-diffusion equation: Here, C H L represents the concentration of H atoms in the normal lattice site (the interstitial site), and C H V and C H 1 represent the concentration of H atoms trapped by vacancies and vacancy clusters and by defects except for them, respectively. The LHS is their time evolution term, and the RHS is the diffusion term. The diffusion coefficient is given by   where the first and second terms in the last line mean the trapping and detrapping processes, respectively, and KðX; TÞ ¼ expðÀX=RTÞ. In these terms, E 1 , N 1 , and h 1 are the H trapping energy, trap site concentration, and H occupancy of the trap sites, respectively, for defects except the vacancies and vacancy clusters. The pre-exponential factors for K(X, T) in the two terms are represented by k 01 and p 01 , respectively. While Eq. [2] is the usual form used so far, we employed the following equation where i is the size of vacancy clusters, namely the number of vacancies composing one vacancy cluster. Vacancies are represented by i ¼ 1: The meaning of E V i ; N V i ; k 0V i ; and p 0V i is the same as in Eq. [2]. Since the time evolution of vacancies and vacancy clusters brings about the change of N V i , we incorporated the effect to the model by adding the term P i @N V i @t h V i to Eq. [3]. Based on the reports [28,29] , we employed the following equations for estimating The equations represent the time evolution of the concentration C V i of vacancies and vacancy clusters whose maximum size is n. The first and the second terms on the RHS of these equations represent the diffusion process and the annihilation process at sink sites such as dislocations, respectively. The terms with the coefficients b þ V 1 þV i and b À V i represent the agglomeration and dissociation processes, respectively. The diffusion coefficient is given by is the sink strength, and C eq V i represents the thermal equilibrium concentration. Table II shows D 0V i and Q V i until i ¼ 6; and we assumed that vacancy clusters > i ¼ 6 are immobile. We used C eq V 1 ¼ expðÀ2:045 Â 10 5 Jmol À1 = RTÞ [29] and assumed C eq V i ¼ 0 for i>1: The reaction constants are given by b þ [28] . N s is the number of available sites, X is the atomic volume, r 0 ¼ 3:3Å , and E b V i is the binding energy of the vacancy cluster V i to a vacancy. In the model, the growth and dissociation of vacancy clusters are represented by the absorption and emission of a single vacancy, respectively. According to the report [30] , the value of E b V i is evaluated up to V 4 by first-principles calculation, and the value for clusters larger than V 4 is given by extrapolation. Table II shows E b V i until i ¼ 6; and for V 5 and V 6 ; the value obtained by extrapolation is shown. We used N s ¼ N L ði À 1Þ=i for i>1 where N L ¼ 5:1 Â 10 29 m À3 is the normal lattice site concentration. As for S V i ; we employed 2:5 Â 10 14 m À2 for all i. The trap site concentration is given by According to the report [32] , when a vacancy traps H atoms, its barrier energy of diffusion increases, so that the jump frequency of the vacancy is several orders of magnitude smaller than the case without trapping H atoms. The similar situation may occur for vacancy clusters. Thus, we assumed that vacancies and vacancy clusters with h V i ! 1 are immobile. In addition, since h V i takes a continuous value in the model, the number of H atoms trapped by vacancies and vacancy clusters, h V i m V i can become a positive real number < 1. In such a case, we regarded the value of h V i m V i as the ratio of vacancies or vacancy clusters that trap a H atom to all vacancies or vacancy clusters. Hence, we assumed that vacancies or vacancy clusters of We evaluated the H trapping energy E V i for vacancy clusters by the molecular static simulation with the EAM potential proposed by Wen et al. [33] . Table III summarizes the value of E V i until V 9 . We ignored trap sites whose H trap energy is < 10:0 kJ mol À1 . According to the table, we can see that m V i increases as the cluster size increases and E V i depends on the number of trapped H atoms. As previously described, the number of trapped H atoms, h V i m V i becomes a positive real number. Thus, we interpolated linearly between two H trapping energies, E V i ðs þ 1Þ and E V i ðsÞ, as follows: where s is a positive integer for the number of trapped H atoms.
We determined the other parameters such as E 1 ; N 1 ; N V i ; k 01 ; p 01 ; k 0V i ; and p 0V i from the experimental thermal desorption spectra shown in Section IV using the procedure described in Section V.

B. Method
We solved Eqs. [1] though [4] simultaneously using the finite difference method in a one-dimensional calculation region. Since the specimen shape is a thin square plate, most of the H atoms are desorbed from the widest surface after diffusing in the thickness direction in which the length is shortest. Thus, the one-dimensional calculation region is a valid approximation. This calculation region is adapted by several reports [34][35][36] . Due to symmetry, we used the half-thickness of the specimen for the calculation region. Also, we chose the initial and boundary conditions reflecting the experimental condition.
First, we applied the model in the previous section to the simulation of the degassing process for the H + strain specimen or the annealing process for the strain specimen, and then we applied it to the simulation of behavior of H atoms, vacancies, and vacancy clusters in the polishing process for both specimens. In the degassing process, as its initial state, we supposed that H atoms and vacancies distribute homogeneously and there is no vacancy cluster. Thus, in the simulation, we set the initial values for the H concentration and vacancies in the one-dimensional calculation region of 2:5 Â 10 À4 m. Furthermore, we assumed the equilibrium state of H atoms between the normal lattice site and trap sites. In the simulation for the annealing process, we set the initial value of H atoms to zero. We performed the simulation isothermally at 303 K for 3600 s under the boundary condition that the concentration of H atoms, vacancies, and vacancy clusters becomes zero on the surface. From the degassing or the annealing simulation, we obtained the distributions of H atoms, vacancies, and vacancy clusters and used them as the initial state of the simulation for the polishing process. Strictly speaking, in the polishing process, the specimen thickness decreases; however, we employed the same one-dimensional calculation region with the degassing or the annealing simulation because the influence of the thickness decrease is minute. We simulated the process at 298 K for 5400 s under the same boundary condition as the case of the degassing or the annealing simulation. Next, we simulated the process for charging the specimen of 3:0 Â 10 À4 m with tracer H atoms. In the simulation, we set a flux J at the surface boundary of the calculation region with length 1:5 Â 10 À4 m, and we held the temperature at 303 K. Then, using the states of H atoms, vacancies, and vacancy clusters, which were  obtained by the H charging simulation, we simulated the heating process of the L-TDS measurement at the heating rate of 60 K h À1 . During the simulation, we counted the amount of desorbed H at a constant interval. Section V describes the determination of the initial values of vacancies in the degassing simulation and of vacancies and H atoms in the annealing simulation and the flux J in the charging simulation. Figure 4 shows the experimental thermal desorption spectra for the strain and H + strain specimens.

IV. EXPERIMENTAL RESULTS
In the figure, the spectra of both specimens have one large peak at about 281 K and a shoulder-shaped peak at a higher temperature. We label these peaks as Peak 1 and Peak 2, respectively. The shape of the spectra is similar to that in the case of the tempered martensitic steel [27] . From the comparison between the strain and H + strain specimens, it is seen that the defects corresponding to Peak 2 are increased by H charging. By fitting a Gaussian curve to each peak of both spectra, we estimated the amount of the desorbed H atoms of each peak, which is shown in Table IV; also, the table shows the total desorbed H atoms. In addition, by measuring spectra with other different heating rates, that is, 30 K h À1 and 120 K h À1 , and then identifying the peak temperature from the Gaussian fitting, we evaluated the detrapping activation energy for each peak of the H + strain specimen from the plot of Choo and Lee [18] shown in Figure 5. We obtained 30:5 kJmol À1 for Peak 1 and 54:9 kJ mol À1 for Peak 2. Subtracting the barrier energy of H diffusion, Q H from these values, the former is roughly close to the H trapping energy of screw dislocations [15,37] , and the latter is close to the energy that vacancies trap H atoms.
Furthermore, to examine the characteristic of Peak 2, we measured the thermal desorption spectra of the H + strain specimen in the case that the degassing process (i.e., aging at 303 K for 1 h) aiming at the removal of H atoms was replaced by the aging process aiming at the annihilation of vacancies. We employed three aging temperatures, 323 K, 373 K, and 423 K, and the aging time was 3600 seconds for the three. Figure 6 shows the results. According to the figure, the H desorption rate of Peak 2 decreases with the aging process, and it almost disappears at ‡ 373 K in the aging temperature. From the results, it is considered that Peak 2 is formed by H atoms released from vacancies and/or vacancy clusters. Peak 1 is also influenced by the aging temperature; in the case with ‡ 323 K, it is revealed that the H desorption rate is almost independent of the aging temperature and the peak temperature is influenced. The defect for Peak 1 is considered in Section V-B.

A. Parameter Determination
First, we determined the simulation parameters that are not already set in Section III from the experimental results shown in the previous section. As the H trapping energy of Peak 1, we used E 1 ¼ 26:6 kJ mol À1 , which is obtained from the difference between the evaluated detrapping activation energy of Peak 1 of the H + strain specimen and the barrier energy of H diffusion, Q H . We chose p 01 ¼ 4:6 Â 10 2 s À1 so that the peak temperature of Peak 1 becomes that of the experimental spectrum and k 01 ¼ 9:0 Â 10 À28 m 3 s À1 from the relation k 01 ¼ p 01 =N L : As for N 1 , according to our several preliminary experiments, the amount of H desorption does not almost change even if the H charging was severer than the current condition. Thus, we can suppose that the trap site of Peak 1 with H trapping energy of E 1 is almost in saturation after charging with tracer H atoms. Hence, from the desorbed H atoms of Peak 1 in Table IV, we assumed that h 1 ¼ 0:95 after the charging of tracer H atoms and estimated N 1 to be 1:1 Â 10 24 m À3 and 1:2 Â 10 24 m À3 for the strain and H + strain specimens. In addition, we chose the flux J for the H charging simulation so that h 1 ¼ 0:95 is satisfied under the experimental condition of charging the specimens with tracer H atoms.
Meanwhile, we determined p 0V i and k 0V i as follows. Although the existence of vacancy clusters can be supposed after charging with tracer H atoms, we considered only vacancies and assumed that h V 1 ¼ 0:3 and then calculated the value of C H;a:c: L is the H concentration in the normal lattice site after charging. The calculated value corresponds to p 0V 1 =k 0V 1 . Even though p 0V 1 =k 0V 1 becomes N L when the H trap energy is a constant, we can suppose that it does not have to be N L since E V 1 depends on h V 1 in the case of vacancies. Thus, we chose p 0V 1 ¼ 5:7 Â 10 5 s À1 so that the peak temperature of Peak 2 approaches the experimental result and then k 0V 1 ¼ 2:7 Â 10 À30 m 3 s À1 is obtained using the calculated value above. We also used the same values with p 0V 1 and k 0V 1 for all the other vacancy clusters. As for the initial value of the concentration of vacancies C init V 1 for the degassing simulation, we chose the value so that the total amount of H atoms trapped by vacancies and vacancy clusters before the heating simulation becomes the experimental value of Peak 2 in Table IV. This determination is possible because the amount of H atoms trapped by defects depends on the amount of the defects. In addition, the size distribution of vacancy clusters before the L-TDS measurement depends on the migration behavior of vacancies in the pre-measurement processes such as degassing, polishing, and charging, and the migration of vacancies is suppressed by their trapped H atoms. Therefore, the initial concentration of H atoms trapped by vacancies in the degassing process influences the thermal desorption spectra. In the degassing simulation, by varying the total H concentration in the calculation region in equilibrium calculation, we determined the concentration of H atoms trapped by vacancies that makes the shape of simulated desorption spectra close to the experimental spectra as the initial value C H;init V 1 . We also determined the initial concentration of H atoms trapped by defects for Peak 1 C H;init 1 from the equilibrium calculation simultaneously.

B. Simulation Results
Herein, we show the simulation results for the pre-measurement process. Figure 7 shows the change of the amount of H atoms trapped by each kind of trap site and of the concentration of vacancies and vacancy clusters in the strain specimen during the annealing, polishing, and charging processes. Figure 8 shows the result for the H + strain specimen where the degassing process was simulated instead of the annealing process.
The initial values of the simulation for the annealing process or the degassing process are as follows:   the strain specimen, and C H;init 1 ¼ 0:25 ppm (0.97), C H;init V 1 ¼ 0:043 ppm (0.34), and C init V 1 ¼ 1:2 Â 10 À6 at% for the H + strain specimen where the value in the parentheses means the occupancy of H atoms for the corresponding trap site. These values were chosen so that the simulated spectra could reproduce the experimental spectra in Figure 4. As for the strain specimen, since no H atoms suppressing vacancy migration are contained before the charging process, as seen in Figure 7(a), vacancy clusters are generated from the beginning of the annealing process because of vacancy migration, as seen in Figure 7(b). However, it is seen that the generation of vacancy clusters is not so large that the decrease in vacancy (V 1 ) concentration is small. In the charging process, the concentration of vacancies and vacancy clusters is almost constant because their migration is suppressed by H atoms. In the case of the H + specimen shown in Figure 8, vacancy clusters are not generated from the beginning to about 0.8 h in the degassing process because of the suppression of vacancy migration by H atoms. After that, since vacancies start to migrate until the end of the polishing process because of the decrease of H atoms, the V 1 concentration decreases slightly, and vacancy clusters are generated rapidly. Compared with the case of the strain specimen, since the initial V 1 concentration is larger, the generation of vacancy clusters is also larger. The tendency in the charging process is similar to that of the strain specimen. Figure 9 shows the amount of H atoms trapped by vacancies, vacancy clusters, and trap sites of Peak 1 and the concentration of vacancies and vacancy clusters. These values are ones before the heating process of TDS. According to the figure, H atoms trapped by V 1 are the largest after Peak 1, and V 5 traps the largest amount of H atoms among vacancy clusters. The tendency of trapped H atoms reflects the tendency of the concentration of vacancies and vacancy clusters. Thus, it means that the amount of trapped H atoms depends strongly on the concentration of vacancies and vacancy clusters rather than the trapping energy. It is also found that the H atoms trapped by V 2 and V 3 are smaller than those for V 4 and V 5 in the strain specimen and those for V 4 to V 6 in the H + strain specimen. This is because the concentration of the former is smaller than that of the latter. The reason why the concentration of V 2 and V 3 is smaller is that their diffusion is faster than that of the vacancies and other clusters. The tendency is almost the same between both specimens though the trapped H atoms and the concentration of  vacancies and vacancy clusters of the strain specimen are smaller than those of the H strain specimen. The total amount of H atoms trapped by vacancies and vacancy clusters is 0.006 and 0.037 ppm for the strain and H + strain specimens, respectively, and the amount of trapped H atoms for Peak 1 is 0.22 and 0.24 ppm for the strain and H + strain specimens, respectively. Note that these values are the same as in the experimental value in Table IV, as intended. Figure 10 shows the thermal desorption spectra simulated in the heating process starting from the distribution of trapped H atoms, vacancies, and vacancy clusters shown in Figure 9. Comparing the simulated spectra to the experimental spectra in Figure 4, we found that the shape of Peaks 1 and 2 is mostly reproduced. However, we can see a spike at about 343 K, and the plateau part and the tail of the high-temperature side of Peak 2 are shorter than those of the experimental spectra. The change of trapped H atoms and of the concentration of vacancies and vacancy clusters in the heating process is shown in Figures 11 and 12. According to the figures, although the amount of trapped H atoms and the concentration of vacancies and vacancy clusters of the H + strain specimen are larger than those of the strain specimen before the heating process, their tendency to change during the heating process is similar. The trapped H atoms start to decrease before 273 K for Peak 1 and just after 273 K for vacancies and vacancy clusters, and by 398 K they become < 10 À8 ppm for all kinds of trap sites. Accompanying the decrease of H atoms trapped by V 1 ; V 1 starts to migrate and vacancy clusters are additionally generated. The generation of V 2 and V 3 is especially more remarkable than that of the others and the decrease of V 2 and V 3 is also so rapid. The rapid decrease causes the rapid release of H atoms trapped by V 2 and V 3 so that the rapid release forms the spike in the simulated spectra at about 343 K in Figure 10. According to Figures 11(a) and 12(a), we observed that the desorption spectra in the temperature range >343 K are formed by H atoms released from vacancy clusters>V 3 . Those vacancy clusters remain slight even beyond the temperature of 373 K and V 1 is also remaining. Figure 13 shows the thermal desorption spectra obtained by simulating the case that the degassing process of the H + strain specimen is replaced by the annealing process.
The figure corresponds to Figure 6. Comparing both figures, we revealed that the decrease of the desorption rate of Peak 2 against the aging temperature was almost reproduced. Thus, the model can properly represent the decrease of vacancies and vacancy clusters for temperature. However, the change of Peak 1 was not simulated in Figure 13. As for Peak 1, the peak temperature in Figure 4 is close to that of the corresponding peak in the spectra of the tempered martensitic steel [27] , and the H trapping energy evaluated in Section IV was close to that of screw dislocations. Thus, it can be considered that Peak 1 is formed by H atoms released mainly from dislocations as reported [27] . Figure 7 indicates that by the aging process at ‡ 323 K, the H trap sites decrease and the H trapping energy changes from the case of 303 K. If it is supposed that the defect for Peak 1 is dislocations, although its decrease by the aging process of high temperature is well known, the increase of its H trapping energy is not so obvious. It might relate to that dislocations absorb vacancies and vacancy clusters as sink sites. However, the mechanism of the change of  dislocations is not so clear that it can be incorporated into the model. In addition, since dislocation trap sites might not change at a temperature of about 281 K, which is the peak temperature of Peak 1, it is considered that their change does not almost influence the result of the TDS measurement. Therefore, we did not consider the change of defects for Peak 1 in the model, so that the change of Peak 1 was not simulated in Figure 13.

VI. DISCUSSION
The thermal desorption spectra simulated for the strain and H + strain specimens in Figure 10 have a spike at about 343 K, which is not seen in the experimental spectra in Figure 3. The spike resulted from the rapid desorption due to the rapid decrease of V 2 and V 3 in Figures 11 and 12. Since their rapid decrease is caused by the fast diffusion of V 2 and V 3 , we considered making the diffusion of both kinds of vacancy clusters slower by enlarging the diffusion barrier energy, Q V 2 and Q V 3 . Figure 14 shows the spectra obtained by the simulation using Q V 2 and Q V 3 , which were increased by 20%. In the figure, we observe that the spike is lower.
According to Figure 15, which shows the change of trapped H atoms and of concentration of vacancies and vacancy clusters, the rapid decrease of the concentration of V 2 and V 3 is mitigated so that the rapid release of the trapped H atoms is also mitigated. From the result, the diffusion of V 2 and V 3 in the specimens may be slower than that of the model.
A similar situation is also mentioned in the report [38] though it discusses the helium diffusion. Therefore, we may need to reconsider the diffusion coefficient that is estimated by the atomistic calculation. Furthermore, it may also be necessary to consider the possibility that solute atoms impede the migration of vacancy clusters as carbon atoms have been reported to make vacancy migration slower [39,40] . In addition, since the spike still remains slightly in Figure 14 even when the diffusion of V 2 and V 3 becomes slow, there is a possibility of other factors slowing the decrease of V 2 and V 3 . Thus, such other factors should be investigated in the future.
Next, we discuss the plateau part and the tail at the high-temperature side of Peak 2, which became shorter than that of the experimental spectra. Since it is   considered that the tail is influenced by vacancy clusters with comparatively large size, we simulated the thermal desorption spectra of the H + strain specimen by varying the maximum size of vacancy clusters, n, and examined the shape of the spectra. Figure 16 shows the magnified view of Peak 2 in the simulated spectra.
According to the figure, we observed that when vacancy clusters > V 4 are included, the tail at the high-temperature side becomes longer than when they are excluded. However, the influence of the cluster size is not so remarkable from V 5 up to V 9 because the concentration of clusters > V 5 becomes gradually smaller than that of V 4 and V 5 , as seen in Figure 9(b). Therefore, even if vacancy clusters > V 9 were included, their influence might not be so large under the current condition. As another possibility, if the ratio of vacancy clusters>V 4 is higher before the charging process, more H atoms can be trapped by such clusters in the charging process and might elongate the tail at the high-temperature side. To realize such a situation at the beginning of the degassing process, vacancy clusters need to already exist instead of just vacancies and grow larger until before the charging process. This situation may be consistent with the report [8] that vacancy clusters whose PAL is < 300 ps can be generated in the steel deformed with H atoms. Furthermore, H atoms trapped by many large clusters may form the plateau part in the spectra. However, we need further theoretical and experimental consideration for examining the amount and size of vacancy clusters before the degassing process.

VII. CONCLUSIONS
We obtained the H thermal desorption spectra of the small-size specimen with 0.3 mm thickness cut out from an iron specimen in which strain was applied and was charged with H atoms concurrently by the L-TDS measurement. Then, we simulated the experimental spectra using the model based on the McNabb-Foster theory, which considered the migration of vacancies and vacancy clusters as well as the growth and dissociation of vacancy clusters. The model considered up to V 9 and employed the parameters based on the atomistic calculations for the behavior of H atoms, vacancies, and vacancy clusters, including H trapping energy that we estimated using molecular static calculations. In addition, we applied the model to the simulation of the behavior of H atoms, vacancies, and vacancy clusters in the L-TDS pre-measurement processes, such as the degassing or annealing, polishing, and charging processes. As a result, we obtained the following.
(1) In the experimental spectra, the H desorption from defects generated by applying strain and charging with H atoms concurrently formed Peak 2 with the shoulder-like shape at temperature after that of Peak 1 corresponding mainly to dislocation. (2) From the experiment for the aging process of the H + strain specimen and the report [27] , it can be considered that the defects corresponding to Peak 2 are vacancies and vacancy clusters. (3) From the simulation for the behavior of H atoms, vacancies, and vacancy clusters in the pre-measurement processes, i.e., the degassing or  annealing, polishing, and charging with tracer H atoms, we obtained the size distribution of vacancies and vacancy clusters in the H + strain specimen before the heating process of TDS, and we observed that the concentration of V 2 and V 3 became smaller than that of V 1 and V 4 -V 6 because of the faster diffusion of V 2 and V 3 . (4) The simulation using the size distribution of vacancies and vacancy clusters as the initial state reproduced, on the whole, the shape of the experimental thermal desorption spectra. However, the spike that appeared at about 343 K, the length of the plateau part and the tail at the high-temperature side of Peak 2 were different from those of the experimental spectra. (5) We found that the spike at about 343 K is caused by the rapid H desorption as a result of the fast diffusion of V 2 and V 3 . In addition, from the dependence of the simulated spectra on the maximum size of vacancy clusters included in the model, it was indicated that if the ratio of vacancy clusters>V 4 were higher in the initial size distribution of the heating simulation, plateau part and tail at the high-temperature side might be improved.
From the last result, the diffusion of V 2 and V 3 should be reconsidered. Also, the amount and size of vacancy clusters that are generated in the specimens by applying strain and charging with H atoms concurrently should be estimated experimentally and theoretically. These results can be significant for the quantitative discussion of the influence of vacancies and vacancy clusters on HE.