Helium metastable species generation in atmospheric pressure RF plasma jets driven by tailored voltage waveforms in mixtures of He and N2

Spatially resolved tunable diode-laser absorption measurements of the absolute densities of He-I (23S1) metastables in a micro atmospheric pressure plasma jet operated in He/N2 and driven by ‘peaks’- and ‘valleys’-type tailored voltage waveforms are presented. The measurements are performed at different nitrogen admixture concentrations and peak-to-peak voltages with waveforms that consist of up to four consecutive harmonics of the fundamental frequency of 13.56 MHz. Comparisons of the measured metastable densities with those obtained from particle-in-cell/Monte Carlo collision simulations show a good quantitative agreement. The density of helium metastables is found to be significantly enhanced by increasing the number of consecutive driving harmonics. Their generation can be further optimized by tuning the peak-to-peak voltage amplitude and the concentration of the reactive gas admixture. These findings are understood based on detailed fundamental insights into the spatio-temporal electron dynamics gained from the simulations, which show that voltage waveform tailoring allows to control the electron energy distribution function to optimize the metastable generation. A high degree of correlation between the metastable creation rate and the electron impact excitation rate from the helium ground state into the He-I ((3s)3S1) level is observed for some conditions which may facilitate an estimation of the metastable densities based on phase resolved optical emission spectroscopy measurements of the 706.5 nm He-I line originating from the above level and metastable density values at proper reference conditions.


Introduction
Non-thermal atmospheric pressure microplasmas have received rapidly growing attention over the last decades due to their relevance for important applications such as plasma medicine (wound healing, cancer treatment), surface modifications, sensors, plasma catalysis, and the removal of volatile organic compounds from exhaust gas streams [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. One of the main advantages of such plasmas is their ability to provide high densities of reactive particle species such as atomic oxygen and/or nitrogen species (RONS) at close to ambient gas temperatures. An important fundamental question (that is, however, also highly relevant for applications) is: how can high reactive particle densities be generated and controlled efficiently as well as selectively?
Microscopic atmospheric pressure plasma jets (µAPPJ) are one of the most promising tools developed for the generation of reactive species. Therefore, these plasma sources have been attracting increasing attention, both from academic and applied perspectives [18][19][20][21][22][23][24]. Such plasma sources are typically operated at a single driving frequency using helium or argon as a feed gas with a small admixture of molecular gases (up to a few percent), e.g. nitrogen, oxygen, water vapour, etc. Due to a lack of fundamental understanding of the mech anisms of plasma generation and control, applications of such discharges are typically optimized empirically. Metastable species, e.g. He-I (2 3 S 1 ) and He-I (2 1 S 1 ), are known to play an important role for plasma generation via Penning ionization and their densities have been measured in single-frequency µAPPJs under selected discharge conditions [25][26][27]. However, concepts to control their density to optimize plasma generation and reactive particle species concentrations are strongly limited.
The main reason why Penning ionization, between helium metastables and reactive gas admixtures, is so important for plasma generation is the high rate of the electron impact excitation from the ground state to the helium metastable levels He-I 2 3 S 1 and He-I 2 1 S 1 . These levels have excitation energies of is 19.8 eV and 20.6 eV, respectively, which is significantly lower than the ionization energy of helium (24.6 eV). Thus, compared to helium ions, helium metastables can be produced more efficiently via electron impact excitation. The energy 'stored' in these metastable states, on the other hand, is high enough to ionize the molecular gas admixtures and, therefore, this (Penning) ionization channel is an important pathway for plasma generation.
For µAPPJs driven by a single frequency voltage waveform at 13.56 MHz and operated in e.g. He/N 2 mixtures, two electron power absorption modes are known: (i) the Ω-mode [28,29], where the ionization/excitation maxima occur inside the plasma bulk during the sheath expansion/collapse phases within the RF period; and (ii) the Penning-mode [29], where ionization/excitation occurs inside the expanded sheaths. Recent investigations show that the ionization dynamics of the Penning-mode has two pathways [29]: (i) the direct channel, where electrons created via the Penning reaction inside the sheath are accelerated by the strong sheath electric field towards the bulk, gain enough energy to cause ionization/excitation, and are multiplied collisionally; and (ii) the indirect channel, where ions, e.g. N + 2 , generated by Penning ionization are accelerated by the sheath electric field towards the adjacent electrode, where they induce the emission of secondary electrons from the surface. Such secondary electrons are accelerated by the sheath electric field towards the bulk, are multiplied collisionally, and contribute to the ionization/ excitation maxima.
At constant nitrogen concentration, increasing the driving voltage ampl itude was found to induce a transition from the Ωto the Penning mode in single frequency µAPPJs operated in He/N 2 . Operation in the Penning mode results in a higher production of helium metastable densities due to a stronger acceleration of electrons generated inside the sheaths, which excite the metastable states [29]. An upper limit for the driving voltage amplitude is defined by a transition of the operation regime of the jet into a 'constricted'-mode, where the power consumption suddenly and drastically increases so that the jet can be damaged by heat [30]. By increasing the molecular gas admixture concentration, the collisionality was found to increase so that electrons cannot gain enough energy to excite helium atoms inside the sheaths. Thus, a transition from the Penning-to the Ω-mode was found to be induced by increasing the reactive gas admixture concentration at constant driving voltage amplitude, since a decrease of the metastable density results in an attenuation of Penning ionization [29].
In µAPPJs operated in 'Penning gas mixtures', Penning ionization is an important destruction mechanism for metastables. Depending on the corresponding reaction rates with various molecular gases, this can cause the helium metastable density to be strongly time modulated within the RF period, e.g. in He + O 2 [27], while in He + N 2 it is only slightly time modulated [29].
In order to expand the operational regime of µAPPJs and to provide a simple concept for effective and selective control of reactive as well as excited species, including metastables and RONS, Voltage waveform tailoring (VWT) was proposed recently [31,32]. The idea is adopted from low pressure discharges [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] and is based on driving the discharge by a voltage waveform that consists of a finite Fourier series of N consecutive harmonics of a fundamental frequency, f 0 . By varying the relative phases and amplitudes of the harmonics different types of waveforms can be generated, e.g. 'sawtooth', 'peaks', and 'valleys'. In accordance with experimental and computational results, it was found that, at atmospheric pressures, 'peaks'/'valleys'-waveforms can induce spatial and temporal asymmetries of the electron dynamics, i.e. the electron power absorption was found to be confined into small regions in space (close to one of the electrodes) and in time (within the fundamental RF period). Results extracted from the particle-in-cell simulations with Monte Carlo treatment of collisions (PIC/MCC) clearly showed that the electron energy distribution function (EEDF) can be controlled efficiently in various regions of interest in space and time within the fundamental RF period by adjusting the driving voltage waveform shape [32]. This provides an opportunity to control the electron driven chemistry in such plasmas by customizing the driving voltage waveform. Technically, this can be done using an arbitrary function generator in combination with a broadband amplifier and a multi-frequency impedance matching unit [41,43]. For applications such impedance matchings are typically important to optimize the power dissipation in the plasma by minimizing the reflected power to save costs. For fundamental research, however, such impedance matchings do not necessarily have to be used as long as the shape of the driving voltage waveform at the electrode is controlled.
Until now experimental investigations of helium metastable densities in µAPPJs have been limited to single frequency discharges [25][26][27]. Here, we present measurements of metastable densities in µAPPJs driven by tailored voltage waveforms in He/N 2 by tunable diode-laser absorption spectr oscopy (TDLAS) for the first time and compare these exper imental findings to results of PIC/MCC simulations performed under identical discharge conditions. We address the effects of the number of consecutive driving harmonics, the driving peak-to-peak voltage, and the reactive gas admixture concentration. We observe that VWT allows us to generate and control metastable densities efficiently. This is understood based on an analysis of the spatio-temporal dynamics of energetic electrons, which can be customized in space and time by VWT. A strong correlation between the metastable creation rate and the electron impact excitation rate from the helium ground state into the He-I (3s) 3 S 1 level is found for some conditions. This might allow us to estimate the metastable densities based on phase resolved optical emission spectr oscopy (PROES) measurements and known metastable density values at proper reference conditions. The paper is structured in the following way: in the second section, the experimental set-up, the data evaluation procedure, and the simulation approach are described briefly. In section 3, the results are presented and discussed. This part is divided into three subsections. First, the helium metastable density and the dynamics of energetic electrons as well as of the helium metastable sources are studied for 'peaks' waveforms as a function of the number of consecutive driving harmonics ( N 4) at a fixed peak-to-peak voltage of 500 V and a fixed N 2 admixture of 0.05%. Then, the effect of the peak-to-peak voltage amplitude on the metastable density and the electron dynamics is investigated at N = 4 and N 2 concentrations of 0.05% as well as 0.1%. The third subsection focuses on the effects of the reactive gas admixture concentration at a fixed number of harmonics and at fixed driving voltage. Finally, conclusions are drawn in section 4.

Experimental set-up
A sketch of the experimental set-up including all diagnostics is shown in figure 1. The measurements are performed using a reference micro plasma jet (COST-jet) developed within the European COST (Cooperation in Science and Technology) Action MP1011 on 'Biomedical Application of Atmospheric Pressure Plasma Technology' [30]. The jet consists of two identical and co-planar parallel electrodes made of stainless steel. The electrode gap is 1 mm. The electrodes are enclosed by two quartz plates confining the plasma volume to 1×1×30 mm 3 . More detailed information about the plasma source, including its design, can be found in [30].
2.1.1. Atmospheric pressure plasma jet driven by tailored voltage waveforms. The discharge is operated in helium (5.0 purity) with nitrogen (5.0 purity) admixture. The gas handling system has been described in [29,32]. The flow rate of helium is fixed at 1000 sccm while the nitrogen flow rate is varied from 0.5 to 4 sccm. This flow rate range of N 2 corresponds to reactive gas concentrations between 0.05% and 0.4% in the He buffer gas. The gas pressure inside the jet is estimated to be 1.02(±0.02)×10 5 Pa.
An arbitrary waveform generator (Keysight 33600A) is used to generate the desired tailored voltage waveform. The signal is amplified by a broadband power amplifier (Vectawave VBA250-400) and then applied to the powered electrode of the COST-jet. Since the level of reflected power is irrelevant for the fundamental studies here, no impedance matching is used and, thus, a high reflected power is present. For the conditions studied in this paper, the applied power is measured in a range of 150-230 W from which more than 99% is reflected backwards. In any application, for which this represents a serious issue, the concepts of recently developed multi-frequency impedance matching networks [41][42][43] can be adopted. We apply the following approach to feed the powered electrode with the desired waveform: the voltage waveform at the electrode is measured using a high voltage probe (Tektronix P6015A with a bandwidth of 75 MHz) and an USB oscilloscope (Picoscope 6402C, 250 MHz bandwidth, 5 Gs s −1 ). The frequency/voltage response of the probe limits the number of applied consecutive harmonics to a maximum of N = 4 for f 0 = 13.56 MHz. This signal is analyzed and corrected for waveform distortions using a LabVIEW software. A corrected signal is then generated by the waveform generator. In this way a corrective feedback loop is used and repeatedly executed until the deviation between the measured and desired phases/amplitudes of all harmonics, used to create the waveform, reaches values of 1%-3%.
During the measurements the peak-to-peak value of the applied voltage waveform, φ pp , is kept constant, to facilitate comparison with PIC/MCC simulations, where the driving voltage is an input parameter. In this work, the maximum value of φ pp is 500 V due to limitations of the broadband amplifier. (Above this threshold, a non-linear response and saturation of the power amplifier occurs.) The values of φ pp used in this work are chosen in order to generate a stable discharge for the range of reactive gas admixture concentrations and for the different numbers of consecutive driving harmonics.
'Peaks' and 'valleys' waveforms are generated by applying a waveform: where φ k are the amplitudes of the individual harmonics of the fundamental frequency and f 0 = 13.56 MHz. The phase angles, θ k , are set to 0 • for all harmonics for 'peaks' waveforms, while the phases of the even harmonics are set to 180° for 'valleys' waveforms. The voltage amplitudes of the harmonics are calculated according to: [33]. Examples of 'peaks'-waveforms calculated for different numbers of driving harmonics are shown in figure 2.
During the measurements the impurity level and the stability of the system are monitored by measuring the time integrated emission spectra from the plasma using a USB grating  spectrometer (Thorlabs CC200/M, having a spectral range of 200-1000 nm and a resolution of 2 nm).

Tunable diode laser absorption spectroscopy.
In order to measure helium metastable densities inside the jet, tunable diode laser absorption spectroscopy (TDLAS) is used. The system is configured to measure the absorption profile of the He-I (2 3 S 1 → 2 3 P J , where J = 1,2) triplet transitions with central wavelengths of λ J=1 = 1083.025 nm and λ J=2 = 1083.034 nm, respectively. The measurements are based on a tunable diode laser source (Toptica, DL100 DFB). The line width of the laser (<10 MHz) is much smaller than the width of the absorption line, which is several GHz due to the strong collisional broadening at atmospheric pressure. The system can scan over two triplet transitions within a range of 14 GHz without any mode-hopping. The repetition rate of the scans is set to 3 Hz.
As shown in figure 1, the laser beam passes through the first 'Orifice 1', where the size of the beam is reduced to 200 µm. Then, a 20 mm thick quartz plate is used to split the laser beam into three parts. The transmitted, highest intensity part of the light passes through 'Orifice 2', where the diameter of the beam is further reduced to 100 µm. This beam passes through the COST-jet, which is placed on a motorized stage (Standa 8MT173V-10, with a step resolution of 1.25 µm) and is detected by the 'Photodiode 1' (Thorlabs DET10C/M, 900-1700 nm). The stage is equipped with a controller (Standa 8SMC5-USB) allowing to set a precise position of the plasma source during the measurements. We use 20 µm steps in this work. The angles at which the beam enters the jet through the quartz planes are 90° with respect to a vector pointing along the electrode gap and 73° relative to the direction of the gas flow through the jet. This configuration is used to avoid that the µAPPJ acts as a plane parallel optical resonator. The absorption length, l, is 1.05 (±0.01) mm. By measuring the laser intensity profile along the gap, the exact electrode gap and the beam size are obtained to be 0.98 (±0.02) mm and 105 (±5) µm, respectively. The detected power of the laser beam is less than 5 µW. The metastable density is expected to be homogeneous along the direction of the gas flow [50]. In order to avoid any influence of plasma formation at the gas inlet as well as any effects of back diffusion of ambient air at the jet nozzle on the measurements of the metastable densities, all results described in this paper are obtained in the middle of the active plasma volume along the direction of the gas flow.
A second, low-intensity beam generated by the beam splitter (20 mm quartz plate) is guided towards a scanning Fabry-Perot interferometer with a built-in photodetector (Toptica FPI 100 with 1 GHz free spectral range) to detect any mode-hoping and to measure the line-width of the absorbing transitions. The third laser beam generated by the beam splitter passes through a hollow cathode lamp and its intensity is measured by the 'Photodiode 2' (Thorlabs DET10A/M, 200-1100 nm). The lamp is used to provide a reference line profile of the specific absorption profiles used in this work. In this way the absolute wavelength of the laser can be determined. The lamp consists of a cylindrical hollow cathode with an inner diameter of 20 mm and a length of 30 mm, placed inside a vacuum envelope equipped with quartz windows. The chamber is filled with helium gas at a pressure of 50 Pa. The lamp current is set between 10 and 20 mA during the measurements.
An example of the three measured signals, i.e. the absorption line profile in the plasma measured by 'Photodiode 1', the Fabry-Perot signal to determine changes of the laser wavelength, and the absorption signal from the hollow cathode lamp measured by 'Photodiode 2' to determine the absolute value of the laser frequency, is shown in figure 3. The data are recorded with a resolution of 345 points per GHz.

Data processing and evaluation procedure
Absolute metastable densities are obtained based on the Beer-Lambert law by measuring the transmittance, T ν : where I(ν) and I 0 (ν) are the intensities of the transmitted radiation as a function of frequency (ν ) with and without the presence of absorbing particle species generated in the plasma, respectively. In order to remove the background noise and light, irradiated from the plasma, four signals are recorded, namely the intensity of the transmitted radiation with (1) plasma and laser switched on, I Lon Pon , (2) plasma on and laser off, I Loff Pon , (3) plasma off and laser on, I Lon Poff , and (4) plasma off and laser off, I Loff Poff . Each of these signals is saved and averaged for 32-128 scans at every position, x, within the electrode gap at which measurements are performed. The measurements are performed automatically using the computer-controlled 'shutter' and 'power amplifier' depicted in figure 1. The population density of metastables, n * , can be derived from the absorption coefficient, k(ν) [51]: where e is the elementary charge, ε 0 is the vacuum permittivity, c is the speed of light in vacuum, m e is the electron mass, f J is the absorption oscillator strength (f J=1 = 0.18 for the 2 3 S 1 → 2 3 P 1 transition and f J=2 = 0.30 for the 2 3 S 1 → 2 3 P 2 transition [52]). F(ν) is a normalised function ( +∞ 0 F(ν)dν = 1) which represents the absorption line shape. Based on equations (2) and (3) the absolute line-averaged density, n * , based on the area, S J , under the line-absorption curve: In our system, the shape of the absorbed line can be described by a convolution of Lorentzian and Gaussian functions, i.e. by a Voigt profile [53]. The Lorentzian part is mainly caused by the pressure broadening mechanism. The influence of other broadening mechanisms such as natural and Stark broadening are negligibly small under the experimental conditions studied here [54,55]. The width (FWHM) of the Gaussian component is due to Doppler broadening [56]. Here, T J is the temperature of the absorbing species, which is assumed to be equal to the gas temperature, k B is the Boltzmann constant, and M He is the mass of a helium atom. Inside the COST-jet the gas temperature is within the range of 300-350 K [57,58]. This results in ∆ν D ≈ 1.72-1.85 GHz. At atmospheric pressure, the Lorentzian width, ∆ν L , is of the order of 10 GHz [27,59].
The triplet transitions at the wavelengths λ J=1 and λ J=2 are only 8 pm or 2.3 GHz apart from each other. Therefore, in contrast to low pressure plasmas, e.g. the hollow cathode discharge used here (see figure 3 (bottom)), at atmospheric pressure these two absorption line profiles overlap significantly and cannot be separated. Based on the arguments presented above, the double Lorentz profile can be used instead of a double Voigt profile [60] as a fit function for the measured absorption spectrum in order to obtain the area below the measured absorption profile, S J . The error of determining S J in this way is below 3%. The corresponding function can be written as the sum of two Lorentz functions for the triplet transition for J = 1 and J = 2: where ∆ν J and ν J are the width and the central position of the corresponding transitions, respectively. f b0 is an offset and is typically equal to the experimental noise level. The number of fitting parameters can be significantly reduced by taking into (4)): An exemplary spectral profile fitted to measured data based on equation (7) (blue dotted line) is shown in figure 3(a). The mode hopping of the laser limits the frequency range used to scan the absorption curve. The fitting procedure is performed automatically after the absorption signal is measured by using a LabVIEW software. We have not observed any significant difference (more than ±0.1 GHz) in the position ν J=2 determined from the fit and from the hollow cathode absorption signal. The obtained width ∆ν J=2 is found to be constant within experimental uncertainties (±0.2 GHz) and independent of the nitrogen admixture concentration, as well as the peak-to-peak voltages used in this experiment. This indicates that the gas temperature inside the jet is approximately constant for all discharge conditions studied. Under discharge conditions, at which the helium metastable densities are low ( 5 × 10 10 cm −3 ), in order to reduce uncertainties, the fitting parameters ∆ν J=2 and ν J=2 are fixed to the values obtained at higher densities. The reproducibility of the measurements was checked over a period of 2 months and the results did not differ by more than 10%-15%, confirming the stability of the experimental system. Based on the assumptions discussed above and including the uncertainties of the used/fitted parameters, we estimate the accuracy of our measured helium metastable densities to be better than 10% in addition to the background noise level that generates an error of ±1.5 × 10 10 cm −3 .

PIC/MCC simulation
The plasma is described computationally by our electrostatic '1d3v' PIC code that incorporates a Monte Carlo treatment of collision processes (PIC/MCC) [61][62][63]. This code has previously been applied to investigations of atmospheric-pressure plasma jets [29,31] and is similar in many respects to that used in studies of atmospheric-pressure nano second pulsed discharges [64]. Details of the discharge model and its computational implementation have already been discussed in some detail in the above references. Therefore, here we only briefly summarize the main features of the model. The model does not aim to account for the details of the complicated plasma chemistry in the jet plasma of the reactive He-N 2 gas mixture. However, a limited set of species and reactions was proven in our previous studies [29,31] to reproduce the main features of the electron dynamics and related ionization and excitation channels in the plasma. Below, we briefly discuss these processes.
The charged particles, that we trace, are electrons, as well as He + , He + 2 , and N + 2 ions. For the electron impact collisions the cross sections are taken from [65] (for e − + He atom collision processes) and from [66] (for e − + N 2 collision processes). (This latter set is largely based on the Siglo cross section set, which is accessible via the LxCat website [67]). The electron + atom/molecule collisions are assumed to result in isotropic scattering. In accordance with this, we use the elastic momentum transfer cross sections for the elastic e − + He and e − + N 2 collisions. In e − + He collisions 50% of the excitation events is assumed to result in the formation of singlet (2 1 S) or triplet (2 3 S) metastable states (He*) either by direct excitation to these levels or via cascades from higherlying states [32,64].
As our studies are limited to small concentrations of nitrogen in the background gas, the collisions of ions are restricted to He atoms as targets. For He + + He elastic collisions an isotropic channel and a backward scattering channel [68] are taken into account. For the collisions of He + 2 ions and N + 2 ions with He atoms the Langevin cross sections of the respective processes are adopted. The ions are created by direct electron impact ionization: via ion conversion: and through Penning reactions: The treatment of the latter two processes is based on the following Monte Carlo approach: whenever He + and He * particles are 'born', random lifetimes are assigned to them based on the reaction rates [69,70] of the respective processes. The given reaction is executed in the simulation after these lifetimes elapsed.
The simulations use different time steps for the various species, optimized according to their collision frequencies.
The smallest time step is required for the tracing of the electrons, due to their extremely high collisionality. Allowing a collision probability P = 1 − exp(−ν∆t) 0.1 (ν being the collision frequency) results in a ∆t e = 4.5 × 10 −14 s upper bound for the time step of electrons. For the ions, significantly longer time steps due to their lower collision frequencies are allowed. For He + , He + 2 , and N + 2 ions, these are, respectively: ∆t He + = 10 ∆t e and ∆t He + 2 = ∆t N + 2 = 100 ∆t e . The buffer gas temperature and pressure are kept constant, at T g = 300 K and p = 1 bar, respectively. The electron reflection probability at the electrodes is set to α = 0.5, which we consider to be a good approximation under the specific conditions covered here [28,29,71]. The ion-induced secondary electron emission coefficients are set as γ He + = 0.2, γ He + 2 = 0.12, and γ N + 2 = 0.07. These values of the γ coefficients are somewhat lower than those used in our earlier studies of the same type of discharge [29,32]. The changes were implemented to ensure a better match (between experiments and simulations) of the domain of stable operation, during our investigations of this plasma source excited by various applied waveforms, at different frequencies and voltage amplitudes. Electron emission at the electrodes due to metastable atoms is neglected, since only a few of these atoms are expected to reach the surfaces due to the decreased diffusion speed at high pressure [72] and due the high probability of the Penning reaction. We use N g = 800 grid points to resolve the gap between the electrodes.

Effect of the number of consecutive driving harmonics
In this section, we study the effects of the number of consecutive harmonics, N, in the 'peaks' type driving voltage waveform (equation (1)) on the helium metastable density profile in the COST jet. Figure 4 shows the measured ((a)-(d)) spatially resolved and time averaged helium metastable densities for a 0.05% N 2 admixture concentration to helium, for φ pp = 500 V, and N = 1 − 4 at a fundamental frequency of 13.56 MHz, along with the corresponding distributions obtained from the PIC simulations ((e)-(h)). These conditions are chosen in order to generate stable discharges and metastable densities that exceed the detection limit of the experimental system (1-2×10 10 cm −3 ).
Our simulations show that the time modulation of the metastable particle densities is very weak (within a few percent). This results from the facts that the main destruction mechanism of the metastables is Penning ionization with nitrogen molecules and the characteristic time of this reaction is longer than one period of the fundamental driving frequency, for the N 2 densities considered here. Therefore, it is sufficient to measure and analyze the temporally averaged metastable densities.
Helium has two metastable species: the He-I 2 1 S 1 singlet metastables and the He-I 2 3 S 1 triplet metastables. In the plasma the population of triplet metastables prevails over that of the singlet metastables because of an efficient conversion channel of superelastic collisions of the singlet metastables with thermal electrons [72,73]: These spin changing collisions are not included in the simulations. Therefore, in figures 4((e)-(h)) we show the sum of the densities of metastable species obtained from the simulation that we believe can be compared with the measured He-I (2 3 S 1 ) triplet metastable densities. Compared to the experiment, the direct results of the simulations ((e)-(h) in figure 4) show higher metastable densities, e.g. a two times higher peak density for N = 4, and much lower metastable densities at x = 0 for N > 1. This can be explained by the fact that the spatial resolution of the PIC data is 1.25 µm, while in the experiment the step size of the measurements is 20 µm and the laser beam diameter is approximately 105 µm. Thus, in order to compare the experimental to the computational data, an averaging procedure is applied. For each measurement position x, that corresponds to the center of the laser beam profile, the PIC data are averaged over a region of interest from x − 52.5 µm to x + 52.5 µ m around x to match the spatial resolution of the experiment (see figures 4((i)-(l)). The measured time averaged helium metastable density profiles are in excellent agreement with these computed profiles for all N after applying this averaging procedure.
The experimental data in figure 4 clearly show that adding more consecutive harmonics to the driving voltage waveform (i) drastically increases the maximum metastable density, from 0.7 × 10 11 cm −3 at N = 1 to 48 × 10 11 cm −3 at N = 4, i.e. by a factor of about 68, and (ii) breaks the spatial symmetry of the density profiles. This is achieved by VWT at constant peak-to-peak voltage and constant reactive gas admixture concentration. This strong effect of VWT on the metastable density profiles can be explained by the effect of the driving voltage waveform on the EEDF in distinct spatio-temporal regions of the discharge [32].
The calculated production rates of helium triplet and singlet metastables via electron impact excitation from the ground state are shown in figures 4(m)-(p). The energy difference between these two states (≈ 0.8 eV) and the different shapes of the respective cross sections cause the production rate of triplets (with a threshold energy of 19.8 eV) to be approximately three times higher than the production rate of singlets (with a threshold energy of 20.6 eV). This ratio changes, however, somewhat with N. As the Penning reaction, that 'destroys' metastables, takes place on a shorter time scale as compared to the time scale of the diffusion of He metastables at atmospheric pressure, we do not consider the latter in the simulations, and consequently the source functions of the metastables are highly correlated with their densities.
As He metastables are generated by energetic electrons, their source (and density) profiles are expected to be correlated with the distribution of the electron impact excitation rate from the helium ground state into other high-lying excited He levels. One of such levels is He-I (3s) 3 S 1 based on which the 706.5 nm emission line (with a threshold excitation energy of 22.7 eV for electron impact excitation from the ground state) is often chosen in experiments to characterise the plasma by phase resolved optical emission spectroscopy (PROES) [29,32]. This method allows determination of the electron impact excitation rate from the ground state into the upper level of the selected transition [29,32,74]. In the following we explore the validity of the expected correlation between the sources of metastables and this excitation rate of the He-I ((3s) 3 S 1 ) level. To facilitate this, we included the computation of the electron impact excitation rate from the helium ground state into the He-I (3s) 3

S 1 level in the simulations. Figures 4(r)-(u)
show the time averaged profiles of this excitation rate. A comparison with the corresponding metastable source functions (figures 4(m)-(p)), in terms of relative magnitudes and spatial distributions, shows indeed a high degree of correlation. The ratio of the maxima of the time averaged excitation rate and the metastable density peaks near the grounded and the powered electrodes are nearly equal for a given set of discharge conditions. This agreement implies the important conclusion that spatially-resolved emission spectroscopy of the He-I 706.5 nm line can be used to estimate the spatial profile of the metastable densities without performing complex TDLAS measurements.
The similarity of the source (excitation rate) of metastables and the excitation rate of the He-I ((3s) 3 S 1 ) level is further explored in figure 5 that shows spatially and temporally resolved distributions of these quantities. According to figure 4(a)), the single frequency scenario (N = 1) is characterized by a spatially symmetric pattern of the metastable densities. This is caused by the symmetric spatio-temporal electron impact excitation dynamics shown in figure 5(a): two excitation maxima with the same intensity are observed during sheath expansion at each electrode at equal distance from the adjacent electrode within each RF period. Similarly, the time averaged He-I (3s) 3 S 1 excitation rate profile (figure 4(r)) as well as the temporally and spatially resolved rate (figure 5(d)) show two identical peaks at similar distances from the adjacent electrode located at the positions of the metastable maxima.
Under these conditions the discharge operates in the Ω-mode, where ionization/excitation occurs mainly in the plasma bulk during sheath expansion and collapse due to the presence of strong bulk electric fields at the times of maximum RF current. The electric field in the bulk is generated to drive the current through the bulk and to ensure current continuity in the presence of many electron-neutral collisions. The field is stronger at the sheath edges compared to the discharge center due to a decrease of the plasma density from the center towards the electrodes.
Inside the sheaths, striations in the excitation patterns near the electrodes (see figure 5(d)) can also be observed [29]. This effect is caused by electrons which after being generated at the electrodes (via secondary electron emission due to ion bombardment), are accelerated towards the bulk by the sheath electric field and cause excitation repeatedly. Based on this electrons gain and lose energy repeatedly and, thus, striations are generated, similar to the Franck-Hertz experiment  [ [75][76][77][78][79]. As mentioned above, we assume that 50% of all excitation events for e − + He atom collision processes form helium metastables. Therefore, the excitation process from the ground state into He-I (3s) 3 S 1 level is not the only process that leads to the generation of helium metastables. Since energy thresholds of the excited levels differ by several eV, the periodic structures, present in the spatio-temporal plot of the excitation rate for each level, have slightly different positions in space and, thus, the striations are not visible in the spatiotemporal plots of the metastable source (see figure 5(a)).
Adding more consecutive harmonics to the 'peaks' driving voltage waveform allows us to customize the dynamics of the boundary sheaths within the fundamental RF period and to break the symmetry of the spatio-temporal excitation dynamics in space and time [31,32]. According to figures 4 and 5 this asymmetry is increased by increasing the number of driving harmonics. This is caused by the fact that adding more harmonics induces a short and fast sheath collapse at the powered electrode, while the sheath is collapsed for a long fraction of the fundamental RF period at the grounded electrode (e.g. figure 5(f)). During the short sheath collapse at the powered electrode a strong electric field is generated to drive a high electron current to this electrode in order to compensate the continuous positive ion flux to this electrode within one RF period [80]. This field accelerates electrons to energies high enough to cause strong excitation/ionization and metastable generation in the vicinity of the powered electrode, while this does not happen at the opposite electrode, i.e. the discharge symmetry is broken by VWT. At higher N, the sheath collapse time is shorter and, thus, a stronger reversed electric field will be generated during the sheath collapse at the powered electrode. This leads to enhanced excitation/ionization and generation of metastables close to the powered electrode, e.g. for N = 4 the peak metastable density is found to increase by almost two orders of magnitude compared to the single frequency case.
As the COST jet is geometrically symmetric, the role of the electrodes can be reversed by changing the driving voltage waveform from 'peaks' to 'valleys' [φ(t) valleys = −φ(t) peaks ], i.e. high metastable densities adjacent to the grounded electrode can be generated if preferred. This symmetry is tested in figure 6(a) that shows the measured spatially resolved absolute densities of helium metastables for 'peaks' and 'valleys' waveforms with N = 4, φ pp = 500 V, and an N 2 concentration of 0.05%. The measured profiles are mirror images of each other with respect to the middle of the gap, x = 0.5 mm. This also confirms that it is sufficient to consider only 'peaks' waveforms in this work, as density profiles for 'valleys' waveforms can be obtained by mirroring the distributions for the 'peaks' waveforms. Figure 6(b) compares the space-and time-averaged helium metastable densities obtained experimentally and from the PIC simulations as a function of the number of harmonics, N, for N 2 concentrations of 0.05% and 1%, and a fixed peakto-peak value of the driving voltage waveform of 500 V. The agreement between the experimental and simulation results is excellent over a wide range of conditions. Figure 7(a) shows the space-and time-averaged metastable density, n * , as a function of N, normalized by its value at N = 1. Additionally, the space-and time-averaged electron impact excitation rate from the ground into the He-I (3s) 3 S 1 state (normalized in the same way as the metastable density) is also shown as a function of N. As we have already discussed, the latter is directly related to the intensity of the plasma emission of the He-I 706.5 nm line [74]. At the low nitrogen admixture concentrations considered here, we find that changing this concentration has a very little effect on the normalised curves. Moreover, when calculating the ratio of the (spaceand time-averaged) metastable density and the excitation rate we arrive at a nearly linear relationship, independently of the nitrogen admixture concentration, as shown in figure 7(b). This finding implies that the increase of the space-and timeaveraged helium metastable density, n * , as a function of N can be estimated by measuring the He-I 706.5 nm line intensity for a given number of consecutive driving harmonics, I N , and by normalizing this intensity by its value at single frequency operation (N = 1), I N=1 , and otherwise identical conditions.
Based on the linear fit to the data (black solid line in figure 7(b)) and the fact that the excitation rate is directly proportional to the He-I 706.5 nm line intensity [29], for such discharges, the following relation is found to be valid within the parameter range studied in this paper: where n * N is the space-and time-averaged metastable density in a discharge driven at N consecutive harmonics and I N is the measured space-and time-averaged intensity of the He-I 706.5 nm line in the presence of N consecutive driving harmonics. n * N=1 and I N=1 are the space-and time-averaged helium metastable density and the space-and time-averaged intensity of the He-I 706.5 nm line at N = 1. In conclusion, this relation allows to determine the space-and time-averaged helium metastable density at any number of consecutive driving harmonics based on measuring the relative emission intensity of the He-I 706.5 nm line and the space-and timeaveraged metastable density in the single frequency reference case. If these data are known, n * N can be determined without doing TDLAS measurements, but based on emission spectroscopy. This is expected to be very useful, since TDLAS measurements are much more complex and difficult compared to emission spectroscopy. The metastable density in the single frequency reference case can be obtained from simulations and/or from reference experiments based on TDLAS.

Effect of the peak-to-peak value of the driving voltage waveform
In this section, we investigate the effect of changing the peakto-peak value of the driving 'peaks'-voltage waveform on the helium metastable density profile in the COST-jet at a constant number of consecutive harmonics (N = 4) of the fundamental driving frequency f 0 = 13.56 MHz. The admixture concentration of N 2 to helium is also kept constant at 0.05% at atmospheric pressure.
The first row of figure 8 shows the measured metastable density profiles for different peak-to-peak values ranging from 350 V to 500 V. This corresponds to the voltage range at which a stable discharge could be generated. The second row shows the corresponding PIC simulation results with a spatial resolution of 1.25 µm, while the third row shows the simulation data spatially averaged to match the resolution of the TDLAS measurements. Excellent agreement between these computational data and the measured metastable density profiles is obtained for all peak-to-peak voltages.
The spatial symmetry of the helium metastable density profile is broken in all cases due to the specific waveform: the metastable density is significantly higher at the powered compared to the grounded electrode. Similar to the findings presented in the previous section, this is caused by the effect of the tailored driving voltage waveform on the spatio-temporal dynamics of energetic electrons traced by the electron impact excitation dynamics shown in figures 9(a)-(c). The highest metastable creation rate (row 4 of figure 8) and the highest electron impact excitation rate from the ground into the He-I (3s) 3 S 1 level (row 5 of figure 8) always occur at the powered electrode during the local sheath collapse phase (see figure 9), when the high energy tail of the EEDF is enhanced at the powered electrode due to the presence of a short sheath collapse that causes an electric field reversal. The spatiotemporal distributions of the metastable creation rate and the electron impact excitation rate from the helium ground state into the He-I (3s) 3 S 1 level are highly correlated in all cases as figure 9 confirms, indicating again that a connection between the metastable density and the plasma emission of the He-I 706.5 nm line may be established.
The absolute values of the maxima of the helium metastable density as well as the space-and time-averaged metastable density increase as a function of the peak-to-peak driving voltage. The latter dependence is found to be nearly exponential as shown in figure 10(a) (please note, the vertical axis has a logarithmic scale). This behaviour is caused by an increase of the plasma density as a function of φ pp due to an enhancement of the electron power absorption at higher driving voltage. An increase of the plasma density leads to a higher electron flux towards the powered electrode, which is required to compensate the ion flux to this electrode on time average. Thus, more electrons are accelerated towards the powered electrode by the reversed electric field during the local sheath collapse and gain energies that are high enough to generate helium metastables. Figure 10(b) shows the ratio of the space-and time-averaged electron impact excitation rate from the ground state into the He-I (3s) 3 S 1 level and the space-and time-averaged helium metastable density as a function of φ pp , normalized by the value of this ratio at φ pp = 500 V. Within the parameter range studied in this paper, the data can well be approximated by a linear fit of the form: where n * φ and I φ is the space-and time-averaged helium metastable density and the emission intensity of the 706.5 nm He line at a distinct peak-to-peak value of the driving voltage waveform, respectively. n 500V and I 500V are the corresponding quantities at a reference peak-to-peak value of the driving voltage waveform of 500 V. This connection between the metastable density and the plasma emission intensity of the 706.5 nm He-I line implies that the latter can be used for an approximation of the metastable density provided that this density is known at a set of reference conditions. Figure 11 presents our experimental and computational results obtained as a function of the nitrogen admixture concentration, at fixed values of the base frequency, f 0 = 13.56 MHz, the number of consecutive harmonics, N = 4, and the peak-topeak value of the driving voltage waveform, φ pp = 500 V. The first row shows the measured metastable density profiles, the second and third rows, respectively, show the corresponding PIC simulation results with a spatial resolution of 1.25 µm and with a decreased resolution, matching the resolution of the experimental data. Similar to the previous sections, we find a good agreement between the experimental and the 'lowresolution' simulation results. The time averaged metastable source functions and the electron impact excitation rate from the helium ground state into the He-I (3s) 3 S 1 level (displayed in rows 4 and 5, respectively) are strongly asymmetric, similarly to the density profiles, due to the effects discussed in section 3.2. We observe a strong depletion of the metastable source with increasing N 2 admixture concentration as an effect of the reduction of the high-energy tail of the EEDF in the presence of lower-threshold-energy excitation channels of the N 2 molecules [32]. This decreased source, together with an increased loss rate due to the Penning process at higher [N 2 ] values, results in a significantly lower metastable density at higher admixture concentrations of N 2 .

Effect of the N 2 admixture concentration
Similarly to our studies of the effects of the number of driving harmonics and of the peak-to-peak value of the driving voltage waveform, we observe here as well in figure 11 a high degree of similarity between the source of the metastable atoms and the excitation rate from the ground state into the He-I (3s) 3 S 1 level. The spatio-temporal behavior of these rates also features very similar patterns as revealed in figure 12. For all cases, the maxima of the rates are observed at the collapsing sheath edge due to the field reversal generated at this location (in space and time). A notable feature is the increase of the sheath width with increasing N 2 admixture concentration due to a depletion of the plasma density as a function of the reactive gas admixture. The depletion of the density of He metastables is quantified in figure 13, at the same conditions as discussed above. The computational and the experimental results show, in a very good quantitative agreement with each other, that the metastable density depends very sensitively on the nitrogen admixture concentration. The data suggest a ∼ [N 2 ] −1 relationship based on a fit to the exper imental data (black solid line in figure 13).

Conclusions
We studied the effects of VWT (by applying 'peaks'-and 'valleys'-type voltage waveforms) on capacitively coupled µ-APPJs based on a synergistic combination of experiments and PIC/MCC simulations. In the experiments we have measured spatially resolved absolute densities of He-I 2 3 S 1 metastables in the plasma by tunable diode-laser absorption spectroscopy. Our measurements have been performed at different nitrogen admixture concentrations and peakto-peak voltages with waveforms that consist of N = 1 to N = 4 consecutive harmonics of the fundamental frequency f 0 = 13.56 MHz. The measured metastable densities and the results obtained from PIC/MCC simulations showed a good quantitative agreement.
The absolute density of helium metastables was found to be significantly enhanced by increasing the number of consecutive driving harmonics as well as by increasing the peak-to-peak value of the driving voltage waveform. The simulations have assisted understanding these observed effects by providing information about the spatio-temporal electron dynamics in the form of spatially and temporally resolved electron impact excitation rates from the ground state into the He-I (3s) 3 S 1 level and helium metastable sources. These data showed that the excitation patterns change from symmetrical (in space and time), i.e. from having maxima at the expanding sheath edge, as well as in the bulk plasma (at N = 1), to highly asymmetric profiles (at high N), where the excitation/ionization maxima are concentrated near the vicinity of the edge of the collapsing sheath at only one electrode, where electric field reversals build up due to the short sheath collapse times caused by the highly asymmetric voltage waveforms applied to the discharge. In this way such tailored voltage waveforms lead to an enhanced high energy tail of the EEDF adjacent to one of the electrodes and cause a strong increase of the local metastable density. Increasing the number of consecutive harmonics used to construct such 'peaks'-and 'valleys'waveforms leads to shorter and steeper sheath collapses at one of the electrodes and, thus, enhances the generation of helium metastables. Increasing the driving peak-to-peak voltage leads to an increase of the plasma density and, thus, to a higher metastable density as well.
A strong depletion of the metastable density was, on the other hand, found by increasing the concentration of the molecular nitrogen gas. This effect is a consequence of a decreased generation of metastable atoms due to a decrease of the high-energy tail of the EEDF (as a consequence of the low-threshold-energy excitation channels in N 2 ) as well as of the increased rate of the Penning reaction destroying He metastables more efficiently at higher [N 2 ].
A high degree of correlation between the metastable creation rate and the electron impact excitation rate from the ground state into the He-I (3s) 3 S 1 level was observed for some conditions. This may facilitate an estimation of the metastable densities based on optical emission spectroscopy measurements of the 706.5 nm He-I line originating from the above level as well as metastable density values at proper reference conditions rather than more complex TDLAS measurements.
The observed effects of tailored voltage waveforms on the helium metastable density in µ-APPJs allow us to realize very high metastable densities and to control these densities conveniently by tuning the driving voltage waveform. As such excited particles play an important role for many applications, e.g. plasma medicine, mass spectrometers, etc, and since existing plasma sources can be upgraded to VWT easily by modifying only the external RF power supply, we expect our results to be highly relevant for technological applications in addition to their fundamental importance.