Rui Qiu

Member, IEEE Wenjun Zhou

Pascal Brault

Member, IEEE Yu Zheng

Gen Li

Han Yang

Han Li

IEEE TRANSACTIONS ON DIELECTRICS AND ELECTRICAL INSULATION A Method to Localize the Axial Position of Partial Discharges in GIL by Gas Sensors

Keywords: insulated transmission lines, partial discharge localization, gas diffusion, temperature rise, convection, gas decomposition, gas sensors

This paper proposes a novel method to realize the partial discharge (PD) localizing function of gas sensors in GIL. Firstly, the diffusion characteristics of decomposition products were studied. The Fuller-Schettler-Giddings model was corrected to obtain decomposition products' diffusion coefficients in high-pressure SF6. The finite element method (FEM) calculated the temperature rise, convection, and diffusion of decomposition products in GIL. With the relationship between the concentration variation and positions, a PD defects location method was proposed. The results indicate that the decomposition products would diffuse rapidly in the radial direction with the convection. Gas diffusion at the top of the GIL was faster than that at the bottom. Hence, the gas sensors were suggested to be set at the top of GIL. The axial diffusion was much slower and mainly affected by the temperature and diffusion coefficient. For a long GIL unit, the diffusion of decomposition products could be simplified to one-dimensional. In diffusion, the delay of the product concentration is negatively correlated with the square of the distance. That quantitative relationship could be used for localizing PD by the concentration increase delay detected by the sensors on both sides of the defect. The localization error would decrease with the increase of diffusion coefficient (rate) and prolonged monitoring (diffusion) time. The accuracy of the localization would also increase when the defect is close to the midpoint of the sensors.

I. INTRODUCTION

AS-insulated transmission lines (GIL) have advantages such as high voltage levels and good mechanical flexibility [START_REF] Koch | Technical aspects of gas insulated transmission lines and application of new insulating gases[END_REF], [START_REF] Tateyama | Conceptual Design and Feasibility Study of Flexible Gas-Insulated Transmission Line Using CF3I Gas Mixture[END_REF], [START_REF] Liang | Electric field reconstruction inside gas insulated transmission line by induced charge tomography[END_REF]. However, partial discharge (PD) is the main problem affecting the stability of GIL [START_REF] Okubo | Partial discharge measurement in long distance SF6 gas insulated transmission line (GIL)[END_REF], [START_REF] Zang | A novel optical localization method for partial discharge source using ANFIS virtual sensors and simulation fingerprint in GIL[END_REF]. The long axial distance of GIL also makes defect localization essential and challenging. As one of the gas analysis technologies, the gas sensors could distinguish the PD defect by monitoring the variation of gas components [START_REF] Ding | Factors affecting PD detection in GIS using a carbon nanotube gas sensor[END_REF], [START_REF] Zhang | Pt-doped TiO2-based sensors for detecting SF6 Decomposition Components[END_REF]. However, same with other gas analysis methods (e.g., gas chromatography [START_REF] Zeng | Evaluating DC partial discharge with SF6 decomposition characteristics[END_REF], infrared absorption spectroscopy [START_REF] Gao | Quantitative detection of multicomponent SF6 decomposition products based on Fourier transform infrared spectroscopy combined with CARS-ELM algorithm[END_REF]), gas sensors still have a shortcoming in localizing PD. That also severely limits the use of gas sensors on the GIL. It is necessary to use other technologies (e.g., traveling wave [START_REF] Li | A revised model for calculating HV cable sheath current under short-circuit fault condition and its application for fault location-Part 1: The revised model[END_REF], optics [START_REF] Zang | A novel optical localization method for partial discharge source using ANFIS virtual sensors and simulation fingerprint in GIL[END_REF], ultra-high-frequency methods [START_REF] Hu | Design of a distributed UHF sensor array system for PD detection and location in substation[END_REF], and ultrasonic methods [START_REF] Avzayesh | A hybrid estimation-based technique for partial discharge localization[END_REF]) to localize the PD. However, some weak PDs are challenging to be monitored by above-mentioned localization methods but could be detected by gas analysis techniques. Therefore, if the PD localization function of gas sensors could be realized, the application and role of gas sensors in GIL insulation diagnosis could be improved.

According to the detection principle, gas sensors could be divided into photoacoustic spectroscopy sensors [START_REF] We | Detection of SF6 decomposition components under partial discharge by photoacoustic spectrometry and its temperature characteristic[END_REF], [START_REF] Cheng | Simultaneous detection of C₂H₂ and CO based on cantilever-enhanced photoacoustic spectroscopy[END_REF], nanotube sensors [START_REF] Zhang | Pt-doped TiO2-based sensors for detecting SF6 Decomposition Components[END_REF], [START_REF] Martin | Detection of SF6 decomposition products generated by DC corona discharge using a carbon nanotube gas sensor[END_REF], etc. Gas sensors were first applied to gas-insulated switchgear (GIS). Due to the compact structure of the GIS chamber, the defect location function of the gas sensor received little attention. However, in the long distance and large volume of GIL chamber, the gas components in different positions vary significantly. The position of the gas sensor arrangement becomes important. The gas sensors could also realize the PD localization by the difference in gas components difference in different positions. That requires the clarification of gas diffusion procession in GIL [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF]. Due to the difficulty in obtaining diffusion coefficients, the current research on the diffusion process of decomposition products in GIL was limited [START_REF] Chen | Numerical calculation and correlative factors analysis on temperature distribution of GIS bus bar based on coupled multi-physics methodology combined with multiple boundary conditions[END_REF], [START_REF] Wu | Contact temperature prediction in three-phase gas-insulated bus bars with the finite-element method[END_REF]. Ming et al. predicted the diffusion coefficients of eight decomposition products by molecular volume and the Fuller-Schettler-Giddings (FSG) model [START_REF] Liu | The computation of diffusion characteristics of decomposition gases in SF6 and SF6/N2 within gas insulated transmission lines[END_REF]. The FSG model has significant errors under high pressure. Therefore, the accuracy of the prediction results needs to be verified. Our team also tested the diffusion process of decomposition gases in GIL [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF]. The diffusion coefficients of SO2F2 and SO2 in SF6 were measured. However, there are diverse decomposition products in SF6. Enriching diffusion coefficients for different decomposition products under high pressure is still necessary.

In addition, the abovementioned research ignored the temperature and convection on gas diffusion in the GIL. In operation, heat would be generated by the current in the GIL conductor. The heat transmission in the gas would generate convection. Then, the temperature and convection would affect the molecular thermal motion and gas transport processes, respectively. Therefore, in order to accurately simulate the diffusion process in GIL, it is also important to study thermal convection. The related experiment and simulation research is adequate and referable. For example, Ham et al. analyzed the temperature rise into the GIL by heat balance equation, finite elements methods (FEM), and temperature rise test, respectively. They found that the fluid dynamics simulation results were closer to the experimental results [START_REF] Ham | Heat transfer in gas-insulated bus bars[END_REF]. The temperature rise test of 550 kV GIL was carried out by

G IEEE TRANSACTIONS ON DIELECTRICS AND ELECTRICAL INSULATION

Minaguchi [START_REF] Minaguchi | Heat transfer characteristics of gas-insulated transmission lines[END_REF]. The temperature rise calculation results could be consistent with the experimental data well. Those researches provide the crucial basis for simulating the diffusion characteristics of decomposition products in GIL with temperature rises and convections. That could also provide the relationship between the variation and position of the components, the basis for localizing PD defects by gas sensors.

In this paper, we tried to determine the diffusion coefficients of some decomposition gases through experiments. Then, the existing FSG model to predict diffusion coefficients could be modified for high-pressure conditions. Furthermore, finite element simulation could analyze the temperature rise, convection, and diffusion in GIL. A novel method for localizing the PD defects in GIL would be proposed based on the diffusion characters from the above results and models.

II. SETUP OF SIMULATION AND EXPERIMENT

As illustrated in Fig. 1, the diffusion process is promoted by two main factors: the concentration gradient and the flow. The diffusion coefficients determined the diffusion rate for different components driven by the same gradients. In addition to the properties of the components, the diffusion coefficient is affected by pressure and temperature. When the temperature rise happens in the gas, with the heat transfer, the convection would also be promoted. Hence, the temperature rises in gas affected the diffusion via flow and concentration gradient. According to the Fig. 1, the diffusion process of decomposition products in GIL could be analyzed via Fick's second law [START_REF] Greenkorn | Momentum, heat, and mass transfer fundamentals[END_REF]:

  2 d pd /SF6 d d C C D v C t        , ( 1
)
where t (h) is the time, Cd (ppm) is the concentration of decomposition products, v is the flow velocity (m/s), and Dpd/SF6 (m 2 /S) is the diffusion coefficient of decomposition products in SF6. The basis for solving (1) lied in obtaining Dpd/SF6 and v.

A. Method to obtain diffusion coefficient

The optimized Fuller-Schettler-Giddings (FSG) model was used to calculate and predict the diffusion coefficient Dd/SF6 [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF]

TT D K k K PP      , (2)
where Kpd/SF6 is the diffusion characteristic coefficient, Kpr/SF6 is the diffusion characteristic coefficient of the reference product, kr is the ratio of Kpd/SF6 to Kpr/SF6 in [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF], [START_REF] Liu | The computation of diffusion characteristics of decomposition gases in SF6 and SF6/N2 within gas insulated transmission lines[END_REF], T (K) is the temperature, and P (Pa) is the absolute gas pressure. T was related to the temperature distribution in GIL. Pressure assumed to be a constant in this paper, is determined by the GIL voltage level. In this paper, SO2F2 was chosen as the reference decomposition product [START_REF] Zeng | Evaluating DC partial discharge with SF6 decomposition characteristics[END_REF]. The diffusion experiments were introduced to determine the Kpr/SF6 of SO2F2. In Fig. 2, the experiment platform consisted of a diffusion chamber and a GC-9760B gas chromatography. The diffusion chamber was equipped with sampling points A and B. They were spaced 1 m apart, which could detect the concentration difference leaded by diffusion. The diffusion source located at 0.5 m from point A.

The chamber and diffusion source were inflated with equal pressure of SF6 and SO2F2, respectively. The pressure was 0.35 MPa, 0.40 MPa, and 0.45 MPa, the diffusion time was 108 h, and the sampling period was 12 h.

B. Temperature rise, convection, and diffusion simulation model in GIL radial direction

Since the axisymmetric structure of GIL, the temperature rises and convection in GIL are mainly in the radial direction. Hence, a two-dimensional radial model was introduced [START_REF] Chen | Numerical calculation and correlative factors analysis on temperature distribution of GIS bus bar based on coupled multi-physics methodology combined with multiple boundary conditions[END_REF], [START_REF] Wu | Contact temperature prediction in three-phase gas-insulated bus bars with the finite-element method[END_REF], [START_REF] Ham | Heat transfer in gas-insulated bus bars[END_REF]. The GIL has two heat sources: the conductor (load current heat) and the shell (induction current heat). The fluid heat transfers would induce temperature gradients and convection.

According to Maxwell's equations, the electromagnetic field in the GIL satisfies:

2 n n s 2 n A ( )A J (conductor) A 0 (shell) jT           (3)
where An (Wb/m) is the magnetic vector perpendicular to the GIS radial section, Js (A/m 2 ) is the current density inside the conductor, ω (rad/s) is the angular frequency, μ (H/m) is the permeability, χ(T) (S/m 2 ) is the material conductivity related to temperature T:

0 ( ) / (1 0.004( 273)) TT     , ( 4
) where χ0 is the material's electrical conductivity at 273 K. For 6063 aluminum alloy used for GIL, χ0 is 55 S/m 2 . The tangential component of the air edge magnetic vector is 0 Wb/m. The heat generation Q of the conductor and shell is

  2 sn ( ) / ( ) d S Q j T T S     JA , (5)
In GIL, three types of heat transfer work: conduction, radiation, and convection. The conduction mainly occurred in the conductor and shell. The radiation appeared on the surfaces of the conductor and shell. The convection primarily exists in SF6 and Air. The conduction process could be solved by [START_REF] Ding | Factors affecting PD detection in GIS using a carbon nanotube gas sensor[END_REF]. Where λs (W•m -1 •K -1 ) is the thermal conductivity of solid material, and Qv (W/m) is the volumetric heat source:

2 vs QT     . (6)
The conservations of mass, momentum, and energy in (7) controlled the heat transfer process by the convention.

     g T g g 2 g g v 0 (mass) ( ) ( ) 1.5 ( ) g (momentum)
(ene rgy)

P T T Q                                    v v v E v v vE v , (7)
where v (m/s) is the gas flow vector at the conductor and shell surfaces, v = 0 m/s; Γ (J•kg -1 •K -1 ) is the heat capacity of the gas material, ρg (kg/m 3 ) is the gas density, λg (W•m -1 •K -1 ) is the gas thermal conductivity, ϑ (Pa•S) is the viscosity. λg and ϑ could be calculated by Sutherland's law [START_REF] Greenkorn | Momentum, heat, and mass transfer fundamentals[END_REF]:

T Su T Su T Su T Su                       (8)
where λ0 and ϑ0 are gas thermal conductivity and viscosity at 273 K, and Su is the Sutherland temperature. The radiative heat transfer between the outer surface of the conductor and the inner surface of the shell could be solved by the Stefan-Boltzmann law:

  44 S ij i j T F T T       , ( 9
)
where σ is the Stefan-Boltzmann coefficient, 5.67×10 -8 W•m - 2 •K -4 ; Ti and Tj are the temperatures of units i and j, ε is the unit emissivity, and Fij is the angular coefficient between units i and j. Furthermore, the radiative heat transfer of the shell's outer surface to the air could be: . [START_REF] Li | A revised model for calculating HV cable sheath current under short-circuit fault condition and its application for fault location-Part 1: The revised model[END_REF] This paper set up a finite element model in Fig. 3 according to the abovementioned equations and a 550 kV GIL. The free triangular mesh was used with a maximum size of 255 mm. The model had a total of 20644 units. The parameters set are elaborated in Table I. Based on the temperature rise simulation results and (2), the DP/SF6 at different radial positions could be obtained. The radial diffusion process was simulated based on the convection results and the diffusion models in [START_REF] Koch | Technical aspects of gas insulated transmission lines and application of new insulating gases[END_REF]. The fluxes on the conductor's outer surface and the shell's inner surface are 0 ppm•m/s. The concentration at the defect Cds(t) was assumed to increase linearly with time.

C. Diffusion simulation model in GIL axial direction

As mentioned, temperature and convection differences are negligible in the GIL axial direction. Hence, the radial temperature rises and convection simulation results could be the constant conditions in the 3-dimensional axial diffusion model in Fig. 4. The length of the GIL model is about 100 m [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF]. The decomposition products flow in from one side and out to the other. The diffusion time is 100 h. A free tetrahedral grid dissection was used. The maximum length of the cells is 10 mm, and there are 2378487 cells in the model. 2), ( 8), and Fig. 6, the Kpr/SF6 could be calculated as 1.521×10 -3 . According to the kr in [START_REF] Martin | Detection of SF6 decomposition products generated by DC corona discharge using a carbon nanotube gas sensor[END_REF], the Kpd/SF6 of 8 specific decomposition products: HF, H2S, SO2, SOF2, SO2F2, S2OF10, CS2, and COS were predicated and listed in Table II. In addition, the Dpd/SF6 of SO2F2 and SO2 in [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF] is 1.113×10 -4 m 2 /s and 1.397×10 -4 m 2 /s, respectively. The prediction value in this paper was 1.119×10 -4 m 2 /s and 1.444×10 -4 m 2 /s. Their relative differences were 0.539% and 3.255%. The prediction values of diffusion coefficients in this paper are close to the experiment one.

III. RESULTS

A. Diffusion experiment results and diffusion coefficients

B. Radial temperature rise and convection

The calculation results were compared with the experimental data at 2.6 kA, 4.4 kA, and 5.6 kA, respectively. Table III demonstrates the comparison results. Their average relative error is 2.75 %. The accuracy of the temperature rise models in GIL is receptible. The simulation result of a 550 kV GIL temperature rise at 5.6

kA is visible in Fig. 7(a). The highest temperature was 86.7 ℃ at the conductor, and the lowest was 48.3 ℃ at the bottom of the enclosure. The temperature distribution in the GIL had a symmetry characteristic, and the gas temperature gradually decreased from top to bottom. Fig. 7(b) illustrates the distribution of convection. The flow velocity at the top of the conductor and shell was the largest. The flow velocity at the bottom of the chamber was low. That could be because the SF6 near the conductor heats up fast. The decreased density caused by heating would make this part of SF6 first overcome gravity to move upwards. When the heated SF6 arrives at the shell, it would slide along the inner surface. In addition, a larger temperature gradient would lead to more convection. Hence, the high-speed area was mainly concentrated on the top of GIL.

There was almost no convection in the low-temperature region at the bottom of GIL. without the temperature rise and convection at 0.85 h, respectively (case, the defect was at a point). The temperature rise and convection would manifestly promote decomposition products to diffuse uniformly. Based on the (1) and ( 2), it could be inferred that a higher temperature rise in GIL could lead to a faster radial diffusion. In addition, the difference in diffusion coefficients did not cause the different diffusion characteristics of decomposition products in Fig. 9(b). It suggests that convection is the primary factor in the radial diffusion, and the temperature rise is the secondary one. That could help to improve the synchronization of the radial monitoring of GIL decomposition products. Fig. 10 illustrates the SO2F2 concentrations on the GIL shell when the defect was at different positions after 2h diffusion. In this paper, the difference between the concentration detected by sensors and that at defects is called the error. The relative error with defects at a, b, α and β were small, only 9.79%, 5.52%, 7.01% and 17.99%, respectively. That may be due to because the strong convection at the top of GIL. Conversely, when the defects were at the bottom c and γ points, where convection was weak and temperature was low, the relative errors increased to 34.71 % and 34.21 %. In the case of defects at c and γ, the anisotropy problem in monitoring results was also noticeable. In summary, GIL temperature rise and convection could diffuse the decomposition products rapidly, simultaneously, and uniformly in the radial direction. Moreover, convection is the primary factor in the radial diffusion. At the top of GIL, with high temperature and strong convection, the difference results could be more accurate and faster. While the temperature rises and convection at the bottom of GIL are weak. The diffusion process is slow, resulting in a large difference between the concentration of decomposition products at the defect and that detected by the gas sensors. Therefore, it was suggested that the gas sensors are set at the top of the GIL.

D. Axial diffusion characteristics of SF6 decomposition products in GIL

The rapid radial diffusion and the large difference in radialaxial dimensions would allow the decomposition products to have diffused in uniform in the radial direction at the beginning of the axial diffusion [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF]. Therefore, this paper ignored the radial distribution of decomposition products at the beginning of axial diffusion. The cross-section at the inlet was a diffusion source (generation rate: 0.1 ppm/h). α, β, γ lines, and αe, βe, γe points were at the top, middle, and bottom of the shell. αe, βe, γe points were 5m away from the inlet. These lines and points were used to investigate concentration distribution and variation characteristics during the axial diffusion. Fig. 11(a) illustrates the concentration distribution in the first 5m on the α line after 100h of diffusion. The concentration decayed exponentially with the increase in the axial distance [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF]: a smaller Ddp/SF6 would lead to more decay. Fig. 11(b) shows an apparent diffusion delay in the concentration variation at αe. The delay was also negatively correlated with Ddp/SF6. The distributions of SO2F2 concentration on α, β, and γ lines are illustrated in Fig. 12. It suggested that, the effect of temperature rise and convection on axial diffusion could be neglected. That was because the convection vector had no axial component. Hence, Ddp/SF6 and temperature mainly controlled the axial diffusion. However, Ddp/SF6 was stable under high pressure. Even the temperature at the top was higher than that at the bottom, see Fig. 7(a), the Ddp/SF6 of SO2F2 at the top of GIL was only 9.98 × 10 -5 m 2 /s, which close to the bottom one of 8.25 × 10 -5 m 2 /s. Therefore, the concentration in the radial direction was almost the same everywhere during the axial diffusion. In addition, the average temperature rise in radial cross-section was used to calculate Ddp/SF6 and the axial diffusion model, i.e., the 3D model was simplified to 1D. Fig. 12 indicates that the difference between the 1D and 3D models was around 0.001 ppm. The axial diffusion process of the decomposition products in GIL could be approximated by a 1D IEEE TRANSACTIONS ON DIELECTRICS AND ELECTRICAL INSULATION model. That also proved that the conclusions and results in [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF] could be referenced. Based on the models in [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF], in a 100 m GIL unit, the decay and delay of the concentration could be:

C C C x t D t x D                   , (12)
where x (m) is the distance, CdD and CdM (ppm) are decomposition product concentrations at defects and gas sensors, t (h) is the time, and η is the decay ratio. The axial diffusion is mainly affected by the diffusion coefficient. The concentration decay and delay caused by axial diffusion are more significant. Therefore, axial diffusion would affect the monitoring results of gas sensors severely. However, that significant difference in concentration variation due to axial diffusion could help to do the axial localization of PD defects.

IV. THE METHOD TO LOCALIZE AXIAL POSITION OF PARTIAL DISCHARGE IN GIL BY GAS SENSORS

The above section gave the decay and delay characteristics of the decomposition products concentration in GIL. That could help to determine the position of PD via the difference in the concentration detected by gas sensors. As shown in Fig. 13, the gas sensors were arranged at 20 m spacing in a 100m GIL unit [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF]. That scheme ensured the presence of at least one sensor on each side of the defect. The GIL unit could be flexibly divided into several 20m-long monitoring segments. The sensor on the left side of the PD was LS, and the right one was RS. When the PD distanced LS and RS were equal, i.e., the defect occurred in the midpoint of the two sensors, the concentration detected by LS and RS varied consistently. In converse, there would be a difference in the variation of the concentrations detected by LS and RS. Based on [START_REF] Avzayesh | A hybrid estimation-based technique for partial discharge localization[END_REF], the distance between the PD and LS and the distance between PD and RS satisfied:

LS LR dp/SF6 RS LS 10 20 20 L t D LL           , ( 13
)
where ΔtLR (h) is the delay between monitored concentrations at LS and RS to reach the same concentration, LLS (m) is the distance between the defect and LS, and LRS is the distance between the defect and RS. LLS + LRS = 20 m (the gas sensors spacing). V. DISCUSSIONS

A. Interference of overheating faults on PD localization

Localized overheating fault may occur with PD together. It could cause the diffusion coefficients to be larger than prediction. To confirm the local overheating, the diffusion coefficient could be determined by introducing data from the third sensor. In this case, the distance between the third sensor and LS or RS was 20 m. According to the simulation model in [START_REF] Qiu | Methods for alleviation of impacts of axial diffusion on decomposition products monitoring in gas-insulated transmission lines[END_REF], the Δt of 20 m diffusion distance could be: Ddp/SF6=0.345/t.

(14) If the measured diffusion coefficient was larger than the predicted one, it could be inferred that GIL have localized overheating faults. When the temperature is hard to obtain for calculating diffusion coefficients, online measurement of diffusion coefficients could also be a alterative solution. However, the presence of a temperature gradient by localized overheating may cause diffusion coefficients change in the axial direction. The localization accuracy may be worse. That will be further studied in our future works.

B. Accuracy of PD localization method

To analyze the accuracy of the method, a series of simulations were done with different diffusion time, PD positions, and diffusion coefficients. The LLS was set as 2.5m, 5m, 7.5m, and 10m, respectively. Considering the symmetric equivalence, this paper did not specifically list the models with LLS over 10m (e.g., the case of LLS = 2.5m, i.e., LRS=17.5 m, is same with that of LLS =17.5m). According to the Table II, Ddp/SF6 was 5×10 -4 m 2 /s, 1×10 -3 m 2 /s, and 5×10 -3 m 2 /s. The generation rate was set as 10 ppm/h. For the different Ddp/SF6, the diffusion time was adjusted from 100h to 10000h. The localization errors are illustrated in Fig. 14. It could be found that the localization error was decreasing due the prolonged diffusion time (see Fig. 14(a) and Fig. 14(d)). In Fig. 14(b), the increase of diffusion coefficients could help to reduce localization errors. That could be explained by the following. [START_REF] We | Detection of SF6 decomposition components under partial discharge by photoacoustic spectrometry and its temperature characteristic[END_REF] was established on the premise that the increase gradient of the concentration was stable. The diffusion delay would not vary during the diffusion process. However, as shown in Fig. 9(b), it would take some time for the concentration growth rate to stable. Until the stabilization, the delay would increase with IEEE TRANSACTIONS ON DIELECTRICS AND ELECTRICAL INSULATION the prolonged diffusion time. In that case, the localization result of ( 13) would be inaccurate due to an unstable diffusion delay. As the diffusion time increases, the delay becomes stable, and the accuracy of (13) gradually increased. A higher diffusion coefficient would accelerate diffusion delay stabilization. Hence, a longer diffusion time and a larger diffusion coefficient would make localization more accurate. That is also the main reason why a larger diffusion coefficient resulted in a shorter required time for the error limitation in Fig. 14(c). In addition, the method has more accuracy for PD defects close the PD is to the midpoint between gas sensors. That is because the closer to the midpoint, the less the difference in concentration detected by two sensors. Then, the unstable delay has less effect on [START_REF] We | Detection of SF6 decomposition components under partial discharge by photoacoustic spectrometry and its temperature characteristic[END_REF]. In fact, Fig. 13 only disclose the theoretical error of the localization methods. There would be more disturbances and factors in the diffusion process in GIL (e.g., axial convection and temperature gradient). In addition, the insulators, especially the supporting one, would obstacle the diffusion process. That requires further simulation and experimental analysis. Our team would also conduct further research on the simulation with more details and the validation experiments on a real GIL in the future. Nevertheless, this paper still could demonstrate the feasibility of localizing PD defects through the gas sensors.

VI. CONCLUSIONS

The modified FSG model could predict the diffusion coefficients of decomposition products in high-pressure SF6. Among the eight typical decomposition products in SF6, HF had the fastest diffusion rate. The diffusion coefficients of H2S, SO2, COS, SOF2, CS2, and SO2F2 were similar. They could be monitored and compared for simultaneous analysis.

Under the effect of temperature rise and convection, the radial diffusion of decomposition products was fast. In contrast, axial diffusion is only affected by temperature and slow. The diffusion rate of decomposition products at the top of GIL was higher than at the bottom. The arranging gas sensors at the top of the GIL was suggested. For the slow axial diffusion, it has more concentration decay and diffusion delay, which lead to gas monitoring accuracy and timeliness problems. According to the simulation results, the diffusion process could be simplified to a 1-dimensional model. The diffusion coefficient could be determined by the average temperature rise. Diffusion delays have a quantitative relationship with the diffusion distance in the form of a quadratic function. With that relationship and the gas sensors' axial array, the defect could be localized by the concentration delay detected by the sensors on both sides of the PD defect. Localization accuracy is positively correlated with diffusion coefficient and diffusion time. In addition, the localization accuracy would increase when the defect is close to the midpoint of the sensors. If the gas diffusion is slow, a longer monitoring time is necessary for high accuracy.

  This work was supported in part by the National Key R&D Program of China 2021YFB2401400 and Science and Technology Project of China Southern Power Grid Co., Ltd. 090000KK52210175. (Corresponding author: Wenjun Zhou).

Fig. 1 .

 1 Fig. 1. The main factors of the diffusion process.

Fig. 2 .

 2 Fig. 2. Diffusion experiments platform.

Fig. 3 .

 3 Fig. 3. GIL temperature rises model settings. TABLE I

Fig. 4 .

 4 Fig. 4. GIL axial diffusion model.

Fig. 5 .

 5 Fig. 5. Simulation and experiment flow chart.

Fig. 5

 5 indicates the SO2F2 concentration at sampling points A and B under different time and pressures. The SO2F2 concentration increased slowly over time. Higher pressure would slow down the diffusion. The diffusion delays ΔtSO2F2 were 2.314 h, 2.663 h, and 3.025 h, respectively, when the pressure was 0.350 MPa, 0.400 MPa, and 0.450 MPa. With the diffusion distance of 1 m, the relationship between the diffusion coefficient and diffusion  IEEE TRANSACTIONS ON DIELECTRICS AND ELECTRICAL INSULATION delay is [16]: Dpd/SF6=2.077×10 -4 /t. (11) Based on (

Fig. 6 .

 6 Fig. 6. The concentration variation of SO2F2 at sampling points A and B under different pressure.

Fig. 7 .

 7 Fig. 7. The distribution of the temperature and convection in 550kV GIL at 5.6 kA. (a) The temperature distribution. (b) The convection distribution, and arrows indicate flow directions.

Fig. 8 .

 8 Fig. 8. The distribution of diffusion coefficients.

Fig. 9 (

 9 Fig. 9(a) and Fig. 9(b) display the diffusion range with and

Fig. 9 .

 9 Fig. 9. The distribution of diffusion coefficients and decomposition products. (a) The diffusion range with the temperature rise and convection. (b) The diffusion range without the temperature rise and convection.

Fig. 10 .

 10 Fig. 10. The SO2F2 concentration on the GIL shell for different defect positions.

Fig. 11 .

 11 Fig. 11. The decay and delay of decomposition products concentration in axial diffusion. (a) Decay. (b) Delay.

Fig. 12 .

 12 Fig. 12. The distribution of SO2F2 concentration on α, β, and γ lines, and the snapshot of axial diffusion process.

  389 10 ( 100) 1.389) /

Fig. 13 .

 13 Fig. 13. An arrangement of sensors in a 100m GIL.

Fig. 14 .

 14 Fig. 14. The localization errors under different conditions. (a) Variations of errors with diffusion time. (b) Effects of diffusion coefficients on the variation of error. (c) The minimum diffusion time for the accuracy above 90%. (d) The errors at different PD positions.

TABLE I MATERIALS

 I PARAMETERS IN THE MODEL

	Parameters (273 K)	SF6	Air	conductor/shell
	ρ (kg/m 3 )	22.82	1.293	2730
	λ0 (10 -2 W•m -1 •K -1 )	1.206	2.44	220(λs)
	ϑ (10 -5 Pa•s)	1.42	1.72	/
	Γ (J/kg•K)	665.18	1005	880
	ε	/	/	0.8/0.9
	Su (K)	243.8	110.4	/

TABLE II DIFFUSION

 II CHARACTERISTIC COEFFICIENTS

	Products	Kpd/SF6	Products	Kpd/SF6
	HF	4.82×10 -3	SO2F2	1.52×10 -3
	H2S	2.62×10 -3	S2OF10	8.98×10 -4
	SO2	1.96×10 -3	CS2	1.60×10 -3
	SOF2	1.69×10 -3	COS	1.90×10 -3

TABLE III COMPARISON

 III BETWEEN THE CALCULATION AND EXPERIMENT RESULTS OF 550 KV GIL TEMPERATURE RISES

				The temperature rises (℃)
	Current (kA)	Top conductor	Bottom of conductor	Top of shell	Bottom of shell
		Exp	Cal	Exp	Cal	Exp	Cal	Exp	Cal
	2.6	41	41	41	41	34	33	31	32
	4.4	57	60	57	59	37	37	32	34
	5.6	84	90	84	89	56	57	46	52

Rui Qiu was born in Harbin, China in 1992. He received the B.Sc. and master degrees from the Harbin University of Science and Technology, Harbin, China in 2015 and 2018, respectively. His research interests include the prediction of gas's dielectric strength, the decomposition and diffusion processes in SF6 and its substitutes.

He is currently pursuing the Ph. Currently, he is a postdoctoral fellow in Wuhan University.

Dr. Zheng is a member of IEEE and works as a secretary for IEEE PES T&D workshop of digital twin for power equipment.


Sanro Mathew

Maeliss Ethis De Corny

Nicolas Chauvet

Laureen Moreaud

Guillaume Laurent

Serge Huant

Erik Dujardin

Gilles Nogues

Guillaume Bachelier

Origin of second harmonic generation in gold plasmonic nanostructures

Keywords: Plasmonics, Nonlinear optics, Second-harmonic generation, Gold, Nanoprisms

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Nonlinear optical processes has been the subject of research for many decades owing to an interest in the pursuit of furthering the fundamental science as well as for its applications. In this context, second-order processes including second harmonic generation (SHG), difference frequency generation (DFG), sum frequency generation (SFG) and spontaneous parametric down conversion (SPDC) has been well studied in bulk media and utilized for frequency conversion as well as for entangled photon pair generation.

1 2 More recently, the advancement in nanofabrication techniques spurred by the growing demand for down scaling components to nanoscale size has led to a renewed interest in the research of nonlinear optics at the nanoscale. Owing to the ability of plasmonic metal nanostructures to confine and enhance electromagnetic field at the sub-wavelength scale, generation of nonlinearity, in particular SHG, from metal nanostructures has therefore been of significant interest.

3 SHG is a process in which an incident wave at frequency ω is converted to an outgoing wave at 2ω. SHG is sensitive to material symmetry. It is forbidden in centrosymmetric materials within the electric dipole approximation. [START_REF] Boyd | Nonlinear Optics[END_REF] Hence, the nonlinearity arises through different mechanisms and its description in metal nanostructures which are centrosymmetric has been the subject of numerous investigations. Briefly, these mechanisms correspond to symmetry breaking at the surface and higher order effects from multipoles within the bulk (magnetic dipole and electric quadrupole). 4 5 The separation of the effects from the different mechanisms has been addressed through several investigations. In the case of thin gold film, surface effects (arising from the normal component of surface polarization) is found to be the dominant contribution to the SHG response, [START_REF] Fu | Surface and bulk contributions to the second-order nonlinear optical response of a gold film[END_REF] while, in the case of split ring resonators, bulk contribution is found to be dominant. [START_REF] Feth | Second-harmonic generation from complementary split-ring resonators[END_REF] Regardless, in general it is widely believed that the contribution from normal surface component always dominate the nonlinear response as was concluded in a recent theoretical work.

9

In this report, we compare experimental and simulated far-field SHG intensity distribution obtained under tightly focused laser beam excitation, and report that the dominant nonlinear contribution in gold nanoprisms arises from the parallel surface component which was considered to play a less significant role in the SHG process.

THEORETICAL FRAMEWORK

Due to the breaking of the inversion symmetry of centrosymmetric materials as in the case of metal interfaces, the conduction electrons experience an asymmetric potential that leads to second order nonlinearity in the induced currents. [START_REF] Rudnick | Second-Harmonic Radiation from Metal Surfaces[END_REF] Thus, for any isotropic material with an inversion symmetry, the second harmonic (SH) surface polarization can be expressed through three components of the susceptibility tensor as:

10 P surf ace (r, 2ω) = 0 χ ⊥⊥⊥ E ⊥ (r, ω)E ⊥ (r, ω)n + 0 χ ⊥ E (r, ω)E ⊥ (r, ω)k + 0 χ ⊥ E (r, ω)E (r, ω)n, (1)
where ⊥ and represents components perpendicular and parallel to the surface and n is the unit vector normal to the surface at point r while k being the unit vector parallel to the surface. For instance, χ ⊥ denotes surface polarization perpendicular to the surface generated by components of electric field parallel to the surface. Although the susceptibility tensor is expressed through three components, the existing models of electron dynamics at the surface based on hydrodynamic model predicts χ ⊥ to be vanishing. [START_REF] Bachelier | Origin of optical second-harmonic generation in spherical gold nanoparticles: Local surface and nonlocal bulk contributions[END_REF] Hence, the surface polarization reduces to:

P surf ace (r, 2ω) = 0 χ ⊥⊥⊥ E ⊥ (r, ω)E ⊥ (r, ω)n + 0 χ ⊥ E (r, ω)E ⊥ (r, ω)k, (2)
Additionally, SH sources arising from the bulk due to higher order effects such as electric quadrupole in the case of isotropic and centrosymmetric material is expressed as:

6 P bulk (r, 2ω) = β E(r, ω)[ ∇ • E(r, ω)] + γ bulk ∇[ E(r, ω) • E(r, ω)] + δ [ E(r, ω) • ∇] E(r, ω), (3)
Within the bulk, the first term of Eq.3 is vanishing whereas the third term corresponding to δ has a nonlinear susceptibility that is an order of magnitude smaller than the second term in the case of single beam SHG. [START_REF] Ethis De Corny | Wavemixing origin and optimization in single and compact aluminum nanoantennas[END_REF] Hence, the bulk polarization can be reduced to:

P bulk (r, 2ω) = γ bulk ∇[ E(r, ω) • E(r, ω)], (4)
The bulk contribution expressed above is often known as non-local bulk because the gradient represents the spatial variation of square of the electric field. The nonlinear susceptibilities in the surface polarization and the material parameter β in the non-local bulk are expressed using Rudnick and Stern parameters as:

4 12 χ ⊥⊥⊥ = - a 4 [ r (ω) -1] e 0 mω 2 , χ ⊥ = - b 2 [ r (ω) -1] e 0 mω 2 , γ bulk = - d 8 [ r (ω) -1] e 0 mω 2 , (5)
where r (ω) is the dielectric function of the material at frequency ω while e and m corresponds to electron charge and mass. The a, b and d are the non-dimensional Rudnick and Stern parameters included to weigh the relative contributions from the different nonlinear sources. Following the work of Bachelier et al., the values of these parameters are taken to be a=0.56-0.25i, d=1 and b=-0.1.

10

As these sources are weighed with different parameters, their separability allows us to evaluate their independent contributions and to simulate their intensity distribution to identify their role in SHG. To this end, it was shown that non-local bulk polarization can be reformulated into an effective surface contribution as:

12 P ef f,⊥ (r, 2ω) = [χ ⊥⊥⊥ -χ ⊥ ] E ⊥ (r, ω) E ⊥ (r, ω) + [χ ⊥ + γ b (2ω) ] E(r, ω) E(r, ω), (6)
Conventionally, the polarization sheet is placed outside and consequently the nonlinear currents are evaluated outside, while, the fundamental fields are evaluated inside the material. As χ ⊥ is vanishing, the effective polarization can then be written as:

P ef f,⊥ (r, 2ω) = χ ⊥⊥⊥ E ⊥ (r, ω) E ⊥ (r, ω) + γ b (2ω) E(r, ω) E(r, ω), (7)
As it is evident in this expression, the nonlocal bulk contribution and surface contributions can thus be separated and this allows us to assign the leading contribution to SHG.

FINITE ELEMENT SIMULATIONS

The separability of the nonlinear contributions thus allows us to simulate them independently to obtain their SHG response maps. The simulations were performed combining analytical methods (via MATLAB) and finite element method (via COMSOL). The simulation was designed to take into account the optical experimental setup that was developed. In a first step, the background field at the fundamental frequency (ω) is evaluated analytically by considering its transmission through a microscope objective, immersion oil and refraction through a glass substrate hosting the nanostructures. Thereafter, using the scattered field formulation, Maxwell's equations are solved to obtain the total field at the fundamental frequency. Then, the nonlinear source currents are computed taking into account the nonlinear susceptibilities given above in the expressions 5-7. In the next step, Maxwell's equations are again solved using these nonlinear sources to compute the harmonic field (2ω) and then the response currents generated by the harmonic field is computed which thereby gives the total current generated in the nanostructure. Finally, the field radiated by the nonlinear currents into the far field is evaluated analytically by taking into account its transmission back through the substrate and the objective lens onto a detector to obtain the SHG response.

NANOFABRICATION AND EXPERIMENTAL PROTOCOL

The protocol used for the fabrication of nanoprisms involves standard E-beam lithography and nanofabrication procedures as well as chemical synthesis for crystalline prisms. The samples were fabricated on a 170µm glass substrate (Zeiss 1.5 coverslip). In a first step, the substrates were cleaned in Nitric acid bath for one minute followed by rinsing using deionized water. Thereafter, substrates immersed in acetone were put in an ultrasound bath for three minutes. In a final step, the substrates were then subjected to an oxygen plasma cleaning for five minutes. Once the cleaning was performed, the substrates were then spin coated with 3% PMMA resist at 4000 rpm/s for 30 seconds and hard-baked for five minutes at 180 • C.

To dissipate charges from the E-beam exposure, the substrates were then spin coated with an electronic conductive resist Elektra92 (ARPC 5090) at 3500 rpm/s for 60 seconds followed by a soft bake at 90 • C for two minutes. Following the E-beam exposure, the conductive resist was removed by rinsing in deionized water for one minute. Thereafter, the exposed resist is developed in Methyl Isobutyl Ketone and Isopropanol (MIBK/IPA) mixture composed in the ratio 1:3 for one minute which is followed by a rapid rinse in isopropanol in a first bath and then two minutes in a second bath of isopropanol as well. In the next step, Gold is deposited using electron gun evaporation of the desired thickness which is finally followed by a lift-off of the unexposed resist, in which, the substrate is immersed in N-methyl-2-pyrrolidone at 80°C for one hour and then rinsed with acetone and isopropanol. On the other hand, chemically synthesized prisms were prepared by using polyvinylpyrrolidone (PVP) and potassium iodide (KI) to reduce gold salt at room temperature. The desired shape, in our case, prisms, are obtained by controlling the pH of the solution.

Additionally the dimension of the structures is controlled by the length of solution resting time. Finally, the prisms deposited at the bottom of solution are collected and dispersed on to a substrate. The characterization of the nanoprisms were performed on an optical setup designed in a way to reach single particle sensitivity i.e to be able to perform single particle measurements. The signal from nanoprisms were collected through an oil immersion objective in reflection mode. For this purpose, a 100x objective of numerical aperture 1.3 is used, thus, providing diffraction limit resolution of 340nm and 170nm corresponding to the fundamental and harmonic wavelengths. The objective and the sample is mounted on a piezo electric stage for longitudinal and lateral scanning. The chromatic aberration induced by the objective is corrected by using motorized telescopes to overlap the focal planes at the excitation and the collection wavelength.

In addition to the motorized telescopes and piezoelectric stage control, a tracking algorithm is developed to automatically optimize the signal acquisition. A pulsed femtosecond Ti:Sapphire oscillator was used as the laser source which delivers 100fs pulses at a repetition rate of 80MHz. The average input power was set at 200µW to avoid excessive damage on the nanoprisms of interest.

RESULTS AND DISCUSSION

In order to probe the origin of nonlinear response from nanoprism antennas, we compare the cartographic SHG maps acquired experimentally to that of the FEM based numerical simulations. [START_REF] Ethis De Corny | Wavemixing origin and optimization in single and compact aluminum nanoantennas[END_REF] The approach used here is to acquire the SHG intensity distribution over the entire nanoprism by simultaneously exciting and collecting the signal from the same positions. In other words, the SHG map does not correspond to a single plasmon mode rather it is a result of different plasmon modes excited as a function of position over the nanoprism. The simulation was therefore developed to account for the exact experimental configuration as briefly described in section 3). Figure 1 shows nanoprisms of size 707nm that was investigated. Chemically Figure 1: 1a is an AFM image of a Gold nanoprism fabricated using E-beam lithography and 1b is an AFM image of a chemically synthesized Gold nanoprism synthesized nanoprisms and lithographic nanoprisms of the same size were studied. The protocol used for the fabrication of these two nanoprisms were presented in the section 4. Figure 2 represents an experimental SHG map on the left acquired on chemically synthesized prisms which are then compared to the simulated individual nonlinear contributions as described in equations ( 5) presented on the right hand side. As the simulation does not take into account the transmission loss through various optical elements, the transmission efficiency of microscope objective and detection efficiency of the photon counting modules, the calculated signal counts represented in the simulated maps are normalized with the highest experimentally obtained intensity counts. Additionally, the red arrow in figure 2 represents polarization of the excitation field normal to the horizontal axis of the prism while blue arrow represents polarization of the detected photons that are normal and parallel to the same axis. In figure 2 one can note that the experimental SHG maps does not change with the change in detection polarization indicated by the blue arrows and there is only a single prominent lobe. In comparison with the simulations, the bulk contribution can be immediately ruled out as the intensity distribution itself is different from what is observed in the experimental map. Whereas, the discrimination between χ s ⊥ and χ s⊥⊥⊥ is more ambiguous. However, noting that the intensity counts of the bulk contribution is dominant, it is evident that there is an apparent contradiction in our interpretation of the comparison made between the SHG maps. To explain this inconsistency, recall that the chemically synthesized nanoprisms are monocrystalline in nature, i.e they are essentially made up of single crystals arranged on top of each other thereby producing a low surface rugosity as it can be seen in the AFM image.

In addition to it, recall that the χ s ⊥ component represents parallel surface contribution generated by electric field components parallel and perpendicular to the surface. This contribution therefore arises from the parallel surface currents. Putting together these two points, rough surfaces displaying irregularities on the surface in the form of pits and bumps could affect the parallel current flow. The rough features could Taking this into account, figure 3 shows a coherent sum of the three nonlinear components by varying the 'b' parameter. Here, it can be seen that the two symmetric lobes diminishes in intensity with b. In the experimental map, one can note the absence of a distinct three lobe pattern lobe which could be due to the characteristic nature of chemically synthesized prism. Regardless, the intensity ratio of the experimental map on the left gives a value of 2.8. This ratio is closer to those given for b = 0.5 in the simulated maps. According to the literature, 4 b = 0.5 corresponds to a flat surface but which diffusely scatters electrons. Increasing the b further close to unity, one can see that contribution from parallel surface component augments the intensity of the bright lobe whereas the two side lobes diminishes in intensity. Therefore, from the intensity ratio, it is evident that the value of b that fits closest to the experimental data is b=-0.5. To decouple the effect of roughness from the flatness of a surface, prisms of the same size as crystalline prisms were fabricated by E-beam lithography. As the nucleation and subsequent layer growth of deposited Gold atoms/molecules in electron gun evaporation does not yield a lattice like structure similar to a crystalline growth, the surface is therefore expected to exhibit significant roughness with increasing thickness. In figure 4 SHG maps obtained with lithographed prisms of varying thickness are compared with crystalline prisms for the same excitation wavelength and polarization shown in the previous figures as well as simulated SHG maps for the larger thickness. As it can be seen, the overall intensity distribution is similar across the three prisms. Additionally, the prism with a larger thickness displays a distinct and comparatively more intense three lobe response which is in agreement with the result obtained from the simulations. This agreement also indicates that the effect of roughness on the response is absent or negligible zince the simulation does not take into account the surface rugosity. One explanation put forth is that, the overall shape or the long range order of the surface of lithograph prism is flat but that contains horizontal features primarily determined by the grain size of Gold due to the process of electron gun deposition. From the point of view of electron scattering, the horizontal distance of a surface feature must be on the order of v f ω to affect the parallel surface current 4 . To put into perspective, for Gold given a wavelength of 1000nm, the value corresponds to 5 Å-10 Å. While the resolution at the angstrom level is inaccessible via AFM, one can note that the grain size of Gold via evaporation is larger than few tens of nanometer which is measurable. Thus, essentially one could argue that the horizontal feature distance is too big to observe the effect from electron scattering and consequently the surface roughness. This, therefore, is in agreement with the conclusion arrived for crystalline prisms while noting that the effect of roughness is not significant in the case of lithographic nanoprisms as it was expected.

Crystalline prism

Lithograph prism Lithograph prism Simulation Thickness: 22nm 22nm 35nm 35nm

Figure 4: Simulated SHG maps and experimental maps corresponding to crystalline and lithographed gold nanoprisms of size 707nm Gold nanoprisms excited at 880nm with polarization perpendicular to the horizontal axis of nanoprisms in the x-y plane.

Conclusion

The generation of SHG in metals is attributed to three main contributions. Among the three, two contributions arise from the surface whereas the third comes from the bulk. Often in the literature, surface contribution, in particular, the one arising from the normal currents are attributed to have the dominant role in SHG because of the breaking of centrosymmetry at the surface resulting in strong electric field gradient field discontinuity. In this report, we compare numerically simulated and experimentally acquired cartographic SHG response maps over a nanoprism to distinguish and determine the dominant SHG source contribution. The map is acquired by simultaneously exciting and collecting different positions on the nanoprism by scanning it under a high numerical aperture microscope objective. The resulting map obtained is therefore an SHG intensity distribution or the efficiency of harmonic response as a function of position over the nanoprism. The experimental configuration is quantitatively accounted for by means of numerical simulations based on finite element method whereby the substrate, high numerical aperture, reflection mode of the experiment are taken into account. In our experiment, chemically synthesized and lithographic nanoprisms of Gold of the same size were investigated. We demonstrate that the SHG source arising from the χ s ⊥ surface contribution is the dominant contribution in the nanoprisms by comparing the response maps of the individual nonlinear SHG source contributions as well as their coherent response. Contrary to the literature, intensity from χ s⊥⊥⊥ is the weakest while that corresponding to the bulk γ has a non-trivial contribution and that from χ s ⊥ remains dominant. Additionally, we verify the effect of roughness on the response by performing the same experiment on a nanoprism of larger thickness since the physical vapor deposition of Gold yields rougher layer with increasing thickness. The approach thus adopted in this investigation can also be extended to investigate the harmonic response from nanoantennas of any shapes or configurations of interest.

BulkFigure 2 : 1 Figure 3 :

 213 Figure 2: Represented on the left are maps obtained on crystalline prisms. On the right side are the simulations corresponding to the three individual nonlinear components. The red arrow represents polarization of the excitation field and blue arrow represents the polarization of the detected photons with respect to the horizontal axis of the prism in the x-y plane. The excitation wavelength is set at 880nm.


domain: spi phys

text:
Jean-Marc Chaix

Renaud Bouchet

Didier Bouvard

Timothée Fabre

Tristan Garnault

Christelle Harnois

Kouami Koutoati

Marie Lachal

Sylvain Marinel

Marlu Cesar Steil

A viewpoint on hot spots in microwave sintering and flash sintering

Keywords: Flash sintering, microwave sintering, hot spots, ceramics

Microwave sintering and flash sintering are field assisted ceramic processing techniques in which heating is totally or partly obtained inside the material itself by a coupling between the ceramic material and the electromagnetic or electrical field. They frequently give rise to the observation of localised temperature excess known as hot spots that can lead to dramatic effects on the sintered parts. This paper explores and discusses the common origin of hot spots in both techniques, their development and the possible ways for limiting their effects.

Introduction

Microwave sintering (MWS) [1]- [START_REF] Chaix | Ref. Module Mater. Sci. Mater. Eng[END_REF] and Flash Sintering (FS) [START_REF] Dancer | [END_REF]- [6] are two techniques that involve the direct coupling between the ceramic powder parts to be heated and an electromagnetic or electrical power source. The effective heat source is therefore inside the bulk of the part and depends on the interaction of the constitutive matter with the electromagnetic field. These techniques are potentially energy efficient for sintering, as the heating a whole furnace up to the sintering temperature is not required, and they are expected to enable fast sintering, due to the limited thermal inertia of the parts. Both techniques however exhibit the problem of hot spots, i.e., the development of local overheating in some regions of the part. It consists in a localization of coupling/heating and is promoted by the high heating rates. Specific additional "non-thermal" effects of the field on the sintering mechanisms and kinetics have been proposed and discussed in a number of papers, as well for FS [7] [8] as for MWS [9]- [12] , while other papers discuss the positive effect of the enabled high heating rates [13] [14] but this discussion is out of the scope of the present paper. The purpose of the present paper is to discuss the common and specific aspects of the hot spot phenomenon in the two techniques and the ways that have been proposed to deal with this problem, in the light of literature data and new results.

FS and MWS techniques FS

Since the early work by Cologna and Raj [15] who first described and named the Flash Sintering technique, FS is obtained by applying a DC or AC field to a ceramic powder compact preheated at a temperature much lower than its usual sintering temperature but high enough for ensuring a non-negligible electric conductivity [START_REF] Dancer | [END_REF]- [6] . Many papers deal with specimens with a dog-bone shape, well suited for direct observation and easy to setup, but classic cylindrical shapes are also used (Figure 1-a). The flash phenomenon is applied using a constant voltage V either when the temperature is stabilized at T0 (isothermal experiments) or at the beginning of a temperature ramp (constant heating rate). In the first case, when suitable (high enough) T0 and V are selected, a first period is observed in which the measured current through the specimen slowly increases (the so-called incubation period), before a dramatic increase in current suddenly occurs (Figure 1-b). This current surge is associated to a strong increase of temperature and a fast shrinkage of the specimen, easily observed in the case of dog-bone samples (see for instance videos by J Francis [START_REF] Francis | Small Flash Movie[END_REF] ). A flash of visible light emission is also observed. Although other interpretations have been proposed for the incubation time, sometimes with very sophisticated physical models [START_REF] Naik | [END_REF] , it is now admitted [18]- [22] that the thermal runaway induced by the Joule effect and the thermally activated electrical conductivity of the ceramics explains the coupled electrical and thermal surge as well as the flash in light. Indeed, if we consider a small thermally homogeneous and isolated sample submitted to a constant electric field 𝐸𝐸, the local Joule power density is

𝑃𝑃 = σ𝐸𝐸 2
(1) where σ is the electrical conductivity and 𝐸𝐸 the electric field, so that the temperature evolution is

d𝑇𝑇 dt = σ𝐸𝐸 2 ρC p (2)
where 𝜌𝜌𝜌𝜌 𝑝𝑝 is the volume heat capacity. When 𝜎𝜎 is thermally activated with an activation energy 𝑄𝑄, this leads to a simple evolution equation:

𝑑𝑑𝑇𝑇 𝑑𝑑𝑑𝑑 = 𝜎𝜎 0 𝐸𝐸 2 𝜌𝜌𝜌𝜌 𝑝𝑝 exp �- 𝑄𝑄 𝑅𝑅𝑇𝑇 � (3)
which can be easily numerically integrated (Figure 1-c) and clearly reproduces the experimental curves. This simple equation also explains the shorter "incubation time" for higher temperature 𝑇𝑇 or higher electric field 𝐸𝐸. It however does not explain the experimentally observed threshold value for the occurrence of the flash phenomenon. This last observation requires to take into account the thermal losses term 𝐿𝐿, that depends on the specimen temperature 𝑇𝑇 and the furnace temperature 𝑇𝑇 0 (typically terms in (𝑇𝑇 -𝑇𝑇 0 ) and (𝑇𝑇 4 -𝑇𝑇 0 4 ))

for conductive/convective and radiative exchanges, respectively). It comes:

𝑑𝑑𝑇𝑇 𝑑𝑑𝑑𝑑 = 𝜎𝜎 0 𝐸𝐸 2 𝜌𝜌𝜌𝜌 𝑝𝑝 exp �- 𝑄𝑄 𝑅𝑅𝑇𝑇 � -𝐿𝐿(𝑇𝑇, 𝑇𝑇 0 ) (4)
The flash phenomenon is only observed when the heat source term exceeds the loss capacity [20][21] [23] , so that the thresholds in 𝑇𝑇 or 𝐸𝐸 correspond to a unique threshold in initial Joule power 𝜎𝜎(𝑇𝑇 0 )𝐸𝐸 2 [24] .

Keeping the applied voltage constant when the flash occurs would lead to an infinite Joule power. In practice, in most experiments the power supply is switched to a constant intensity control, when the intensity reaches a predefined value 𝐼𝐼. The "classic" flash sintering technique therefore includes two steps, known as voltage control and current control stages, respectively (Figure 1-b). The intensity control step is generally maintained for several tenths of seconds, in order to improve the final density, as the densification during the short flash is significant but far from reaching full density. It is to be noticed that the switch intensity control, performed by an electronic device, is not instantaneous, so that uncontrolled overcurrent/overheating is unavoidable and systematic. evolution of temperature and current during a FS experiment calculated from Equation 3 using data for 3YSZ from Todd et al. [18] (c).

MW

Microwave sintering [1]- [START_REF] Chaix | Ref. Module Mater. Sci. Mater. Eng[END_REF] is operated in a microwave cavity, to which microwaves produced by a generator, generally a magnetron and more recently a solid state generator, are brought through a coaxial or rectangular waveguide. Two major types of cavities are used. The single mode cavity is the semi closed section of a waveguide, whose size can be adjusted to approach resonant conditions. In this case, the electric field distribution corresponds to stable stationary waves (Figure 2-a) with defined high field nodes and antinodes. In such cavities, specimens with a size smaller than the half wavelength can be placed in almost constant electric fields areas. The second type of cavity, called multimode cavity, is a large one, generally a parallelepiped, in which multiple reflections of the waves on the sides lead to multiple propagation modes and to a complex and almost unpredictable distribution of the field (Figure 2-b). Several or large specimens can be placed at the same time in the multimode cavity; they are often placed on rotating plates in the purpose of obtaining a homogeneous time average of the field around it.

In this technique, the specimen can be directly heated by the interaction with the field, in most cases the electric field, the power density 𝑃𝑃 being:

𝑃𝑃 = 2𝜋𝜋𝜋𝜋𝜀𝜀 0 𝜀𝜀 𝑒𝑒 ′′ 𝐸𝐸 2 2 ( 5
)
where 𝜀𝜀 0 is the dielectric permittivity, 𝜋𝜋 the wave frequency and 𝜀𝜀 𝑒𝑒 ′′ the imaginary part of the effective dielectric permittivity of the specimen. 𝜀𝜀 𝑒𝑒 ′′ is often written 𝜀𝜀 𝑒𝑒 ′′ = 𝜀𝜀 𝑒𝑒 ′ tan 𝛿𝛿 , where 𝜀𝜀 𝑒𝑒 ′ is the real part of the relative permittivity and tan 𝛿𝛿 the so-called loss factor. In Equation 5, 𝐸𝐸 is the electric field inside the material: it is in practice different from the externally applied field [START_REF] Chaix | Ref. Module Mater. Sci. Mater. Eng[END_REF] , but this point will not be discussed in the present paper.

In the case of microwave sintering, as high temperatures are required and the cavity walls are cold, the specimen is generally embedded in a thermally insulated box constituted by a highly porous fibrous or porous material, which is supposed to have low interaction with the field [3][25] . A major point for the present paper is that the effective dielectric permittivity, which includes a contribution 𝜎𝜎 2𝜋𝜋𝜋𝜋𝜀𝜀 0 of the electrical conductivity 𝜎𝜎, strongly increases with increasing temperature, almost as a thermally activated property, in particular when this increase is due to the electrical conductivity [3][10] . This temperature dependency is an issue for the control of the heating rates during microwave sintering [26]- [28] . In the same way as in Joule heating, the microwave heating is therefore subject to thermal runaway phenomena, as demonstrated by Roussy et al. [29] using a simple model leading to practical criteria in the case of a second degree polynomial dependency of permittivity with temperature, or by Spotz et al. [30] in the case of alumina.

A direct link with FS was recently made by Bykov et al. [31] who suggested to use the thermal runaway in microwave heating to perform microwave flash sintering on the basis of an efficient control of the microwave input power.

Hot spots/localization of heating

Basic instability

The tendency to produce local overheating is a common feature to MWS and FS that reflects the intrinsic instability of the heating by direct coupling when the coupling parameter (𝜀𝜀 𝑒𝑒 ′′ or σ, respectively) increases with temperature. A similar equation defines the dissipated power density in both cases (Equation 1 and 5):

𝑃𝑃 ∝ γ(𝑇𝑇) 𝐸𝐸 2 (6)
where the coupling property γ(𝑇𝑇) of the material has a positive temperature coefficient 𝑑𝑑γ 𝑑𝑑𝑇𝑇 . The instability can be understood using a simplified calculation. If we consider a small volume δ𝑉𝑉 of the sample with a small perturbation of temperature δ𝑇𝑇 with respect to its surroundings (Figure 3), this volume will exhibit a power excess δ𝑃𝑃 with respect to its surroundings:

δ𝑃𝑃 ∝ 𝑑𝑑γ 𝑑𝑑𝑇𝑇 𝐸𝐸 2 δ𝑇𝑇 δ𝑉𝑉 (7)
while the outwards heat loss δ𝐿𝐿 by heat diffusion (heat conductivity κ) along distance δ𝑙𝑙 through the surface δ𝑆𝑆 of the small volume is expected to dissipate this temperature excess:

δ𝐿𝐿 ∝ κ δ𝑇𝑇 δ𝑙𝑙 δ𝑆𝑆 ( 8
)
Where δ𝑙𝑙 is the heat diffusion distance. The temperature perturbation will be damped when (δ𝑃𝑃 -δ𝐿𝐿) < 0 and amplified if (δ𝑃𝑃 -δ𝐿𝐿) > 0. In the case of positive temperature coefficient, (δ𝑃𝑃 -δ𝐿𝐿) can be positive and its sign is determined by the amplitude of the dimensionless power ratio

δ𝑃𝑃 δ𝐿𝐿 ∝ 𝑑𝑑γ 𝑑𝑑𝑑𝑑 𝐸𝐸 2 κ δ𝑉𝑉δ𝑙𝑙 δ𝑆𝑆 (9)
As both δ𝑉𝑉 δ𝑆𝑆 and δ𝑙𝑙 scale with the size 𝜆𝜆 of the small volume, the dimensionless power to loss ratio Π can be written:

Π = 𝑑𝑑γ 𝑑𝑑𝑑𝑑 𝐸𝐸 2 κ 𝜆𝜆 2 (10)
Π determines the stability of the system. For large values of this ratio, any small over heating will be amplified, leading to heat localization and the formation of a hot spot. This oversimplified calculation helps understanding the conditions that favour the development of instabilities: high fields, i.e., high power and high heating rates ; high positive temperature coefficients 𝑑𝑑γ 𝑑𝑑𝑇𝑇 ; low heat conductivity κ, i.e., low heat dissipation. Equation 10 also means that the perturbation size 𝜆𝜆 must also be large enough, which of course requires that the sample itself is large enough. In any given processing conditions of a given material a large sample must systematically be unstable with respect to hot spots. These simple conclusions are valid for FS and MWS as well.

In the case of FS, Dong [32] proposed a more classic perturbation analysis of stability and obtained a critical size in the case of Joule heating in a 2D geometry (perpendicular to the 𝐸𝐸 field axis), with thermally activated conductivity 𝜎𝜎 = 𝜎𝜎 0 exp (-𝑄𝑄 𝑅𝑅𝑇𝑇 ) (a typing error is present in equation 10 and 11 in Dong's original paper [32] , where 𝜆𝜆 should be read 𝜆𝜆 2 ; this error is reproduced in some citing papers).

𝜆𝜆 𝑐𝑐 = � 2κ 𝐸𝐸 2 𝜎𝜎 0 𝑄𝑄 𝑅𝑅𝑑𝑑 0 2 exp (- 𝑄𝑄 𝑅𝑅𝑑𝑑 0 ) � 1/2 (11)
that could be written, in the present notations:

𝜆𝜆 𝑐𝑐 = � 2κ 𝐸𝐸 2 𝑑𝑑γ 𝑑𝑑𝑑𝑑 � 1/2 (12)
This is typically equivalent to setting that the system is unstable when Π > Π 𝑐𝑐 = 4. Dong [32] estimated 𝜆𝜆 𝑐𝑐 to typically 0.8-1.6 mm around 1000°C for YSZ, but it is to be noticed that in the case of thermally activated conductivity, 𝑑𝑑γ 𝑑𝑑𝑇𝑇 strongly increases with temperature, so that even smaller samples could be the seat of instabilities when temperature increases during FS. In the case of MW heating of alumina, Spotz et al. [30] also discussed a critical specimen size through the use of the Biot number. The intrinsically unstable character of microwave heating and the development of hot spots was also studied and discussed on the basis of computer simulations [33] . ; local melting of the thermal insulation around the specimen due to localised heating of the insulating material (c). [START_REF] Garnault | Frittage Par Micro-Ondes de Céramiques Mises En Forme Par Impression 3D[END_REF] Development of the localisation

In the case of microwave heating, the instability generally develops as a hot spot observed on the surface of the specimen (Figure 4-a) and the local overheating can lead to the fracture of the specimen (Figure 4-b) [28] . Hot spots can also develop inside the specimen, as shown in experiments as well as in numerical calculations [28][33] . They can also develop inside the insulating material placed around the specimen (Figure 4-c) in which the coupling with microwave is expected to be low but the very low thermal conductivity favours the development of hot spots [28] .

In FS experiments, when the local thermal instability turns into localized thermal runaway, the localization develops along a preferential current path (high temperature = high conductivity) which results in a "hot line" rather than a hot spot [35][36] . The process therefore leads to strongly inhomogeneous samples (Figure 5): they are almost not densified but a region around a preferential current path that shows densification and grain growth, with sometimes melted zones [35][37] , cavities [6] and even fractures or as-drilled holes [START_REF] Becker | [END_REF] . [START_REF] Bichaud | Frittage « flash » de Céramiques Sous Courant Alternatif[END_REF] ; top and cross section of a LVP (Li3V2(PO4)3) pellet showing the localisation of densification and fracture (b).

Jerby et al. have studied the development of the thermal instability into a "localized microwave heating" by experiments and numerical models and qualified the phenomenon as an illustration of the "Matthew effect" [START_REF] Jerby | [END_REF] , in reference to the Matthew apostle question "why the rich get richer and the poor get poorer?", an effect whose applications go far beyond thermal runaway [40] . In FS as well as in MWS, when the field is kept constant, the power is increasing with temperature, so that the higher the temperature, the higher the temperature increase. This does not only lead to an increase of the temperature difference (Figure 6-a), but also to the shrinking of the hot region as shown on normalised curves (Figure 6-b).

It can be noticed that high temperature localisation leads to locally faster densification due to sintering. As higher density means higher electrical conductivity and higher dielectric permittivity, sintering phenomena tend to enhance the localisation. The simulations presented in the present paper are performed using the electrical conductivity for YSZ proposed by Todd et al. [18] , i.e., σ = 9302 exp(171kJ.mol -1 /RT) S/cm -1 .

Is there a limit to this temperature concentration increase? A stationary temperature profile may be reached at which the heat loss equals the electromagnetic power source. Hewitt et al. [21] have shown that in the case of FS, for a long cylindrical specimen of radius 𝑟𝑟 with side heat exchange with ambient, a steady state is always expected in the fixed current situation, while in the fixed voltage configuration a steady state is only obtained when the dimensionless "strength of the ohmic heating"

𝐸𝐸 2 𝜎𝜎 0 𝑄𝑄 𝑅𝑅𝑑𝑑 0 2 κ
𝑟𝑟 2 (using the same notations as above) is small enough. It can be noticed that this dimensionless parameter is almost the same as the above mentioned power to loss ratio Π in Equation 10, here applied to the specimen characteristic size. This is consistent as the stability analysis and the existence of a stationary state are indeed not really different problems. Anyway, even when a steady state is mathematically expected, events such as local melting or fracture of the material may occur before the steady state is reached.

Sources of instability nucleation

Thermal exchanges

The bulk heating in both MW and FS naturally produce a thermal gradient often mentioned as "inverted gradient" with respect to the usual gradient obtained by surface heating in a classic furnace with hot walls. This inverted gradient provides an inhomogeneous temperature profile at the scale of the specimen size, suitable to trigger an instability in any large enough specimen. In the case of FS, these thermal gradients were observed in early modelling results in FS [41][42] and clearly analysed as the source of the development of instabilities in a 1D axisymmetric model [18] . Other papers only obtain limited thermal gradients in calculations, which can be related to the small thickness and width of the dog-bone shaped modelled samples [43] . In MWS, although field penetration limitations have to be taken into account in some cases, simulations clearly showed the development of such bulk hot spots [28][33] [44] . Inverted gradients are a major source for nucleating centred hot spots in both FS (Figure 6-c) and MWS. Off-centre hot spots are however frequently observed (Figure 5), revealing that other sources for nucleating hot spots may be present.

Field inhomogeneity

Defects in the field distribution are other sources of instability nucleation. The uncontrolled field non uniformity in multimode MW cavities (Figure 2) is one of these sources, often discussed in review papers [9][45] . The specimen is itself a source of inhomogeneity as the field distribution outside and inside results from the interaction of the microwave applied field and the specimen [28][44] [46] . The field, and therefore the dissipated power depend on the specimen shape and position in the cavity and may exhibit strong inhomogeneity even at the beginning of heating (Figure 7). In FS, the field results from the application of the voltage at contacts. When the system evolution departs from symmetry, the initial perturbation is not centred: for instance in a dogbone shaped specimen in horizontal position hanged by the electrical wire (Figure 8-a), the slight off-centring is enough to induce an upward deviation from the current flow, so that the non-axially developed instability leads to a curved sample [19] . The quality of contacts is also a significant potential source of defects leading to localization (Figure 8-b): different contact configurations may lead to a variety of final defects in the FS processed parts [START_REF] Becker | [END_REF] .

a b

Composition inhomogeneity

Composition inhomogeneity is the third typical source of localisation nuclei. In MWS, the presence of particles of a highly absorbing material in a weakly coupling matrix is a cause of hot spots [47] . It must be noticed that such hot spots can in some cases be used as sintering aids ("internal susceptors") [48] . Local composition homogeneity is also mentioned to promote hot spots in other domains of microwave heating such as food defrosting [49] In FS, local composition inhomogeneity induces preferential current paths that lead to the development of hot spots (hot lines) [50] . Even small composition variations, inducing small conductivity variations are enough to trigger preferential paths (Figure 9).

Figure 9.

Example of inhomogeneous sample in FS: simulation of localisation in a zirconia specimen containing a local volume with slightly higher density, leading to 5% or 10% higher conductivity.

Preventing and using hot spots/localization Thermal and power management

As the intrinsic instability tends to develop from any source of localisation, except for very small specimens, the first strategy consists in suppressing, limiting or compensating the sources of nucleation.

The very first idea is to reduce the thermal exchanges by using thermal insulation. It was tested with some success in the case of FS [51] using different types of insulation on small cross section (1.9X3mm 2 ) dog-bone shapes specimens. The second idea is to balance the natural thermal gradient by an adapted application of current. One example is the "current path management" [52] that was tested on numerical models and experiments in the case of a flat dog-bone specimen: two connections instead of one were used at each end, so that two privileged current paths, symmetric with respect to the dog-bone axis, partially balance the natural tendency for axial concentration of the current and heat production. With the same purpose of balancing the localisation tendency of current, Molina-Molina et al. [53] proposed to use a three-phase electrical power supply and three (or more) contacts to apply a rotating electric field ("Multi-Phase Flash Sintering") that is expected to prevent privileged current path to develop in the specimen.

In the case of MWS, thermal insulation is systematically used but it does not allow preventing localisation. For this purpose, the commonly used technique is the so-called hybrid microwave heating: parts of materials with a significant ability to couple with microwave and a weak evolution of the permittivity with temperature are placed close to or around the specimen to be heated/sintered. They act as "susceptors" that both transform the MW energy into heat (then classically transmitted to the specimen) and limit the field applied to the specimen itself.

Hybrid microwave heating appears as a good compromise [28][46] , and various strategies have been developed for optimizing susceptor configurations with respect to the properties of materials to be processed [25][54] .

Process parameters management

The second strategy is to limit or delay the development of hot spots so that high enough density is achieved before the localisation phenomena reach unacceptable levels.

Very high heating rates being conditions that promote instability, the first idea is to use a controlled heating rate, as already mentioned. This is efficient in MWS when associated to hybrid heating [27][28] . In FS, controlled current ramps used instead the usual voltage control before a 10 s current control showed efficient limitation of the current spike and favourably compared to both classic FS in the case of 3YSZ [55] . With a similar purpose, a "step-wise" technique, in which a series of progressive flashes with increasing fixed current short stages is applied, was successfully applied to sintering SOFC electrolytes [56] . In FS, the fixed current step limits the dissipated power in the specimen but does not by itself stop localisation: the power limitation may help reaching a steady state but does not guarantee it. In their experiments on BaTiO3 performed by applying a constant heating rate at fixed field E M'Peko et al. [57] observed that current localisation occurs when the limiting current I is high. They plotted the observations in "processing maps", i.e., plots of experimental results in the plane (field E, current density I). The plots showed an almost constant current density limit between "safe" and "localisation", but in these experiments the furnace temperature and electric field were strongly correlated. Several papers have reported results about such plots in the case of constant furnace temperature on different materials: 3YSZ [58] , ceria [59] , BiFeO3 [60] . All these results are schematically summarised in Figure 10. As expected from previous papers no flash occurs below a critical field value Ec. Localisation is observed for high field and currents, i.e., above the solid limit curve. The density of the flashed sintered specimens is low at low values of E and increases with increasing values of E and I, leading to possible optimal processing conditions between the dashed and full curves of Figure 10. In a first order evaluation, the experimental limit should be close to constant power hyperbolas 𝐸𝐸𝐼𝐼 = 𝑃𝑃 * , the product 𝐸𝐸𝐼𝐼 being the maximum value of the power provided to the specimen during the FS experiment, at the transition from fixed field 𝐸𝐸 to fixed current 𝐼𝐼. These processing maps are useful but only partly capture the experimental conditions, as they are drawn for a given value of the furnace temperature T0 and a given hold time in the fixed current step, as pointed out by Trombin and Raj [58] . As previously noticed, higher furnace temperature is expected to decrease the thermal gradient that favours the hot spot formation, while increasing the hold time might favour the densification (and grain growth) as well as the development of the hot spot, so that complementary maps using the plane (current density, holding time) should also be used [6][58] . This kind of consideration was also used for developing a "continuous flash sintering" technique in which, after the flash phenomenon is initiated, a long specimen is moved at an optimal velocity between two electrodes [61] .

Using hot spots

Hot spots and localisation phenomena can also be considered as an opportunity, when the hotspot nucleation can be controlled. Jerby et al. [62][63] have developed the concept of localized microwave-heating intensification and identified applications such as drilling of glass, thermite reaction ignition, or local sintering of metal powder. Riminesi & Olmi [START_REF] Riminesi | [END_REF] proposed to use localized microwave heating against biotic agent growth within wooden artefacts and on the surface of stone artefacts.

Conclusion

Hot spots in FS and MWS result from the thermal runaway due the intrinsic instability of heating with an electrical or electromagnetic field when the coupling between the field and the material, i.e., the heat source, increases with increasing temperature. This is the case in most ceramics. The localisation of the thermal runaway results from the competition between local heat source and local heat evacuation. It develops as the amplification of any localised temperature or electromagnetic power excess: any localisation nucleus develops when it is larger than a critical size that depends on the material and processing conditions. This amplification phenomenon is enhanced by high field and high temperature coefficient 𝑑𝑑γ 𝑑𝑑𝑇𝑇 of the coupling parameter γ (electrical conductivity σ and imaginary part of the permittivity 𝜀𝜀 𝑒𝑒 ′′ , respectively) and by low values of the thermal conductivity. In practice, high fields are required to get high heating rates.

A major and systematic source from which the localisation develops is the inverted temperature gradient of the part inherent to these heating techniques. Any default of the specimen or the process conditions with local character can also initiate the phenomenon: local composition or field variations, inhomogeneous electrical contact can be the practical source of hot spots.

Two types of approaches have been proposed in the literature to prevent or limit hot spots. The first one consists in limiting the nucleation, e.g. by trying to suppress thermal gradients or by using geometries and field applicators that are expected to balance the thermal gradients.

The second approach is illustrated by the "processing maps", that aim at identifying process conditions in which optimal results can be obtained. These maps evidence the compromises to be found between the target properties (density, homogeneity) and the process parameters.

Although localisation is almost unavoidable, both types of approaches contribute to obtain acceptable conditions for the densification of homogeneous ceramic parts while keeping the expected advantages of field assisted techniques, namely fastness and energy saving.

Figure 1 .

 1 Figure 1. Flash sintering (FS): dog-bone and cylinder specimens (a); two steps experiments with fixed voltage V (fixed field E) before the flash event then fixed current Imax for time tm (b)evolution of temperature and current during a FS experiment calculated from Equation 3 using data for 3YSZ from Todd et al.[18] (c).

Figure 2 .

 2 Figure 2. Geometry and electric field in the two major devices used in microwave sintering: single mode (a) and multimode(b) cavities.

Figure 3 .

 3 Figure 3. Small volume δV of length λ corresponding temperature profile (red line) and local approximate temperature gradient (red dotted line)

Figure 4 .

 4 Figure 4. Hots spots during MWS: hot spot on the specimen surface observed by IR camera (a); local fracturing of a specimen during sintering due to local thermal stress stresses induced by hot spot (b); local melting of the thermal insulation around the specimen due to localised heating of the insulating material (c).[START_REF] Garnault | Frittage Par Micro-Ondes de Céramiques Mises En Forme Par Impression 3D[END_REF]

Figure 5 .

 5 Figure 5. Current flow concentration and "hot lines" in FS experiments: mostly non densified alumina-zirconia composite pellet showing densification and grain growth and even melted eutectic along the flow axis (a)[START_REF] Bichaud | Frittage « flash » de Céramiques Sous Courant Alternatif[END_REF] ; top and cross section of a LVP (Li3V2(PO4)3) pellet showing the localisation of densification and fracture (b).

Figure 6 .

 6 Figure 6. Initiation and development of localisation event: typical evolution of a representative physical entity (temperature T, local power P, current density I) (a) and of the normalised value showing the narrowing of the hot region (b) and numerical example of a FS cylindrical homogeneous specimen submitted to axial current and side cooling (c). The simulations presented in the present paper are performed using the electrical conductivity for YSZ proposed

Figure 7 .

 7 Figure 7. Examples of inhomogeneous electric field application field in MWS in a single mode cavity: field in and around an empty conical specimen (a); power density in cylindrical and conic centimetre specimens slightly off-centred from a field node (b). Simulations are performed for a 1 cm size zirconia specimen at room temperature.

Figure 8 .

 8 Figure 8. Examples of inhomogeneous electric field application field in FS: off-axis contacts in a dog-bone-like specimen (a); numerical simulation showing the development of the instability in a pellet with perfect bottom contact and localised top contact (b).

Figure 10 .

 10 Figure 10. Schematic representation of processing map in the case of fixed temperature FS experiments with 2 steps: fixed field E then fixed current I.

Acknowledgements

This work was supported by the French National Research Agency through ANR-17 CE08 0021 and ANR-20-CE05-0040-02 projects and by the Centre of Excellence of Multifunctional Architectured Materials "CEMAM".