On the use of recombination rate coefficients in hydrogen transport calculations

The commonly accepted picture for the uptake of hydrogen isotopes (HIs) from the gas phase across the surface into a metal with an endothermic heat of solution for HIs is that of dissociation followed by thermalisation in a chemisorbed surface state and finally overcoming a surface barrier to enter the metal bulk where the HIs occupy interstitial solute sites. To leave the metal bulk the HIs first transition to the chemisorbed surface state from which they then enter gas phase by recombining into a diatomic molecule. This model is generally attributed to the work of Pick and Sonnenberg from 1985. They clearly distinguish surface states and subsurface solute sites where the recombination flux is proportional to the square of the concentration of chemisorbed atoms due the diatomic nature of this Langmuir–Hinshelwood process. In an effort to compare their extended model with an earlier surface model by Waelbroeck, which uses an expression for the recombination flux proportional to the square of the sub-surface interstitial solute concentration, they derive an effective recombination coefficient. However, also with the so-derived Pick and Sonnenberg recombination coefficient, the Waelbroeck model is only applicable under certain conditions. But, due to its simplicity, it is often used in boundary conditions of diffusion trapping type calculations, generally ignoring whether or not these conditions are met. This paper will use the full Pick and Sonnenberg model implemented in the TESSIM-X code and in simplified algebraic approximations, to show the limits of applicability of the Waelbroeck–Ansatz in modelling hydrogen transport in metals foreseen for the first wall of magnetic confinement fusion devices.


Introduction
The uptake of hydrogen isotopes (HIs) into metals has been studied extensively in literature and different models have been proposed and compared [1][2][3][4][5][6][7]. In [6] different models for the uptake have been compared and among those models, the one * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
proposed by Pick and Sonnenberg in 1985 [7] is currently the accepted picture for uptake from the gas phase and release of retained HIs across the metal surface. They describe the uptake from the gas phase into a metal with an endothermic heat of solution for HIs as a two step process: first the HI molecules dissociate and chemisorb in sites at the surface where they thermalise with the surface. Then, the chemisorbed HI-atoms overcome a barrier (comparable to the heat of solution ΔE Sol ) in a thermally activated Arrhenius process and enter the bulk where they occupy solute sites through which the surface couples to the continuum bulk diffusion. The reverse process of leaving the metal bulk across the surface is, therefore, also a two step process: first the solute HI-atoms, in an Arrhenius process, overcome a barrier (comparable to the activation energy of diffusion ΔE Diff ) to enter the chemisorbed surface state. Then, they recombine in a Langmuir-Hinshelwood process to an HImolecule and desorb into the gas phase. This model involves two HI concentrations: the concentration CA of chemisorbed atoms at the surface and the concentration CS0 of HI in subsurface solute sites. Since usually ΔE Sol > ΔE Diff this process is asymmetric: the transition rate constant from surface to bulk is much smaller than the transition rate constant from bulk to surface. The Pick-Sonnenberg model has been investigated in literature (e.g. [8,9]) and expressions for surface limits and steady-state solutions were derived. However, due to its complexity it has not been used extensively in diffusion trapping modelling for the interpretation of laboratory experiments on retention and permeation or to make predictions for future fusion devices. But newly developed codes such as MHIMS [10][11][12] or HIIPC [13] and, as will be described in this publication, TESSIM-X [14,15] feature implementations of the full Pick-Sonnenberg model. In an earlier model by Waelbroeck et al [5] this asymmetry in the transition rates is ignored and the release to the gas phase is assumed to be proportional to CS0 2 based on the hypothesis that CA ∝ CS0. This hypothesis alleviates the need to treat the chemisorbed surface states, thus greatly simplifying the model description. The Waelbroeck model has been used extensively in literature to interpret gas-driven permeation experiments, e.g. [16,17], ion-driven retention measurements e.g. [18] and to derive surface limits, e.g. [19]. In the model by Waelbroeck the release from the surface is described by Γ WB Rec ∝ K r CS0 2 and the governing parameter K r is called 'recombination coefficient' despite the fact that it does not really describe recombination in the sense of the actual Langmuir-Hinshelwood process occurring from the chemisorbed surface states. Pick and Sonnenberg investigated the meaning of K r in view of their more-extended surface model and derived an expression for K r that is commonly referred to as the Pick-Sonnenberg recombination coefficient K PS r . Due to its simplicity the combination of Waelbroeck's Ansatz for the 'recombination' flux Γ WB Rec from the surface to the gas phase together with Pick-Sonnenberg's expression for the 'recombination rate' coefficient K PS r has been used extensively in literature to derive boundary conditions for diffusion trapping type calculations, e.g. [18,20]. However, the limits of this simplified description of the underlying surface processes are hardly ever touched upon. This paper will compare algebraic steady-state solutions of the Pick and Sonnenberg surface model with results obtained from the simplified Waelbroeck-Ansatz. Both of these steadystate solutions ignore surface saturation effects to yield simple algebraic expressions. To investigate the influence of surface saturation the algebraic results will then be compared to the final state of fully time-dependent simulations of the Pick and Sonnenberg surface model including site-blocking terms, as implemented in the diffusion trapping code TESSIM-X. The paper will show that even in absence of surface saturation effects, the Waelbroeck-Ansatz is not able to describe the different limits arising from the Pick and Sonnenberg surface model. It will also show that the K r parameter used in the Waelbroeck-Ansatz is generally not a material parameter that can be used to model surface limits, but is a parameter specific to an experiment that should not be used for extrapolating laboratory results on retention and permeation to future fusion experiments, in particular if ion-implantation is involved.

Surface model description
In figure 1 the energy landscape of the Pick-Sonnenberg [7] surface model as implemented in TESSIM-X is shown. At the chemisorbed surface state, with concentration CA of chemisorbed HIs, the following fluxes in units of (m −2 s −1 ) are balanced in steady-state: the influx Γ Diss of diatomic HI gas with molecule mass m HI at pressure P HI and temperature T: (1) The out-flux Γ Rec of diatomic HI molecules after recombination in a Langmuir-Hinshelwood process (see e.g. [21]): The flux Γ Bulk of HI atoms from the chemisorbed state to the subsurface solute sites: To bulk transition rate s −1 .
The flux Γ Surf of HI atoms from the subsurface solute sites to the chemisorbed state: The areal density of host material surface atoms is computed from the number density of the host lattice density ρ(m −3 ) and the monolayer width . The concentration of chemisorbed HIs CA is the unit-less fraction of sites occupied per surface atom and ranges from 0 to η, with η the number of sites per surface atom. Similarly, the subsurface solute concentration CS0 is the unit-less fraction of sites occupied per bulk atom. Site blocking terms like (1 − CA/η) are only introduced for entering the chemisorbed surface state, no site blocking is introduced for entering the subsurface solute, since the diffusion trapping concept (see e.g. [22]) is derived and only valid for low site occupancy in the bulk.
At the subsurface solute location with concentration CS0 the following fluxes in units of (m −2 s −1 ) are balanced in steady-state: Γ Bulk (equation (3)), Γ Surf (equation (4)) and Γ Diff : In TESSIM-X the two flux balances are written as ordinary differential equations (ODEs) by invoking material conservation (equation (6)) at each of the two boundary surfaces: the left-or inlet-side and the right-or outlet-side.
In total these flux balances add four additional ODEs to be solved in TESSIM-X, two on each side: CS0 is the boundary value of the solute diffusion profile which is discretised on a grid in order to obtain finite difference approximations of the spatial derivatives needed in the bulk solute diffusion equation (see e.g. equation (1) in [23]). Thus, its cell width Δx in equation (6) must match the grid spacing at the surface Δx Surf to assure mass conservation in the numeric solution.
To derive simple algebraic steady-state results for CS0 from the above equations they are re-written without the site blocking terms and by collecting all non-concentration dependent terms into single pre-factors Ω X .
In steady-state in the absence of temperature gradients, the solute gradient is constant, i.e. the solute profile is linear and thus the diffusion flux can be written as The Ω Diff pre-factor for Γ Diff thus also contains the system dimension ΔX Grad .
In their paper Pick and Sonnenberg derive their 'recombination coefficient' K PS r from a surface flux balance on the inlet side of a gas-loading experiment using Waelbroeck's Ansatz: Γ Diss − 2 × K PS r (ρCS0) 2 = Γ Diff ≡ 0. They set Γ Diff ≡ 0 in steady-state, which is accurate for a gas-loading experiment when the entire sample volume is homogeneously filled. They then invoke Sieverts' law for solution of a diatomic molecule to obtain the expected steady-state value of CS0 = S(T) √ P and then solve for K PS r . Using equation (9) this yields for the effective recombination coefficient K PS r in Waelbroeck's Ansatz: From this derivation three prerequisites for the applicability of the Waelbroeck's Ansatz in combination Pick-Sonnenbergs 'recombination coefficient' can be deduced: firstly, the surface must always be in steady-state. Secondly, the surface concentration CA must always be small compared to η, i.e. the surface must not saturate. Thirdly, the diffusion flux from the surface into the bulk or vice versa must be small compared to the influx. While the first two prerequisites are often quoted in literature, the third prerequisite, Γ Diff ≡ 0, is often neglected. Also it should be noted that the derivation of equation (11) only describes the released flux of HIs off a surface that is required such that the net uptake of HIs from a gas phase influx fulfils Sieverts' law. In the original publication by Waelbroeck et al [5] the parabolic dependence of Γ WB Rec on CS0 stems from the assumption that in steady-state CA is directly proportional to CS0 (called 'surface-sub-surface equilibrium' SSE hypothesis). The applicability of the SSE hypothesis will be discussed in detail in section 3.2. However, looking at the derivation of equation (11) an alternative interpretation is that the parabolic dependence of the effective recombination flux Γ WB Rec in Waelbroeck's Ansatz on CS0 does not stem from Langmuir-Hinshelwood type recombination but from the pressure dependence in Sieverts' law. Thus the Waelbroeck Ansatz describes a required effective recombination flux Γ WB Rec to counter the influx from dissociation Γ Diss of an ambient gas phase to yield a steady-state value of CS0 as required by Sievert's law for the given ambient gas pressure. It does NOT describe the Langmuir-Hinshelwood type recombinative release from the surface due to a flux from the bulk to the surface, but this is exactly what it is mainly used for in boundary conditions in diffusion trapping calculations when ion beam/plasma loading is modelled. This fact will become even more apparent in the next section where the Waelbroeck Ansatz in combination with Pick-Sonnenbergs 'recombination coefficient' (from now on referred to as the WB-model) is compared with steady-state results obtained from equations (7)-(10) which describe the Pick-Sonnenberg model without site blocking (from now on referred to as the PS-model). The treatment is very similar to [8] with the exception that Sieverts' law will be enforced by the choices of rate parameters Ω X in equation (9) which allows to identify the similarities and discrepancies between the PS-and WB-model.
In steady-state, the surface can limit the uptake/release of HIs into/from the bulk only in the presence of sources or sinks in the bulk that lead to diffusive fluxes to or away from the surface respectively (i.e. Γ Diff = 0). Without such sources or sinks the surface follows Sieverts' law and is in equilibrium with the ambient gas atmosphere. The diffusion fluxes driven by the sources and sinks in steady-state are linear gradients pivoted by CS0 at the surface, thus, if the transport across the surface is not fast enough, they increase or deplete CS0. Therefore, the surface limits manifest themselves in values of CS0 different from the not-surface-limited case: for gas-loading the non-surface-limited value of CS0 is given by Sieverts' law CS0 ≡ S(T) × √ P. For bulk sources without ambient gas pressure (e.g. ion beam loading) the non-surface-limited value of CS0 = 0 and the HI release is only limited by diffusion. In the absence of strong sources or sinks in the bulk, the transport across the surface does not limit the emission/uptake of HIs and, therefore, there is no need for either the WBor PS-models in boundary conditions for diffusion trapping modelling and simple diffusion-limited transport or Sieverts' law suffices to describe the boundary value CS0.

Enforcing Sieverts' law
Sieverts' law describes the uptake of HIs into a surface in equilibrium with an ambient gas atmosphere at pressure P and temperature T. The uptake is described by the solute concentration as function of temperature and gas pressure as is a material constant that must not depend on the details of the energy landscape in figure 1 but only on the heat of solution ΔE Sol . As is also shown in Pick and Sonnerberg's paper [7], this means that the activation energies and pre-factors in the rate constants cannot be chosen arbitrarily. In the absence of diffusion fluxes to/from the surface, the steady-state solution of equation (12) for CS0 must yield Sieverts' law as in equation (13).
Since the values of the Ω X factors in equation (9) always have to be chosen such that equation (13) is fulfiled, equation (13) can be used to compute the value of one of the Ω X factors and thus reduce the number of free parameters of the surface model. The obvious choice is to compute the value of Ω Bulk from equation (13) since the transition from surface to bulk is most affected by the heat of solution as can be seen in figure 1. Also, Ω Bulk is experimentally hardly accessible whereas reasonable estimates for the other rate parameters Ω X can be derived (see sections 3.6 and 3.7). Therefore, for the remainder of the paper Ω Bulk is defined as in equation (14) unless specified otherwise.

Fixed diffusion flux from bulk to surface
To better show what the WB-and PS-model describe, it is instructive to investigate how they treat the presence of a fixed diffusion flux Γ Bulk Diff from or to the bulk across the left inlet surface. Of course in a more realistic case Γ Bulk Diff would vary as CS0 varies, since from equation (10) follows Γ Diff ∝ ΔCS, but for simplicity a constant Γ Bulk Diff is assumed here in this initial example. For Γ Bulk Diff > 0 (flux to the right) it acts as a sink and for Γ Bulk Diff < 0 (flux to the left) it acts as a source for HIs at the surface. The pressure on the inlet side is P HI . In the PS-model, the flux balance at the surface in steady-state contains two equations for the two unknowns CS0 and CA as in equation (15) 0 In the WB-model, the flux balance contains a single equation for a single variable CS0 as in equation (16) 0 Solving equation (15) for CS0 and CA yields: The result for CA in equation (17) shows that the sum of all fluxes to or away from the surface (Γ Bulk Diff from/to the bulk, Γ Diss from the gas phase) must be balanced via recombinative losses from the chemisorbed surface states. The result for CS0 PS in equation (17) shows that CS0 is determined by the net transport of HIs from the bulk (Γ Bulk Diff /Ω Surf ) and from the chemisorbed surface states (Γ Bulk /Ω Surf ) to the subsurface solute sites. Since, as explained above, surface limits manifest themselves in the value of CS0, this means that surface limits are due to these two net transport contributions. Either can limit the transport across the surface and, therefore, both need to be taken into account in a surface model.
By enforcing Sieverts' law one removes the dependence on Ω Rec in the solution for CS0 in equation (18). Solving equation (16) for CS0 and using the definition of K PS r from equation (11) yields: The result for CS0 for the PS-model after enforcing Sieverts' law allows to see the similarities between the PS-and the WB-model: the results for both models in equations (18) and (20) share the S(T) Diff Ω Diss + P HI term. By using the solution for CA from equation (17) together with equation (14) in Γ Bulk from equation (9), this term can again be identified as Γ Bulk /Ω Surf . As before this term describes the net transition of HIs from the chemisorbed surface states into the subsurface solute sites. However, the WB solution lacks the other contribution Γ Bulk Diff /Ω Surf which is the net transition of HIs from the bulk to the subsurface solute sites. This difference shows that the WB-model does not describe the release from the bulk, but only describes the transition from surface to bulk and thus is only equivalent to the PS-model if the transition from bulk to surface controlled by Ω Surf is not the rate-limiting process. This raises the concern whether it makes sense to model the release of HIs from the bulk (e.g. during a thermal desorption experiment) for surface limited systems with the WB approach using some 'custom' choice of K r since it cannot describe all processes contributing to a surface limit properly. Also using a 'custom' choice of K r to reproduce CS0 PS , by equating equations (18) and (19), means that K r would have to become dependent on Γ Bulk Diff . This means that custom choices of K r are not material parameters but experiment specific parameters that cannot generally be used for predictions. In the work by Waelbroeck in [4,5] the conditions for cases where the release from the surface can be described by K r CS0 2 are discussed at length (see sections IX and X in [5]). The conclusion is that an equilibrium between CS0 and CA must exists such that CS0 ∝ CA which is called the SSE hypothesis. From equation (17) it is evident that this SSE hypothesis is equivalent to Γ Bulk Diff /Ω Surf ≈ 0. This means that either the diffusion flux Γ Bulk Diff is zero or that the transition rate constant Ω Surf into the chemisorbed surface state is very large, i.e. no barrier ΔE Surf to enter the surface exists. The fact that the SSE hypothesis breaks down for large values of ΔE Surf > ΔE Diff is also stated in section IX in [5], however, even in the absence of a barrier, Ω Surf will still be finite and not infinite. Thus the SSE hypothesis breaks down also for large values of Γ Bulk Diff which was not considered at that time. This again shows the important prerequisite: Γ Diff ≈ 0 for applying the Waelbroeck Ansatz as pointed out in the introduction. It should be noted that ignoring Γ Diff in the PS-model flux balance equation (15) leads to results for CS0 and CA that support the SSE hypothesis independent on Ω Surf and thus suggest an equivalence of the PS-and WBmodel (see e.g. derivation of equations (5) and (6) in [24]). As the following sections will show the PS-and WB-model are generally not equivalent for Γ Diff = 0.
In the next examples the models will be compared for more realistic systems with different surface limits to map out where the WB-model can actually be applied and where the more advanced PS-model is required.

Permeation: outlet side
The steady-state permeation flux on the outlet side of an ion or gas-driven permeation experiment can be described by a fixed source solute concentration on the inlet side CS Src that results in a diffusion flux towards the outlet side due to the gradient CS0−CS Src ΔX , where ΔX is the thickness of the sample. The balance equations for the outlet side are the same as in equations (15) and (16) with Γ Diss = 0 since P HI = 0 on the outlet side and that the diffusion flux now properly depends on CS0 as Γ Bulk Diff = −Ω Diff × CS0 − CS Src . Solving equation (15) under these conditions for CS0 now yields the result in equation (21) Similarly, solving equation (16) with Γ Diss = 0 and the CS0 dependent diffusion flux yields equation (22) with the same definition of β as in equation (21) CS0 WB = β The first interesting feature of this comparison is that the PSresult in equation (21) for lim α→∞ becomes identical to the WB result in equation (22). This again highlights the fact that the WB-model does not describe the transition from the bulk to the surface correctly, it is only applicable when the transition to the surface is much faster than diffusion to the surface: Ω Surf Ω Diff , i.e. the surface transition is not the rate-limiting step.
In a system with a diffusion-limited outlet side, the release of solute HIs at the surface is fast enough such that no solute accumulates and CS0 ≡ 0. For a given CS Src , this yields the highest permeating flux. If the release through the outlet side is not fast enough, the solute piles up at the subsurface solute location resulting in 0 CS0 CS Src . Therefore, the ratio of = CS0/CS Src is a measure of how strong the outlet side limits the permeation flux. Using the PS-result in equation (21) one obtains equation (23) for this concentration ratio .
For given material parameters α and β the critical source concentration strength CS Src Crit , defined in equation (24), controls whether the permeation is flux is limited ( ≈ 1) by the outlet surface or not ( ≈ 0). For Θ 1, becomes ≈ 1 whereas for Θ 1, becomes ≈ 1/ (1 + α). So the system is not surface limited only if both, a strong source and a large value of α 1 are present, i.e. when the transition from CS0 to CA is not rate-limiting and the solute does not pile up at the outlet side.
Performing the same analysis with the WB result for CS0 in equation (22) one obtains equation (25) for the concentration ratio = CS0/CS Src Again the results for WB-equation (25) and PS-model equation (23) become identical for lim α→∞ . But in contrast to the PS-model, the WB-model always suggest a non-surfacelimited system for strong sources Θ * 1 since it cannot model surface limits due to finite values of α.
While is a mathematically convenient measure for the presence of surface limits, it is experimentally hardly accessible. Therefore, the expression for the solute concentration gradient ΔCS = CS0 − CS Src , driving the measurable permeation flux Γ Perm = −Ω Diff ΔCS, as function of Θ is given for the WB-in equation (27) and the PS-model in equation (26).
Both models yield the same scaling of Γ Perm with the source strength: for the surface limited case Θ 1, Γ Perm is ∝ Θ 2 . For Θ 1 both models yield Γ Perm ∝ Θ, but only the WBmodel then suggests a non-surface-limited system, i.e. ≈ 0.
In the PS-model a linear scaling of Γ Perm with Θ only means a non-surface-limited system, if also α 1. Otherwise, the system is surface limited by the pile up of solute in CS0 on the outlet side due to the small transition rate from CS0 to CA.
As long as α is large the WB-and PS-models yield very similar results. Assuming that the transition from the subsurface solute to the chemisorbed surface state is comparable to a diffusive jump of one monolayer width a 0 , Ω Surf can be calculated as in equation (28) Ω Surf from equation (28) together with the definition for Ω Diff in equation (10) yields α = 6ΔX/a 0 . Since the typical sample thickness in a permeation experiment ΔX (μm to mm) is a 0 (10 −10 m), α usually is very large (10 4 to 10 7 ). Thus, the WB and PS-model should be indistinguishable on the outlet side unless a massive barrier exists that blocks the transition to the surface, bridging the 4-7 orders of magnitude between Ω Surf and Ω Diff . However, if the gradient length ΔX in Ω Diff becomes small, e.g. when CS Src is located close to the surface, then α becomes small and the two models should diverge significantly. This will be investigated in the next sub-section where the WB-and PS-model are compared at the inlet side of an ion-driven permeation experiment.

Inlet side, ion-driven permeation
In figure 2 the solute profile on the inlet side of an ion-driven permeation experiment in steady-state is shown following the concept in [16,17].
As explained above the gradient length (i.e. the sample thickness ΔX) driving the permeation flux Γ Perm is much larger than the mean projected range R P where most of the implanted ions are deposited. Therefore, Γ Perm Γ SURF Diff and it can be neglected in the flux balance at the surface. This means that in steady-state the diffusion flux Γ SURF Diff at depth R P must be equal to the influx from the ion beam Γ Ion which defines C Max in figure 2 as in equation (29).
The rest of the flux balance (equation (30)) at the surface for the PS-model is equivalent to equation (15) with P HI = 0 on the in-and outlet side, such that the uptake from the gas phase is negligible (Γ Diss = 0).
Similarly, the flux balance for the inlet side of ion-driven permeation applying the WB-model yields equation (31) 0 Assuming a diffusion-limited outlet side (zero solute concentration) the permeation flux is given by Γ Perm = Ω Diff C Max . Solving equations (30) and (31) for CS0 and using the result for C Max from equation (29) the permeation flux can be written as equation (32) for the PS-and equation (33) for the WB-model.
Comparing the result in equation (32) with the temperature independent result for diffusion-limited ion-driven permeation Γ DL Perm in equation (34) (i.e. CS0 ≡ 0 at in-and outlet side), one finds that for ω + λ 1 √ Γ Ion 1 the system becomes surface limited and Γ Perm is larger than in the diffusion-limited case. For ω 1 the system is surface limited independent on the ion flux, whereas for ω < 1 a critical ion flux, Γ Crit Ion , can be defined as in equation (35) below which the system becomes surface limited: The ω parameter is the inverse of the α parameter in the previous section dealing with the outlet side. Small values of ω mean that the rate constant for transition from the subsurface solute to the chemisorbed surface state is faster than rate constant for diffusion from the implantation depth R P . For small values of ω the transition from the subsurface solute to the chemisorbed state is not rate-limiting and again the WB-model becomes identical to the PS-model. For ω < 1 and Γ Ion Γ Crit Ion , Γ Perm scales as ∝ √ Γ Ion for ω 1 and/or Γ Ion Γ Crit Ion , Γ Perm scales as ∝ Γ Ion . Therefore, a linear scaling of Γ Perm with Γ Ion does not mean the system is not surface limited, it might be limited due to ω 1 which can only be understood within the PS-model.
Estimating values for the ω parameter using the same assumptions for Ω Surf in equation (28) as before for the α parameter yields ω ≈ a 0 / (6R P ) which for typical ion energies ≈ 100 (eV) yields ω ≈ 0.01. This means that small differences in the activation energies for diffusion ΔE Diff and the transition to the surface ΔE Surf can result in ω 1 depending on temperature and, thus, introduce a surface limit that cannot be described by the WB-model. This will be demonstrated in the numerical examples in sections 3.6 and 3.7.

Inlet side, gas-driven permeation
For the inlet side of a gas-driven permeation experiment the flux balance is given in equation (36) for the PS-and in equation (37) for the WB-model The inlet side is considered surface limited if CS0 is less than the value predicted by Sieverts' law S(T) √ P because the diffusion flux from the surface into the bulk is strong and the transition from the gas phase to the subsurface solute is too weak. Solving equation (36) and in equation (37) and using the result to calculate the ratio CS0/ S(T) √ P yields the results in equation (38) for PS-and in equation (39) for the WB-model.
Again the models in equations (38) and (39) become identical for slow diffusion and/or fast transition from subsurface solute to the chemisorbed surface state, i.e. for γ ≈ 0.
Assuming zero pressure on the outlet side of the gas-driven permeation system and, thus, no solute on the outlet side the permeation flux is ∝ CS0 yielding equation (40) for PS-and in equation (41) for the WB-model.
From equations (38) and (40) for the PS-model a critical pressure P Crit = ξ/(γ + 1) 2 can be defined. For pressures P P Crit , Γ Perm ∝ P and CS0/ S(T) √ P PS,WB < 1, i.e. the system is surface limited on the inlet side. This linear dependence of Γ Perm on pressure is a commonly accepted sign of surface-limited permeation in interpreting experiments using the WB-model [19]. Usually experimentalists increase the pressure until Γ Perm ∝ √ P and then the WB-model suggests that the experiment is not surface limited anymore. Indeed for P P Crit both the WB and PS-model yield Γ Perm ∝ √ P, but only the WB-model yields CS0/S(T) √ P WB ≈ 1, i.e. suggesting that the system not surface limited. For the PS-model for P P Crit the value of CS0/S(T) √ P PS ≈ 1/ (1 + γ) depends on the value of γ. If γ 1 then CS0/S(T) √ P PS < 1 even for P P Crit suggesting that despite the ∝ √ P scaling of Γ Perm the system might still be surface limited.
Since γ is the inverse of the α parameter derived for the outlet side in section 3.3 an estimate for the numerical value of γ can be derived. Based on the estimates for α of (10 4 to 10 7 ), in the absence of the surface barrier (i.e. ΔE Surf ≈ ΔE Diff ), the expected value for γ is very small and thus for non-surfacelimited system both the WB and the PS-model should again be indistinguishable.

Example: ion-driven permeation through tungsten
So far a comparison of the PS-with the WB-model has shown that they are only equivalent if the transition from subsurface solute sites to the chemisorbed surface state is not the rate-limiting process. The question remains for which cases this transition is actually critical, i.e. when is the WB-model not applicable to describe the experimental conditions properly. For the different sample sides and load scenarios three parameters were derived and values thereof were defined which determine if the transition from subsurface solute sites to the chemisorbed surface state is the rate-limiting process: outlet side α ≈ 0, ion-driven inlet side ω 1 and gasdriven inlet side γ 1. From the foregoing discussion on the values of α and γ the outlet side and the gas-driven inlet side should generally be well described by the WBmodel except for very strong transport barriers with ΔE Surf > ΔE Diff . However, for ion-driven experiments, due to the steep near surface gradients, the condition ω > 1 cannot be generally excluded. To determine numerical values for the ω and λ parameters in equations (32) and (33) for ion-beam loading of pure W, literature values for the solubility and solute diffusion coefficient from [25,26] are used. These two references suggest essentially the same solubility, but they differ in the activation energies for diffusion ΔE Diff with values of 0.27(eV) and 0.39(eV), respectively. In this numerical example D(T) = D 0 exp −ΔE Diff /K B T is taken from [25] and Ω Surf = δ 0 6D 0 a 2 0 × exp −ΔE Surf /K B T is calculated following equation (28). The value of ΔE Surf is thereby varied from ΔE Surf = ΔE Diff from [25] to ΔE Surf = 0.5(eV) which is roughly a factor of two larger than what [25] suggests for ΔE Diff . From equation (10) follows Ω SURF Diff = ρD(T) R P . Assuming implantation of D at 100 (eV) and impact along the surface normal, a mean projected range R P ≈ 3.8(nm) is calculated by SDTrimSP [27,28]. Following equation (9), the value for Ω Diss is the product of the gas kinetic factor Φ(m HI T) and the sticking probability Ξ(T ). Since it is generally assumed [7] that sticking on clean surfaces is high (i.e. Ξ ≈ 1) Ω Diss = 2Φ(m HI T) is used. Based on these assumptions the following temperature-dependent expressions for ω and λ are obtained: To investigates the influence of a rate-limiting step from the subsurface solute to the chemisorbed surface state controlled by ω, the ratio of permeation fluxes in equations (32) and (33) to the permeation flux for diffusion-limited conditions equation (34) is calculated leading to equation (44).
In figure 3 the temperature evolution of ω and λ for W is shown for different values of ΔE Surf . For low temperatures and high ( 0.4(eV)) values of ΔE Surf follows ω 1 and, thus, the release out of the inlet surface is limited in the PS-model which cannot be modelled in the WB approach. The value for λ does not depend on ΔE Surf since under the assumptions above its temperature dependence stems mainly from the solubility S(T ) and to a lesser extent also from the gas kinetic factor.
In figures 4 and 5 the ratio of the WB-and PS-model based permeation flux to the diffusion-limited permeation flux from equation (44) is compared for different values of ΔE Surf and Γ Ion , respectively.
In both figures the WB-model cannot describe the lowtemperature surface limit due to a finite transition rate from the subsurface to the surface controlled by the ω parameter. This leads to an underestimation of the permeation flux at low temperatures in figure 4 by a factor of 3 for ΔE Surf = 0.4(eV) and more than an order of magnitude for ΔE Surf = 0.5(eV) compared with the PS-model. At very high temperatures 1000 K, when the solubility increases, both models predict a strong surface limit due the increased back transition Γ Bulk from the chemisorbed surface state controlled by the λ parameter. As can be seen in figure 5 this surface limit is strongly dependent on the ion flux Γ Ion at high temperatures: The higher Γ Ion , the more the diffusion flux (which needs to compensate it and is, thus, equal to Γ Ion ) from the ion deposition location  R P to the subsurface dominates over the flux Γ Bulk from the surface: thus, the higher Γ Ion the lower the influence of the surface on the permeation flux at these high temperatures. At low temperature λ is too small to lead to a flux dependence in equation (44).
It should be noted that these calculations are based on the algebraic results and therefore do not include the effect of surface saturation which might affect the high flux cases. However, at the high temperatures of 1000-2000 (K), where this surface limit appears, saturation of the chemisorbed state is unlikely even for the highest flux of 10 24 (m −2 s −1 ) considered here.

Example: measurement of K r for tungsten
An often cited publication describing the measurement of the 'recombination rate' coefficient K r for use in the WB-model is by Anderl et al [29]. In this ion-driven permeation experiment the recombination coefficient was calculated from the measured permeation and re-emission flux using diffusion and solubility data from [26]. The equation number 2 in [29] used to calculate K r follows from the flux balance in equation (45).
Solving the flux balance in equation (45) for C Max , CS0 and K r yields equation (46), identical to the result in [29] In the derivation of equation (45) the re-emission flux is assumed to be equal to the incident ion flux in steady-state. This is justified, as the permeation flux is orders of magnitude smaller, even for the thin 25 μm foil used in the experiment in [29]. In [29] the re-emission flux was measured, but the values are not given, therefore, for the comparison here, this assumption needs to be made. Using the result for the Γ Perm for the PS-and WB-models from equations (32) and (33) allows to calculate the expected value of K Anderl r for the two models which can be compared with the experimental result from [29].
Using the WB Ansatz in the equation for K Anderl r leads finally to equation (48) which is identical the Pick-Sonnenberg version of the 'recombination coefficient' from equation (11), as expected. In contrast, using the PS-model leads to a fluxdependent expression for K Anderl r in equation (47) for non zero values of ω. This again highlights the fact that Waelbroeck's 'recombination coefficient' is not a material parameter, but depends on the experimental conditions due to its dependence on Γ Ion in the presence of surface barriers controlled by the ω parameter.
In figures 6 and 7 the different solutions for K Anderl r are compared for different values of ΔE Surf and Γ Ion respectively, also shown is the fit to the experimental data K Lit r = K 0 exp −E r /K B T value using the 'all combined' values for K 0 = 3.2 × 10 −15 (m 4 s −1 ) and E r = 1.16(eV) from table 2 in [29]. In both figures 6 and 7 the literature result K Lit r is in fundamental disagreement with the WB-model based value  K WB r : while K WB r decreases with temperature, K Lit r increases with temperature. However, comparing K Lit r with K PS r for the different values of ΔE Surf in figure 6 and different values of Γ Ion in figure 7, a match can be found. This suggests that the discrepancy between the WB-model based K Anderl r with the experiment could be explained by a finite transition rate from subsurface solute to the chemisorbed surface state due to finite, non-zero values of ω.
An even better fit of K Lit r could be found if in Ω Surf not just ΔE Surf but also the pre-factor ν BS thus −E r = −1.16(eV) = 2ΔE Sol − ΔE Diss . Using ΔE Sol ≈ 1.0(eV) from [26] this yields a required barrier for adsorption of ΔE Diss ≈ 3(eV). From the discussion in [7,24] the general assumption is that for clean surfaces ΔE Diss ≈ 0 which is supported by [30]. Therefore, such a high barrier for adsorption/sticking seems unlikely, but could be due to impurities. On the ion bombarded inlet side steady sputter cleaning of the surface should suppress impurities, but still their influence cannot be ruled out.

Comparison to TESSIM-X calculations
To investigate the effect of surface saturation in an iondriven permeation experiment the full version of the Pick-Sonnenberg model, including site blocking terms (see equations (1)-(4)), as implemented in TESSIM-X, is used. The permeation through a 30 μm thick W foil is modelled using the same choice of experimental and material parameters as in the numerical example in section 3.6. The transition from subsurface solute to the chemisorbed surface state uses the same activation energy as the solute diffusion coefficient (ΔE Surf = ΔE Diff ) and is, thus, not a rate-limiting factor, still this transition rate is finite in the PS-model. Surface saturation in the Pick-Sonnenberg model can only occur when the concentration CA of HIs in the chemisorbed surface state approached its maximum value η, the number of surface sites per W atom. For simplicity we assume η ≡ 1 and choose the pre-factor χ 0 and activation energy ΔE Rec for the Langmuir-Hinshelwood recombination coefficient χ(T ) in equation (2) based on the work in [21]. In the simulation the W foil is exposed to a linear D-ion flux ramp from 10 16 to 10 23 (m −2 s −1 ) at constant temperature of 600 K. The high temperature is used such that the system is in steady-state at all times and can be directly compared to the algebraic results from the PS-and WB-model. The activation energy ΔE Rec is varied from 0.6 to 1.2(eV), the same range of values as given in [21] for different surface states on W. For the TESSIM-X calculation to be identical to the algebraic treatment in section 3.4, only the inlet side is modelled using the Pick-Sonnenberg model, the outlet side is modelled as a diffusion-limited surface leading to CS0 = 0 = CA at the outlet side.
In figure 8 the permeation flux during the linear Γ Ion ramp is shown for different choices of ΔE Rec . Since the algebraic results for CS0 and thus for Γ Perm do not depend on Ω Rec in the PS-or WB-model, they do not vary with ΔE Rec , they only vary with Γ Ion . As long as Γ Ion is low and/or ΔE Rec is low, the surface does not saturate and the full Pick-Sonnenberg model in the TESSIM-X calculation perfectly matches the algebraic results from the PS-and WB-model. However, for ΔE Rec = 1.2(eV) and Γ Ion ≈ 10 22 (m −2 s −1 ) the TESSIM-X permeation flux increases drastically, because the inlet side saturates, blocking the uptake of HIs from the bulk and thus the solute diffusing towards the surface piles up leading to an increase of CS0. As, according to equation (29), the permeation flux is proportional to C Max , which is pivoted by CS0 and thus increases as CS0 increases, Γ Perm increases (see section 3.4).
In figure 9 the subsurface solute concentration during the linear Γ Ion ramp is shown for different choices of ΔE Rec . Again, as before for Γ Perm , the results for the PS-or WB-model they do not vary with ΔE Rec but with Γ Ion only. For nonsaturated surfaces the PS-model perfectly matches the CS0 result from TESSIM-X and as the surface starts approaches saturation for ΔE Rec 1(eV), the PS-model starts to deviate slightly, as expected. In contrast the WB-model is totally off, basically predicting CS0 ≈ 0 under all conditions. The reason is again that the WB-model cannot describe the finite transition rate from subsurface solute to chemisorbed surface, even when it is not limiting the permeation rate. However, in the present example the C Max value driving the permeation flux, is much larger than CS0 because it is dominated by the Γ Ion /Ω SURF Diff term in equation (29), resulting in CS0 C Max at the inlet side and the release of HIs being limited by diffusion. Thus, the failure of the WB-model to correctly calculate CS0 remains hidden when only looking at Γ Perm . This is the reason why the failure of the WB-model is generally not detected: during ion beam loading the finite transition rate across the surfaces can be hidden by the strong ion-implantation bulk source, unless the surface saturates and blocks or when ΔE Surf > ΔE Diff . This is also visible in both the PS-model and the TESSIM-X result: for ΔE Surf = ΔE Diff and non-saturating surfaces, they also yield CS0 C Max , albeit at CS0 levels much higher than the WB-model. This also means that, without saturation/blocking and ΔE Surf ΔE Diff , surfaces react essentially like diffusionlimited surfaces to the diffusion out-flux of implanted HIs. Therefore, neither the PS-nor the WB-model are often needed to describe the inlet side of ion-driven experiments in diffusion trapping calculations.
The solution of the PS-model for the different loading conditions also produces an expression for the chemisorbed In figure 10 the chemisorbed surface concentration CA as calculated by TESSIM-X and the PS-model are compared. Again the PS-model perfectly matches the TESSIM-X result until CA approaches η, then the TESSIM-X solution yields the expected saturation values whereas the PS-model keeps increasing.
In the absence of a barrier from the subsurface solute to the chemisorbed surface (ΔE Surf ΔE Diff ) and when the surface does not saturate both the PS-and WB-model reproduce the TESSIM-X permeation flux results in this example calculation. While the PS-model actually matches the CA and CS0 underlying the observable permeation flux, the WB-model matches the permeation flux simply because its wrong result for CS0 does not matter under these conditions. While the presented PS-model cannot handle saturation effects, the solution it provides for CA allows to detect surface saturation and thus it can be used to detect the limits of its applicability.
It should be noted that Sievert's law is also enforced in the TESSIM-X calculations, thus, changing ΔE Rec also affects ΔE Bulk in equation (3) such that it compensates the effect of the significant changes in CA in figure 10 on the value of CS0 in figure 9. The fact that CS0 1 and hardly changes with ΔE Rec is the reason why the permeation flux in figure 8 is also hardly affected by ΔE Rec unless saturation occurs. Without saturation or surface blocking recombination cannot easily affect permeation because increasing ΔE Rec also means increasing ΔE Bulk (see equation (14)) to maintain Sievert's law, thus partly compensating the increase in Γ Bulk from the higher values of CA. As shown by the algebraic model results in figure 5 only at high temperatures where the increase in ΔE Bulk is compensated by temperature, the reduced recombination leads to a significant increase of both CA and CS0 and thus affects permeation.

Conclusion
The comparison of the Pick-Sonnenberg model with the Waelbroeck-Ansatz has shown that in the absence of surface barriers for the transition from subsurface solute to chemisorbed surface sites and/or low diffusion fluxes from the bulk towards the surface, the two approaches lead to identical permeation flux results. This result is in line with the work by Waelbroeck in [5]. However, the presence of impurities at the surface might introduce surface barriers that reduce the finite transition rate from CS0 to CA even further. Even in the absence of impurities, the break of crystal symmetry at the surface should lead to a barrier ΔE Surf > ΔE Diff as is suggested in recent DFT calculations in [31].
In general, the Waelbroeck-Ansatz only describes the uptake from the surface into the bulk and, therefore, it can only describe surface limits due to either a too weak or too strong transport from surface to bulk: for gas-loading the transport from surface to bulk is too weak if the supply of HIs from the gas phase across the surface cannot compensate the permeation flux. In contrast, for HI volume sources, from ion beam loading or de-trapping during TDS, the transport back from surface to bulk is too strong, if it leads to pile up of solute in CS0 due to limited release via recombination from CA. As shown in the numerical example in section 3.6 this back-reaction surface limit is only important at high temperatures where S(T ) is large. What the Waelbroeck-Ansatz cannot describe is the effect of high diffusion fluxes from the bulk towards the surface in combination with the finite transition rate from CS0 to CA and this is commonly neglected when using it in boundary conditions for diffusion trapping calculations. In particular for ion implantation, from a plasma or a laboratory ion beam, at the moderate temperatures commonly encountered in experiments, strong diffusion fluxes towards the surface exist, due to the close proximity of the mean ion implantation range to the surface. These high diffusion fluxes can then lead to surface limits that cannot be described with the Waelbroeck-Ansatz as was show in the example calculations for ion-beam-driven permeation in W and the determination of K r following Anderl's work in [29].
Trying to model such limits using the Waelbroeck-Ansatz leads to an ion-flux dependence of K r . Therefore, databases of K r values that have been compiled need to be treated with caution when used in ion-beam-loading experiments: K r is not a material parameter but an experiment-dependent parameter and must not be used to extrapolate to the higher ion fluxes expected in fusion devices compared with laboratory experiments. Generally, the Waelbroeck-Ansatz with values of K r measured from gas-driven permeation can only be applied to gas-phase loading, which is what is was derived for. Also, it should be noted that in diffusion-trapping modelling one should first try to use non-surface-limited boundary conditions (Sieverts' law or diffusion-limited transport) before introducing additional input parameters that reduce the likelihood of the solution given the experimental data in light of Occam's Razor. This is particularly true for the inlet side of ion-driven experiments. In table 1 the main results w.r.t. to the applicability of the Waelbroeck-Ansatz are summarised for the different loading conditions.
The comparison of the full version of the Pick-Sonnenberg model in TESSIM-X, including blocking terms for CA, with the simple algebraic PS-and WB-models shows that at very high fluxes and high barriers for the Langmuir-Hinshelwood recombination rate, the surface can saturate, making the simple models no longer applicable. For non-saturated surfaces the algebraic PS-model perfectly matches the numeric TESSIM-X results for CS0 and CA whereas the WB-model fails even for non-limiting but finite transition rates from subsurface solute to chemisorbed surface states. The comparison has also shown that this failure of the WB-model to correctly treat the transition of diffusion fluxes from the bulk across the surface is often not recognised because the pile-up of CS0 due to the finite transition rate is dwarfed by the concentrations introduced by the implantation source. This makes the WB-model equivalent to diffusion-limited modelling of the inlet side of ion-driven permeation experiments unless the surface is limited by the strong back-reaction from the chemisorbed surface state. However, as stated before, these back-reaction surface limits require high S(T ) values which for, e.g. W requires very high temperatures (see figure 5) exceeding typical values, at least from laboratory experiments.
The simplified algebraic steady-state version of the Pick-Sonnenberg model using the fluxes from equation (9) in this paper is similar to that in [8] with the exception that Sieverts' law is intrinsically enforced in all results which allows to see the similarities between the PS-and WBmodel. Also, this interestingly removes any dependence on Ω Rec from the solutions for CS0 or the permeation flux Γ Perm through the sample. This makes this implementation of the Table 1. Summary of load cases and conditions where the Waelbroeck-Ansatz is applicable and when surface limits are not present. The definitions of the Ω X parameters are given in equations (9) and (10) and further definitions are given in equations (1)-(4).

Case
Critical source strength WB-applicable for NOT surface limited for Pick-Sonnenberg model an interesting alternative to the Waelbroeck-Ansatz since it also only involves CS0 and only requires an estimate for Ω Surf as a single additional parameter. The other input parameters, S(T ) and Ω Diss , are also needed for the Waelbroeck-Ansatz with the Pick-Sonnenberg model for K PS r . Thus, using this model would not add numerical complexity or introduce additional solution variables, but would be a more general approach to handling surface limits in diffusion-trapping calculations.