T N Tran

I O Golosnoy
email: i.golosnoy@soton.ac.uk

P L Lewin

G E Georghiou

Numerical modelling of negative discharges in air with experimental validation

Keywords: corona discharge, surface discharge, finite element simulation, needle plane, Trichel pulses

Axisymmetric finite element models have been developed for the simulation of negative discharges in air without and with presence of dielectrics. The models are based on the hydrodynamic drift-diffusion approximation. A set of continuity equations accounting for the movement, generation and loss of charge carriers (electrons, positive and negative ions) are coupled with Poisson's equation to take into account the effect of space and surface charges on the electric field. The model of a negative corona discharge (without dielectric barriers) in a needle--plane geometry is analysed first. The results obtained show good agreement with experimental observations for various Trichel pulse characteristics. With dielectric barriers introduced into the discharge system, the surface discharge exhibits some similarities and differences to the corona case. The model studies the dynamics of volume charge generation, electric field variations and charge accumulation over the dielectric surface. The predicted surface charge density is consistent with experimental results obtained from the Pockels experiment in terms of distribution form and magnitude.

1.

Introduction Non-thermal plasmas produced from gas discharges have received a lot of attention in recent years due to their broad applications and practical issues [START_REF] Georghiou | Numerical modelling at atmospheric pressure leading to plasma production[END_REF][START_REF] Kogelschatz | Dielectric-barrier discharges: their history, discharge physics, and industrial applications[END_REF]. In the absence of a dielectric, corona discharges have been utilised for the treatment of gaseous pollutants, surface treatment, lighting and electrostatic painting [START_REF] Georghiou | Numerical modelling at atmospheric pressure leading to plasma production[END_REF][START_REF] Penetrante | Non-thermal plasma techniques for pollution control[END_REF][START_REF] Stoffels | Plasma needle: a non-destructive atmospheric plasma source for fine surface treatment of (bio)materials[END_REF]. Corona discharges can also lead to serious problems such as arcs in RF processing or power losses in transmission lines [START_REF] Georghiou | Characterisation of point-plane corona in air at radio frequency using a FE-FCT method[END_REF]. In the presence of a dielectric(s), the discharges are often regarded as dielectric barrier discharges (DBD) or silent discharges. At atmospheric pressure, DBDs have found their applications in large-scale ozone production, plasma display panels, thin film deposition, actuators for flow control, CO 2 lasers and excimer ultra-violet light sources [START_REF] Georghiou | Numerical modelling at atmospheric pressure leading to plasma production[END_REF][START_REF] Lagmich | Model description of surface dielectric barrier discharges for flow control[END_REF]. Additionally surface charge accumulation at triple-junctions in gas-insulated switchgear leading to breakdown of insulation system has also been a topic of research for a while [START_REF] Murooka | Nanosecond surface discharge and charge density evaluation, part i: review and experiments[END_REF].

Depending on the polarity of the applied voltage, the plasma structure varies from filamentary to diffuse patterns. The latter is produced from a negative applied voltage and is often preferred for certain applications due to its stability, uniformity and efficiency [START_REF] Georghiou | Numerical modelling at atmospheric pressure leading to plasma production[END_REF]. The negative corona discharge has been well studied experimentally and the most distinctive feature is often regarded to as Trichel pulses [START_REF] Trichel | The mechanism of the negative point to plane corona near onset[END_REF]. Numerical modelling work on this subject has mostly been limited to 1D-1.5D models [START_REF] Morrow | Theory of negative corona in oxygen[END_REF][START_REF] Espel | Numerical modelling of negative corona discharge in SF 6 at atmophseric presure[END_REF][START_REF] Soria-Hoyo | Particle-in-cell simulation of Trichel pulses in pure oxygen[END_REF]. Regarding dielectric barrier discharge, most simulation work have concentrated on the positive discharge in a needle--plane geometry while attention on the negative polarity is not as extensive. Thus the aim of this paper is to develop 2D axisymmetric finite element models for the negative corona and barrier discharges and validate the model parameters by comparing with experimental data. The hydrodynamic drift-diffusion approach is utilised which consists of coupled continuity equations for the charge transport and Poisson's equation for the electric field.

Experimental studies

Figure 1 shows the schematic diagram of the experimental setup used to measure the Trichel current pulses in air at atmospheric pressure. A needle electrode of 250 µm radius of curvature was utilised and connected to a negative dc high voltage supply. The needle was placed 3.3 mm away from the plane electrode which is made of stainless steel and has a radius of 10 cm and thickness of 5 mm. The plane electrode was grounded via a shielded low inductance 75 Ω shunt resistor. The voltage across the shunt resistor is transmitted to a digital oscilloscope via a 75 Ω impedance co-axial cable and terminated using a 75 Ω termination to eliminate reflections. The entire experiment was housed in a Faraday cage to minimise external interferences. With regard to measurements of surface charge accumulation over solid dielectrics, various techniques have been developed over the years including dust figures, Lichtenberg figures, electrostatic probes and the Pockels technique [START_REF] Murooka | Nanosecond surface discharge and charge density evaluation, part i: review and experiments[END_REF][START_REF] Sam | Surface discharge measurements of polymeric materials[END_REF]. Out of these methods, the electro-optic Pockels technique is the only method that offers quantitative analysis of surface discharge dynamics. The schematic diagram of the Pockels experiment is shown in figure 2. The principle of the method is based on the Pockels characteristic of some materials such as Bisthmuth Silicon Oxide (BSO). Under the electric field caused by surface charge, the refractive indexes of the BSO crystal along two perpendicular principal axes (so-called fast and slow axes) are altered. The change in the refractive indexes will lead to a linear phase retardation of the incident light. The change in phase is linearly proportional to the internal electric field. By measuring the light intensity of the reflected light, the phase retardation can be determined and hence the surface charge density. Further details of the experiment can be found in [START_REF] Sam | Surface discharge measurements of polymeric materials[END_REF]. This method has been used to measure surface charge density produced by applying a negative square voltage waveform to a 25 µm radius of curvature needle electrode placed 500 µm away from the Pockels cell.

3.

Model description

Governing equations

Based on the well-known hydrodynamic drift-diffusion model [START_REF] Georghiou | Numerical modelling at atmospheric pressure leading to plasma production[END_REF], the set of continuity equations (to account for the movement, generation and loss of electrons, positive and negative ions and for the development of the space charge) coupled with Poisson's equation (to include the effect of space charge on the electric field) is given as

N N W N W N t N ∇ - ∇ - - - = ∂ ∂ . β η α (1) ( ) p p p p ep p e np p n e e p N D W N N N N N W N t N ∇ - ∇ - - - = ∂ ∂ . β β α (2) ( ) n n n n np p n e e n N D W N N N W N t N ∇ - ∇ - - = ∂ ∂ . β η (3) ( ) ( ) 0 . 0 r p e n e V N N N ∇ ε ∇ + - - = ε (4)
where t is time, e the electronic charge, ε 0 the vacuum permittivity and ε r the relative permittivity of the medium, V the electric potential; N e , N p and N n the electron, positive and negative ion concentrations respectively; e W , p W and n W the electron, positive and negative ion drift velocities; α, η, β, D e , D p and D n the ionisation, attachment, recombination and diffusion coefficients for electrons and positive and negative ions respectively. The swarm parameters, which are functions of the local electric field, V E -∇ = , were extracted from the approximations published in [START_REF] Kang | Numerical study on influences of barrier arrangements on dielectric barrier discharge characteristics[END_REF] and summarised in [START_REF] Georghiou | Numerical modelling at atmospheric pressure leading to plasma production[END_REF] with an addition of the recombination coefficient equal to 2×10 -13 (m 3 .s -1 ) as commonly used in the literature [START_REF] Georghiou | Numerical modelling at atmospheric pressure leading to plasma production[END_REF].

The plasma kinetic model is simplified in equations 1-4; only three sorts of charged particles: electrons, positive and abstract negative ions are taken into account. This model is, however, adequate as shown by many previous publications (see references in the review [START_REF] Georghiou | Numerical modelling at atmospheric pressure leading to plasma production[END_REF]). Particularly in the case of negative ions, their density is at least one order of magnitude lower compared to the positive ion and electron densities [START_REF] Georghiou | The theory of short-gap breakdown of needle point-plane gaps in air using finite-difference and finite-element methods[END_REF]. The complex ionisation, attachment and recombination kinetics were hence not considered as this would add considerably to the computational burden without affecting the results significantly. If other cases where the negative ions have a significant contribution to the development of the discharge were to be considered then the issue of complex plasma kinetics should be included (discharges in SF 6 , for example).

In the presence of a dielectric, zero current through the barrier and charge accumulation at the surface were assumed. The surface charge density can be calculated by integrating the normal component of charged particle fluxes at the surface as expressed in equation 5

where the charge fluxes through the surface

i i i i norm norm J e N W = r r , , , i e p n =
. After being deposited on the surface, the charge is to be stationary; any charge motion over the dielectric is neglected. It has long been established that secondary emission processes are the main mechanism sustaining negative discharges [START_REF] Trichel | The mechanism of the negative point to plane corona near onset[END_REF][START_REF] Morrow | Theory of negative corona in oxygen[END_REF]. These electron emission processes at the cathode include the Fowler-Norheim process (electrons emitted due to a high electric field at the cathode), photo-emission (electrons produced by photons generated in discharge plasma) and ion-impact emission (electrons released as positive ions strike the cathode surface). Previous studies [START_REF] Soria-Hoyo | Particle-in-cell simulation of Trichel pulses in pure oxygen[END_REF][START_REF] Gupta | Theory of step on leading edge of negative corona current pulse[END_REF] argue that ion-impact is the key mechanism sustaining the discharge. For this reason, the only secondary process that has been taken into account in this work is the ionimpact. By using the expression given in equation 6, this process may represent others by varying the ion-impact secondary emission coefficient, γ.

input flux of electrons

Different values of γ have been used in the literature and they often fall in the range from 10 -3 to 10 -2 [START_REF] Morrow | Theory of negative corona in oxygen[END_REF][START_REF] Espel | Numerical modelling of negative corona discharge in SF 6 at atmophseric presure[END_REF][START_REF] Soria-Hoyo | Particle-in-cell simulation of Trichel pulses in pure oxygen[END_REF][START_REF] Gupta | Theory of step on leading edge of negative corona current pulse[END_REF][START_REF] Georghiou | The effect of photoemission on the streamer development and propagation in short uniform gaps[END_REF].

Computational domain

The simulation geometries utilised in this work are shown in figure 3. In the negative corona discharge model, a hyperboloid-shaped needle electrode was assumed of 250 µm radius of curvature. The distance between the needle and the ground electrode was set at 3.3 mm. These were to match with the dimensions of the experimental setup. In the negative surface discharge model, a system of two dielectrics was modelled to resemble the Pockels cell. A 160 µm thick BSO layer (ε r = 56) was set on top of a 800 µm thick glass (ε r = 3). Both layers share the same radius of 10 mm. The BSO surface was set 500 µm away from the needle electrode. A different electrode of 25 µm radius of curvature was used to facilitate a higher electric field.

Discharge parameters

In order to compute the discharge current, two different methods have been utilised. The first approach is based on the energy conservation law which was detailed in [START_REF] O'sullivan | A model for the initiation and propagation of electrical streamers in transformer oil and transformer oil based nanofluids[END_REF] where the discharge current, I V , is calculated in a similar manner to [START_REF] Morrow | The discharge current induced by the motion of charged particles in time-dependent electric fields; Sato's equation extended[END_REF].

1 V D a V V d I QdV P dV V dt     = +     ∫ ∫ (7)
where V a is the external applied voltage, , ,

D i i i e p n P E e N W = ⋅ ∑ ur uur and 2 0 1 2 r Q E = ε ε ur .
Equation 7 is an extension of the well known Sato's expression [START_REF] Sato | Discharge current induced by the motion of charged particles[END_REF] to time-dependent voltage [START_REF] Morrow | The discharge current induced by the motion of charged particles in time-dependent electric fields; Sato's equation extended[END_REF]. Two expressions are detailed in [START_REF] Morrow | The discharge current induced by the motion of charged particles in time-dependent electric fields; Sato's equation extended[END_REF]; the first one uses the total field as equation 7 whereas the second one utilises the Laplacian field only. Here it is more convenient to use the first expression since all data required is readily obtained from the model. The other method to calculate the current is to investigate the charge transport across domain boundaries. The electric field is calculated in quasi-static approximation in the model. The fixed applied potential requires surface charges to be built on the electrode surfaces. The time variation of the total surface charge at the cathode is supplied by the external circuit (current I S ) or by the plasma fluxes, J terms in equation 8. The current is hence calculated by integrating the charged particle fluxes through the cathode surface plus the changes of the induced surface charge at the cathode.

( )

  = + + ⋅ + σ ⋅ ⋅     ∫ ur r ur r r r (8)
where the needle surface charge density is computed from Gauss' theorem ( ) ( ) is assumed due to metal electrode). It is worth noting that the current calculated at the ground electrode via equation 8 is consistent with that calculated from the cathode surface. This is as expected since no charges should be accumulated in the external power source.

0 ca r c m n E E σ = -⋅ ε ε -

Details of numerical algorithm

The Convection and Diffusion application mode in COMSOL was utilised to solve the continuity equations [START_REF] Comsol Inc | COMSOL Multiphysics simulation environment. Version 3.5a manual[END_REF]. Streamline artificial diffusion (anisotropic) was added to avoid oscillations and overshoots. COMSOL facilitates various forms of artificial diffusion when solving convection and diffusion problems [START_REF] Zienkiewicz | The finite element method for fluid dynamics[END_REF][START_REF] Codina | Comparison of some finite element methods for solving the diffusionconvection-reaction equation[END_REF]. Isotropic diffusion introduces numerical viscosity in all directions, streamline diffusion is in the direction of drift velocity only and crosswind diffusion is orthogonal to the streamline direction [START_REF] Zienkiewicz | The finite element method for fluid dynamics[END_REF]. Detailed comparisons of different artificial diffusion types have been presented in [START_REF] Codina | Comparison of some finite element methods for solving the diffusionconvection-reaction equation[END_REF]. It has been confirmed that streamline anisotropic diffusion gives relatively good results and reasonable run time efficiency on the test problems. Crosswind diffusion provides better accuracy for the test systems but runs 10 times slower. The streamline anisotropic diffusion has proved effective to solve discharge problems in gas and liquid [START_REF] O'sullivan | A model for the initiation and propagation of electrical streamers in transformer oil and transformer oil based nanofluids[END_REF][START_REF] Lee | Finite-element analysis of corona discharge onset in air with artificial diffusion scheme and under Fowler-Nordheim electron emission[END_REF] although other techniques for transport problems can be used [START_REF] Georghiou | An improved finite-element fluxcorrected transport algorithm[END_REF].

In addition, the Electrostatics mode was employed to solve Poisson's equation. The surface charge density variation over time, equation 5, was calculated by the integration of an additional differential equation. The time derivative in equations 5, 7 and 8 is linked with the time variation of the plasma electric field. The variable-order variable-step-size backward differentiation formulas (BDF) are used by COMSOL for time differentiation. (BDF is a further development of the multistep methods [START_REF] Gear | Numerical initial value problems in ordinary differential equations[END_REF].) It makes the current calculations consistent with time discretisation of the transport equations 1-3. Space integration is performed using the standard Gaussian quadrature (4 th order) built in COMSOL [START_REF] Comsol Inc | COMSOL Multiphysics simulation environment. Version 3.5a manual[END_REF].

Integration is undertaken on the same mesh used for the convection and diffusion problem.

Time steps in this model are automatically determined by COMSOL. The choice of time step is driven by the specified relative (0.01) and absolute (0.001) tolerances. The tolerances control the error in each integration step in a way that the step is accepted if the estimated local errors satisfy both tolerances [START_REF] Comsol Inc | COMSOL Multiphysics simulation environment. Version 3.5a manual[END_REF]. The time increment is chosen automatically according to the algorithm summarised in [START_REF] Gear | Numerical initial value problems in ordinary differential equations[END_REF]. As a result, the time steps vary around 10 -10 s during the pulse formation and increase to ~10 -9 s during the decay stage. Different tolerances have been tested by comparing the solutions with that obtained using fixed time steps (10 -12 s). It was found that the relative tolerance values of 0.01 and 0.001 (absolute tolerances were set 10 times smaller in both cases) give similar solutions (5% difference) to that calculated from fixed time steps. The use of larger relative (0.1) and absolute (0.01) tolerances resulted in a 20% difference hence was discarded.

4.

Negative corona discharge Experimental studies on this subject revealed a well-known phenomenon known as Trichel pulses [START_REF] Lama | Systematic study of the electrical characteristics of the "Trichel" current pulses from negative needle-to-plane coronas[END_REF]. This phenomenon occurs when a negative corona discharge takes place in an electronegative gas and the measured discharge current shows a regular pulse train. Important Trichel pulse parameters include the individual pulse shape, their magnitude and frequency of occurrence. Regarding the modelling side of the phenomenon, most of the work have been undertaken in 1D-1.5D (continuity equations are solved in 1D along the symmetry axis while Poisson's equation is solved using the effectively 2D disc method) [START_REF] Morrow | Theory of negative corona in oxygen[END_REF][START_REF] Espel | Numerical modelling of negative corona discharge in SF 6 at atmophseric presure[END_REF][START_REF] Soria-Hoyo | Particle-in-cell simulation of Trichel pulses in pure oxygen[END_REF][START_REF] Gupta | Theory of step on leading edge of negative corona current pulse[END_REF]. Reports on full 2D simulations are rather limited [START_REF] Akishev | Numerical simulations of Trichel pulses in a negative corona in air[END_REF]. Thus the objectives of this work are to enhance the understanding of the negative corona discharge phenomenon by the full 2D axisymmetric model, to get the model input parameters validated by experimental data and to set a benchmark to further proceed with the negative surface discharge model.

Boundary and initial conditions

The simulation geometry of this model is shown in figure 1(a) and the corresponding boundary conditions are summarised in table 1.

( ) 0 = ∇ - ⋅ p p N D n Concentration N n = 0 Potential V = -5500 (V) Boundaries 3&4 Convective flux ( ) 0 = ∇ - ⋅ e e N D n Convective flux ( ) 0 = ∇ - ⋅ p p N D n Convective flux ( ) 0 = ∇ - ⋅ n n N D n Zero charge/ Symmetry ( ) 0 0 = ⋅ E n r ε ε Boundary 5 Convective flux ( ) 0 = ∇ - ⋅ e e N D n Concentration N p = 0 Convective flux ( ) 0 = ∇ - ⋅ n n N D n Ground V = 0
In this model, a step voltage of -5.5kV was applied to the needle. A Gaussian distribution of neutral plasma (seed electrons and positive ions) with a peak value of 10 16 m -3 and width of 25 µm was set at the needle tip according to the following equation [START_REF] Tran | Two dimensional studies of Trichel pulses in air using the finite element method[END_REF]

        - - - - × = 2 0 2 0 2 0 2 0 max , 2 ) ( 2 ) ( exp s z z s r r N N p e (9)
where N max = 10 16 m -3 , r 0 = 0 µm, z 0 corresponds to the needle tip position and s 0 = 25 µm. The initial condition has been confirmed to only speed up the pulse formation and not to change the discharge characteristics.

The secondary emission coefficient γ was set to 0.004; as will be shown later this value gives the best agreement with experiment. During the development of the model, it was realised that major discharge activities occur in the region close to the cathode. For this reason, meshing around this region proves very important. A fine mesh of 1 µm in size was created in the volume 100 µm away from the needle axially and radially. Further away from the ionisation region, a coarser mesh was used. A total number of one million degrees of freedom was solved for the process time duration of 30 µs.

Simulation results

Details of the first discharge current pulse and of the pulse train are shown in figure 4. The current pulses have a magnitude of around 2 mA and frequency of occurrence of 142 kHz. It is interesting to note that the first current pulse magnitude is not larger than the subsequent ones as opposed to that observed when the code was used for a shorter gap [START_REF] Tran | Two dimensional studies of Trichel pulses in air using the finite element method[END_REF]. This can be explained by the fact that for a longer gas gap, the frequency of occurrence is lower hence space charge particles have sufficient time to clear away from the gap [START_REF] Lama | Systematic study of the electrical characteristics of the "Trichel" current pulses from negative needle-to-plane coronas[END_REF]. As a result, the subsequent current pulses occur with similar conditions to the first one. The model suggests that the first pulse only has a single peak as opposed to double peaks claimed by other simulation work [START_REF] Morrow | Theory of negative corona in oxygen[END_REF][START_REF] Gupta | Theory of step on leading edge of negative corona current pulse[END_REF]. The rise time of the first pulse is around 10 ns and the pulse falls to its 50% peak value in approximately 30 ns. These values have also been reported in [START_REF] Akishev | Numerical simulations of Trichel pulses in a negative corona in air[END_REF]. It can be seen that the two current calculation methods give very similar results. Volume integration takes longer to calculate as compared to surface integration. In the following sections, the method of surface integration was used to compute the results for that reason.

Different stages (I, II, III and IV) during the formation of the first Trichel pulse are marked in figure 4(a) and are discussed below. To support the analyses, the axial distribution of electron, positive and negative ion and space charge densities are plotted in figure 5. In addition, 2D snapshots of the distribution of charged particle density are shown in figure 6 at critical times. The variation of the electric field over time is shown in figure 7 along the symmetry axis and in figure 8 as a 2D distribution.

Electron multiplication (from I to II: 0-300 ns)

During this phase the seed electrons move away from the cathode under the influence of the electric field. As they travel, electrons gain high velocity and create additional electrons and positive ions through impact ionisation. The highest concentration of positive ions and electrons are found near the needle (figure 5 -300 ns curve) where the electric field is strong. As the electrons travel further away from the needle, attachment process becomes dominant in the low electric field region. This leads to a lower concentration of electrons and an increase of negative ion density over time in this region. An overall positive space charge occupies the volume near the needle (0.35 mm) while a slightly negative space charge fills the rest of the gap, figure 5(d). Since the created positive ions strike the surface of the needle, secondary electrons are released according to equation 6. This cloud of electrons then repeats the same cycle as the previous one hence the charged particle densities keep increasing. Evidently, it is this secondary process that sustains the negative discharge. The typical 2D charge density snapshot at 300 ns in figure 6(a) clearly demonstrates that positive ions are sandwiched between electrons and the cathode. The initial applied electric field is not heavily distorted during this phase until point II, figures 7 and 8, due to the low order of charge density (~10 17 m -3 ). Plasma sheath formation (from II to III: 300-348 ns) At this phase, the formation of positive ions and electrons has reached the critical concentration (order of 10 20 m -3 ) where the space charge field strongly distorts the initial applied field. The formation of a dense plasma near the cathode increases the field on the cathode while decreases the field within the plasma (figures 7 and 8 at 330 ns and 348 ns).

The movement of charged particles during this phase can be clearly observed in figure 5.

Until 348 ns, the peaks of both positive ion and electron densities increase and move towards the cathode. This process occurs as a result of the electric field enhancement at the needle by the plasma and has been experimentally identified as the motion of an ionising wave towards the cathode [START_REF] Cernak | Experimental confirmation of positivestreamer-like mechanism for negative corona current pulse rise[END_REF]. The negative ion density in this phase increases accordingly however it is still one order less than the positive counterpart. The axial space charge density plot, figure 5(d), clearly demonstrates an abrupt rise in concentration and fast movement of the charge particles during this phase. The positive space charge region is now reduced to 0.05 mm away from the needle followed by an increase of the negative charge concentration up to a distance of 0.3 mm. The dense plasma formation can also be observed in the 2D charge density distribution, figure 6(b). At 348 ns the peak concentration of positive ions moves closest to the cathode (approximately 10 µm away) and reaches its maximum density (2.6 × 10 20 m -3 ). This time instance also corresponds to the peak of the current waveform (point III) where the volume of positive ion concentration is reduced to its minimum. This creates a thin sheath on the cathode. High gradients of concentrations require application of a very fine mesh near the needle. It is also interesting to note that the peaks of both positive ion and electron densities are not on the symmetry axis but situated 50 µm away from the axis, and dense plasma actually forms a charged torus ring rather than a disk. At 348 ns, the electric field at the needle is approximately 6 times higher than the initial applied field whereas the field within the plasma reaches almost zero (figure 7). The ionisation region (where ionisation dominates attachment) is hence greatly reduced.

Decay process (from III to IV: 348-500 ns)

It occurs as soon as the plasma sheath formation is complete. The ionisation region has been reduced to its minimum and it is unable to sustain further discharge. The positive ion and electron densities start to drop while the negative ion density starts increasing due to attachment (figure 5, curves 450 and 500 ns). At 500 ns, the negative ion density is even higher than the electron density and of the same order as the positive charge density, figure 5(c). The space charge plot shows a decrease in magnitude of charge densities in this phase, figure 5(d). The volume occupied by positive space charge, having shrunk in the previous phase, now expands out to 0.1 mm axially. During this phase the electric field tends to return to the initial field distribution (curve 500 ns, figures 7 and 8). The next Trichel pulse occurs after the space charge particles have been gradually cleared away from the cathode by the electric field.

Influence of the ion-impact secondary emission coefficient

In the negative corona discharge model, the ion-impact secondary emission coefficient is the only parameter that is tuneable assuming the swarm parameters for air is well established. Effects of γ on the formation of Trichel pulses can be studied by varying the coefficient while keeping all other parameters unchanged as specified in the model. The first current pulse magnitude and the frequency of occurrence were calculated and plotted as a function of γ values (figure 9). As can be observed, the first Trichel pulse peak value gradually increases from 1.86 to 1.98 mA as γ increases from 0.002 to 0.006. Similarly, the frequency increases approximately linearly from 70 kHz to 195 kHz. The averaged charge per pulse has also been calculated for each γ value and it was confirmed that the charge also increases linearly from 0.09 to 0.12 nC over this range of γ. These effects can be explained by noting that an increase of γ leads to a more intensive electron production and more frequent discharge activities. Current (A)

1.84x10 -3 1.86x10 -3 1.88x10 -3 1.90x10 -3 1.92x10 -3 1.94x10 -3 1.96x10 -3 1.98x10 -3 2.00x10 -3
Frequency Current magnitude

Influence of the applied voltage

As a validation test for the model, the influence of the applied voltage on the pulse frequency was investigated. According to experimental studies [START_REF] Lama | Systematic study of the electrical characteristics of the "Trichel" current pulses from negative needle-to-plane coronas[END_REF], the relationship between the applied voltage and pulse frequency can be expressed as

2 0 1 ) ( 1 rS V V V K T F - = = ( 10
)
where F is the frequency (kHz), T the time period between pulses (ms), K 1 a constant depending on the electrode conditions and the surrounding environment, V the applied voltage (kV), V 0 = -2.3 (kV) the threshold voltage, r the radius of curvature (mm) and S the gap spacing (mm). From this formula, the plot of the frequency F against V(V-V 0 ) should demonstrate a straight line with a gradient of K 1 /(rS 2 ).

The model was run for a range of applied voltage (from -5.4 to -5.7 kV) whilst keeping all other parameters constant (γ = 0.004). Figure 10 shows the dependence of the predicted frequency F on V(V-V 0 ). The plot clearly shows the linear relationship with a gradient of 38 kHz kV -2 . The increase of the pulse frequency at higher voltages can be understood as follows: as the voltage increases, the electric field is stronger giving rise to a more intensive ionisation and a larger drift velocity. As a consequence the dense plasma forms at shorter times, charged particles are cleared away quicker and the next discharge pulse occurs sooner.

V(V-V 0 ) (kV 2 )

Experimental results

Using the experimental setup described in section 2, the Trichel pulse results obtained at -5.5 kV applied voltage and 3.3 mm gap spacing were recorded and shown in figure 11.

The Trichel pulses captured have a frequency of around 144 kHz which is very close to the simulation results obtained with a similar electrode geometry and the secondary emission coefficient set at 0.004. Details of the first experimental current pulse are superimposed onto the simulation result and are shown in figure 11(b). As can be seen, the pulse magnitude is around 3.5 mA which is higher than the simulation magnitude. Although the rising front of the simulation current waveform is slightly longer than that of the experimental curve, the model seems to predict the right decaying tail. The observed discrepancy is probably due to the influence of the power supply circuitry on the current which is not included into the model. The stochastic appearance of seed electrons in the neighbourhood of the needle and the exclusion of other less important secondary processes may also be accountable. The artificial diffusion used in the model may influence the current magnitude but hardly affect the decaying tail. The injected charge results calculated by integrating the current waveforms are shown in figure 11(c). An excellent agreement between the experimental and simulation injected charge is met. Each current pulse injects approximately 0.11 nC into the system. In order to verify the empirical relation in equation 9, the experiment was repeated with a range of applied voltage (from -5.4 to -5.7 kV) and the corresponding frequency was recorded. Figure 10 shows the observed frequency F against V(V-V 0 ). As can be seen, the experiment confirmed the linear relationship between the two parameters with a gradient of approximately 32 kHz kV -2 and a good agreement between the experimental and simulation results can be observed with a maximum prediction error of 15%.

5.

Negative surface discharge Having had the negative corona discharge model validated by comparison with the published data and experiment, the model can be used to solve the problem when dielectric barriers are present. Previous experimental observations of a negative surface discharge suggest that the discharge produces a diffuse and circular/ring distribution of charge [START_REF] Manabe | Formation mechanism of surface corona on dielectric plates under negative impulse voltage in atmospheirc air[END_REF][START_REF] Celestin | Patterns of plasma filaments propagating on a dielectric surface[END_REF]. The earlier model proposed in [START_REF] Manabe | Formation mechanism of surface corona on dielectric plates under negative impulse voltage in atmospheirc air[END_REF] is based on calculations along the strongest field line and can only predict the approximate size of the charge pattern. By modifying the full 2D model presented in section 3, a dynamic simulation of charge accumulation at the dielectric surface can be achieved.

Boundary and initial conditions

In this simulation, the insulation system has been introduced into the gas gap with the details given in figure 3(b). The simulation domain is similar to the corona discharge case with different boundary conditions at dielectric surfaces -boundaries 5-7 in figure 3(b). Most of the initial conditions were kept unchanged from the corona discharge model except the applied voltage which was set to -2.5 kV. Poisson's equation was solved over all 3 subdomains (air gap, BSO and glass) whilst the continuity equations were solved in the gas subdomain only. The same secondary emission coefficient as previously (γ = 0.004) was used in this model. The upper boundary of the BSO layer was set to free convective flux for electrons and negative ions to let them go through and the fluxes of these charged particles are used to calculate surface charge density at this surface (B5). The flux of positive ions at this boundary, in contrast, was set to zero as it is believed that no positive ions from the dielectric can be released into the gas. The new set of boundary conditions in this model is summarised in table 2. The same as in Table 1 The same as in Table 1 The same as in Table 1 The same as in Table 1 except Potential V = -2500 (V) for boundary 2 Boundary 5

Convective flux

( ) 0 = ∇ - ⋅ e e N D n Insulation/ Symmetry 0 ) ( = ∇ - ⋅ p p p p N D N W n Convective flux ( ) 0 = ∇ - ⋅ n n N D n Surface charge ( ) s D D n ρ = - ⋅ 2 1 Boundary 6 - - - Continuity ( ) 0 2 1 = - ⋅ D D n Boundary 7 - - - Ground V=0
The field in the gas gap is enhanced by the presence of the insulators and a much lower voltage in comparison with the corona discharge model is needed to be applied at the needle electrode. During the model implementation it has been found that a much finer mesh than that used in the corona discharge case is required in order to obtain a stable algorithm. The mesh size in the region close to cathode was hence reduced to 0.5 µm.

Simulation results

With such a gap spacing and electric field, the dielectric barrier discharge undergoes similar discharge events as in the corona discharge case. The significant difference is that only one discharge current pulse forms in this case. The predicted current waveform is shown in figure 12. The pulse is characterised by a fast rising front and a slower decaying tail. On the decaying curve, there exist some oscillations. The current magnitude is just above 1.7 mA, the initial pulse width is around 100 ns and the tail lasts approximately 1 µs. 15 respectively show the 2D charged particle density, electric field and surface charge density distributions at selected times. Similar to the previous analysis in section 4, the Electron Multiplication phase (0-50 ns) is characterised by the low concentration of positive ions in between the needle electrode and a lower concentration of electrons, figure 13(a). The initial field distribution is not distorted during this phase. Figure 14(a) clearly demonstrates the undisturbed equipotential lines in the gas gap. The surface charge accumulation during this phase is diminutive (figure 15, curve 50 ns). This phase is then followed by the Plasma Sheath Formation period (50-62 ns) during which the positive and electron densities increase in magnitude but decrease in size, figure 13(b). Both charged species move towards the cathode as a consequence of the field distortion, figure 14(b). The field distribution plot indicates an increase in magnitude of the field near the needle tip and an increase of the field in the bulk of the insulators due to the surface charge. At the instant the current waveform reaches its peak (62 ns), the surface charge density has a peak value of 0.05 nCmm -2 , figure 15.

The Decay process (62-2000 ns) starts after the current has reached its peak. During this phase, positive ion and electron densities decrease but negative ion density increases in magnitude. These negatively charged particles are then swept away towards the dielectric surface due to the electric field. At the beginning of the Decay phase, the increase rate of surface charge density is the largest (figure 15, curves 62-500 ns). This is a result of the movement of the dense cloud of the negative ions, which has a high density just after the Plasma Sheath Formation stage. At the end of the Decay stage, the surface charge density seems to reach its saturation level. The increase rate is much smaller than initially and after 2 µs the charge density acquires its peak at around 0.37 nCmm -2 . During the Decay process, the surface charge not only increases in magnitude but also expands radially. This can be explained by the effects of the surface charge field. Initially the charge deposits inside a circle of ~0.22 mm radius, figure 15 (curve 50 ns). As the dielectric surface gradually accumulates the negative charge and the charge density has reached a certain level, the surface charge field tends to oppose the Poisson's field and pushes the incoming ions away from the symmetry axis. At the end of the current pulse, the charge accumulates over a ring of ~0.7 mm radius. The drop of charge density at the centre is evident since the beginning of the Plasma Sheath Formation phase. The feature is linked with the dense plasma position being slightly offset from the symmetry axis, figure 13(b). It is worth noting here that this argument would not be applicable to the positive discharge as in that case electrons are produced at the streamer head by photoionisation.

Effort has been made to investigate the origin of such a ring-shaped density distribution since in the case of positive streamer with photoionisation as the driving mechanism these off-axis structures were attributed to numerical instabilities [START_REF] Pancheshnyi | Comment on "The role of photoionization in positive streamer dynamics[END_REF]. This effect was believed to originate from the incorrect representation of 3D phenomena in 2D cylindrical symmetry. To verify the proposed algorithm, space mesh refinement has been performed and smaller time steps have also been examined. The mesh refinement indicated the convergence of the algorithm. In addition, the artificial diffusion coefficient was varied significantly but the results hardly changed. Also a uniformly distributed initial distribution, instead of the Gaussian condition (equation 9) for positive ions and electrons was used and the off-axis behaviour was still observed. Effects of boundary conditions have also been analysed. The boundary condition (equation 6) was artificially modified by including the tangential component of ion outlet velocity. It was shown that the algorithm is not sensitive to the boundary condition implementation. In addition, the negative corona model has been modified to negative plane-plane geometry (infinite radius of curvature needle) to analyse effects of geometry on the observed off-axis behaviour. The off-axis effect was not present in this plane-plane test. As a result, it is quite possible that the observed behaviour may be a genuine phenomenon and it should be subjected to further discussions in the research community.

At the end of the current pulse (2 µs -figure 14(c)), the field lines in the air gap tend to return to the Laplacian distribution while the bulk field is still heavily distorted by surface charge. This field distortion appears in the region where the electrons form during the Electron Multiplication phase and this distortion prevents the development of subsequent Trichel pulses. The next pulse can be generated only after the surface charge has dissipated from the BSO. The Pockels experiment suggests that, in general, it takes 1-1.5 seconds for surface charge to decay.

Influence of the ion-impact secondary emission coefficient

In the case of the negative corona discharge (section 4), the Trichel pulse frequency and the first pulse magnitude linearly increase with the ion-impact secondary emission coefficient. As there is only one discharge pulse for the barrier discharge, this part focuses only on the effects of γ coefficient on the pulse magnitude, the surface charge density distribution and the total surface charge (integration of charge density over the dielectric surface). In comparison with the corona discharge case, in this part a wider range of secondary emission coefficients is considered. Figure 16 shows the changes in the peak current magnitude and total surface charge as γ varies from 0.002 to 0.012. It is evident that the current magnitude increases linearly with γ and the variations in the current magnitude are much larger than those calculated for the corona discharge (figure 9). The total surface charge also increases but seems to reach a saturation level as γ reaches the high values. The charge density distribution at the end of the discharge current pulse calculated for these values of γ is shown in figure 17.

Generally, an increase in γ value results in a rise in charge density and a small radial expansion. As γ increases, more charged particles are produced near the cathode and these particles will move towards the dielectric surface at the end of the current pulse. The charge density peak seems to move nearer to the symmetry axis as γ increases. At high values of γ, the increases in the charge density seem to drop due to the field suppression created by the surface charge. Excessive injection of electrons from the needle reduces the total field but does not result in additional charge accumulation above a particular threshold.

Experimental results

The surface charge density distribution generated from a negative 2.5 kV square voltage waveform was measured using the Pockels technique and is shown in figure 18. The discharge current waveform, figure 18(a), demonstrates a current pulse of 35 mA magnitude, initial pulse width of around 10 ns and the tail of approximately 70 ns duration. Although these numbers are different from the simulation results, the injected charge seems to be of the same order. The measurement indicates an amount of 0.4 nC charge injected into the system compared with 0.21 nC from the simulation. The experiment suggests 20 times shorter duration but also 20 times stronger current magnitude than predicted. The long tail in the predicted current waveform may stem from the artificial numerical diffusion (section 3.4) but there is another possible source of uncertainty. The electric field is calculated by quasi-static Poisson's equation and effects of the displacement current are neglected. Since the duration of the real discharge is quite short and the BSO has a high permittivity the quasi-static approximation could bring significant errors.

The measured two-dimensional charge density distribution and a typical radial charge density are shown in figures 18 (b) and (c). The charge is distributed over the diameter of around 2.1 mm with the peak charge density of 0.6 nCmm -2 . A drop in the charge density distribution at the symmetry axis position is noticeable. The predicted radius of the injected charge region is 30% less than that measured experimentally. The agreement is better for higher γ values (~0.01). Taking into account the uncertainties in γ and differences in the experimental and modelled needle shapes it can be concluded that the set of equations 1-4 could be effectively used to analyse dielectric barrier discharges. • At the Electron Multiplication stage of the discharge the peaks of both positive ion and electron densities are not on the symmetry axis but situated at some distance away from the axis and the dense plasma actually forms a charged torus ring rather than a charged ellipsoid. Possible effects of numerical instabilities and implementation of boundary conditions were investigated in detail. It was concluded that the algorithm is stable and converges as the mesh size and time steps are reduced.

It was suggested that the observed behaviour could be a result of imposed 2D axial symmetry which might not accurately describe the whole problem and full 3D modelling is needed to achieve perfect results. Nevertheless even 2D studies allow describing key features of the dielectric barrier discharges observed experimentally.

• It results in ring-like charge distribution on the surface of the dielectric barrier with a noticeable drop in the charge density at the symmetry axis position. Although such behaviour looks artificial, the charge distribution measured by the Pockels technique supports arguments for the observed off-axis effect. The problem should be subjected to further discussions in the research community. • Subsequent pulses in the Trichel sequence appear when the plasma particles are cleared away from the air gap and the field returns to its initial distribution. • In the case of the dielectric barrier only one discharge current pulse is produced instead of a pulse sequence as for the corona discharge case. The accumulated surface charge distorts the electric field, depresses electrons multiplication and prevents the formation of subsequent pulses. • It was found that the model predictions and experimental results agree well when the ion-impact secondary emission coefficient γ = 0.004 is used. • The magnitude of the current pulse depends linearly on γ but at high values of γ the total surface charge reaches a saturation level.

Figure 1 .

 1 Figure 1. Schematic diagram of the negative corona discharge experiment.

Figure 2 .

 2 Figure 2. Schematic diagram of the Pockels experiment.

3 .

 3 (a) Negative corona discharge (b) Negative dielectric barrier discharge Figure Simulation geometry of the negative discharge models.

  field vectors in the plasma over the surface and inside the cathode respectively,

Figure 4 .

 4 Figure 4. Predicted discharge current waveform of (a) the first Trichel pulse and (b) Trichel pulse sequence.

5 .

 5 (a) Electron density distribution along the symmetry axis (b) Positive ion density distribution along the symmetry axis (c) Negative ion density distribution along the symmetry axis (d) Space charge density distribution along the symmetry axis Figure Predicted distributions of charged particle densities along the symmetry axis.

  (a) 300 ns (II) (b) 348 ns (III) (c) 500 ns (IV) Figure 6. Predicted 2D distribution of charged particles densities at critical times. Surface plots are in log scale and contour plots are in linear scale.

Figure 7 .Figure 8 .

 78 Figure 7. Predicted electric field variations with time along the symmetry axis.

Figure 9 .

 9 Figure 9. Trichel pulse characteristics predicted for different secondary emission coefficients.

Figure 10 .

 10 Figure 10. Measured and predicted frequency F as a function of V(V-V 0 ).

11 .

 11 (a) Experimental Trichel pulse sequence (b) Details of the first current pulse (c) Injected charge waveforms Figure Trichel pulse results obtained from the experiment and the negative corona discharge model.

Figure 12 .

 12 Figure 12. Predicted discharge current from the negative surface discharge model.

Figures 13, 14

 14 Figures 13,[START_REF] Georghiou | The theory of short-gap breakdown of needle point-plane gaps in air using finite-difference and finite-element methods[END_REF] and 15 respectively show the 2D charged particle density, electric field and surface charge density distributions at selected times. Similar to the previous analysis in section[START_REF] Stoffels | Plasma needle: a non-destructive atmospheric plasma source for fine surface treatment of (bio)materials[END_REF], the Electron Multiplication phase (0-50 ns) is characterised by the low concentration of positive ions in between the needle electrode and a lower concentration of electrons, figure13(a). The initial field distribution is not distorted during this phase. Figure14(a) clearly demonstrates the undisturbed equipotential lines in the gas gap. The surface charge accumulation during this phase is diminutive (figure15, curve 50 ns). This phase is then followed by the Plasma Sheath Formation period (50-62 ns) during which the positive and electron densities increase in magnitude but decrease in size, figure13(b). Both charged species move towards the cathode as a consequence of the field distortion, figure14(b). The field distribution plot indicates an increase in magnitude of the field near the needle tip and an increase of the field in the bulk of the insulators due to the surface charge. At the instant the current waveform reaches its peak (62 ns), the surface charge density has a peak value of 0.05 nCmm -2 , figure15.The Decay process (62-2000 ns) starts after the current has reached its peak. During this phase, positive ion and electron densities decrease but negative ion density increases in magnitude. These negatively charged particles are then swept away towards the dielectric surface due to the electric field. At the beginning of the Decay phase, the increase rate of surface charge density is the largest (figure15, curves 62-500 ns). This is a result of the movement of the dense cloud of the negative ions, which has a high density just after the Plasma Sheath Formation stage. At the end of the Decay stage, the surface charge density seems to reach its saturation level. The increase rate is much smaller than initially and after 2 µs the charge density acquires its peak at around 0.37 nCmm -2 . During the Decay process, the surface charge not only increases in magnitude but also expands radially. This can be explained by the effects of the surface charge field. Initially the charge deposits inside a circle of ~0.22 mm radius, figure15(curve 50 ns). As the dielectric surface gradually accumulates the negative charge and the charge density has reached a certain level, the surface charge field tends to oppose the Poisson's field and pushes the incoming ions away from the symmetry

  Figures 13,[START_REF] Georghiou | The theory of short-gap breakdown of needle point-plane gaps in air using finite-difference and finite-element methods[END_REF] and 15 respectively show the 2D charged particle density, electric field and surface charge density distributions at selected times. Similar to the previous analysis in section[START_REF] Stoffels | Plasma needle: a non-destructive atmospheric plasma source for fine surface treatment of (bio)materials[END_REF], the Electron Multiplication phase (0-50 ns) is characterised by the low concentration of positive ions in between the needle electrode and a lower concentration of electrons, figure13(a). The initial field distribution is not distorted during this phase. Figure14(a) clearly demonstrates the undisturbed equipotential lines in the gas gap. The surface charge accumulation during this phase is diminutive (figure15, curve 50 ns). This phase is then followed by the Plasma Sheath Formation period (50-62 ns) during which the positive and electron densities increase in magnitude but decrease in size, figure13(b). Both charged species move towards the cathode as a consequence of the field distortion, figure14(b). The field distribution plot indicates an increase in magnitude of the field near the needle tip and an increase of the field in the bulk of the insulators due to the surface charge. At the instant the current waveform reaches its peak (62 ns), the surface charge density has a peak value of 0.05 nCmm -2 , figure15.The Decay process (62-2000 ns) starts after the current has reached its peak. During this phase, positive ion and electron densities decrease but negative ion density increases in magnitude. These negatively charged particles are then swept away towards the dielectric surface due to the electric field. At the beginning of the Decay phase, the increase rate of surface charge density is the largest (figure15, curves 62-500 ns). This is a result of the movement of the dense cloud of the negative ions, which has a high density just after the Plasma Sheath Formation stage. At the end of the Decay stage, the surface charge density seems to reach its saturation level. The increase rate is much smaller than initially and after 2 µs the charge density acquires its peak at around 0.37 nCmm -2 . During the Decay process, the surface charge not only increases in magnitude but also expands radially. This can be explained by the effects of the surface charge field. Initially the charge deposits inside a circle of ~0.22 mm radius, figure15(curve 50 ns). As the dielectric surface gradually accumulates the negative charge and the charge density has reached a certain level, the surface charge field tends to oppose the Poisson's field and pushes the incoming ions away from the symmetry

Figure 13 .

 13 Predicted 2D distribution of charged particles densities at critical times for the dielectric barrier discharge. Surface plots are in log scale and contour plots are in linear scale. (a) 0 ns (b) 62 ns (c) 2000 ns

Figure 14 .

 14 Figure 14. Predicted 2D electric field distribution at critical times for the dielectric barrier discharge.

Figure 15 .

 15 Figure 15. Predicted surface charge density distribution over time.

Figure 16 .

 16 Figure 16. Predicted discharge current pulse magnitude and total surface charge against γ.

Figure 17 .

 17 Figure 17. Predicted surface charge density obtained from different secondary emission coefficients.

  (a) Current and injected charge waveforms (b) 2D surface charge density distribution (c) Radial charge density distribution Figure 18. Experimental negative surface discharge results.

  discharges in air without and with dielectrics were analysed by 2D axisymmetric finite element models. The model input parameters were verified by experiments. It was shown that the model successfully describes both the localised nature of the discharges and their global characteristics. The observations can be summarised as follows:

Table 1 .

 1 Boundary conditions for the negative corona discharge model.

	Application	Convection and	Convection and
	mode	diffusion	diffusion
		N e	N p

Table 2 .

 2 Boundary conditions for the negative surface discharge model.

	Application	Convection and	Convection and	Convection and	Electrostatics
	mode	diffusion	diffusion	diffusion
		N e	N p	N n	V
	Boundaries
	1,2,3&4

7.


Soroush Nazarpour

MIND-IN Albert Cirera

Variation of adhesive force in the interface of Pd and SrTiO3 as a consequence of residual stresses

Keywords: Oscillating Indentation, Nanoscratch Test, Friction Force Microscopy, dislocation theory

Initially, Pd thin films were deposited over a hard substrate using electron beam physical vapor deposition. The growth and the surface roughness of the films were analyzed and their effects upon conventional indentation test were discussed.

Afterwards, an experimental method is described which can measure the critical fracture force in thin films using oscillating indentation. Initially, repetitive contacts at a single point with the purpose of identifying the fracture time provide the fracture force versus fracture time plot. Non-linear curve fitting of the data reveals the theoretical fracture force by a single indentation, which is called critical fracture force. Arguments are put forth to show the relation between piling up height and applied force. Discrepancies were observed in the plot of the ratio between total indentation depth and pilling up height versus applied force when higher loads than critical fracture force were applied.

Discrepancies appear as a result of indenting the substrate. Nanoscratch test facilitated the possibility to measure adhesion strength and adhesion energy of the films considering measured critical fracture force as maximum applied force. The relation between residual compressive stresses, adhesion strength, plastic deformation, and pilling up area was discussed using dislocation theories. Indentation with high applied load leaves behind large plastic deformation and reduces the accuracy and reliability of the test results. Hence, lower loads (in the order of nN) were applied using Atomic Force Microscopy in the friction mode. Pulling off force was mapped in each thickness Confidential: not for distribution.

Introduction

The successful performance and reliability of metallic thin films is often limited by their mechanical integrity and their adhesion to the substrate. Therefore, many applications of thin films and coatings require knowledge of their mechanical properties. Generally, thin metal films on rigid substrates sustain high stresses during growth process. The importance of stress has not decreased in all the years of intense investigation and it is still one of the main reasons for malfunction of thin film applications. Severe damages such as cracking [1], buckling [2], lifting [3], curling [4], and peeling [5] has been reported as a consequence of residual stresses in thin films. In fact, residual stress develops in thin films when the spacing of two neighboring atoms deviates from local equilibrium spacing. Numerous reports could be found in the literature describing the mechanism behind stress generation in thin films [6][7][8][9][10]. It is well established that residual stresses affect upon adhesive properties of the film [START_REF] Mattox | Handbook of Physical Vapor Deposition (PVD) Processing[END_REF]. Hence these stresses should be considered while measuring thin film mechanical properties. One powerful method for measuring thin film mechanical properties is nanoindentation [START_REF] Marshall | [END_REF][13][14][15].

Marshal and Evans [START_REF] Marshall | [END_REF] introduced a model in which an indentation creates a plastic impression on the film surface. The residual strain in the plastic impression causes the film to buckle upwards in order to minimize the strain energy, leaving a circular crack at the interface beneath the indentation, which results in a circular blister on the film. In their model, if the accumulated compressive stress in the film plus the stress associated with indentation is greater than the critical stress required for buckling, then the film will buckle from the interface and a ring crack forms in the interface and expands until the driving force for crack extension drops below a critical value [16]. However, this method is highly impractical for very thin films, particularly ones that are soft or have a very strong interface with the substrate. In fact, the strain energy in the Marshal and Evans [START_REF] Marshall | [END_REF] model basically is too small to induce delamination. In order to solve this problem, various researchers [17][18][19][20] have deposited thick overlayers capable of accumulating large amount of strain energy over the interface. In our previous work [21], a new method was introduced based on the oscillating indentation to accurately measure the critical fracture force in Pd thin films. In this article, the critical fracture force has been described as a critical force that indenter fractures the film and reaches to the interface without indenting the substrate. This method indents generally at low indentation loads to prevent the plastic deformation while repetitive contact of the indenter with the surface using oscillating indentation modifies the final contact area.

Here, it is tried to apply this method to measure the adhesive properties of Pd thin film.

Nanoindentation with a scratch test module was applied to measure critical stress and adhesion energy. During the scratch test, the load was ramped from 0N until the critical fracture force in each thickness measured by oscillating nanoindentation. The scratch length was set as 1000µm at a moving velocity of 50µm/s. Moreover, in order to overcome against plastic deformation such as piling up and sinking in area, the measured critical fracture by oscillating nanoindentation was converted into the critical applied load in Atomic Force Microscopy (AFM) using geometrical ratio between Berkovich indenter and AFM tip.

Experimental Method

The materials used in this study were 101B grade (Crystal GmbH) both sides polished (100) oriented SrTiO 3 (STO) substrate and 99.99% pure palladium. Pd and STO have been selected for this study since Pd/STO shows the behavior of a soft film over hard substrate. Moreover, application of Pd as electrode material and STO as insulator is widespread in electronic devices such as planoparallel capacitors. delamination of electrode metallic layers results in malfunction of these devices which might be attributed to residual stress. This is the reason to study the effect of residual stress upon interface properties of Pd and STO. Electron beam evaporator was used for palladium evaporation. Graphite crucibles were used for resistive evaporation. The evaporation chamber was evacuated with mechanical and turbo molecular vacuum pumps until the base pressure of the chamber was 3.1×10 -7 torr [22][23]. Palladium was evaporated at 1.2 kV accelerating potential and 140 mA emission intensity. Three thickness of Pd film (50, 100, and 150 nm) were deposited over STO substrate at room temperature. The ex situ film thickness was measured by a Dektak 3030 profilometer (Veeco). The crystalline structure was determined by X-Ray Diffraction (XRD) analyses in -2 scan by a 4-circle X-ray diffractometer with Cu-K radiation (MRD PHILIPS). Topographic images were recorded with Digital Instruments Multimode-Veeco AFM in tapping mode ( = 300 kHz). The nanoindentation tests were performed by a MML NanoTest TM nanoindenter with a Berkovich indenter (face angle 65.27º), whose tip radius is 100 nm.

From each indentation cycle, hardness, elastic module, and Poisson's ratio values were determined from the shape of unloading curve by means of Oliver and Pharr method [24]. Average elastic module and hardness were calculated by repeating indenting measurement with various forces at different points of the sample. Scanning electron microscopy (SEM) (ESEM Quanta 200 FEI) has been utilized for measuring the contact area after nanoindentation. Afterwards, repetitive contacts at a single point were obtained by oscillating the specimen against indenter that is mounted on a freely swinging pendulum. The process has been described in [21]. Afterwards, stress was measured in -tilt mode from 0º to 180º (2 = 38º -42º) in order to plot the d-spacing versus Sin 2 . Afterwards, by considering the measured Poisson ratio and elastic module, residual stress was measured for each Pd film [25,26]. To reduce the plastic deformation (piling up and sinking in regions), AFM MFP 3D (Asylum) was used to indent Pd films applying very low loads. Measurements are performed in buffer medium PBS using cantilever SNL (Veeco) with nominal spring constant equal to 0.06 nN/nm and a monolithic diamond tip (ND-DYIRS-4) with around 100nm tip radius.

Cantilever spring constant was calculated again using the thermal noise method and 102.8pN/nm was obtained. Moreover, AFM force curves are detected by pushing the AFM tip onto a substrate and recording the deflection and the stroke of the cantilever.

Pulling off force was mapped in a square with a 10×10 grid in 5×5 µm 2 area for each thickness of Pd film when maximum vertical force did not exceed 15nN. The average pulling off force was studied in each thickness of Pd films.

Results and Discussion

Oscillating nanoindentation

Surface roughness is a very important issue while measuring the adhesive properties of a film using nanoindentation. Since the contact area is measured indirectly from the depth of penetration, the natural roughness of the real surfaces causes errors in determination of the area of contact between the indenter and the substrate [27][28].

Hence, as much film is smoother, the accuracy of the indentation results is higher.

Previously it has been shown [22] that thin Pd films transfer the morphology of the underneath layer or substrate toward the final surface. Hence, surface roughness of Pd films depends upon substrate roughness as well as deposition conditions. Therefore, both sides polished (100) oriented STO with the root-mean-square (RMS) roughness less than 3nm has been selected as substrate. Low roughness of STO results in negligible contribution of the substrate roughness in the surface roughness of Pd films. Furthermore, deposition condition thoroughly modifies the surface quality of the films.

It has been found experimentally that at vacuum level less than 5×10 -5 Torr, and low accelerating potential and emission intensity (low deposition rate) reduces the surface roughness. In fact low vacuum level avoids the collision of the ad-atoms with the atmosphere in the chamber and more homogeneous thin films are expected.

Accelerating potential and emission intensity defines the size of the impinging atoms.

Basically smaller ad-atoms provide lower surface roughness. After impinging the adatoms over the substrate, surface diffusion controls the movement of ad-atoms over surface. At low substrate temperature (room temperature), surface mobility of ad-atoms is limited. Therefore, process such as coalescence which increases the surface roughness does not or not completely occurs. During impinging and growth, residual stresses generate in the films. Dislocations, defects and stacking faults reside in the lattice as a consequence of these stresses. Normally, deposition temperature by enhancing the surface mobility of the ad-atoms avoids of these lattice imperfections. Thus it could be expected that deposited Pd films at room temperature exhibits smooth surface and contains small grain size which are surrounded by grain boundaries with high density of dislocations. increasing the thickness of the films, surface roughness would be increased; RMS roughness of the films did not exceed 30nm. Hence, possible errors that arise due to surface roughness of the films are negligible during measurement of the contact area using nanoindentation. Contact area is a vital parameter while measuring hardness and elastic module of the films. Normally mechanical properties such as hardness and elastic module are measured from conventional force-indentation depth curves, as it would appear for a loading and unloading cycle into an elastic-plastic material. In fact, the ways in which the plasticity affects the interpretation of elastic unloading data can be dealt with the shape of the perturbed surface in the elastic analysis. Tabor et al. [29][30] demonstrated that the shape of entire unloading curve and the total amount of recovered displacement can be accurately related to the elastic modulus and the size of impact compression for conical indenter. However, there would be some important parameters which we should be aware of, such as high temperature stability and isolation of vibration which are the essential environment requirements during conducting the low force in order to measure the mechanical properties of thinner films. Oliver and Pharr method [24]. This method is a commonly used method for obtaining hardness and elastic module of the films. However there are some shortcomings in this method. First, at very low loads, there is a considerable offset error in the loadindentation depth curve. Therefore, hardness and modulus tend to decrease at very low loads due to effect of the elastic-plastic deformations [31]. However, by applying higher loads, the substrate effects must be considered. STO hardness has been measured equal to 8.9 ± 1.1 GPa. While measuring the hardness, we assumed that STO hardness is 9 GPa. During hardness measurement of soft films over hard substrate, since the film is extremely soft compared to the substrate and it yields at much smaller loads in comparison with the substrate, the film itself accommodates all plastic deformation and only when the indenter is close enough to the interface, substrate starts to yield [32][33].

This reveals the importance of measuring the critical fracture force during indentation [21]. Another potential problem is that this method was deduced from a purely elastic contact model for monolithic material. During indentation, only sinking in is assumed occurring on material surface and therefore, depth of contact is always less than real depth. However, in the case of soft thin films, when piling up occurs, depth of contact is higher than real depth. Hence this method is not accurate for plastic material where material plastically piles up around the indenter. These problems are the origin of the error bars in fig 3. Furthermore, metal films are tolerant of defects in their grain boundaries, whose influence can be blunted by plastic deformation effects. Hence, variation of hardness versus thickness of metal films is small. In fact, by increasing the thickness, ad-atom surface mobility enhances and hence, the density of dislocation in grain boundaries decreases. This is the reason that at lower thickness, higher density of GB dislocations results in higher hardness. Normally, hardness versus thickness of metallic films decreases until reaching to the bulk hardness.

Going beyond the determination of elastic modulus and hardness, a new method is introduced to measure fracture force at the interface of the soft thin films and substrate using repetitive contacts at a single point by oscillating the specimen against indenter [21]. Oscillating indentation with different applied loads was performed on the deposited Pd films and indentation depth versus time was drawn (not shown). The results reveal the fatigue crack growth and an abrupt film delamination. The fatigue crack growth was detected since the diamond indenter tip begins to move away from the surface caused by expanding the crack volume as evidence by a decrease in the indentation depth. Fig. 4(a) represents the fracture force versus fracture time. It could be seen that at lower loads, fracture time increases. Theoretically, at higher indentation time, hardening occurs in the film which results in the asymptotic shape of the fracture force versus fracture time plot [21] [21] that the ratio between h t and h i is constant by increasing the applied load until a critical value. Indeed, this critical value is critical fracture force. Inconsistency could be detected at loads higher than critical one which is due to the substrate deformation irrespective of the tip rounding. Fig 5(c) points out this variation as well. Simply, this could be a way to test the accuracy of the measured critical fracture force using piling up heights.

Effect of residual stress on the interface adhesion strength

Basically, large percentage of the plastic deformation in soft thin films is associated with piling up region and this area expands at higher loads. Randal [34] showed that at low loads, little pile up is present, although at higher loads, pile up area increases and a substrate relaxation occurs when indentation depth goes beyond the film thickness. In fact, applied load generates stresses around indentation area which will be released by piling up phenomenon. Lee et al. [35] showed a dependency between pile up height and applied stress into the specimen. Since thin films normally undergo residual stresses, pile up height might be dependent to the residual stresses in the films. Fig 6 represents the residual stresses in the deposited thin films measured using X-ray diffractions. In fact, stress was measured from the slope of the d-spacing versus Sin 2 plot. Afterwards, by considering the measured Poisson´s ratio and elastic module using nanoindentation, residual stress was measured for each Pd film [36]. Underwood [37] showed the effect of residual stresses upon the shift in the shape of the pattern around indentation. The pattern around indentation consists of piling up and sinking in area. Therefore, any mechanical properties measured using nanoindentation, are dependent to the residual stresses in the film. Marshal and Evans [START_REF] Marshall | [END_REF] showed that residual stresses in the films modify the crack extension force depending on the buckling of the films. Therefore, the effect of the residual stress upon the critical stress for the delamination ( c ) and the fracture energy release rate (G c ) for interface delamination (adhesion energy) is studied.

Nanoindentation with a scratch test module was applied to measure critical stress and adhesion energy. The applied force for interface delamination, defined as the critical loads (P c ) and the critical scratch track width (d c ) were measured from the scratch morphology as shown in fig 7. Critical stress for the delamination ( c ) (adhesion strength) was measured using equation (1) as follow [38][39][40]:

Eq. (1)
where µ is the measured friction coefficient of indenter sliding given by nanoscratch tester about 0.025 and f is the Paisson´s ratio measured in each thickness using nanoindentation. Thereafter, adhesion energy was measured using equation (2).

Eq. (2)

where t is the thickness of the film and E f is the elastic module measured by nanoindentation. Fig 8 represents the variation of the adhesion strength and adhesion energy versus residual stress in the films. It can be seen that at higher compressive residual stress and higher thickness of Pd film, adhesion strength and adhesion energy increased. However, since the contribution of the residual stress mixed with the contribution of the thickness, separation of the two parameters is extremely difficult.

Huang et al. [41] showed that residual compressive stress enhanced the adhesive properties of SiN/glass. Hence, an increase in residual compressive stress and thickness enhances the adhesive properties. Basically, compressive residual stress generates during thin film deposition. During indentation, these stresses blunt the cracks´tip and suppress their propagation. This is why typical force-displacement curves of the films which undergo compressive stress shift toward lower displacement. In fact, compressive residual stress can aid the nucleation of the dislocations during deposition. Elastic mismatch stress in epitaxial films provides a driving force for threading dislocations that extend from free surface of the film to the interface between film and substrate [42][43][44]. The propagation of threading dislocations in the film layer results in misfit dislocation lines along the interface. This process relaxes the elastic mismatch strain and increases the total length of the dislocation [45] when the thickness of the film is higher than Matthews´critical thickness [42]. During indentation with the forces lower than critical fracture force, the density of the misfit dislocation increases which leads to an increase in the adhesion strength of the films. Indeed, indentation may lead to appearance of perfect Lomer-Cottrell dislocation which results in interlocking between misfit dislocations in the interface between film and substrate. Hence, adhesion strength increases due to the interface hardening. There would be two possible solutions avoiding these dislocations. First solution is already discussed by applying critical fracture force onto the films with the fast ramping load. This can avoid the interface dislocation due to the fact that fast indentation avoids the relaxation in the interface, particularly when thin film ruptures under the applying force. Second solution would be applying very low loads using small tips. This is nowadays feasible using Atomic Force Microscopy in friction mode.

Pulling off mapping by Atomic Force Microscopy

As mentioned before, using AFM provides the possibility to conduct force measurement with low applying forces. Before conducting force measurement, it is necessary to accurately determine cantilever spring constant. Various methods are described in the literature for calibration of AFM cantilevers [46][47][48][49][START_REF] Tortonese | Proc. SPIE: Micromachining and Imaging[END_REF]. In our experiments the so-called thermal noise method was used, as it was found to deliver most reliable results. Spring of the piling up region corresponds to plastic deformation around the indentation area. It could be seen that maximum height of pile up area decreased from 70 nm until 60 nm by reducing the applied force from 250nN to 150nN. However, further decreasing the applied force from 150 nN to 50nN almost suppressed the appearance of the pile up area. It coincides by the fact that still clear indentation area could be seen which reveals that the film underwent plastic deformation. The surface after indentation with 50nN (f) is smoother in comparison with the indentation at higher forces (d and e). Hence, in agreement with the variation of the pile up height induced by indentation, it could be seen that in low indentation depth and low applying loads, pile up height is negligible. It should be noted that by further decreasing the indentation load (<1nN) indentation area disappears due to the elastic behaviour of the film. Hence, measurement of the exact indentation area would not be possible. Moreover, applying indentation with low forces avoids of tip rounding. Thereafter, to investigate the adhesive properties of the film, pull off test was done over a 10×10 grid (5 µm ×5µm area) of each Pd film. Each grid has been selected equal to 0.5µm×0.5µm. It should be noted that it might be possible to decrease the grade size, but the attention should be made about plastic deformation around each indentation. Fig 9 (f) shows that plastic deformation are after indentation with 50nN is a sphere with around 100nm radius. This area on the surface is even smaller at lower applying loads (<5nN), although in the area far from the surface (deeper area) it is not possible to calculate it. It should be considered that the distance between indentation points should not be lower than the plastically deformed area.

Therefore, in order to enhance the accuracy, the distance between indentation points was selected about 0.5µm. Indeed, when the AFM is retraced, the adhesion force causes a kink in the force curve. Pulling off force varied from 4.23×10 -10 N until 6.96×10 -10 N. It should be noted that the indentation on the area around hillocks entirely represents high pulling off force. This might be in connection with the presence of hillocks. It has been previously showed that hillocks appear on the surface of the metallic film as a result of stress relaxation [START_REF] Nazarpour | [END_REF]. It could be concluded that the area around a hillock undergoes lower residual stress.

Therefore, the density of dislocation is lower than not-relaxed area. This in turn results in an area with lower hardness. When an area exhibits softer properties, while retracing the AFM tip, film sticks to the tip which results in a larger pile up area around indentation imprint. Therefore, pulling off force increases in this area.

Conclusion

In this study, Pd thin films were grown using electron beam physical vapor deposition. Surface roughness and crystal structure of the deposited film were studied and their effects upon accuracy of the force-displacement curves were discussed. Moreover, the mechanical properties and interface adhesion energy of 50nm, 100nm, and 150nm Pd thin film were measured using oscillating nanoindentation. Critical fracture force was introduced as a critical force to reach into the interface without impacting the substrate. It was measured in each thickness of Pd films and their connection with the pile up area was drawn. It was found that at higher applied forces than critical fracture force, the plot of h t /h i versus applied force represents irregular variations as a result of indenting the substrate. Nanoscratch test assists with the measurement of the adhesion strength and adhesion energy of Pd films. It has been found that compressive residual stress enhances the adhesive properties of the films. The compressive residual stress can aid the nucleation of the dislocations during deposition. The propagation of threading dislocations in the film layer results in misfit dislocation lines along the interface. This process relaxes the elastic mismatch strain and increases the total length of the dislocation in the interface. During indentation with the forces lower than critical fracture force, the density of the misfit dislocation increases which leads to an increase in the adhesion strength of the films. However, applying forces equal to critical fracture force lead to large plastic deformation and pile up area around the indentation area. This plastic deformation may reduce the accuracy of the indentation process. Therefore, Atomic Force Microscopy was used to apply low loads leave behind a small plastically deformed area. Pulling off force was mapped for each Pd film using applied forces not more than 15nN. It is found that the area around hillocks represents higher pulling off force. This is due to the fact that presence of hillocks releases the residual stresses of the area around it. Hence, AFM tip indents a softer area which results in higher pulling off during retracing. It could be concluded that during nanoindentation, applying forces equal to the critical fracture force can enhances the accuracy of the results. Furthermore, applying low loads using AFM perfectly reduces the plastic deformation around the indentation area. exhibits the pulling off map of (a). It can be seen that pulling off force varied from 4.23×10 -10 N until 6.96×10 -10 N.

  Fig 1 represents highly textured Pd film toward (111) crystal direction. Basically, face centered cubic thin film materials tend to growth toward (111) preferred orientation. Fig 2 shows the smooth surface of 50nm Pd thin film. Although by

Fig 3

 3 Fig 3 represents hardness and elastic module of the deposited Pd films calculated using

  . The crossover of the plots in fig 4(a) and Y-axis identifies critical fracture force. Theoretically, critical fracture force is a theoretical fracture force that reaches the indenter to the interface between film and substrate in the first stage of indentation. Fig. 4(b) represents the critical fracture force versus thickness of Pd thin film. As shown in fig. 4(b), critical fracture force increases at higher thickness of Pd thin film. Fig 5(a) shows SEM image of 150nm Pd film which is indented with measured critical fracture force (1.3mN load for 150nm Pd film (fig 4(b))). Pilling up area could be detected around the indentation area. Fig 5(b) reveals the height profile along the line which is drawn in fig 5(a). h t is the total depth of indentation and h i is piling up height. It has been shown previously

  constant has been measured equal to 102.8 pN/nm. Fig 9 shows the indentation area using AFM diamond tip over 50nm Pd films. Fig 9 (a), (b), and (c) correspond to the indentation with 250nN, 150nN, and 50nN, respectively. Fig 9 (d), (e), and (f) are the counter plots of the fig 9 (a), (b), and (c), respectively. It could be seen that pile up area around the indentation area became smaller at lower applying loads. The height and area

Fig 10 (

 10 Fig 10(a) shows the grids in an area in the 50nm Pd film. Crosses in fig 10(a) represent

  Fig 11 represents the average pulling off forces of each thickness of Pd films versus indentation load. Apparently, after a critical applied load, discrepancy appears in plot of the pulling off force versus load. This irregular variation of the points (showed by the squares in fig 11) might be related to the plastic deformation in the films. Basically, as much indentation load increases, sticking increases during retracing. Consequently, pile up area generates around the indentation and its height becomes larger as much the indentation load increases. In principal, other parameters such as tip rounding, porosity, and generation of the stacking faults during indentation might impact on the accuracy of the results. Ability in applying loads even in the range of pN is a significant advantage of using Atomic Force Microscopy. Further study is essential to clarify the main reason for observed discrepancies in fig 11.

Figure 1 .

 1 Figure 1. X-ray diffractions of 50nm, 100nm, and 150nm deposited Pd thin films over STO substrate at room temperature. Single (111) diffraction of Pd film shows that deposited films are highly textured toward (111) preferential orientation.

Figure 2 .

 2 Figure 2. AFM topographic 2D image (a) and 3D (b) image of the 50 nm Pd film represents a smooth surface with RMS roughness less than 10nm.

Figure 3 .

 3 Figure 3. Measured hardness and elastic module using nanoindentation versus thickness of Pd film

Figure 4 .

 4 Figure 4. Fracture force versus fracture time for different thickness of Pd thin film (a) and critical fracture force versus thickness of Pd thin film (b) measured using oscillating nanoindentation [21].

Figure 5 .

 5 Figure 5. SEM image of 150nm Pd film indented with 1.3mN (a) and height profile (b) along the line which is drawn over (a).

Figure 6 .

 6 Figure 6. Residual stress measurement using variation of d-spacing of (111) crystal plane. Large error bars are due to the 2 position of (111) crystal diffraction (around 40°). Basically, accuracy of the measurement enhances when the planes at higher 2 positions ( 180) were studied.

Figure 7 .

 7 Figure 7. SEM surface morphology along the scratch track on 150 nm Pd film facilitating the measurement of the critical scratch track width.

Figure 8 .

 8 Figure 8. Variation of the adhesion strength and adhesion energy versus residual stress in the films. It can be seen that at higher compressive residual stress and higher thickness of Pd film, adhesion strength and adhesion energy increases.

Figure 9 .

 9 Figure 9. AFM topographic image of the indentation areas using AFM diamond tip over 50nm Pd films. (a), (b) and (c) correspond to the indentation with 250nN, 150nN, and 50nN, and. (d), (e), and (f) are the counter plots of (a), (b), and (c), respectively.

Figure 10 .

 10 Figure 10. AFM topographic image of 5µm×5µm area of 50nm Pd film (a). This area is divided to 10×10 grids when the crosses in each grid represent the indentation area. (b)

Figure 11 .

 11 Figure 11. The average pulling off forces of each thickness of Pd films versus indentation load.

FigureFigure 11 (

 11 Figure 01 (Fig 1.tif)