Benjamin Vincent
email: benjamin.vincent@cnrs-orleans.fr

Sedina Tsikata

Stéphane Mazouffre

Incoherent Thomson scattering measurements of electron properties in a conventional and magnetically-shielded Hall thruster

Reliable measurements of electron properties are key for the validation of simulations of Hall thrusters and similar plasma sources. In this paper, a recently-developed incoherent Thomson scattering diagnostic is successfully applied to electron property measurement in a low-power Hall thruster discharge. Investigations of electron properties along radial and azimuthal directions in both conventional and magnetically-shielded architectures are performed. The high diagnostic sensitivity gives access to electron property measurements in a plasma environment with densities as low as 10 16 m -3 . Electron temperatures reaching several tens of eV and drift velocities on the order of 10 6 m s -1 are measured in the exit plane region. The impact of the magnetic field (intensity and direction) and discharge voltage are investigated.

Introduction

A Hall thruster is an electric propulsion device first conceived and developed in the mid-1960s in the Soviet Union. This thruster was operated for the first time in space aboard the USSR Meteor satellite in 1972 [1][START_REF] Kim | 30th International Electric Propulsion Conference[END_REF][START_REF] Mazouffre | [END_REF][4] and is today widely used for station-keeping on telecommunications satellites [START_REF] Goebel | Fundamentals of Electric Propulasion[END_REF][START_REF] Charles | [END_REF][START_REF] Hoskins | The 33rd International Electric Propulsion Conference[END_REF]. Interest in such technologies is driven by the significant fuel mass economies achievable in comparison to chemical propulsion-based systems. In spite of the long flight heritage of these devices, an understanding of their underlying physics remains poor: questions such as the nature of instabilities occurring over different temporal and spatial scales [START_REF] Choueiri | [END_REF], interactions between the thruster plasma and external cathode plasma [9], the origins of non-Maxwellian electron velocity distribution functions [10,11], the nature of plasmawall interaction and the consequences of secondary electron emission [12][13][14][15], all remain to be fully understood. Early simulation efforts have shed light on aspects such as Hall thruster instabilities and particle dynamics, illustrated for example in Refs. [START_REF] Hirakawa | Numerical simulation of plasma particle behavior in a hall thruster 32[END_REF][START_REF] Adam | [END_REF][18][START_REF] Hofer R R | Collection of Technical Papers[END_REF] and newer modeling efforts have been undertaken in several research teams in recent years, including in the teams of Refs. 4, 20-24, in order to advance physical understanding and ultimately, make possible fullypredictive discharge modeling. However, all codes must be fully validated through comparison with accurate experimental measurements for electron properties.

Electron property measurement in Hall thrusters has so far been extensively performed with Langmuir probes, including in Refs. [25][START_REF] Shastry | Near-Wall Plasma Properties and EEDF Measurements of a 6-kW Hall Thruster 45[END_REF][START_REF] Dannenmayer | [END_REF][28], to name only a few. Optical emission spectroscopy has also been used as an alternative method to estimate electron temperature in thrusters in works including Refs. 29-32. Such techniques are adequate for far plume measurements, however, their application in the acceleration region of the Hall thruster plasma has limitations [33]. The expected high electron densities and temperatures and the presence of strong magnetic fields in this region complicate signal measurement and interpretation. More specifically, for Langmuir probes, the measurement duration must be minimized with fast translation stages in order to avoid sputtering, and the helical trajectories of magnetized electrons require more complex probe theories to interpret the signals acquired. For optical emission spectroscopy, the use of line intensities in collisional-radiative models to determine electron temperatures is associated with inaccuracies for both low and high electron temperatures [30,[START_REF] Sommerville | An Optical Diagnostic for Xenon Hall Thrusters Including Metastable Contributions 42[END_REF].

Thomson scattering applied to such measurements offers the possibility of unambiguous and perturbationfree measurement of electron properties. The challenge associated with this technique is achieving a high enough diagnostic sensitivity for measurements in a low-density plasma, such as that of the Hall thruster. In the coherent regime, Thomson scattering has already proven to be an efficient technique to measure electron density fluctuations in such plasmas, seen in the implementation of the PRAXIS diagnostic [START_REF] Tsikata | Small-scale electron density fluctuations in the Hall thruster, investigated by collective light scattering[END_REF]. In the incoherent regime, Thomson scattering allows measurement of thermal properties of electrons. This technique has already been successfully used in several low-pressure plasma sources to estimate electron density and temperature, including in Refs. [START_REF] Irisawa | [END_REF][36][37][38][39][40][41] or even the electron drift velocity in Ref. 42.

Recent studies by Washeleski [START_REF] Washeleski | Laser Thomson scattering measurements of electron temperature and density in a Hall-effect plasma[END_REF][START_REF] Washeleski | [END_REF] have shown that this technique can be used in Hall thrusters, but the low signal levels and large stray light intensities make its application challenging. In our previous work, a sufficiently sensitive diagnostic named THETIS (THomson scattering Experiments for low-Temperature Ion Sources) was developed and successfully tested on a cathode discharge [45] and a planar magnetron [46]. In this paper, the results from the implementation of this diagnostic on a low-power Hall thruster are presented.

In Section 2, a description of the main features of the diagnostic and the data analysis procedures used are presented. Results of experimental investigations on two thruster architectures under different discharge conditions are discussed in Section 3. A summary of the main results is given in Section 4.

Experimental schemes

Overview of the experimental set-up 2.1.1. Test facility

Experiments were performed in the NExET (New Experiments in Electric Thrusters) vacuum chamber at ICARE. This chamber is a stainless steel cylinder 0.7 m in diameter and 1.7 m in length, equipped with a primary dry pump, turbomolecular pump, and cryogenic pump. The pressure was measured with a Leybold Penningvac PTR 90 gauge positioned at the top of the chamber, approximately 1 m downstream of the gas injection position .The pumping assembly achieves a base pressure of 10 -4 Pa and operating pressure below 2 × 10 -2 Pa with a xenon mass flow rate up to 2.5 mg s -1 . Background pressure effects are not negligible [47] for these operating conditions. Ionization of the ambient gas is expected to contribute to higher discharge currents than at lower chamber pressures, for the same applied voltage. Slightly higher electron densities (due to ambient ionization) and lower electron temperatures (due to increased electron-neutral collisions) are also expected under our operating conditions. Such effects have been observed for previous studies investigating facility effects, for example, in Ref. [START_REF] Walker | Effects of Facility Backpressure on the Performance and Plume of a Hall Thruster[END_REF].

A c c e p t e d M a n u s c r i p t

The discharge voltage is regulated with an SM1500 Delta Elektronika DC power supply. Discharge current is measured through a 0.5-1.0 W Stangenes current transformer probe and monitored with a LeCroy HDO6104 digital oscilloscope.

Plasma sources

Incoherent Thomson Scattering (ITS) measurements were performed on two 200 W Hall thrusters of the ISCT200 (ICARE Small Customizable Thruster 200W ) type: the ISCT200-GEO and ISCT200-MS. Both thrusters were equipped with BNSiO 2 channels and have the same channel width and mean diameter, a geometry referred as "2S 0 " in previous publications [START_REF] Lejeune | Impact of the channel width on Hall thruster discharge properties and performances IEPC[END_REF]. The so-called 2S 0 geometry corresponds to a configuration with a channel width to diameter ratio twice as large as the usual value. In this paper, thrusters were operated with xenon injected through a porous ceramic gas distributor placed behind a stainless steel ring anode, as seen in Fig. 2.

The ISCT200-GEO has been extensively studied at the ICARE laboratory [START_REF] Mazouffre | [END_REF]. It has a standard magnetic field topology, achieved using an assembly of permanent magnets, with the maximum field intensity at the center of the channel and inside the exit plane. The magnetic field intensity is changed by varying the number of magnets of the internal and external magnetic circuits. Three peak magnetic fields were used in this study: B low (15 mT), B mid (25 mT) and B high (31.5 mT). The ISCT200-GEO is shown in Figure 1a. The typical magnetic profile of this thruster is illustrated in the upper half of Figure 2.

The ISCT200-MS has been studied in recent years [51,52] and has a magnetic circuit designed for "magnetic shielding", a concept originally developed by researchers at NASA-JPL [START_REF] Mikellides | Magnetic Shielding of the Acceleration Channel Walls in a Long[END_REF] and recently extended to miniature Hall thrusters (in the 300-500 W range) by Conversano and colleagues [START_REF] Conversano | [END_REF]. This configuration is designed to reduce channel erosion [START_REF] De Grys | Demonstration of 10,400 Hours of Operation on 4.5 kW Qualification Model[END_REF] and the maximum magnetic field position is slightly shifted downstream of the exit plane.

In this work, measurements with a peak magnetic field of 15 mT (B low ) are discussed. The ISCT200-MS is shown in Figure 1b and the typical magnetic field line shape of this shielded thruster is illustrated in the lower half of Figure 2.

Both thrusters were operated with a 5 A MIREAtype hollow cathode equipped with a LaB 6 emissive element heated with 130 W of power and run with a xenon flow rate of 0.2 mg s -1 .

The laser beam used to probe the plasma was kept fixed. The thrusters were translated along the x-axis (Figure 3) for axial exploration of electron properties. Displacements along the y-axis and z-axis were used to change the scattering configuration.

Description of the ITS diagnostic

A detailed description of the recently-developed ITS diagnostic named THETIS and the procedure for data analysis can be found in a published article [45]. Only a brief description is made here.

Optical elements

A schematic of the THETIS diagnostic transmission and detection branches is shown in Figure 3.

The transmission branch produces the scattering signal. This is drawn with a dark green line in Figure 3. It starts with a frequency-doubled, 10 Hz pulsed Q-A c c e p t e d M a n u s c r i p t switch Nd:YAG laser delivering pulses of 430 mJ (pulse width 6 ns). Mirrors and lenses conduct and focus the laser beam into the observation volume of the detection branch. The laser beam is finally sent to a largeaperture beam dump to reduce back reflections.

The detection branch is used to analyze the Thomson scattering signal and is indicated with a light green line in Figure 3. It starts at the observation volume position, where the light from the transmission branch has been focused. From this volume, the Thomson scattering radiation and other stray radiation is emitted. This signal is focused with a collection lens onto the entrance of an optical fiber bundle, composed of a 5 × 3 array of 300 µm-diameter fibers. In the configuration of these experiments, the collection lens referred as L1 is a 100 mm-diameter plano-convex lens with a 200 mm focal length positioned at the top of the vacuum chamber 0.564 m from the observation volume (magnification factor of 0.55). The fiber bundle transmits the collected light to the optical table, where light from a 1 × 15 array of the fiber stack is collimated with a 50 mm-diameter plano-convex lens of 125 mm focal length. The stray light from reflections and Rayleigh scattering are then attenuated with a 25 mm × 25 mm VBG-NF (Volume Bragg Grating Notch Filter). The remaining light is focused onto the 1 mm-wide entrance slit of a spectrometer using a 300 mm focal length lens to fit the f/9.7 numerical aperture of the spectrometer. Light is then typically dispersed by the spectrometer with a 600 lines mm -1 grating suited for measurement of the high electron temperatures (tens of eV) expected at the Hall thruster exit plane. Finally, the dispersed light is detected with an ICCD camera with a 10 ns gate, synchronized with the Q-switch signal from the laser.

Scattering configuration

With this detection branch configuration, the observation volume is 0.3 mm in diameter and 2.7 mm in length. Plasma properties measured with the diagnostic correspond to an average over this volume.

Electron properties along an observation wave vector k are probed with Thomson scattering diagnostics. As represented in Figure 4, the observation wave vector direction in this work is defined by the vector equation k = k i -k s . k i corresponds to the wave vector of the incident laser beam and k s corresponds to the wave vector of the light entering the detection branch.

In these experiments, the thruster vertical position was chosen so that the laser beam passes at a height r mean × cos(45°) above the midplane (where r mean is the mean channel diameter). With this configuration, a lateral translation of the thruster allows the alignment of the observation wave vector along either the radial or azimuthal direction, if the observation volume is positioned in front of the channel mean diameter.

By convention, a positive spectral shift of the Thomson spectrum corresponds to a positive electron drift velocity. With the data analysis convention used, a positive drift velocity measured along the azimuthal direction corresponds to electrons drifting in the anticlockwise direction (i.e. along k). Along the radial direction, a positive drift will correspond to an electron moving towards the internal pole of the thruster.

Data analysis

For calibration of the transmission factor of the detection branch, the same procedure as that discussed in Ref. 45 is followed. The vacuum chamber was filled with nitrogen to a pressure of 1000 Pa and the signal from Raman scattering accumulated over 6000 laser pulses. Due to the high electron temperatures expected with the Hall thruster, the Raman calibration, as well as most of Thomson spectrum acquisitions, were performed with the 600 lines mm -1 grating and a central wavelength of 532 nm. The spectrometer slit width was set to 1 mm. With this grating and slit width the Boltzmann distribution of the vibrational states, and not the individual Raman lines, is resolved.

In our recent work, the Thomson scattering signal analyzed is obtained from the subtraction of two consecutive spectra types accumulated over 6000 laser pulses. The first record ("Thomson"), is obtained with the plasma and laser on. During this record, the Thomson scattering, plasma emission, electronic noise and the remaining stray light from Rayleigh and reflections are measured. The second record, ("background"), is obtained with the plasma on but the laser off. During this record, only the plasma emission and electronic noise signals are measured. For both plasma sources studied to date (the hollow cathode and planar magnetron) the subtraction of these two records were enough to isolate Thomson scattering from plasma emission lines.

For analysis of the Thomson spectra in the present work, the data analysis procedure has been modified with respect to the procedure described in Ref. 45 to be more robust to potential distortion from plasma emission lines. Due to the high electron temperatures, the spectral range covered by the Thomson scattering signal is large. Inside the 28.6 nm spectral range observed, several intense emission lines from the Hall thruster plasma are present. Even during the 10 ns observation gate of the detector, some emission lines remain at amplitudes at least one order of magnitude above the maximum Thomson scattering signal. The subtraction of the plasma emission record can lead to dips in the final spectrum (for example, near 530 nm and 542 nm in the Thomson spectrum of Figure 5), or leave stray peaks (for example near 526 nm in the A slow drift in the plasma conditions due to thermalization of the thruster could lead to a drift of plasma emission intensity. This slight difference occurring between the "Thomson" and "background" records would appear significant in comparison to the weak Thomson scattering signal. Another possible reason is that during "Thomson" records, the emission lines are measured in the presence of the intense electric field of the laser. This electric field may induce significant Stark splitting and shifting [START_REF] Ernst | [END_REF]57]. During the "background" records such splitting is absent, the recorded emission lines would not be identical to those obtained during "Thomson" records.

V + 45° k s k k i k k s k i
To limit distortion by the emission lines, the Thomson spectra are analyzed in the following way before fitting. Negative values of intensity (counts) below three times the standard deviation of the electronic noise level from the iCCD (σ elec noise ) can only be attributed to distortion from plasma emission lines following subtraction. Intensities of the acquired Thomson spectra below this critical value were truncated at this value (3 × σ elec noise = -2.5 counts with the experimental parameters used during the campaign). To mitigate the positive distortion of the Thomson spectrum (i.e., cases where the residual emission lines add up positively to the Thomson signal), a Savitzky-Golay filter was applied to the Thomson spectrum. This low-pass filter reduces the fast rises of intensity due to plasma emission lines, while leaving unchanged the slowly-varying Thomson spectral shape. This filter was applied independently on the left and right sides of the Thomson spectrum, separated by the rejected central spectral range around 532 nm. Data obtained from this preliminary data processing is used as input to the algorithm for Thomson spectral analysis presented in Ref. 45. A Gaussian function convoluted with the instrument function is used as a fitting function for the Thomson spectrum. From the Gaussian standard deviation, an electron temperature is determined, while the shift of the Gaussian mean from the laser wavelength is used to determine the bulk electron drift velocity. Finally, the integrated intensity of the Thomson spectrum gives an estimation of the electron density averaged inside the observation volume.

Similar to the procedure presented Ref. 45, uncertainties were estimated for each electron parameter. Uncertainties were systematically calculated from the quadratic mean of the propagation of the main sources of uncertainty. The principal sources of this uncertainty were the uncertainty in the value of the differential cross-section for Raman scattering of N 2 (≈8 %), the laser energy fluctuations (<1 %) and the quality of the fit to the data, determined from the co-variance matrix obtained during fitting, and dependent on the signal-to-noise ratio of the fitted data.

The spectrum example provided in Figure 6 shows a case obtained for lower electron temperatures, for which both the 600 lines mm -1 grating and a less-dispersive 2400 lines mm -1 grating were used to constitute the spectra. Origins of the weakly non-Gaussian spectrum shape will be discussed in a forthcoming paper.

Results and discussion

Preliminary remarks

To determine whether the electron properties are isotropic, measurements along both the radial and azimuthal directions were made by switching between the two scattering configurations presented in Section 2.1.3. In the subsequent plots, values of electron properties obtained along the azimuthal direction are plotted with blue triangles and radial observations plotted with yellow circles. While the observation of non-isotropic electron temperature or drift velocity has a physical significance, this is not the case for electron density, which is expected to be the same in both directions. Differences in the electron density measured at the azimuthal and radial positions may be due to azimuthal plasma non-uniformity (possibly influenced by the cathode position [58]), spatial inhomogeneities of the magnetic field or gas feeding [START_REF] Kawashima | Effects of Azimuthal Nonuniformity on the Hall Thruster Discharge IEPC[END_REF]). Electron beam fluorescence [START_REF] Jarrige | Gas in a Hall Effect Thruster using Electron Beam[END_REF] might be a way to verify the azimuthal uniformity of the neutral gas injection in the future. A slight displacement (below 0.3 mm) of the measurement position as the thruster is translated laterally to A c c e p t e d M a n u s c r i p t switch between the two scattering configurations may also account for small differences in measured density between the azimuthal and radial positions.

As mentioned previously, for reliable Thomson spectrum measurement, an average over 6000 laser pulses was necessary to increase the signal-to-noise ratio. For the presented study, the Thomson scattering spectrum is measured at random phases in the discharge current oscillation and is therefore expected to be representative of the mean electron properties. Phase-resolved measurement of electron properties will be performed in forthcoming studies, and has indeed already been performed in another of our studies on pulsed plasmas [46].

Axial exploration of electron properties was performed within the axial range 2-40 mm in the azimuthal direction and 3-40 mm in the radial direction with respect to the exit plane. The closest approach to the exit plane for both directions is that for which stray light reflections can be sufficiently attenuated for the Thomson signal to be unobscured. The difference in the closest approach for the two directions is due to the different levels of light reflected off the central pole which enter the detection branch. While with the radial scattering configuration the inner pole behind the observation volume is directly illuminated, it is not the case with the azimuthal scattering configuration. The furthest axial distance is that for which the plasma density remains sufficient for a signal to be detected. As shown with the Thomson spectrum example in Figure 5, a Gaussian curve can be fitted to the data with confidence despite the low electron density.

Influence of the magnetic field

Magnetic field intensity

Axial profiles of electron properties with the ISCT200-GEO for various magnetic field intensities were determined.

For each magnetic field, the discharge power was kept constant at 285 W, with the discharge voltage and current of respectively 250 V and 1.14 A. The discharge current was kept constant using an adjustment of the xenon mass flow rate to 1.18 mg s -1 , 1.60 mg s -1 and 1.38 mg s -1 for, respectively, the B low , B mid and B high magnetic field intensities. Results from these explorations are shown in Figure 7.

As shown in the first row of Figure 7, the electron density has only a weak dependence on the magnetic field intensity. The different values of electron density remains between 4.5 × 10 16 m -3 and 0.8 × 10 16 m -3 from 2 mm to 40 mm from the exit plane. The electron density is lowest with the B low configuration (possibly attributable to less efficient confinement) and the highest electron density is observed with the nominal B mid magnetic field.

In Figure 7, the electron temperatures of the B mid configuration decrease monotonically from approximately 40 eV near the exit plane to 10 eV at 40 mm downstream, for both azimuthal and radial directions. These values are on the same order of magnitude as recent PIC simulations predictions [START_REF] Boeuf | [END_REF].

In the B low configuration, the electron temperatures are higher than for the B mid case, and azimuthal temperatures exceed the radial temperatures. There are some possible reasons for such anistropy in temperature. A weaker magnetic trap, provided a similar electric field amplitude is achieved for the B low and B mid configurations, would lead to higher E × B drift velocities. We note that the E × B drift is the driver of azimuthal instabilities like the electron cyclotron drift instability [62], confirmed in thruster simulations to contribute to electron heating in both azimuthal and radial directions (with a stronger contribution in the azimuthal direction) [14,63]. Radial electron temperatures may also be lower due to radial losses of energetic electrons to the walls and the injection of cooler wallemitted electrons into the plasma. This latter effect might, however, be more significant inside the thruster channel. In order to fully understand the causes of the anisotropy, complementary measurements of (i) the electric field axial profiles using laser-induced fluorescence, and (ii) instabilities implicated in heating using coherent Thomson scattering would be necessary. Investigations of the influence of geometrical effects on anisotropy (varying channel widths at fixed mean channel diameters) could also be instructive.

The reproducibility of the electron properties presented in this paper has been confirmed by a few successive measurements performed under identical plasma conditions, which produce identical spectra in amplitude and width.

Variations from one measurement to another are due generally to variations in the level of the emission lines.

The time-varying discharge current and the associated power spectrum density for the different magnetic fields are shown in Figure 8 and Figure 9. The integral of the power spectrum density over frequency between 10 3 Hz and 2 × 10 7 Hz gives 9.8 × 10 -3 A 2 , 3.6 × 10 -3 A 2 and 81 × 10 -3 A 2 for respectively the B low , B mid and B high configurations. For these configurations, the electron density decreases as the power contained in discharge current oscillations increases. The larger discharge current oscillations observed in the B low and B high configurations may be linked to axial plasma oscillations (such as transit time oscillations [START_REF] Boeuf | [END_REF]) which may contribute to a larger anisotropy in the electron temperatures and may be associated with poorer electron confinement. This study on the correlation between the power density spectrum of the discharge current and the presence of anisotropic electron temperature values should be repeated in order to verify the statistical significance of this observation. An order of magnitude estimate of the expected electron density can be determined using the expression for the xenon mass flow rate:

ṁn = 1 η prop ṁi = 1 η prop n i m i Av i (1)
where ṁn and ṁi are the mass flow rate of xenon atoms and ions, respectively, through a section A. η prop is the propellant utilization efficiency. m i and v i correspond to, respectively, the mass and velocity of xenon ions while n i is the density. Under the hypothesis of quasineutrality this equation can be used to estimate n e :

n e = η prop ṁn m i Av i (2)
At the position 10 mm downstream of the exit plane, we may assume that most of the ions have been accelerated. Using previously-estimated values of the parameters of Equation 2 with the same thruster at similar discharge conditions (η prop ≈80 % [START_REF] Grimaud | Magnetic shielding topology applied to low power Hall thrusters[END_REF], v i ≈16 km s -1 [START_REF] Bourgeois | Influence de la topologie magnétique[END_REF] and a channel surface A expanded by a 45°plume divergence angle) we obtain an estimated electron density of 10 17 m -3 at 10 mm from the exit plane. The slightly lower values measured at this position (a few 10 16 m -3 in Figure 7) may be attributed to lower ionization efficiency or higher thruster divergence than that used to estimate the density.

Magnetic field configuration

The difference in electron properties between the standard (ISCT200-GEO) and magnetically-shielded configurations (ISCT200-MS) was investigated with similar discharge conditions. To allow a direct comparison, the thrusters were operated at identical discharge power levels. The ISCT200-GEO was operated at 250 V, 1.18 mg s -1 in the B low configuration, with the magnetic field directed towards the outer pole. These conditions gave a discharge power of 285 W. The shielded ISCT200-MS was operated at 250 V, 1.05 mg s -1 for a discharge power of 285 W. The corresponding magnetic field was directed towards the inner pole.

The comparisons between the two architectures are shown in Figure 10. As seen in the first row of Figure 10, higher electron densities are reached with the ISCT200-MS, up to 4 × 10 16 m -3 . This difference with the ISCT200-GEO, operated at the same discharge conditions, can be explained by the downstream shift of the ionization and acceleration position [51] of the shielded ISCT200-MS.

Electron temperature is plotted on the second row of Figure 10. Much higher electron temperatures are found with the ISCT200-MS Hall thruster, reaching as high as 75 eV. Recent simulations also report much higher electron temperatures with a magneticallyshielded thruster than with a standard magnetic configuration, a feature associated with a downstream shift of the axial profile [START_REF] Garrigues | [END_REF]. Nevertheless, absolute values of electron temperature from Ref. 66 are lower: only 30 eV is reached with the standard configuration, while and 53 eV is attained with the MS configuration. The difference between the electron temperatures from simulations and experiments under similar conditions is yet to be understood, but it is possible that the contribution of plasma instabilities -not captured in the simulations -may contribute to electron heating to some degree. The effect of instabilities such as the electron cyclotron drift instability and the modified two-stream instability in electron heating has been demonstrated in several Hall thruster simulations in recent years, including in Refs. 14, 17, 24, 67, 68. Drift velocities for the two architectures are presented in the third row of Figure 10. In the magnetic shielding configuration, the spectral shift of the Thomson spectra at positions near the exit plane is much more significant than for the standard architecture and the measured drift velocity can be better interpreted.

For previous measurements on a planar magnetron [46], we reported sensitivity for electron drift velocity estimation as low as the velocity corresponding to the Doppler shift associated to a spectral width of one pixel (at the corresponding wavelength). This sensitivity would correspond to approximately 100 km s -1 with the 600 lines mm -1 grating and the 8-pixel wavelength binning used for this paper.

Nevertheless, the Thomson scattering signal obtained from the Hall thruster plasma is more challenging to analyze because of the lower signal-to-noise ratio and distortion of the Thomson spectrum by plasma emission lines. For this reason, electron drift velocities below 500 km s -1 obtained in this campaign were not considered for data discussions. This lower limit corresponds to the drift velocity associated with the Doppler broadening of the FWHM of the instrument profile. In the figures where the estimated drift velocity is shown, a green strip is drawn in order to indicate this threshold. The electron drift velocity direction and magnitude along the azimuthal direction are expected to be easier to predict because the E × B drift should be one of the dominant contributions in the acceleration region. Given the magnetic and electric field orientations in the acceleration region, the azimuthal drift velocity of electrons is expected to be anti-clockwise with the ISCT200-GEO (facing the ion beam) and clockwise for the ISCT200-MS. Along the radial direction, it is more difficult to predict the drift velocity trend as other contributions, such as curvature and gradient drifts, may be on a similar order. With the ISCT200-GEO, the electron drift velocity along the azimuthal direction is systematically positive with only a few values above the threshold of 500 km s -1 . As expected, these positive values correspond to electrons drifting in the clockwise direction. Along the radial direction, the drift velocity is weaker and oscillates around zero.

With the ISCT200-MS the azimuthal drift velocity is significantly negative close to the thruster exit plane. These negative values are due to the reversed magnetic field direction and are consistent with a clockwise E× B drift. Such significant values of electron drift velocity are observable only with the ISCT200-MS because the acceleration region, where the electric field is strong, is shifted downstream, rendering it accessible to the ITS probing beam. As with the standard configuration, electron drift velocities in the radial direction are comparatively low.

Magnetic field direction

To gain further insight into the electron drift velocity values, investigations were performed on the same ISCT200-MS thruster with the magnetic field directed towards the outer, and then inner, poles. Similar discharge conditions used for comparison with the ISCT200-GEO in Section 3.2.2 were used. With A c c e p t e d M a n u s c r i p t the magnetic field directed towards the outer pole, a discharge voltage of 250 V, with 1.14 A discharge current for a xenon flow rate of 1.2 mg s -1 , was used. With the magnetic field directed towards the inner pole, the same discharge power was used, at 250 V with a xenon flow rate of 1.05 mg s -1 . This 15 % increase in the required propellant flow rate suggests the influence of the magnetic field direction on the thruster performance is not negligible.

As shown in the two first rows of Figure 11 electron temperature and density profiles are only weakly affected by the magnetic field direction. In contrast, the electron drift velocity is significantly affected. Along with both the azimuthal and radial directions, the electron drift velocity switches in sign when the direction of the magnetic field is switched. When the magnetic field direction is reversed, the magnetic field gradient and curvature are expected to remain the same, with only the E × B contribution to the electron drift velocity expected to change.

The fact that the measured drift is of nearly the same magnitude regardless of the magnetic field direction suggests that the E × B drift dominates the diamagnetic and curvature drifts near the exit plane. This observation might differ at distances further downstream [69]. The maximum electron drift velocity can be used to estimate the maximum electric field in the acceleration region of the thruster. For a magnetic field of approximately 15 mT where electron drift velocity reaches its maximum value of 1800 km s -1 , the electric field would be 27 kV m -1 . LIF measurements performed on the same thruster at lower discharge voltage (200 V) and 1 mg s -1 of xenon flow gave a slightly higher maximum value (38 kV m -1 ) in Ref. 70.

Discussion on the high electron temperatures observed

The results from these experiments show electron temperatures easily reaching 60 eV for certain conditions, a result which does not match any typical probe measurements on large and small thrusters, where authors have previously measured temperatures in the range of few tens of eV in the near-field region, for example, in Refs. 25, 28 and in Ref. 71, where for a similar discharge voltage in a larger thruster, the electron temperature is only 30 eV in the acceleration region.

We first consider examples of temperatures in the 60 eV range, using spectra obtained on the ISCT200-MS 2 mm from the exit plane with the same discharge condition as Section 3.2.3. These spectra are shown in Figure 12 and concern magnetic fields directed towards and away from the inner pole. The corresponding electron properties, deduced from Gaussian fits to these spectra, are shown in the Figure legend.

The following observations may be made: (i) As discussed in Section 3.2.3, the presence of a strong drift shifts the center of the spectrum away from the laser frequency, and a switch in the magnetic field direction (from directed outwards to inwards) also switches the direction of this shift. This is clearly visible in Figure 12. (ii) In addition to the spectral shift, the spectra are asymmetric. A long tail is visible along the spectral direction associated with the electron drift velocity, while the opposite side of the spectrum shows a fast decrease of signal intensity. These characteristics are consistent with non-Maxwellian EVDFs, and therefore, the presence of electrons which are not in thermal equilibrium. While a Gaussian fit may not be strictly relevant, it can still reasonably fit one side of the spectrum, and the electron temperature is estimated from such a fit. The spectra shown, though much noisier than the example provided in Figure 5 and Figure 6, are still sufficiently well-resolved to show their wide spectral range. (iii) Lastly, we consider challenges associated with the analysis of these spectra. The high-electron temperature, low-electron density regime of the cases shown results in low-amplitude, broad Thomson spectra. These spectral features, combined with the presence of several residual emission lines, contribute to the noisiness of the data. In spite of this, the large spectral width of the cases shown is unambiguous. What contributes to this relatively large width? We may speculate that this observation is due to oscillatory discharge behavior, averaged out during the 10-minute acquisition time for each spectrum. Although we have sought to limit these investigations to regimes in which the discharge current standard deviation is low, this does not rule out the possible influence of high-frequency plasma instabilities [START_REF] Adam | [END_REF] which may generate electric field fluctuations. The accelerating potential is not truly stationary. This may create oscillations in the bulk distribution which, when averaged, would appear as unusually broad spectra, leading to apparently higher electron temperatures than would be expected from the average applied potential. To avoid this, it would be necessary to obtain a highly time-resolved measurement (ideally, a single shot measurement with shorter laser pulse duration), at present beyond reach for these low plasma densities. In addition to this effect, high-frequency instabilities have been demonstrated to contribute to true electron heating in thrusters. This dynamic heating is discussed in a number of works, including in Ref. 67, where electron temperatures reaching the 100 eV range have been mentioned. It is worth and radial (points in orange) directions for the shielded ISCT200-MS thruster. The discharge conditions are fixed to at 250 V and 1.14 A for both magnetic field directions: towards the outer pole (left column) and towards the inner pole (right column).

noting that these simulations do not fully account for the energy dissipation mechanisms for electrons. Nevertheless, a full understanding of the contribution of plasma oscillations to the experimentally-measured spectral widths may be achievable through numerical modeling, to be examined in future work. Further insights into electron behavior may be obtained from the calculation of the electron velocity distribution function from the Thomson spectrum profile [72]. However, this would require a higher signal-to-noise ratio, as discussed in Ref. 45, obtained for example under discharge conditions with higher electron densities.

Influence of the discharge voltage

Lastly, a study of the influence of the discharge voltage on the electron properties was performed, at a fixed position 2 mm from the exit plane. Observations were made along the azimuthal direction. The discharge voltage was varied from 112 V to 300 V at constant discharge power of 300 W, with an adjustment of the xenon mass flow rate from 2.3 mg s -1 to 0.9 mg s -1 . The investigations were performed on the ISCT200-MS with the magnetic field directed toward the external pole. At the axial position probed, electron drift velocity is expected to be dominated by the E × B contribution, giving an anti-clockwise rotation of electrons.

As seen in Figure 13, the electron density is significantly lower at higher discharge voltages. This trend might be expected from the decrease in xenon mass flow rate imposed in order to maintain the discharge power at 300 W when the voltage is increased.

Provided that the accelerated ions have undergone most of the voltage drop at the position probed, their axial velocity can be expressed as v i ≈ qU d /m i , where U d is the discharge voltage and q is the elementary charge (with only singly-charged ions considered).

Substituting this expression of the velocity into Equation 2, under the hypothesis of a constant ionization efficiency (η prop ) and divergence angle (i.e., A is constant), the product n e √ U d will be proportional to the xenon mass flow rate ṁn .

A constant discharge power is equivalent to a fixed value of the product I d U d , and equivalent to a constant value of ṁn U d , since I d ∝ ṁn . The constant discharge power can also be expressed as n e U 3/2 d , using the substitution of n e √ U d for the mass flow rate in ṁn U d . We can refer to the constant discharge power as c. The value of c is slightly lower at 300 V than at 112 V, whereas it would be expected to be constant. This observation is likely consistent with a slight axial shift in the accelerating electric field profile (meaning that the measurement, conducted at a fixed axial position, misses the peak density), or potentially, an increase in the propellant efficiency [START_REF] Grimaud | Design and characterization of a 200W Hall thruster in "magnetic shielding[END_REF] as the voltage is increased. Measurement of the full axial profiles of the electron properties at several different voltages would allow for more detailed interpretations. The relatively large electron densities obtained with the ISCT200-MS at low discharge voltages facilitate the analysis of the Thomson signal. Near the exit plane, quantitative deviations from the typical Gaussian form of the Thomson spectra were also observed at low voltages and will be analyzed and discussed in a forthcoming paper.

The variation in electron temperature shown in Figure 13 is non-monotonic. At low discharge voltage, electron temperature first rises linearly from 8 eV (see Figure 6) to 13 eV before a fast increase to 25 eV at 217 V, followed by a slow decrease to around 15 eV at higher voltages.

Previous studies [START_REF] Raitses | Characterization of plasma in a Hall thruster operated at high discharge voltage 41st[END_REF] using probes investigated the influence of channel width and discharge voltage on the maximum electron temperature. It was reported that regardless of the channel width, a linear increase of the maximum electron temperature was observed at low discharge voltages, with saturation in temperature A c c e p t e d M a n u s c r i p t observed around 60 eV at large discharge voltages.

The trend of the electron temperature can best be interpreted by considering the effect of varying voltage on the position of the peak electric field. For the ISCT200-MS, there is a tendency for the electric field get broader and shift further outwards as the discharge voltage increases, at least for the voltage range (200 to 300 V) considered in Ref. 51. This trend contrasts with some other observations, including on a larger shielded thruster [START_REF] Hargus | [END_REF]76], where the electric field instead shifts upstream as voltage increases. If the electron temperature is correlated with drift velocity, as expected [START_REF] Garrigues | [END_REF], then the peak in the temperature may appear to be situated at 217 V because for this voltage, the peak electric field coincides with the observation volume position 2 mm from the exit plane. The subsequent decrease in the electron temperature at higher voltages would merely be due to the peak electric field moving further downstream of the observation volume position. This explanation would be valid even if the electric field instead shifts upstream as voltage increases; the key point is that a voltage change produces shifts in the electric field not necessarily captured by a fixed observation volume position. A more complete experiment would involve investigations of the electron properties at several axial positions for each voltage.

As shown in Section 3.2.3, at a position 2 mm from the exit plane, the electron drift velocity is dominated by the E × B contribution, and thus the electric field is proportional to the azimuthal drift velocity. This accounts for why the trend in electron drift velocity shown in Figure 13 is similar to that of the electron temperature, which, as we have discussed, depends on the electric field. As with the electron temperature, the drift velocity profile is a reflection of the progressive axial shift of the peak electric field into, and then out of, the fixed observation volume.

The temporal profiles and power density spectra of the discharge current for a number of discharge voltages are shown in Figure 14 and Figure 15. The discharge current main oscillations are in the 2 × 10 5 Hz range at higher discharge voltage and may be consistent with transit time oscillations. At lower discharge voltages, oscillations in the discharge current are more concentrated in the 4 × 10 4 Hz (breathing mode) range. At the lowest discharge voltage of 112 V, current oscillations are significant with a peakto-peak amplitude as large as the mean current, with rotating spokes expected to appear [START_REF] Mazouffre | Plasma Sources Science and Technology[END_REF][START_REF] Matyash | [END_REF]. With this discharge mode, the power spectra of the oscillations cover a broadband in the 10 4 Hz range, but with welldefined oscillations, synchronization of the Thomson signal acquisition with different points in the discharge current period appears achievable, giving access to time-resolved measurements of the electron properties.

At an intermediate discharge voltage, near 217 V, there is a minimum in the power density integrated over the entire frequency range, yet the measured electron temperature is maximum for this voltage. As discussed previously, it is not yet possible to rule out the contributions of high-frequency instabilities (not visible in the discharge current measurement) to the apparent high electron temperatures. Analysis of the discharge current frequency content can only capture oscillations in the ion (and to a smaller degree, electron) current.

The frequency of the discharge current oscillations depends on the charged particle flow and does not necessarily reflect the actual frequency of instabilities involved in electron heating. Azimuthal instabilities in the MHz frequency range, for example, such as those identified in experiments in Refs. 79 and 34, could potentially appear in the discharge current frequency spectrum, provided (i) they produce ion or electron current fluctuations at the same frequencies as the instabilities themselves, and that (ii) these fluctuations in current are large enough to be identifiable in the frequency spectrum of the discharge current. However, given the complex way in which azimuthal instabilities contribute to electron heating and transport 14, 21, the discharge current frequency content alone is not a sufficiently sensitive diagnostic for their presence and contribution to the electron dynamics.

Overall, complex modifications of the discharge properties when the discharge voltage is varied are expected.

More detailed studies (experimental, possibly with coherent Thomson scattering, and numerical, possibly with particle-in-cell simulations) of the contributions of different instabilities to electron heating are required.

In addition, the previous points regarding the shift of the electric field with voltage and the possible influence on the measurement can be clarified using complementary laser-induced fluorescence measurements of the electric field profile.

The discharge conditions at 249 V (300 W) investigated during the voltage variation are similar to those at 250 V (285 W) investigated during the axial exploration discussed in Section 3.2.2 and Section 3.2.3. However, at 249 V the azimuthal electron temperature at 2 mm (circled in Figure 13) is only about 14 eV, far below the 40 eV azimuthal temperature at the same position at 250 V (circled in Figure 11).

To explain this difference, we note first that the 249 V (300 W) discharge conditions correspond to a 5 % increase of the discharge power (from 285 W to 300 W). Given the 200 W nominal power of the thruster, it is not surprising that far from this power, small excursions in the discharge conditions induce changes to the averaged plasma properties such as electron Figure 16: Power density spectra of the discharge current of the ISCT200-MS thruster associated with two similar discharge conditions. The 250 V (285 W) conditions were used for axial explorations in Section 3.2.2 and Section 3.2.3, while the 249 V (300 W) conditions were obtained during the discharge voltage exploration (at a fixed power) of this section.

A c c e p t e d M a n u s c r i p t temperature. This point can most easily be illustrated by considering the discharge current oscillations for the 249 V (300 W) and 250 V (285 W) conditions.

Figure 16 shows the power density spectrum of the discharge current for the two discharge conditions. In the 10 kHz to 100 kHz frequency range, pronounced oscillations are present in the discharge current at 249 V (300 W) which are absent or attenuated at 250 V (285 W). The discharge conditions, though similar in voltage, show differing behavior which could partly account for the large difference in the measured electron temperatures for these conditions.

Conclusion

In this paper, implementations of the recentlydeveloped Thomson scattering diagnostic known as THETIS to Hall thruster investigations were presented. These studies focused on the determination of electron properties in two architectures of a 200 W thruster: a standard configuration and unshielded configuration. Electron properties -temperature, density, and drift velocity -were measured along both the azimuthal and radial directions.

Investigations of electron properties on the standard Hall thruster with different magnetic field intensities revealed axial profiles of electron densities in the 10 16 m -3 range and temperatures around 25 eV. Electron properties measured on the magneticallyshielded equivalent thruster confirmed (i) the axial shift of the plasma in comparison to the standard architecture, seen in higher densities downstream of the exit plane, and (ii) large and significant non-isotropic electron drift velocities, up to 1800 km s -1 near the thruster exit plane. Reversal of the magnetic field was found to produce a reversal of the dominant E × B drift component at the exit plane. Discharge voltage explorations performed at a fixed axial position and discharge power revealed electron properties consistent with a shift in the electric field profile with voltage, however, complementary measurements of the axial electric field profile at different voltages are required to fully understand the influence and degree of such shifts. Large electron temperatures are apparent in many cases (exceeding what might be expected from the average applied potential), an observation which may be linked to plasma oscillations and electron heating by plasma instabilities. Combined studies (experimental and numerical) of plasma instabilities and their possible contribution to these observations will be key.

Figure 1 :Figure 2 :

 12 Figure 1: Hall thruster models used for ITS investigations

Figure 3 :

 3 Figure 3: Schematic view of the THETIS diagnostic setup for Hall thruster investigations (not to scale). The transmission branch is indicated with a dark green line and the detection branch with a light green line. L0, L1, L2 and L3 indicate lenses; A0, A1, A2, and A3 indicate apertures; VBG indicates the notch filter. The plasma volume is shown in blue.

Figure 4 :

 4 Figure 4: Scattering configurations (not to scale) used for investigation of electron properties along the azimuthal (left) and radial (right) directions. Both views are presented facing the thruster plume (in blue). The laser beam is in green and the detection volume in orange.

Figure 5 :

 5 Figure 5: Thomson spectrum example obtained along the azimuthal direction with the ISCT200-GEO thruster. The corresponding electron density and temperature are ((2.6 ± 0.2) × 10 16 m -3 and (20 ± 2) eV respectively. Discharge parameters are a xenon flow rate of 1.6 mg s -1 , discharge voltage of 250 V and a discharge current of 1.13 A. This spectrum is obtained for an axial distance from the exit plane of 4 mm and peak magnetic field of 25 mT at the exit plane.

Figure 6 :

 6 Figure 6: Thomson spectrum example obtained along the azimuthal direction with the ISCT200-MS thruster. The corresponding electron density, temperature and drift velocity are (1.3 ± 0.1) × 10 17 m -3 , (3.7 ± 0.3) eV and(93 ± 12) km s -1 , respectively. Discharge parameters are a xenon flow rate of 2.3 mg s -1 , discharge voltage of 112 V and a discharge current of 2.67 A. This spectrum is obtained for an axial distance from the exit plane of 2 mm and a peak magnetic field of 15 mT at the exit plane.

Figure 7 :Figure 8 :

 78 Figure 7: Comparison of electron density (first row) and electron temperature (second row) obtained with the ISCT200-GEO thruster at 250 V and 1.14 A for different magnetic field intensities. From left to right, the columns correspond to the B low , B mid and B high configurations.

Figure 9 :

 9 Figure 9: Power density spectra of the discharge current for the B low , B mid and B high magnetic field intensities, for the conditions of Figure 8.

0Figure 10 :

 10 Figure 10: Comparison of electron density, temperature and and drift velocity along both azimuthal and radial directions. The discharge conditions are fixed to at 250 V and 1.14 A both magnetic field configurations: standard ISCT200-GEO (left column) and unshielded ISCT200-MS (right column)

Figure 11 :

 11 Figure11: Comparison of electron density, temperature and drift velocity along both azimuthal (points in blue) and radial (points in orange) directions for the shielded ISCT200-MS thruster. The discharge conditions are fixed to at 250 V and 1.14 A for both magnetic field directions: towards the outer pole (left column) and towards the inner pole (right column).

Figure 12 :

 12 Figure12: Thomson specta obtained along the azimuthal direction 2 mm from the exit plane of the ISCT200-MS thruster, for 250 V discharge voltage and 1.14 A discharge current (same conditions as Section 3.2.3). Top: magnetic field directed towards the outer magnetic pole; n e =(3.8 ± 0.3) × 10 16 m -3 , T e =(41 ± 4) eV and v e,drif t =(1000 ± 100) km s -1 . Bottom: magnetic field directed towards the inner magnetic pole; n e =(3.5 ± 0.2) × 10 16 m -3 , T e =(58 ± 7) eV and v e,drif t =(-950 ± 130) km s -1 .

Figure 13 :

 13 Figure13: Electron density, temperature and drift velocity along the azimuthal direction 2 mm from the exit plane for the ISCT200-MS thruster. Various discharge voltages were investigated at a fixed power of 300 W.

Figure 14 :

 14 Figure 14: AC component of the discharge current for different discharge voltages for the ISCT200-MS thruster. The DC component of the current is 1.14 A.

Figure 15 :

 15 Figure 15: Power density spectra of the discharge current for different discharge voltages for the ISCT200-MS thruster, for the same conditions as Figure 14.

  A c c e p t e d M a n u s c r i p t A c c e p t e d M a n u s c r i p t A c c e p t e d M a n u s c r i p t

Acknowledgments

This work was supported by the French Space Agency (CNES) under the grant R & T R-S18/PF-0005-105. The authors would like to thank C. Boniface, T. Minea, J. Fils, J. Schwindling, O. Tuske, and G. Largeau.


Amjad Mouhaidali

Damien Tromeur-Dervout

Olivier Chadebec

Jean-Michel Guichon

Sebastien Silvant

J-M Guichon

Electromagnetic Transient Analysis of Transmission line based on rational Krylov approximation

Keywords: Electromagnetic transient, frequency dependency, rational krylov approximation, Vector fitting, transmission line model, passivity

This paper shows a comparison between Vector fitting and rational Krylov fitting techniques for the determination of rational models concerning the fitting accuracy, the computational performances and the model order. Primarily, the mathematics behind the second technique are presented. It should be noted that rational Krylov fitting have never been used in transmission line modeling. A new procedure is proposed to use rational Krylov fitting instead of vector fitting in the universal line model (ULM). Furthermore, it is demonstrated that this procedure has several advantages over the traditional one. Two illustrative examples involving a transmission system are presented for validation of the new procedure.

I. INTRODUCTION

In the frequency domain and for electromagnetic transient studies, rational approximation is used to reproduce a given frequency behavior. This approximation leads to recursive convolution and fast computations in the time domain [START_REF] Semlyen | Fast and accurate switching transient calculations on transmission lines with ground return using recursive convolutions[END_REF]. Examples of applications are the modeling of transmission lines and cables [START_REF] Morched | A universal model for accurate calculation of electromagnetic transients on overhead lines and underground cables[END_REF], power transformers at high frequencies and the Frequency-Dependent Network Equivalents (FDNE) [START_REF] Gustavsen | Simulation of transmission line transients using vector fitting and modal decomposition[END_REF]. A general rational representation of a frequency dependent function f (s) has the form:

f (s) ≈ r(s) = N m=1 c m s -q m + d (1)
where s = jω in rad/s, the model order N; the poles q m and residues c m are either real or come in complex conjugate pairs; d is a constant. In the case of cable modeling, the matrix f corresponds to the frequency dependent characteristic admittance Y c (s) or the propagation function H(s).

The vector fitting VFIT [START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF], [START_REF] Gustavsen | Improving the pole relocating properties of vector fitting[END_REF] technique is widely used by the engineering community for the achievement of rational models. VFIT has a relatively simple formulation, high computational efficiency and delivers high model accuracy. Recently, VFIT was combined with frequency partitioning and model order reduction via balanced realizations to approximate FDNE [START_REF] Ramirez | Vector fitting-based calculation of frequency-dependent network equivalents by frequency partitioning and model-order reduction[END_REF] or transmission lines and cables [START_REF] Cervantes Martinez | Accurate and reduced order identification of propagation function for electromagnetic transient analysis of cables[END_REF]. A novel approach is proposed in [START_REF] Morales | A loewner/mpm-vf combined rational fitting approach[END_REF], where the Matrix Pencil Method (MPM) and the Loewner Matrix (LM) are incorporated into the VFIT technique.

A new rational approximation method has emerged recently and it is based on the Rational Kylov fitting RKFIT [START_REF] Berljafa | Rational Krylov Decompositions: Theory and Applications[END_REF]. This method is based on Rational Arnoldi Decomposition (RAD) for poles relocation. The used orthonormal RADs deliver a well-conditioned basis for the numerical linear algebra problems, and a more robust and faster convergent algorithm. The main difference between VFIT and RKFIT is the poles relocation method. Both methods solve a Least Square (LS) problem at each iteration, however VFIT uses a projection space represented in the partial fraction basis and RKFIT uses a projection space represented in discrete-orthogonal rational functions.

In this paper, the advantages of RKFIT are highlighted regarding the automatic determination of the model order (N ) with a built-in model order reduction technique, the fitting accuracy and the pole relocation convergence. It is proposed to use RKFIT in the implementation of the Universal Line Model (ULM) [START_REF] Morched | A universal model for accurate calculation of electromagnetic transients on overhead lines and underground cables[END_REF] instead of VFIT.

Unfortunately, RKFIT cannot guarantee the passivity of the model. If passivity violation is detected, post fitting passivity algorithms are necessary. Various passivity enforcement technique are presented in [START_REF] Grivet | Passive macromodeling: Theory and applications[END_REF].

This paper is organized as follows. The rational Krylov algorithm is presented in Section II. Section III assess the numerical performances through different examples. Section IV demonstrate the application of the new approach for a single overhead line and an underground cable system.

II. RATIONAL KRYLOV FITTING ALGORITHM

A. Mathematical background of RKFIT

RKFIT is an iterative algorithm for solving nonlinear weighted rational least square problems [START_REF] Berljafa | Rational Krylov Decompositions: Theory and Applications[END_REF]. Let's consider The following matrices and vector:

• A = diag(s i ) ∈ C Ns,Ns , where N s is the number of frequency samples • F = diag(f i ) ∈ C Ns,Ns with f i = f (s i ) • b ∈ C Ns ,
where b is the starting vector. The algorithm tries to approximate F b as r(A)b, where r is a rational function of the form:

r = p m+k q m (2)
where p m+k and q m are polynomials of degree less or equal to m + k and m respectively such that the relative error ε er is minimal:

ε er = ||F b -r(A)b|| 2 ||F b|| 2 → min (3)
The RKFIT algorithm is extended to approximate a family of matrices all sharing a common denominator q m [START_REF] Berljafa | Rational Krylov Decompositions: Theory and Applications[END_REF]. This feature fits all the element of Y c (s) with a common pole set. The rational Krylov fitting is a fitting method based on Krylov space. We will provide a brief overview of the Krylov methods. A more exhaustive analysis can be found in [START_REF] Liesen | Krylov subspace methods: principles and analysis[END_REF], [START_REF] Saad | Numerical methods for large eigenvalue problems: revised edition[END_REF]. Given the matrix A and the vector b, the corresponding order-m Krylov subspace is defined as:

K m+1 (A, b) def = span {b, Ab, . . . , A m b} (4)
The starting point in the RKFIT algorithm is to consider a non zero polynomial q m ∈ P m that has no roots in Λ(A), the set of eigenvalues of A. The rational Krylov space of order m associated with (A, b, q m ) is defined as:

Q m+1 (A, b, q m ) def = q m (A) -1 K m+1 (A, b) (5)
To extract information out of the rational Krylov space, the first step is to construct the rational Arnoldi decomposition of the space Q m+1 (A, b, q m ) as follows:

AV m+1 K m = V m+1 H m (6)
where the columns of V m+1 are an orthogonal basis of Q m+1 (A, b, q m ). K m and H m are unreduced upper-Hessemberg (m+1)×m pencils such that {h j+1,j /k j+1,j } m j=1

are the root of q m . The rational approximation associated with the rational Arnoldi decomposition in (6) satisfies:

r(A)b = V m+1 K m f K † m H m V m+1 K m † b (7)

B. Pole relocation process

The pole relocation process aims to find new poles that would provide a better approximation to F b from previous ones. Two linear spaces in C Nc are used, where the search space is S = R(V m+1 ) = Q m+1 (A, b, q m ) and the target space T = R( Vm+k+1 ) = Q m+k+1 (A, b, q m ). Let's consider the orthogonal projection onto T as P T = Vm+k+1 V * m+k+1 . Then we can write that r(A)b should be the orthogonal projection of F b onto the target space T :

||F b -r(A)b|| 2 = ||F b -Vm+k+1 V * m+k+1 F b|| 2 = ||(I -P T )F b|| 2 (8)
Then we search for vector v ∈ S such that:

v = argmin v∈S ||v||2=1 ||(I -P T )F v|| 2 (9)
Defining the matrix S = F V m+1 -Vm+k+1 ( V * m+k+1 F V m+1 ), a solution is given by v = V m+1 ĉ, where ĉ is the right singular vector of S associated to the smallest singular value σ min . Let a nonzero polynomial qm ∈ P m that has no roots in Λ(A), then:

Q m+1 (A, b, q m ) = Q m+1 (A, qm (A)q m (A) -1 b, qm ) (10)
The pole relocation consists in finding the polynomial qm such that:

v = V m+1 c = qm (A)q m (A) -1 b (11)

C. Partial fraction basis calculation

The approximant r is calculated via a coefficient vector c ∈ C m+k+1 such that:

c = V * m+k+1 (F b) ||b|| 2 (12)
The computation of lead to the evaluation of the relative error and compared it to the desired approximation tolerance. The full development of the rational Krylov fitting process is proposed in [START_REF] Berljafa | Rational Krylov Decompositions: Theory and Applications[END_REF], [START_REF] Berljafa | Generalized rational krylov decompositions with an application to rational approximation[END_REF]. The high-level description of RKFIT is presented in Fig 1 .   1: Takes initial q m 2: repeat 3:

Set search space S := Q m+1 (A, b, q m ) 4:

Set target space T := K m+k+1 (A, q m (A) -1 b)

5: Find v = argmin v∈S ||v||2=1 ||(I -P T )F v|| 2 6: Let qm ∈ P m such that v = qm (A)q m (A) -1 b 7:
Set q m = qm 8: until stopping criteria satisfied 9: construct approximent r

III. ACCURACY COMPARISON OF RKFIT AND VFIT

In [START_REF] Victor | Algorithms for the rational approximation of matrix-valued functions[END_REF] a selection of algorithms for rational approximation are discussed. Table I recapitulates the main RMS error values for various example and different model order studied in [START_REF] Victor | Algorithms for the rational approximation of matrix-valued functions[END_REF]. It was shown that RKFIT outperforms VFIT for rational approximation regarding the approximation accuracy. In this section, we will focus on the main differences of the two methods. First the model order determination, then the accuracy and pole convergence with weight, after that the choice of the starting poles and at the end the system passivity. An example from the literature used as benchmark for transmission line studies will be fitted with both methods. This case study is taken from [START_REF] Gustavsen | Improving the pole relocating properties of vector fitting[END_REF] and represent the 6 × 6 admittance matrix of the distribution network of Fig. 2 measured from nodes A and B. The magnitude of the elements of the admittance matrix are plotted in Fig. 3. We choose to fit only the first column of the admittance matrix. This case is challenging due to multiple resonance peaks. The fitting accuracy is examined by comparing the fitting RMS error calculated as:

RM S error = Ns k=1 Nc n=1 |f f itted n (s k ) -f accurate n (s k )| 2 N s N c ( 13
)
The Pole relocation process can be monitored by computing so-called Hausdorff distance between the pole sets at two consecutive iterations [START_REF] Grivet | Passive macromodeling: Theory and applications[END_REF]. The simulations are performed using a 16 GB RAM computer, i5-9600 3.10 GHz-processor with MATLAB 2019a. The element of the vector are fitted using the relaxed nontriviality constraint of VFIT, provided in the function vectfit3.m, and the unstable poles are forced to be stable by flipping them.

A. Model order determination and reduction technique

In the first iteration of RKFIT, the model order is chosen to be smaller than the number of frequency samples and high enough to approximate the function with a low deviation error. When RKFIT founds an adequately good approximation, it reduces the model order and tries not to degrade the approximation accuracy. The admittance matrix is fitted using RKFIT with reduction technique. The approximation error is presented in Figure 4 for two error tolerances. The maximal model order is fixed to 70 and the starting poles is chosen to be infinite. RKFIT achieve an RMS error below the tolerance at iterations 2 and 5, after which the degree reduction takes place. The denominator degree is reduced to N = 33 and N = 49 for tolerances T 1 and T 2 respectively. After 1 or 2 iterations of pole relocation, RKFIT reaches the fixed tolerance and the algorithm is stopped.

B. Fitting accuracy and pole convergence

In order to compare the fitting accuracy, The model order is fixed for RKFIT and VFIT. The fitting is done with three weighting schemes (no weight, strong inverse magnitude weight, and weak inverse magnitude weight). The weight is common for all the vectors elements. The initial poles are chosen to be complex and linearly spaced. The RMS error obtained by fitting the admittance matrix through both methods for different fitting orders and weighting schemes is given in Table II. This table validates that RKFIT outperforms VFIT for different weighting schemes and orders regarding fitting accuracy except with the case with weak inverse magnitude weight and high model order. However, VFIT is slightly faster.

The convergence of the pole relocation process as function of the iteration count is shown in Fig. 5. The model order is N = 50 and no weight scheme is applied. The Haussdorf distance is presented on the left side of the graph and the RMS error on the right side. It is chosen to restrict the lowest distance to 10 -8 for illustration purpose. This graph shows that RKFIT converges faster than VFIT. After 5 iterations RKFIT reaches the set of poles that delivers the smallest RMS error and the distance becomes zero. For RKFIT, it is remarked the correlation between the distance and the RMS error. When the distance becomes zero, the RMS is at the lowest value. However for VFIT, the RMS error reaches the lowest value when the distance starts to decrease linearly. This observation can be used as the fitting stopping criteria. For the same previous example, if the model order is increased to N = 100, VFIT does not converge. Fig. 6 shows the convergence of the pole relocation process. RKFIT converges within 4 iterations while VFIT diverges, the Haussdorf distance rises linearly and the RMS error oscillates.

C. Choice of the starting poles

One of the main advantages of RKFIT is that the user requires no a priori knowledge of the initial poles location and distribution. Poles can be chosen to be infinite and RKFIT will relocate them. Fig. 7 shows the convergence of the pole relocations for different initial pole distributions. This figure demonstrates that RKFIT converges faster than VFIT for different starting poles. It should be noted that the infinite value is fixed by the ieee arithmetic representation for positive infinity in [START_REF]Ieee standard for binary floating-point arithmetic[END_REF].

D. passivity

For a multi-port electrical network, the passivity is guaranteed if and only if its transfer admittance matrix Y n (s) has a purely positive real part through the whole frequency range [START_REF] Grivet | Passive macromodeling: Theory and applications[END_REF]. Y n (s) is positive definite if and only if its Hermitian part Y H (s) has only positive eigenvalues:

eig(Y H (s)) = eig (Y n (s) + Y * n (s)) (14)
where Y T n (s) is the conjugate transpose of Y n . The admittance matrix is fitted with different model order to assess the passivity violation for RKFIT and VFIT. The model's passivity proprieties are assessed using the half singularity test developed in [START_REF] Semlyen | A half-size singularity test matrix for fast and reliable passivity assessment of rational models[END_REF] via the main function driver RPdriver.m in [START_REF] Gustavsen | Matrix fitting toolbox. Users guide and reference[END_REF]. Table III illustrates the passivity assessment for different model order approximation. It is shown that RKFIT can deliver passive model for relatively optimal model order, However, for high order models passivity violation appears. Despite the fact that passivity violation cannot be anticipated, several passivity enforcement algorithm can solve this issue.

IV. UNIVERSAL LINE MODEL BASED ON RKFIT

In the ULM [START_REF] Morched | A universal model for accurate calculation of electromagnetic transients on overhead lines and underground cables[END_REF], the propagation function H and the characteristic admittance Y c are approximated using VFIT [START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF]. The rational approximations are used for time domain implementation in EMTP-type simulation tools [START_REF] Hermann | Electromagnetic Transients Program: Reference Manual:(EMTP theory book)[END_REF]. As outlined in the previous section, RKFIT gives more accurate results than VFIT. Accordingly, a fitting procedure for Y c and H is proposed in this part using RKFIT instead of VFIT.

The proposed fitting approach has several advantages over the fitting approach implemented in the EMT-softwares:

1) The proposed approach delivers accurate fitting results with low RMS error and guarantees convergence of the algorithm.

2) The model order is automatically calculated by the RKFIT algorithm. The user only defines the RMS fitting error tolerance.

3) The user requires no a priori knowledge of the pole location. The poles can be equal to the ieee arithmetic representation for positive infinity, represented in Matlab by the value Inf.

A. Rational approximation in the ULM with RKFIT 1) Fitting the characteristic admittance: The characteristic admittance Y c is easily fitted since it doesn't have resonance peaks. All the elements of Y c are fitted with the same set of poles expressed in the following form:

Y c (s) ≈ Y c f itted = N m=1 c m s -q m + D (15)
2) Fitting the propagation function: The approximation of the propagation function is more difficult since its elements contain propagation delays. The modal contributions of H are calculated and a proper delay is assigned to each mode. The delay is calculated using the optimization technique proposed in [START_REF] Gustavsen | Optimal time delay extraction for transmission line modeling[END_REF]. The criterion proposed in [START_REF] Kocar | New procedure for computation of time delays in propagation function fitting for transient modeling of cables[END_REF] is used to group repetitive or close delays. Time delays are removed from the modal components and then fitted using RKFIT. The propagation function is then fitted using the poles from the modal approximation. The only remaining unknowns are the residues, and H is approximated as:

H(s) ≈ H f itted = K i=1 Ni m=1 R m,i s -a m,i e (-sτi) ( 16
)
where K is the number of grouped modes, τ i are the collapsed time delay, R m,i is the matrix of residues calculated for the corresponding poles a m,i of the grouped modes and H f itted stands for the rational approximation of H considering the grouped modes. Time domain implementation is based on the equation developed in [START_REF] Gustavsen | Avoiding numerical instabilities in the universal line model by a two-segment interpolation scheme[END_REF] to reduce integration and interpolations errors.

B. Numerical Results

The proposed method is used with different transmission line configurations and the calculated results are compared with the corresponding results from VFIT. In all cases, the characteristic admittance and the propagation function matrices are calculated for the frequency range from 0.01Hz up to 1M Hz with a sampling rate of logarithmically spaced 10 points/decade.

1) Overhead Transmission Line: Let's consider a 30km single conductor over-head line over a lossy ground as shown in Fig. 8 and presented in [START_REF] Gustavsen | Optimal time delay extraction for transmission line modeling[END_REF]. For the fitting of Y c and H, the initial poles are chosen to be infinite and the maximum allowable RMS-error is 10 -3 . As for the propagation function, the time delay is τ = 1.039×10 -4 sec. The magnitude of H is plotted in Fig 11 . The use of the proposed approach reduces the RMS error compared to VFIT and it can be seen on the magnitude plote specially at high frequencies. 2) Underground cable system : Let's consider the cable system in Fig. 12 with parameters given in Table IV. The resulting order for the fitting of Y c is N = 9 for RKFIT and N = 12 for VFIT. Fig. 13 shows the rms error as a function of the iteration number. It is seen that RKFIT gives a faster convergence and a lower order of approximation, although the final result is quite similar. The time delays are calculated and presented in Table V. Three time delays are nearly equal. Then the six modes of propagation are lumped into four distinct modes, and then subjected to rational approximation. The propagation function H is fitted in the phase domain. Fig. 14 shows the deviation of the resulting fitting of the elements of H for VFIT and RKFIT. Clearly, a better accuracy is obtained using RKFIT than VFIT. The maximum residue/pole ratio for RKFIT is 355.5 and for VFIT is 42.27. The maximum residue/pole is lower than the critical threshold of 1000 [START_REF] Kocar | Improvement of numerical stability for the computation of transients in lines and cables[END_REF]. The transient simulation is performed using the obtained model via RKFIT. The step voltage excitation is presented in Fig. 15, where the first core is energized with a unit step voltage at t = 0. Fig. 16 shows the simulated receiving end voltage. The response is seen to be stable for the RKFIT model. The results are compared to accurate results obtained using the Numerical Inverse Laplace Transformation [START_REF] Moreno | Implementation of the numerical laplace transform: A review task force on frequency domain methods for emt studies, working group on modeling and analysis of system transients using digital simulation, general systems subcommittee, ieee power engineering society[END_REF].

V. DISCUSSION Till now, there are no comprehensive convergences analyses for VFIT and RKFIT. However, an example from [START_REF] Berljafa | Rational Krylov Decompositions: Theory and Applications[END_REF], [START_REF] Lefteriu | On the convergence of the vector-fitting algorithm[END_REF] was constructed where VFIT is repellent. Despite this example and Figure 6 in this work, VFIT is being successfully used for various problems. Nonetheless, RKFIT is based on a numerically well-defined computational scheme and have multiple advantages compared to VFIT. The proposed method determines the model order by itself and the reduction technique is implemented to reduce the model order if the fixed rms error allows it. The order reduction technique is not a post processing technique. It is one of the method proprieties and it is mathematically proven by theorem 6.1 in [START_REF] Berljafa | Rational Krylov Decompositions: Theory and Applications[END_REF]. The final set of poles calculated by RKFIT is mathematically proven to be the optimal solution [START_REF] Berljafa | Rational Krylov Decompositions: Theory and Applications[END_REF]. We observe this in the Haussdorff distance when it becomes equal to zero. However for VFIT, the Haussdorf distance decreases continuously and that means that the poles are slightly varying. The main drawback of RKFIT as VFIT is that they don't guarantee the passivity of the obtained model.

VI. CONCLUSION A new method for the rational approximation is studied and compared to well-known vector fitting. The proposed method is applied to transmission line modeling and specially in the universal line model. It was shown that this method delivers more accurate results than VFIT in the frequency domain and stable results in the time domain. The distinct advantage of the proposed method is the ease of the use by the engineers regarding the model order and the starting poles. Jean-Michel Guichon received the Diploma degree in electrical engineering and the Ph.D. degree from Grenoble Institute of Technology, Grenoble, France, in 1998 and 2001, respectively. He is an Associate Professor at the Universite Grenoble Alpes. He is has been with the Grenoble Electrical Engineering Laboratory (G2ELab), Grenoble, France, since 2003, working in the field of electromagnetic modeling for power electronics. He is focused in the development of numerical tools to compute the electromagnetic field and the electrical behavior in low frequency (without propagation).

Sebastien Silvant graduated from the Ecole Polytechnique in Paris and from the EPFL in Lausanne in 2008. After leading innovation in high voltage circuit breakers, managing teams dedicated to numerical simulation software development, Sbastien has been responsible of a R&D group at Supergrid Institute, with the objective of developing models and applying real-time simulation to the control and protection of power systems and converters with a focus on HVDC.

Fig. 1 .

 1 Fig. 1. High-level description of RKFIT

Fig. 2 .

 2 Fig. 2. Power system distribution system

Fig. 3 .

 3 Fig. 3. Magnitude of admittance matrix of the power system

Fig. 4 .

 4 Fig. 4. Model order reduction with RKFIT

Fig. 5 .

 5 Fig. 5. Convergence of pole relocation process by RKFIT against VFIT

Fig. 6 .

 6 Fig. 6. Convergence fo pole relocation process by RKFIT and VFIT

Fig. 7 .

 7 Fig. 7. Convergence of poles relocation process with different starting poles

Fig. 8 .

 8 Fig. 8. single conductor overhead line

Fig. 9 .Fig. 10 .

 910 Fig. 9. Fitting of Yc using RKFIT and VFIT The magnitude function of Y c is shown in Fig. 9. Fig 10 reports the RMS-error from the application of both methods.It is seen that RKFIT and VFIT reach the same RMS-error, however RKFIT converges faster. It is worth mentioning that both methods have the same number of poles at the end.

Fig. 11 .

 11 Fig. 11. Magnitude of the propagation function

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Cross section of three phase 145kV underground cable

Fig. 14 .

 14 Fig. 14. Propagation function H

Fig. 15 .

 15 Fig. 15. Line energization by step voltage.

Fig. 16 .

 16 Fig. 16. Time domain simulation of unit step voltage response at the cable receiving end.

Amjad

  Mouhaidali has an electrical engineering diploma and the Ph.D. degree from Grenoble Alpes University (France) in 2015 and 2019, respectively. Currently, He is a research engineer based at Super-Grid institute, particularly in the HVDC cable systems and junctions research program. His research areas include the modeling of frequency-dependent and thermal effects in cable systems. Damien Tromeur-Dervout received the Ph.D. degree in applied mathematics from University of Paris VI and ONERA, France , in 1993. Since 1994, He is with University Lyon 1, working on research and development activities related to numerical analysis of methods for high performances computing. In september 2015, he joined supergrid institute, program 1: Supergrid Architecture & Systems, for developping numerical methods for the simulation of electrical networks . Olivier Chadebec was born in 1973, Sens, France. He received the Diploma in Electrical Engineering from the Grenoble Institute of Technology and the Ph.D. degrees from the Univ. Grenoble Alpes, France, in 1997 and 2001, respectively. He is currently CNRS senior researcher (directeur de recherche CNRS) and leads the "Models, Methods and Methodologies Applied to Electrical Engineering research group of G2ELab, Univ. Grenoble Alpes, France. He has been a visiting scientist at the Technology Centre of the Federal University of Santa Catarina, Brasil, in 2012-2013. His research activities focus on computational electromagnetics (finite element and integral methods), equivalent magnetic sources identification by solving inverse problems and low magnetic fields measurements. He has co-authored more than 250 papers published in international journals and conference proceedings.

TABLE I RMS

 I ERROR VALUES FOR ALL TESTED EXAMPLES IN[START_REF] Victor | Algorithms for the rational approximation of matrix-valued functions[END_REF]

	Example	Model order	RMS error
			VFIT	RKFIT
	Matrix fitting	10	1.139 × 10 -2	6.246 × 10 -3
	Toolbox [15]	20	6.069 × 10 -3	1.233 × 10 -3
	CD player SLICOT	10	1.683 × 10 1	3.675 × 10 0
	benchmark [16]	20	1.674 × 10 1	9.061 × 10 -3
	ISS model SLICOT	10	6.729 × 10 -4	8.735 × 10 -5
	benchmark [16]	20	3.016 × 10 -4	1.507 × 10 -5
	Nonlinear eigenvalue,	10	1.476 × 10 -4	2.292 × 10 -10
	buckling example [17]	20	1.270 × 10 -11	2.389 × 10 -12
	ISS model with	10	5.557 × 10 -1	4.702 × 10 -2
	additive normally	20	6.280 × 10 -2	3.798 × 10 -2
	distributed noise [16]

TABLE II RELATIVE

 II APPROXIMATION ERROR OF VFIT AND RKFIT

		Weight	Order	RMS error	CPU
		No weight	50 70	1.58 × 10 -4 1.55 × 10 -8	0.87s 1.11s
	VFIT	Strong inverse	50	1.00 × 10 -3	0.85s
		magnitude weight	70	4.43 × 10 -7	1.14s
		Weak inverse	50	2.65 × 10 -4	0.78s
		magnitude weight	70	1.19 × 10 -8	1.16s
		No weight	50 70	9.88 × 10 -5 1.35 × 10 -8	1.04s 1.43s
	RKFIT	Strong inverse	50	3.18 × 10 -4	1.00s
		magnitude weight	70	3.08 × 10 -8	1.45s
		Weak inverse	50	2.02 × 10 -4	1.06s
		magnitude weight	70	1.51 × 10 -8	1.44s

  Table VI presents the fitting rms error of the modal propagation function using VFIT and RKFIT. The Table VI shows that a more accurate fitting is obtained with RKFIT using less number poles than VFIT.