Investigations into penetration depth profiles of hydrogenic species in beryllium plasma-facing components via molecular dynamics simulations

During the operation of nuclear fusion reactors, plasma-facing components lining the reactor vessel are continually bombarded by plasma species. The penetration and subsequent trapping of these bombarding plasma ions has implications for component damage as well as in-vessel inventory. Accurately predicting the expected ion penetration depth profiles at a range of plasma ion and surface temperatures typical of fusion reactor operating conditions will inform the scrape-off layer design to limit particle radiation damage and tritium trapping in order to prolong the lifetime of the plasma-facing components and satisfy the DT fuel cycle requirements. By defining a statistical distribution for ion penetration depth and describing the evolution of its parameters across the fusion parameter space of interest, the expected ion deposition depth profiles can be calculated for any subset of ion and surface temperature ranges as needed. Molecular dynamics simulations were used to study the bombardment of beryllium lattices with surface temperatures of up to 1100 K by 5 eV–150 eV deuterium and tritium ions, and the resulting ion penetration depths were investigated. The distributions of two penetration depth quantities, considered from the perspectives of lattice damage and hydrogen retention are defined and their distribution parameter dependence on surface and ion temperature is identified. The expected positive correlation between penetration depth and ion temperature is observed, where the non-linear relationship between these quantities indicates the expected form of the velocity dependence of nuclear stopping power at low bombardment energies. Isotope effects on the distributions are also investigated, with results suggesting that heavier ions have comparably lower mobility within the sample and will generally accumulate closer to the surface. A short study on ion deposition rates is also performed; a non-linear increase of deposition rate with increasing bombarding ion energy has been observed, and evidence of a weak positive surface temperature correlation has been noted.

During the operation of nuclear fusion reactors, plasma-facing components lining the reactor vessel are continually bombarded by plasma species.The penetration and subsequent trapping of these bombarding plasma ions has implications for component damage as well as in-vessel inventory.Accurately predicting the expected ion penetration depth profiles at a range of plasma ion and surface temperatures typical of fusion reactor operating conditions will inform the scrape-off layer design to limit particle radiation damage and tritium trapping in order to prolong the lifetime of the plasma-facing components and satisfy the DT fuel cycle requirements.By defining a statistical distribution for ion penetration depth and describing the evolution of its parameters across the fusion parameter space of interest, the expected ion deposition depth profiles can be calculated for any subset of ion and surface temperature ranges as needed.Molecular dynamics simulations were used to study the bombardment of beryllium lattices with surface temperatures of up to 1100 K by 5 eV-150 eV deuterium and tritium ions, and the resulting ion penetration depths were investigated.The distributions of two penetration depth quantities, considered from the perspectives of lattice damage and hydrogen retention are defined and their distribution parameter dependence on surface and ion temperature is identified.The expected positive correlation between penetration depth and ion temperature is observed, where the non-linear relationship between these quantities indicates the expected form of the velocity dependence of nuclear stopping power at low bombardment energies.Isotope effects on the distributions are also investigated, with results suggesting that heavier ions have comparably lower mobility within the sample and will generally accumulate closer to the surface.A short study on ion deposition rates is also performed; a non-linear increase of deposition rate with increasing bombarding ion energy has been observed, and evidence of a weak positive surface temperature correlation has been noted.

Introduction
Due to beryllium's low atomic number, good thermal properties and in-situ repairability [1,2], it has been identified as the material of choice for the plasma-facing components (PFCs) of many operational and planned future nuclear fusion reactors, whether as a bulk material or as a coating [3].During reactor operation, the PFCs are continuously bombarded by plasma species with energies of ∼100 eV [4], and up to ∼1 keV in some scenarios [5,6].This interaction between PFCs and the edge plasma results in the deposition and occasional implantation of plasma species onto the surface or into the bulk of the PFCs; accordingly, determining the deposition rates and penetration depth profiles at fusion-relevant plasma ion energies and PFC temperatures is vital for understanding contributions to PFC damage as well as the trapping and retention behavior of plasma species in the PFC [7,8].
The growing availability of high-performance computing resources has fueled the use of increasingly complex molecular dynamics (MD) models for studying plasmasurface interactions under fusion reactor-like conditions.This is evidenced by the recent developments of many new MD interatomic potentials purpose-built for fusion research [9][10][11][12].MD systems utilizing these novel potentials have been used to investigate the properties and behavior of tungstenbased alloys (containing tantalum [10], rhenium [11], or zirconium-carbide [12]) as possible candidate materials for the PFCs of future reactors, as well as to study the intermetallic beryllium-tungsten structures expected to form in ITER [9].
In this work, an experimentally validated MD model of a rough beryllium surface, the details of which are discussed in section 2, is used to estimate the penetration depth profiles of deuterium and tritium ions with energies of up to 150 eV impinging on beryllium PFCs with temperatures up to 1100 K.The distributions of two distinct penetration depth quantities are reconstructed: the maximal penetration depth (MPD), considered from the perspective of assessing lattice damage, and the final penetration depth (FPD) as considered from a hydrogen retention perspective.Further details regarding the definitions of these quantities and the statistical significance of the presented distributions are also included in section 2.
Despite the availability of existing methods for estimating the stopping power of materials and computing the resulting ion penetration ranges, namely the binary collision approximation (BCA) [13] method utilized by the popular collection of software packages SRIM [14], MD has been demonstrated to be significantly more suitable for studying sub-keV surface interactions such as those involved in this study [15][16][17].Additionally, the flexibility of MD allows for the modeling of more imperfect and subsequently more realistic systems, such as those where surface reconstruction has taken place, or where preparatory bombardments were used to artificially roughen the surface and induce damage to the bulk material.
In the following sections, the effects of plasma ion energy and beryllium surface temperature on individual penetration depth distribution parameters are quantified; the presented relationships, coupled with the tabulated data available in the appendices, will allow the reader to reconstruct custom penetration depth distributions based on their parameter ranges of interest.Section 3 additionally explores the isotope effects of bombarding beryllium by tritium, and the resulting changes to the penetration depth distributions that should be expected.The significance of the presented results, accompanied by deuterium ion deposition rates, are then discussed in section 4. Lastly, the observed relationships are compared against trends derived from data computed using the SRIM package [14], and conclusions are made regarding possible paths of research that can follow from this work in the future.

MD model
The MD model, implemented using the LAMMPS Molecular Dynamics MD Simulator [18], simulates 9786 beryllium atoms (standard atomic weight of 9.0122 g mol −1 [19]) organized in a hexagonal close-packed (HCP) lattice with the parameters a = b = 228.58pm and c = 358.43pm [20].The lattice, illustrated in figure 1, fills a 3-dimensional region spanning 40 Å by 40 Å in the periodic x and y dimensions and 50 Å in the non-periodic z dimension, thereby simulating the upper section of a monocrystalline beryllium PFC tile with an infinitely large surface and a thickness of 50 Å, oriented such that the ABAB planes are stacked along the z-axis.This system size was chosen to balance computational cost against finite-size effects arising from the use of periodic boundary conditions, with a focus towards a larger tile thickness that ultimately limits the highest lattice bombardment energy that can be investigated.Each atom in the lattice was assigned a random initial velocity sampled from a uniformly-distributed ensemble of velocities representing the lattice temperature T Be .Interactions between all atoms in the system were computed using a 3-body Tersoff potential [21], the parameters of which were adapted from an existing beryllium hydride analytical bond-order potential (ABOP) developed for studies of plasma-wall interactions in fusion reactors [22].This potential was chosen as it was, at the time of writing, the only well-established potential describing the interactions between beryllium and hydrogen that was intended for simulations of non-equilibrium processes and validated at a wide range of beryllium lattice temperatures.A flowchart illustrating the order of steps (grouped into 'stages') performed by the MD model (right), accompanied by a render of the MD system during the production stage (left).Rendered using OVITO [23], the image shows the beryllium lattice (red) bombarded by a deuterium ion (blue).The roughness of the surface, a deuterium atom absorbed into the bulk lattice, and structures resembling crystal grain boundaries on the exposed xz plane should be noted.In the flowchart, T H refers to the ion temperature of the hydrogenic species being investigated.

Equilibration stage
In order to induce surface reconstruction in the upper section of the bulk beryllium lattice, as well as to remove any artefacts arising from initial atom positions near the boundary conditions, the system was stepped forward in time and allowed to dissipate excess energy until an equilibrium state was reached.Time integration of the system was performed using the Verlet algorithm on Nosé-Hoover style equations of motion which sample from the canonical NVT ensemble.In order to reduce computation time while maintaining high accuracy in the initial equilibration stages, the timestep dt with which the integrator updates particle positions and velocities was set to be dynamic, with a lower limit of 1 attosecond and an average dt of 2-4 attoseconds depending on the lattice temperature once the system was near equilibrium.During this equilibration stage, the system progressed for a minimum of 200 picoseconds (although often requiring as long as 500 picoseconds), and was considered to be equilibrated once the average change in the total energy per beryllium atom was no more than 1 meV in the last 100 picoseconds of simulation time.

Surface preparation stage
During the surface preparation stage, as visualized in figure 1, the newly-formed 'perfect' beryllium surface is exposed to hydrogen plasma (in the form of repeated surface bombardment by deuterium or tritium ions), in an effort to induce damage to the surface and upper bulk sections of the beryllium tile.A deuterium or tritium ion (isotope masses of 2.0141 g mol −1 and 3.0160 g mol −1 [24], respectively) was created at a random position 25 Å above the surface sufficiently away from its influence.At the modeled energies and proximity to the surface, the 'ion' can be assumed to be uncharged as a result of the The energy E per Be atom in 8 equilibrated trial lattices of periodic lengths L and a fixed non-periodic length of 50 Å were computed (black) and fitted to an exponential decay curve (green solid).The asymptote (green dashed) shows the ideal periodic length to be around L = 44.70Å.
Auger neutralization process [25].The ion was then assigned a velocity with a magnitude equivalent to a plasma ion temperature T H of 50 eV directed towards the surface with a maximum approach angle of 27 • to the normal of the surface; as ions get accelerated towards the PFC while in the plasma sheath region, the maximum surface approach angle was chosen by limiting the individual velocity angular components in x and y to 20 • [26] (bombardment events with shallower surface approach angles are increasingly infrequent and more likely to result in surface reflection regardless).To simulate the energy transfer between the bombarding ion and the beryllium lattice, all atoms in the system during this stage operate under a constantenergy NVE regime, with only the bottom 5 Å of the lattice still sampling from a constant-temperature NVT ensemble and thereby acting as a heat sink.Upon impact with the lattice, the ion either reflects off the surface, adsorbs onto the surface, or is absorbed into the bulk of the beryllium tile (the latter case can be seen in figure 1).This excess energy propagates through the lattice into the lower region where it is gradually dissipated, eventually returning the system to a new equilibrium state.Repeated surface bombardments induce damage to the surface via the sputtering of surface species or adsorption of plasma species to the surface, as well as damage to the top layers of the bulk beryllium lattice in the form of trapped hydrogen atoms, lattice vacancies or dislocations.
Under fusion-like conditions, bombardments of the 16 nm 2 of surface area exposed by this model occur significantly less frequently than can be feasibly simulated within the timescales of this model.However, it can be assumed that the system spends most of its time between successive bombardments in an equilibrium state, which does not need to be simulated.As a result, it is satisfactory to start the next surface bombardment immediately after an equilibrium state from the previous bombardment is reached.For the system presented here, the beryllium slab has been bombarded 20 times with a reequilibration time of 25 picoseconds between each bombardment.The state of the 'prepared' system (including all trapped hydrogen atoms) was then saved, and all successive bombardments during the production phase started from this state and reset back to this state once adequate data from each bombardment were collected.A total of five individual systemsone for each investigated temperature-were prepared using this method.

Validation stage
In order to ensure that the system is not sensitive to the choice of internal model parameters or the conditions that were used to prepare it, a number of additional 'validation' systems were tested.To determine the optimal lattice size while maximizing the length of the non-periodic z-dimension and minimizing the total number of atoms in the system (a number which contributes significantly to the system's overall computational cost), 8 systems with a fixed non-periodic z length of 50 Å and varying periodic x and y lengths between 5 Å and 40 Å were equilibrated, and the final energy per beryllium atom in each lattice was measured.As shown in figure 2, an exponential decay curve of the form E = a + b exp (−L/c) was fitted to the computed average energy E per Be atom in each trial system with a periodic length L (where a, b, and c are the fitted parameters), and the ideal periodic size of a system with a fixed non-periodic length of 50 Å was computed to be L = 44.70Å approaching the curve's asymptote to 2 decimal places.As such, a periodic length L = 40 Å was selected as being sufficiently close to the thermodynamic limit.
Additional investigations were performed to evaluate the effects of the thermal damping time t damp used as a part of the Nosé-Hoover NVT thermostat, as well as the bombardment energy at which the surface of each system was prepared.Two test systems with 0.1× and 10× the nominal t damp value of 100dt (approximately 403 attoseconds considering the average dynamic timestep dt = 4.03 attoseconds at equilibrium), as well as one system with a surface preparation bombardment energy of 200 eV (4× the nominal value of 50 eV) were equilibrated and their surfaces prepared; each system was then bombarded by a number of 5-50 eV deuterium ions as described in the surface preparation subsection, and the sputtering behavior and penetration depth profiles were compared.To compare the sputtering behavior, the sputtering yield at each bombardment energy was fitted to a universal sputtering relation [27], and the fitting parameters Q and threshold energy E th were found to be comparable to within one standard deviation of each other for all test systems.To compare the penetration depth profiles, the technique described in the following production and data analysis subsection was used, and no significant change in behavior between the test systems was observed.
As systems with different lattice temperatures were used throughout this study, the effect of lattice temperature on the density of the beryllium tile was also considered to play a role in the bombarding ion penetration behavior.In order to investigate the thermal expansion of the lattice, a radial distribution function (RDF) was computed for 5 equilibrated systems with lattice temperatures between 300 K and 1100 K, and the average interatomic nearest-neighbor (NN) separation corresponding to the first RDF peak position was estimated for each system.An increase in NN separation with temperature has been observed, from the base value of (2.2566 ± 0.0003) Å at T Be = 300 K to (2.273 ± 0.001) Å at T Be = 500 K.However, as shown in figure 3, there is no evidence of further lattice expansion at higher temperatures, presumably as a result of the use of periodic boundary conditions and a fixed system size restricting this behavior.A possible workaround for this is discussed in the conclusion.

Production and data analysis stages
An incident ion impacting a beryllium surface will either reflect off the surface, adsorb onto the surface, or penetrate the surface and enter the bulk of the lattice.In the case of reflection, the ion is not in contact with the lattice and thus penetration depth is undefined; however, the rate of reflected ions can be used to compute the deposition rate, a quantity that is useful when paired with penetration depth distribution profiles for estimating the absolute values of trapped ion inventory or bombardment-induced damage.In the cases of adsorption or absorption, the bombarding ion penetration depth can be calculated with respect to the exact z-position of the lattice top surface.
To compute the surface z-position of a perfect HCP lattice with its ABAB planes stacked along the z-axis, a histogram of atom z-positions near the surface can be compiled that will appear as a set of normal distribution probability density functions, the peak of each corresponding to the z-position of an ABAB plane; the maximum of the peak with the highest zvalue will represent the uppermost surface plane.In imperfect scenarios where surface reconstruction has occurred, atoms have been deposited on the surface, or the lattice planes are at a slight angle with respect to the z-axis, a generalized normal distribution [28] with a s > 2 'flat top' will be a better fit to the histogram peaks.Additionally, should the layers be skewed sufficiently and the individual histogram peaks overlap significantly, a convolution of two (or more) such distributions will aid the fitting process; this method of surface identification is visualized in figure 4.
With the knowledge of the lattice surface position, two definitions of penetration depth (also illustrated in figure 5) can be established: the MPD, defined as the maximum depth that a bombarding ion reaches below the lattice surface, which considers penetration depth from the perspective of lattice damage; and FPD, defined as the bombarding ion's final depth below the surface after it has transferred all of its excess energy to the lattice (sufficient 'settle time' has elapsed), and which instead considers penetration depth from the perspective of hydrogen retention.
Outgassing effects present themselves as the gradual ascent of the trapped ion towards the lattice surface.The rate of the ion's ascent, calculated as the gradient of the ion's zposition after settle time, can be used as an indication of its expected residence time before it reaches the surface and  Graphical representations of the two defined penetration depth quantities, calculated using a sample trajectory of a 150 eV deuterium ion bombarding a 300 K beryllium surface.The magenta dashed line represents the maximum penetration depth (MPD), while the green line represents the final penetration depth (FPD).An outgassing rate is computed by calculating the ion's rate of ascent (if any) from the settle time at which it is expected to have dissipated all of its excess energy (dashed red line); this rate is used to estimate the bombarding ion's corrected FPD at the time of impact (dashed blue line).subsequently desorbs.The FPD is also corrected for these outgassing effects by using the gradient to compute the ion's expected penetration depth at the time of impact.Therefore, unlike the MPD which describes a maximum value irrespective of simulation time, the FPD is always measured at the time of impact.This correction is graphically represented in figure 5.It should be noted, however, that such outgassing effects are expected to be negligible at the simulated timescales.
Throughout this study, a statistically significant number of lattice bombardments was carried out for each investigated lattice temperature and plasma ion temperature (100 bombardments for the majority of points in the explored parameter space, with up to 300 bombardments for T Be = 300 K), and the bombarding ion's MPD and FPD was recorded for each bombardment.As outlined in the previous subsection, each bombardment involves creating a hydrogen ion randomly above the lattice surface, assigning it a velocity with a known magnitude and a random but constrained surface approach angle, and allowing the ion to impact the beryllium surface.Each bombardment, however, starts from the same state, and only the bombarding ion's initial position and surface approach angle is randomly varied.As a result, for each point in the parameter space, a distribution of MPD and FPD values is collected.The Moyal distribution [29], initially proposed for modeling the energy loss of fast charged particles by ionization of the target medium and used throughout high-energy physics [30], has been found to fit well to both penetration depth quantities.
In order to study the effects of surface temperature and bombarding ion energy on the two penetration depth distributions (as well as to compare between them), Moyal probability density functions (PDFs) were reconstructed from the collected penetration depth data for both the MPD and FPD by minimization of the negative log-likelihood function.Being a two-parameter distribution, as described by equation (1), the shape of any Moyal PDF can be described by its location parameter µ and scale parameter σ; these parameters are used throughout the following sections to more concisely describe the change in the shape of the distribution throughout the investigated fusion parameter space.In terms of the FPD, as an example, an increase in the location parameter µ indicates that bombarding ions tend to settle deeper below the surface, while an increase in the scale parameter σ indicates that the ions' final depth tends to be more varied; similar logic can also be applied to the MPD.The tabulated parameter values can also be useful to anyone who wishes to reconstruct their own expected penetration depth distribution profile based on their parameter ranges of interest.

Results
Five deuterium lattices with surface temperatures T Be between 300 K and 1100 K were equilibrated and prepared as outlined in the previous section.Each lattice was then repeatedly bombarded by deuterium ions with energies representing plasma ion temperatures T D between 5 eV and 150 eV.Between 80 and 300 unique bombardments were carried out for each of the 5 surface temperatures and 24 bombarding ion energies, and the resulting MPD and FPD data were fitted to Moyal PDFs for comparison and further analysis.A summary of these results can be seen in figure 6; in subfigures (a) and (b), a graphical representation of the evolution of the shape and peak position of the MPD and FPD distributions with increasing T D can be seen, while subfigures (c) and (d) describe the behavior of this evolution in terms of the distribution parameters.By keeping the surface temperature T Be constant and increasing the deuterium bombarding energy T D , a shift in the peaks of both the MPD and FPD distributions towards greater depths below the surface can be observed, as illustrated in figures 6(a) and (b).This result, as expected, is due to the fact that more energetic particles are able to travel deeper below the surface into the bulk of the beryllium lattice before dissipating their energy.By graphing the location parameter µ of both distributions against T D as shown in figure 6(c), it can be noted that this tendency towards greater depth at higher energies appears to follow a relationship described by equation ( 2), where a, b and c are fitting parameters.During the fitting procedure, a single value of the parameter c was shared between the MPD and FPD curves for each of the distribution parameters; see table 1 for details of the fitted parameter values.
Comparing µ between MPD and FPD in figure 6(c), it can be noted that MPD locations have a consistently higher value at the same T D ; indeed, this behavior is also visible in figures 6(a) and (b), where all MPD PDFs are consistently shifted to the right when compared to FPD PDFs at the same T D (as visualized by the dashed lines).This is due to the fact that, as MPD is a measure of maximal penetration depth, no FPD value can by definition be greater than the corresponding MPD value for a single bombardment.However, this behavior also suggests that bombarding ions consistently backscatter towards the surface once they have sufficiently low kinetic energy.Secondly, it can be noted at low T D , the FPD location is negative; this agrees with the observed behavior that at low bombardment energies, the bombarding ions will adsorb and rest on the surface rather than be absorbed into the bulk of the lattice.However, MPD µ always remains positive, suggesting that even at T D conditions as low as 5 eV, the bombarding ion will perturb the surface and reflect off lower lattice layers before returning to rest on surface.A final observation that can be made relates to the nonlinearity of the µ(T D ) relationship, which is a result of the fact that the stopping power experienced by the bombarding ion once inside the bulk of the lattice increases with the ion's velocity in the investigated range of T D energies [31]; this is discussed further in the following section.
Comparing instead the scale parameter σ between MPD and FPD as shown in figure 6(d), the parameter is observed to increase nonlinearly with increasing T D for both penetration depth metrics; this is apparent as a widening of the distribution PDF, as visualized in figures 6(a) and (b).This widening is a result of the fact that higher bombardment energies allow the bombarding ions to initially penetrate deeper into the lattice, but also increase the total number of collisions necessary for the ion to dissipate its energy, thereby increasing the likelihood that the ion will also be backscattered and redirected towards the surface.An increase in T D should therefore be thought of as allowing access to a larger range of penetration depths, rather than simply resulting in increased penetration.It can also be noted that the distribution parameter σ is generally higher for FPDs across all investigated T D , or rather, FPD PDFs are generally wider than MPD PDFs at the same T D ; this behavior indicates that, between individual bombardments at the same bombarding ion energy, the range of maximum depths reached by the ions is more consistent than their range of FPDs.
The effect of beryllium surface temperature on all penetration depth measurements appears to be negligible.As a result, all results quoted in this section are computed using aggregated data from all simulated surface temperatures.Possible causes of this behavior are discussed in section 4.

Deposition rates
As discussed previously in section 2, only a fraction of ions bombarding the beryllium lattice will remain absorbed in the bulk lattice or adsorbed on the surface; some ions will reflect off the surface immediately, and some may backscatter due to collisions with atoms in the bulk lattice while remaining energetic enough to reach the surface and desorb.Although not the primary target of investigation during this study, the total number of ions that remain deposited in or on the sample after dissipating all of their excess energy has been measured for 5 lattice temperatures T Be .These deposition rates, as shown in figure 7, appear to increase linearly with T D in the investigated Figure 6.Summary of the compiled MPD and FPD distribution profiles.Subfigures (a) and (b), using data from the T Be = 300 K lattice, illustrate the widening and shifting of the distributions towards the right with increasing T D , while subfigures (c) and (d) describe the evolution of the distribution shapes with T D using the distribution parameters µ and σ.Vertical dashed lines are included to aid in the comparison between MPD and FPD peak positions in (a) and (b), respectively; this behavior, where MPD peaks (being a measure of a maximum) are consistently found at greater depths below the surface, can also be seen in subfigure (c).
Table 1.Fitted parameters a, b and c, as well as their corresponding one standard deviation errors of fit, describing the relationship given by equation (2) between the MPD and FPD distribution parameters and the bombarding deuterium energy T D .Only a single parameter c was fitted for each of the distribution parameters µ and σ.Due to negligible surface temperature effects as well as to improve the statistical significance of the discussed relationship, penetration depth data from all investigated lattice temperatures were aggregated for this fit.range of deuterium plasma ion temperatures; with increasing T D , bombarding ions are more likely to penetrate the surface and subsequently dissipate their energy while surrounded by atoms in the bulk lattice, giving rise to this positive trend.However, this behavior also suggests the necessary existence of a plateau in deposition rates at significantly higher values of T D , as discussed further in section 4.
A weak temperature dependence can also be noted in figure 7, highlighted by the fitting of a straight line with a fixed gradient to the measured deposition rates of each investigated lattice temperature.A small increase in deposition rate with increased lattice temperature is visible; this behavior, if confirmed, is probably caused by the decreased lattice surface density at higher lattice temperatures, which results in an increased probability of surface penetration.However, due to the noisy nature of the collected data in this energy range, this trend should only be considered indicative at this stage.

Isotope effects
As early future fusion plasmas are expected to consist of a mixture of both deuterium and tritium, this penetration depth study has been extended to also involve a single 300 K lattice bombarded by tritium ions with energies ranging between 5 eV and 150 eV.The results are plotted and compared against deuterium data in figure 8, and the corresponding fitted parameter values are available in table 2. Primarily, it can be noted that the tritium MPD and FPD distribution location parameter µ grows more slowly with increasing bombarding ion energy as compared to deuterium; tritium ions can therefore be expected to accumulate closer to the beryllium's surface than deuterium at comparable plasma ion temperatures.This effect is theorized to be the result of tritium experiencing stronger stopping power due to its mass as it participates in collisions with other lattice ions in the bulk.Similar behavior is visible for the distribution scale parameter σ, which increases at a slower Deposition rates obtained by analyzing all bombardment trajectories from previous penetration depth studies.The rates are fitted to a straight line with a common gradient and a varying intercept.In this energy region, deposition rates appear to increase by 8 ± 1% with every 100 eV of T D , and have a baseline value of 59 ± 1% for the 300 K lattice and 61 ± 1% for the 1100 K lattice.on MPD and FPD distributions, respectively.Distributions for both deuterium (solid lines) and tritium (dashed lines) are included for comparison, with dashed vertical lines marking tritium MPD and FPD peaks for ease of comparison between the two subfigures.Subfigures (c) and (d) show the effects of bombarding ion temperatures on the individual penetration depth distribution parameters.Markers for parameter values at individual bombarding ion energies are included for tritium only.The 10-25 eV region of high noise apparent in subfigure (d) is due to the large density of data points collected in that energy range (with the aim of estimating the sputtering threshold energy [27], see section 2.3 for more information).rate compared to deuterium, resulting in taller tritium MPD and FPD PDFs compared to deuterium at similar bombarding ion energies.Due to the decreased mass ratio between the bombarding tritium ion and the beryllium lattice atom, more efficient momentum exchange results in the tritium ion experiencing fewer collisions than deuterium before dissipating its excess energy, decreasing its overall range of mobility in the sample.
Table 2. Fitted parameters a, b and c with their corresponding one standard deviation errors of fit, describing the relationship given by equation ( 2) (where T D has been replaced by T T ) between the penetration depth distribution parameters µ and σ, and the tritium ion energy T T .As in table 1, only a single parameter c was fitted for each of the distribution parameters.However, unlike the aggregated deuterium dataset, the data points plotted and fitted here represent only a single 300 K beryllium surface bombarded by tritium.

Discussion
As discussed in section 2, the non-periodic height or z-size of the beryllium lattice has been set to 50 Å in order to make it possible to perform a sufficiently large number of bombardments to derive statistically significant conclusions.This size, however, imposes a maximum limit on the bombarding ion energy at which penetration depth can be studied; ions with a high enough energy will simply pass through the bulk of the lattice and exit via the lower 'surface' -a state for which penetration depth is undefined despite the fact that neither reflection nor backscattering has occurred; such ions are still considered to be deposited, as the modeled z-size represents only the uppermost section of a much thicker beryllium tile that exists below the virtual 'lower surface'.Despite this, as deposition rates are primarily influenced by interactions near the surface, they can be investigated for an extended range of bombarding ion energies.It should be noted, however, that deposition rates estimated in this way will be overestimated as high-energy bombarding ions will often pass through the surface without having a chance to backscatter.With this in mind, additional deuterium bombardments with energies up to 1 keV were performed and the deposition rates were computed, as shown in figure 9. Due to the aforementioned weak beryllium temperature dependence (see figure 7), deposition rate data were aggregated over all surface temperatures and fitted to a curve of the form y = a − b exp(cx), where a, b and c are fitting parameters.It can be noted that at this extended bombarding energy range, the relationship no longer appears linear, and instead plateaus at 83 ± 1%.To verify this behavior, a reference deposition rate dataset was also computed using TRIM [14] (stopping power version SRIM-2008, using the 'surface sputtering' damage model on a 1 µm thick Be layer with default density and displacement energy values); 99 999 ion bombardments were performed at 25 deuterium ion energies up to 1 keV and at 6 surface approach angles between normal surface incidence and 25 • .This reference dataset, although fitting well to a curve of the same form, predicts significantly higher deposition rates of up to 98 ± 3% and is a result of TRIM operating outside of the energy range at which the BCA algorithm is effective [13].
The stopping power experienced by a bombarding ion traveling through a medium refers to the average energy it dissipates per unit path length traversed in that medium; the energy can be dissipated via inelastic collisions between bound electrons (electronic stopping power) or via elastic collisions between the ion and nearby lattice nuclei (nuclear stopping Comparison between the TRIM ion range parameter (computed from 99 999 bombardments by deuterium ions up to 150 eV at a range of surface approach angles up to 25 • normal to the surface) against the maximum likelihood penetration depths for the MPD and FPD distributions.The proposed penetration depth dependence on bombarding ion energy appears to fit well to both the model and reference datasets.power).For bombarding ions with a mass comparable to target lattice atoms, and especially at sub-keV bombarding energies, nuclear stopping power is the primary pathway for energy dissipation in the lattice; indeed, for deuterium ions bombarding beryllium with energies up to eV, electronic stopping power can be considered to be negligible [32].All MD models emulate nuclear stopping power via the interatomic potential, and 'frictional' forces may be added to simulate electronic stopping power.The MD model in this study, utilizing solely a Tersoff potential, implements no electronic stopping power effects.Any observed isotope effects on stopping power, such as those seen in figure 8 where tritium appears to experience greater stopping power than deuterium at comparable ion energies, are therefore solely a result of the decreased mass ratio between the bombarding ion and target lattice atom.Despite this, it should be noted that isotope effects on electronic stopping exist and are applicable at sufficiently high energies and mass ratios [33].
As a final consideration, in order to verify the form of the relationship between the bombarding ion energy and the MPD/FPD distribution parameters calculated in the previous section, a comparison can again be made to reference data generated using TRIM.Reference ion range and straggle data were extracted from a subset of the TRIM dataset described at the beginning of this section (bombarding ion energies up to 150 eV).For easier comparison between the distributions outputted by TRIM and the Moyal distributions used in this study, the most likely penetration depths (corresponding to the PDF peaks) for both the MPD and FPD distributions were compared against the TRIM 'ion range' parameter, as shown in figure 10.The reference TRIM datasets indeed fits well to the proposed relationship; the magnitude of the penetration depth overestimation by TRIM is inconsequential as it is known to be erroneous-TRIM predicts an unfeasible average penetration depth of 4 Å at a bombarding deuterium energy of only 5 eV.

Conclusion
In this work, a model beryllium lattice was bombarded by deuterium and tritium ions, and the effects of bombarding ion energy and beryllium surface temperature on bombarding ion penetration depth profiles and ion deposition rates were investigated.Two measures of penetration depth have been defined: the MPD that considers ion bombardment from the perspective of PFC damage by particle radiation, and the FPD that, paired with the knowledge of deposition and outgassing rates, can be used for estimating the trapping and retention of hydrogen species in beryllium PFCs.The Moyal distribution, described by its location parameter µ and scale parameter σ, has been found to fit well to both the MPD and FPD distributions.The expected increase of both penetration depth measures with bombarding ion energy has been observed, while neither the MPD nor the FPD appears to be affected by beryllium surface temperature.The relationship between the penetration depth distribution parameters and the bombarding ion energy has been shown to follow the form y = a − b exp (cx), where y is the MPD or FPD µ or σ distribution parameter, x represents the deuterium or tritium ion temperature, and a, b and c are fitted parameters whose values are tabulated in tables 1 and 2 up to deuterium and tritium energies of 150 eV, respectively.The presented data allows for the reconstruction of penetration depth profiles for any combination or range of D/T ion energies up to 150 eV without the need for additional MD simulations.
The increase in ion deposition rate with bombarding ion energy also appears to follow the aforementioned form y = a − b exp (cx) up to deuterium energies of 1 keV.The form of this relationship is in agreement with the reference dataset computing using TRIM.A small but consistent increase in deposition rate has also been noted with increasing beryllium surface temperature, although modifications to the model will be needed to verify this behavior; the effects of beryllium temperature are mainly expected to manifest via thermal lattice expansion resulting in decreased lattice density, but this behavior was limited to the surface as a result of the periodic boundary conditions.Future models may be able to overcome this limitation by using a lattice created with modified lattice parameters that already account for the expected thermal lattice expansion; alternatively, using an NPT ensemble (or a similar regime that modifies the system volume as needed) during system equilibration should be considered, if feasible.
Several scenarios have also been identified that could potentially modify the deuterium and tritium penetration depth distributions in ways that the current model does not yet account for; as an example, thin layers of beryllium oxide were observed to form readily over beryllium PFCs even at low oxygen pressures [34].Aside from modifying the surface composition, these oxide layers were also noted to promote deuterium trapping [35].Neutron-induced damage to beryllium PFCs is also suspected to increase tritium retention [36].
Additionally, the results presented in this study open pathways for further research; the orientation of the beryllium lattice can, for example, have an effect on the observed penetration depth profiles.Exposing particular crystallographic orientations to the surface may allow bombarding ions to channel between the crystal planes and experience significantly reduced stopping power.Inversely, alternative lattice orientations could prove beneficial by preferentially backscattering bombarding ions.In any case, analysis of existing bombarding trajectory data with a focus on ion-induced collision cascades could provide additional information regarding the effects of ion energy and beryllium temperature on damage propagation.Lastly, investigations into whether the Moyal distribution and its discussed parameter evolution can be extended to penetration depth distributions of alternative ion-target pairs with a similar mass ratio may highlight a new method for rapidly estimating penetration depth profiles of low-energy ion bombardments.

Appendix A. Deuterium MPD, FPD and deposition rate data
The tables below summarize the data obtained from the analysis of all deuterium bombardments, sorted by the deuterium plasma temperature T D and grouped by the beryllium surface temperature T Be .The total number of bombardments performed by T D ions on a T Be lattice is given by N total .For computation of the deposition rate, three values are provided: N sorb , N refl and N trans , representing the total number of adsorbed and absorbed ions (ion remains in lattice), reflected ions (ion leaves via top surface), or transmitted ions (ion leaves via bottom surface), respectively, such that N total = N sorb + N refl + N trans .For the reconstruction of penetration depth distributions, the Moyal distribution location µ and scale σ parameters are provided for both the maximal penetration depth (MPD) and final penetration depth (FPD).

Appendix B. Tritium MPD, FPD and deposition rate data
The table below summarizes the data obtained from the analysis of all tritium bombardments, sorted by the tritium plasma temperature T T .See appendix A for a description of the tabulated values.

Figure 1 .
Figure 1.A flowchart illustrating the order of steps (grouped into 'stages') performed by the MD model (right), accompanied by a render of the MD system during the production stage (left).Rendered using OVITO[23], the image shows the beryllium lattice (red) bombarded by a deuterium ion (blue).The roughness of the surface, a deuterium atom absorbed into the bulk lattice, and structures resembling crystal grain boundaries on the exposed xz plane should be noted.In the flowchart, T H refers to the ion temperature of the hydrogenic species being investigated.

Figure 2 .
Figure 2.The energy E per Be atom in 8 equilibrated trial lattices of periodic lengths L and a fixed non-periodic length of 50 Å were computed (black) and fitted to an exponential decay curve (green solid).The asymptote (green dashed) shows the ideal periodic length to be around L = 44.70Å.

Figure 3 .
Figure 3. Radial distribution functions computed for the coldest (T Be = 300 K) and hottest (T Be = 1100 K) equilibrated beryllium lattice.Each peak corresponds to the increased probability of finding an atom situated a distance r away from any other atom in the lattice, with the first peak corresponding to the nearest-neighbor (NN) distance.A normally distributed probability density function has been fitted to each NN peak, estimating the NN separation to be (2.2566 ± 0.0003) Å and (2.264 ± 0.003) Å for the 300 K and 1100 K lattice, respectively.

Figure 4 .
Figure 4.A histogram of atom z-positions within 2 Å of the predicted surface position at z = 100 Å, compiled from 250 successive lattice snapshots of an equilibrated system with a lattice temperature of 300 K.A convolution of two generalized normal probability density functions was fitted to the histogram (black solid line), with the rightmost peak corresponding to the surface position (green dashed line).

Figure 5 .
Figure 5. Graphical representations of the two defined penetration depth quantities, calculated using a sample trajectory of a 150 eV deuterium ion bombarding a 300 K beryllium surface.The magenta dashed line represents the maximum penetration depth (MPD), while the green line represents the final penetration depth (FPD).An outgassing rate is computed by calculating the ion's rate of ascent (if any) from the settle time at which it is expected to have dissipated all of its excess energy (dashed red line); this rate is used to estimate the bombarding ion's corrected FPD at the time of impact (dashed blue line).

Figure 7 .
Figure7.Deposition rates obtained by analyzing all bombardment trajectories from previous penetration depth studies.The rates are fitted to a straight line with a common gradient and a varying intercept.In this energy region, deposition rates appear to increase by 8 ± 1% with every 100 eV of T D , and have a baseline value of 59 ± 1% for the 300 K lattice and 61 ± 1% for the 1100 K lattice.

Figure 8 .
Figure 8.Comparison between deuterium and tritium penetration depth data.Subfigures (a) and (b) show the effects of bombarding ionon MPD and FPD distributions, respectively.Distributions for both deuterium (solid lines) and tritium (dashed lines) are included for comparison, with dashed vertical lines marking tritium MPD and FPD peaks for ease of comparison between the two subfigures.Subfigures (c) and (d) show the effects of bombarding ion temperatures on the individual penetration depth distribution parameters.Markers for parameter values at individual bombarding ion energies are included for tritium only.The 10-25 eV region of high noise apparent in subfigure (d) is due to the large density of data points collected in that energy range (with the aim of estimating the sputtering threshold energy[27], see section 2.3 for more information).

Figure 9 .
Figure 9.Extended deposition rate data obtained by performing additional MD bombardments (black) and aggregating over all beryllium surface temperatures.Reference data (cyan) were generated using the TRIM software package and the SRIM-2008 stopping power database.Data from both sources were fitted to curves of the form y = a − b exp(cx), where a, b and c are fitting parameters.As BCA results should only be treated qualitatively at low energies[13], the overestimated ion deposition rates predicted by TRIM still agree with MD model results on the form of the relationship.

Figure 10 .
Figure10.Comparison between the TRIM ion range parameter (computed from 99 999 bombardments by deuterium ions up to 150 eV at a range of surface approach angles up to 25 • normal to the surface) against the maximum likelihood penetration depths for the MPD and FPD distributions.The proposed penetration depth dependence on bombarding ion energy appears to fit well to both the model and reference datasets.

Table 3 .
T Be = 300 K bombarded by deuterium ions.

Table 4 .
T Be = 500 K bombarded by deuterium ions.

Table 5 .
T Be = 700 K bombarded by deuterium ions.

Table 6 .
T Be = 900 K bombarded by deuterium ions.

Table 7 .
T Be = 1100 K bombarded by deuterium ions.

Table 8 .
T Be = 300 K bombarded by tritium ions.