E Faudot

A Cherukulappurath Mana

F Brochard

S Heuraux

Effect of the electrode/wall area ratio on the Plasma potential in discharge and Tokamak plasmas

The plasma potential is a key parameter in plasma discharge or fusion plasma to control plasma-wall interaction or ExB drift. The magnitude of the plasma potential depends on the overall energy transmitted to the plasma via DC or RF devices such as electrodes or antennas. Knowing the plasma potential from the exciting source is useful to prevent high energetic fluxes to the wall or to improve the plasma confinement or expalining shear velocity. The aim of the present model is to calculate this plasma potential whth respect to a DC or RF source in a magnetized or unmagnetized plasma. This double saturated probe model takes into account the electron saturation current and is able to derive the plasma potential as a function of the electrode/wall area ratio for a DC or RF discharge in a Helium/Argon plasma with or without a magnetic field. The results of the model are compared to Aanesland's model in the case of the unmagnetized capacitive sheath, and PIC simulations. The magnetized model is applied to a plasma column with a perpendicular capacitive current in an RF discharge. It appears that the plasma potential can increase to almost the RF potential value at a low wall/electrode area ratio (A wall /A el lower than 5) while the same potential is collapsing as soon as the area ratio (perpendicular over the parallel current area) is higher than the electron/ion saturation current ratio. This is directly due to the saturation of the electron current, preventing the plasma potential from following the imposed RF potential by the electrode, so that the maximum value can be as lower as the floating potential. The perpendicular current involved is mainly a conduction current modeled as a resistive collisional current. In fusion plasma, the maximum plasma potential can rise to higher values than in plasma discharges but the collapse of the potential still occurs for long collisional biased flux tubes.

I. INTRODUCTION

Most experimental plasma devices are using a local electric field strong enough to ionize a low-pressure gas. This electric field can be obtained either with DC (Direct Current) or RF (Radio-frequency) biased electrodes 1 (plasma discharges, thrusters) or RF antennas 2 (fusion plasma, tokamaks, stellarators). In case of RF electrode, this one can be capacitively coupled or directly coupled leading to very different self-biasing of the plasma and electrode 3 . This last coupling is involved in the present work. RF electrodes or antenna induce electrostatic and electromagnetic waves able to bias the quasi-neutral plasma resulting in plasma potential variations, which is a key parameter to understand the behavior of the plasma in either magnetized or unmagnetized devices. For example, a 3D (3 dimensions) map of the plasma potential 4 in a magnetized device allows to understand the ExB drift 5 and the heat fluxes driving hot spot phenomena on ICRH (Ion Cyclotron Resonant Heating) antenna structures 6,7 , and also arcing phenomena 8,9 . Another example is the radial potential profile in the SOL (scrape-off layer) of tokamaks which can trigger better confinement regimes both in linear devices 10 and in tokamaks (H-mode) 11,12 .

The plasma potential magnitude depends strongly on the electrode-plasma interaction and plasma-wall (wall=ground) interaction as well. In other words, it depends on the sheath potentials between the electrode, a) Electronic mail: eric.faudot@univ-lorraine.fr DOI: 10.1109/TPS.2022.3156208 the plasma, and the wall. The sheath physics 13 has to be extended to RF sheath physics 1,14,15 taking into account capacitive currents 16 in unmagnetized or magnetized plasma devices. These currents can change the sheath potential according to their I-V characteristic 17 . As a consequence, the plasma potential can be deduced from the balance of the sheath potential and currents all over the internal geometry of the device and in-vessel pieces. In the most simple model, a 1D particle flux is set to deduce the I-V characteristic of the sheath as it has been widely explained in the Mott-Smith-Langmuir paper 18 . This I-V characteristic can be expressed in planar, cylindrical, and spherical geometries, revealing the strong dependence of the boundaries at the plasma edges. In the special case of planar geometry, the I-V curve is characterized by the classical exponential function used in most modelings, reaching a plateau starting when the potential (of the electrode or a Langmuir probe) is higher than the plasma potential. This plateau corresponds to the saturation of the electron current, and also to the reversal of the sheath potential because electrons have to be accelerated in the sheath 19 . These negative sheaths or electron sheath, first suggested by Langmuir 20 , in which electron density is higher than ion density prevent the plasma potential to rise up to hundreds of volts, even if the RF magnitude of the source supplies thousands volts. Then the understanding of the saturation mechanism leading to the mitigation of the plasma potential is very important to prevent high energetic plasma fuxes accelerated in the sheaths and flowing along ExB drifts. According to the area ratio of the wall to the electrode, this saturation will occur or not. This effect has been measured experimentally 21 but no model was able to explain the evolution of the plasma potential with respect to area ratio and RF exciting potential whatever the sheath main population (electronic or ionic). For exemple in the edge of tokamak plasma, most of sheath can be electronically saturated close to the RF antennas for long enough flux tube connected to the antenna. This saturation is in this case strongly dependant on the magnetic connection length and can be rather complex in the edge. In plasma discharges, this saturation can explain why it is very difficult to increase the plasma potential at several hundreds of volts, despite the powerfulness of the source.

This saturation of the plasma potential usually occurs at high electrode/wall area ratios and has been investigated starting from Lieberman's model [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF] based on the fraction of electrons lost during one RF period 15,19 and for high non ambipolar flow [START_REF] Balruud | [END_REF] . In most modeling 21 the plasma potential is deduced from already saturated sheath (electronic or ionic), so that the complete oscillation of the sheath during one RF period can not be calculated. But this saturation can also be calculated more self-consistently from a "double saturated probe" model (DSP) based upon a continuous sheath characteristics, using the "tanh" function. This is the big novelRF cas ty which allows to calculate the real time plasma potential during one RF oscillation, and to extend the model to a magnetized plasma column, which is also a new possibility of the model. On top of that, the effect of different currents, conduction and/or displacement currents can be investigated thanks to this new modeling. The main goal is to calculate the plasma potential in front of a direct coupling RF biased electrode, magnetized or not, which would be really useful for plasma discharge community, but also for plasma/wall interaction in the Scrape off layer of tokamaks, especially in the vicinity of ICRH antennas.

The nonlinear resolution of this model is presented in the first part in DC and RF cases with only conduction currents in the sheath. In a second part, the model is extended to capacitive RF sheaths. This kind of modeling is close to Aanesland's one. In the final 3rd part, the model is adapted to a magnetized plasma column allowing perpendicular RF currents possibly associated with collisions. Such models have been performed in the past 16 but without taking into account the saturation of electron currents. Here we complete the study of these plasma potentials in the RF environment. In addition, the 3 models are compared to PIC (Particle-In-Cell) simulations which have been performed on a GPU ensuring a very high number of super-particles (up to 10 8 ).

II. SATURATED PLASMA POTENTIAL

The classical I-V characteristics for a 1D planar probe are described in the Mott-Smith-Langmuir paper 18 . This is the simplest way to model the electron saturation is to consider a 1D planar probe, receiving only a longitudinal plasma flux. In reality this is a 3D flux depending on the sheath expansion and then the saturation is not a straight horizontal line but can be given by an OML theory which is quite complicated to introduce in our time dependant model. Replacing this I-V characteristic by a saturated continuous one as plotted in figure 1 using a tanh function, allows solving the current conservation all over the electron and ion saturation regions. This I-V characteristic is given by

I(V ) = I i - I e 2 1 + tanh - eV k B T e + 1 (1)
with I i = j i S the ion saturation current, j i the ion current density, I e = j e S the electron saturation current, j e the electron current density, S the collected current area, V is the potential gap within the sheath, k B the Boltzmann constant, e the elementary charge and T e the electron temperature. This double saturation characteristic implies that no sheath extension would increase the ion and electron saturation currents with respect to the sheath potential (according to Child Langmuir law or OML Theory 18 ). This crude assumption of the model may yield some brutal transition that would not be observed in measurements, nevertheless, the model can take into account the electron saturation current continuously all over the potential domain as plotted in figure 1.

Figure 2 shows an example of a plasma square box with 2 electrodes connected to the same RF or DC source. This configuration allows to decrease the minimum area ratio S 2 /S 1 (equal to 1), with S 1 and S 2 the electrode and wall areas respectively, and has been chosen for the PIC simulations performed afterward to correlate the analytical results. In a 2D configuration, the area ratio is a length ratio L 1 /(2L -L 1 ). This kind of asymmetric plasma discharge can be modeled as a simple electric circuit (figure 3) in which sheaths are represented as diodes, because of their nonlinear I-V characteristics. Due to asymmetric areas of the electrode and grounded wall, each sheath does not drive the same amount of current which is proportional to their respective areas.

A. "Double saturated probe" (DSP) model

Applying the current conservation to this asymmetric double saturated probe model, with φ = eVp kB Te with V p the plasma potential in volts, and the tanh function for electron current yields :

S 1 j i - j e 2 (1 + tanh(1 -φ + φ bias )) + (2) S 2 j i - j e 2 (1 + tanh(1 -φ)) = 0
where φ bias is the potential applied to the electrode, it can be equal to φ dc = eV dc kB Te or φ rf = eV rf kB Te . The ion saturation current is supposed constant all over the potential domain. This is certainly a crude assumption but it simplifies the calculations and does not induce any significant error as long as j e ≫ j i . This model takes into account the electron saturation region so that the sheath can reverse as soon as the applied potential is strongly positive. Applying a DC potential to the electrode, the plasma potential φ can be deduced directly :

φ = 1 + 1 2 φ dc + 1 2 ln S 1 + S 2 -A S 1 + S 2 + A -(3) 1 2 asinh A. cosh (φ dc ) + (S 2 -S 1 ) sinh (φ dc ) (S 1 + S 2 ) 2 -A 2
with A = (S 1 + S 2 )(2j i /j e -1).

On the contrary, to obtain the plasma potential from the RF potential, it must be averaged over one RF period as shown in equation below.

φ = 1 + 1 2 ln S 1 + S 2 -A S 1 + S 2 + A -(4) 1 2 asinh A. cosh(φ rf (t)) + (S 2 -S 1 ) sinh(φ rf (t)) (S 1 + S 2 ) 2 -A 2

B. PIC code description

The PIC code has been written in Matlab using CUDA instructions implemented with Matlab GPU commands. The advantage is a very short code because many advanced commands are already implemented (to search particles in each grid cell for example, or to solve a linear system or to use Fourier transform). The drawback is that the code is slower than a compiled C++ CUDA code. For example, a 2048x2048 matrix multiplication 10% faster with a CUDA code compared to the GPU matlab function. On condition that optimized GPU commands are used, and avoiding loops and preferring vectorized instructions, the productivity gain is real as it could have been using Python language with its CUDA libraries.

Particle motion is achieved using a leap frog method, allowing to "center" the electric field between the initial and final position of the superparticle. The Lorentz force is analytically integrated and written in matrix multiplicators to obtain the exact position even if the timestep is larger than the cyclotron period. The other option would have been to use the well-known Boris algorithm, but it appears the exact method written through matrices is as fast as the Boris method on the GPU. A plasma source injects particles in a square or rectangular box to maintain either the ion density constant or the electron density, depending on the loss rate of each of them. The global plasma temperature tends to decrease with time because the fastest particles are lost the firsts, "cooling" especially the electron distribution function. This effect could be canceled with a heating mechanism but is not yet implemented. 2 electrodes are connected directly (direct coupling) to the RF generator. All other boundaries are connected to the ground according to figure 2. The configuration where a blocking capacitor is put between the RF generator and the electrode (capacitive coupling) [START_REF] Raizer | Radiofrequency Capacitive Discharges[END_REF] is not investigated here because it does not bias the plasma with recpect to the ground contrary to a direct coupling.

Here the goal is to study if the plasma potential can be as high as the RF or DC source. This version of the code includes magnetized plasma motion with an electrostatic Poisson solver only, neither collision nor ionization. The Poisson solver is written using a finite difference Fourier transform method using fft2 functions in 2 dimensions. This method apprears to be the fastest for grid size larger than 256x256 because it can be directly run on the GPU and also because it requires only N 2 log(N ) operations (with N the number of gridpoints) . For smaller grid size, the best method is the decomposition method based on specific LU decomposition for sparse matrix, allowing to use much lower memory than a Gauss method. But this optimized decomposition method can only be used on the CPU. This code has been compared to xoopic [START_REF] Verboncoeur | An object-oriented electromagnetic pic code[END_REF] in unmagnetized and magnetized plasma to be validated. In both cases, the results on density, potential, and time dynamics (plasma frequency) are quite close (less than 1% discrepancy).

On an RTX 6000 GPU, each simulation lasts 5 min for 5000 time steps, which is pretty faster than those executed on a single CPU, between 10 to 50 times faster than on a Core i7 7820HQ. 5000 time steps are enough to reach a quasi-steady state plasma potential. If the duration of the simulation is too long, the cooling effect due to fast particle loss can sensitively alter the result. In PIC simulation fastest particles in the Maxwellian velocity distribution leave the simulation box first while the plasma source is a complete Maxwellian distribution instead of a Maxwellian tail replacing only the fast particles lost [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF] . A way to avoid this cooling is either to inject the same fast particles or heat the plasma with another artificial way.

C. PIC and DSP model comparison

The simulation box is a 5 by 5 cm square as shown in figure 2. The plasma parameters are as follows : n e = 10 16 m -3 , T e = T i = 3 eV, V dc = 100 V. PIC simulations have been performed using 2 [START_REF] Verboncoeur | An object-oriented electromagnetic pic code[END_REF] super-particles in a 500x500 grid steps, which means a space step of 0.1 mm smaller than the Debye length (λ De = 0.13 mm). The time step is 0.1 ns which is lower than 1/ω pe . According to the algorithm, the time step does not need to be lower than 1/ω ce because the cyclotron motion is calculated from the analytic integration of the motion equation with Lorentz force (corresponding to implicit Tajima's method [START_REF] Tajima | Computational Plasma Physics[END_REF] ) and is valid at any time, on condition the Larmor radius remains smaller than a grid step, to see the same electric field during the gyrations (which is the case for electrons). Of course, the ions motion being much slower their characteristic time steps are always valid.

In figure 4, one can see a comparison between the model of equation 3 and PIC simulations for a direct coupling of the electrode facing the plasma.

PIC curves do not fit the theoretical curves based on the "tanh" current conservation within a static area sheath equal to electrode area (S 1 = cst in figure 4). The transition from the unsaturated to the saturated regime, which corresponds to an area ratio S 2 /S 1 higher than the saturation current ratio j e /j i is very sharp in the theoretical model while it is smoother in the PIC simulations. This is due to the plateau of I-V characteristics in which no sheath expansion occurs while this phenomenon is naturally simulated in the PIC simulation. The sheath thickness is close to 0 when the applied voltage at the electrode is equal to the plasma potential, this happens at small area ratios. But as soon as the plasma potential decreases the sheath potential gap (= φ-φ dc ) is reversed and increases (in absolute value) so that the sheath thickness increases as well, resulting in the increase of the electron saturation current passing through the sheath. As a consequence the ratio j e S 1sh /j i S 2 increases, retarding the strong drop of the plasma potential. The static electrode area S 1 has to be replaced by an area depending in the sheath size S 1sh = S 1 + βλ sh . Replacing S 1 by S 1sh in equation 3 and solving the new nonlinear equation numerically, the collapse is much smoother and almost matches the PIC curves (figure 4). For a good fitting β has been set to 1.5 in Argon, and 5 in Helium. The sheath thickness is simply deduced from the Child-Langmuir law

λ sh = αλ De |φ -φ dc | 3 4
, with α ≈ 1, λ De the Debye length and φ -φ dc the potential drop in the sheath. So taking into account the sheath expansion in the analytical model results in a more realistic collapsing transition, more consistent with experimental measurements.

In the first part of the curves, at low area ratios, the slope is much higher in PIC simulations. This can be explained by a small heating of the plasma due to the relatively high applied potential compared to the initial electron temperature. At a small area ratio, the ion and electron temperatures are likely to increase at the beginning of the simulation before reaching lower steady-state values but still higher (by 20%) than the initial loading temperature because the sheath thickness can not be neglected compared to the plasma box size. This heating is not present in the analytical model in which the temperature is relatively low (3 eV) corresponding to a steady state discharge plasma.

In the end, the main mechanism consisting in lowering the plasma potential while the area ratio is increasing is demonstrated. In addition, it is clear that the transition area ratio for which the plasma potential starts to fall down matches the saturation current ratio j e /j i , which simply corresponds to j e S 1 /j i S 2 > 1. If the DSP model provides an analytical expression of the plasma potential for a RF or DC source as a function of the electrode/wall area ratio and can be used as a first calculation to find out the evolution of the plasma potential, the too sharp collapsing region needs the introduction of a self-consistent sheath expansion to be closer to the real experiments.

III. ASYMMETRIC DOUBLE PROBE WITH CAPACITIVE SHEATH IN UNMAGNETIZED PLASMAS

One considers now an asymmetric double probe model with capacitive sheaths [START_REF] Metze | [END_REF] . These sheaths can drive nonlinear conduction currents as shown in the I-V characteristics but also displacement currents due to their capacitance. These capacitances c 1 and c 2 are normalized to k B T e /e because φ = eV /k B T e .

c 1 dφ rf dt -(c 1 + c 2 ) dφ dt = j i (S 1 + S 2 ) -(5) j e 2 (S 1 + S 2 + S 1 tanh(1 -φ + φ rf ) + S 2 tanh(1 -φ))
with

c 1 = ǫ 0 ǫ sh S 1 λ sh .k B T e /e ( 6
)
c 2 = ǫ 0 ǫ sh S 2 λ sh .k B T e /e ( 7
)
λ sh = αλ De eV k B T e 3 4 (8)
with 1 < α < 3 (see next section for more details [START_REF] Chabert | Physics of Radiofrequency Plasmas[END_REF][START_REF] Faudot | [END_REF] ), V = kB Teφ e the time averaged potential drop in the sheath and λ sh the sheath thickness given by the Child-Langmuir law. ǫ sh = 1 -ω 2 pe ω 2 is the relative permittivity of the sheath in the cold plasma approximation. An RF sheath can be considered as electronically empty when the sheath potential is much higher than the electron temperature, as a consequence, the sheath relative permittivity ǫ sh can be considered equal to 1 15 .

Renormalizing the last expression to j i (S 1 + S 2 ) and replacing the sheath capacitances c 1 and c 1 + c 2 by c rf and c respectively :

c rf dφ rf dt -c dφ dt = 1 - (9
)

j e 2j i 1 + S 1 S 1 + S 2 tanh(1 -φ + φ rf ) + S 2 S 1 + S 2 tanh(1 -φ)
No analytical solution can be found for this equation so that it needs to be solved numerically.

A. Comparison to Aanesland's model

Let's remind Aanesland's formula and assumptions to calculate the plasma potential from RF potential V rf and surface ratio S 2 /S 1 between the wall and the active RF electrode respectively 19 .

cos -1 V V rf = π j e j i -exp - eV k B T e S 2 S 1 (10)
with

j e j i = 0.61 2πm e m i ( 11
)
Figure 6 depicts the comparison of Aanesland's and the DSP models for the plasma parameters given in table I. Both models are quite close at lower density (n 0 = 10 15 m -3 instead of 10 16 ) and for Argon gas, which corresponds to a higher capacitive regime for currents within the sheath. As a matter of fact, conduction currents are proportional to the density, while displacement currents are assumed to be inversely proportional to the sheath thickness, depending on the square root of the density. In Aanesland's model, the maximum potential is the RF amplitude V rf (here 100 V) instead of the average over one RF period in the DSP model. This unsaturated (no displacement current) average plasma potential is given by Godyak 14 : V p = V rf /π + V f l , with V f l = T e ln(j e /j i ) the floating potential of the plasma. This result is retrieved for the Helium plot at S 2 /S 1 = 1 with the DSP model. The helium conduction currents are relatively larger than Argon current compared to displacement currents, which explains why there is no increase of the plasma potential before the collapsing transition as it is the case for the Argon case. The higher the capacitive current, the higher the averaged potential which tends to the RF amplitude in saturated regime (high capacitive currents). Aanesland's model can not calculate this first increasing step, while the collapsing transition is close to the DSP model.

B. PIC comparison

DSP model, Aanesland's model, and PIC simulations are compared for the parameters given in table I. Again PIC simulation parameters are the same as for the DC case, 5x5 cm box, 500x500 grid steps, etc... To better fit the PIC curves, the α parameter in the DSP model have been set to 0.3 instead of 0.8, 50/27, 2 5/4 , or even 3 as a function of the author's model or simulations 15,[START_REF] Chabert | Physics of Radiofrequency Plasmas[END_REF][START_REF] Faudot | [END_REF][31] ). The correcting factor for the sheath capacitance has also been calculated by Lieberman 15 as equal to 2.452 for the first RF Harmonic. The question of the RF varying capacitance of the sheath is not straightforward, when the sheath oscillates, its capacitance is much higher at low RF potential and the Child law is only valid at high RF potential (compared to electron temperature in eV). So if we find here a capacitance 3 times higher than expected, it is not far from Lieberman's correction factor. All other parameters have not been changed.

Each plot corresponding to each model appears in figures 7 and 8 for Helium and Argon plasma respectively. The plasma potential V p at the beginning of the plot (S 2 /S 1 = 1) should theoretically start at V rf /π + V f l , which means around 45 V. This value appears to be higher in both PIC and DSP simulations, but the plasma potential continues to increase up to 80 V in Helium and almost 95 V in Argon, before collapsing at area ratios higher than 5 and 10 for both species.

PIC results are quite close to the DSP model except at higher area ratios for which the plasma potential is higher. As for the DC biasing, the sheath thickness changes the real area ratio because the sheath is growing as soon as the plasma potential is lower than the applied RF potential, which also explains the lower decrease rate of the potential in the PIC curve. This effect is even stronger when λ sh is higher than L 1 the electrode width, which means for high area ratios S 2 /S 1 . The same tendencies are observed in the Argon plot. The effect of the sheath extension is less important than in the DC case, because the RF oscillations of the sheath already imply a Child Langmuir thickness whose time averaged is almost constant all over the area domain. This correction is then only visible at higher area ratios. Equation 9 is solved without the correction of the sheath thickness in the electrode area S 1 because the non linear sheath thickness makes the computation much more unstable.

In a near future, this model could be compared to experimental measurements in such asymmetric discharges in Argon 3,19 or Helium, to discuss its experimental validity.

IV. SYMMETRIC DOUBLE PROBE WITH MAGNETIC FIELD

A similar electric model can be built in the case of an RF magnetized discharge. One considers a plasma column connected on both sides to the RF electrode and the grounded wall respectively as shown in figure 9. Each sheath at both sides presents the same area so that the double probe model is symmetric except that perpendicular currents are allowed and must be added to the current conservation equation as follows : Perpendicular current means currents perpendicular to the magnetic field direction. They can be either conduction current (with the conductivity σ) and are modelled by σE = σ(-∇V ) = σ(V 0 -V )/d at the first order. V 0 is the plasma potential outside of the biased plasma column, and V at the center of the plasma column. In this model, V 0 is supposed to be equal to 0 because it is actually much smaller than the applied RF potential, that is why the conduction current can be rewritten as γφ in equation 14. d is the typical scale length of the potential gradient perpendicular to the magnetic field. It can be either a sheath thickness if particles involved in perpendicular current are not well magnetized, or ρ the Larmor radius for magnetized particles. Perpendicular currents can also be polarization currents due to RF oscillations and then are expressed using an equivalent perpendicular capacitance c ⊥ normalized to k B T /e so that c ⊥ = ǫ 0 ǫ ⊥ S ⊥ λ sh .(k B T e /e). Referring to the sketch in figure 10, and ignoring the sheath capacitance c 1 and c 2 as well as the perpendicular inductance M , the equation to be solved can be written as follows : 13)

j + j ⊥ = 0 ( 12
)
-S ⊥ γφ -c ⊥ dφ dt = (
S [2j i -j e /2 (2 + tanh(1 -φ + φ rf ) + tanh(1 -φ))]
Here perpendicular currents are collected over the perpendicular area, corresponding to the lateral boundaries of the biased plasma column, while parallel currents are collected over the parallel area corresponding to 2 times the electrode area. Renormalizing to j i S ⊥ , and rewriting c ⊥ and γ :

-γφ -c ⊥ dφ dt = (14) S S ⊥ 2 - j e 2j i (2 + tanh(1 -φ + φ rf ) + tanh(1 -φ))
The time-averaged perpendicular permittivity used in the perpendicular capacitance can be given by the magnetized cold plasma permittivity tensor, ǫ ⊥ = 1-

ω 2 pi ω 2 -ω 2 ci . If ω ≪ ω ci then ǫ ⊥ ≈ 1 + ω 2 pi ω 2 ci
which is a low RF frequency domain in which capacitive currents occur as presented in equation 14.

If ω ≫ ω pi then ǫ ⊥ ≈ 1 - ω 2 pi ω 2
which is a high RF frequency domain applying in small plasma discharges devices in which magnetic field hardly exceeds 100 mT. In that case ǫ ⊥ is negative, which means the current behaviour is inductive, and then equation 14 should be rewritten as follows, now neglecting the capacitive current :

- d dt (γφ + i L ) = (15) S S ⊥ d dt 2 - j e 2j i (2 + tanh(1 -φ + φ rf ) + tanh(1 -φ)) ( 16
)
with φ = M diL dt and M the normalized perpendicular inductance of the plasma column as shown in figure 10.

ci ν ω ci c ⊥ λ 2 De Csλ sh or ω -1 ci ω -1 ci -γ dφ dt - φ M = ( 17
) S S ⊥ d dt 2 - j e 2j i (2 + tanh(1 -φ + φ rf ) + tanh(1 -φ))
This equation has not successfully been solved, because the time derivative of the exponential term on the RHS is numerically unstable. Nevertheless assuming the behavior of capacitive or inductive current is the same on the plasma potential averaged amplitude, the equation 14 can be solved with both capacitive and inductive parameters presented in table II. The last assumption being very hazardous, PIC simulations have been performed to confirm the similarity of both capacitive and inductive behaviors. But before we first analyze the consequence of capacitive currents alone and resistive currents alone on the time average plasma potential.

A. Capacitive (or inductive) currents only

Setting γ to 0 leads to consider only capacitive currents perpendicular to the magnetic field in the model. In that specific case, one can see in figure 11, that these currents tend to increase the plasma potential up to the maximum value equal to the floating potential plus the RF amplitude 32 . The fact that these currents are phase shifted by π/2 prevents the saturation to occur because the maximum value of the RF potential is reached when capacitive currents are 0 and thus can not reverse the sheath potential. This means that a 90°phased shift current can not make the plasma potential decrease with respect to the area ratio as presented in figure 11. In the present model, only capacitive currents are solved, not inductive currents. Indeed the combination of these opposite phased shift currents can yield in-phase currents, as resistive ones, in the case of resonance for example. This resonance occurs in ICRH heating in which the RF frequency is close to the ion cyclotron frequency. An LC series resonance can also be excited considering the association of the sheath capacitance and plasma column inductance. In that case, RF in-phase currents can reverse the sheath potential and consequently make the plasma potential saturate to a lower value when the area ratio increases. Unfortunately, the model with all kinds of currents stays unsolved as well as the only inductive one. One way to reproduce resonant RF currents effects is simply to consider a strong enough resistive perpendicular current.

B. DC or resistive currents only

Now setting the capacitance to 0, only the DC (or inphase) currents are taken into account. These currents are resistive in case of plasma discharges, collisions slow down the ion motion, making the current to decrease with ν the ion-neutral frequency. On the contrary, in a well-magnetized plasma, such as in the SOL of a tokamak, these collisions increase the ion mobility, and thus the current perpendicular to the magnetic field. This opposite behaviour can be seen in figures 12 and 13 or 14.

These DC currents decrease the plasma potential, by reversing the sheath potential during a fraction of the RF period, ion currents becoming higher than the electron saturation currents collected by the RF electrode. The higher these DC currents, the faster is the collapse of the plasma potential. On the contrary, when these currents are low enough, no saturation occurs and the maximum value is equal to the averaged rectified RF potential, approximately V rf /π.

In the tokamak case, figure 12, a higher collision rate means higher ion mobility and then a sooner collapse of the plasma potential. In such a magnetized plasma the ratio ν/ω ci is always lower than 1 and ions are considered as well magnetized. In typical plasma in the SOL of tokamaks, this transition takes place at 50 < S ⊥ /S < 100 for a collision frequency of 1 MHz. Long flux tubes connected to ICRH antenna (several meters) are expected to be much less biased than the plasma in the near field of the antenna. In this example the RF potential is 100 V, that's why the initial averaged plasma potential (for S ⊥ /S = 1) is equal to V rf /π, approximately 35 V. In front of typical ICRH antenna this RF potential can be as high as 1000 V 32,[START_REF]22th Topical Conference on Radio-Frequency Power in Plasmas[END_REF] , while the measured probe potential far from the antenna is of the order of 200 V 17,[START_REF] Kubic | Measurement of sheath potential in RF-biased flux tubes using a retarding field analyzer in Tore Supra tokamak[END_REF] , which seems to be consistent with the present model. On the plasma discharge side, figures 13 and 14 compare the plasma potential in Helium and Argon plasma respectively. It appears that the collapse of the plasma potential is triggered at a lower collision rate in Helium because the ratio j e /j i is more than 3 times lower than for Argon. Here the ratio ν/ω ci is always higher than 1 even in Helium, that's why ions are considered as unmagnetized.

A high ion-neutral collision frequency ν means a low perpendicular current, and then no reverse of the sheath potential, no collapse of the plasma potential which stays constant all over the computed range.

On the contrary, a low frequency means higher perpendicular current and a strong collapse (at ν = 1 MHz) of the plasma potential at an area ratio as low as 10 or 20 in Helium and Argon. This means that an increase of the RF potential does not always drive an increase of the plasma potential due to this saturation.

C. Both currents

The association of both resistive and capacitive currents yields both behaviors on the plasma potential evolution. First an increase in the first part of the plot and next a decrease due to higher perpendicular ion DC currents. According to the DC conductivity and capacitance, the plasma potential can either only decrease (low capacitive currents) or only increase or firstly increased and then decreased as presented in figure 15 and 16. All behaviors can be observed, making the understanding of plasma potential even more complex. In the plasma discharge case, i.e. figure 15, the double behavior appears at higher collisionality when the resistive perpendicular current is low enough compared to capacitive currents, and high enough to next make decrease the potential. In the same way, the tokamak case shown in figure 16 exhibits the same behavior except the first transition region is shorter, appearing for lower area ratios (S ⊥ /S <5). The collapse region is not as brutal as in plasma discharge at higher collision- ality, i.e. higher ion current. But in tokamaks, other perpendicular currents are allowed (inertia currents, viscous currents [START_REF] Rozhansky | [END_REF] , and even resonant ion cyclotron currents in the resonant layer of the ICRH antenna), which means the curve corresponding to ν = 10 7 Hz is more likely to be observed, meaning that for relatively long flux tubes connected to the RF source, the plasma potential can be more than 5 times lower than at the source point.

D. PIC comparison

This last section aims at validating previous analytical models with PIC simulations. They have been run in a 2D simulation box as presented in figure 9. The RF biased flux tube appears along the magnetic field and allows perpendicular capacitive or inductive currents, but resistive currents are not taken into account in the PIC simulation because Monte-Carlo collisions are not implemented. The area ratio S ⊥ /S is now equal to 2L/L 1 . In that case, the RF electrode mainly drains electron current while the other side and lateral boundaries of the biased plasma column drain ion currents. The plasma parameters are the same as in table II and the PIC simulation parameters are a little bit different than in previous sections. Now the simulation box is rectangular, 2.5 cm along x (magnetic field direction) and 5 cm along y. The grid is now composed of 512 by 1024 cells and the superparticle number is 2 [START_REF] Raizer | Radiofrequency Capacitive Discharges[END_REF] . The time step remains the same 5.10 -11 s. 460 seconds are needed to perform the 12000 time steps, corresponding to 12 RF periods. The goal here is to demonstrate that we do not need resistive currents to see the double behavior of the plasma potential. In figures 17, 18 and 19, it appears that the 2 behaviours can be seen especially in the tokamak plasma case (red curve), at higher magnetic field (B = 3 T). In the plasma discharge case, the increase region is smaller, meaning that capacitive or inductive currents are too high to let the plasma potential increase. The other effect appearing clearly in PIC simulations (see figures 20 and 21 for Helium gas) is the depletion of the RF biased flux tube 4,36,37 , due to electron parallel pumping and to ion perpendicular ejection. In figures 20 and 21, the average density profile along y-direction has been plotted for 3 electrode sizes (0.4, 0.12 and 0.006 the total width). It can be seen that the density in front of the electrode can be as low as a tenth of the initial density for B=0.1 T. This also means that the saturation of the electron current is reached much sooner, explaining the strong collapse of the plasma potential. At higher magnetic fields, i.e. tokamak case, the ion currents are smaller and due to smaller Larmor radius the biased plasma column is only depleted at the edges (figure 21), which means that the depletion is not as strong as in the low B case, leading to a smoother collapse of the plasma potential at higher area ratios; a high area ratio meaning thin and long-biased flux tube. Moreover, the area ratio for which the plasma potential collapses seems to be not very dependant on the ion mass in the PIC simulation. One can see that for Argon the collapsing transition occurs at just a little bit higher area ratio than for Helium or Deuterium (instead of a factor 3), but at lower magnetic fields (0.1 T) the transition occurs at area ratios more than 2 times lower than at higher magnetic field (3 T) for the same species. So the perpendicular current depends more on the magnetic field than on ion mass. As a perspective, this effect could be validated in our experimental magnetized plasma device Aline 38 for the plasma discharge case.

Finally, the PIC results are compared to DSP analytical model and nicely fit together as shown in figures 17, 18 and 19, on conditions that parameters α, κ and ν are well chosen. Here κ is a correction parameter in the perpendicular capacitance of the plasma column, c ⊥ = 1 ωci is replaced by c ⊥ = κ ωci . This correction factor is of the same nature as α for the parallel sheath capacitance. This factor takes into account the geometric shape of the sheath which depends non-linearly on the plasma potential and density oscillations and the distortion of the supposed Maxwellian distribution for electrons. The tokamak and discharge simulations are fit with the corresponding parameters (conduction and capacitive coefficient) presented in table II.

V. CONCLUSION

The classical double probe model able to calculate the plasma potential of an asymmetric plasma discharge has been modified to take into account the electron saturation current. This can result in a collapse of the plasma potential for the highest asymmetric discharges. The asymmetry is characterized by the wall/electrode ratio from 1 (symmetric) to more than 100. The new model is called the double saturated probe (DSP) model because it takes into account both ion and electron saturation currents. To do so, the exponential is replaced by a hyperbolic tangent. The simplest expression of the model has been solved analytically in RF and DC discharges and shows a very brutal fall of the plasma potential as soon as the area ratio is higher than the electron over ion saturation current ratio. These results have been compared to PIC simulations in which the collapsing region is not as sharp and is more extended over higher ratios. But solving the DSP model with a modified electrode area depending on the sheath thickness allows to almost fit the PIC curves.

Next, the model has been enriched with capacitive sheaths, adding capacitive currents within the sheath. The model has been solved numerically and shows the same results as Aanesland's model in the case of a high capacitive sheath. Again the collapse of the plasma potential occurs at the corresponding electron over ion saturation current ratio in Argon or Helium. The main difference with Aanesland's model appears in the first part of the plot, for low area ratio (lower than 5), in which the plasma potential first rises, starting from V rf /π to a maximum value which can be close to the RF potential before collapsing at higher area ratio.

Finally, the same model has been adapted to a magnetized plasma column with perpendicular capacitive currents. The same behavior is observed in plasma discharge and tokamak plasma : first, an increase of the potential and then the collapse at a higher ratio. This occurs at area ratios depending more on the magnetic field strength than on the ion species mass, contrary to the unmagnetized model. In the magnetized case, the area ratio is defined as the perpendicular current area (the area of the plasma column) over the parallel current area (the electrode area). It is like an asymmetric discharge and then the plasma potential for long plasma columns can be much lower than the RF source amplitude which can be more than 5 times lower. This study is useful to understand that it is impossible to bias the plasma at very high potential without saturation due to current asymmetries. It also shows that in tokamak plasmas the high RF potential induced by the slow wave of ICRH antennas is strongly mitigated far from the antenna.

Figure 1 .

 1 Figure 1. 1D planar probe I-V characteristic and the corresponding "tanh" I-V characteristic.

Figure 2 .

 2 Figure 2. Sketch of the 2D plasma box used in the PIC simulation.

Figure 3 .

 3 Figure 3. Electric model of an asymmetric plasma discharge.

Figure 4 .

 4 Figure 4. Plasma potential (in Volts) as a function of the ground/electrode area ratio S2/S1. DSP model (dashed and dotted lines) is compared to PIC simulation results (solid lines with markers).

Figure 5 .

 5 Figure 5. Electric equivalent circuit of the asymmetric double probe capacitive model

Figure 6 .

 6 Figure 6. Plasma potential with respect to the wall vs RF electrode ratio, respectively S2/S1. The applied RF amplitude is 100 V.

Figure 7 .

 7 Figure 7. Helium plasma potential Vp with respect to wall/electrode area ratio in an unmagnetized plasma for PIC, Aanesland's and DSP evaluations.

Figure 8 .

 8 Figure8.Argon plasma potential Vp with respect to wall/electrode area ratio in an unmagnetized plasma for PIC, Aanesland's and DSP evaluations.

Figure 9 .Figure 10 .

 910 Figure 9. Sketch of the 2D magnetized plasma box with an applied RF voltage

Figure 11 .

 11 Figure 11. Plasma potential with respect to perpendicular area over parallel area, respectively S ⊥ /S . 3 different cases are presented (dicharge plasma in He and Ar, and SOL tokamak plasma). The applied RF amplitude is 100 V. Only capacitive currents are taken into account.

Figure 12 .Figure 13 .Figure 14 .

 121314 Figure12. Plasma potential with respect to perpendicular/parallel area ratio, respectively S ⊥ /S in a Tokamak plasma. The applied RF amplitude is 100 V. Only DC collisional (ion-neutral) currents are taken into account for a SOL tokamak plasma.

Figure 15 .

 15 Figure 15. Plasma potential as function of the perpendicular/parallel area ratio for several ion-neutral collision frequencies ν in a Helium plasma discharge.

Figure 16 .

 16 Figure 16. Plasma potential as function of the perpendicular/parallel area ratio for several ion-neutral collision frequencies ν in a tokamak plasma.

Figure 17 .

 17 Figure 17. PIC simulations of the plasma potential as function of the perpendicular/parallel area ratio in discharge and tokamak Argon plasmas.

6 Figure 18 . 1 Figure 19 . 5 . 10 7

 618119510 Figure 18. PIC simulations of the plasma potential as function of the perpendicular/parallel area ratio in discharge and tokamak Helium plasmas.

Figure 20 .Figure 21 .

 2021 Figure 20. Plasma density at low magnetic field (B=0.1 T) for Helium.

Table I .

 I Typical plasma discharge parameters used in figures 6, 7 and 8.

	gas	He, Ar
	n0	10 16 m -3
	Te	3 eV
	V rf	100 V
	S2/S1	1 → 250
	RF freq.	20 MHz
	λ

sh 10 -4 → 10 -2 m c λ 2 De Csλ sh c rf S 1 λ 2 De (S 1 +S 2 )Csλ sh

Table II .

 II Typical plasma parameters for a RF discharge or tokamak plasma in the "Scrape off layer" (SOL).

		RF discharge Tokamak
	gas	He, Ar	D
	n0 (m -3 )	10 16	10 18
	Te (eV)	3	20
	B (T)	0.1	3
	ν (Hz)	10 6 -10 8 10 4 -10 6
	V rf (V)	100	"
	S2/S1	1 → 100	"
	RF freq.	20 MHz	"
	perp. scale λ sh or ρci	ρci
	γ	Cs νλ sh or ν ω

VI. ACKNOWLEDGEMENTS

The authors would like to thank Laurent Marot, Paul Hiret, Rodrigo Antunes and Roland Steiner from Basel university for their usefull suggestions and discussions. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work is also supported by the French National Research Agency (ANR) under project SHEAR ANR-19-CE30-033-01.

VII. DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.


domain: sde spi

text:
Dimassi

H Mhamdi

K Azaiez

T Fiorido

R Benabderrahmane Zaghouani

J L Lazzari

M Bendahan

W Dimassi

Rabia Benabderrahmane
email: rabia.benabderrahmane@gmail.com

Room temperature NO2 gas sensor based on stain-etched porous silicon: towards a low-cost gas sensor integrated on silicon

Keywords: Nano-Porous silicon, Macro-porous silicon, stain-etching, NO2 detection, room temperature

Room temperature NO2 gas sensor based on stain-etched porous silicon: Towards a low-cost gas sensor integrated on silicon

Introduction

Several materials with large specific surface areas, such as polymer membranes [START_REF] Bai | Gas Sensors Based on Conducting Polymers[END_REF], carbon nanotubes [START_REF] Slobodian | Multi-wall carbon nanotube networks as potential resistive gas sensors for organic vapor detection[END_REF], and metal oxides [START_REF] Sertel | Development of MgO: TiO2 thin films for gas sensor applications[END_REF][START_REF] Oliveira | Enhancement of acetylene gas sensing properties for ZnObased gas sensor produced by plasma immersion ion implantation and deposition[END_REF][START_REF] Barbosa | Investigation of electronic and chemical sensitization effects promoted by Pt and Pd nanoparticles on single-crystalline SnO nanobelt-based gas sensors[END_REF] have been studied as prospective candidates for future gas sensors. Actually, the most used gas sensors are based on metal oxides whose principle is based generally on chemisorption interactions which are oxydo-reduction reactions based on the exchange of electrons between the gas and the sensitive layer. The reactions take place on the surface on specific adsorption sites thermodynamically favorable. These reactions depend on the gas type and adsorption sites. It has been shown that oxygen is necessary to catalyze the sensing reaction in which the gas reacts generally with the adsorbed oxygen molecules. Currently, in the most developed sensors, the sensitive layer is nanostructured, which considerably increases the number of adsorption sites. In this work, one way of improvement consists of increasing the sensor specific surface which can be achieved by using a porous silicon substrate. Porous silicon (PS), an interesting material with large specific surface area, is easily produced by chemical [START_REF] Mhamdi | Study of n-WO3/p-porous silicon structures for gas-sensing applications[END_REF][START_REF] Azaiez | Enhancement of porous silicon photoluminescence property by lithium chloride treatment[END_REF] or electrochemical [START_REF] Azaiez | Enhancement photoluminescence property of porous silicon treated with bismuth (III)[END_REF] processes. It is indeed a very variable material who's optical, electrical and morphological properties could be modulated by porosity and thickness. It presents a high specific surface up to 1500 m 2 /cm -3 . Chemical functionalization of the porous silicon surface shows the potential for developing a variety of gas sensors [START_REF] Mhamdi | Study of n-WO3/p-porous silicon structures for gas-sensing applications[END_REF]. Recently, PS gas sensor has been investigated especially for NO2 detection [START_REF] Hui-Qing | The light-enhanced NO2 sensing properties of porous silicon gas sensors at room temperature[END_REF][START_REF] Bdaiwi | Fabrication of Gas Sensor Device for H2 and NO2 from Porous Silicon[END_REF]. Although PS sensors show good sensitivity when exposed to NO2, some problems such as selectivity, repeatability, response-recovery characteristics and long-term stability need to be considered. In the literature, some research works have reported incomplete recovery, and irreversible behaviors limiting the PS integration in final gas-sensing devices [START_REF] Gaburro | Role of microstructure in porous silicon gas sensor for NO2[END_REF][START_REF] Massera | Improvement of stability and recovery time inporous-silicon-based NO2, sensor[END_REF]. In addition, the use of multilayer structures with different morphologies based on porous silicon (different pore sizes, porosity) obtained on the same substrate has been investigated for NO2 gas sensor [START_REF] Bolotov | Formation of Multilayer Structures with Integrated Membranes Based on Porous Silicon[END_REF]. In other reports, porous silicon membranes are investigated in gas sensors using hybrid and composite materials leading to the selectivity enhancement [START_REF] Bolotov | Formation of Multilayer Structures with Integrated Membranes Based on Porous Silicon[END_REF]. The PS with ordered pore channels may influence considerably the gas-sensing properties, including reversibility and response-recovery time. Moreover, porous silicon sensing properties could be ameliorated by catalyst deposition and modifying PS surface [START_REF] Galstyan | Investigations of hydrogen sensors made of porous silicon[END_REF].

However, the porous silicon presents unstable Si-H bonds on its surface which are rapidly degraded under ambient conditions leading to surface oxidation and therefore the interface states and dangling bonds presence.

The formation of SiOx layer is certainly playing an important role on gas detection.

In this work, we present an original study focused on the use of two different types of porous silicon (PS) obtained by a simple and low-cost method as a gas senor at room temperature: nano-porous silicon (nano-PS) and macro-porous silicon (macro-PS). The PS gas sensing properties and response/recovery characteristics at room temperature are investigated and discussed.

Experimental details

Samples preparation :

Nano porous (nano-PS) and macro porous (macro-PS) silicon samples were synthesized via stain etching technique during 2 min using etching solution (HF, HNO3, H2O) = (1, 2, 3) in volumetric ratio (Fig. 1(a)). P-type and N-type silicon substrates with a resistivity of ̴ 1-3 Ω.cm and a thickness of 200 µm were used to elaborate nano-PS and macro-PS samples respectively. To elaborate macro-PS samples, light is used during the etching process in order to create holes necessary for the etching process.

Before PS elaboration, the substrates have followed different cleaning steps to eliminate organic greases by immerging samples in boiling acetone during 10 min then in ethanol during 5 min and deionized water. Finally, the native oxide was etched in 10 % acid hydrofluoric during 1 min, rinsed in deionized water and dried with nitrogen. After porous silicon formation, aluminum (Al) electrodes are deposited on PS surface by thermal evaporation with a thickness of about 100 nm in a parallel configuration.

Physicochemical samples characterization:

The formed porous silicon was studied using a Fourier Transform Infrared (FTIR) Nicolet MAGNA-IR 560 Spectrometer. The PS morphology was investigated by a scanning electron microscope (SEM). The elemental analysis of PS surface was performed by the X-ray Dispersive Energy (EDX) coupled with SEM. Effective minority carrier lifetime (τeff) was studied using a microwave induced photoconductive decay technique (μW-PCD). The measurements are performed by Semilab WT-2000 system using a 904 nm laser excitation source.

Gas sensing measurement:

The gas-sensing performances of the PS sensors were tested in a dynamic gas sensing testing system, as shown in Fig. 1 (b). The sensor was placed on a heating plate connected to a temperature controller to adjust the operating temperature. By measuring the resistance change over time between two electrodes under gas, the sensor dynamic response curves were collected. In this work, the gas sensor response is defined as: Ra/Rg (When sensor is exposed to an oxidizing gas) (Eq. 1-a)

Or Rg/Ra (When sensor is exposed to a reducing gas) (Eq. 1-b)

Where Ra and Rg are the sensor resistances measured in dry air and under tested gas, respectively. The PS sensors are tested toward exposed gas at room temperature (RT) (25 ± 2 °C). The response and recovery times are also calculated. The response time is measured as the time necessary for the resistance to reach 90 % of the equilibrium value after gas exposure. The recovery time is obtained as the time required reaching 10 % above initial resistance after gas removal. The samples were exposed to different gas as NO2, NH3, O3 and CO. Ozone is generated from dry air, using a system based on a UV lamp. The dry air passes through a quartz tube illuminated by the UV lamp. Exposure of oxygen molecules to UV rays leads to the generation of ozone molecules.

Results and discussion:

Porous silicon characterization:

For silicon nanostructuring, there are two possible approaches: the top-down and the bottom-up approaches. In this work, porous silicon is obtained by the top-Down approach consisting in miniaturizing the silicon substrate until obtaining nanometric structures by etching the initial substrate in an acidic solution. The stain-etching process is composed of two chemical reactions occurring simultaneously: the silicon surface oxidation by the nitric acid and the silicon oxide etching by hydrofluoric acid. Porous silicon is characterized by the morphology, the porosity and the thickness. Porosity is defined as the fraction of voids within the porous structure [START_REF] Yang | Quantifying biofilm structure using image analysis[END_REF]. The quantification of the material loss in the etching process is used as a gravimetric measurement to evaluate the PS porosity and thickness [START_REF] Zaghouani | Study of WO3-decorated porous silicon and Al2O3-ALD encapsulation[END_REF]. The following equations are used to obtain the average porosity (P) (Eq. 2). The sample was weighted before porosification (m1), after porosification (m2), and after porous silicon dissolution in 1 wt.% NaOH solution (m3).  is the density of silicon and S is the etched area. The porosity of the elaborated porous silicon is about 38% and 66 % for nano-PS and macro-PS respectively. Fig. 2 shows the nano-PS surface morphology obtained by SEM microscopy. The top surface SEM images (Fig. 2(a-b)) shows the presence of pores with nanometric sizes less than 10 nm confirming the synthesis of nano porous structure. The PS layer thickness estimated thanks to the cross-section SEM image (Fig. 2(c)) is about 0.64 µm. Fig. 3(a-b)

shows the macro-PS surface morphology. We distinguish the presence of pores with sizes more than 50 nm proving the formation of macro-porous silicon layer with a thickness of 1.6 µm (Fig. 3(c)). The X-ray Dispersive Energy (EDX) coupled to the SEM performed on the nano-PS porous silicon sample shows that no metallic impurities related to the synthesis method were detected (Fig. 2(d)). The presence of an oxygen peak indicates the surface oxidation. In fact, the hydrogenation via HF and the formation of Si-H bonds is only stable during low durations. Stored in ambient atmosphere, the Si-H bonds are rapidly degraded leading to the porous surface oxidation and therefore the presence of dangling bonds which could play an important role in gas sensing mechanism. In order to confirm this observation, FTIR spectra of freshly prepared and oxidized PS samples were performed (Fig. 4) showing the vibration bands attributed to SiHx and Si-O-Si stretching mode at 2112 cm -1 and 1092 cm -1 respectively, the SiH2 scissors mode (907 cm -1 ), the Si-Si stretching mode (628 cm -1 ), SiHn and O-SiHx wagging mode at 672.7 cm -1 861 cm -1 respectively [START_REF] Alaya | Enhancement of physical properties of stain-etched porous silicon by integration of WO3 nanoparticles[END_REF]. Comparing the two spectra, we notice the reduction of the SiHx peak at 2112 cm -1 against the increase in Si-O peak intensity confirming the PS surface oxidation.

5(c)

) and macro-PS sample (Fig. 5(d)). As we can notice, the lifetime in freshly prepared nano-PS layer is homogeneous and greater than that in silicon thanks to the surface hydrogenation during the etching process in HF/HNO3 solution. For the macro-PS sample, a slight enhancement in τeff is also obtained. We note also that the average τeff measured in the macro-PS (13µs) is greater than that of the nano-PS sample (1.9 µs) attributed to the quality of the initial silicon substrate used.

Gas sensing properties:

The gas sensing properties of the nano-PS and macro-PS sensors towards 4 ppm of NO2 gas were studied at room temperature in dry air. Fig. 6(a) shows the variation of the electrical resistance as a function of time with a temperature increasing from 25 to 50 °C. The results show that the resistance decreases with increasing the temperature proving the semiconductor behavior of the PS sensor. The nano-PS gas sensor response as a function of time at room temperature is shown in Fig. 6 (b). This sensor has a response about 1.5 at ambient temperature, as well as reasonable response (45 s) and recovery times (67 s). The gas sensor based on nano-PS is therefore operational in ambient temperature. The recovery time is considered interesting thanks to the low thickness of the PS layer (about 640 nm). Fig. 6 (c) shows the good repeatability of the nano-PS sensor for three dynamic response cycles for an exposure time of 1 min and for two dynamic cycles for an exposure time of 3 min. For each test, complete recovery of initial resistance is observed. This indicates that porous silicon has a chemical stability when it comes to NO2, the gas particles are able to be adsorbed, and then desorbed rapidly at room temperature once the experiment is over. The gas selectivity of porous silicon sensors is also studied. Fig. 6(d) shows the sensor response under different gas (4 ppm NO2, 4 ppm NH3, 200 ppb O3 and 4 ppm CO) at room temperature. As noted, the best response is recorded to NO2 however the sensor response towards NH3 and O3 is very low. When exposed to CO, the sensor does not display any sensitivity. The CO gas is still presenting problems of detection as few research works have reported on CO sensing in particular at room temperature. These results show that NO2 is the most compatible gas for this sensor. Table 2 recapitulates the different characteristics of the nano-PS sensor.

For the different detected gas, the nano-PS sensor exhibits a rapid response proving its efficiency while the highest recovery time is obtained when it is exposed to NH3.

The gas sensing properties of the macro-PS sensors towards 4 ppm of NO2 gas were studied at room temperature. Fig. 7 shows the macro-PS gas sensor response as a function of time. The sensor had a response of 2.35 at room temperature, much higher than nano-PS sensor response. The response and recovery times of this sensor were determinated and the values were about 64 s and 45 min, respectively. This sensor shows a very slow recovery characteristic attributed to the high thickness of the PS layer (about 1.6 µm). The selectivity of macro-porous silicon sensors is also studied. Fig. 7(b) shows the sensor response under different gas (4 ppm NO2, 4 ppm NH3, 200 ppb O3 and 4 ppm CO) at room temperature. The macro-PS sensor exhibited the important gas sensing characteristic for NO2 gas. Table 3 recapitulates the different characteristics of the macro-PS sensor.

We can consider that NO2 gas is the most compatible gas for this sensor.

Compared to the nano-PS sensor results, the increase in sensitivity for the macro-PS sensor is attributed to the obtained morphology (Fig. 2 and Fig. 3). The nano-PS presents nano-sized pores less than 10 nm however the macro-PS has macro-pores. Besides, the macro-PS sample has a higher porosity justifying that a higher quantity of gas molecules is adsorbed in the macro-PS than in the nano-PS sensor. Comparing the recovery times measured for the two sensors, we note that the gas desorption takes much more time with macro-PS sensor attributed to the porous layer thickness (1.6 µm for macro-PS against 0.64 µm for the nano-PS).

Sensing mechanism:

The response of most gas sensors is due to the changes in carrier concentration, usually induced by adsorptiondesorption reactions on the sensing materials surface. In recent years, many researchers have investigated the NO2-sensing mechanism of porous silicon [START_REF] Gaburro | Opposite effects of NO2 on electrical injection in porous silicon gas sensors[END_REF][START_REF] Boarino | Towards a deeper comprehension of the interaction mechanisms between mesoporous silicon and NO2[END_REF]. It is found that the PS belongs to the surface-controlled type as both surface states and oxygen adsorption play important roles in the NO2 detection [START_REF] Sun | Microstructure, electrical and gas sensing properties of meso-porous silicon and macro-porous silicon[END_REF]. NO2 acts as an oxidizing gas, which will capture the electrons from the sensing material. In the first step, the atmospheric oxygen molecules, adsorbed onto the PS surface in the form of O -and O2 -depending on the operational temperature, are extracting electrons from the PS conduction band ((Eq. 7)-(Eq. 9)). In the case of the n-type porous silicon sensor, it leads to the formation of electron-depleted space-charge layer under the PS surface [START_REF] Ali | Effects of surface passivation in porous silicon as H2 gas sensor[END_REF], and then the resistance increase.

O2 (gas)→ O2 (ads)

(Eq. 7)

O2 (ads) +é→ O2 -(ads), T<100 °C (Eq. 8)

O2 -(ads) +é→ 2O -(ads), 100 °C<T<300 °C (Eq. 9)

In the second step, when the PS sensor is exposed to NO2, the gas molecules may be directly adsorbed onto the PS surface by trapping electrons (Eq. 10), and interact with the chemisorbed oxygen (such as [START_REF] Gaburro | Role of microstructure in porous silicon gas sensor for NO2[END_REF] and (Eq. 12)) [START_REF] Ghimbeu | Electrostatic spray deposited zinc oxide films for gas sensor applications[END_REF].

- O and - 2 O ) ((Eq.

NO2 (gas) + é → NO2 -(ads)

(Eq. 10) NO2 (gas) +O2 -(ads) + 2é → NO2 -(ads) + 2O -(ads) , T<100 °C (Eq. 11) NO2 (gas) +O -(ads) + 2é → NO2 -(ads) + O 2-(ads) , 100 °C<T<300 °C (Eq. 12)

The reaction with NO2 will increase further the electron depletion layer width and potential barriers, and then increase the sensor resistance [START_REF] Stievenard | Are electrical properties of an aluminium-porous silicon junction governed by dangling bonds[END_REF]. Consequently, NO2 molecules capture the electrons from PS surface, leading to the decrease of free electron concentration [START_REF] Konstantinova | Influence of NO2 molecule adsorption on free charge carriers and spin centers in porous silicon[END_REF] leading to the sensor resistance increase. The PS electrical conductivity is significantly modulated depending upon the presence or the absence of NO2, explaining the resistance variations. Fig. 8 presents the NO2-sensing mechanism diagram of the n-type porous silicon at room temperature. Initially, the sensing material presents a flat band structure (Fig 8 (a)). At room temperature, the oxygen molecules are chemisorbed in the form of O2 -, leading to the formation of a depletion layer (Fig. 8(b)).

This formed layer would be raised as a result of NO2 gas molecules interaction with the chemisorbed oxygen (Fig. 8 (c)).

However, when n-type PS is exposed to a reducing gas (electron donor) such as ammonia (NH3), the gas molecules adsorbed on the surface are interacting with the chemisorbed oxygen molecules. This reaction leads to the release of electrons trapped in the surface of oxygen and their injection into the porous silicon conduction band [START_REF] Mhlongo | A highly responsive NH3 sensor based on Pd-loaded ZnO nanoparticles prepared via a chemical precipitation approach[END_REF]. Thereby, the PS resistivity is decreasing. The same explanation could be conducted with p-type porous silicon.

Conclusion :

In summary, two types of porous silicon (nano-PS and macro-PS) were prepared for NO2 detection through simple chemical method. It was found that the prepared PS is promising for low-concentration NO2 sensing at room temperature. Using the same technique and synthesis parameters, the two prepared sensors exhibit different behaviors. The substrate doping, and accordingly the obtained morphology, should be accounted for the PS gas sensors. In this study, we obtained responses towards low concentration of NO2 at room temperature so the operational temperature is no longer an issue when using porous silicon as a sensing material. Therefore with porous silicon, the heater is no longer an essential component of the system.

Figure Caption

Fig. 5

 5 Fig. 5 presents the 2D mapping image and the corresponding histogram with 250 µm resolution showing the minority carrier lifetime in p-silicon substrate (Fig. 5(a)), nano-PS sample (Fig. 5(b)), n-silicon substarte (Fig.

Figure 1 :

 1 Figure 1: (a) Schematic illustration of the adopted porous silicon synthesis process. (b) The schematic diagram of the static gas-sensing test system.

Figure 2 :

 2 Figure 2: SEM images of nano-PS sample: (a-b) top surface and (c): cross section images (d) EDX diagram of porous silicon.

Figure 3 :

 3 Figure 3: SEM images of macro-PS sample: (a-b) top surface and (c): cross section images.

Figure 4 :

 4 Figure 4: FTIR absorption spectrum of porous silicon PS sample.

Figure 5 :

 5 Figure 5: The 2D mapping image of the minority carrier lifetime in (a) p-silicon, (b) Nano-porous silicon (nano-PS), (c) n-silicon and (d) Macro-porous silicon (macro-PS).

Figure 6 :

 6 Figure 6: (a) Resistance variation of nano-PS as a function of time under dry air with simultaneous increase of temperature from 25 °C to 50 °C (b) Dynamic response and (c) Resistance variation of the nano-PS sensor exposed to 4 ppm of NO2 (d) Selectivity towards NO2, NH3, O3 and CO at room temperature.

Figure 7 :

 7 Figure 7: (a) Dynamic response of macro-PS as a function of time under dry air, (b) Selectivity of macro-PS towards NO2, NH3, O3 and CO at room temperature.

Figure 8 :

 8 Figure 8: The schematic NO2-sensing mechanism of n-type porous silicon.

  Figure 4

Figure 5 Figure 6

 56 Figure 5

TableTable 1 :

 1 Comparative of morphological and electronic parameters of nano-PS and macro-PS samples.

Table 2 :

 2 The gas sensing properties of the nano-PS sensor towards different gases at room temperature.

Table 3 :

 3 The gas sensing properties of the macro-PS sensor towards different gases at room temperature.

Acknowledgements:

This work was supported by the bilateral Tuniso-French PHC-UTIQUE project (Code CMCU: 20G1108-Code Campus France: 44146SM). The authors wish to thank Damien Chaudanson and Alexandre Altié for SEM measurements done at the Centre Interdisciplinaire de Nanoscience de Marseille (CINaM).

Conflict of Interest Statement:

The authors declare that they have no conflict of interest.