Surface recombination in Pyrex in oxygen DC glow discharges: mesoscopic modelling and comparison with experiments

Surface recombination in an oxygen DC glow discharge in a Pyrex (borosilicate glass) tube is studied via mesoscopic modelling and comparison with measurements of recombination probability. A total of 106 experimental conditions are assessed, with discharge current varying between 10 and 40 mA, pressure values ranging between 0.75 and 10 Torr, and fixed outer wall temperatures ( Tw ) of −20, 5, 25 and 50∘C . The model includes O+O and O+O2 surface recombination reactions and a Tw dependent desorption frequency. The model is validated for all the 106 studied conditions and intends to have predictive capabilities. The analysis of the simulation results highlights that for Tw=−20∘C and Tw=5∘C the dominant recombination mechanisms involve physisorbed oxygen atoms ( OF ) in Langmuir–Hinshelwood (L-H) recombination OF+OF and in Eley–Rideal (E-R) recombination O2+OF , while for Tw=25∘C and Tw=50∘C processes involving chemisorbed oxygen atoms ( OS ) in E-R O+OS and L-H OF+OS also play a relevant role. A discussion is taken on the relevant recombination mechanisms and on ozone wall production, with relevance for higher pressure regimes.


Introduction
The interaction between low-temperature plasmas and surfaces is of paramount importance in various scientific and technological domains.Understanding it is crucial to advance materials science, biomedicine and nanotechnology, potentially helping to optimize plasma-based processes such as plasma etching, surface modification, thin film deposition, plasma-tissue treatment and atmospheric entry.In essence, delving into the interdisciplinary field of plasma-surface interactions not only enhances our fundamental understanding of the interdependence between two media in direct contact but also paves the way for ground-breaking technological innovations with far-reaching implications.
Low-temperature plasmas are known to modify solid surfaces by providing charges and reactive and energetic gaseous species (Oehrlein 1997, Kersten et al 2001, Vanraes et al 2021).However, the densities of active gaseous species in the plasma are also affected by surface kinetics.As such, a detailed understanding of this kinetics is often required to control plasma properties and design plasma processes (Kushner 2009, Murphy and Park 2017, Bogaerts et al 2022).Indeed, self-consistent macroscopic models, either fluid or kinetic, require an accurate description of the interaction of the plasma with the surface.However, the surface interaction data strongly depend on working conditions and are often derived from fitting experimental measurements obtained for a particular condition, which limits the predictive capabilities of models dependent on plasma-surface interactions.
Oxygen plasmas are widely used in material processing applications and in biomedical applications at high pressure.In particular, atomic oxygen is a key reactive species for many plasma processing applications.It is also a main species of interest for product separation in plasmas for electrochemical conversion at intermediate and high pressures (Rohnke et al 2004, Chen et al 2020, Pandiyan et al 2022).The main mechanism of production of O atoms in pure oxygen plasmas is electron impact dissociation of O 2 .At pressures below tens of Torr, three-body collisions in the plasma volume have low probability and thus O atoms are predominantly lost by surface recombination back into O 2 (Gousset et al 1987, Macko et al 2004, Annušová et al 2018, Dias et al 2023).The probability for this latter process, γ O , defined as the ratio between the number of atoms that successfully recombine on the surface and all the atoms that strike the wall, is highly dependent on the type of surface and on working conditions of pressure, current, gas temperature and wall temperature.Hence, atomic oxygen surface recombination is a very complex process and has been widely studied in plasma reactors for different conditions and surfaces.These studies are listed and compared in the recent review by Paul et al (2023).The recombination is very effective (above 10%) on metallic and catalytic surfaces (Morgan and Schiff 1964, Mozetic and Cvelbar 2007, Mozetic et al 2015), less effective on crystals and semi-catalytic surfaces and even less (below 1%) on inert materials like silica-based or glass-based materials such as Pyrex and quartz.These materials are common in plasma reactors and several studies have measured the dependence of O recombination on surface temperature and other parameters for Pyrex (Pagnon et al 1995, Macko et al 2004, Booth et al 2019, 2020), for quartz (Kim and Boudart 1991, Balat-Pichelin et al 2003, Cartry et al 2006, Lopaev et al 2011a) and similar materials (Rakhimova et al 2009, Krištof et al 2012).Moreover, atomic oxygen recombination in silica-like materials is of paramount importance in defining the heat transfer to the surface of space vehicles in atmospheric re-entry conditions (Kovalev et al 2005, Kovalev and Kolesnikov 2005, Bourdon and Bultel 2008).
The same recombination process can be studied theoretically and numerically (Cacciatore et al 1999, Rutigliano andCacciatore 2016).Plasma-surface interactions take place in multiple time scales.The shortest time-scale, of the order of fs, is associated with the vibrational motion of adsorbed species (Graves andBrault 2009, Neyts andBrault 2017).The description of elementary steps of adsorption, desorption, diffusion, and reactions on the surface on the atomic level occurs in time-scales in the range ns-ms (Guerra and Marinov 2016).However, experiments often involve an evolving surface dynamics during several minutes.Mesoscopic models are a solution to bridge the gap between the detailed atomistic description and the macroscopic scales of real systems (Alves et al 2018).These are coarse-grained models describing surfaces and their occupation by surface site densities and surface species densities.In these models, the system evolution is described by a series of reactions with prescribed rates.Mesoscopic models can adopt a stochastic (Guerra and Loureiro 2004, Guerra and Marinov 2016, Marinov et al 2017) or a deterministic description of the system, whose mathematical resolution can be either analytical (Guerra 2007, Rakhimova et al 2009, Lopaev et al 2011a, Booth et al 2019) or numerical (Gordiets et al 1996, Gordiets and Ferreira 1998, Cartry et al 2000).They require as input the energy barriers for each elementary step and other physical parameters that depend on the atomistic interactions described by more microscopic models, and thus are not fully self-consistent.However, mesoscopic models can effectively provide an interesting physical insight of the surface kinetics processes, receiving input from the gas phase chemistry and be used to describe both the steady-state and the dynamics of the system up to any time-scale.
Despite the wide variety of models and experimental conditions, there is a widely accepted view of atom surface recombination processes (Seward andJumper 1991, Cacciatore andRutigliano 2009).It can generally occur via two mechanisms: Eley-Rideal (E-R) processes where a particle from the gas phase directly recombines with an adsorbed (chemisorbed or physisorbed) particle, which are thus dependent on the adsorbed density to the first order; Langmuir-Hinshelwood (L-H) processes where a physisorbed particle diffuses along the surface and recombines with another adsorbed particle, that are dependent on the adsorbed density to the second order.In the case of oxygen plasmas, if E-R and L-H processes involve one oxygen atom O and one oxygen molecule O 2 , they can give origin to ozone O 3 , another interesting species for applications such as air and water treatment.These surface processes were evidenced on silica-like materials by experiments by Janssen and Tuzson (2010) and by Marinov et al (2013a) and studied numerically by Lopaev et al (2011aLopaev et al ( , 2011b)), by Marinov et al (2013b) and by Guerra et al (2019).More recently, the ozone kinetics, including its surface production, was investigated by Booth et al (2023) and by Meyer et al (2023) on Pyrex (borosilicate glass).Concerning O atom recombination as a whole, several experimental and modelling works have assessed the relative importance and the characteristic parameters of E-R and L-H mechanisms in oxygen plasmas, as this evaluation can change the predictive accuracy of mesoscopic models (Macko et al 2004, Guerra 2007, Guerra and Marinov 2016, Booth et al 2019, Guerra et al 2019).Booth et al (2019) recently studied the recombination of O atoms on the surface of a Pyrex (borosilicate glass, SiO 2 + B 2 O 3 ) tube containing a pure O 2 DC glow discharge over the pressure range 0.2-10 Torr and the current range 10-40 mA, for fixed outer wall temperatures of 5 • C and 50 • C. The atom loss rate was measured through optical emission actinometry with current modulation and the atom recombination probability γ O was deduced from it.The value of γ O passes through a minimum at approximately 0.75 Torr.At pressures below this minimum, γ O decreases with pressure and increases strongly with the discharge current.This effect was attributed by Booth et al (2019) to incident ions and fast neutrals with sufficient kinetic energy to clean or chemically modify the surface, generating new adsorption sites.This idea was further developed by Afonso et al (2024).At pressures above 0.75 Torr, Booth et al (2019) found an approximate Arrhenius behaviour of γ O with the incident atom temperature.γ O falls below the Arrhenius law when the temperature and the O flux are lower.The authors explained these results by an E-R mechanism with incident O atoms recombining with both chemisorbed and more weakly bonded physisorbed atoms on the surface.Booth et al (2019) proposed a phenomenological E-R model that explains the observed behaviour while depending only on a minimum set of parameters that are adjusted to match a large set of experimental measurements.
Despite considering a low desorption frequency and highlighting the role of physisorbed atoms in O recombination, Booth et al (2019) neglect L-H mechanisms, as well as the impact of O 2 and of O 3 production on the surface.Indeed, the model by Booth et al (2019) may correspond to an effective way of capturing to some extent a more complex physical picture.However, a more accurate description of O atom surface recombination on Pyrex may be required to clarify the dominant mechanisms.Such a description should result in modelling with predictive capabilities, being able to foresee the influence of plasma-surface interactions in different systems where these play an important role.The works in the PhD thesis by Morillo-Candás (2019) and in the MSc thesis by Silveira (2023) employed a mesoscopic model to assess the influence of different surface mechanisms (those suggested by Booth et al (2019) and others) on O atom recombination in the same conditions studied by Booth et al (2019).Parametric studies were done to feed the investigation of the importance of different mechanisms.This work builds on those studies to develop a more complete description of surface kinetics of O atom recombination in Pyrex in oxygen discharges.Despite the large number of studies of O recombination on Pyrex, the results are rather inconsistent, showing a large irreproducibility attributed to different experimental conditions (e.g.different recombination probabilities are found by Booth et al (2019) and Booth et al (2020)).As such, a model validated against consistent and reproducible measurements is still lacking, that can provide a complete description and evaluation of recombination mechanisms in a widely used material in plasma reactors such as Pyrex.
In this work, mesoscopic modelling of surface processes is employed to assess the interaction between an oxygen DC glow discharge and a Pyrex surface in the conditions studied by Booth et al (2019) and in additional conditions.The experimental conditions are presented in section 2. Two experimental data sets are taken and described.One from Booth et al (2019)  The mesoscopic surface kinetics model used in this work is described in section 3, together with three different reaction sets.The simplest one is the set by Booth et al (2019), effective at describing experimental results with a minimum set of parameters but neglecting likely relevant processes.Then, the base reaction set employed in previous works (Guerra andMarinov 2016, Marinov et al 2017) is recalled and further developed into an upgraded reaction set providing a more complete description of surface kinetics while keeping the simplicity of a mesoscopic model.The upgraded model considers the kinetics of impacting O 2 and of O 3 surface production and assumes a desorption frequency dependent on wall temperature.The results of the upgraded model are validated against the two sets of experimental data in section 4.1 and the relevance of different recombination mechanisms is analysed in section 4.2.These are further discussed, together with the surface mechanisms for ozone production, via literature comparisons, in section 4.3.Finally, the conclusions are summarized in section 5.

Discharge set-up and conditions
The discharge conditions studied experimentally and numerically by Booth et al (2019Booth et al ( , 2020Booth et al ( , 2022Booth et al ( , 2023)), Western et al (2020), Viegas et al (2023)and Dias et al (2023) are assessed in this work.The experimental set-up is described in detail by Booth et al (2019Booth et al ( , 2020Booth et al ( , 2022Booth et al ( , 2023))and a simplified scheme is presented in figure 1, adapted from the work by Booth et al (2019).It consists of a DC glow discharge in O 2 , ignited in a Pyrex (borosilicate glass) tube of 2 mm thickness, 1.0 cm inner radius and 56 cm length.The discharge is ignited by electrodes located in side-arms separated by 52.5 cm.The discharge current is varied between 10 mA and 40 mA and the pressure values are changed within the interval between 0.2 Torr (0.27 mbar) and 10 Torr (13.33 mbar).In this work, only the high pressure regime identified by Booth et al (2019) for pressures of 0.75 Torr (1 mbar) and above is addressed.In these conditions, ions are assumed to have low energies and thus to not modify the surface.The surface kinetics of atomic oxygen recombination in Pyrex for pressures below 0.75 Torr, conditions in which the ion kinetic energy can have an effect on surfaces, was recently addressed by Booth et al (2019) and by Afonso et al (2024).The cylindrical tube outer surface is kept at a constant temperature variable between −20 • C and 50 • C.This is guaranteed by a water/ethanol mixture flowing through an outer envelope and connected to a thermostatic bath.The temperature drop across the Pyrex tube wall is considered to be negligible, of less than 2 K (Booth et al 2019).The gas flow rate is kept low, between 3 sccm and 10 sccm, as to ensure that the gas residence time (>1 s) is longer than the lifetime of all the active species in the discharge.The leak rate of air into the system is small, below 0.015 sccm, corresponding in the worst case to less than 0.4% N 2 in the mixture.Such a N 2 concentration in O 2 gas is not expected to have any noticeable effect neither on discharge parameters nor on the kinetics of oxygen plasma species.As a result, an axially uniform plasma column is studied with constant gas composition (almost pure O 2 ) and cylindrical symmetry, in interaction with the cylindrical Pyrex vessel.Two experimental data sets are addressed in this work, which provide the input conditions for the surface models and the experimental outputs to compare with simulations and study surface recombination.

Experimental data from Booth et al (2019): data set 1
Data set 1 was obtained by Booth et al (2019) and shared by the authors.It contains a total of 56 different conditions: 28 for wall temperature T w = 5 • C and 28 for T w = 50 • C.These include pressure values of 0.75, 1, 1.5, 2, 3, 5, 7.5 and 10 Torr and currents of 10, 20, 30, and 40 mA.For each condition, the gas temperature was measured from the rotational structure of the O 2 (b) optical emission at 762 nm.It was determined both radially averaged (T g,av ) and radially resolved, with 1 mm resolution.The gas temperature near the wall (T nw ) was determined at 1 mm from the wall from High Resolution Two-photon Absorption Laser-Induced Fluorescence (HR TALIF) measurements of the Doppler profile of O atoms.Booth et al (2019) found from measurements and 1D radially-resolved calculations that T nw is related to the other temperatures: T nw ≃ T w + 0.28 (T g,av − T w ) . (1) The atomic oxygen fraction was also measured for every condition via actinometry.This technique was also used to measure the temporal variations of O atom fraction during discharge modulation, from which is obtained the effective loss frequency of O atoms ν O loss .From that measurement, the wall recombination coefficient of O (γ O ) is also obtained, employing the expression (Booth et al 2019): where L av is the average volume loss rate of O, R tube is the tube radius of 1 cm, [O] av is the average O atom density and [O] nw is its near-wall value.v th,O is the thermal velocity of O calculated as: where k B is the Boltzmann constant and M O is the mass of the O atom.Booth et al (2019) calculated [O]av [O]nw from 1D model results of T g (r) and the assumption that the O atom density profile follows the total gas density, determined from the ideal gas law, which is reported to induce up to a 30% decrease of γ O with respect to the case considering [O] av = [O] nw .L av was computed by considering only the three body recombination of two atoms with rate coefficient k = 3.34 × 10 −30 • (1/T g,av ) • e −170/Tg,av cm 6 • s −1 , which results in L av being reported as less than 10% of ν O loss .The resulting γ O values were shared by Booth et al (2019) with the current authors and are represented in figure 2 for the two different wall temperatures.

Experimental data at different wall temperatures: data set 2
Data set 2 was obtained by Dias (2019) in a similar set-up as data set 1, with the same dimensions and characteristics, but in a different tube and measurement campaign.It contains a total of 50 conditions: 11 for T w = −20 • C, 13 for T w = 5 • C, 13 for T w = 25 • C and 13 for T w = 50 • C.These include pressure values of 0.8, 1, 1.5, 2, 3, 5 and 7.5 Torr and currents of 20 and 40 mA.The gas temperature was measured in the same way as by Booth et al (2019) for T w = 5 • C and T w = 50 • C. For the remaining wall temperatures, T g,av was calculated assuming a linear dependence of T g,av with T w .T nw was then obtained from equation (1).The atomic oxygen density was measured via actinometry and its temporal variation was obtained from the temporal evolution of line emission intensity ratios, which led to obtaining the effective loss frequency of O atoms ν O loss .2, for the different wall temperatures, together with those of data set 1. Considering the significant amount of data and of pressure and current conditions, in this work the γ O values are represented, in either logarithmic or linear scale, for different wall temperatures as function of the inverse of the gas temperature near the wall, as γ O = f(300/T nw ).We believe this way of plotting the data is the one that best represents the direct dependences of γ O .We should notice that T nw always increases with T w , pressure and current in the studied range of conditions (Booth et al 2019, 2022, Dias 2019).
Figure 2 shows not only the decrease of γ O with T w and T nw , but also the similarity between γ O from data sets 1 and 2 for T w = 5 • C and T w = 50 • C, even though the data was obtained in different tubes.Nevertheless, the data exhibit an uncertainty in reproducibility of up to 39% between the two sets.We should notice that this difference between data sets 1 and 2 is much lower than that reported by Booth et al (2020) where for similar conditions ν O loss can be more than double of the one measured by Booth et al (2019), which is attributed to different Pyrex tube histories (Cartry et al 1999).In the remaining figures in this work, an errorbar is represented to take into account the experimental uncertainty of γ O .Taking into account that the uncertainty of measurement is low and the uncertainty in the conversion of ν loss O to γ O is of the order of L av /ν O loss ≃ 5%, we consider that the experimental uncertainty is mostly due to reproducibility.To guide the eye, we add a 30% error bar to γ O measured in experiments.

Surface kinetics model
The mechanisms taking place at the surface are effectively taken into account via a mesoscopic (not atomistic) description (Kim and Boudart 1991, Guerra 2007, Guerra and Marinov 2016, Marinov et al 2017).In this model, the Pyrex surface is considered fully covered by a series of adsorption sites, which can reversibly or irreversibly hold O atoms and O 2 molecules incoming from the O 2 plasma, where dissociation takes place.Reversible adsorption is associated with physisorption, at surface sites F V of total density [F], in which the bond between the particle and the surface is due to van der Waals forces.Due to their relatively low energy bonds, physisorbed atoms O F or molecules O 2,F can diffuse along the surface or desorb back to the gas phase.Physisorption can occur at any site on the surface.Conversely, chemisorption identifies irreversible adsorption, where a true chemical bond (covalent or ionic bond) exists between the adsorbed atom and the surface.It takes place on surface defects S V of total density [S], where the bonds with O are promoted.Due to the high energy of the bond and the lack of sufficiently high temperatures to promote thermal desorption below ≃2000 K (Kim and Boudart 1991), it is assumed that a chemisorbed atom O S can only leave the surface site through recombination with another O atom.In reality, there may be a distribution of bonding energies for species adsorbed at the surface (Marinov et al 2017).As a simplification, we consider only two types of adsorbed species: physisorbed with zero energy bond and chemisorbed with relatively high energy bond.In this work, different sets of reactions and parameters are taken, which are summarized in the next subsections.The considered processes generally include physisorption, thermal desorption from physisorption sites, chemisorption, E-R recombination between a gas-phase and an adsorbed species, chemisorption via surface diffusion of physisorbed atoms, and L-H recombination between two adsorbed species.
In this work the deterministic method is numerically employed for different reaction schemes, as explained by Marinov et al (2017) and recently recalled by Afonso et al (2024).Through the deterministic method, a system of rate balance equations is solved for the surface densities n i of every species i.For reactions j with only one surface reactant k (e.g.R1-R5 in table 1), the model solves: where c R i,j and c L i,j are the stoichiometric coefficients of species i in reaction j on the right-hand side (production) and left-hand (destruction) side of the reaction equation, respectively, S j is the corresponding source or loss term and r j is the rate of reaction j in units of site −1 • s −1 .For reactions j where there is one surface reactant O F and one surface reactant k (e.g.R6-R8 in table 3), the rate balance equations are solved as: where is the fractional coverage of O F .Then: Finally, the recombination probability is calculated in this model as in the works by Guerra (2007) and Marinov et al (2017), as the quotient from dividing the surface recombination rate (written as in equations ( 4) and ( 5), including the multiplication by the number of recombining O atoms per recombination reaction) by the incoming flux of atoms to the surface, such that: where ϕ O is the thermal flux of incoming O atoms from the gas phase to the surface.
In the works by Guerra and Loureiro (2004), Guerra and Marinov (2016), Marinov et al (2017)and Afonso et al (2024), another method is used to solve the system of surface reactions, the Kinetic Monte Carlo (KMC) method.Instead of dealing with average surface coverages the KMC method follows every particle and every site on the surface.Therefore, using KMC, the state of the system can be described as a list of all active sites in the simulation domain together with their occupancies, evolving in time according to the listed reactions and their probabilities.As explained in the works by Guerra andMarinov (Guerra andMarinov 2016, Marinov et al 2017), KMC allows to obtain more accurate results than the deterministic method that takes approximations concerning the diffusion of physisorbed species.In this work, KMC is not employed, since the studied conditions are similar to those in the works by Marinov et al (2017) and by Afonso et al (2024), where it was confirmed that the KMC and deterministic results are close, within a 20% difference.
The input required for the model from the experimental data sets are pressure, T g , T w and [O] av .T nw is calculated from equation ( 1) and where N is the gas density calculated from the ideal gas law.The thermal fluxes of incoming O and O 2 particles from the gas phase to the surface are calculated as: where v th are the thermal velocities calculated as in equation ( 3).As noticed in the work by Viegas et al (2023), near the wall it may be more accurate to assume Nevertheless, it was verified that the impact of changing that assumption is very subtle in the conditions under study.The following subsections describe the different reaction sets and kinetic data employed in the surface kinetics model.

Reaction set from Booth et al (2019)
Based on the experimental results of data set 1 (see figure 2), Booth et al (2019) analysed that, for each value of T w , for pressure above 2 Torr, the logarithm of γ O falls linearly with 300/T nw , which constitutes an Arrhenius behaviour.
According to Booth et al (2019), this suggests that recombination is dominated by E-R processes that depend on the energy of incoming atoms, and thus on T nw , and not by L-H mechanisms involving diffusing physisorbed atoms, whose rates only depend on T w , and thus whose γ O contribution should be constant for a fixed T w .Booth et al (2019) further noticed that the data points for low T nw (pressures below 2 Torr) do not lie on the Arrhenius line and instead there is a 'knee' in the Arrhenius plot (in figure 2, more visible for data set 1 than for data set 2).Moreover, they assessed that the decrease of γ O with T w for the same T nw implies a varying occupation of sites with T w .Both factors suggest an involvement of physisorbed atoms.Based on this evaluation, Booth et al (2019) suggested an effective reaction set describing the atomic oxygen recombination in a simple way.Their aim was to propose a phenomenological but robust surface reaction model with a small number of adjustable parameters, validated against well-defined measurements (data set 1), thus minimizing the arbitrariness of their values.The list of proposed reactions is presented in table 1, along with the corresponding rates.The description of parameters is given in table 2. The list includes physisorption (R1), desorption (R2), chemisorption (R3) and E-R recombination with chemisorbed (R4) and physisorbed (R5) atoms.In this work, the system is solved analytically but also employing the mathematical framework previously described (equations ( 4)-( 6)).For the approximate analytical solution, equations ( 19) and ( 26) from the work by Booth et al (2019) are used: It should be noticed that the multiplication by R/2 in equation ( 11) is not presented in equation ( 26) in the work by Booth et al (2019) but was deemed necessary and added in this work.Concerning the numerical solution, the process rates are taken from Booth et al (2019) but written differently, in a compatible way with our mathematical formulation.For the choice of parameters, Booth et al (2019) considered two extreme cases regarding the activation energy for recombination: in case 1 the activation energies for E-R recombination are the same for chemisorbed and physisorbed O atoms (E S ER = E F ER ), which would mean that the activation energy depends only on the interaction between the adsorbed atom and the gas phase atom, and is independent of the nature of the bond between Table 1.List of surface reactions in the model from Booth et al (2019).V and S, respectively, are the volume and the surface of the reactor (V/S = R tube /2).[O]nw = [O]av • Tg,av/Tnw is the estimated density of O atoms near the wall.R is the ideal gas constant, Tw is the wall temperature and Tnw is the near-wall temperature.the adsorbed atom and the surface; in case 2 the E-R recombination with physisorbed atoms is considered to have zero activation energy (E F ER = 0), which corresponds to considering that the activation energy depends on the bond between the adsorbed atom and the surface, which is much lower for physisorbed than chemisorbed atoms (Kim and Boudart 1991, Guerra 2007, Marinov et al 2017).A clarification on the nature of the activation energy for recombination should resort to DFT and/or molecular dynamics calculations and falls outside the scope of this work.For each case, the obtained analytical expressions of γ O were used by Booth et al (2019) to fit the experimental results of data set 1. The recombination barriers, the desorption barrier and the recombination reaction probability were obtained from those fits.The parameters used in case 1 in the model by Booth et al (2019) are listed in table 2.
The numerical calculations employing the surface kinetics model in this work and the reactions, rates and parameters of tables 1 and 2 are compared to the approximate analytical solution and the experimental results in figure 3 Booth et al (2019), employing the analytical calculations (figures 11 and 12 in that work).As shown by Booth et al (2019), the calculated γ O obtained from 'case 2' parameters would be similar to the ones shown here ('case 1'), although with different contributions to recombination from E-R with chemisorbed and E-R with physisorbed O atoms.Booth et al (2019) found the fitting parameters of table 1 by considering γ O = ∑ j S j,O+O /ϕ O (equations ( 18) and ( 25) in that work), instead of the loss probability γ O = 2 ∑ j S j,O+O /ϕ O (equation ( 7) in this work).For that reason, there is necessarily a factor 2 difference between our model simulations employing the given parameters and the experimental γ O .In figure 3 the simulation results are divided by 2 to provide an appropriate comparison.It should be noticed that the same result is obtained by dividing the pre-exponential factor for recombination by 2, which is what might have been found from the fits by Booth et al (2019) if equation ( 7) had been used to define γ O .The accuracy of these statements is confirmed in figure 3, where the match between the analytical and numerical solutions is shown.This match means that the analytical expressions given by Booth et al (2019) are very closely retrieved by taking in our deterministic numerical model the mechanisms suggested in that paper with the rate coefficients proposed.
Figure 3 shows that the phenomenological model by Booth et al (2019) effectively correctly describes the experimental trends of γ O , with the curvature in the γ O = f(300/T nw ) dependence evincing a key role of physisorbed O atoms, which is accounted in a simplified way through O + O F recombination.However, the actual picture may be more complex and this may correspond to an effective way of describing a sequence of elementary steps (Booth et al 2019).In this work, we intend to provide a description of surface kinetics that accurately depicts the more complex physical picture of O recombination, not fully represented in the model by Booth et al (2019).Indeed, if all the expected recombination processes were included, namely L-H recombination, the Table 3. List of surface reactions in the base model (Marinov et al 2017).ϕ O is the thermal flux of gas-phase O (O(g)) to the wall (equation ( 8)), R is the ideal gas constant, Tw is the wall temperature and Tnw is the near-wall temperature.

Nbr.
Reaction simulated γ O in figure 3 would be up to two orders of magnitude higher, as shown by Morillo-Candás (2019) (figure 8.25 in that work).This is due especially to the low value of ν d (10 13 s −1 ) taken by Booth et al (2019), that would guarantee a high population of physisorbed atoms (between 10 −3 and 3 × 10 −2 in the calculations of figure 3).To provide such a more accurate description, we start by presenting our base model, employed in previous works, and then we build on it towards an upgraded model.

Base reaction set
The base model considered in this work is the one by Marinov et al (2017), based on the works by Kim and Boudart (1991), Guerra (2007)and Guerra and Marinov (2016) ) is not part of this base model.On the one hand, the bond between physisorbed atoms and the surface is relatively weak, so that recombination involving a physisorbed atom can be relatively easy.On the other hand, if the kinetic energy of the atoms arriving to the surface is relatively high, physisorbed atoms may be removed from the surface by the impingement of gas phase atoms, without recombining, or the incoming atoms may escape the potential well and proceed towards the closest O S .None of these processes are taken into account in the base model, with the atoms arriving from the gas phase to an occupied physisorption site being assumed to simply reflect (Guerra 2007).
The rates considered in table 3 are explained in the work by Guerra (2007).However, in this work, as in the works by Booth et al (2019) and Guerra et al (2019), the kinetic energy of incoming particles from the gas phase is considered, by taking T nw instead of T w in the rates involving those particles  (adsorption and E-R recombination).The parameters in table 4 were picked in the work by Guerra and Marinov (2016), based on previous studies of O recombination in SiO 2 surfaces (Kim and Boudart 1991, Seward and Jumper 1991, Gordiets and Ferreira 1998, Cartry et al 2000, Guerra 2007, Lopaev et al 2011a) and on the agreement with previous experiments (shown in figure 7 in that work).That work includes a discussion of the uncertainties of the chosen parameters, through comparisons with those considered in other works.The different ways to obtain surface parameters were also discussed in the work by Marinov et al (2017), pointing out that these can be obtained from ab-initio calculations, direct measurement, tuning to match experimental measurements and 'guesstimation'.It was concluded that even in such a relatively simple system consisting of only two chemical elements, oxygen and silicon, a significant uncertainty still remains, with the recombination mechanisms themselves not being fully understood.As such, this base model, as other mesoscopic models, remains an effective way to account for complex phenomena only through a few reactions and parameters.In figure 4, the base model simulation results of γ O are compared to the data set 1 experimental measurements, for T w = 5 • C and T w = 50 • C.
The contributions of the three different recombination mechanisms are also presented.Figure 4(a) shows that, at T w = 5 • C, the base model simulation results are close to the experiments, with the dominant recombination mechanism being L-H recombination between physisorbed atoms.This mechanism provides a visible curvature, in agreement with experiments.Nevertheless, it is visible that at T w = 5 • C and high T nw (300/T nw < 0.9) the steepness of the simulation trend and the experimental trend are different, which suggests there should be a stronger role of the O+O S recombination mechanism in those conditions.Despite the agreement at T w = 5 • C, figure 4(b) shows that the simulation results are farther from reproducing the measurements at T w = 50 • C. In this case, the simulated γ O is lower than the experimental results, up to a factor 4, and is given mostly by O+O S and O F +O S recombination mechanisms, which cannot provide the curvature observed in experiments.This analysis suggests that recombination processes providing the curvature, such as O F +O F or O+O F (see figure 3), should have a more prominent role than that attributed in the base model.Following the discussion by Guerra and Marinov (2016), we should notice that, of all the parameters of table 4, Table 5. List of surface reactions in the upgraded model.ϕ O is the thermal flux of gas-phase O (O(g)) to the wall (equation ( 8)), ϕ O2 is the thermal flux of gas-phase O 2 (O 2 (g)) to the wall (equation ( 9)), R is the ideal gas constant, Tw is the wall temperature and Tnw is the near-wall temperature.Reactions R1-R8 are the same as presented in tables 1 and 3, while reactions R9-R15 are new.

Nbr.
Reaction Rate (site the pre-exponential factor of the desorption frequency ν d is the most uncertain parameter, with values in literature between 10 13 s −1 and 10 16 s −1 , and may depend directly on T w (Ibach et al 1980).These insights are valuable to develop an upgraded model that matches the experimental results through a more complete description of surface kinetics while keeping the simplicity of a mesoscopic model.

Upgraded reaction set
The reaction set developed in this work is based on the one presented by Guerra et al (2019) and on parametric studies by Morillo-Candás (2019) and Silveira (2023) and will be further presented in this section.It is summarized by the list of reactions and rates in table 5 and the list of parameters in table 6.
In table 5, reactions R1-R8 are the same as presented earlier in tables 1 and 3, while reactions R9-R15 are new.
Starting from the base reaction set, the work by Guerra et al (2019) further considers O 2 physisorption (R9) and thermal desorption (R10) and the recombination between O 2,F and O F (R14,R15), to account for gaseous O 3 formation on the surface in the simplest possible way.The desorption energy for O 2 is taken in that work as 17.5 kJ • mol −1 , close to the value of 18.3 kJ • mol −1 suggested by Lopaev et al (2011a) by matching experimental measurements in quartz (also a silicalike surface).This value is significantly lower than that considered for O, of 30 kJ • mol −1 , also close to the value suggested by Lopaev et al (2011a) of 31.6 kJ • mol −1 .It should be noticed that processes R14 and R15 are distinguished since we consider that in each case it is the first reactant that moves towards the second reactant, having to overcome the bonding energy of the second reactant to the surface.In this work, all the recombination processes involving weakly bonded atoms and no chemisorbed atoms are considered to have no energy barrier.For that reason, R14 and R15 end up having the same rate.Moreover, the diffusion and desorption frequencies are considered to be the same for O F and O 2,F , as is the diffusion energy barrier.It is considered, as by Guerra et al (2019), that molecular chemisorption (into O 2,S ) cannot take place because the distance between two neighbouring chemisorption sites is statistically too long to allow the simultaneous chemisorption at two surface sites.The chemisorption of one atom of an incoming molecule and the ejection of the other atom to the gas phase is also neglected.
Apart from the O+O S and O F +O S recombination processes, the upgraded reaction set considers further recombination mechanisms involving physisorbed species.The possibility of L-H recombination between two physisorbed atoms (R8) was considered in the work by Rakhimova et al (2009) and was discussed and included already in the works by Guerra and Loureiro (2004), Guerra and Marinov (2016), Marinov et al (2017), Alves et al (2018), Guerra et al (2019), where it was found to be one of the dominant mechanisms.E-R recombination involving weakly bonded atoms is not included in the base reaction set of table 3, but is put in evidence in the works by Lopaev et al (2011aLopaev et al ( , 2011b))and by Booth et al (2019).This mechanism (R5) is included in the upgraded set for the sake of completion.E-R recombination mechanisms resulting in ozone formation (R11 and R13) are also considered.These mechanisms were proposed by Lopaev et al (2011aLopaev et al ( , 2011b) ) and by Marinov et al (2013b) and in those works constitute some of the main mechanisms for O surface recombination and efficient paths for ozone production in the pressure range between 10 and 50 Torr.They allowed to explain the high ozone concentrations observed in experiments that could not be explained only by gas-phase processes.The rates of these processes are written in this work following an analogy with other E-R rates, considering probability 1 of recombination interactions.Recently, the works by Meyer et al (2023) and by Booth et al (2023)  Besides the added reactions and rates, the upgraded model includes changes in parameters (table 6) with respect to the base model (table 4).Indeed, in table 6, the additional parameters and the modified parameters are signaled by different marks.Some parameters were changed to adjust the simulation results to the experimental data at T w = 5 • C (see figure 4(a)).Following the analysis of figure 4, the influence of O+O S recombination at high T nw is increased by slightly decreasing the recombination barrier of this process (E S ER ) from 17.5 kJ • mol −1 to 15.0 kJ • mol −1 .In order to counteract the increase of the simulated γ O , the barrier of O F +O S recombination (E S LH ) is slightly increased from 17.5 kJ • mol −1 to 20.0 kJ • mol −1 , with effects mostly at low T nw .Following the discussion of parameters by Guerra and Marinov (2016), these values are reasonable.They are within the intervals, or lie close to the values in literature of 15.8 kJ • mol −1 (Lopaev et al 2011a(Lopaev et al , 2011b)), 17 kJ • mol −1 (Kim and Boudart 1991), 25 kJ • mol −1 (Kovalev et al 2005) and 25.5 kJ • mol −1 (Guerra 2007).Moreover, they remain close to the diffusion barrier E D (15 kJ • mol −1 ) (Lopaev et al 2011a).While in the base model E S  LH = E S ER , this is no longer the case in the upgraded model where there is a 5 kJ • mol −1 difference between them.
It is reasonable to consider that an extra diffusion barrier of this relatively small magnitude may exist in the vicinity of a chemisorption site due to local structure modifications, as also considered by Afonso et al (2024) in the case of metastable chemisorption sites (in that case an extra barrier of 15 kJ • mol −1 was taken).Apart from these small changes, the main change in the upgraded model concerns the parameter with largest uncertainty in literature, the pre-exponential factor of the desorption frequency ν d (Guerra and Marinov 2016).By varying this parameter with wall temperature, an agreement between simulations and experiments is obtained for every case in both experimental data sets, as shown in section 4.1, without changing any of the other parameters in table 6 (most of them taken directly from table 4 of the base model).
The desorption frequency pre-exponential factor ν d was discussed by Ibach et al (1980) for a well characterized system of CO interaction on Ni(111).It was noticed that transitionstate theory predicts a linear dependence with T w for ν d for the desorption of a single atom out of a one dimensional potential: (Marschall and MacLean 2011).This justifies the usual choice of ν d = 10 13 s −1 (Lopaev et al 2011a, Booth et al 2019).However, by addressing the equilibrium between the gas phase and the adsorbed phase, Ibach et al (1980) found that ν d should have values closer to 10 15 s −1 , as used in the seminal work by Kim and Boudart (1991), rather than 10 13 s −1 , and it should be dependent on T w and adsorption site coverage.Booth et al (2023) also argued that taking into account collisions between gaseous particles and the surface can increase ν d from 10 13 s −1 to 10 15 s −1 .Ibach et al (1980) further presented experimental data demonstrating an exponential increase of ν d with the coverage, which, in its turn, decreases with T w .The same effect was reported with CO on Ru(001) surfaces (Pfnür et al 1978).According to Ibach et al (1980), such an increase of ν d with coverage is indicative of repulsive interactions between the adsorbed species, that facilitate desorption or avoid adsorption.
In this work, we estimate ν d by matching the simulation results with experimental data set 2 and then confirming the good agreement also with data set 1. We tried to correlate the estimated values of ν d with the coverage of physisorption sites (which depend not only on T w but also on discharge parameters), but no good fit or clear dependence of ν d on coverage was found.Nevertheless, we find an exponential growth of the estimated ν d with the median coverage for each T w , in agreement with Ibach et al (1980), which can also be described as an exponential decay of ν d with T w .Indeed, the approximate values of ν d that match experimental data set 2 are 2.0 × 10 15 s −1 for T w = −20 • C (253.15 K), 9.0 × 10 14 s −1 for T w = 5 • C (278.15 K), 4.5 × 10 14 s −1 for T w = 25 • C (298.15 K) and 3.0 × 10 14 s −1 for T w = 50 • C (323.15 K).These values are within the limits discussed by Guerra and Marinov (2016) (10 13 s −1 and 10 16 s −1 ) and are perfectly fitted by the expression represented in figure 5: The fit is given by equation ( 13).
In this way, by varying the most uncertain and temperature dependent parameter, ν d , we find an effective way to account for more complex, multi-step kinetics phenomena absent from mesoscopic models.It should be noticed that it is experimentally difficult to distinguish differences in ν d from differences in the desorption barrier energies E O d and E O2 d and that the same effect could be obtained by varying E O d and E O2 d with the wall temperature (or surface coverage), instead of changing ν d .Indeed, by fitting experimental results, Tait et al (2006) proposed decreasing desorption barriers of alkanes/hydrocarbons on different surfaces with surface coverage, for small coverages, indicating the presence of adsorption sites with different adsorption energies.This would mean increasing desorption barriers with T w , which is coherent with the results in this work.The variation of ν d is preferred in this work with respect to a possible variation of E O d and E O2 d because ν d is more uncertain (Guerra and Marinov 2016).Further research is required to clarify this point.
To highlight the role of the different changes between the base and upgraded models, figure 6 shows the comparison with experimental data set 1 of the simulation results of the base model, the base model with ν d (T w ) and the upgraded model.The comparison demonstrates that the consideration of the wall temperature dependent pre-exponential factor of the desorption frequency is essential to bring the γ O results to the appropriate values at each T w , but the remaining changes and additions are still necessary to obtain an adequate evolution of γ O with T nw .The relevance of each recombination mechanism of the upgraded model and its influence on the shape of the γ O (300/T nw ) curve is further shown and explained in section 4.2.

Validation of the upgraded model against experimental results
The results of all the models presented in the previous section were already compared to experimental data set 1 in figures 3, 4 and 6. Figure 7 shows the recombination probability γ O as function of the reverse of T nw for the conditions of data set 2. The experimental results are compared to the calculations of the analytical model by Booth et al (2019), already shown to provide the same results as the corresponding numerical model, to the simulations of the base model and to those of the upgraded model, all described in section 3.
Figure 7 shows that the simulation results employing the upgraded model match very well, mostly within the ±30% error bars, the measurements of data set 2. This may not look too surprising since the pre-exponential factor of the desorption frequency ν d (T w ) is chosen for each T w to match these measurements.Nevertheless, the results show a remarkable consistency of the model for a significant range of T nw (corresponding to different conditions of current and pressure) and of T w , just by assuming a smooth T w dependency of a single parameter, ν d (T w ), with all the other parameters remaining based on literature and close to the base reaction set.The consistency and agreement with experiments is also maintained when changing the experimental conditions, in this case from those of data set 1 (figure 6) to those of data set 2. For both data sets, it is shown that the base model presents a worse agreement with experiments than the upgraded model, calculating too high loss probability γ O at low T w and too low γ O at high T w .Also the dependency of γ O from T nw is more accurately described by the upgraded model than the base model, as it is visible that the steepness of the decrease of γ O with 300/T nw is generally higher in experiments than in the base model simulations.
Finally, figure 7 shows that the model by Booth et al (2019) presents a good agreement with the measurements of data set 2 for T w = 5 • C and for T w = 50 • C, as in the case of data set 1.
Moreover, the agreement within the errorbars is confirmed for T w = 25 • C, the case in between the previous two wall temperatures.Nevertheless, the model by Booth et al (2019) overestimates γ O for T w = −20 • C, which points to a limitation of the model, which has only two recombination reactions and no T w dependence in any parameter, to being extended to varying conditions.

Recombination mechanisms
The validated upgraded model is then used to analyse through which mechanisms is atomic oxygen recombined back into molecular form in the Pyrex wall for different conditions.This is done by assessing the contributions of the different recombination reactions in table 5 to the calculation of the recombination probability γ O .The comparison of the different contributions is shown firstly for the conditions of data set 1 in figure 8.Besides the total γ O from experiments and simulations, the main contributions to γ O for each wall temperature case are displayed.These are the 4 most relevant recombination mechanisms in each case among the 5 following processes: From the O recombination processes of table 5 (R4, R5, R7, R8, R11, R13, R14 and R15), we should notice that R13-15 are negligible due to the O 2,F fraction of physisorption sites (between 2 × 10 −7 and 3 × 10 −6 ) being one to two orders of magnitude lower than the O F fraction (between 8 × 10 −6 and 4 × 10 −5 ), according to the upgraded model in the conditions of data set 1.This takes place in spite of O 2 having between 4 and 27 times higher number density in the plasma than O, mostly due to the higher surface desorption barrier considered for O (30 kJ • mol −1 ) than for O 2 (17.5 kJ • mol −1 ), which is determinant for the loss of physisorbed species.As a result, R11 is the dominant reaction for O 3 wall production, as pointed out by Booth et al (2023).
Figure 8(a) shows that for low wall temperature, T w = 5 • C, the main recombination mechanisms according to the model  are L-H with O F + O F (R8) and E-R with O 2 + O F (R11), which are also the main processes inducing the curvature in the results, represented as log(γ O ) vs300/T nw .The curvature is the result of the O F fraction tendentiously increasing with T nw while the particle fluxes from the gas phase have a complex dependency of T nw (whose rise increases the thermal velocity), pressure and current, since the dissociation degree tends to decrease with pressure and increase with current, within It should be noticed that the contribution of E-R with O + O S (R4) depends only on T nw , which increases with T w , with that process becoming relevant for high values of T nw , in agreement with literature (Macko et al 2004, Guerra 2007).The remaining contributions do not present such a sharp increase with T nw due to not having a T nw dependent rate or not depending on an increasing flux of atomic oxygen with T nw .As a result, they become less relevant with increasing T nw , when compared to E-R with O + O S (R4).Moreover, the contribution of L-H recombination with O F + O S (R7) increases with T w , following the dependence of the recombination rate.In the opposite direction, the relevance of E-R with O 2 + O F (R11), of E-R with O + O F (R5) and of L-H with O F + O F (R8) decreases with T w due to the higher desorption rate of physisorbed atoms.
The analysis of the contributions of the different recombination mechanisms is continued in figure 9, where the same comparisons are presented as in figure 8  The dependence of the different contributions to recombination on T w is further shown in figure 10, where γ O simulated with the upgraded model in some conditions of data set 2 is presented as function of T w , along with the contributions of the five major recombination mechanisms.Two cases of discharge current (20 mA and 40 mA) and pressure (0.8 Torr and 3 Torr) are considered, in a total of 4 subfigures.To assist the interpretation of the results, figure 11 presents the surface site coverages for the same conditions.The surface fraction of O S is not shown because it always occupies close to 100% of the chemisorption sites, i.e. a 1.5 × 10 −2 fraction of the total number of sites.It is clearly visible from figure 10 that the contributions of E-R with O + O S (R4) and L-H with O F + O S (R7) increase with T w , in the first case due to the increasing flux and in the second case due to the increasing rate.They are thus expected to be the dominant recombination mechanisms for T w higher than 50 • C and, in the case of O + O S , to keep increasing with T w (Guerra and Marinov 2016) 10 that increasing pressure generally induces higher recombination probabilities, in agreement with the higher coverage of O F presented in figure 11.The flux of O 2 responsible for the O 2 + O F recombination producing ozone also increases with pressure.The dependence on current is less clear in figure 10, although there is also a general increase of γ O with current.One aspect that can be noticed when increasing current from 20 mA to 40 mA is that the contribution from O + O F approaches that of O 2 + O F , due to the increasing dissociation degree in the plasma.

Discussion on recombination mechanisms and ozone production
The measured and simulated γ O are within the expected range for Pyrex for the studied conditions, between 10 −4 and 10 −2 (Pagnon et al 1995, Gordiets and Ferreira 1998, Macko et al 2004).In the work by Macko et al (2004), where γ O was measured in Pyrex in the afterglow, E-R mechanisms appear to be dominant for wall temperatures above 250 K (≃ − 20 • C), due to a linear evolution of log(γ O ) vs1000/T w , without distinction between T w and T nw .These results were described by calculations by Guerra (2007).In this work, we find almost linear dependences of log(γ O ) with 300/T nw for wall temperatures between 253 K and 323 K, while having both E-R and L-H mechanisms being significant for recombination, which may not be contradictory with the results by Macko et al (2004).In fact, L-H was suggested to be dominant or as important as E-R for O recombination on silica-like materials in the range of T w between 250 K and 500 K in the works by Kim and Boudart (1991), Gordiets et al (1996), Gordiets and Ferreira (1998), Cacciatore et al (1999), Guerra and Marinov (2016), Marinov et al (2017), Guerra et al (2019).Those studies assume E-R recombination to take place between O and O S .In this work the picture is more complex, with E-R mechanisms O 2 +O F and O+O F being as significant or even more than O+O S for T w up to 298.15 K (25 • C).Within this description, the linear dependence of log(γ O ) with 300/T nw can be obtained from both E-R and L-H mechanisms, while L-H with O F +O F is shown to be fundamental to obtain a curvature at low T nw .We should add that O+O F is not a dominant recombination mechanism in the upgraded model, being consistently less significant than O 2 +O F and O F +O F recombination mechanisms.This is not the case employing the reaction set by Booth et al (2019), where the desorption frequency is low and O+O F may be an effective way to consider other recombination mechanisms involving O F .Indeed, employing the upgraded model, the estimated average number of hops of a physisorbed atom during its lifetime on the surface (calculated as in equation ( 10) in the paper by Marinov et al (2017)) lies between 6 and 10 in the conditions under study.This means that long distance diffusion of O F may not have an important role and thus the differences between considering an effective O+O F recombination mechanism and a more detailed scheme may be more conceptual than practical.
The surface kinetics simulations results should also be discussed in terms of coherence with the knowledge of the plasma under study, namely in what concerns ozone wall production (via R11 and R13-15).Janssen and Tuzson (2010) studied the isotopic compositions of ozone in an oxygen discharge in Pyrex at pressures between 1.75 and 12.5 Torr and concluded that the atom recombination coefficient into ozone γ O→O3 should be of the order of 4 × 10 −3 at room temperature and slightly higher at lower temperatures.This value was considered in the recent work by Meyer et al (2023).The model results by Lopaev et al (2011aLopaev et al ( , 2011b) also reported γ O→O3 in an oxygen glow discharge in quartz increasing almost linearly with pressure between 5 and 50 Torr up to 6 × 10 −3 (approximately 5 × 10 −4 at 5 Torr and 10 −3 at 10 Torr), via the O+O 2,F (R13) and O 2 +O F (R11) processes.The analysis of measurements by Marinov et al (2013a) in oxygen postdischarge after short pulses found that at room temperature and 5 Torr, ozone production should account for approximately 25% of the oxygen atom loss rate on the surface of non-treated silica fibres.Nevertheless, the combined experimental and modelling study by Marinov et al (2013b) on the same pulse post-discharge in the pressure range of 1-5 Torr for discharge currents between 40 and 120 mA concluded that ozone formation at the wall does not contribute significantly to the total ozone production.In that work, the total γ O was assumed to be 2 × 10 −4 , lower in post-discharge than under direct plasma exposure, as is typically the case in literature.Upper limits were suggested by Marinov et al (2013b) for γ O→O3 as 10 −5 at 1 Torr and 10 −4 (half of γ O ) at 5 Torr.Simulations of similar post-discharge conditions by Guerra et al (2019) found slightly lower γ O→O3 , within the upper limits presented by Marinov et al (2013b).Under conditions similar to those of this work, for 40 mA, 20 • C wall temperature and pressures between 0.5 and 4 Torr, Booth et al (2023) reported simulation results in the early afterglow of O 3 surface production fluxes of 2 × 10 18 − 1.5 × 10 19 m −2 • s −1 and O fluxes to the wall of the order of 1 − 5 × 10 22 m −2 • s −1 , which results in γ O→O3 of the order of 10 −4 .The simulation results of γ O→O3 presented in this work clearly increase with pressure and O 2 density, as in the literature, and they decrease with current and T w , due to dissociation and desorption, respectively.They have values between 10 −5 and 10 −3 , within the ranges found in the discharge studies by Janssen and Tuzson (2010) and by Lopaev et al (2011aLopaev et al ( , 2011b) ) and in the post-discharge study by Booth et al (2023), but above the upper limit of 10 −4 indicated for post-discharge by Marinov et al (2013b).
Another way to verify O 3 wall production is to assess its self-consistency with simulation results of O 3 volume losses.This will be done more consistently in the near future by coupling the surface model to the global model presented by Dias et al (2023).However, a first verification can be done here that the ozone surface production rate, obtained mostly via reaction R11 but also through R13-15, is consistent with the expected ozone loss rate in the plasma.The ozone production rate can be calculated as 2R tube , where γ O→O3 = γ 11 + γ 13 + γ 14 + γ 15 , v O,th is the atomic oxygen thermal velocity at temperature T nw and R tube is the reactor tube inner radius of 1 cm.For data set 1, for T w = 50 • C, ν O→O3 lies between 1 s −1 and 13 s −1 , within the range of O loss frequency measured by Booth et al (2019).Therefore, S O3 is limited between 3 × 10 21 m −3 • s −1 (low pressure 0.75 Torr) and 9 × 10 22 m −3 • s −1 (high pressure 10 Torr).In the simulation results of the zero-dimensional chemical kinetics model by Viegas et al (2023) (without O 3 surface production or loss) in the conditions of that work, addressing the same discharge as in this work for 30 mA discharge current and T w = 50 • C, the calculated total loss rate of O 3 is 1.2 × 10 21 m −3 • s −1 at 0.75 Torr and 8.5 × 10 22 m −3 • s −1 at 10 Torr, slightly lower than the production rate in this work.However, in the work by Viegas et al (2023) the O 3 density from zero-dimensional results is shown to be underestimated with respect to the O 3 density from 1D results by up to a factor 3. A higher O 3 density would of course result in a higher loss rate, coherent with the production rate in this work.As such, the O 3 production rate found in this work is reasonable and the associated mechanism (R11) may be a valuable complement to the reaction set by Dias et al (2023).

Conclusions
This work addresses surface kinetics in an oxygen DC glow discharge in a Pyrex (borosilicate glass) tube with 1 cm inner radius at intermediate pressure.Oxygen surface recombination is studied via mesoscopic modelling and comparison with experiments.Two experimental sets of conditions and of measured data are assessed: one previously published and a new one.They compose a total of 106 experimental conditions with discharge current varying between 10 and 40 mA, pressure values ranging between 0.75 and 10 Torr, and fixed outer wall temperatures (T w ) of −20 • C, 5 • C, 25 • C and 50 • C. The lower pressure conditions addressed recently by Afonso et al (2024), in which work ion-induced reversible surface modification was put forward, are left out of this study.
Three reaction sets for mesoscopic surface modelling are presented and compared: one based on the work by Booth et al (2019), one by Marinov et al (2017) and an upgraded reaction set developed in this work.It is shown that the analytical model by Booth et al (2019) can be reproduced numerically and its results of recombination probability (γ O ) agree well with measurements at T w = 5, 25 and 50 • C but overestimate measurements at T w = −20 • C. The model developed in this work is based on the mesoscopic model by Marinov et al (2017) but adds further reactions to it and is inspired by the parametric studies by Morillo-Candás (2019) and Silveira (2023).The upgraded model considers the pre-exponential factor of the O and O 2 desorption frequency as being exponentially dependent on the inverse of the wall temperature.This is based on insights by Ibach et al (1980) and attributed to repulsive interactions between adsorbed species, whose population decreases with wall temperature.After these upgrades, the model is validated by a remarkable agreement of γ O with experiments for the 106 studied conditions of wall temperature, pressure, current and near-wall temperature.This validation effort is a great display of versatility and of the expected predictability when changing conditions.It provides confidence in results and allows to use the model to deepen the knowledge in such complex reaction kinetics.Despite the good results and the indications that the pre-exponential factor of the desorption frequency may depend on the wall temperature, this consideration may as well be an effective way of capturing more complex kinetics, such as multi-step processes, varying adsorption and desorption energy barriers and varying distributions of reactivity of surface sites.The chosen formulation of the upgraded model allows to provide a more complete description of the plasmasurface interactions, while keeping the simplicity of a mesoscopic model.
The As such, the more complex picture introduced in this work allows to describe experimental results of O recombination probability through a combination of both E-R and L-H mechanisms with different reactants.Finally, the ozone wall production rate predicted in this work for Pyrex is considered reasonable and is shown to be relevant not only for O recombination but also for the global ozone production rate in plasmas in the studied conditions.The ozone wall production rate clearly increases with pressure, with relevance also for higher pressure discharge regimes.In the future, we intend to couple the surface kinetics model to a plasma global model and extend it to other pressure regimes and gas mixtures.
The same formalism developed in this work can be used to study recombination for other surfaces and atoms, both using a deterministic method as presented here or a Kinetic Monte Carlo model.However, each case has its one specificities, that needs to be taken into account.First, the appropriate surface parameters (active surface site density, adsorption energies, diffusion barriers and frequencies, etc) have to be considered, depending as well on the metallic or dielectric nature of the surface.Obtaining these parameters cannot be achieved solely within a mesoscopic model and requires fundamental experimental studies, such as the ones reported in this work, and/or the support from ab initio calculations (Marinov et al 2017).An assessment of the need to consider different elementary processes is also required.For instance, physisorption on top of chemisorption sites, dissociative chemisorption (Nave et al 2014, Kroes 2021), abstraction (Sholl 1997, Khanom et al 2003) or a distribution of reactivity among the adsorption sites (Donnelly et al 2011, Marinov 2019) may have to be included.Nevertheless, once the relevant mechanisms have been identified, the extension of the current model to describe any of them is straightforward.
at wall temperatures of 5 • C and 50 • C. Another one from the MSc thesis by Dias (2019) in similar conditions at wall temperatures of −20 • C, 5 • C, 25 • C and 50 • C.

Figure 1 .
Figure 1.Scheme of the experimental set-up.Adapted from Booth et al (2019).© IOP Publishing Ltd.All rights reserved.
, for T w = 5 • C (figure 3(a)) and for T w = 50 • C (figure 3(b)).The contributions of the two E-R recombination mechanisms are also represented.The same type of figure was shown by

Figure 3 .
Figure 3. O atom wall loss probability in Pyrex as function of the reverse of the near-wall temperature, for (a) Tw = 5 • C and (b) Tw = 50 • C. Results from experimental measurements from data set 1 (Booth et al 2019) and from analytical and numerical calculations employing the reaction set from Booth et al (2019) (case 1 in that work).The numerical simulation results are divided by 2 for reasons explained in the text.The calculations presented include the total γ O (× and +) and the contributions of the different recombination mechanisms.

Figure 4 .
Figure 4. O atom wall loss probability in Pyrex as function of the reverse of the near-wall temperature, for (a) Tw = 5 • C and (b) Tw = 50 • C. Results from experimental measurements from data set 1 (Booth et al 2019) and from simulations employing the base model.The simulation results presented include the total γ O (•) and the contributions of the different recombination mechanisms.

Figure 5 .
Figure 5. Desorption frequency as function of the wall temperature.ν d estimated by matching simulation results of γ O with the measurements of data set 2 for Tw = −20 • C, 5 • C, 25 • C, 50 • C.The fit is given by equation (13).

Figure 6 .
Figure 6.O atom wall loss probability in Pyrex as function of the reverse of the near-wall temperature, for (a) Tw = 5 • C and (b) Tw = 50 • C. Results from experimental measurements from data set 1 (Booth et al 2019) and from simulations employing the base model, the base model with modified ν d (Tw) and the upgraded model.

Figure 7 .
Figure 7. O atom wall loss probability in Pyrex as function of the reverse of the near-wall temperature, for (a) Tw = −20 • C, (b) Tw = 5 • C, (c) Tw = 25 • C and (d) Tw = 50 • C. Results from experimental measurements from data set 2 and from the three different models presented in section 3.

Figure 8 .
Figure 8. O atom wall loss probability in Pyrex as function of the reverse of the near-wall temperature, for (a) Tw = 5 • C and (b) Tw = 50 • C. Results from experimental measurements from data set 1 (Booth et al 2019) and from simulations employing the upgraded model.The simulation results presented include the total γ O (•) and the contributions of the 4 most relevant recombination mechanisms.

Figure 9 .
Figure 9. O atom wall loss probability in Pyrex as function of the reverse of the near-wall temperature, for (a) Tw = −20 • C, (b) Tw = 5 • C, (c) Tw = 25 • C and (d) Tw = 50 • C. Results from experimental measurements from data set 2 and from simulations employing the upgraded model.The simulation results presented include the total γ O (•) and the contributions of the 4 most relevant recombination mechanisms.
but for data set 2, thus extending the wall temperature range of the analysis.The previous assessment is coherent for the different conditions of data set 2. We can conclude from the figures that, according to the model presented in this work and in the studied conditions, L-H with O F + O F (R8) and E-R with O 2 + O F (R11) are the dominant oxygen recombination mechanisms for T w = −20 • C and T w = 5 • C, while for T w = 25 • C and T w = 50 • C the contributions of E-R with O + O S (R4) and L-H with O F + O S (R7) also play a relevant role.

Figure 10 .
Figure 10.O atom wall loss probability in Pyrex as function of the wall temperature for the conditions of data set 2, for (a) I = 20 mA and p = 0.8 Torr, (b) I = 20 mA and p = 3 Torr, (c) I = 40 mA and p = 0.8 Torr and (d) I = 40 mA and p = 3 Torr.Results from simulations employing the upgraded model.The simulation results presented include the total γ O (•) and the contributions of the 5 most relevant recombination mechanisms.

Figure 11 .
Figure 11.Fractions of physisorption site occupation as function of the wall temperature for the conditions of data set 2 and figure 10.Occupation by (a) O F and (b) O 2,F .Results from simulations employing the upgraded model.
analysis of the upgraded model results highlights that L-H recombination with O F + O F and E-R recombination with O 2 + O F are the dominant oxygen recombination mechanisms in Pyrex for T w = −20 • C and T w = 5 • C. For T w = 25 • C and T w = 50 • C the contributions of E-R with O + O S and L-H with O F + O S also play a relevant role for O surface recombination.
In this work, γ O for data set 2 is calculated from equation (2), by taking the approximation [O] av /[O] nw = 1, and calculating L av as was done by Booth et al (2019) for data set 1. It should be noticed that in the work by Viegas

Table 2 .
Booth et al (2019) in the model fromBooth et al (2019)(case 1 in that work).It is found in this work that P 0

Table 6 .
List of parameters in the upgraded model.To distinguish new parameters from those already present in the base model, the mark ( * ) signals parameters modified in the upgraded model and the mark (#) signals parameters that are introduced as new in the upgraded model.
. The contributions of the mechanisms involving O F and not O S (O + O F , O 2 + O F and O F + O F ) decrease with T w due to the decreasing fraction of O F , shown in figure 11.These mechanisms are expected to keep being dominant for T w lower than −20 • C. It is also noticeable in figure