Chemical kinetics in an atmospheric pressure helium plasma containing humidity

Atmospheric pressure plasmas are sources of biologically active oxygen and nitrogen species, which makes them potentially suitable for the use as biomedical devices. Here, experiments and simulations are combined to investigate the formation of the key reactive oxygen species, atomic oxygen (O) and hydroxyl radicals (OH), in a radio-frequency driven atmospheric pressure plasma jet operated in humidified helium. Vacuum ultra-violet high-resolution Fourier-transform absorption spectroscopy and ultra-violet broad-band absorption spectroscopy are used to measure absolute densities of O and OH. These densities increase with increasing H2O content in the feed gas, and approach saturation values at higher admixtures on the order of 3 10 cm 3 for OH and 3 10 cm 3 for O. Experimental results are used to benchmark densities obtained from zero-dimensional plasma chemical kinetics simulations, which reveal the dominant formation pathways. At low humidity content, O is formed from OH by proton transfer to H2O, which also initiates the formation of large cluster ions. At higher humidity content, O is created by reactions between OH radicals, and lost by recombination with OH. OH is produced mainly from H2O + by proton transfer to H2O and by electron impact dissociation of H2O. It is lost by reactions with other OH molecules to form either H2O + O or H2O2. Formation pathways change as a function of humidity content and position in the plasma channel. The understanding of the chemical kinetics of O and OH gained in this work will help in the development of plasma tailoring strategies to optimise their densities in applications.


Introduction
The interaction of non-thermal atmospheric pressure plasmas (APPs) with biological matter and their potential applications as biomedical devices [1][2][3][4] are currently a topic of significant interest. APPs have been shown to be effective in many different areas of biomedicine, such as sterilization, 5-7 cancer treatment, [8][9][10][11][12] and wound healing, [13][14][15] and have recently been identified as potential triggers of beneficial immune responses. 16 First trials on patients confirm the effectiveness of APPs. 13,17,18 APPs may offer advantages compared to conventional therapeutics due to their typically small dimensions, offering the possibility of locally confined treatment, low production cost, and the potential to tailor sources for specific applications.
A key question in plasma interactions with biological matter is the role played by plasma produced reactive species (RS). RS are known to interact with cells and their membranes, and often serve as signaling agents in cell metabolism. 19,20 They can also cause severe damage to cells at high concentrations. 19,21 For APPs to fulfil their potential in any biomedical application, a full characterization of the sources used to produce them is necessary, including the quantification of RS produced. Reactive oxygen and nitrogen species (RONS), such as atomic oxygen and nitrogen (O and N), ozone (O 3 ), excited states of molecular oxygen (e.g. O 2 (a 1 D)), or nitric oxides, have previously been quantified both experimentally and numerically in O 2 and N 2 containing plasmas. [22][23][24][25][26][27][28] Here, the production of RS in an enclosed APP operating in helium with small contents of humidity is investigated. Water is typically present in the direct vicinity of biological material, and can easily enter the gas phase via evaporation. Therefore, RS produced from water vapor can be created during the treatment of the material when plasmas are applied. Water is also usually present as a feed gas impurity. 29 Therefore, the investigation of RS directly produced from water vapor, such as O and the hydroxyl radical OH, is of interest for biomedical applications. These species can act as precursors for longer-lived species such as hydrogen peroxide (H 2 O 2 ), an important signaling agent in cells, 19,30 and O 3 . In high concentrations both of these species can have toxic effects on biological material.
The quantification of RS in APPs represents a challenge for diagnostics based on optical emission from excited states, since the plasma emission is strongly quenched by the ambient gas due to the high pressure. Laser Induced Fluorescence (LIF) and Two-photon Absorption LIF (TALIF) have been previously used to detect species such as O and OH produced from water vapor. [31][32][33][34] However, in order to accurately predict the effect of quenching using these techniques at atmospheric pressure, the densities of all potential quenching particles are needed. This is increasingly challenging in complex gas mixtures and in regions with gradual gas mixing, like the plasma effluent. These techniques also rely on quenching rate coefficients for investigated species with all possible quenchers, which for some cases, particularly quenching involving water molecules, are only poorly known. The implementation of faster laser systems such as picosecond or femtosecond lasers [35][36][37] can help to quantify the effect of these quenching processes. In addition to accounting for the effects of quenching to obtain absolute density measurements using TALIF, an additional calibration measurement involving a gas with a known quantity is typically needed. An alternative diagnostic technique, which is independent of collisional quenching, is mass spectrometry. This technique has recently been used to detect RS such as OH and H 2 O 2 produced from water vapor in the plasma effluent. 38 Similar to LIF and TALIF, this technique requires a calibration measurement to obtain absolute species densities. Mass spectrometry has also been used to detect high order protonated water clusters 39,40 produced in APPs. Species deposited in a liquid by plasma treatment are sometimes investigated by means of absorption spectroscopy in the liquid phase, and electron paramagnetic resonance spectroscopy. 41,42 However, to calculate gas-phase densities from liquid-phase densities, usually a calibration is required.
An established optical diagnostic technique for the quantification of OH in the gas-phase is ultra-violet (UV) Absorption Spectroscopy (AS), 34,[43][44][45] which is independent of collisional quenching and does not require an additional calibration measurement. However, measuring ground state densities of atomic species produced in water-containing plasmas, such as O, is challenging using AS since the energy gaps between the ground and excited states of the atoms are large. Therefore, the required excitation wavelengths typically lie in the vacuum ultra-violet (VUV) spectral range, which is strongly absorbed by air. However, atomic species, in particular O and N, have previously been quantified in an APP using synchrotron radiation and a spectrometer with an ultra-high spectral resolution, so-called VUV high-resolution Fourier-Transform Absorption Spectroscopy (VUV-FTAS). 22,23 In this work, we combine VUV-FTAS and UV broad-band AS (UV-BBAS) to determine absolute densities of gas-phase O and OH in a radio-frequency APP jet operated in helium (He) for different values of humidity up to 1.3%. We combine the experimental investigations with zero-dimensional, plug-flow plasma simulations to model the chemical kinetics in the source. These models are commonly used to study properties of atmospheric pressure plasmas. [46][47][48][49][50][51] The role of humidity on the plasma chemistry in APPs has been subject of numerical investigations in the past, 48,[52][53][54] that have established the baseline understanding of these systems. We build upon these prior works by comparison of modeling results to experiments performed for the same conditions. Species densities are measured and simulated mainly in the plasma bulk, to validate the reaction mechanism. The resulting reaction mechanism can then be used to investigate important formation pathways for different RS, and to predict additional species densities that are difficult to measure. In many applications, reactive species exit the plasma source into ambient air where the chemical kinetics will differ from the active plasma region. This transition is not investigated here, but the validated reaction mechanism constructed in this work will act as a base to be built upon for future studies in this area.

Atmospheric pressure plasma jet
The production of atomic oxygen (O), and hydroxyl radicals (OH), in an atmospheric pressure plasma jet (APPJ) operating in 5 slm helium (either grade 4.6 with 25 ppm N 2 and 7 ppm O 2 impurities, or grade 5 with 3 ppm H 2 O and 2 ppm O 2 impurities) with water vapor (H 2 O) admixtures is investigated. The plasma jet used in this work is shown in Fig. 1 and is the same as described by Dedrick et al. 23 The jet has a plane-parallel electrode configuration. One electrode with an area of 2.4 Â 0.86 cm 2 is powered by a sinusoidal voltage at a frequency of 13.56 MHz while the other electrode, which is the housing of the source, is grounded. Powered and grounded electrodes are separated by 0.1 cm, keeping the critical dimensions and operating parameters close to the 'COST Reference Microplasma Jet', 55 but with a smaller surfaceto-volume ratio (22 cm À1 here instead of 40 cm À1 for the 'COST Reference Microplasma Jet'). An impedance matching network unit (L-configuration) is used to optimise the power coupled into the plasma. The applied voltage across the gap is monitored using a high-voltage probe. A list with the equipment for the power coupling and voltage monitoring can be found in Appendix B.
We generally conduct experiments at a fixed generator power. For measurements of OH, we set the generator power close to the arcing point of the plasma in pure helium, which is independent of the generator power when using different generators and which typically occurs around 520 AE 10 V pp . Starting at this point also maximises the measurement range with respect to water content. The water content is varied while keeping the settings on the generator constant. The implications of this on the coupled power to the plasma will be further discussed in Section 2.4. At higher voltages, the plasma tends to extend around the powered electrode when operated in pure He, and transitions from a homogeneous glow-like discharge into a constricted ''arc'' mode at the electrode edges, which can damage the source.
For the O measurements, the source is operated in a vacuum vessel, with limited options for visually identifying the plasma mode through a transparent vacuum flange. In this case measurements are carried out at a lower voltage of 470 V pp to avoid the 'constricted mode' and related damage of the electrodes.
For the spatially resolved measurement of OH (presented in Section 4.1), a different plasma source, as described elsewhere, 22 has been used. It utilises the same design concept, i.e. the same gap size of 1 mm, but a slightly larger electrode area of (1.1 Â 3) cm 2 , compared to the source described earlier.
Since the surface-to-volume ratios for both designs are very similar, we assume the RS production to be comparable under similar operating conditions (gas flow, power density).
Water vapor is admixed into the He flow using two mass flow controllers and a homemade bubbler, which consists of a domed glass adapter (Biallec GmbH) clamped to a KF40 flange. Two stainless steel pipes welded to the flange provide gas in-and out-lets. Both mass flow controllers are fed with dry He, while the outlet flow of one controller passes through the bubbler before being mixed with the other. By changing the ratio of the humidified to the dry He flow, the water vapor content of the total gas flow can be regulated. 34,43,45 With the humidity level in the He flow leaving the bubbler being saturated (see below), the vapor pressure p vap H 2 O can be calculated using the semi-empirical expression given below. 56 The total amount of water in the vapor phase can then be calculated using the vapor pressure of H 2 O (in bar) and the flow rate of the He through the bubbler F bubbler He , as described in ref. 34 where T w is the water temperature in 1C. To check that the He flow exiting the bubbler is saturated with H 2 O, the weight loss of the bubbler due to evaporation of the water is measured for different He flow rates as a function of time. The absolute water concentration in the He flow is given by the ratio of the two former quantities. The results are shown in Fig. 2. A small systematic drop of the measured water concentration with increasing He flow is observed. This may reflect a temperature drop of the water inside the bubbler, which is not temperature controlled, either because of an increased evaporation rate, or fluctuations in the laboratory room temperature, since measurements were taken over several days. The averaged result for the water concentration of (16.9 AE 2.0) g m À3 corresponds to a water temperature of (20 AE 2.0) 1C assuming full saturation. Since this temperature represents the typical 'room temperature' at which the measurements were taken, we can confirm that the He flow out of the bubbler is saturated with water vapor. The uncertainty of 2 1C would lead to an uncertainty of approximately 14% in the calculation of the water vapor content in the gas.

VUV high-resolution Fourier-transform absorption spectroscopy
Absolute line-averaged O atom ground state densities are measured at the DESIRS beamline at the synchrotron SOLEIL, 57 with its unique ultra-high resolution VUV Fourier-Transform spectrometer 58 able to cover the complete VUV spectral range down to 40 nm with a resolving power (l/Dl) of up to 10 6 . The atomic oxygen transition O(2p 4 3 P J=2 -3s 3 S 1 ) is investigated in this work.
The measurement and analysis procedure is described in detail elsewhere. 22 The spectrometer yields a transmission spectrum S T , which includes the convolution of the plasma transmission T (accounting for Doppler and pressure broadening of the corresponding spectral line profile) with the sinc-shaped instrumental function Fðs 0 Þ ¼ sinðpðs 0 À sÞÞ pðs 0 À sÞ of the FT spectrometer where s is the wavenumber of the transition and S 0 the reference spectrum without absorber. Absolute densities are obtained from the transmission spectra using Beer-Lambert's law where A(s) is the absorbance, l the length of the absorbing medium (here defined by the width of the electrode as 8.6 mm), and k the absorption coefficient, which includes the ground state density, the spectral line ''Voigt'' profile, which takes into account the Doppler and pressure broadening of the spectral line, the statistical weights g J for the different states and the transition probabilities. For the evaluation of these transmission spectra, the different broadening mechanisms are taken into account as fixed values during the fitting process. The instrumental broadening is set to Ds I = 0.87 cm À1 as described elsewhere. 22 The Doppler width Ds D = 0.24 cm À1 is calculated for a gas temperature T g = 304 K. This value is determined as the OH rotational temperature from the absorption measurements by fitting the OH(X 2 P i , u 0 = 0) -OH(A 2 S + , u 00 = 0) rotational transitions using a spectral simulation (see results discussed in Section 4.2 and shown in Fig. 9). The detailed working principle of the simulation will be described in the next section. Rotational temperatures were obtained for different water contents ranging from 0.1 to 1.3%. A standard deviation of 2.2 K shows that the gas temperature stays fairly constant within the investigated range of H 2 O admixtures. Finally, the pressure broadening is determined as Ds L = 0.37 cm À1 from an average of several automated fits to the data using the previously specified values for Ds D and Ds I . This value is in reasonable agreement with Ds L = (0.46 AE 0.03) cm À1 for He measured by Marinov et al. 59 in the 'COST Reference Microplasma Jet' using Doppler-free TALIF, albeit for a different optical transition. A typical transmission spectrum is presented in Fig. 3. We only evaluate the strongest J = 2 transition of the O(2p 4 3 P J=2 ) triplet from the fine structure split ground state to the first electronically excited state because of the low signal-to-noise ratio of the weaker J = 0, 1 transitions. In order to estimate the total ground state density n O ¼ P J¼0À2 n J , the Boltzmann factor is applied, where E J is the energy of the state and k B is the Boltzmann constant. The main uncertainties in this technique lie in the estimated absorption length (uncertainty of 5%), and the accuracy of the transition probability A ik (r3% 60 ), which are included in the expression for the absorption coefficient in eqn (4). A change of the gas temperature within 10 K influences the Boltzmann factor calculated using eqn (5) by less than 2%. We therefore estimate the systematic error in all VUV-FTAS measurements presented here to be within 10%.

UV broad-band absorption spectroscopy
Absolute OH densities are measured in the same plasma source using UV-BBAS using two different experimental setups to ensure reproducibility. The experimental setup UV-BBAS I is presented in Fig. 4(a). Light from an ultra-stable broad-band plasma lamp (Energetiq EQ-99) is guided through the middle of the plasma channel and focused on the entrance slit of a 320 mm spectrograph (Isoplane SCT320) with a 2400 grooves per mm grating. Spectra are recorded using a photodiode array detector (Hamamatsu S-3904). The setup is described in detail elsewhere. 61 The second setup (UV-BBAS II), which is shown in Fig. 4(b), comprises several different components, mainly a UV LED (UVTOP-305-FW-TO18, Roithner Lasertechnik GmbH) as light source and a CCD camera (Andor Newton 940) in combination with a spectrometer (Andor SR-500i) as detector. For the UV-BBAS II setup, the plasma is mounted on an automated x-z stage, allowing for spatially resolved measurements in the plasma channel. The experimental setup is described in detail elsewhere. 24 To calculate the absorbance in eqn (4), four signals are required: plasma on and light source on (I P,L ), plasma on only (I P ), light source on only (I L ) and a background with both plasma and light source off (I 0 ). Each signal is integrated over a time period of 50 ms, with a plasma stabilisation time of 4 s beforehand. A schematic showing this sequence is shown in Fig. 5. The plasma transmission T in eqn (4) is calculated as TðsÞ ¼ I P;L ðsÞ À I P ðsÞ I L ðsÞ À I 0 ðsÞ An example spectrum of the OH absorbance A is shown in Fig. 6. Using two setups, an admixture range of 200-13 000 ppm humidity content is investigated. Measured OH rotational absorbance spectra of the transition OH(X 2 P i , u 0 = 0) -OH(A 2 S + , u 00 = 0) are fitted using a spectral simulation in order to obtain absolute OH(X 2 P i , u 0 = 0) densities. The fitting programme is based on a calculation of the Einstein coefficients and wavelengths for the individual transitions within the investigated rotation band, as described by Dieke and Crosswhite. 62 Based on the selection rules for the total angular momentum J = L AE S and the angular momentum L (without electron spin S = 1/2 for the OH radical), relative intensities are calculated for 12 possible branches, using expressions derived by Earls. 63 An experimental value for the radiative lifetime for a rotationless upper state F 1 ( J 00 = 0.5) has been determined as 0.688 ms 64 (here, F 1 donates the doublet component of the upper state with J = L + 1/2, in accordance with Diecke and Crosswhite 62 ). Therefore, all calculated relative Einstein coefficients can be normalised to this value.
Our calculated values are in good agreement with those from Goldman and Gillis. 65 As in Dilecce et al., 45 the spectral fitting includes an instrumental function, whose width represents the spectral resolution of the spectrometer, which depends on the pixel size of the detector array, the optical grating and the width of the spectrograph's entrance slit. We assume the instrumental function to be Gaussian. Examples of measured and simulated absorbance spectra are shown in Fig. 6. Here, the instrumental width is 56 pm (UV-BBAS I) or 34 pm (UV-BBAS II), which is much larger than the Doppler (Dl D (304 K) = 0.098 cm À1 = 0.93 pm) and pressure broadening (estimated as Dl P (1 atm) = 0.07 cm À1 = 0.66 pm, as in ref. 66). The fitting programme is also used to calculate OH rotational temperatures.
The main systematic uncertainties of UV-BBAS lie in the estimation of the absorption length (5%), and the accuracy of the calculated Einstein coefficients, which we estimate here to be within 10%. For the absorbance measured with the UV-BBAS II setup (featuring the LED), the standard deviation of the noise is in the order of 3 Â 10 À4 , which places a lower limit on the measurable OH density at 3.6 Â 10 13 cm À3 . For the UV-BBAS I setup (featuring the ultra-stable light source), the noise level of the measured absorbance is typically an order of magnitude lower, and therefore disregarded in the uncertainty estimation. The combination of the systematic error and a statistical error of 7% is shown as error bars in the results that follow.

Determination of plasma power
For accurate comparison between simulation and experiment the rf power dissipated in the plasma is a particularly important. Experimentally, this so-called plasma power is measured by determining current, voltage and phase shift using current (Ion Physics Corp. CM-100-L 1 V/A) and voltage probes (PMK-14KVAC).   The probes are installed between the impedance matching unit and the plasma source. The time averaged power P is given by where U and I are the voltage and current amplitudes, respectively, and j is the phase shift between the two. Parasitic power losses, e.g. into the plasma source or the rf cables, are accounted for by measuring the power deposited in the system without a gas flow, so that the ignition of the plasma is inhibited. The subtraction method is then used 67,68 for a given current to determine the plasma power P d (I 2 ) = P on (I 2 ) À P off (I 2 ) The net power P d is the difference between the power measured with and without plasma, P on and P off , respectively. For a given plasma volume V plasma , the corresponding plasma power per unit volume p d is given by The instrumental phase shift of the measurement system (probes, BNC signaling cables, and digital oscilloscope) is determined using a variable air capacitor with known phase shift (MFJ 282-2018-1). For the calibration measurement, the plasma source with its rf cable to the matching box is replaced by this capacitor.
Current and voltage waveforms are recorded by a fast oscilloscope (LeCroy WaveSurfer10, 10 GS per s sample rate). The voltage and current amplitude as well as the corresponding phase shift are determined by a Fourier analysis of this data. P d was found to be approximately constant (within 15%) as a function of feed gas water content at a constant generator power and matching settings. The average of P d over several different water contents is used as an input for the simulations over the whole range of water content. The average value of P d was determined as 2.8 W (E14 W cm À3 ) for the UV-BBAS measurements of OH (at approximately 510 V pp ), and 2.1 W (E10 W cm À3 ) for the measurement of O using VUV-FTAS (at approximately 470 V pp ). These values are used as the input for the simulations, unless otherwise stated.
Power measurements are carried out separately from the density measurements using two power generators: coaxial RFG-150-13 (150 W maximum output power, same model as used for the OH measurements using the UV-BBAS II setup), and coaxial RFG-50-13 (50 W maximum output power, smaller range of powers for a better stability). We find a similar average power using these two setups, and that the power stays constant as a function of water content within one measurement set with a standard deviation of all points below 5%. We estimate a total uncertainty of 15% from repetitive measurements. These variations are small enough to not significantly influence measured species densities, particularly for OH, which we found to be only weakly dependent on applied voltage, and therefore power (increase of about 40% when the voltage is increased from 490 to 850 V pp , not shown here).

Model description and reaction mechanism
To better understand the dynamics of reactive species in cold atmospheric pressure plasmas, zero-dimensional plasma chemical kinetics simulations (global models) are often used. 51 In this work, experimental results are compared to those obtained using the GlobalKin code as described elsewhere. 50 GlobalKin solves the continuity equation for mass conservation for both charged and neutral species, taking into account particle production and loss through gas phase reactions and interactions with surfaces where N denotes the number density of heavy particles i, S V the surface to volume ratio, L D the diffusion length, D the diffusion coefficient, g the surface sticking coefficient, f the return fraction of species from walls, and S i the source term, which accounts for gas phase production and losses. In addition, the electron energy conservation equation is solved to calculate the electron temperature in the plasma by taking into account the balance of power input and loss of electron energy due to elastic and inelastic collisions with heavy particles d dt where n e is the number density of electrons, T e the electron temperature, m e and M i the electron and heavy particle masses, respectively, n mi the electron collision frequency, k the reaction rate coefficient and De l the electron energy gain/loss through inelastic collisions. GlobalKin also incorporates a two-term approximation Boltzmann solver, which updates the electron energy distribution function during the simulation, and calculates electron impact rate coefficients, using electron impact cross sections as an input. From the electron energy distribution function electron transport coefficients are also determined for the use in the continuity equation.
In this work we apply a temporally constant power deposition corresponding to the time averaged power measured in the experiment. For rf APPs the electron heating is strongly modulated in time, leading to a power and electron impact rate coefficients, that vary during the rf cycle. 69,70 This effect is not captured in our model. However, Lazzaroni et al., 70 investigated the differences between a conventional global model, using a time averaged power deposition, and one that takes into account time-varying power deposition within the rf period. For their case, using a He/O 2 reaction mechanism, the densities of neutral species calculated by the modified model (O, O 3 , O*, O 2 *) were typically within a factor of 2 of those from the conventional model. The trends in the results of the two models were similar. Therefore, we expect that neglecting the time-varying power deposition in our model will only lead to a quantitative difference in the results, while the trends should remain valid.
Gas temperatures are self-consistently calculated using the GlobalKin code using 50 d dt Here, N g is the gas density, c p the specific heat of the gas, T g the gas temperature, DH i the change of enthalpy for reactions with rate R i , k the thermal conductivity of the gas, and T s the surface temperature of the reactor wall. Therefore, GlobalKin balances gas heating via electron collisions (first term on right hand side), chemical reactions (second term), and heat exchange with surrounding walls (third term). Here, we assume T 0 g = 295 K (room temperature) as the initial temperature of the gas before entering the plasma channel. Coupled powers are typically small in this work, therefore it is assumed that the reactor wall is not significantly heated and T s is set to 295 K. This is in good agreement with previous observations, 24 where the electrode temperature was measured using an infrared thermometer in a very similar plasma configuration under a variation of plasma power.
The model incorporates 43 species and 390 reactions. Table 1 contains the species in the mechanism. The plasma reaction mechanism is in Appendix A (Tables 5-8). At the surfaces, it is assumed that most neutral and negatively charged species (except electrons) do not react, while positive ions are neutralised with a probability of 1. The species assumed to react differently are listed in Table 4 in the Appendix A. A detailed discussion of the role of surface interactions in a similar simulation system is given elsewhere. 71 In this work, the model is solved for a channel length of 2.4 cm, with a gas flow rate of 5 slm, which corresponds to a gas velocity of about 11 m s À1 . Using a pseudo-1D plug flow, temporally computed densities are converted into spatially dependent quantities. Where species densities are presented as a function of humidity content, densities are extracted from the simulation at the axial centre of the source (at 1.2 cm), which is the position where measurements were made.
From the plasma dimensions, the diffusion length L D , a necessary parameter for determining diffusion losses of particles, is calculated as 72 for a plasma with rectangular cross section (x Â y) and length l. For the plasma source used here, L D = 0.0316 cm. (This is larger than the 'COST Reference Microplasma Jet' with L D = 0.0225 cm.)

Pathway analysis
The PumpKin software 73 is used to identify the production and destruction pathways for the selected neutral species. The reaction pathway of a species of interest results from (a) analysing the elementary reactions that contribute directly or in subsequence to the formation of this species and selecting the significant ones only, (b) algebraically summing up the formal notations of these reactions, and (c) eliminating shorterlived species, here the electrons and some ions, to end up with a simplified net reaction. The short-lived species are defined as those with a lifetime shorter than a lifetime set by the user, which we denote as t p . Note, that this so-called net reaction should not be mistaken as an elementary reaction, since it is specific to the choice of eliminated species. This approach is particularly useful for understanding the production and destruction of species which are formed via complex reaction pathways involving a chain of elementary reactions, as opposed to simply one or two. Among the neutral species considered in this work He* and He 2 * metastables have the shortest effective lifetime. For the pathway analysis, we therefore choose t p to be slightly shorter than the lifetime of these species, which is in accordance with previous work. 52 Table 2 shows the strong dependence of the  simulated lifetime on the humidity admixture for specific plasma conditions. These findings support the conclusion of Niemi et al., 46 that the metastable character of these helium species at high pressure is significantly reduced in the presence of small admixtures or even impurity levels of molecular gases through Penning ionization under atmospheric pressure conditions.

OH densities along plasma channel
The density of OH along the plasma channel for an intermediate humidity content of 5400 ppm and a plasma power density of 18 W cm À3 is shown in Fig. 7. Experimental results show a rapid increase of the OH density over the first 2 mm of the channel. With increasing distance from the gas inlet the density stays approximately constant between 3.5 Â 10 14 cm 3 and 4.0 Â 10 14 cm 3 up to the end of the channel. A similar trend is also observed in the simulation. Absolute simulated and experimental densities agree well within around 25%, which is likely within the combined uncertainties of experimental data (as shown as error bars in Fig. 7) and simulations. Uncertainties in the simulation results would most likely occur due to uncertainties in used reaction rate coefficients, and considered reaction pathways, and were shown to be within a factor 10 for a He/O 2 reaction mechanism under similar plasma conditions. 74 To gain insight into the dynamics of OH formation, a pathway analysis is performed for the three regions highlighted in Fig. 7, which correspond to the fast build-up of OH at the entrance of the plasma channel (0-0.2 cm), a steady-state region (2-2.5 cm) and the decay of OH in the plasma effluent (3.3-3.5 cm). The dominant production and consumption pathways for OH, averaged over each region, are shown in Fig. 8.
At the entrance of the discharge channel (0-0.2 cm), the gas consists mainly of the initial feed gas mixture plus some rapidly forming species such as ions and electrons. Therefore, the main production reaction for OH is through electron impact with water vapor, either via dissociation or dissociative attachment e + H 2 O -OH + H + e (60%) (14) e + H 2 O -OH + H + e (9%) The products of eqn (15) are ground state atomic hydrogen and OH in its excited OH(A) state, while the products of eqn (14) are both in their ground states. The percentage contribution for each reaction to the total production of OH is shown in brackets.  Another production mechanism for OH is through the formation, and subsequent destruction, of charged water clusters, as was previously identified by Ding and Lieberman. 52 The formation of these clusters is typically a multi-step process. For positive clusters, this process usually begins through ionisation of H 2 O, either through electron impact or Penning ionisation with He*. These water ions then collide with water molecules to form the cluster ion H + (H 2 O), which accumulates additional water molecules through a series of reactions: Here, the reaction below the solid line represents the net reaction. Similar processes occur for negative ion clusters OH À (H 2 O) n , which are included in the reaction mechanism for n r 3.
The main consumption pathway for OH in the first 0.2 cm of the jet is the formation of H 2 O 2 , and recombination to water through The rapid increase of OH density within the first two millimeters raises the question, if averaging the pathways over this region is a valid analysis. For both the production and consumption pathways, the contribution of each reaction does not change significantly, if evaluated for separate points within the first 0.2 cm instead of averaging over this region. However, the ratio of the total rate of production and consumption changes significantly, leading to the increase in OH density over this region. Further from the gas inlet, the rates of production and consumption equalise leading to an equilibrium OH density.
In the quasi steady state region (2-2.5 cm), the previous pathways still dominate. However, additional species with intermediate lifetimes build up along the channel and begin to play a role in the formation of OH. For example, hydroperoxy radicals (HO 2 ) promote production of OH through reactions with H In the afterglow region (3.3-3.5 cm), a rapid decay of OH occurs both in experiment and simulation, as shown in Fig. 7. Short lived species such as ions and electrons recombine rapidly, while metastable species like He* and He 2 * are consumed through Penning ionization with water before reaching this region. Therefore, the chemistry in the plasma effluent is dominated by intermediate and long lived neutral species, where OH is produced mainly through reactions between H and longer-lived neutral species: In this region, consumption occurs at a higher rate than production, leading to a decrease in the OH density, with reactions (22) and (18) (collisions with H 2 O 2 and OH) dominating.

OH densities under varying humidity content
The density of OH measured by UV-BBAS in the centre of the plasma channel (at position 1.2 cm) as a function of the H 2 O content in the feed gas is shown in Fig. 9(a). The OH density increases sub-linearly with increasing H 2 O content, as previously observed. 34,43 Absolute densities obtained in this work also agree well with results obtained by others. 34,43 Absolute OH densities measured with the two different experimental setups agree well within the uncertainties in each measurement. The simulated OH densities at different feed gas humidity contents are shown in Fig. 9(a). In general, good agreement in the trends of experimental and simulation results is observed. Absolute OH densities agree particularly well at low H 2 O contents o2000 ppm. Towards higher H 2 O contents, simulated densities are higher than those measured experimentally. The largest difference is a factor 1.8 at the highest H 2 O content, which is reasonable agreement given the previously mentioned uncertainties.
OH(X) rotational temperatures and gas temperatures calculated using GlobalKin are shown in Fig. 9(b), and found to be in good quantitative agreement with each other. In both experiment and simulations, temperatures stay fairly constant with increasing water content. While in the simulation, a very small decrease of the gas temperature is observed, the experimental data is more scattered, and a clear trend cannot be observed taking into account the uncertainties of the measurement.
The main production and consumption pathways for OH at different H 2 O admixtures are shown in Fig. 10 (17)). At any stage of the clustering process, the clusters can be destroyed by dissociative recombination with electrons Towards higher H 2 O contents, this pathway is gradually replaced by direct electron impact dissociation or dissociative electron attachment with H 2 O (eqn (14) and (16)). OH is mainly consumed by reactions with other OH radicals (eqn (18) and (19)) and O Towards higher water admixtures, the contributions of these reactions to the consumption of OH decrease slightly, and reactions of OH with H, and more slowly forming species such as H 2 O 2 and HO 2 become more important.
In both experiment and simulation, OH densities increase rapidly with increasing H 2 O at low H 2 O content, and less rapidly at high H 2 O content. The transition between these two regimes occurs at lower H 2 O content (around 2000 ppm) in the experiment compared to the simulation (around 3000 ppm). This leads to the increasing discrepancy between simulation and experiment at higher H 2 O contents where the experimental OH densities saturate and the simulated OH densities continue to slowly increase. The reason for this transition is investigated by looking at the most important formation pathways for OH, which are production by electron impact dissociation and dissociative attachment of H 2 O (eqn (14)-(16)), and consumption via reactions with OH to form H 2 O 2 and H 2 O (eqn (18) and (19)). As the gas temperature remains relatively constant with changing water content, rate coefficients for consumption of OH also stay approximately constant. The reaction rate coefficient for production of OH is dependent on the electron temperature T e , and the electron density n e . Fig. 11(a) shows the two quantities as a function of humidity content. T e and n e show opposite trends with increasing H 2 O admixture. T e , which is calculated from balancing the electron energy sources and losses (see eqn (11)), increases with increasing H 2 O content due to increasing electron energy losses in inelastic collisions with water molecules. For constant power input, the increased electron energy losses and T e are balanced by a decrease in n e with increasing H 2 O content.
The effect of these changes on the total rate coefficient for dissociation k diss = k 14 + k 15 + k 16 and the dissociation frequency R = k diss n e is shown in Fig. 11(b). Due to the variation in T e , k diss increases with increasing humidity content, exhibiting a  View Article Online similar trend to the OH densities shown in Fig. 9. Here, the transition from a fast to a slow increase also occurs around 3000 ppm. The dissociation frequency R exhibits a peak at this humidity content, which represents the optimum between the increasing T e and decreasing n e . Thus, the origin of the transition between fast and slow increase in OH density with increasing humidity content is a result of the transition between an increasing dissociation frequency below 3000 ppm H 2 O to a decreasing dissociation frequency above 3000 ppm H 2 O. Overall, the dissociation rate R Â n H 2 O , and therefore the OH density, increases over the whole range due to the increasing value of n H 2 O . Based on this discussion, the differing H 2 O contents at which the transition occurs in the experiment and simulation may indicate that the rate of electron energy loss with increasing H 2 O content is misrepresented in the simulation. Another reason for the discrepancy between the simulated and measured trend in OH densities at higher water contents might be due to an additional consumption mechanism for OH, which is not taken into account in this work, such as the population of vibrationally excited states, which would also scale with T e .

O densities as a function of humidity content
The absolute O density measured by VUV-FTAS in the centre of the discharge (at x = 1.2 cm, triangles) as a function of H 2 O content in the feed gas is shown in Fig. 12 Simulated O densities are around a factor 2 lower than those obtained experimentally. The measured densities agree well with previous measurements in a similar source using twophoton absorption laser induced fluorescence. 34 A possible explanation for the difference in absolute O densities and trends between the experiment and simulation may be limitations in the global model, particularly the accuracy of the rate coefficients used, as discussed earlier. O is not directly produced from H 2 O due to electron or heavy particle impact dissociation in significant amounts at the electron temperature of interest. As a result, O must be formed in a process taking at least two steps, meaning that the uncertainties in multiple rate coefficients will play a role in determining the uncertainty in the simulated O density. As a result, the simulated O density is likely to have a larger uncertainty than the simulated OH density, whose dominant formation occurs directly from electron collisions with H 2 O. As shown in Fig. 13, the dominant production mechanism of O is via recombination of two OH molecules to form H 2 O and O. At lower H 2 O Fig. 11 (a) Electron density n e and temperature T e , and (b) combined rate coefficient k diss for electron impact dissociation and dissociative attachment and dissociation frequency R = k diss n e . Conditions are 5 slm total He flow and 14 W cm À3 plasma power. contents, O is also formed through processes involving positive ion water clusters: With increasing H 2 O admixture the formation of O 2 is also increased. As a result, electron impact dissociation of O 2 becomes a more important production pathway for O: Where the numbers in brackets represent the electron energy thresholds for these reactions. O is mainly consumed by reactions involving OH forming O 2 and H (eqn (27)).
The measured and simulated O densities show an increased discrepancy towards smaller H 2 O admixtures o1000 ppm. A possible explanation for this might lie in the presence of unintentional air impurities in the experiment, which have been found previously to be able to influence the chemical kinetics in atmospheric pressure plasmas. [75][76][77] For the measurement of O, we use helium with a purity level of 99.999%, whereas the main impurities are H 2 O (3 ppm) and O 2 (2 ppm). Additional small impurities could arise from residual gases in the feed gas line. Simulations for two different non-zero O 2 impurity concentrations in the order of typical O 2 impurities originating from the feed gas supply are shown in Fig. 12

Numerical investigation of the production of longer-lived species
The OH density reaches a steady-state value in the simulation well before the end of the plasma channel in both simulation and experiment. Particularly at higher H 2 O content, this is a result of OH being primarily produced by direct electron impact dissociation of H 2 O in a one step process (eqn (14)) and is consumed in interactions with other OH molecules. Atomic hydrogen behaves similarly, being produced mainly by electron impact dissociation and consumed via interactions with surfaces, 71 both of which occur relatively quickly. However, other species do not reach a steady state within the length of the plasma channel, and instead continuously increase in density up to the outlet of the plasma source. This is particularly true for slowly forming, long-lived species such as O 2 , H 2 O 2 and H 2 , as shown in Fig. 14. This finding suggests that the length of the plasma source, or the gas flow rate, and therefore the residence time of the gas, can be used to control the ratio of different species densities by taking advantage of the different timescales required for them to reach steady-state.
First, we will discuss the formation of O in more detail. The density of O does not reach a steady-state value in the simulation within the plasma channel for most investigated conditions using a He-H 2 O gas mixture. Long timescales for simulations of atmospheric pressure He-H 2 O plasmas to reach steady-state have also been found by others. 78 This is in contrast to the case where similar sources are operated in He-O 2 mixtures. 79 In the work described in ref. 79, O densities approach steady-state towards the end of the plasma channel of the AAPPJ. In Fig. 14, O densities are increasing sharply within the first few millimetres of the channel, and then at a lower rate up to the end of the channel. Therefore, O densities follow a similar dependence as the OH densities also shown in Fig. 14. This is not surprising when considering that both the dominant production and consumption pathways are related to OH, i.e. production by reactions of two OH molecules to form H 2 O and O  (27)). The fact that O is still building up within the channel, while OH approaches a steady-state value, is due to the continuous build-up of O 2 in the channel, also shown in Fig. 14. Electron impact dissociation of O 2 (eqn (29) and (30)) provides an additional formation mechanism for O further into the channel, although eqn (19) and (27) are still the dominant production and consumption pathways for O. Overall, this leads to a slow increase of the O density while the O 2 density continues to increase.
The formation of species that reach steady-state on timescales longer than the residence time in the discharge channel are usually comprised of a complex multi-step processes. As an example, we demonstrate the dominant pathways for formation of O 2 , which is an important precursor for the formation of excited states of O 2 , such as O 2 (a 1 D). Since O 2 is a slowly forming species, we look at dominant production and consumption pathways for a longer timescale t p than those previously given in Table 2. The time scale of interest in the simulation is chosen so that only He, H 2 O, O 2 , O 2 (a 1 D), H 2 , and H 2 O 2 are treated as long-lived species, in accordance with previous studies. 52 The computational lifetimes of the shortest-lived species of these six are listed in Table 3 for different H 2 O contents.
The two main net production reactions for the formation of molecular oxygen are found to be Many different pathways are possible in order to obtain these net reactions. A few examples of these pathways are Note that eqn (33) and (34)  The cluster ions produced are consumed by dissociative recombination with electrons, and H, which is formed in that process, is lost by surface recombination. Other pathways exists that involve the formation of clusters, but are not explicitly discussed here.  Eqn (35) has a different net production reaction than the others. In this case, OH produced from electron impact dissociation reacts with H 2 O 2 to form reactive HO 2 , which then forms O 2 during reactions with OH.

Conclusions
In this work, the chemical kinetics in an rf atmospheric pressure plasma with humidity are investigated using experimental and numerical techniques. By using this combination, computed species densities are benchmarked against experimental densities. The simulations are then used to reveal the dominant formation pathways of species of interest, here OH and O, and longer lived species such as O 2 , which is an important precursor for the formation of its excited states, such as O 2 (a 1 D). This work provides a detailed understanding of chemical kinetics in the active plasma. In many applications, reactive species will exit the plasma source and transit into an effluent region where they will mix with ambient air, where their chemical kinetics will differ. While this is not considered in this work, the results presented here provide a basis to be built on in future work to understand reactive species kinetics in this transition region.
Absolute number densities of O and OH are determined experimentally using VUV high-resolution Fourier-transform absorption spectroscopy, UV broad-band absorption spectroscopy, and numerically by using the 0-D plasma chemical kinetics code GlobalKin.
Absolute OH densities and formation pathways are investigated as a function of position in the discharge. Three different regions can be identified i.e. (a) a strong increase of OH density in the first few millimeters of the plasma channel, (b) a quasi steady-state region, and (c) a rapid drop of OH density in the plasma effluent region. During the fast increase and steady-state regions, OH is mainly produced via fast processes such as electron impact dissociation of H 2 O, and consumed predominantly via reactions with other OH molecules to form H 2 O 2 or H 2 O. These relatively simple chemical kinetics make it possible for OH to reach an equilibrium value within the plasma channel.
Other species, whose densities have not been measured, are investigated numerically as a function of position in the plasma channel. Simulation results show that the H density approaches a steady-state value within the plasma channel, similarly to OH as discussed previously, as it is mostly formed directly via electron impact dissociation of water, and consumed at surfaces to form stable H 2 . However, most other species generated in the He-H 2 O plasma studied in this work do not reach a steady-state value within the length of the plasma channel due to more complex formation mechanisms. This has been shown using O 2 as an example. Therefore, the length of the plasma source could be used as a control parameter to tune the chemical composition of the gas at the end of the plasma jet for applications.
Both OH and O densities are also investigated as a function of the humidity content in the He feed gas. It is found, both in experiments and simulations, that O and OH densities increase non-linearly with increasing feed gas humidity, offering the possibility of tailoring reactive species densities by changing the feed gas composition.
The maximum OH density is on the order of 3-4 Â 10 14 cm À3 (13-17 ppm). It is found that at very low water content, OH is mainly produced via reactions between H 2 O + and water molecules to form OH and protonated water clusters of the form H + Á (H 2 O) n , while electron impact dissociation of H 2 O becomes an increasingly important production pathway with increasing water content. The main loss channel for OH at all H 2 O contents is recombination to form H 2 O 2 .
The maximum O density on the other hand is found to be in the order of 3 Â 10 13 cm À3 (1.3 ppm). Recombination of two OH molecules is the most important production process for O at all H 2 O contents, while at very low water content, OH is also strongly produced via reactions between OH + and water molecules to form O and protonated water clusters. Since the dominant destruction pathway of O is recombination with OH to form O 2 and H, the formation of O is strongly coupled to the OH density in the gas flow. At higher H 2 O concentrations, electron impact dissociation of accumulated O 2 can also contribute to the production of O. It is also found that towards low H 2 O content, production of O from air impurities in the ppm range originating from the feed gas can increase the O density via direct electron impact dissociation of O 2 . Towards higher H 2 O admixtures, this effect becomes less significant due to increased production via collisions involving OH. Therefore, larger amounts of purposely admixed molecules lead to a better control of the plasma properties and reactive species than operating the source with small or no intentional admixtures.

Conflicts of interest
There are no conflicts to declare. Table 4 shows wall recombination coefficients and return species for the simulations used in this work.  Table 5 shows electron impact reactions used in this work. Reaction rate coefficients are either taken from the literature, or calculated by the GlobalKin two-term Boltzmann equation solver.

A. Reaction mechanism
For the latter, reaction rate coefficients are indicated as f (E), and collisional cross sections are taken from the indicated literature. Electron impact cross-sections are taken from several databases, for    84 Although not all of the reactions from these databases are included in the plasma-chemical reaction mechanism shown here, they are still accounted for in the Boltzmann solver calculation for the electron energy distribution function and electron transport coefficients. Any other approach to obtain reaction rate coefficients is denoted by footnotes. Table 6 shows reaction rate coefficients for ion-ion recombination processes. It is generally known that ion-ion recombination processes can occur both as two-or three-body processes, depending on the gas pressure. Two-body reaction rate coefficients for several different gases have been obtained in ref. 85 at low pressure, and found to be in the order of 10 À13 m 3 s À1 or lower. Taking into account the He density at atmospheric pressure and at 315 K, and the rate coefficient for three-body ion-ion recombination proposed by Kossyi, 86 the effective two-body reaction rate amounts to a value in the order of 10 À12 m 3 s À1 . Due to this higher effective rate coefficient under our conditions, we only include three-body ion-ion recombination rate coefficients in this work. These reactions are found to be particularly important for the destruction of the higher-mass water clusters, which are abundant at higher H 2 O admixtures. Similar observations have been made by Liu et al., 48 who, after an analysis of the robustness of their chemistry model, only included a few ion-ion recombination reactions in their simplified models, a large fraction of which were three-body recombination processes for the collisions of higher mass cluster ions. We also found that under our conditions, ion-ion recombination between positive He ions and negative ions are negligible due to the rapidly decreasing He ion density with increasing water content, which again is in accordance with the findings of Liu et al., 48 and the fact that He ions undergo charge exchange reactions with most neutral species due to their high ionisation potential. Table 7 shows reaction rate coefficients for collisions between ions and neutrals. In this table a number of three-body processes are included. Three-body processes are typically characterised by a pressure dependence. The nature of these reactions mean that this pressure dependence normally takes the form of a curve exhibiting low-and high-pressure limits. In the low-pressure limit the effective rate coefficient (i.e. the three-body rate coefficient multiplied by the   third body density) is linear with the third body density. In the high pressure limit the effective rate coefficient is independent of the density of the third body. In the region between the two limits the effective rate coefficient is non-linear with the third body density.
For a number of reactions, this transition region occurs around atmospheric pressure, therefore effective rate coefficients must be calculated using available knowledge of the high and low pressure limits. The coefficients which have been explicitly calculated for  87 using an exponential fit, and using constant values n = 16, B = 5000, and k L = 10 À24 . f Value is listed as a lower limit in reference. g Estimated branching ratio. h Third body is H 2 O in reference. atmospheric pressure are marked as ''effective'' in Tables 7 and  8. Among these reactions is the formation and destruction of protonated water clusters H + Á(H 2 O) n , where the rate coefficients are given by Sieck et al. 87 Here, the expressions given in ref. 87 are used to calculate the effective rate coefficients for these reactions under our plasma operating conditions (atmospheric pressure, T g = 280-350 K). Results are fitted with an Arrhenius expression, where possible, in order to keep the temperature dependence for these reactions, since the formation of cluster ions is highly temperature dependent. 87 The rate coefficients for the formation of the two highest order clusters taken into account in this work are estimated by extrapolating the coefficients k 300 0 and A given by Sieck et al. using an exponential fit, and using constant values n = 16, B = 5000, and k L = 10 À24 (see Sieck et al. 87 for further description of these coefficients). Table 8 shows reaction rate coefficients for collisions between ions and neutrals. Similar to the description of ion-neutral reactions in Table 7, a number reactions rates in Table 8, including neutral reactions of several O and H containing species, are specified as ''effective''. Data to calculate rate coefficients for these reactions has generally been taken from the IUPAC chemical kinetics database. 88 For the calculation of ''effective'' decay rates, and generally for three-body processes, we have multiplied the three-body rate coefficient from the respective sources by a background gas dependent efficiency factor, if available, where the collider gas is different from He in the reference. This accounts for the fact that He is a less effective quencher compared to other to the more abundant protonated water clusters under all investigated conditions, and the fact that in this work, the focus lies on the investigation of the neutral particle dynamics.

B. List of equipment
The equipment used for the investigations in this work is listed in Table 9. 88 l a In m 3 s À1 and m 6 s À1 for two-and three-body collisions, respectively. b Value in an upper limit in reference. c Estimated value in reference. d Estimated branching ratio. e Branching ratios taken from Sanders. 168 f Third body is Ar instead of He in reference. The gas efficiency factor is assumed to be 1. g Third body is Ar instead of He in reference. The gas efficiency factor is assumed to be 0.65. This factor is calculated by dividing reaction rate coefficients for He and Ar as background gases for the same reaction measured by Zellner et al. 169 h Effective rate coefficients calculated from pressure dependent rates for 1 atm and fitted by an Arrhenius expression in the temperature range 280-350 K. i Third body is N 2 instead of He in reference. The gas efficiency factor is assumed to be 0.43. This factor is calculated by dividing reaction rate coefficients for He and N 2 as background gases for the same reaction measured by Hsu et al. 170 j Recommended rate coefficient in reference is for N 2 background gas instead of He. We apply a gas efficiency factor of 0.41 to the low-pressure limit reaction rate coefficient to account for this. This factor is calculated by dividing the room temperature rate coefficient from the given reference for He background gas (measured by Forster et al. 171 ) by the recommended value (measured by Fulle et al. 172 ). k Third body is Ar instead of He in reference. The gas efficiency factor is assumed to be 0.77. This factor is calculated by dividing reaction rate coefficients for He and Ar as background gases for the same reaction measured by Campbell and Thrush. 169 l Third body is N 2 instead of He in reference. The gas efficiency factor is assumed to be 0.61. This factor is calculated by dividing reaction rate coefficients for He and N 2 as background gases for the same reaction measured by Lin and Leu. 173