H Huang
email: he.huang@laas.fr

A Boyer

S Ben Dhia

B Vrignon

Susceptibility Analysis of an Operational Amplifier Using On-Chip Measurement

Keywords: integrated circuit, on-chip measurement, conducted susceptibility I

This paper presents an electromagnetic immunity study of a simple operational amplifier by using the on-chip measurement. With the technology of on-chip non-invasive sensor, the internal inaccessible signal (voltage/current) can be obtained accurately in time domain. This approach grants a good insight in the internal transient response caused by the external electromagnetic interference. The validity of the on-chip measurement results is discussed comparing with the off-chip measurement and simulation.

INTRODUCTION

These last years, the susceptibility of integrated circuit (IC) has become an important issue for all circuit classes. To reduce the redesign time, the need for prediction of the risks of noncompliance during the design phase has become critical for IC manufacturers [START_REF] Ramdani | The electromagnetic compatibility of integrated circuits -past, present and future[END_REF]. By using a circuit simulator with the available models of components (passive components, ICs, parasitic parts, transmission lines …), the circuit functional failure caused by the electromagnetic interferences (EMI) could be predicted. From the simulation results, the design of circuit can be easily optimized. Although the simulation models are already very accurate, it's hard to consider and reconstruct all the elements of the real chips, especially the parasitic resistances, inductances and capacitances of interconnects, I/O pads, substrate or package. These elements play an important role in the filtering of high frequency interference. Several simulation results (e.g. I/O voltages) could be verified by off-chip measurements, like S parameters and probing testing [START_REF] Vrignon | On-chip sampling sensors for high frequency signals measurement: evolution and improvements[END_REF]. However, verifying the exactness of the simulated voltage profiles of internal inaccessible circuit nodes requires on-chip measurement methods.

As a recent technology began in the early 2000s, the onchip measurement has progressed a lot in the recent years [3] [4]. In [START_REF] Vrignon | On-chip sampling sensors for high frequency signals measurement: evolution and improvements[END_REF], the authors have demonstrated the precise results of the on-chip measurement for high frequency signals in time domain, like crosstalk, delay, and supply bounce. These measurement results could be verified by the simulation. In [START_REF] Deobarro | Use of on-chip sampling sensor to evaluate conducted RF disturbances propagated inside an integrated circuit[END_REF], the first use of on-chip sensor for the characterization of the coupling of external disturbances to IC power supply pin is presented. The comparison between the on-chip and off-chip measurements in [START_REF] Ben Dhia | On-Chip Noise Sensor for Integrated Circuit Susceptibility Investigations[END_REF] reveals the advantage of on-chip measurement in high frequency (more than 50 MHz). The differences between the two measurements is related to the power distribution network (PDN) which represents the physical interconnects of component, parasitic capacitances and the package pins.

As a common electronic block, failures of the operational amplifiers (OPA) caused by electromagnetic interference is widely discussed and analyzed by experimental measurement, simulation modeling and mathematical demonstrations [START_REF] Graffi | EMI-induced failures in integrated circuit operational amplifiers[END_REF] [8] [START_REF] Edout | EMC of Analog Integrated Circuits[END_REF]. The coupling of EMI on OPA inputs and power supply leads to a distortion of the output voltage, especially the generation of an offset. The origin of the problem lies in the distortion of the current which circulates in the differential pair. The characterization of these currents is essential to validate the exactness of a prediction model of the susceptibility to EMI, but this task remains difficult as the differential pairs are inaccessible.

This paper aims at presenting the susceptibility analysis of a simple operational amplifier based on the use of on-chip sensors. These sensors provide a precise time domain characterization of the voltage of internal nodes. It can help to understand the reaction of the amplifier to EMI, specially the currents of the differential pair, and to validate or improve the circuit model. Section II details the on-chip measurement method and the experimental set-up. Experimental results are presented and compared with the first simulation results in Section III. Finally, the conclusion and future work are provided in Section IV.

II. ON-CHIP SENSOR AND MEASUREMENT SETUP

A. On-chip sensor

The architecture and the principle of signal reconstruction of the on-chip sensor used in this test have been presented in [START_REF] Ben Dhia | On-Chip Noise Sensor for Integrated Circuit Susceptibility Investigations[END_REF] and are described in Fig. 1. An on-chip sample and hold (S/H) cell probes directly the voltage signal. The high input impedance of the probe ensures that the measurement is not invasive. The sampling command is delayed by a delay cell. The delay is set according to the test frequency, in order to get a good time resolution of the reconstructed waveform.

With this method, a very high virtual sampling frequency can be achieved without the constraints of hardware bandwidth. According to [START_REF] Vrignon | On-chip sampling sensors for high frequency signals measurement: evolution and improvements[END_REF], the bandwidth of sensor in our test (90 nm technique) is up to 14 GHz, which is sufficient for the measurements in this study (up to 2 GHz).

B. Device under test

The amplifier under test is embedded in a test chip designed in Freescale ® CMOS 90-nm process. It is a simple Miller operational amplifier with a P-channel input stage, in order to accept a common mode-voltage down to zero. As illustrated in Fig. 2, the input stage comprises 3 main parts: a differential pair of PMOS (M2 and M3), a bias current source of PMOS (M1) controlled by an external bias voltage (V bias ), and the active load (M4 and M5). The output terminal (V out ) is an amplified copy of the voltage difference between the two input terminals (V + and V -).

Two 500 Ω current-sensing resistors are placed in both branches of the differential pair and four on-chip sensors are connected to the resistor terminals. The measurement of the voltages across each resistor provides the drain current of each transistor of the differential pair by using ( 1) and [START_REF] Vrignon | On-chip sampling sensors for high frequency signals measurement: evolution and improvements[END_REF]. The sensor measurements have to be synchronized to a same time references to ensure a correct extraction of the current.

      m im im m R t V t V t i 2   (1)       p ip ip p R t V t V t i 2   (2)
The measurement repeatability of the on-chip sensor for low voltage (0 to 3.75V) is estimated to be ±4 mV, adding with the accuracy of external equipment given as ±3 mV, the accuracy of the low voltage sensor used in this device is ±7 mV. Since two voltage sensors operate simultaneously for the current measurement, the accuracy of current i m and i p is evaluated to ±28 μA (±7 mV * 2 / 500 Ω). In this case study, the current of the differential pair is defined by the bias current, the addition of two current-sensing resistors has a little influence of the function of amplifier. That is why the 500 Ω resistors are chosen which we can get a good accuracy and a little voltage drop. If the value of the added current-sensing resistor is too important to modify the normal activity of the circuit, like the supply current measurement noted in [START_REF] Vrignon | On-chip sampling sensors for high frequency signals measurement: evolution and improvements[END_REF], a small resistor like 1.5 ohm in this case could be used not to be too invasive, but on the other hand, the current measured will be less accurate.

V-V+

A fifth on-chip sensor is also installed at the output terminal of the operational amplifier. This output voltage can be measured in the PCB card by an oscilloscope and a highfrequency probe, after a transmission via package, socket, PCB track and connector, which represent an important interconnect that will filter the timing profile of high frequency signals. This off-chip measurement could be used to verify the on-chip sensor result for DC or low frequency test.

C. Experimental set-up

As described in Fig. 3, the operational amplifier is in voltage-follower configuration in our test, because this topology maximizes the susceptibility to EMI [START_REF] Masetti | Failures induced on analog integrated circuits from conveyed electromagnetic interferences: A review[END_REF]. Normally, if there is no disturbance injected in the input terminal of voltage-follower, the output will track the input voltage. A conducted EMI voltage is superimposed to a constant voltage equal to 1.

III. EXPERIMENTAL RESULTS AND ANALYSIS

A. Off-chip and on-chip measurements Among the five sensor test points, only the output voltage of amplifier could be compared with an external probe testing. EMI with a given amplitude A emi is injected at various frequencies F emi on the non-inverting input V + . As shown in Fig. 4, a negative DC offset can be observed in the output by both external and on-chip measurements. Normally the AC part of the EMI-induced distortion is not very harmful because it can be filtered easily. However, the DC offset is nearly not possible to eliminate.

Although both measurement methods give almost the same result of DC offset, the waveforms are not exactly the same. An oscillation with frequency about 300 MHz is completely filtered by the interconnector between the internal V out terminal in chip and the probe test point in the PCB card (Port 2 in Fig. 3). This phenomenon is clearly visible on the on-chip measurements.

B. On-chip current measurement

The on-chip measurement of current illustrates a good repetition, and a small dispersion between the different devices under test. Fig. 5 illustrates the two currents of the differential pair. The EMI injected in the input is the same as the case (a) of the V out measurement (Fig. 4). The current waveforms reveal the distortions induced when the EMI is injected in the input. When the input voltage is higher than V out , in the minus side the transistor is saturated, and the transistor of the plus side is blocked. The internal ringing exists also in the waveform of currents, but this oscillation is not as obvious as in the voltage measurement. The rapid peak could be related to the sensor

C. Modeling and simulation

As illustrated in Fig. 6, the modeling of system is constructed from the existing transistor net list. Besides, a first simple passive distribution network (PDN) model is built from a two-port S parameter characterization between port 1 and port 2 in the PCB card, as noted in Fig. 3. This simple PDN model fits with the measured impedance profile up to 100 MHz. On-going modeling works are done to extend the validity range of this model. All the simulations are performed with Cadence. This preliminary susceptibility simulation is limited to 100 MHz owing to the precision of PDN model. As shown in Fig. 8, the simulation demonstrates the activity saturation/blockage of differential currents, and the current level of simulation fits well with the measurement. However, these is a gap of saturation/blockage time of transistor between the simulation and the measurement, this difference may relate with the internal parasitic elements which are not modeled enough accurately in the simulation. Both the DC and susceptibility simulations validate the results of on-chip current measurement.

Vin

Figure 1 .

 1 Figure 1. Architecture and the sequential equivalent-time sampling principle of the synchronous on-chip sensor [10]

Figure 2 .

 2 Figure 2. Operational amplifier stucture and the positions of on-chip sensors

Figure 3 .

 3 Figure 3. Experimental setup for the susceptibility analysis of the amplifier

Figure 4 .

 4 Figure 4. Comparison between external probe measurement and on-chip sensor measurement of Vout at same EMI amplitude Aemi = 1V and different EMI frequency Femi: (a) 20 MHz; (b) 50 MHz; (c) 100 MHz; (d) 2 GHz

Figure 5 .

 5 Figure 5. The current of the differential pair of amplifier measured by the on-chip sensor with EMI (Femi = 20 MHz, Aemi = 1 V)

Figure 6 .Figure 7 .

 67 Figure 6. Modeling of the operational amplifier in voltage-follower configuration avec a preliminary PDN model Fig. 7 illustrates the evolution of the internal current of the minus side of the differential pair for different DC values of input where no interference is injected. The simulation verifies the current value obtained by the on-chip measurement.

Figure 8 .

 8 Figure 8. The comparison of the current in the differential pair between onchip sensor and simulation with EMI: (a) Femi = 20 MHz, Aemi = 1 V; (b) Femi = 50 MHz, Aemi = 1 V

IV.

CONCLUSION

This paper has presented an on-chip measurement method used in the susceptibility analysis of a simple operational amplifier. With the on-chip measurement, several important internal inaccessible signals could be observed, like the current flowing in a differential pair of an operational amplifier. Measuring these currents is very useful to understand the susceptibility of the circuit to electromagnetic interferences and to verify the modeling. According to the comparison with the off-chip measurement, the on-chip measurement reveals a good precision in the characterization of high frequency signals. A preliminary model of the operational amplifier has been proposed. The simulation results fit partially with on-chip measurements. Parasitic elements will be introduced in the model to extend its validity at higher


H Huang
email: he.huang@laas.fr

A Boyer

S Ben Dhia

The detection of counterfeit integrated circuit by the use of electromagnetic fingerprint

Keywords: integrated circuit, counterfeit detection, electromagnetic emission I

This paper proposes a preliminary demonstration of an alternative method of detection of counterfeit integrated devices, based on the « electromagnetic fingerprint » obtained by electromagnetic emission measurements. The principles of the methodology are explained. Two case studies are presented and two data analysis methods are discussed.

INTRODUCTION

A counterfeit electronic is a device whose material, performance, or characteristics are knowingly misrepresented by the vendor, supplier, distributor, or manufacturer [START_REF] Livingston | Avoiding Counterfeit Electronic Components[END_REF]. This use of electronic counterfeits can cause reduced performance of circuits, such as instability of clock frequency, operating life decrease, a lower storage memory space, or failure of the whole system. In recent years, there is a growing number of reported incidents related to counterfeit integrated circuits (ICs) [START_REF]Defense Industrial Base Assessment: Counterfeit Electronics[END_REF]. The most commonly counterfeited semiconductor types are analog ICs, microprocessors and memories [START_REF]Defense Industrial Base Assessment: Counterfeit Electronics[END_REF]. The Alliance for Gray Market and Counterfeit Abatement estimates that nearly 10 percent of technological products sold in the global market are counterfeit, which represents about 100 billion dollars loss for the electronic companies every year [START_REF] Pecht | Bogus: Electronic manufacturing and consumers confront a rising of counterfeit electronics[END_REF].

Three major techniques exist to produce counterfeit circuits: re-marking components as a higher grade and more expensive chip, re-packaging old devices up to non-qualified components, and duplicating counterfeits through inferior parts or materials [START_REF] Bastia | Next Generation Technologies to Combat Counterfeiting of electronic components[END_REF]. The struggle against this problem relies in a better management of the supply chain and the market [START_REF] Mcfadden | Supply Chain Risk Mitigation for IT Electronics[END_REF], the development of advanced authentication methods (e.g. watermarking [START_REF] Koushanfar | Can EDA combat the rise of electronic counterfeiting?[END_REF]) and serialization technologies (e.g. 2-D bar code [START_REF] Bastia | Next Generation Technologies to Combat Counterfeiting of electronic components[END_REF]). Besides, industry always looks for non-destructive, rapid and cheap tests. A large number of detection techniques exist and are already used by industry and detection laboratories: visual inspection of the appearance of component (texture, mold mark, pitch form…) in order to detect raw counterfeit devices, decapsulation to verify the die layout, material analysis methods, like fluorescent X-ray or C-mode scanning acoustic microscope for package analysis, and different levels of electrical testing (V/I characteristics, ESD test, operation life test…) [START_REF] Hewett | Methods used in the detection of counterfeit electronic components[END_REF]. The electrical tests are adapted to detect elaborated counterfeit devices. They consist in obtaining electrical signals of IC pads such as voltage/current/impedance profile and comparing the results to expected figures given by datasheets, or verifying the response to given vectors. However, this type of methods has some drawbacks. Firstly, the characterization of IC pads provides only information about the status of interfaces. Secondly, a large number of vectors and a good knowledge of circuit operation are required to perform electrical test.

Although there are lots of methods to detect counterfeit devices, most of them are destructive, and none can cover 100% of counterfeit types. Moreover, the detection methods have to evolve because the counterfeiting techniques may adapt to them. Furthermore, we need a large numbers of detection techniques to deal with the large number of counterfeiting possibilities.

The electromagnetic emission (EME), also called "electromagnetic fingerprint" in this article, is a contactless side channel related to the IC transient activity. It depends on numerous circuit parameters such as technology, placement and routing, embedded code, internal filtering, packaging, temperature, aging [8] [9]... Any modification of one of these parameters can lead to a significant change of the electromagnetic emission. This principle emerges as a new idea to detect counterfeit ICs [START_REF] Gross | Detecting counterfeit electronic components using EMI telemetric fingerprints[END_REF] [START_REF] Keller | System and method for physically detecting counterfeit electronics[END_REF]. However, no demonstration of the feasibility of the method has been shown in literature. This paper aims at presenting a first study about the feasibility to distinguish authentic and different devices through electromagnetic emission measurements.

II. DESCRIPTION OF THE DETECTION METHODOLOGY

A. Principles of the method

A counterfeit component or circuit exhibits a different or degraded operation compared to the nominal behavior of the original device. The electromagnetic (EM) fingerprint is a parasitic electromagnetic signal produced by IC internal activity, measured in conducted and/or radiated mode, in given experimental conditions. The device under test could be a small active component, a complicated electronic system or whatever electronic active device which can generate the electromagnetic emission.

The detection test consists in comparing the EM fingerprint measured from a suspect test device and a reference fingerprint from good devices. The devices under test must be powered and set in a given configuration which induces at least a transient activity or current consumption. If all the measurements are done in similar conditions, any significant difference of the electromagnetic emission will indicate a difference between the devices. However, only one known good device is not enough to extract a reference fingerprint due to measurement errors and process dispersion. It is better that the reference fingerprint could include statistical information about the EM fingerprint of authentic devices, which can be extracted by the measurements of a number of known good devices.

B. Electromagnetic measurements

The electromagnetic noise produced by the circuit activity can be coupled according to three modes: conducted mode, farfield radiated mode or near-field mode. The two last modes ensure contactless measurement. In addition, the near-field mode allows localization of the source of electromagnetic [START_REF] Slattery | Measuring the electric and magnetic near fields in VLSI devices[END_REF]. Proven methods for characterizing electromagnetic emissions from ICs between 150 KHz to 1 GHz are proposed by IEC 61967 standard [START_REF]Integrated circuits -Measurement of electromagnetic emissions, 150 kHz to 1 GHz[END_REF], which is extensively used by IC manufacturers for electromagnetic compatibility qualification. The concepts, requirements and advices given by this standard can be followed. Besides, all the EME measurements do not damage the device under test, and the type of measurement could be chosen according to the device and the requirement of detecting application.

The measurement can be done in time domain with oscilloscopes or frequency domain with narrowband receivers such as spectrum analyzer (e.g. in this study a MS2667C portable spectrum analyzer of Anritsu is used). The measurement in the frequency domain outweighs the one in time domain by having a lower noise floor. Moreover, as transient current induced by circuit activity has an impulse nature, its spectrum covers a large frequency range. Therefore, in our two case studies more differences can be detected in frequency domain than in time domain. No matter which measurement method is chosen, the same measurement conditions for the different devices should be guaranteed.

C. Analysis method

The analysis of the emission spectrum measured in frequency domain is discussed in this paper. The amplitude measured at a finite number of frequencies M are stored in a vector E(k), k = [1;M]. If the circuit activity is periodic, the amplitude can be extracted at harmonic frequencies. The reference fingerprint is noted E REF (k) and the fingerprint of the suspect device under test (DUT) is noted E DUT (k). The reference fingerprint E REF (k) contains the average emission level measured at frequency k over the N reference samples. Only the points whose amplitude is larger than the noise floor are taken into account to minimize measurement errors. According to the work of [START_REF]These de binghong[END_REF], the distribution of their emission level around this average level is supposed to be normal. The dispersion around the average level of E REF (k) in frequency k linked to measurement uncertainties and process differences is given in term of standard deviation σ REF (k), which can be obtained by the measurements of N reference samples of certain times (e.g. 5 times per reference sample).

Two basic statistical criteria are proposed for the analysis of the EM fingerprint between the reference and suspect devices. These two statistical methodologies are very simple and common for estimating the relationships among variables. The first criterion is called z-score (also called standard score) [START_REF] Livingston | Avoiding Counterfeit Electronic Components[END_REF]. For the different frequencies, this estimator gives the difference between the suspect device and reference fingerprints, divided by the standard deviation of the reference fingerprint σ REF .

For each frequency k, this figure provides an indication about the probability of differences between the suspect device and reference emission levels, so a z-core in frequency k close to 0 means that the emission levels of the two devices are similar in this frequency.

( ) ( ) ( ) ( ) k k E k E k z REF REF DUT σ - = (1)
The second and alternative criterion is the determination coefficient R² [START_REF]Defense Industrial Base Assessment: Counterfeit Electronics[END_REF]. In regression analysis, this topic relates to goodness of fit. It is computed over a group of frequencies (from frequency i to frequency j). This coefficient determines whether a linear relation exists between suspect device and reference device fingerprints. The interval of the value of R 2 is [0, 1], a determination coefficient close to 1 indicates a strong linear relation between suspect device and reference fingerprints, and vice versa.

( ) ( ) ( ) ( ) ( ) ( ) 2 : : 2
: , : :         ⋅ = j i E j i E DUT REF REF REF j i E j i E Cov j i R σ σ (2)
Where, Cov is the covariance of the reference vector and the DUT vector, and σ is the standard deviation of these two vectors respectively. Like all the other statistical methodologies, these two criteria that we choose in this study cannot represent all the relationships between the DUT and the reference, but as the preliminary analysis, they are effectual to identifier the differences related to the counterfeits.

III. CASY STUDY I : DISTINCTION BETWEEN COMPONENTS WITH TECHNOLOGICAL DIFFERENCES

A. Case study presentation

As buying certified counterfeit circuits is uncertain, the proposed study is done on "simulated counterfeit devices", i.e. with known small technological differences or stressed devices representing two large kinds of counterfeits in the electronic product market.

A digital circuit has been designed in CMOS 0.25 µm process. As resumed in Table I, the reference electromagnetic fingerprint is obtained from a group of 9 authentic circuits designed in the same technology. The reference devices are called Core0(ref). Nine other components are set to be the suspect components under test. One component called Auth0 is designed with the same technology as the device of the reference group and is randomly chosen to verify whether its EM fingerprint can be identified as similar to the reference fingerprint. Eight other components under test used a similar technology which differs from the amount of distributed capacitor and substrate isolation. They are called Core1 (from Core1-1 to Core1-8). Their electromagnetic emissions are measured and compared to the reference fingerprint in order to verify whether they can be separated from the reference devices.

B. Comparison of emission level

Conducted emission test have been chosen in this case study. The transient current that returns to the ground is measured by the use of a 1 Ω resistor probe, as the 1 Ω direct coupling method detailed in the standard 61967-4 [START_REF]Integrated circuits -Measurement of electromagnetic emissions, 150 kHz to 1 GHz[END_REF]. Since the emission becomes negligible above 300 MHz, the measurements are performed from 150 KHz up to 300 MHz. The fundamental frequency of emission is 4 MHz, Fig. 2 compares the reference fingerprint with the EM emission of one Core1 and Auth0 devices. The general shapes of the three spectra are similar and, in spite of small differences between these different devices, it is difficult to conclude about the differences between the three samples. Only statistical analysis can provide a more precise conclusion by revealing differences or correlations between EM fingerprints.

C. Analysis of the differences

The two statistical criteria ( 1) and ( 2) are applied on the different suspect samples with the reference devices. Fig. 3 compares the z-score of one Core1 sample and Auth0. The result reveals that a huge difference exists between the fingerprint of Core1-3 and Core0 over all the considered frequency range. The highest z-score reaches more than 60 times the standard deviation of the reference fingerprint S REF , and the average difference is 8.75 times S REF . The probability for this Core1 device to have a similar EM fingerprint as the reference devices is very small. In contrast, the difference between Auth0 fingerprint and the reference fingerprint does not exceed 2 times S REF and its average value is about one standard deviation. It is unlikely that Auth0 and reference devices have different EM fingerprints. Fig. 4 presents a scatterplot between the emission levels of the reference and the tested devices (Auth0 and Core1-3). A clear linear relation exists between Auth0 and reference group emission levels. This is underlined by a determination coefficient very close to 1. The small differences are due to measurement errors and process dispersion. In contrast, the linear relation between Core1 and reference levels is not so obvious, as demonstrated by the lower value of R². Analyses on the other devices provide similar results, as shown in Fig. 5. The average z-scores and determination coefficients measured with each Core1 devices show large differences with reference devices, which indicate that their EM fingerprints are not correlated to the reference fingerprint due to their internal differences. Auth0 EM fingerprint presents similar statistical properties to the reference fingerprint. This result reveals that design differences such as those brought in Core1 induces significant changes in the EM fingerprint. These differences can be detected from electromagnetic measurements between circuits.

A. Case study presentation

Recent publications have shown that IC aging provokes a change of electromagnetic emission profile [START_REF] Ben Dhia | Characterization of the Electromagnetic Modelling drifts of a nanoscale IC after Accelerated Life Tests[END_REF]. Degradation mechanisms accelerated by harsh conditions such as Negative Bias Temperature Instability or Hot Carrier Injection induce a spread of the transient current due to IC activity. One source of counterfeit devices is related to the reuse of non-qualified and recycled old devices. This part aims at demonstrating that the comparison of the EM fingerprints of a group of reference components with a same but aged device can provide a detection method.

In this case study, Core0 components are reused, except that they have been submitted to stress conditions according to High Temperature Operating Life (HTOL) test [15]. During 408 hours, the samples have been powered under high temperature conditions (150°c). Before and after the accelerated-life test, the same conducted electromagnetic emission test of the previous case study is employed. The reference group is constituted by the fresh components. After aging the components are still functional and a little variation of the emission level is observed, as shown in Fig. 6.

B. Analysis of differences

The two statistical criteria defined in II.C have been applied on the measured emission results. Fig. 7 presents the summary of statistical analysis done on the fresh and aged Core0 samples. The average z-score shows that the average emission level differences between aged components and reference group are at least two times larger than the fresh samples. Moreover, the linearity of the relation between the reference fingerprint and the devices fingerprint of tends to degrade after the accelerated aging. Although these differences are less significant than these in the last case study, the result of this study reveals that the aging of circuit can induce detectable changes in the EM fingerprint.

V. CONCLUSION

This paper has proposed a preliminary demonstration of an alternative method dedicated to the analysis of traceability of integrated devices, based on a measurement of the « electromagnetic fingerprint » of a circuit. The presented results have shown that significant changes in the electromagnetic fingerprint of a circuit can be measured when some modifications of the original design are provided or when it is submitted to an accelerated-aging test. Conducted measurements have been leaded to extract the electromagnetic fingerprint in the case studies of this paper, and several parallel tests applied on other circuits with radiated or near-field measurement methods have provided positive results. More studies are required to evaluate the robustness of the method and to compare the advantages and drawbacks with other detection techniques. Also, improvements have to be made to increase the speed of measurements (parallel tests, rapid analyzer, compromise between measurement accuracy, noise floor and measurement time).

Fig. 1 Figure 1 .

 11 Figure 1. Description of the detection methodology

Figure 2 .

 2 Figure 2. Comparison between reference (Core0), Auth0 and Core1-3 fingerprints

Figure 3 .

 3 Figure 3. Comparison of z-score between Auth0 and Core1-3 fingerprint

Figure 4 .

 4 Figure 4. Scatterplot of Auth0 (a) and Core1-3 (b) emission levels vs.reference emission level

Figure 5 .

 5 Figure 5. Comparison of statistical analysis results between reference, Auth0 and Core1 devices: Average Z-score (a), Determination coefficient (b)

Figure 6 .

 6 Figure 6. Comparison between reference (Core0), Core0-1 (fresh) and Core0-1 (aged) fingerprints

Figure 7 .

 7 Figure 7. Comparison of statistical analysis results between fresh and aged Core0 devices: Average Z-score (a), Determination coefficient (b)

TABLE I .

 I IDENTIFICATION OF DEVICES OF CASY STUDY I

	Group name	Reference group		Suspect group
	Identity	Core0 (ref)	Auth0	Core1-1 ~ Core1-8
	Number	9	1	8
	Technology	Technology Core0	Technology Core1


Francis Filbet
email: filbet@math.univ-lyon1.fr

Chang Yang
email: yangchang@hit.edu.cn

MIXED SEMI-LAGRANGIAN/FINITE DIFFERENCE METHODS FOR PLASMA SIMULATIONS

Keywords: Cartesian mesh, semi-Lagrangian method, Hermite WENO reconstruction, guidingcenter, drift-kinetic model. 2000 Mathematics Subject Classication. 65M08, 65M25, 78A35

In this paper, we present an efficient algorithm for the long time behavior of plasma simulations. We will focus on 4D drift-kinetic model, where the plasma's motion occurs in the plane perpendicular to the magnetic field and can be governed by the 2D guiding-center model.

Hermite WENO reconstructions, already proposed in [25], are applied for solving the Vlasov equation. Here we consider an arbitrary computational domain with an appropriate numerical method for the treatment of boundary conditions.

Then we apply this algorithm for plasma turbulence simulations. We first solve the 2D guiding-center model in a D-shape domain and investigate the numerical stability of the steady state. Then, the 4D drift-kinetic model is studied with a mixed method, i.e. the semi-Lagrangian method in linear phase and finite difference method during the nonlinear phase. Numerical results show that the mixed method is efficient and accurate in linear phase and it is much stable during the nonlinear phase. Moreover, in practice it has better conservation properties.

Introduction

In the context of strongly magnetized plasma simulations, the motion of particles is confined around the magnetic field lines; the frequency of this cyclotron motion is faster than the frequencies of interest. Therefore, the physical system can be reduced from the 6D Vlasov-Maxwell system to a four or five dimensional model by averaging over the gyroradius of charged particles (See for a review [START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF][START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF]). In this paper we focus on 4D drift-kinetic model, where the movement of the plasma in the plane perpendicular to the magnetic field can be governed by the guiding-center model. This reduced model could help us to investigate plasma turbulence problems with an acceptable computational time. More especially, using the 2D guiding-center model, we can focus directly on the difficulties of treatment of boundary conditions on arbitrary computational domain. Moreover, these reduced models have the conservative properties, which can be used as criterion to evaluate the good performance of numerical methods.

To develop accurate and stable numerical techniques for plasma turbulence (4D drift kinetic, 5D gyrokinetic and 6D kinetic models) is one of our objectives. In [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-Lagrangian code[END_REF] several numerical solvers have been developed using an Eulerian formulation for gyro-kinetic models. However, spurious oscillations often appear in the nonlinear phase when small structures occur and it is difficult to distinguish physical and numerical oscillations. Moreover, for these models semi-Lagrangian methods are no more conservative, hence the long time behavior of the numerical solution may become unsuitable. At contrast, a class of numerical methods based on the Hermite interpolation [START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF], together with a weighted essentially non-oscillatory (WENO) reconstruction can be applied either to semi-Lagrangian or to finite difference methods [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially-non-oscillatory reconstruction for Vlasov equations[END_REF].

On the one hand, the semi-Lagrangian methods are very efficient and fast but do not conserve mass in an arbitrary grid. On the other hand, the finite difference methods are mass conservative and more stable for long term simulations but have a restrictive CFL constraint. We thus propose a mixed method to carry forward the advantages of each method, i.e. the semi-Lagrangian method in linear phase where the solution is relatively smooth, the finite difference method during the nonlinear phase where a lot of small structures appear. We will apply the mixed method to the 4D drift-kinetic model to verify its efficiency.

The numerical resolution of 4D drift-kinetic model in a cylinder has been already performed via polar coordinates [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]. However, polar coordinates are not always suitable since artificial singular points appear in the computational domain coming from the change of variable. Moreover, for more complicated geometry, it is not straightforward to perform the appropriate change of variables. From these considerations, we focus here on the discretization of transport models on a Cartesian mesh and develop a suitable technique as in [START_REF] Filbet | An inverse Lax-Wendroff method for boundary conditions of Boltzmann equations[END_REF] to treat boundary conditions. This method is based on the inverse Lax-Wendroff method [START_REF] Filbet | An inverse Lax-Wendroff method for boundary conditions of Boltzmann equations[END_REF]. On the other hand, to compute the electric potential from the Poisson equation, we apply an extrapolation technique [START_REF] Filbet | Numerical Simulations of Kinetic Models for Chemotaxis[END_REF]. To test the numerical algorithm, the guiding center model will be solved on a D-shape domain [START_REF] Miller | Noncircular, finite aspect ratio, local equilibrium model[END_REF]. We will first compute a steady state solution by solving numerically a nonlinear Poisson equation and then investigate its stability properties.

The paper is organized as follows : in Section 2, we present the derivation of the 4D drift-kinetic and the 2D guiding-center models, and their conservative properties. Then in Section 3, we recall the Hermite WENO reconstructions developed in [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially-non-oscillatory reconstruction for Vlasov equations[END_REF] for solving the Vlasov equations, and the treatment of boundary conditions corresponding to the Poisson equation. In Section 4, we numerically compute a steady state solution for the 2D guidingcenter model in a D-shape domain. Then, we study its stability properties by performing numerical simulations with a perturbed steady state as initial data. Finally, we numerically solve the 4D drift-kinetic model with the mixed methods, and present the numerical results of the conservative properties of the 4D drift-kinetic model and the evolution of the distribution function of density.

Mathematical models

The Vlasov equation for the distribution function f in standard form in standard notation is

(2.1) ∂f ∂t + v ⋅ ∇ x f + e m E + v × B c = 0,
where t ∈ R + is the time variable, x ∈ Ω ⊂ R 3 is the space variable, v ∈ R 3 is the velocity variable, m is the particle mass, e is its charge, E is the electric field and B is the magnetic field. We assume the electric field is computed by E = -∇φ, where φ is electric potential whereas the magnetic field is uniform B = B e z , where e z stands for the unit vector in the toroidal direction. Moreover, we assume that f is vanishing at infinity of velocity field and periodic boundary condition is taken in z direction.

To derive the drift kinetic model, we start to perform a change of variable according to the drift direction E × B; that is

(2.2) w = w w = v -U ε v , with U = c E × B B 2 .
Next, we decompose E into components along B and perpendicular to B :

E = E + E e z .
Substituting this expression in (2.2), it yields

E + v × B c = E e z + ε c w × B.
Then we introduce a new distribution function g, such that

g(t, x, w) = f (t, x, v)
for which we get that it is solution to

∂g ∂t + (U + εw ) ⋅ ∇ x g + w ∂ z g + 1 ε U T ∇ x U ∇ w g + e m E ∂ w g + w T ∇ x U ∇ w g + e cm w × B ⋅ ∇ w g = 0.
Finally, we integrate the previous equation in velocity field w ∈ R 2

∂ f ε ∂t + U ⋅ ∇ x f ε + w ∂ z f ε + e m E ∂ w f ε = -ε R 2 w ⋅ ∇ x g dw ,
where f ε = ∫ R 2 g dw . By passing formally to the limit ε → 0, we obtain the drift-kinetic model

(2.3) ∂ f ∂t + U ⋅ ∇ x f + w ∂ z f + e m E ∂ w f = 0.
On the other hand, the self-consistent potential φ is solution to the quasi-neutrality equation [START_REF] Crouseilles | An Isogeometric Analysis approach for the study of the gyrokinetic quasi-neutrality equation[END_REF] (2.4)

-∇ ⋅ ρ 0 Bω c ∇ φ + e ρ 0 T e (φ -φ) = ρ -ρ 0 ,
where ω c = eB m i is the ion cyclotron frequency, and T e and ρ 0 are the electron temperature and density profiles respectively which only depend on x . The ion density profile is given by

ρ(x, t) = R f (x, v , t)dv
and φ represents the average on the magnetic field lines, that is,

φ = 1 L z Lz 0 φdz,
with L z the length in the z variable. Dropping the ¯over the various quantities and replacing w by v, the non-dimensional form of the drift-kinetic model can be written as

(2.5) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂f ∂t + U ⋅ ∇ x f + v ∂ z f + E ∂ v f = 0, U = E × B B 2 , E = -∇φ, -∇ ⋅ ρ 0 (x ) B ∇ φ + ρ 0 (x ) T e (x ) (φ -φ) = ρ -ρ 0 .
The following proposition shows some properties of the drift-kinetic model when we ignore the effect of boundary conditions: Proposition 2.1. Let us consider Ω = R 3 or the torus. Then the drift-kinetic model (2.5) verifies the following properties :

(1) If f is smooth, we have the maximum principle

0 ≤ f (t, x, v ) ≤ max x,v (f (0, x, v )), t ≥ 0.
(

) L p norm conservation, for 1 ≤ p ≤ ∞ d dt R Ω f (t, x, v ) p dx dv = 0, t ≥ 0. 2
(3) Kinetic entropy conservation

d dt R Ω f ln f dx dv = 0, t ≥ 0.
(4) Energy conservation

(2.6) d dt R Ω (f -f M )v 2 dxdv + Ω φ(ρ -ρ 0 )dx = 0.
For practical applications, this model has to be supplemented with suitable boundary conditions when considering a domain of the form

Ω = {(x, y, z) ∈ R 3 ∶ (x, y) ∈ D, 0 ≤ z ≤ L z },
with D a two dimensional domain. We assume that the electric potential is vanishing at the boundary ∂D

(2.7) φ(x) = 0, x ∈ ∂D × [0, L z ],
and the distribution function is given by

(2.8) f (x , z, v ) = f M (x , z, v ), x ∈ ∂D × [0, L z ], with U(x) ⋅ n x ≥ 0,
where f M is a stationary solution to the drift-kinetic model (2.5). Furthermore, periodic boundary conditions are assumed for the potential and the distribution function in the z direction

φ(x , 0) = φ(x , L z ), f (x , 0, v ) = f (x , L z , v ), x ∈ D, v ∈ R.
Finally, we can derive the guiding centre model from (2.5) by integrating the equation with respect to (z, v ). We get that the reduced density ρ ∶ x ↦ R is solution to the guiding centre system of equations (2.9)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ ρ ∂t + U ⋅ ∇ x ρ = 0, U = E × B B 2 , E = -∇φ, -∇ ⋅ ρ 0 (x ) B ∇ φ = ρ -ρ 0 .
In this case, the solution ρ satisfies the following properties Proposition 2.2. Let us consider the two dimensional domain D = R 2 or the torus. Then the guiding centre model (2.9) verifies the following properties :

(1) If ρ is smooth, we have the maximum principle

0 ≤ ρ(t, x ) ≤ max x (ρ(0, x ).
(

) L p norm conservation, for 1 ≤ p ≤ ∞ d dt D ρ(t, x ) p dx = 0. 2
(3) Energy conservation

(2.10) d dt D ρ φ dx = 0.
For practical applications, we assume that the electric potential is vanishing at the boundary ∂D (2.11) φ(t, x ) = 0, x ∈ ∂D, t ≥ 0.

Numerical schemes

In this section we present the hybrid method based on the Hermite WENO reconstruction already proposed in [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially-non-oscillatory reconstruction for Vlasov equations[END_REF]. On the one hand, we apply a semi-Lagrangian method for a general transport equation written in a non conservative form. On the other hand, we apply a finite difference method, which enforces the conservation of mass when the equation is written in the conservative form. These methods are coupled with the inverse Lax-Wendroff procedure to discretize accurately boundary conditions in an arbitrary 2D geometry. Finally in the subsection 3.3, we discretize the Poisson equation for the electrical potential (2.9).

3.1.

Hermite WENO reconstruction for semi-Lagrangian methods. We briefly remind the high order Hermite interpolation coupled with a weight essentially non-oscillatory (HWENO) reconstruction for semi-Lagrangian methods. The semi-Lagrangian method becomes a classical method for the numerical solution of the Vlasov equation because of its high accuracy and its small dissipation [START_REF] Cheng | The integration of the Vlasov equation in configuration space[END_REF][START_REF] Sonnendrücker | The semi-Lagrangian method for the numerical resolution of Vlasov equation[END_REF]. For a given s ∈ R + , the differential system

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ dX dt = A(t, X), X(s) = x,
is associated to the transport equation

(3.1) ∂f ∂t + A(t, X) ∇ X f = 0.
We denote its solution by X(t; s, x). The backward semi-Lagrangian method is decomposed into two steps for computing the function f n+1 at time t n+1 from the function f n at time t n :

(1) For each mesh point x i of phase space, compute the backward characteristic X(t n ; t n+1 , x i ), the value of the characteristic at time t n who is equal to x i at time t n+1 . (2) As the function f of transport equation verifies

f n+1 (x i ) = f n (X(t n ; t n+1 , x i )),
we obtain the value of f n+1 (x i ) by computing f n (X(t n ; t n+1 , x i )) by interpolation, since X(t n ; t n+1 , x i ) is not usually a mesh point.

We apply a third order Hermite interpolation coupled with a weighted essentially nonoscillatory procedure, such that it is accurate for smooth solutions and it removes spurious oscillations around discontinuities or high frequencies which cannot be solved on a fixed mesh. Consider a uniform mesh (x i ) i of the computational domain and assume that the values of the distribution function (f i ) i and its derivative (f ′ i ) i are known at the grid points. We define two quadratic polynomials in the interval I i :

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ h l (x) = f i + f i+1 -f i ∆x (x -x i ) + (f i+1 -f i ) -∆xf ′ i ∆x 2 (x -x i )(x -x i+1 ), h r (x) = f i + f i+1 -f i ∆x (x -x i ) + ∆xf ′ i+1 -(f i+1 -f i ) ∆x 2 (x -x i )(x -x i+1 ).
The polynomial of degree 2 h l verifies

h l (x i ) = f i , h l (x i+1 ) = f i+1 , h ′ l (x i ) = f ′ i , while h r verifies h r (x i ) = f i , h r (x i+1 ) = f i+1 , h ′ r (x i+1 ) = f ′ i+1 .
The idea of WENO reconstruction is now to use the cubic polynomial when function f is smooth, otherwise, we use the less oscillatory polynomial of degree 2 between h l or h r . Thus, let us define H 3 as follows

H 3 (x) = w l (x) h l (x) + w r (x) h r (x),
where w l and w r are WENO weights. To determine these WENO weights, we follow the strategy given in [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF] and first define smoothness indicators by integration of the first and second derivatives of h l and h r on the interval I i :

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ β l = x i+1 x i ∆x(h ′ l ) 2 + ∆x 3 (h ′′ l ) 2 dx = (f i -f i+1 ) 2 + 13 3 ((f i+1 -f i ) -∆xf ′ i ) 2 , β r = x i+1 x i ∆x(h ′ r ) 2 + ∆x 3 (h ′′ r ) 2 dx = (f i -f i+1 ) 2 + 13 3 (∆xf ′ i+1 -(f i+1 -f i )) 2 .
Then we set w l and w r as

w l (x) = α l (x) α l (x) + α r (x) and w r (x) = 1 -w l (x),
where

α l (x) = c l (x) (ε + β l ) 2 and α r (x) = c r (x) (ε + β r ) 2 .
where c l = (x i+1 -x) ∆x, c r = 1 -c l and ε = 10 -6 to avoid the denominator to be zero.

Observe that when the function f is smooth, the difference between β l and β r becomes small and the weights w l (x) ≈ c l (x) and w r (x) ≈ c r (x). Otherwise, when the smoothness indicator β s , s = l, r blows-up, then the parameter α s and the weight w s goes to zero, which yields

(3.2) w l (x) ≈ 1, w r (x) ≈ 0 or w l (x) ≈ 0, w r (x) ≈ 1.
Finally, let us mention that the following fourth-order centred finite difference formula is used to approximate the first derivative at the grid point x i

(3.3) f ′ i = 1 12∆x (f i-2 -8f i-1 + 8f i+1 -f i+2 ).

Hermite WENO reconstruction for conservative finite difference methods.

When the velocity A is not constant in (3.1), the semi-Lagrangian method is not conservative even when divA = 0, hence mass is no longer conserved and the long time behavior of the numerical solution may be wrong even for small time steps. Therefore, high order conservative methods may be more appropriate even if they are restricted by a CFL condition.

In this section, we extend Hermite WENO reconstruction for computing numerical flux of finite difference method. Suppose that {f i } 1≤i≤N is approximation of f (x i ). We look for { fi+1 2 } 0≤i≤N such that the flux difference approximates the derivative f ′ (x) to k-th order accuracy :

fi+1 2 -fi-1 2 ∆x = f ′ (x) + O(∆x k ).
To approximate the flux fi+1 2 , we define a piecewise polynomial G such that it is exactly known on a set of points

x = x k+1 2 , k = i -l, . . . , i + r. G(x k+1 2 ) = G k+1 2 = ∆x k j=-∞ f j , i -l ≤ k ≤ i + r.
Thus, given the point values {f i }, we can compute G(x) by an interpolation method and therefore deduce the numerical flux by

(3.1) fi+1 2 = dG dx x=x i+1 2 .
Now to interpolate the function G(x), we apply a high order Hermite WENO scheme and outline the procedure of reconstruction only for the fifth order accuracy case.

The aim is to construct an approximation of the flux f - i+1 2 by the Hermite polynomial of degree five together with a WENO reconstruction from point values {f i } :

(1) We construct the Hermite polynomial H 5 such that

H 5 (x i+j+1 2 ) = G i+j+1 2 , j = -2, -1, 0, 1, H ′ 5 (x i+j+1 2 ) = G ′ i+j+1 2 , j = -2, 1, (2)
We construct cubic reconstruction polynomials H l (x), H c (x), H r (x) such that :

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ H l (x i+j+1 2 ) = G i+j+1 2 , j = -2, -1, 0, H ′ l (x i-3 2 ) = G ′ i-3 2 , H c (x i+j+1 2 ) = G i+j+1 2 , j = -2, -1, 0, 1, H r (x i+j+1 2 ) = G i+j+1 2 , j = -1, 0, 1, H ′ r (x i+3 2 ) = G ′ i+3 2 ,
where G ′ i+1 2 is the sixth order centered approximation of first derivative

G ′ i+1 2 = 1 60 [(u i+3 + u i-2 ) -8(u i+2 + u i-1 ) + 37(u i+1 + u i )].
Let us denote by h l (x), h c (x), h r (x), h 5 (x) the first derivatives of H l (x), H c (x), H r (x), H 5 (x) respectively. By evaluating h l (x), h c (x), h r (x), h 5 (x) at x = x i+1 2 , we obtain

h 5 (x i+1 2 ) = -8f i-1 + 19f i + 19f i+1 + 3G ′ i-3 2 -6G ′ i+3 2 27 and ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ h l (x i+1 2 ) = -2 f i-1 + 2f i + G ′ i-3 2 , h c (x i+1 2 ) = -f i-1 + 5 f i + 2 f i+1 6 , h r (x i+1 2 ) = f i + 5 f i+1 -2 G ′ i+3 2
4 .

(3) We evaluate the smoothness indicators β l , β c , β r , which measure the smoothness of

h l (x), h c (x), h r (x) on the cell [x i , x i+1 ] β l = x i+1 x i ∆x(h ′ l (x)) 2 + ∆x 3 (h ′′ l (x)) 2 dx = l 2 1 + 3 l 1 l 2 + 75 16 l 2 2 , with l 1 = u i -u i-1 , l 2 = -3u i-1 + u i + 2G ′ i-3 2 , β c = x i+1 x i ∆x(h ′ c (x)) 2 + ∆x 3 (h ′′ c (x)) 2 dx = c 2 1 + 2 c 1 c 2 + 25 12 c 2 2 , with c 1 = u i -u i-1 , c 2 = u i-1 -2u i + u i+1 , β r = x i+1 x i ∆x(h ′ r (x)) 2 + ∆x 3 (h ′′ r (x)) 2 dx = r 2 1 + 39 16 r 2 2 , with r 1 = u i+1 -u i , r 2 = u i -3u i+1 + 2G ′ i+3 2 .
(4) We compute the nonlinear weights based on the smoothness indicators

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ w l = α l α l + α c + α r , α l = c l (ε + β l ) 2 , w c = α c α l + α c + α r , α c = c c (ε + β c ) 2 , w r = α r α l + α c + α r , α r = c r (ε + β r ) 2 ,
where the coefficients c l = 1 9, c c = 4 9, c r = 4 9 are chosen to get fifth order accuracy for smooth solutions and the parameter ε = 10 -6 avoids the blow-up of α k , k = {l, c, r}.

(5) The flux f - i+1 2 is then computed as

f - i+1 2 = w l h l (x i+1 2 ) + w c h c (x i+1 2 ) + w r h r (x i+1 2
). The reconstruction to f + i+1 2 is mirror symmetric with respect to x i+1 2 of the above procedure. To discretize the Laplacian operator ∆ x φ near the physical boundary, some points of the usual five points finite difference formula can be located outside of interior domain. For instance, Figure 1 illustrates the discretization stencil for ∆ x φ at the point (x i , y j ). We notice that the point x g = (x i , y j-1 ) is located outside of interior domain. Let us denote the approximation of φ at the point x g by φ i,j-1 . Thus φ i,j-1 should be extrapolated from the interior domain.

We extrapolate φ i,j-1 on the normal direction n

(3.2) φ i,j-1 = wp φ(x p ) + wh φ(x h ) + w2h φ(x 2h ),
where x p is the cross point of the normal n and the physical boundary D. The points x h and x 2h are equal spacing on the normal n, i.e. h = x p -x h = x h -x 2h , with h = min(∆x, ∆y), ∆x, ∆y are the space steps in the directions x and y respectively. Moreover, wp , wh , w2h are the extrapolation weights depending on the position of x g , x p , x h and x 2h . In (3.2), φ(x p ) is given by the boundary condition (2.11), whereas φ(x h ), φ(x 2h ) should be determined by interpolation.

x

y n ◾ ◾ ◾ ◾ ◾ ◾ ◾ • x g • • • • • x p P * 0 ◆ • • • • • P * 1 ◆ • • • • • ◆ P * 2 ⊡ i -2 i -1 i i + 1 i + 2 j -2 j -1 j j + 1 j + 2 ⊙ x h ⊙ x 2h
Figure 1. Spatially two-dimensional Cartesian mesh. • is interior point, ◾ is ghost point, is the point at the boundary, is the point for extrapolation, the dashed line is the boundary.

For this, we first construct an interpolation stencil E, composed of grid points of D. For instance, in Figure 1, the inward normal n intersects the grid lines y = y j , y j+1 , y j+2 at points P * 0 , P * 1 , P * 2 . Then we choose the three nearest points of the cross point P * l , l = 0, 1, 2, in each line, i.e. marked by a large circle. From these nine points, we can build a Lagrange polynomial q 2 (x) ∈ Q 2 (R 2 ). Therefore, we evaluate the polynomial q 2 (x) at x h and x 2h , i.e.

φ(x h ) = 8 =0 w h, φ(x ), φ(x 2h ) = 8 =0 w 2h, φ(x ),
with x ∈ E. We thus have that φ i,j-1 is approximated from the interior domain.

However, in some cases, we can not find a stencil of nine interior points. For instance, when the interior domain has small acute angle sharp, the normal n can not have three cross points P * l , l = 0, 1, 2 in interior domain, or we can not have three nearest points of the cross point P * l , l = 0, 1, 2, in each line. In this case, we alternatively use a first degree polynomial q 1 (x) with a four points stencil or even a zero degree polynomial q 0 (x) with an one point stencil. We can similarly construct the four points stencil or the one point stencil as above.

Numerical simulations

In this section, we present numerical simulations. We first consider the 2D guiding-center model in a D-shaped domain. The steady state solution and perturbed solution are studied. Then we perform the ion turbulence instability simulation with the 4D Drift-Kinetic model in a cylinder domain.

4.1. Test 1 : Steady state solution for the guiding center model in a D-shaped domain. We consider now the 2D guiding-center model in a D-shaped domain Ω presented in Section IV of [START_REF] Miller | Noncircular, finite aspect ratio, local equilibrium model[END_REF] and depicted in Figure 2(a). The mapping X from curvilinear coordinates ξ = (ξ 1 , ξ 2 ) to physical coordinates x = (x 1 , x 2 ) is given by x = 1.7 + [0.074(2ξ 1 -1) + 0.536] cos[2πξ 2 + arcsin(0.416) sin(2πξ 2 )],

y = 1.66[0.074(2ξ 1 -1) + 0.536] sin(2πξ 2 ), for -231 74 ≤ ξ 1 ≤ 1, 0 ≤ ξ 2 ≤ 1.
We now search a steady state solution for the guiding-center model in the D-shaped domain Ω. We first notice that ρ(φ), a function of φ, is a solution of the guiding-center equation (2.9). Then it remains to solve a nonlinear elliptic equation In the previous proposition, the existence of a positive solution is obtained by Schauder's fixed point theorem [START_REF] Evans | Partial Differential Equations[END_REF], while the uniqueness is straight. We now solve (4.1) in Ω. Figure 2(b) illustrates that the boundary ∂Ω is embedded in Cartesian mesh. Thus the numerical scheme presented in section 3.3 should be applied, which gives us a nonlinear system for φ. Then by applying a Newton method to this nonlinear system, we obtain a steady state solution φ 0 of (4.1) shown in Figure 3(a). Then by injecting φ 0 into the guiding-center equation (2.9), we get the steady state density ρ0 (see Figure 3(b)). Moreover, we plot the velocity field of steady state solution in Figure 4. By comparing the streamline in Figure 4 and the constant line of coordinates ξ = (ξ 1 , ξ 2 ), it is interesting to notice that these lines don't coincide, especially in the core of the D-shaped domain. The velocity is fast near the edge of the D-shaped domain but much slow in the core. Now we verify that (φ 0 , ρ0 ) is the steady state solution by long time simulation. That is to take the pair (φ 0 , ρ0 ) as an initial solution for the guiding-center model (2.9), then to compare the difference between (φ(t), ρ(t)) and (φ 0 , ρ0 ). We will measure these differences by a relative error as

(4.1) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ -∇ ⋅ ρ 0 B ∇ φ = ρ(φ) -ρ 0 in Ω, φ = 0 on ∂Ω.
E(u(t)) = u(t) -u 0 1 u 0 1 , u = φ, ρ.
The finite difference method with HWENO reconstruction is used for solving guiding-center equation (2.9). The time step is taken to be ∆t = 0.001. Figure 5 presents the relative errors of the potential φ and the density ρ. We observe that the solution remains steady for long time simulation with a relative error of magnitude of 10 -4 .

4.2.

Test 2 : Numerical simulation of the guiding center model in a D shape. Now we still consider the previous initial data (φ 0 , ρ0 ) which is a stationary solution of the guiding-center model, but perturb it of magnitude of ε.

In Figure 4, we have noticed that the streamline is different with respect to the constant line of coordinates (ξ 1 , ξ 2 ). On the other hand, we denote Then by the definition of streamline, we have

U = ⎛ ⎝ U x U y ⎞ ⎠ = ⎛ ⎜ ⎝ -∂φ ∂y ∂φ ∂x ⎞ ⎟ ⎠ . 1 1
U y -U x dy dx = 0, which implies d dx φ(x, y(x)) = ∂φ ∂x + ∂φ ∂x dy dx = U y -U x dy dx = 0.
Thus, φ(x, y) = const represents a streamline, i.e. the steady state function ρ0 revolves along the isoline of potential φ 0 . In this test case, we perturb the function ρ0 along the streamline, that is

ρ = ρ0 (1 + ε cos(2πkξ 2 ) exp(-2 φ 0 -φ p 2 ε 4 ),
with φ p = -0.1, k = 5 and ε = 0.1. Figure 6 illustrates the evolution of density governed by the guiding-center model. We present the difference between the perturbed density and the steady state density, i.e. δρ(t) = ρ(t) -ρ0 . We observe that the difference of density δρ revolves, and small filaments appear at time t = 200. Until the time t = 300, we can clearly identify the filaments.

4.3. Test 3 : Drift-kinetic model for ion turbulence simulation. In this subsection, we reproduce the ion turbulence simulation [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]. This simulation has been realized by different methods but in cylindrical coordinates [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Samtaney | Numerical aspects of drift kinetic turbulence: ill-posedness, regularization and a priori estimates of sub-grid-scale terms[END_REF]. We will perform the simulation in Cartesian mesh with the numerical methods presented in section 3.

The discretization of the Drift-Kinetic model can be developed very similarly as the one for the guiding-center model. Here, we present some principle discretization steps.

The Vlasov equation of system (2.5) can be split into three equations :

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂f ∂t + U ⋅ ∇ x f = 0, ∂f ∂t + v ∂ z f = 0, ∂f ∂t + E ∂ v f = 0.
Thus when applying the Semi-Lagrangian method, we use the Strang splitting method [START_REF] Strang | On the construction and comparison of difference schemes[END_REF] for time discretization, and a second order predictor-corrector method for searching the characteristic foot. The Semi-Lagrangian method is suitable for linear phase simulation, because it doesn't have CFL constraint and it is very accurate. However, it doesn't preserve well the conservation properties of physical models [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially-non-oscillatory reconstruction for Vlasov equations[END_REF] in nonlinear phase. Therefore, we should use the conservative finite difference method presented in section 3.2, where the 4th order Runge-Kutta method is used for time discretization. The criterion to pass from the Semi-Lagrangian to the finite difference methods is as follows

(4.2) Ω R [f (t n ) -f (t n-1 )] dvdx > h 3 ,
where h is the smallest space step. The quasi-neutrality equation of the system (2.5) is a three-dimensional elliptic problem. Obviously, a direct resolution of this equation is very costly. However, we notice that the diffusion term

ρ 0 (x ) B
only depends on x , and the second term of the quasi-neutrality equation is vanishing by taking average in z-direction. Thus averaging the the quasi-neutrality equation in z-direction, we get a 2D average equation

(4.3) -∇ ⋅ ρ 0 (x ) B ∇ φ = ρ -ρ 0 in D.
Taking difference between the quasi-neutrality equation and the average equation, it yields a fluctuation equation :

(4.4) -∇ ⋅ ρ 0 (x ) B ∇ φ ′ + ρ 0 (x ) T e (x ) φ ′ = ρ ′ -ρ in Ω = D × [0, L z ],
with φ ′ = φ -φ. Note that the terms

ρ 0 (x ) B and ρ 0 (x )
Te(x ) are independent on z. Thus the fluctuation equation (4.4) can be solve slice by slice in z-direction. Moreover, the Dirichlet boundary conditions can also be decomposed for the average equation (4.3) and the fluctuation equation (4.4) as follows

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ φ = 0, ∀x ∈ ∂D, φ ′ = 0, ∀x ∈ ∂D × [0, L z ].
Therefore, the resolution for Poisson equation in Section 3.3 can be applied for (4.3) and (4.4). Now we perform the ion turbulence simulation. The plasma is initialized by exciting a single ion temperature gradient (ITG) model (m, n) (where m is a poloidal mode and n is a toroidal mode). The distribution function is thus considered at the initial time as the sum of an equilibrium and a perturbed part: f = f eq + δf . The equilibrium part f eq is chosen as a local Maxwellian

f eq (r, v ) = n 0 (r) (2πT i (r)) 1 2 exp ⎛ ⎝ - v 2 2T i (r) ⎞ ⎠ ,
while the perturbation δf is determined as

δf = f eq ε exp - (r -r p ) 2 δr cos 2πn L z + mθ ,
where the profiles n 0 (r), T i (r) and T e (r) satisfy ∂ r P (r) P (r) = -κ P cosh -2 r -r p δr P , for P = n 0 , T i and T e , together with the normalization Let us first compare the different discretization methods. The Semi-Lagrangian methods with cubic Hermite reconstruction and the HWENO reconstruction are used to solve the 4D Drift-Kinetic model, with small time step such that the CFL number is small than 1. Then we compare the numerical results with the one obtained by the mixed Semi-Lagrangian/finite difference method, where large time step (CFL> 1) is used for the Semi-Lagrangian method in linear phase and small time step (CFL< 1) is used in nonlinear phase. We emphasize that the Semi-Lagrangian method switches to the finite difference method automatically by the criterion (4.2).

In Figure 7, we summarize relative errors of the conservation laws for the Drift-Kinetic model for the different methods. We notice that these three methods have almost the same results when t < 3000, while these results differ significantly when t ≥ 3000. We thus denote the linear phase for t < 3000 and the nonlinear phase for t ≥ 3000.

The Semi-Lagrangian methods can not conserve well the mass in the nonlinear phase, while the finite difference method conserve exactly the mass (see Figure 7(a)). Then from Figures 7(b), 7(c), we observe that the Semi-Lagrangian method with Hermite reconstruction loses completely the conservation properties for L 2 norm and entropy for long time simulation, since it involves too much spurious oscillation. At contrast, the Semi-Lagrangian method and the finite difference method with HWENO reconstruction work much better. Finally, we see the mixed Semi-Lagrangian/finite difference method has better energy conservation property than the Semi-Lagrangian methods in the nonlinear phase (see Figure 7(d)). Therefore, the mixed Semi-Lagrangian/finite difference method is better than the Semi-Lagrangian methods for long time ion turbulence simulation. Moreover, since the Semi-Lagrangian method is used in linear phase, thus our mixed method is more efficient than the pure finite difference method.

We next investigate the L 2 norm and entropy convergence properties of the mixed Semi-Lagrangian/finite difference method by refining mesh size in different directions. The results calculated with mesh size n x × n y × n z × n v = 64 × 64 × 32 × 65 is used as a reference solution. In Figure 8, we observe that a significant improvement is obtained by refining in x, y directions. The refinement in z direction doesn't improve the results. The ones obtained by refining in v direction is slightly better than the reference solution. Figure 9 presents the refinement results for energy conservations. We see again the refinement in x, y directions improves best the energy conservations. At last, we present the evolution of the distribution function during ion turbulence simulation. We first notice in Figure 10 that the instability develops exponentially in the linear phase, where the growth of instability is measured by the quantity φ(t, x, y, z) 2 dx dy dz, where (x, y) ∈ {(x, y) ∶ x 2 + y 2 = r 2 p }. Then the instability reaches a saturation point, which corresponds to the starting point of the nonlinear phase.

In Figure 11, we show the evolution of distribution function f at v = 0. We see that till the time t ≤ 2000 the instability can not be identified very clearly, that is why the Semi-Lagrangian method can be applied in the linear phase. At time t = 3000, we reaches the saturation point, and five vortices are developed. These vortices rotate and create small filaments. At time t = 4000, more small structures appear in the distribution function. At this moment, the Semi-Lagrangian method can not conserve well the invariant quantities, while the conservative finite difference method performs much better. Finally, the instability continues to develop small structures of the distribution function till the mesh size, and we attain to a relatively steady state.

Conclusion and perspective

In this paper, we have presented an efficient algorithm for long term plasma simulations. We first derive the 4D drift-kinetic and the 2D guiding-center models, and present their conservative properties. The Hermite WENO reconstructions are applied for solving the Vlasov equations, which was proved to be robust [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially-non-oscillatory reconstruction for Vlasov equations[END_REF] in computational performance. Moreover, to adapt the arbitrary computational domain, we discretize the models on Cartesian meshes, and the special numerical methods for the boundary conditions, as the inverse Lax-Wendroff method for the Vlasov equation [START_REF] Filbet | An inverse Lax-Wendroff method for boundary conditions of Boltzmann equations[END_REF] and the extrapolation method for the poisson equation [START_REF] Filbet | Numerical Simulations of Kinetic Models for Chemotaxis[END_REF], are proposed. Next, we solve the guiding-center model on a D-shape domain. A steady state solution (φ 0 , ρ0 ) is found numerically. Then we perturb the steady state density ρ0 along the streamline, and use this perturbed density ρ as the initial condition for the guiding-center model. We observe that the difference of density δρ = ρ -ρ0 revolves, and the filaments appear for long term simulation.

Finally, we simulate the 4D drift-kinetic model with the mixed methods, i.e. the semi-Lagrangian method in linear phase and finite difference method during the nonlinear phase. Numerical results show that the mixed method is efficient and accurate in linear phase and it is much stable during the nonlinear phase. Moreover, it preserves well the conservative properties. We thus conclude that our mixed method is efficient for realistic and high dimensional plasma turbulence simulations.

3. 3 .

 3 Discretization of the Poisson equation (2.9)-(2.11). We use a classical five points finite difference approximation to discretize the Poisson equation (2.9)-(2.11). So it remains to treat the Dirichlet boundary conditions on ∂D.

Figure 2 .

 2 Figure 2. D-shaped domain. (a) Constant lines in coordinates ξ = (ξ 1 , ξ 2 ); (b) D-shaped domain embedded in Cartesian mesh.

For a suitable functionProposition 4 . 1 .

 41 ρ, we have a unique solution of equation (4.1) : Consider the equation (4.1) with ρ 0 = 1, B = 1, ρ(φ) = e -φ -1. The function ρ is bounded on C = {φ ∈ H 1 0 (Ω) ∶ φ ≥ -ln 2} ⊂ H 1 0 (Ω). Then (4.1) has a unique weak solution in C.

Figure 3 .

 3 Figure 3. A steady state solution of the guiding-center model (2.9) in Dshaped domain. Mesh size is n x × n y = 240 × 440.

Figure 4 .Figure 5 .

 45 Figure 4. Streamline of the velocity field corresponding to steady state solution of guiding-center model (2.9) in D-shaped domain.

Figure 6 .

 6 Figure 6. Instability simulation for guiding-center model in D-shaped domain. The difference between the perturbed density and the steady state density is presented, i.e. δρ(t) = ρ(t) -ρ0 .

rmax r min n 0

 0 (r) dr = r max -r min , T i (r p ) = T e (r p ) = 1. This gives the formulas P (r) = C P exp -κ P δr P tanh r -r p δr P , where C T i = C Te = 1 and C n 0 = , we choose the following parameters r min = 0, r max = 14.5, κ n 0 = 0.055, κ T i = κ Te = 0.27586, δr T i = δr Te = δrn 0 2 = 1.45, ε = 10 -6 , n = 1, m = 5, L = 1506.759067, v max = 8, r p = rmax+r min 2 , δr = 4δrn 0 δr T i .

Figure 7 .

 7 Figure 7. Comparison of different reconstruction methods for the conservation laws for the Drift-Kinetic model. The mesh size is n x × n y × n z × n v = 128 × 128 × 32 × 65.

Figure 8 .

 8 Figure 8. Refinement of mesh size in different directions for L 2 norm and entropy conservations of the Drift-Kinetic model. Mixed Semi-Lagrangian/finite difference method is used.

Figure 9 .Figure 10 .

 910 Figure 9. Refinement of mesh size in different directions for the energy conservations of the Drift-Kinetic model. Mixed Semi-Lagrangian/finite difference method is used.

Figure 11 .

 11 Figure 11. Evolution of ion turbulence. The distribution function is shown for the velocity v = 0. The mesh size is n x = n y = 128, n z = 32, n v = 65. Mixed Semi-Lagrangian/finite difference method is used.

Acknowledgment Both authors are partially supported by the European Research Council ERC Starting Grant 2009, project 239983-NuSiKiMo and the Inria project Kaliffe. Chang YANG is also supported by Fundamental Research Funds for the Central Universities and Program for Innovation Research of Science in Harbin Institute of Technology (PIRS OF HIT T201402)