Vacuum ultraviolet radiation in and from an atmospheric pressure plasma source

This work presents a computational study of the atmospheric pressure plasma of the COST jet source in He + 0.1% N2 buffer gas, incorporating the treatment of the VUV resonant radiation from He atoms at a wavelength of 58.4 nm. The plasma is driven by single- and multi-frequency waveforms at various peak-to-peak voltages between 400 V and 685 V. The simulations, which include the photoemission of electrons from the electrodes as well as the photoionization of the N2 molecules indicate that these processes, induced by the VUV radiation, have relatively little influence on the discharge characteristics due to other more efficient charged particle production mechanisms, e.g., Penning ionization. While the plasma characteristics are computed along one spatial dimension perpendicular to the electrodes, the tracing of the VUV photons is executed in the real 3D geometry of the plasma source to allow the computation of the VUV photon flux leaving the device through the orifice. Photon fluxes up to ≈5×1015 cm−2 s−1, corresponding to an energy flux of ≈18 mW cm−2 are found at multi-frequency operation.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. an inert gas (He or Ar) and a molecular gas or molecular gases (N 2 , O 2 , H 2 O). Due to the relatively low power consumption of these sources the heating of the gas is normally insignificant making these devices applicable for the treatment of temperature-sensitive substances, including biological tissues.
The applications of APPJs [23] are mainly based on the efficient production of reactive oxygen and nitrogen species (RONS) in the gas mixture plasma, see, e.g. [24][25][26], therefore optimisation of the conditions (e.g., gas composition, driving voltage amplitude and waveform) at which a high production rate of RONS can be achieved, is essential [27][28][29]. This requires both experimental investigations that reveal the charged particle and excitation/ionization dynamics and provide information about the densities of radicals, as well as modelling of the discharge plasma. In the past years, a series of combined experimental and particle simulation studies [28][29][30] has been conducted to characterise the physics and chemistry of the COST reference jet [31] that is a standardized type APPJ source (COST = European Cooperation in Science and Technology). These computational studies have used a simplified plasma chemistry model as compared to detailed chemical kinetics modelling investigations [32,33], but contributed significantly to the understanding of different electron heating modes, mode transitions, and the generation of radicals under simple and tailored voltage waveform (TVW) radio frequency (RF) excitation of the discharges. These studies [28][29][30], and many others, have however, disregarded the effect of the presence of the intensive short-wavelength radiation in the plasma, on the discharge characteristics.
At typical conditions in APPJs, employing a high concentration of the 'carrier' gas (Ar or He) and a low concentration of molecular gas(es) a significant portion of the electron-impact excitation processes populates the lower-lying excited states of the noble gas atoms. These states are either metastable states (from which dipole radiation is forbidden to the ground state, but their presence contributes to the sustainment of the plasma via Penning ionization) or resonant states (from which shortwavelength, vacuum ultraviolet (VUV) radiation is emitted). In case of Ar, the VUV resonant lines have wavelengths around λ ≈ 106 nm, while for He, λ ≈ 58 nm. Such short-wavelength radiation can have an effect both on the discharge physics and, if emitted, can be either the basis of certain applications, while it can deteriorate others. For the former case (i.e. when the VUV radiation serves as the basis of a discharge application), sterilization with high-frequency discharges and their afterglows [34,35], the decomposition of various hazardous organic molecules [36], reactive oxygen species generation in bio-relevant liquids [37], and the creation of excited OH radicals at the conditions applicable to plasma-based cancer treatment [38] can be quoted. As a negative example, the VUV radiation damage of wafers in inductively coupled plasmas (quantified by measuring the hole currents generated in SiO 2 films) [39] is mentioned.
As the energy of the VUV resonance photons matches the energy difference between an excited state and the ground state of the noble gas atoms, this radiation can be efficiently reabsorbed by the ground-state atoms that are present with a high density at elevated pressures. Absorption events are followed by emissions, further absorptions, etc, leading to the trapping of the radiation. This peculiar way of transport has been investigated in a number of studies, of which only a few can be mentioned here. In [40], the VUV radiation spectrum from a He discharge operated at 40 Torr between plane-parallel electrodes was measured. The discharge was created by ∼kV pulses with μs duration. Strong emission on the 2 1 P → 1 1 S resonant transition was observed, with a time delay in the order of 100 ns, which can be explained by the long escape time of the VUV photons due to their significant trapping. Accompanying simulations have also considered the formation of He * 2 excimers and demonstrated the contribution of the photoionization of N 2 molecules to the propagation of ionization waves in a higher-pressure (250 Torr) He-N 2 discharge. Photoionization by the resonant photons was also identified to be an important factor in the propagation of ionization waves in the investigations reported in [41]. In [42], the importance of the electron photoemission from the discharge tube wall on the propagation of the ionization wavefront was pointed out. It is interesting to note that the VUV radiation in the plasma may even affect the operation of probes via influencing their surface temperature and thermal balance [43].
The VUV radiation transport in an Ar plasma jet has been investigated in [44]. An accurate solution to this problem was found to be essential for a correct prediction of the spatial distribution of the excited atoms. The effect of the outflow of the excited atoms in the effluent was also considered in this work. A detailed spectroscopic investigation of an APPJ (COST plasma jet) in the VUV domain was reported in [45]. The source was operated in He or Ar or He/Ar mixture, with small amounts of molecular gas admixtures. The measurements were conducted with a He-filled VUV spectrometer, which precluded observation of the resonance radiation of He, but allowed the observation of strong radiation from He * 2 excimers. The short wavelength radiation from an Ar plasma jet has been investigated in [46]. The effects of different wavelength domains were probed by applying appropriate filters for various wavelength ranges. This study concluded that the UV/VUV radiation from the plasma may significantly influence the plasma-liquid chemistry in the case of aqueous targets. Strong radiation originating from the Ar * 2 excimer as well as from N, H, and O atomic lines was observed in the study of an argon plasma jet [47]. As in this work the plasma radiation was coupled to the spectrometer via a MgF 2 window (with a cut off wavelength of 115 nm) the resonant lines were not possible to measure. For the case of biological targets (e.g. skin) the 'hard' VUV radiation was found to be blocked by biomolecular layers, supporting the safe use of such devices in human treatment [48]. A recent work [49] presented direct measurements of the sub-110 nm radiation from an atmospheric pressure plasma jet. In the experiments, the discharges were established in pure He and Ar gases, as well as in mixtures including small amounts of oxygen or nitrogen. Besides the strong resonant lines of He and Ar, intensive spectral features of nitrogen and oxygen also showed up in the measured spectra even in the case of operation with pure buffer gases due to the contact of the plasma with the ambient air. In the same setup, the VUV spectra of air discharges have also been recorded. The works cited above have provided evidence that in some cases the VUV radiation may have a significant effect on the discharge properties.
In this work, (i) we investigate the importance of the processes driven by the resonant VUV photons of He (photoionization and photoelectron emission) in the discharge in the COST APPJ and (ii) compute the flux of the VUV photons leaving this plasma source via its orifice, at different operating conditions. These investigations are based on an particle-incell/Monte Carlo collisions (PIC/MCC) code into which the tracing of the resonant VUV photons is incorporated. The paper is organised as follows. In section 2.

Model geometry
Here, we consider the COST plasma jet [31], of which the simplified model geometry is shown in figure 1. The cross section of the plasma source is H × D = 1 mm × 1 mm, and the length of the channel is L = 30 mm. The electrodes are situated at the top and the bottom of the cell, while from both the left and right sides the plasma is confined by dielectric (quartz) plates. The plasma is created by a RF voltage applied at the powered electrode. The quartz plates allow easy access to the system for optical diagnostics with a spatial resolution both along the direction of the gas flow and along the direction between the electrodes. The feed gas (in our case, He with a small amount of N 2 at atmospheric pressure) flows through the system (usually at a rate of 1 standard liter per minute in the experiments). This flow of the gas within the plasma channel does not influence the plasma characteristics, however, carries radicals in the effluent outside the plasma source. On the way towards the sample to be treated further chemical reactions between these radicals and the constituents of the ambient atmosphere can also take place, which are, however, not modelled here.
In our model, the x axis is oriented perpendicular to the electrodes (situated at x = 0 and x = H), the y axis is parallel with the plasma channel, with y = 0 being fixed at the position of the orifice (exit plane). The z axis is perpendicular to the quartz plates (situated at z = 0 and z = D). In the simulations of the discharge, variation of the plasma characteristics along the x coordinate only is accounted for in the 1D3V PIC/MCC code. The y and z directions are not resolved, the system is considered to be infinite along these axes. This simplification limits the accuracy of the computations, e.g., due to neglecting the formation of sheaths at the quartz plates. Nonetheless, earlier works based on the same approach (see, e.g., [28][29][30]) resulted in computational results being in good agreement with experimental data, regarding a number of discharge characteristics, as mentioned in section 1. When computing the probability of the escape of VUV photons from the system, the full 3D geometry of the source as shown in figure 1 is considered on the other hand. This approach is necessary to provide estimates for the escaping photon flux, as in the case of a 1D treatment the photons could hit only the electrodes. More details will be given below.

PIC/MCC simulation of the plasma including VUV photons
The PIC/MCC method belongs to the particle-mesh techniques. While charged particles move in space, their densities as well as the electric potential distribution are computed at points of a grid. The potential (in electrostatic simulations) is defined via the Poisson equation that takes into account the space charge generated by the charged particles (electrons and ions), as well as the potentials at the electrodes, which serve as boundary conditions. The densities of the charged particles are computed, the Poisson equation is solved, and particles are moved after short time steps according to the actual value of the electric field at their positions. In the meantime, particles may reach the electrodes, where various processes may take place as defined by the surface model, and particles may as well collide within the gas phase with the atoms/molecules of the background gas. The simulations work with superparticles this way significantly reducing the number of particles to be traced. The ratio of the number of real particles and computational particles (i.e., the superparticles) is called the superparticle weight, W, for more details, see, e.g. [50].
The charged particles considered in the simulations are electrons, He + ions, He + 2 ions, N + 2 ions, and N + 4 ions. The 1D3V (i.e. spatially one-dimensional) PIC/MCC code employed here is the same as previously used in studies of plasma jets [28,30], with the extensions that (i) N + 4 ions (formed via a conversion process from N + 2 ions) are also considered as charged particles, (ii) their two body recombination with electrons (which is believed to be the dominant recombination process based on the rates given in [51]) is also taken into account, and (iii) the tracing of the He resonant VUV photons is incorporated into the simulation code, including the photoemission of electrons from the electrodes and the photoionization of N 2 molecules. According to the spatially 1D approach, the finite extent of the system in the y and z directions (see figure 1) is not taken into account, except for the calculation of the transport of the VUV photons.
For electron-He atom collisions and electron-N 2 molecule collisions we use the same set of processes and corresponding cross sections as in [28,30,52]. The collisions with He atoms include elastic scattering, excitation, as well as ionization. The electron-N 2 molecule collisions include elastic scattering, rotational, vibrational and electronic excitation, and ionization. All types of electron-atom/molecule collisions are assumed to result in isotropic scattering. Correspondingly, the elastic momentum transfer cross section is used for elastic collisions. The heavy particle processes (that involve excited states and ionic species), recombination, and radiative processes are listed in table 1.
Reaction 1 (Penning ionization) is an essential ionization channel for the sustainment of the plasma. Besides this reaction, He + and N + 2 ions are also created in electron-impact VUV photon + N 2 → N + 2 + e − Photoionization [58] ionization processes. These ions participate in three-body conversion reactions, which convert He + to He + 2 (process 3) and N + 2 to N + 4 (process 6). These reactions are rather significant at atmospheric pressure and largely determine the composition of the plasma. As the main constituent of the background gas is He, in the transport of the various ionic species He atoms are considered only as collision targets. For He + ions, the approach of [54] is followed and the elastic scattering is approximated with contributions of an isotropic channel and a backward scattering (charge exchange) channel. For the other ions, He + 2 , N + 2 , and N + 4 (processes 4, 5, and 7, respectively) the Langevin cross sections for the elastic collisions with He atoms are computed.
The density of the plasma is constrained by two-body recombination, process 8, which has a rate coefficient k rec = 2 × 10 −6 (300/T e ) 1/2 cm −3 s −1 , with the electron temperature, T e given in units of Kelvins. This process has been identified to be more significant than the three-body reactions assisted by electrons and heavy neutrals [51].
In the model, half of the excitation events to singlet He levels is assumed to result in an atom in the 2 1 P state either directly or via cascades from higher singlet levels [52]. This state decays to the He ground state via the emission of a resonant photon that can be reabsorbed in the gas (processes 9 & 10). Due to their high energy, VUV photons can also ionize the N 2 component of the background gas: the cross section for this process is ≈2.5 × 10 −17 cm 2 at the wavelength of He resonant radiation [58]. In the plasma source considered here nitrogen is present in the buffer gas with a concentration of 0.1%, i.e. the density of N 2 molecules is 2.4 × 10 16 cm −3 , resulting in a photoionization mean free path of λ pi ≈ 1.67 cm. While this value is much greater than the electrode separation, due to the 'diffusive' motion of the VUV photons their path can be in the range of λ pi , i.e., process 11 (of table 1) can be important in the present system. More information about the treatment of the radiation will be given later.
Time is discretised in the simulation and between consecutive collisions the charged particles move along trajectories that result from the solution of their equations of motion, using the leapfrog scheme. The collision probability after each Δt time step is computed according to where ν proj,tot is the total collision frequency of the given projectile. As mentioned above, for the ion species a single type of target species (He) is considered. Thus, for the He + , He + 2 , N + 2 , and N + 4 ions: where n He is the He atom density, σ proj,He is the total cross section for collisions between the given ion and He atoms, ε com is the centre-of-mass energy of the colliding particles, and g is their relative velocity. In the case of ionic projectile species we take into account the thermal motion of the background gas particles via choosing random partners from the Maxwell-Boltzmann distribution of the background atoms in possible collision events. (As the flow of the buffer gas is much slower than the thermal motion a zero mean velocity is assumed when sampling from this distribution.) For electrons, we use the cold gas approximation, i.e., take g to be equal to the projectile (electron) velocity. The electrons can, however, collide with both He and N 2 particles of the background gas, i.e., for the electrons where n N 2 is the density of N 2 molecules, ε e and v e are, respectively, the energy and the velocity of the electron. In the simulations, a collision is executed when the condition R 01 < P c is met after a Δt time step, where R 01 is a random number having a uniform distribution over the unit interval. (Note, that repeated occurrences of R 01 below always correspond to independent random numbers.) In order to achieve accurate results, the parameter settings in PIC/MCC simulations should obey the stability and accuracy conditions, see e.g., [59,60]. These conditions include (i) the need for a high temporal resolution (small enough Δt) to resolve plasma oscillations, (ii) the requirement for a high enough spatial resolution, to resolve the Debye length, (iii) the Courant-Friedrichs-Lewy condition that requires that particles should not fly a longer distance during a time step than the grid division, and (iv) the need to minimise 'missed' collisions by choosing time steps sufficiently small to keep the collision probabilities at small values (P c 0.1). The first of these conditions depends on the time step for electrons, while the second on the number of the divisions of the spatial grid. Both of these together control the fulfillment of the third condition, while the last one depends on the time steps (at given cross sections). For the system considered here, the most demanding condition on the time steps, especially in the case of the electrons, is imposed by the last condition, due to the high density of the background gas. To ensure accurate computations, the electron time step is set to Δt e = 4.5 × 10 −14 s. In order to optimize the run time of the simulations we use different time steps for the ion species: Δt He + = 10Δt e and Δt He + = 100Δt e . The fulfillment of the accuracy and stability conditions ((i)-(iii), above) requires a spatial grid with points between 200 and 500 in the simulations. With these settings, and the time steps specified above, these conditions are fulfilled for the domains of parameters studied here.
As mentioned above, a significant fraction of the excitation processes of He atoms leads to the population of the 2 1 P state, from which resonant VUV photons are emitted. Resonance radiation from this state largely dominates the radiation from higher, n 1 P (n > 2) states [61], therefore radiation from higher P states is not taken into account in the present model. Here, photons are treated as discrete (classical) particles and are treated using a Monte Carlo approach [61] (also described in details in [52]). They are followed from their creation, through their propagation, repeated absorption and re-emission events, until they either reach the boundary of the region of interest or get lost in photoionisation processes. As already mentioned, the propagation of the VUV photons is followed in the real 3D geometry of the plasma source (see figure 1). Whenever an excitation event to the He 2 1 P state at a position x 0 and at time t 0 occurs in the 1D PIC/MCC simulation of the plasma, in the photon tracing module a VUV photon is scheduled to appear at a 3D position (x 0 , y 0 , z 0 ) and at a (later) time t = t 0 − τ 0 ln(1 − R 01 ), where y 0 = R 01 L and z 0 = R 01 D (see figure 1), and τ 0 is the lifetime of the 2 1 P state.
The central wavelength of the 2 1 P → 1 1 S transition is λ 0 = 58.4334 nm. Unperturbed, stationary atoms emit photons near this wavelength, as defined by the natural broadening of the spectral line caused by the finite lifetime of the excited state. In a plasma environment the radiation of excited atoms may be perturbed by other particles, which is accounted for, e.g., by pressure broadening. Moreover, the thermal motion of the emitters alters the radiation wavelength due to the Doppler effect. The natural and pressure broadening cause a Lorentz profile, whereas the Doppler broadening results in a Gaussian profile. These mechanisms combine to yield a Voigt profile [62].
Considering a single photon, the emission wavelength is sampled randomly [61] from the Lorentz-broadened lineshape as: where Δλ L is the total homogeneous (natural and pressure) broadening [52]. In the algorithm, this wavelength is altered by the Doppler shift; first, a random velocity vector, v A , for the (radiating) atom is sampled from the Maxwell-Boltzmann distribution of the background atoms. (A stationary background gas is assumed here as well like in the case of ion-atom/molecule collisions.) Second, a random direction, represented by a unit vectorê ph , is assigned to the radiated photon. Finally, the photon is emitted with a wavelength [61]: The probability of the absorption of the photons is computed using where Δt ph is the time step for the photon tracing, c is the speed of the light, n 0 is the background gas density, and σ phabs is the photoabsorption cross section: where g 1 and g 2 are the statistical weights of the respective levels, A 21 is the Einstein coefficient for spontaneous emission, and V is the Voigt-profile (computed using the method given in [63]). The cross section (7) and the probability of absorption during a given time step are computed by taking a randomly chosen (potentially absorbing) atom from the background gas with velocity v B , and using the wavelength that takes into account the Doppler shift caused by the motion of the potentially absorbing atom. The time step used here is Δt ph = 4.5 × 10 −18 s. Following each time step, the probability given by equation (6) is compared to a random number R 01 , and an absorption event takes place if R 01 < P c,vuv is fulfilled. Upon absorption, a random lifetime to the excited state is assigned. Emission events, of which the details were given above, are initiated when this lifetime is passed. At the electrode surfaces, secondary electron emission and elastic electron reflection are taken into account [29], with the following coefficients: γ(He + ) = 0.3, γ(He + 2 ) = 0.2, γ(N + 2 ) = 0.1, γ(N + 4 ) = 0.1, and γ(VUV) = 0.1, as well as r = 0.6. These values (except γ(VUV)) were used in previous simulations of the same plasma source [28,29].
The discharge is excited by a RF voltage, with either of the two excitation waveforms: • A simple cosine waveform, V(t) = V 1 cos(ωt), where ω = 2π f 1 and f 1 = 13.56 MHz; • A 'TVW', V(t) = 4 k=1 V k cos(kωt + θ k ), which consists of four consecutive harmonics of the base frequency f 1 . V k and θ k are, respectively, the amplitudes and phases of the individual harmonics. With θ 1 = θ 3 = 0 and θ 2 = θ 4 = π, the above form generates a valleys-type waveform that has long, relatively flat positive parts and a short, high magnitude negative part (see e.g. figure 2 of reference [64]).
The amplitudes of the harmonics appearing in both waveforms are set according to a specified peak-to-peak voltage, V pp , as φ k = 2V pp (N − k + 1)/(N + 1) 2 with N = 1 for the single-frequency case and N = 4 for the multi-frequency case. These cases will also be referred to '1F' and '4F' cases, respectively.
The convergence of the simulations (which is monitored by the stabilization of the number of superparticles, except of their statistical fluctuations) takes several 100 RF cycles. Data for analysis are collected over further 100 RF cycles using 10 copies of the converged state system (i.e., data are collected from 1000 RF cycles). Due to the high computational demand of the code the number of superparticles of the dominant ionic species (being N + 4 ) is kept at 5 × 10 4 ± 10%. All 'active' particles in the simulations have the same superparticle weight; these include all the charged species, the atoms in the resonant and metastable states, as well as the VUV photons. The number of superparticles for a given species can be revealed by comparing the density of this species with that of N + 4 . With the settings specified above, the simulation of a single RF cycle takes about one hour on a fast CPU.

Results and discussion
Below, we illustrate the discharge characteristics obtained from the computational model described in the previous section. The results are presented for (i) N = 1 with voltages V pp = 400 V and 685 V and (ii) N = 4 with voltages V pp = 400 V and 475 V. The lower voltages are close to the minimum maintaining voltages of the plasma, the higher values still belong to the stable operation regime, according to experiments [29]. The effects of some of the individual elementary processes will also be illustrated during the discussion. Figure 2 shows the potential distribution in the discharge for single-frequency (N = 1, (a1)) and multi-frequency (N = 4, (b1)) operation for the same voltage of V pp = 400 V (for an easier comparison), with excitation waveforms specified in section 2.2. The small panels at the bottom of this figure (a2) and (b2) give a graphical representation of the time-dependence of the discharge voltage, φ(t), which, for N = 1 is identical to the excitation voltage, φ(t) ≡ V(t), while for N = 4 it is modified by the DC self-bias voltage (η), as φ(t) = V(t) + η. In this case, η = 28 V is generated by the electrical asymmetry effect [65][66][67]. A notable difference between the φ(x, t) distributions between the N = 1 and N = 4 cases is, that the distribution is highly asymmetric for the multi-frequency case: the sheath forms at the powered electrode during about one fifth of the (fundamental) RF period only [68]. Correspondingly, a sheath at the grounded electrode is found during most of the RF period. At the powered electrode, thus, both sheath formation and sheath collapse are rather fast; this has consequences on the distributions of the excitation and ionization rates, as will be demonstrated later.
Next, we analyse the composition of the plasma and the density distributions of the charged species. As, apart from the electrons, the density of these species does not exhibit significant temporal dynamics, only the temporally averaged distributions are shown in figure 3. The first observation is that N + 4 ions are the dominant charged species in the plasma for all the conditions. For the multi-frequency ('4F') cases (see the right column of figure 3) electrons are present in the central region of the plasma with same density as the N + 4 ions, however, for the single-frequency ('1F') cases, especially for the lower voltage case of V pp = 400 V (panel (a1)) the temporally averaged electron density is lower at all positions as compared to the density of the N + 4 ions. The plasma at these conditions operates in the 'non-neutral regime', see [69].
The second most important ionic species appears to be N + 2 , with at least an order of magnitude lower density as compared to that of N + 4 . It is interesting to note that while the N + 2 density peaks at the center of the plasma at N = 1 and V pp = 400 V, a minimum at or near the center of the gap is found at the other conditions, as a consequence of the N + 2 → N + 4 conversion, which is especially efficient when the motion of N + 2 ions is primarily diffusive (i.e., in regions with low electric field). He + and He + 2 ions represent relatively minor contributions to the positive charge density. Except for the low-voltage 1F case these ions are present mainly near the electrodes. The He + ions created by electron impact rapidly convert to He + 2 ions, therefore the density of the latter dominates among the helium ionic species. Even though the secondary electron emission coefficient of the helium ionic species is higher than that of nitrogen ionic species, N + 4 ions dominate the ejection of electrons at the electrode surfaces. The peculiarities of the density distributions can be understood based on the analysis of the of the rates of elementary processes determining their balance. This discussion, presented next, is however, limited to the cases with the higher peak-to-peak voltages. Figure 4 shows the spatio-temporal distributions of the rates of the electron impact ionization of He, as well as the electron impact ionization, Penning ionization, and photoionization of N 2 . Figure 5 depicts similar data for the ion conversion processes He + → He + 2 and N + 2 → N + 4 , as well as the loss of N + 4 due to recombination with electrons. The data are shown for the N = 1, V pp = 685 V and the N = 4, V pp = 475 V cases, in the left and right columns of these figures, respectively.
In the N = 1 (V pp = 685 V) case, electron impact ionization of He takes place within the sheath regions where electrons emitted from the electrodes and created in Penning ionization can be accelerated to high energies. This 'shape' of the source function of He + , together with fact that these ions have a limited diffusion length due the ion conversion process (3) explains why these ions are only present near the electrodes (as seen in figure 3(a2)). Electron impact ionization of N 2 occurs mainly near the edges of the expanding sheaths and also near the electrodes at the collapse of the adjacent sheaths. The peak of the electron impact ionization rate of N 2 is only somewhat lower than that of He, although N 2 is present with a concentration of 0.1% only. The reason for this is the significantly  lower ionization potential of N 2 as compared to that of He. Excitation of He atoms is much more prominent than their ionization, this results in a significant Penning ionization of N 2 molecules, that actually dominates the overall ionization rate (see figure 4(a3)). The rate of this process does not exhibit any temporal dynamics as the lifetime of He metastable atoms associated with this reaction exceeds the RF period. Photoionization of N 2 molecules, although with a low rate, spreads completely in space and time as the radiation fills the whole volume of the plasma source.
The rate of the conversion of He + to He + 2 (see panel (a1) of figure 5) shows the same patterns as the creation of He + but it  is shifted in time by about 0.15T, corresponding the characteristic lifetime of He + ions, defined by this conversion process. Their short lifetime results in pronounced temporal dynamics. The motion of He + 2 ions is dominated by drift towards the electrodes, only a small fraction of them diffuses towards the centre of the plasma, creating a low density population there (see figure 3(a2)). The N + 2 → N + 4 reaction, similarly to the main creation channel of N + 2 (Penning ionization), shows no temporal dynamics. Its spatial distribution resembles that of the Penning ionization rate. The recombination rate (figure 5(a3)) is appreciable only at times and locations, where electrons are present in high density and preferably with low energy; observe the correlation between S N + 4 ,rec (x, t) in figure 5(a3) and the electric field distribution, E(x, t), in figure 6(a). As the recombination rate is significantly lower than the creation rate of N + 4 , its role is relatively minor at these conditions.
Consequently, their impact into the electrodes represents the major loss channel for N + 4 . The spatio-temporal distributions of the reaction rates for the N = 4 (V pp = 475 V) case, shown in the right columns of figures 4 and 5 are more complex due to the TVW excitation of the discharge. The ionization rate of He is almost an order of magnitude higher as compared to the single-frequency case (compare panels (a1) and (b1) of figure 4). Similarly to the 1F case, intensive ionization occurs within the sheath regions that build up for a short time at the powered electrode and for a long fraction of the RF period at the grounded electrode. Ionization also takes place near the edge of the expanding sheath (at t/T ≈ 0.45, x/H ≈ 0.2 at the powered electrode and at t/T ≈ 0.65, x/H ≈ 0.7 at the grounded electrode), as well at sheath collapse at the grounded electrode (at t/T ≈ 0.45, x/H ≈ 0.85). As the duration of this collapse is short as compared to the RF period, electrons have to be accelerated by an electric field reversal towards the grounded electrode (see figure 6(b)) to ensure the flux compensation of oppositely charged species on time average. The patterns of the ionization rate of N 2 (figure 4(b2)) follow closely those of He, although production of N + 2 ions within the sheath regions is less important, the maxima of the rates are also comparable. Penning ionization (figure 4(b3)) strongly dominates ion production, however, due the dominance of inelastic collisions between electrons and He atoms in the vicinity of the grounded electrode the highest S N 2,Pen occurs in this region as well. Photoionization (figure 4(b4)) becomes spatially non-uniform at N = 4, its rate is, however, much less modulated in space as compared to the rate of Penning ionization.
Following the source of He + ions their conversion to He + 2 is confined to the spatial regions near the electrodes. N + 2 ions are converted to N + 4 with the highest rate in the vicinity of the grounded electrode. As this region belongs to the sheath domain for most of the RF period there is also significant loss of the N + 4 ions towards the electrodes. Recombination loss of N + 4 is pronounced in the bulk region of the plasma where electron density is the highest and electron energy is the lowest. It is interesting to note that the recombination rate decreases at times when the powered sheath is expanded, as during this period of time an electric field appears in the bulk, as seen in figure 6(b), that accelerates electrons to higher energies, creating less favorable conditions for recombination. Figure 7 shows the source function of the He VUV photon generation (the excitation rate of the He 2 1 P level), the density of the He 2 1 P level, as well as the density of metastable He atoms for the conditions studied above. For the N = 1 case, the VUV radiation has a source within the sheath regions and near the edges of the sheaths at expansion and collapse, with values around 10 17 cm −3 s −1 , see figure 7(a1). As explained in section 2.2, electron impact excitation of He atoms is the primary source of the VUV radiation, and the VUV photons propagate in a way that they are several times reabsorbed by the ground state He atoms and are re-emitted. This peculiarity of the transport of the radiation creates He atoms in the 2 1 P in the whole domain with more or less uniform distribution as figure 7(a2) reveals. The density of the He atoms in this excited state is in the order of 10 10 cm −3 . For a comparison, the density of He metastable atoms is shown in figure 7(a3); the metastable density appears to be about an order of magnitude higher as compared to that of the 2 1 P state.
As compared to the 1F case, a significantly higher excitation rate as well as 2 1 P and metastable atom densities are found in the 4F scenario (see the right column of figure 7). The highest production rate is in the order of 10 19 cm −3 s −1 near the grounded electrode, at times of sheath collapse. Production of the 2 1 P state is also appreciable at the same time at the powered side of system, and within the expanded sheaths at both electrodes. The highest density of He atoms in the 2 1 P state is found in the vicinity of the grounded electrode as a consequence of the strong excitation rate in that region. The density of metastable atoms ( figure 7(b3)) is also about an order of magnitude higher in this case, as compared to the 1F scenario.
The importance of processes related to the VUV photons: (i) secondary electrons from the electrodes due to photoemission as well as (ii) their contribution to the ionization in the plasma via the photoionization of the N 2 molecules, can be evaluated by comparing the simulation results presented so far (which include these effects) with result obtained by disregarding these processes. Figure 8 aids this evaluation by showing the density distribution of the N + 4 ions for various conditions. The data reveal that the density distributions change only marginally (by a few percent at most) when the effects of VUV photons are omitted. As an additional test, a set of simulations was also carried out by neglecting all surface processes (i.e., secondary electron emission due to all species as well as the elastic reflection of electrons). The decrease of the density of the dominant ions is between 0% and ≈20%, depending on the conditions. This comparison concludes that (i) surface processes have generally limited effects on the plasma density and (ii) the resonant He VUV radiation has very little effect on the plasma density. This latter conclusion is somewhat surprising considering the effectiveness of this radiation to enhance the number of charged particles. Other reaction pathways-especially the Penning ionization-are, however, much more important in the ionization and thus minimise the effects of the VUV radiation. Table 2 gives the results of our computations for the flux of VUV photons that leave the plasma source via the orifice, for the various discharge conditions. For N = 1, the escaping flux is ≈10 14 cm −2 s −1 for the lower voltage (400 V pp ) case and increases by a factor of about 23 when V pp is increased to 685 V. For N = 4, the photon flux reaches values between 10 15 and 10 16 photons per 1 cm 2 per second. The further propagation of these photons towards any surface in contact with the effluent is largely determined by the conditions, especially by the mixing of the effluent with the ambient gas. Table 2 also gives the VUV energy flux for the cases investigated. The highest value is found for the 4F case with V pp = 475 V, in this case the energy flux (at the wavelength of 58.4 nm) is 17.7 mW cm −2 .
These computed photon fluxes, unfortunately, cannot directly be compared to experimental results for exactly the same plasma source and operating conditions due to the lack of measured values. Therefore, we quote here some experimental data, which were obtained using absolutely calibrated spectrometers, for other plasma sources. While these sources are different, it is interesting to note that the fluxes are rather similar to those obtained here. In [70], e.g., VUV fluxes of 10 15 -10 16 cm −2 s −1 were obtained in in situ measurements of radiation from low-pressure microwave-produced plasma in Ar/O 2 gas mixtures, while VUV fluxes up to ∼10 17 cm −2 s −1 on Ar resonance lines were measured in a 600 W ICP source at pressures of ≈10 mTorr [71]. In an ICP source as well, VUV fluxes of up to ∼10 16 cm −2 s −1 were found at 25-400 W power and at pressures between 1 and 50 mTorr [72]. VUV energy fluxes in the order of a few mW cm −2 , corresponding to  photon fluxes of ∼10 15 cm −2 s −1 were measured at an electrically floating substrate in a low-pressure 13.56 MHz radiofrequency plasma reactor used for polymer surface treatments in [73].

Summary
In this work, a computational PIC/MCC study of the COST jet plasma in a buffer gas consisting of He with a 0.1% admixture of N 2 , including the treatment of the VUV resonant radiation of He has been presented. The variation of the plasma parameters was computed along one spatial dimension, whereas the VUV photons were traced in the 3D geometry of the COST jet to allow the calculation of the flux of these photons leaving the plasma source via the orifice. The discharge conditions covered various driving voltage waveforms, composed of N = 1 and 4 RF harmonics, and peak-to-peak voltages within the domain of the stable operation of the plasma source. As processes induced by the VUV photons the simulations have incorporated the photoemission of electrons from the electrodes and the photoionization of N 2 molecules. For the operating conditions covered in this study, the results indicated that these processes influence only marginally the plasma characteristics, e.g., charged particle densities. This can be explained by the much higher rate of the Penning ionization process between He metastable atoms and N 2 molecules and by the relatively little importance of surface processes, in general, in the sustainment of the plasma.
VUV photons (with a wavelength of 58.4 nm) with a flux in the range of ≈10 14 − 5 × 10 15 cm −2 s −1 have been found to leave the plasma source via the orifice, corresponding to an energy flux between ≈0.38 mW cm −2 and ≈18 mW cm −2 .
The propagation and the effects of the VUV radiation outside the plasma source, as well as its effects on the effluent, or on any closely placed target, which largely depend on the ambient environment, have not been addressed in this work, but may be topics of future research.
The author thanks Jan Benedikt for a useful discussion on the VUV spectroscopy of plasma jets and Julian Schulze for his comments on the manuscript.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.