Alfio E Vinci
email: alfio.vinci@cnrs-orleans.fr

Marco R Inchingolo

Stéphane Mazouffre
email: stephane.mazouffre@cnrs-orleans.fr

Jaume Navarro-Cavallé

Ion dynamics in the magnetic nozzle of a waveguide ECR thruster via laser-induced fluorescence spectroscopy

Keywords: magnetic nozzle, laser-induced fluorescence, ECRT, electric propulsion

Xenon ion velocity is mapped in the magnetic nozzle of a circular waveguide electron cyclotron resonance thruster operating at 5.8 GHz by means of laser-induced fluorescence spectroscopy in the near-infrared spectral range. An array of thruster operational parameters are explored to investigate the influence on the acceleration profile and terminal ion velocity. Owing to several mechanisms which broaden the measured spectra, e.g. Paschen-Back/Zeeman effect, inference of the most probable velocity along with the axial kinetic temperature requires full lineshape modeling, especially in the near-field plume and inside the source. Ions are effectively accelerated along the magnetic nozzle, reaching up to 12 000 m s -1

for the lowest neutral pressure tested. A relatively large axial kinetic temperature is observed, typically in the order of 5000 K, which can be attributed to an extended ionization region that overlaps with the acceleration region.

Introduction

A new class of in-space propulsion based on electrodeless discharges is being researched to extend the mission range typically covered by mature technologies, e.g. gridded ion thrusters and Hall thrusters. These devices embed a magnetic nozzle (MN) to enhance the momentum of the ejected ions and come with the promise of increased durability, simpler satellite integration and a current-free plume, relieving from the need of a neutralizer and allowing the extensive use of alternative propellants [1][2][3]. However, these devices are still far from the performance levels achieved by classical technologies [4][START_REF] Takahashi | Helicon-type radiofrequency plasma thrusters and magnetic plasma nozzles Reviews of Modern Plasma Physics[END_REF][START_REF] Takahashi | [END_REF].

The Electron Cyclotron Resonance Thruster (ECRT) is of central interest in this thruster class, as it exhibits interesting performance metrics [7]. The working principle of an ECRT relies on the electron cyclotron resonance to efficiently transfer power from an electromagnetic wave (typically in the microwave range) to the electrons, rising their perpendicular temperature. This phenomenon occurs locally at the resonance layer, where the condition ω ce = ω is met, with ω ce = eB/m e , the electron cyclotron angular frequency and ω the frequency of the externally applied electromagnetic wave [START_REF] Stix | Waves in plasmas[END_REF][START_REF] Bittencourt | Fundamentals of plasma physics[END_REF]. The conversion of electron thermal energy into ion kinetic energy is then conveyed by the plasma expansion in the magnetic nozzle [START_REF] Ahedo | [END_REF]. The discharge and acceleration mechanisms thus require the existence of an externally applied magnetic field.

Two different microwave injection designs have been reported in the literature for ECRTs, the coaxial [7,[START_REF] Jarrige | Characterization of a coaxial ecr plasma thruster 44th[END_REF][START_REF] Correyero | [END_REF][13] and the truly electrodeless circular waveguide design [START_REF] Jahn | Physics of Electric Propulsion[END_REF][START_REF] Miller | [END_REF][16][START_REF] Crimi | Microwave driven magnetic plasma accelerator studies (cyclops) Tech. rep[END_REF][START_REF] Sercel | An experimental and theoretical study of the ECR plasma engine[END_REF]. In the first case, a central conductor is present in the plasma chamber which can be seen as an extension of the coaxial cable delivering the microwaves from the source to the thruster. This component is subject to strong erosion and heat loads and it is absent in the waveguide version. The way microwave power is coupled to the plasma appears to be a relevant factor for electron heating and consequently ion acceleration.

Lower performance levels have been reported for the waveguide configuration [13] for similar operating conditions (input microwave power of 25 W at 2.45 GHz, 1 sccm of Xenon) and thruster dimensions. A better understanding of the reasons behind the performance discrepancy is of great importance since it could lead to improved design strategies. For this purpose, mutual efforts coming from both simulations and experiments are necessary. Thus far, because of the limitations set by intrusive diagnostics, little is known about the ion dynamics and plasma potential distribution inside the discharge chamber region of a waveguide ECRT, the only information coming from non-self-consistent plasma simulations [START_REF] Inchingolo M R | Hybrid pic-fluid simulation of a waveguide ecr magnetic nozzle plasma thruster Space Propulsion Conference[END_REF]. In addition, although the 83 MN physics is common to other similar devices, 84 ion acceleration profiles along the nozzle have never 85 been reported for a waveguide thruster, being this 86 information of key importance for validating future 87 MN models and plasma simulations. 108

The rest of this study is organized as follows: the 109 experimental arrangement is detailed in Section 2; after 110 commenting on the recorded spectra and discussing 111 the procedure to infer the velocity profiles through 112 lineshape modeling, the results are presented in Section 113 3; ultimately, conclusions are drawn in Section 4. The microwave frequency used is 5.8 GHz, allowing for smaller waveguide dimensions than the usual 2.45 GHz frequency, and consequently lowering the size of the device, its mass flow rate and power requirements. A higher frequency comes with the need of a stronger field for reaching the resonance condition (B res = 0.207 T at 5.8 GHz). The field generated by the PM alone suffices the ECR condition inside the PC.

The PC is made of stainless steel. Its diameter and length are D T = 36 mm and L T = 20 mm. The PC is separated from the aluminium circular waveguide at the back by a quartz window sealed with a FFKM oring. Xenon gas is injected radially in the PC from twelve, 1 mm diameter holes located at z = -18 mm on the lateral PC walls.

Microwave power is delivered from the source to the thruster through a waveguide transmission line (TL). More information on the TL used for this work can be found in reference [START_REF] Inchingolo | Submitted to Plasma Sources Science and Technology[END_REF] and [START_REF] Inchingolo | Direct thrust measurements of a circular waveguide electron cyclotron resonance thruster 37[END_REF]. The thruster waveguide and consequently the PC is electrically floating with respect to the rest of the TL. This is made possible by a vacuum-gap DC-block consisting of atoms and ions. In previous works [START_REF] Mazouffre | [END_REF][23][24], the ii) accurately measure the laser wavelength by way of a calibrated wavemeter whose absolute accuracy is 80 MHz (≈ 60 m s -1 ). Modulation of the laser beam power is achieved using a mechanical chopper at ∼ 500 Hz.

The laser is coupled to a 50 µm multi-mode optical fiber which transports the beam from the optical bench to the inside of the vacuum facility. Here a custom collimator is employed to shine the laser beam towards the backplate of the ECRT in the direction of the thruster axis z, see Figure 1. At the output of the collimator, the laser beam diameter is about 4 mm, whereas the power density is maintained to ∼ 2 mW mm -2 , which ensures a weak saturation effect. Note that, in the described configuration, k • v < 0, where k is the laser wavevector and v the ion velocity. The collimator is kept at a safe distance (>20 cm) from the thruster exit plane in order to minimize plume perturbations and avoid overheating and damage of optics. The fluorescence light is collected using a bi-convex lens with 60 mm focal length focusing on a 200 µm core diameter optical fiber. In order to resolve the ion velocity at z < 0, the detection branch is tilted by 60 deg with respect to the thruster axis, cf. Figure 1. Both the excitation and detection branches are mechanically linked to a linear translation stage, whose spatial resolution is 1 mm and accuracy ≈ 300 µm. Any relative displacement between the two optical branches is inhibited. This enables points to be probed in the range of -1 cm < z < 20 cm along the thruster axis.

The 200 µm optical fiber transports the fluorescence light to a monochromator which isolates the selected line from the whole spectrum. Afterwards, a photomultiplier tube converts the light signal into voltage, which is fed into a lock-in amplifier operating at the laser modulation frequency to distinguish the signal from the natural plasma emission. A signal deconvolution is not performed in obtaining the lineshapes.

Xenon ions are probed in the 5d 2 [4] 7/2 metastable level. The selected optical transition is detailed in Table 1, together with the excitation and fluorescence wavelengths [23,[25][26][27][START_REF] Kramida | NIST Atomic Spectra Database (ver. 5.8)[END_REF].

Results and discussion

Observed spectra

Fluorescence spectra are resolved along the thruster axis in a region that extends from z < 0 up to several centimeters downstream the exit plane. The farthest distance being probed is limited by the signalto-noise ratio of the recorded spectra and it depends on velocity. Following the discussion reported in [24], the 279 raw spectra recorded near z = 0 do not image the 280 local ion velocity distribution function (IVDF). The relatively strong magnetic field inside the source and in the near-field plume determines the splitting of the energy levels. This phenomenon leads to a complex lineshape as exemplified in Figure 4. Other LIF studies on ECR thrusters [START_REF] Correyero | [END_REF][START_REF] Wachs | [END_REF] did not observe such a prominent effect on the experimental profiles. This fact has been attributed to the weaker magnetic field within the probed domain [START_REF] Wachs | [END_REF]. In the current experiment, it is found that departure from a quasi-Gaussian profile occurs at field strengths greater than ∼ 30 mT. As a result, modeling of the fluorescence profiles needs to be addressed for the accurate determination of the ion velocity.

It is worth mentioning that the presence of multiple populations ions, e.g. resulting from chargeexchange (CEX) collisions, cannot justify the profiles recorded in this experiment. In fact, this hypothesis would not provide an explanation for the absence of slow ions farther downstream the exit plane.

Lineshape modeling

The purpose of this section is to provide a broad overview of the theory that has been implemented to reconstruct the fluorescence profiles of interest for this experiment. A more detailed treatment can be found in [24] and the references therein included.

The spectrum is modeled accounting for seven stable xenon isotopes having a natural abundance >1 % [30]. Each isotope leads to the so called isotopic shift due to mass and volume differences in the nucleus, whereas those isotopes with odd atomic number ( 129 Xe and 131 Xe) possess a nonzero nuclear spin I and therefore lead to the hyperfine structure (HFS). On top of that, the external magnetic field further splits each spectral line into several components. This effect is accounted for either in the weak or strong field approximation. The former is called Zeeman effect, while the latter corresponds to the Paschen-Back effect. The choice of relying either on the weak or strong field approximation depends on the local strength of the external magnetic field in the context of the hyperfine interaction [START_REF] Woodgate | Elementary Atomic Structure[END_REF]. In practice, the complete modeling procedure consists in:

(i) finding every allowed optical transition;

(ii) quantifying the associated detuning frequency;

(iii) determining the relative intensity of each line;

(iv) accounting for Doppler broadening and Doppler shift.

First, the upper and lower energy levels of the probed transition are solved by determination of the quantum numbers of relevance. For the even isotopes, the electronic angular momentum J and the magnetic M J are required. In the weak field regime, the total angular momentum F and the magnetic M F quantum numbers are needed to identify the hyperfine sub-levels [25,26]. However, in the strong field case, F and M F are undefined quantum numbers therefore the set (I, J, M I , M J ) needs to be used instead. Afterwards, the allowed transitions are set by the values of ∆M F,J . The latter are selected by accounting for the laser polarization with respect to the local direction of magnetic field. In this experiment, the laser beam is parallel to the magnetic field, therefore only the σ components (∆M F,J = ±1) are allowed. The selection rule ∆M I = 0 is additionally considered whenever the strong field condition is met. These selection rules lead to 12 lines for the each even isotope, 36 for 129 Xe and 108 for 131 Xe.

The energy shift of each transition is computed as the sum of three terms: the isotopic shift, the contribution of the hyperfine structure, the contribution of the magnetic field. The values of the isotopic shifts of interest are tabulated in [START_REF] Smith | Diode Laser-Induced Fluorescence of Xenon Ion Velocity Distributions[END_REF]. It is chosen to refer these shifts to 132 Xe since it is the most abundant isotope. The hyperfine structure term is nonzero for the odd isotopes only and it is function of the hyperfine constants A and B [25,[START_REF] Woodgate | Elementary Atomic Structure[END_REF][START_REF] Sargsyan | [END_REF]34].

The magnetic field contribution is a function of the magnetic quantum number M F,J,I , the strength of the magnetic field and the Landé g-factor. Complete expressions are found in [START_REF] Woodgate | Elementary Atomic Structure[END_REF] and the required constants are provided in [START_REF] Kramida | NIST Atomic Spectra Database (ver. 5.8)[END_REF]35].

At this stage, the complete spectrum is built by determining the relative intensity of each line.

The formulation is reported in [START_REF] Woodgate | Elementary Atomic Structure[END_REF] for each ∆M F,J component ‡. These intensities require to be weighted accounting for natural abundance of each isotope [30] and hyperfine components tabulated in [36] as a function of I, ∆J and ∆F . Doppler broadening of each transition is addressed by applying a Gaussian profile exhibiting identical width and frequency shift to each line. Here the laser beam profile is considered as a Dirac delta function due to its bandwidth (< 1 MHz), and the natural width (5 MHz) is neglected in comparison to the Doppler broadening full width at half maximum (in the order of 1 GHz at typical kinetic temperatures) [26].

Eventually, the theoretical fluorescence profile stems from the sum of all individual Gaussian profiles.

The width of each Gaussian component is function of the ion kinetic temperature T i , which is a free parameter in this analysis. It is worth mentioning that the broadening linked to T i represents the spread in the axial ion velocity. As an example, Figure 5 shows ‡ The results are obtained by direct substitution of |JM J ⟩ and |F M F ⟩ in the equations reported in [START_REF] Woodgate | Elementary Atomic Structure[END_REF]. The extent to which these intensities are approximations was not explored here. red evolution with T i . The highlighted point at T i = 4900 K identifies the value that best fits the spectrum recorded at z = 10 mm, whereas T i that minimizes χ2 red of the whole series equals 4700 K. These data are relative to the set of profiles recorded at:

P F = 200 W, ṁ = 2 sccm, Ic = 5 A.
two experimental profiles along with a few theoretical 385 lineshapes computed using arbitrary values of T i for 386 the sake of comparison. In this case, one can observe 387 a good agreement with the experimental data when 388 T i ≈ 5000 K. For a more consistent determination 389 of T i , the algorithm relies on the computation of 390 the reduced-χ 2 metric (χ 2 red ), which quantifies the 391 goodness of fit given the experimental and theoretical 392 lineshapes. For a fixed set of thruster operating 393 parameters, χ 2 red is computed for each spectrum and 394 then averaged over all spectra giving χ2 red . By that 395 means, T i is chosen with a resolution of 100 K as the 396 value that outputs a value of χ2 red closer to 1. An 397 example of how χ 2 red evolves as a function of T i is 398 provided in Figure 6. In doing so, it is intrinsically assumed that T i does not change throughout the plasma expansion. However simplistic, this procedure did not show the need to consider different values of T i as the axial position z varies. In fact, it shall be shown later that the theoretical and experimental profiles are generally found to be in good agreement.

A different approach would consist in computing the ion kinetic temperature that provides the best fit for each individual measured lineshape, so as to yield a function in the form of T i (z) for a fixed set of thruster operation parameters.

For further completeness, an implementation of this kind has been performed within the fitting algorithm. However, as a general remark, no clear trend of T i associated with the axial coordinate is found. Thus, a relatively small (≲ 10%) difference in T i along z is most probably related to concurrent uncertainties, e.g. noise level, quantum constants, robustness of the fitting algorithm, thruster stability at fixed operation parameters, departure from equilibrium. Therefore, this more general analysis suggests that the proposed use of T i (z) = T i can be regarded as a valid approach.

As a final step, once the Doppler broadened line is generated, the most probable ion velocity is estimated from the Doppler shift equation [23]. This shift is equal to the off-set in the frequency axis that gives the best match between the experimental and modeled lineshapes with a precision of 10 MHz (≈ 8 m s -1 ).

The error associated with the ion velocity computation results from several contributing factor, such as the uncertainty on hyperfine constants and isotopic shifts in input, the quality of the raw data, the robustness of the fitting algorithm. By exploring the full range of uncertainties of the input constants to the model (i.e. isotopic shifts and hyperfine constants) and different filtering parameters on the noisy signals, it is estimated that the typical uncertainty in the velocity value is below the absolute accuracy of the wavemeter. show the theoretical and experimental LIF profiles at z = 10 mm and z = 100 mm when the additional coil of the thruster is powered with 0 A, 5 A and 10 A, respectively. Also in these cases, the modeled lineshape is in close agreement with the measurements. A kinetic temperature between 4700 K and 6400 K is found to fit these data, leading to χ2 red ≃ 5 -20, see Figure 6.

It shall be pointed out that the typical T i observed in this experiment is notably larger that what reported in previous studies concerning magnetic nozzle dynamics [24], where T i < 1000 K were found.

At such low plasma pressures typical of the source in analysis here, ion heating due to electron impact collisions cannot explain such high temperatures.

Thermalization is rather a relevant phenomenon for other operating regimes related to VASIMR and MPD thrusters. Conversely, in this case, a large value of T i in the axial direction, i.e. a dispersion in the axial velocity, is rationally attributed to an extended ionization region that overlaps with the acceleration region. This process then would lead to an anisotropy in the ion kinetic temperature between parallel and perpendicular directions. The same phenomenon has been often observed in Hall thrusters, where the spread in axial velocity essentially matches the spread in potential across which ions are produced and it reads up to a few tens of eV [START_REF] Mazouffre | [END_REF]23,[START_REF] Szabo | Fully kinetic numerical modeling of a plasma thruster[END_REF]. However, the dispersion here is noticeably lower than what is typically found in Hall thrusters due to the presence of the magnetic nozzle which enhances the axial 478 electron conductivity, whereby significant potential 479 drops within short distances are unrealistic.

In 480 principle, dispersion in ion velocity could also originate 481 from thruster unstable and/or oscillatory behaviors 482 during a laser frequency tuning scan. The occurrence 483 of these transients was continuosly monitored by 484 measuring the microwave reflected power and V T and 485 data were acquired only during steady state operation. 486

Velocity profiles

487

A parametric approach is followed in performing the 488 present experiment. Measurements are clustered into 489 three main blocks, each of which analyzes the effect 490 due to the variation of a single external parameter: i) 491 an electromagnet current series; ii) a mass flow rate 492 series; iii) an input power series. Results in terms of 493 most probable axial velocity of the ions are respectively 494 reported in Figure 11, 12 and 13. These values are 495 computed from the Doppler shift of the fluorescence 496 spectra as discussed in Section 3.2.

497

The electromagnet current series, see Figure 11, 498 is recorded using constant values of P F = 200 W and 499 ṁ = 2 sccm, whereas the current driven in the exter-500 nal electromagnet is varied between I c = 0, 5, 10 A. In-501 creasing the electromagnet current has three main ef-502 fects: i) the magnetic field strength in the PC in-503 creases; ii) the ECR layer shifts downstream; iii) the 504 expansion rate of the MN decreases or, in other terms, 505 the nozzle becomes longer. It is observed that, at 506 z = -10 mm, the ions are mostly directed towards the 507 thruster backplate (v z ≈ -500 m s -1 ), where they will 508 recombine. This suggests that the peak in plasma po-509 tential is positioned between -10 < z < 0 mm. This 510 result holds regardless of I c . Most of the momentum 511 gain occurs within 5 cm downstream the exit plane, 512 yet at different rates depending on I c . A higher I c is 513 associated with a lower acceleration, which is compat-514 ible with a lower axial variation rate of the magnetic 515 field intensity. Farther downstream instead, any dif-516 ference linked to the shape of the MN is smoothed 517 out and every configuration leads to a nearly identi-518 cal axial velocity profile. This suggests that varying 519 I c does not significantly affect the total plasma poten-520 tial drop in this experiment. A rather constant V T 521 has been observed between the three working points 522 (≈ 70 -80 V), qualitatively confirming that the condi-523 tion in the source were similar. However, it is expected 524 that a non-negligible difference among the various MN 525 shapes might be preserved in terms of plume diver-526 gence, as visually noticed in Figure 3. This aspect, 527 though analyzed in [START_REF] Inchingolo | Submitted to Plasma Sources Science and Technology[END_REF], requires further investigation 528 and experimental data.

529

The mass flow rate series, cf. The farthest distance at which the ion velocity could 532 be resolved is limited by the LIF signal-to-noise ratio.

533

The latter is deteriorated at larger mass flow rates mean free path for electron-neutral and ion-neutral collisions. The first effect leads to a reduced electron temperature, hence a smaller potential drop along the expansion. On the other hand, the electron-neutral collisions contribute to electron cooling inside the source therefore lowering the internal energy available for ion acceleration through the MN, whereas ionneutral collisions represent a drag force component which further reduces the ion kinetic energy [START_REF] Collard | [END_REF].

The input power series, whose results are reported in Figure 13, is recorded at z = 10 cm at constant value of ṁ = 2 sccm and I c = 0 A. Within the explored range of P F , the ion axial velocity along the MN axis shows an approximately linear growth with the microwave input power, although it is expected to saturate as the thrust efficiency decays for higher powers, as shown by direct thrust measurements performed on this device [START_REF] Inchingolo | Direct thrust measurements of a circular waveguide electron cyclotron resonance thruster 37[END_REF]. As opposed to the case shown in Figure 12, increasing the power at constant ṁ rises the available energy per electron, therefore T e increases along with the plasma potential drop. It is noted that the value of v z for P F = 200 W measured in the power series differs by roughly 10 % from the same point shown in Figure 11 and Figure 12. This discrepancy is attributed to a slightly different condition of the thruster operation, since the two measurements are performed over different time frames. As a matter of fact, the reading of the floating potential at the source walls differs by about 15 V from one day to the next. This qualitatively indicates that the conditions in the source were slightly different and most probably a lower plasma potential drop was achieved, giving rise to a difference in the measured ion velocity.

Conclusion

The dynamics of xenon ions has been spatially resolved by way of laser-induced fluorescence spectroscopy along the axis of a circular waveguide electron cyclotron resonance thruster.

In the near-field plume and inside the source chamber, inference of ion velocity and parallel ion temperature requires intensive lineshape modeling efforts. One must account for several broadening mechanisms, namely isotopic shift, hyperfine structure and Paschen-Back/Zeeman effect. This is especially true for ECRTs operating at relatively high frequencies which imply the need for more intense magnetic fields to meet the resonance condition.

Findings show that ions are considerably accelerated within a few centimeters downstream the exit plane, although a non-negligible amount of kinetic energy is gained in the far-field region. When the divergence of the magnetic nozzle is reduced, the accelera-

114 2 .Figure 1 Figure 1 .

 211 Figure1shows a schematic view of the experimental 117 arrangement in the vacuum chamber facility at the 118 University of Carlos III de Madrid (UC3M). The 119 thruster unit is shown at the top and a detailed 120 description is available in[START_REF] Inchingolo | Submitted to Plasma Sources Science and Technology[END_REF][START_REF] Inchingolo | Direct thrust measurements of a circular waveguide electron cyclotron resonance thruster 37[END_REF]. It mainly comprises 121 a circular waveguide constituting the plasma chamber 122 (PC) and a magnetic field generator. The latter, is 123 made up by two separate components: a permanent 124 magnet (PM) and an electromagnet. The PM consists 125 of a radially magnetized ring, which is made of 126 twelve Sm-Co YXG-32 magnets glued together and 127 enclosed in an aluminium box. The brown solid line 128 in Figure2(a) depicts the magnetic field strength 129 produced by the PM along the thruster axis. On the 130 other hand, the electromagnet can be used to adapt the 131 shape of the MN and displace the resonance position. 132 The effects of the electromagnet on the magnetic field 133 strength and divergence are shown in figures 2(a)-(b) 134 and 3.

Figure 2 .Figure 3 .

 23 Figure 2. (a) Measured axial profiles of the ECRT magnetic field amplitude, (b) simulated 2-D maps of the MN highlighting the ECR surface as a function of Ic. z = 0 refers to the exit plane of the discharge tube. The dash-dotted line indicates the resonance at 5.8 GHz.

194

  measurement technique was described in details along 195 with a description of the optical bench used in this 196 study. The probed species is excited using a laser 197 beam produced by an amplified tunable single-mode 198 laser diode in the near-infrared spectral range. The 199 diode can deliver up to 600 mW in the 810 -840 nm 200 wavelength range. In the present experiment, the 201 laser remains mod-hop free over a frequency tuning 202 range of about 5 GHz. By splitting the primary 203 laser beam into multiple beams, it is possible to: i) 204 monitor continuously the laser mode and detect any mode hop by means of a Fabry-Pérot interferometer;

z

  = 15.0 cm 834.71 834.72 834.73 834.74 834.75 834.76 834.LIF signal [a.u.]

Figure 4 .

 4 Figure 4. LIF spectra resolved along z. Wavelength values are expressed for air. Thruster operation parameters are: P F = 200 W, ṁ = 2 sccm, Ic = 0 A.

Figure 5 .

 5 Figure 5. Comparison of experimental (squares) versus modeled (solid lines) lineshapes using arbitrary values for T i . The experimenal data points are Doppler shifted to fit the model profiles. The bar plot images the allowed optical transitions. Thruster operation parameters are: P F = 200 W, ṁ = 2 sccm, Ic = 5 A.

Figure 6 .

 6 Figure 6. Instance of χ 2red evolution with T i . The highlighted point at T i = 4900 K identifies the value that best fits the spectrum recorded at z = 10 mm, whereas T i that minimizes χ2 red of the whole series equals 4700 K. These data are relative to the set of profiles recorded at: P F = 200 W, ṁ = 2 sccm, Ic = 5 A.

Figure 7 .

 7 Figure 7. Experimental (squares) versus modeled (line) lineshape of Xe II. The experimental data points are Doppler shifted to fit the model profile. Thruster operation parameters are: P F = 200 W, ṁ = 1 sccm, Ic = 0 A.

Figure 8 .

 8 Figure 8. Experimental (squares) versus modeled (line) lineshape of Xe II. The experimental data points are Doppler shifted to fit the model profile. Thruster operation parameters are: P F = 200 W, ṁ = 2 sccm, Ic = 0 A.

Figure 9 .

 9 Figure 9. Experimental (squares) versus modeled (line) lineshape of Xe II. The experimental data points are Doppler shifted to fit the model profile. Thruster operation parameters are: P F = 200 W, ṁ = 2 sccm, Ic = 5 A.

Figure 10 .

 10 Figure 10. Experimental (squares) versus modeled (line) lineshape of Xe II. The experimental data points are Doppler shifted to fit the model profile. Thruster operation parameters are: P F = 200 W, ṁ = 2 sccm, Ic = 10 A.

Figure 7 to

 7 Figure 7 to Figure 10 provide different instances of model outcomes compared to experimental spectra.

Figure 7

 7 Figure 7 reports a series of axial scans which cover a broad region spanning from the thruster exit plane up to 20 cm downstream. The model adequately predicts the fluorescence lineshape for all of these very different conditions, i.e. from strong magnetic field at z = 0 up to fairly negligible magnetic field at z = 200 mm. For this series of spectra, a kinetic temperature of 3900 K represents the best fitting value, yielding χ2 red ≃ 2. Figure 8, 9 and 10

Figure 12 ,Figure 11 .Figure 12 .Figure 13 .

 12111213 Figure 11. Most probable axial velocity profile of Xe II for different electromagnet currents. Constant thruster operation parameters are: P F = 200 W, ṁ = 2 sccm.

  In this context, 88 this work examines the ion dynamics in the plume 89 of a waveguide ECRT via near-infrared laser-induced 90 fluorescence (LIF) spectroscopy. The LIF technique 91 has the advantage of being weakly intrusive and 92 capable of high spatial resolution when compared to 93 electrostatic diagnostics, allowing to probe regions 94 otherwise inaccessible without affecting the expansion 95 itself. However, in MN-based discharges, the relatively 96 intense external magnetic field can hinder the easy 97 inference of the ion velocity from the shift of the line 98 center. This is particularly true inside the source and 99 in the proximity of the exit plane, where the field might 100 lead to visible Zeeman and Paschen-Back effects. An 101 accurate estimation of the Doppler shift in this region 102 requires a full modelling of the fluorescence profile 103 by accounting for all broadening mechanisms. As a 104 further result of the lineshape modelling, the axial ion 105

kinetic temperature can be also inferred from a fit of 106 the LIF profiles, offering insights on the ionization-107 acceleration region in the discharge chamber.

Table 1 .

 1 Xe II probed optical transition.

		Transition	Wavelength (air)
	Xe II	5d 2 [4] 7/2 → 6p 2 [3] • 5/2 6p 2 [3] • 5/2 → 6s 2 [2] 3/2	834.724 nm 541.915 nm
	of two waveguides aligned and separated by a distance

This work has been supported by the PROMETEO

634 project, funded by the Comunidad de Madrid, under 635 Grant reference Y2018/NMT-4750 PROMETEO-CM, 636 and it was also carried out in the frame of the 637 Santander Chair of Excellence from the Carlos 638 III University in Madrid granted to Dr. S. 639 Mazouffre in 2020. Additional support came 640 from the ESPEOS project, funded by the Agencia 641 Estatal de Investigación (Spanish National Research 642 Agency), under Grant number PID2019-108034RB-643 I00/AEI/10.13039/501100011033).


Denis Carnier
email: denis.carnier@kuleuven.be

Arthur Correnson
email: arthur.correnson@cispa.de

Christopher Mcnally
email: mcnallyc@mit.edu

Youssef Moawad
email: y.moawad.1@research.gla.ac.uk

Goose: an OCaml environment for quantum computing

In this presentation, we showcase Goose: an OCaml library to model, simulate and compile low-level quantum programs. Goose is designed as a playground with an emphasis on extensibility and accessibility to non-experts. The library is compatible with the Open-QASM standard, and targets a variety of backends with a minimalistic, circuit-based IR.

Introduction

Quantum computing [START_REF] Nielsen | Quantum Computation and Quantum Information: 10th Anniversary Edition[END_REF] is an emerging model of computation that exploits non-classical effects like superposition and entanglement to achieve algorithmic speedups. A radical break from classical computing, the model presents new challenges for programming languages research.

In particular, researchers are interested in more expressive ways to design and develop quantum algorithms, and the ability to leverage existing tooling for classical programming languages in the construction of quantum programming languages (QPLs) [START_REF]Qiskit: An open-source framework for quantum computing[END_REF][START_REF][END_REF][START_REF] Smith | A practical quantum instruction set architecture[END_REF][START_REF] Bichsel | Silq: A high-level quantum language with safe uncomputation and intuitive semantics[END_REF][START_REF] Green | Quipper: A Scalable Quantum Programming Language[END_REF][START_REF] Altenkirch | A functional quantum programming language[END_REF].

OCaml is a mature ecosystem, widely used for programming languages research. However, to the best of our knowledge, it has a scarcity of tools for quantum programming. We therefore present Goose 1 : an OCaml library to model, simulate and compile quantum programs. The library is designed to act as a playground to experiment with quantum programming in the OCaml ecosystem. Figure 1: A high-level overview of Goose's modular pipeline. 1 Goose is available at https://qgoose.github.io.

Figure 1 gives a high-level overview of Goose. The library is built around an intermediate representation (IR) for quantum circuits. This IR can be used as an entry point to, for instance, build simulators, compilers, and static analyzers for quantum programs. Goose includes a full state-based simulator to interpret the IR, designed as a functor parametrized by a linear algebra module. This design choice allows one to implement many tools simply by specializing the simulator. We demonstrate it by deriving both a symbolic simulator and an emitter to C. Further tools could be implemented, for example optimizers and compilers for other languages.

To ensure compatibility with other existing tools related to quantum computing, we provide a frontend for OpenQASM 2.0, a standard low-level quantum circuit representation [START_REF] Cross | Open Quantum Assembly Language[END_REF].

3 An example

|0⟩ H |0⟩ |0⟩ Figure 2: 3-qubit GHZ circuit.
In Figure 2, we provide an example of a GHZ [START_REF] Nielsen | Quantum Computation and Quantum Information: 10th Anniversary Edition[END_REF] circuit entangling three qubits. Listing 3a shows the corresponding Goose IR representation. This is implemented with a list comprehension that specifies the CX gates, with an initial H gate. In Listing 3b, we display the C code that was generated based on the results of the symbolic simulator. cadd, cmul, and csub are C functions for complex number arithmetic. Finally, s and o are arrays for storing the input and output state, respectively. (1, 0, 0, 0, 0, 0, 0, 1) T

Conclusion and future work

We presented Goose: an open source library for quantum computing in the OCaml ecosystem. Currently, the library's front-end supports a subset of the OpenQASM 2.0 circuit representation, a standard for interoperability between quantum computing tools. Additionally, Goose implements a generic simulator for quantum circuits. We specialized this simulator to both a symbolic simulator and a C emitter.

In the future, we want to connect Goose to other existing projects, for instance Twist [10], a QPL whose existing OCaml-based interpreter calls out to C++ quantum simulation libraries, or VOQC [START_REF] Hietala | A verified optimizer for Quantum circuits[END_REF], a formally verified compiler for quantum circuits that can be extracted to OCaml.

  Explicit construction of a GHZ circuit in the Goose IR. cfloat *s = (cfloat*)malloc(N*sizeof(cfloat)); cfloat *o = (cfloat*)malloc(N*sizeof(cfloat)); for (int i = 0; i < N; i++) s[i] = (cfloat) {0.0, 0.0}; o[0] = cmul(SQRT1_2,cadd(s[0],s[1])); o[1] = cmul(SQRT1_2,csub(s[6],s[7])); o[2] = cmul(SQRT1_2,cadd(s[2],s[3])); o[3] = cmul(SQRT1_2,csub(s[4],s[5])); o[4] = cmul(SQRT1_2,cadd(s[4],s[5])); o[5] = cmul(SQRT1_2,csub(s[2],s[3])); o[6] = cmul(SQRT1_2,cadd(s[6],s[7])); o[7] = cmul(SQRT1_2,csub(s[0],s[1])); (b) Sample of generated C code to simulate this circuit for n = 3. Each line computes one amplitude in the output state: 1 √ 2

Acknowledgements

We would like to thank the anonymous reviewers who have helped to improve this paper through their criticism and suggestions. Moreover, we especially thank Paulette Koronkevich for inspiring the name Goose. C.M. was supported by the CQE-LPS Doc Bedard Fellowship. This research is co-funded by a PhD studentship from the College of Science and Engineering at the University of Glasgow, by the Research Fund KU Leuven, by the Flemish Research Programme Cybersecurity and by the European Union (ERC, UniversalContracts, 101040088). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.


Ramon Estalella-Rodríguez

Carlos Olalla
email: carlos.olalla@urv.cat

Angel Cid-Pastor
email: angel.cid@urv.cat

Isabelle Queinnec
email: queinnec@laas.fr

Numerical Methods for the Robust Stability and Performance of Power Conversion Systems

One of the open problems in the design of power converters is related to the selection of components, and the effects that these choices may have, beyond steady-state properties, in the dynamical behaviour of the system. This paper describes a numerical method that may help in finding appropriate components. The method is based on Lyapunov functions that are efficiently manipulated to derive linear matrix inequality (LMI) based conditions. With the proposed approach, the domain of component values that verify given dynamical properties can be estimated. The resulting parametric domain can then be used by the practicing engineer to optimise other factors, such as weight, cost or volume. The advantages of the proposed approach have been illustrated with the design of two different input filters for power converters operating as constant power loads.

I. INTRODUCTION

The design of power filters is an important part of the design of switching power converters. These converters often require input filters to reduce current and voltage ripple and to comply with conducted electromagnetic interference (EMI) regulations. Since input filters are typically designed on top of existing converters, designers must consider the possible interactions that may arise. Such interactions can not only make a switching converter lose some of the dynamic performance but even make it unstable.

The basis for what can be considered a conventional design can be seen on [START_REF] Erickson | Fundamentals of Power Electronics 2nd Edition[END_REF]- [START_REF]Design Techniques for Preventing Input Filter Oscillations in Switched-Mode Regulators[END_REF]. This method allows to select the parameters of the filter by using simplifications of the Nyquist stability criteria. Although the method is straightforward to apply, the simplifications might lead to conservative results.

This paper takes another route by using Lyapunov-based methods. Such approaches allow to ensure stability and performance properties of a dynamical system by solving linear matrix inequalities (LMIs). Such numerical tools should help This work has been partially supported by project PID2020-120151RB-I00 financed by MCIN/ AEI /10.13039/501100011033. It has also been partially supported by University Rovira i Virgili through the Martí-Franquès Research grants Programme with reference 2019PMF-PIPF-26. to find regions in the space of the aforementioned parameters of the filter, ensuring design specifications. These specifications may include attenuation, stability, decay rate or damping ratio, to name some. Within the resulting region of parameters, the designer can choose appropriate components, according to other restrictions such as weight, cost, volume, fitting to commercial values, etc.

The paper follows the ideas in [START_REF] Ramos | An LMI approach to compute robust stability domains for uncertain linear systems[END_REF] to manipulate parameterdependent Lyapunov functions and derives some original LMI conditions that are useful in the definition of the constraints. Although the proposed methods could be used to design different parts of a power converter, these ideas are illustrated with two examples of input filter design for power conversion systems behaving as constant power loads.

The paper is organized as follows. First, the basis for the robust stability and performance is introduced in Section II. Section III describes a bisection-based method to enlarge the parametric space that satisfies the stability and performance constraints. Section IV shows the case example of input filter design and evaluates the robustness of the proposed approach. Finally, conclusions are given in Section V.

II. NUMERICAL METHODS FOR ROBUST STABILITY AND PERFORMANCE

The investigation of the domain of stability that we propose is based on a convex combination of Lyapunov functions, following the ideas in [START_REF] Ramos | An LMI approach to compute robust stability domains for uncertain linear systems[END_REF]. Besides of stability, transient performance is addressed by robust pole placement.

A. Robust Stability

Consider the following continuous-time linear system

ẋ(t) = Ax(t) (1)
It has been proven that A is asymptotically stable if and only if the LMIs

P > 0; A ′ P + PA < -I (2)
or equivalently, with

W = P -1 W > 0; AW + WA ′ < -I, (3)
admit a feasible solution. In ( 2) and (3), I is the identity matrix of appropriate dimensions and A ′ is the transpose of the state matrix A. The bold notation in P and W indicates that these are variables that must be found to prove stability. In the case where A is not precisely known, and as long as the unknown parameters appear linearly, it is possible to define A as a polytope that contains the uncertain domain of parameters. The polytope can be modeled as the convex combination of its N vertices A j as follows

A(p) = N j=1 λ j A j ; N j=1 λ j = 1; λ j ≥ 0. ( 4
)
where p represents the vector of parameters that are uncertain. If there exists a common Lyapunov matrix W that verifies (3) for all the vertices of the polytope (4), the stability of the uncertain matrix A can be assured. Note that only by verification of the vertices, the properties hold in the entire domain of p. Also, this is only a sufficient condition for robust stability and such approach often leads to conservative results.

For less conservative results, the set of LMIs

W j A ′ j + A j W j < -I; j = 1, ..., N W k A ′ j + A j W k + W j A ′ k + A k W j < 2 N -1 I j = 1, ..., N -1; k = j + 1, ..., N (5)
may be used. In [START_REF] Leite | Robust pole location by parameter dependent state feedback control[END_REF], W j corresponds to a convex combination of positive definite matrices. The resulting Lyapunov function is defined as in (4), such that W(p) depends on p.

W(p) = N j=1 λ j W j ; N j=1 λ j = 1; λ j ≥ 0. (6)
Proof of the sufficiency of these LMIs is given on [START_REF] Ramos | An LMI approach to compute robust stability domains for uncertain linear systems[END_REF]. The resulting parameter-dependent Lyapunov function provides some degree of flexibility, and therefore the results are less conservative than using a unique W for all vertices.

B. Robust Performance Constraints

The dynamical characteristics of the system to be designed can be specified by pole placement techniques. First, it is possible to relate some properties to a circle in the complex plane centered at the origin with radius ρ [START_REF] Leite | Robust pole location by parameter dependent state feedback control[END_REF]. Other characteristics can be related to a conic sector in the left-half plane with an angle equal to θ [START_REF] Chilali | H∞ Design with Pole Place Constraints: An LMI Approach[END_REF]. Figure 1 shows these conditions in the complex plane.

1) Circle of Radius ρ: The maximum natural frequency of oscillation of the poles of a system can be constrained to ρ. The LMI conditions that verify the location of the system poles in a generic circular region have been derived in [START_REF] Leite | Robust pole location by parameter dependent state feedback control[END_REF]. The LMIs for the particular case where the circle is centered at 0 + 0j are shown in Eqs. [START_REF] Erickson | Optimal single resistor damping of input filters[END_REF] to [START_REF] Mosek Aps | The MOSEK optimization toolbox for MATLAB[END_REF], where A rj = A j -ρI. 2) Conic Sector of Angle θ: As shown in [START_REF] Chilali | H∞ Design with Pole Place Constraints: An LMI Approach[END_REF], the dynamic system (1) has its poles inside a conic sector of angle θ if and only if there exists W > 0 such that

sinθ(AW + WA ′ ) cosθ(AW -WA ′ ) cosθ(WA ′ -AW) sinθ(AW + WA ′ ) < 0 (10)
combining this LMI with the sufficient condition for parameter-dependent Lyapunov functions described in [START_REF] Ramos | An LMI approach to compute robust stability domains for uncertain linear systems[END_REF], two novel LMIs (11) and ( 12) have been derived. These LMIs consider parameter-dependent Lyapunov functions [START_REF] Chilali | H∞ Design with Pole Place Constraints: An LMI Approach[END_REF] and are one of the contributions of the paper.

C. Parameter Region

As stated in the introduction, the proposed algorithm does not return a single parameter value, but rather a region or a domain of values. As in [START_REF] Ramos | An LMI approach to compute robust stability domains for uncertain linear systems[END_REF], [START_REF] Leite | Robust pole location by parameter dependent state feedback control[END_REF], consider the following continuous-time linear system:

ẋ(t) = A 0 + m i=1 α i E i x(t) (13)
where x ∈ R n , A 0 ∈ R n×n is the nominal matrix of the system, E i ∈ R n×n (i = 1, ..., m) are given matrices representing the perturbation directions and α i ∈ R (i = 1, ..., m) are scalar values defining the amount of perturbations allowed. The parameter dependent matrix in (13) is equivalent to the representation in [START_REF] Ramos | An LMI approach to compute robust stability domains for uncertain linear systems[END_REF]. Note that the sum of A 0 with the different E i matrices multiplied by the scalar values in the vector α i results in the vertices of the domain of A(p). It is worth to remark, however, that this approach allows to consider the influence of the parameters in p separately.

Each perturbation direction is related to a vector α i , which contains two values (lower and upper) that are connected to the range of possible parameter values. Thus, the vector is defined as

α i := [-α i , α i ].
The aim is to maximise the domains [-α i , α i ] ∀i = 1, . . . , m while ensuring stability and pole placement conditions.

A stable initial A 0 is assumed. In the case where pole placement is required, also an initial A 0 with its poles inside the region is assumed.

The problem is written such that each perturbation direction can be related to an uncertain parameter of the system. The range of values for each parameter depends on: the initial values (set on A 0 ), the scale of the parameters (set on E i ) and

A rj W j + W j A ′ rj A rj W j W j A ′ rj -ρW j ] < -I 0 0 0 ; j = 1, ..., N (7)
  A rj W j + W j A ′ rj + A rj W k +A rk W j + W j A ′ rk + W k A ′ rj A rj W j + A rj W k + A rk W j W j A ′ rj + W j A ′ rk + W k A ′ rj -ρ(2W j + W k )   < 1 (N -1) 2 I 0 0 0 j = 1, ..., N ; k = 1, ..., N ; k ̸ = j (8)         A rj W k + A rk W j + W j A ′ rk +W k A ′ rj + A rj W l + A rl W j +W j A ′ rl + W l A ′ rj + A rl W k +A rk W l + W l A ′ rk + W k A ′ rl A rk W j + A rj W k + A rl W j +A rj W l + A rl W k + A rk W l W j A ′ rk + W k A ′ rj + W j A ′ rl +W l A ′ rj + W k A ′ rl + W l A ′ rk -2ρ(W j + W k + W l )         < 6 (N -1) 2 I 0 0 0 j = 1, ..., N -2 ; k = j + 1, ..., N -1 ; l = k + 1, ..., N (9)
sinθ(A j W j + W j A ′ j ) cosθ(A j W j -W j A ′ j ) cosθ(-A j W j + W j A ′ j ) sinθ(A j W j + W j A ′ j ) < -I ; j = 1, ..., N (11)
sinθ(A j W k + A k W j + W j A ′ k + W k A ′ j ) cosθ(A j W k + A k W j -W j A ′ k -W k A ′ j ) cosθ(-A j W k -A k W j + W j A ′ k + W k A ′ j ) sinθ(A j W k + A k W j + W j A ′ k + W k A ′ j ) < 2 N -1 I j = 1, ..., N -1 ; k = j, ..., N (12)
the limits of α i , which are found by checking the feasibility of the LMI conditions. A problem particular to each system consists in finding the state matrix A and defining the perturbation matrices E in a way that it is possible to overcome the usually non linear relationships between the elements in the original matrix A and the parameters of the system. Perturbation matrices E i must be defined carefully, since the same parameter can appear in several positions of the matrix.

Considering the developments shown in this section, we derive the following proposition.

Proposition 1. For a given uncertain system (4), if there exist positive definite matrices W i , ∀i = 1, . . . , N , such that the conditions (7), ( 8), ( 9), ( 11) and (12) are satisfied, then the stability of matrix A(p) and its pole-placement performance in the subset formed by the circle of radius ρ and the conic sector of angle θ is guaranteed.

Proof 1. Directly follows from the elements given in Sections II-A and II-B.

III. BISECTION ALGORITHM

Proposition 1 allows to verify if a range of parameters in matrix A satisfy the conditions mentioned above. However the objective of the paper is to find this range of parameters in matrix A. Finding how large the perturbation amount can be is the job of the bisection algorithm.

Given an initial feasible A 0 matrix and perturbation directions E i , the bisection algorithm finds a feasible range of -α i , α i ∀i = 1, ..., m. An additional boolean vector ψ i is used in order to define the desired "direction of movement" for each α i . The vector is defined as ψ i := [ψ i , ψ i ], ∀i = 1, ..., m. As an example, for ψ i = [0, 1], the search is restricted to positive values of α i , and the corresponding vector

-α i , α i is [0, α i ].
The directionality can be useful to only increase or to only decrease some particular parameter through α i . As an example in the context of power electronics, the search for smaller parameters may translate to smaller components and therefore lower costs. Conversely, the search for larger parameters may be related to parasitic elements, in order to assess the influence of these undesired factors in the dynamic behaviour.

The search of the feasible domain is carried out in two steps. First, the maximum independent non-symmetric domain for each individual perturbation direction α i is found. This step is noted as independent because only one perturbation direction is tackled at a time. It is non-symmetric since the distances between A 0 and the limits of the stability domain are not necessarily symmetric i.e. maximum and minimum α i values -α i , α i do not necessarily have the same magnitude. Figure 2a shows the idea of these domains for case in which N = 2. The domain is found by bisection, i.e.: testing the feasibility or the unfeasibility of the LMI conditions.

The resulting feasible domains are used as a scaling parameter to obtain the simultaneous non-symmetric domain. Using the previous domain as a scaling parameter helps in dealing with the case when there is a big discrepancy in the possible range of movement between directions. Simultaneous means that all the E i are used simultaneously to create the polytope to be tested by the LMIs. Again, this domain is found by bisection. Now the bisection only works on a single scalar value of α and the individual α i are obtained as a product of α and the scaling parameter. Figure 2b shows an example of how the final polytope may look like. The algorithm described above can be posed as follows:

Algorithm 1.

1. Define the problem description: matrices A 0 , E i , ψ i and requirements ρ max , θ max .

2. Obtain the independent non-symmetric domains. ∀ i = 1, ..., m if ψ i = true, max(α i ) subject to LMIs (7), ( 8), ( 9), (11), (12). else α i = 0. if ψ i = true, max(α i ) subject to LMIs (7), ( 8), ( 9), (11), (12). else α i = 0. 3. Obtain the simultaneous non-symmetric domain.

max(α) subject to LMIs (7), ( 8), ( 9), (11), (12). 4. Each individual α i is obtained by α and the scaling parameter.

This final domain is used to obtain the range of feasible parameters of the system. Any combination of the α i obtained is feasible. This is also true for the parameters of the system if they are independent with respect to each other.

IV. CASE EXAMPLE: ROBUST FILTER DESIGN

The basis for what can be considered conventional filter design can be seen on [START_REF] Erickson | Fundamentals of Power Electronics 2nd Edition[END_REF]- [START_REF]Design Techniques for Preventing Input Filter Oscillations in Switched-Mode Regulators[END_REF]. These conventional methods

L b C f R f L f R f P i L f i L b v C f (b) LC filter with parallel R f and L b branch, loaded by a CPL.
Fig. 3: Generic and simplified models rely on simplifications of the stability criteria that might lead to conservative results. Figure 3a shows a typical representation of the design problem, a power converter with a control loop and an input filter. In Figure 3b, the converter is simplified as a constant power load (CPL), and the filter is a third order case. The tools described above are used to perform input filter design for two different filters.

Reducing the size of the filter is one of the objectives, which is in accordance with decreasing the capacitance and the inductance of the elements in the filter. A feasible oversized filter can be used as a starting point for the search of the feasible domain of parameters. As illustrated in Figure 4 the feasible domain depends on the initial system A 0 .

Therefore, a recursive approach that runs Algorithm 1 and finds new matrices A dnew is proposed. This iterative process will be finished when the range of the filter parameters (P i ) obtained is lower than some given threshold ∆ pi . A final bisection search will be performed, with the objective of finding the largest possible region of parameters satisfying the dynamical properties. This final search will start at the most favourable point found and go in the opposite direction in order to find a large parameter region. The algorithm used is defined on Algorithm 2. Algorithm 2. 1. Set the specifications of the problem. 2. Check A 0 feasibility with α i = 0. Otherwise an error is returned. 3. Run Algorithm 1 to obtain the α i . 4. Compute A dnew with an average value of α i . 5. If (P iold -P inew > ∆ i ) → A 0 = A dnew and go back to 3, else obtain A 0 with the maximum values of α i . 6. Run Algorithm 1 with ψ i = 1 -ψ i . 7. Compute the final parameter region A 0 and α i .

A. Case Example: LC Filter with parallel L-R damping branch

In this filter topology, shown in Figure 3b, the aim is to reduce the filter parameters while being able to obtain a parameter region from which a designer can extract viable filter parameters. The damping is given by a L-R branch in parallel with the filter inductance. The losses are low since the DC power does not involve the damping resistor. Such system has the following state space representation A matrix:

A =    0 0 -1 L f 0 -R f L b -1 L b 1 C f 1 C f β C f    , with x =   i L f i L b v C f   (14)
as the state vector, and β = P V 2 , P being the set power and V the voltage on the CPL. The following perturbation matrices will be defined:

E 1 =   0 0 -1 L fsc 0 0 0 0 0 0   ; E 2 =   0 0 0 0 0 0 1 C fsc 1 C fsc β C fsc   E 3 =   0 0 0 0 -R f 0 L bsc -1 L bsc 0 0 0   ( 15
)
This means that the uncertain model will have N = 2 3 = 8 vertices. The suffix " sc " is short for "scaling". All the parameters with this suffix are used in order to keep the perturbation values as close as possible to each other. The designed filter parameters can be obtained by

A des = A 0 + k i=1 α i E i ( 16
)
where A des and A 0 have the same structure as A. The initial parameters of the filter in A 0 are noted with the suffix " 0 ", whereas the range of parameters that has been found is noted with the suffix " d ". In the case at hand, each of the parameters in the filter only depends on one particular α i and E i (they are independent from each other), thus they can be found as follows:

L f d = (-A dj (1, 3)) -1 = 1 L f0 + α 1 1 L fsc -1 (17)
C f d = (A dj (3, 1)) -1 = 1 C f0 + α 2 1 C fsc -1 (18)
L b d = (-A dj (2, 3)) -1 = 1 L b0 + α 3 1 L bsc -1 (19)
For this particular filter, only one search is performed, implementing the three perturbation directions simultaneously. Note that, since α 1 , α 2 and α 3 are vectors of dimension 2, a range for each parameter L f , C f and L b is found. The damping resistance R d is obtained as the optimum value for damping the filter according to [START_REF] Erickson | Optimal single resistor damping of input filters[END_REF]. The search has been programmed such that after each run of the bisection algorithm a new R d is obtained which will be used to update the A 0 matrix. For the numerical example, performance conditions have been set to: attenuation of at least -0 dB at the switching frequency of f sw = 250 kHz and a damping ratio of at least ζ = 0.3.

The attenuation is enforced with the use of ρ, since the poles of the system are in agreement with the magnitude asymptotes of the input-to-output current transfer function. Limiting the frequency of the poles sets a minimum distance between them and the switching frequency f sw , at which a specific attenuation is desired. For a given attenuation atten, and a slope of the high frequency asymptote of the filter dB/Oct, ρ is found as follows:

ρ = 2πf sw 2 p ; p = atten dB/Oct (20)
Then, the minimum damping ratio is enforced through the conic region of angle θ. For a given minimum damping ratio ζ, θ is found as follows:

θ = cos -1 (ζ). (21)
Table I shows numerical values tested for this type of filter, including the initial values from which Algorithm 2 has started, the resulting range of filter parameters and finally a possible filter option. Note that the initial system is well oversized. The methods described above are able to reduce the overall size of the components and provide the designer with a range of parameters. From that range, the design engineer can choose the most suitable values according to parameters other than stability or dynamical performance. It can be seen that all of the parameters have a range that allows to easily find commercially available components. The huge upper range of L f indicates that in practice there is not limit on how large that inductor can be.

L f C f R d C b P R L C f v C f v C b i L f
Figure 5 shows the frequency response of the final filter (as well as the maximum and minimum values). Attenuation at f sw is higher than the specifications. It is also possible to check that the damping is higher than the specifications.

B. Case Example: LC Filter with RC Parallel Damping

The schematic for this type of filter can be seen on Figure 6. For this case the state matrix of the system is

A =    -R L L f -1 L f 0 1 C f β C f -1 R d C f 1 R d C f 0 1 R d C b -1 R d C b    , with x =   i L f (t) v C f (t) v C b (t)   (22)
as the state vector, and β = P V 2 . The perturbation matrices are defined as:

E 1 =   -R L L fsc -1 L fsc 0 0 0 0 0 0 0   , E 2 =   0 0 0 1 C fsc β C fsc - 1 R d C fsc 1 R d C fsc 0 0 0   , E 3 =   0 0 0 0 -1 R dsc C f 1 R dsc C f 0 0 0   , E 4 =   0 0 0 0 0 0 0 1 R dsc C bsc -1 R dsc C bsc   . (23)
In order to show the versatility of the method, we employ a different approach in this case. Algorithm 2 is employed twice, for two different sets of parameters of the filter. The first time, the algorithm finds a range for the values of L f and C f (E 1 and E 2 ). The second time, the algorithm looks for the maximum admissible range of values of R d and C b (E 3 and E 4 ). Since only two perturbation matrices are applied at the same time the resulting uncertain model will have N = 2 2 = 4 vertices. The range of parameters is obtained by: In practice, such a large upper bound means that there is no limit in choosing a large capacitor Cb. Notice that since the parameters R d and C b depend on each other (27) even though the alpha region is an orthogonal plolytope the filter parameter region is not. In order to show a range of parameters that are all compatible with each other, the region previously mentioned has been conservatively constrained to be orthogonal. It could be the case that this simplification limits the values too much, in that case it would be necessary to look at the non-orthogonal region.

L f d = -A d (1, 2) -1 = - 1 L 0 + α 1 1 L sc -1 (24)
C f d = A d (2, 1) -1 = 1 C 0 + α 2 1 C sc -1 (25)
R d d = (C f A d (2, 3)) -1 = 1 R d0 1 + α 3 1 R dsc -1 (26)
C b d = (R d d A d (3, 2)) -1 = R d d R d0 C b0 + α 4 R d d R dsc C bsc
Figure 7 shows the magnitude frequency response of the possible and chosen filter values. It can be seen that the

C. Influence of Initial Conditions

One of the possible problems of the proposed approach is its sensitivity to different initial conditions in A 0 . Although there is no guarantee that the proposed method will achieve good results in all possible cases, the results have been found to present reasonable robustness to changes in A 0 . In order to illustrate this point, a few tests are shown. The tests compare the results in the previous section with results obtained with different A 0 .

The example from Section IV-A has been recalculated with smaller and larger initial parameters. The results are given on Table IV. First, it can be seen how the lower bound for L f , C f and L b is relatively similar. Also, for different initial conditions the final ranges of values present differences, but they still allow the designer to choose similar or even smaller components.

The case in Section IV-B has been tested against initial conditions with values 40 % higher or lower than those in the original A 0 . The resulting ranges can be seen in Table V. The resulting region still allows the designer to choose the same final values. In addition, it can be seen how the case that allows the lowest bound in R d is also the case that exhibits a higher lower bound in C b . This illustrates the trade-off between the different ranges in the components.

D. Numerical Complexity

The method proposed in the paper employs LMIs that are efficiently solved by interior-point methods. The complexity of the problem can be computed in terms of the number of parameters to be evaluated m and the size of the Lyapunov matrices, which is linked to the number of states in A, i.e. n. Additionally, since the algorithms rely on a feasibility test, the precision tolerance of the bisections also plays a role.

In the examples demonstrated above, conventional computers at the time of write-up of the paper, using MATLAB with YALMIP [START_REF] Löfberg | YALMIP[END_REF] and MOSEK [START_REF] Mosek Aps | The MOSEK optimization toolbox for MATLAB[END_REF] can solve Algorithm 1 with running times in the order of seconds up to one or two tens of seconds. Then, Algorithm 2 calls Algorithm 1 around ten V. CONCLUSIONS This paper presented a numerical approach to physical parameter selection by means of uncertain models. LMIs have been used to implement performance conditions which force a minimum attenuation and damping of the system. A modified bisection algorithm has been used in order to find a region of parameters compatible with the performance conditions. It has been shown that it is possible to start with large parameters and iterate to find parameter regions with smaller overall parameters until certain tolerance conditions are met.

It is important to keep in mind that this approach may lead to conservative results and thus it does not guarantee the most optimised values. Nevertheless it gives the designer a tool to quickly estimate a feasible set of parameters.

A couple of design examples have been shown. Two different filter topologies have been used in order to implement the numerical methods described. Any filter value defined inside the found parameter region can be chosen according to other design criteria such as cost, size or other arbitrary criteria.

Fig. 1 :

 1 Fig. 1: Region defined with the pole placement LMIs.

  Maximum independent non-symmetric domain example.

  Simultaneous non-symmetric domain example.

Fig. 2 :

 2 Fig. 2: Bisection algorithm domains examples

  Generic filter plus closed loop converter.

Fig. 4 :

 4 Fig. 4: Example of the domain obtained by the bisection algorithm with A 0 and the domain obtained with A dnew .

Fig. 5 :

 5 Fig. 5: Input current -to -output current io iin frequency response of the LCRL filter. The solid line is for the chosen values and the dashed lines represent the extreme vertices of the polytope, i.e. the frequency response with the maximum and the minimum values of the parameters.

Fig. 6 :

 6 Fig. 6: LC filter with parallel RC branch.

  and III show the results of the example design. The performance conditions are an attenuation of -80 dB at the frequency of 250 kHz and a damping ratio of ζ = 0.3. As it can be seen in the table, the initial values of the filter are oversized to L = 470 µH, C = 800 µF, R d = 1 Ω, and C b = 2000 µF. The first search returns a range of values from which the designer can choose: L ∈ 73.2 -1890 µH and C ∈ 59.3 -162 µF. Once these components are chosen to be L = 220 µH and C = 100 µF, the second search returns another range of values for R d and C b : R d ∈ 0.67 -1.5 Ω and C b ∈ 458 -410e9 µF.

Fig. 7 :

 7 Fig. 7: Input current -to -output current io iin frequency response of the filter.

TABLE I :

 I LC with parallel RL branch numerical example

TABLE II :

 II Filter values after L f and C f iteration

	Param.	First Selection	Result. Range	Final Selection
	L f /µH	220	220	220
	C f /µF	100	100	100
	R d /Ω	1	0.67 -1.5	0.7
	C b /µF	2000	458 -410e9	940(2x470)

TABLE III :

 III Filter values after R d and C b iteration attenuation at f sw = 250 kHz is better than the specification of -80 dB set by the design constraints. The damping constraints are also met.

TABLE IV :

 IV Comparison of different initial points for the LC filter with parallel RL branch times per design cycle which gives a rough estimate of how long the whole procedure can take.

TABLE V :

 V Filter values after L f and C f iteration, with double initial values