Judy Najnudel

Thomas Hélie

David Roze

Henri Boutin

Simulation of an ondes Martenot circuit

Keywords: Electronic music instrument, non-linear, powerbalanced, simulation

come

I. INTRODUCTION

A S the audio industry is evolving from analog to digital, the preservation of analog machines and instrumentsand with them, that of the craft of generations of engineersis proving critical [START_REF] Davies | The preservation of electronic musical instruments[END_REF], [START_REF] Bonardi | The preservation, emulation, migration, and virtualization of live electronics for performing arts: An overview of musical and technical issues[END_REF]. The concern has been especially raised about the ondes Martenot; invented in 1928 by Maurice Martenot, this pioneer electronic musical instrument [START_REF] Leipp | Les ondes Martenot[END_REF] is no longer produced in its original form and some of its components are now obsolete, such as specific triode vacuum tubes. A solution consists of modeling its circuit and building a real-time virtual instrument, so that the community of composers, musicians and musicologists may have access to virtual facsimiles.

Several techniques are available to build virtual analog instruments (for a review, see [START_REF] Pakarinen | A review of digital techniques for modeling vacuum-tube guitar amplifiers[END_REF]- [START_REF] Välimäki | Virtual analog effects[END_REF]). They include Modified Nodal Analysis [START_REF] Ho | The modified nodal approach to network analysis[END_REF], [START_REF] Wedepohl | Modified nodal analysis: An essential addition to electrical circuit theory and analysis[END_REF], the Nodal DK method [START_REF] Holters | Physical modelling of a wah-wah effect pedal as a case study for application of the nodal dk method to circuits with variable parts[END_REF] (for audio applications, see e.g. [START_REF] Yeh | Automated physical modeling of nonlinear audio circuits for real-time audio effects-part i: Theoretical development[END_REF], [START_REF] Yeh | Automated physical modeling of nonlinear audio circuits for real-time audio effects-part ii: BJT and vacuum tube examples[END_REF]), Wave Digital Filters [START_REF] Fettweis | Wave digital filters: Theory and practice[END_REF] (see e.g. [START_REF] Bilbao | Wave and scattering methods for numerical simulation[END_REF]- [START_REF] Bogason | Modeling time-varying reactances using wave digital filters[END_REF]). Real-time simulations of circuits with vacuum tubes have been derived from such methods [START_REF] Angelo | New family of wave-digital triode models[END_REF], [START_REF] Pakarinen | Wave digital modeling of the output chain of a vacuum-tube amplifier[END_REF] or from time-continuous state-space representations combined with numerical schemes [START_REF] Dempwolf | Discretization of parametric analog circuits for real-time simulations[END_REF]- [START_REF] Macak | Simulation of Fender type guitar preamp using approximation and state space model[END_REF]. The approach used in this paper is based on a state-space representation, which satisfies the power balance of the physical system structured into conservative, dissipative and external parts, known as Port-Hamiltonian Systems (referred to as PHS) [START_REF] Maschke | An intrinsic Hamiltonian formulation of network dynamics: Non-standard Poisson structures and gyrators[END_REF], [START_REF] Van Der Schaft | Port-Hamiltonian systems theory: An introductory overview[END_REF]. This formulation can be combined with numerical methods that preserve the power balance structure (and passivity) in the discrete-time domain for both linear and non-linear systems. This has proved to be relevant for simulations of audio (electronic or multi-physical) systems [START_REF] Lopes | Explicit second-order accurate method for the passive guaranteed simulation of port-Hamiltonian systems[END_REF]- [START_REF]Passive simulation of the nonlinear port-Hamiltonian modeling of a Rhodes piano[END_REF]. In a previous paper [START_REF] Najnudel | Simulation of the ondes Martenot ribbon-controlled oscillator using energy-balanced modeling of nonlinear time-varying electronic components[END_REF], the first stage of the ondes Martenot circuit (a variable oscillator) has been modeled and simulated using this method, with a special focus on electro-mechanical power balance. This paper addresses the full modeling and simulation of the circuit of the ondes Martenot No. 169, which is composed of 5 coupled stages. Moreover, analyzing numerical experiments allows some model simplifications, which in turn grant access to real-time. This paper is organized as follows. Section II describes the history and the composition of the ondes Martenot as well as the way the instrument is played. In section III, the PHS formalism is briefly recalled and illustrated with a simple example. In section IV, each isolated stage is first separately modeled as a PHS and the coupling of the stages is then addressed. In section V, the physical parameters for simulation are examined and set. Section VI is devoted to off-line powerbalanced simulations of the complete circuit. Finally, in section VII, simplifications are proposed based on the analysis of numerical experiments. They lead to a reduced order PHS from which a real-time audio plug-in is built for the sound synthesis.

II. HISTORY, ORGANOLOGY, CIRCUIT AND PRINCIPLE

A. History

The advent of electrical engineering marked a disruption in many industrial domains, amongst them musical instrument making. Providing both new timbres and new playing modalities, it gave the 20th century instrument corpus its singularity. Nonetheless, instrument making is not immune to technological obsolescence, and only few of these innovative instruments are still played today. This is the case of the ondes Martenot; only 264 ondes were built, but around 1500 pieces were composed for it [START_REF] Laurendeau | luthier de l'électronique[END_REF]. As they were handmade and subjected to successive enhancements, each ondes is unique and identified by a reference number. The considerable task of cataloguing and evaluating the condition of the remaining instruments has been achieved by the Musée de la Musique in Paris [START_REF] Ramel | Conservation and restoration of electroacoustic musical instruments at the Musée de la Musique, Paris[END_REF]. It appears that for some of the ondes, making them playable again would necessitate irreversible interventions. This poses a big deontological dilemma from the preservation point of view and argues in favor of building facsimiles.

B. Organology

This paper focuses on the ondes No. 169, manufactured in 1937, and kept in the Musée de la Musique in Paris. This specimen has a great heritage value, particularly because of its electronic circuit with vacuum triodes. Also, it is kept in playing conditions and is fully documented [START_REF] Courrier | Analyse de fonctionnement onde 169[END_REF]. It is based on heterodyne processing, like the Theremin [START_REF] Glinsky | Theremin: Ether music and espionage[END_REF] and the Trautonium [START_REF] Couprie | Oskar sala[END_REF]. It generates a harmonic signal which feeds specific loudspeakers (called diffuseurs, Fig. 1a). The player controls the pitch with a sliding ribbon (continuous pitch change) or-in later versions of the instrument-with a mobile keyboard (discrete pitch change and vibratos, Fig. 1b). The loudness is controlled with an intensity key, by pressing a small bag containing a mixture of conductive and insulating 1a on the right, from [START_REF] Laurendeau | luthier de l'électronique[END_REF]), circuit, ribbon, keyboard (Fig. 1b), intensity key with the powder bag underneath circled in red (Fig. 1c).

powders. This key operates as a rheostat [START_REF] Quartier | Intensity Key of the Ondes Martenot: An Early Mechanical Haptic Device[END_REF] (Fig. 1c). Thanks to this intensity key, the instrument's expressiveness has been praised since its creation. However, the lack of robustness of its circuit (subject to drifting, particularly in older models) is dreaded by musicians (called ondists) who are accustomed to perform repairs on the fly before concerts.

C. Circuit and operation

The technique of heterodyning was developed originally for wireless telegraphy [START_REF] Belrose | Reginald Aubrey Fessenden and the birth of wireless telephony[END_REF], in order to shift high frequency signals into the audio domain. In the ondes Martenot specifically, two oscillators generate quasi-sinusoidal voltages at high frequency; one is set to 80 kHz and the other, variable below 80 kHz, is controlled by the player. Their sum is an amplitude-modulated signal, whose envelope is detected with a triode vacuum tube in series with an RC circuit. The envelope fundamental frequency, equal to the difference between the oscillators' frequencies, lies in the audible frequency range. The resulting voltage is amplified through two successive triode vacuum tubes, which increase the harmonic distortion present in the signal, due to their non-linear characteristics.

Finally, the diffuseur, selected by the player, converts the electrical waveform into sound and in turn modifies its spectral content. The complete circuit is structured into five stages devoted to specific functionalities (see Fig. 2a), namely : (1) fixedfrequency oscillator, (2) variable-frequency oscillator, (3) demodulator, (4) preamplifier, and (5) power amplifier. These stages are connected through transformers that introduce coupling: this will be naturally taken into account in the global PHS model (addressed in section IV-C). However, in first step, these stages are presented separately, ignoring coupling (no load on the secondary winding, 2b-2d). Qualitative specificities of these stages are stated below.

1) Oscillators: In both oscillators, a voltage noise V start at the triode grid is amplified by the triode vacuum tube and filtered by a LC circuit, before feeding the grid through a transformer of ratio P . Each transformer has two windings on the secondary; the second secondaries (for oscillators 1 and 2) are connected in series to form the demodulator voltage input.

2) Demodulator: The voltage input, equal to the sum of the oscillators' outputs, can be assimilated to an amplitudemodulated sinewave: cos(Φ) + cos(Φ -φ m ) = 2 cos(Φφ m /2)cos(φ m /2), where F = Φ/2π = 80 kHz is the fixed oscillator frequency, f m = φm /2π is the target frequency, and F -f m is the variable oscillator frequency. The demodulator is composed of a triode vacuum tube in series with an RC circuit (R 4 and C 21 on Fig. 2c). As the triode grid polarization is close to zero, the qualitative behavior of the triode between the grid and the cathode is that of a diode: grid current only appears for positive grid voltages (Fig. 3). The triode vacuum tube is loaded with capacitors in parallel with a transformer. This load forms a resonant bandpass RLC filter (Fig. 2c): the resistance R p models the transformer losses, L p is the primary inductance, C dem is the equivalent capacitor of C 9 and C 10 in parallel. Its resonance frequency lies in the audible range. The stage amplification in a common cathode pattern increases with the triode load. As a consequence, the oscillators' frequencies and their intermodulation products are filtered out before amplification.

3) Preamplifier and power amplifier: Both stages rely on the triode vacuum tube amplification properties in a common cathode pattern [START_REF] Cohen | Modélisation, analyse et identification de circuits non linéaires: Application aux amplificateurs guitare à lampes pour la simulation en temps réel[END_REF].

III. POWER-BALANCED MODELING

A. Framework

Detailed presentations of PHS are available in [START_REF] Van Der Schaft | Port-Hamiltonian systems theory: An introductory overview[END_REF]. Here we rely on a differential-algebraic formulation adapted to multiphysical systems [START_REF] Falaize | Passive guaranteed simulation of analog audio circuits: A port-Hamiltonian approach[END_REF], which allows the representation of a dynamical system as a network of 1) components storing energy E as a regular definite nonnegative function of their state x (E = H(x) ≥ 0), 2) instantaneous dissipative components for which efforts e w are related to flows w through a constitutive law z that dissipates power (e w = z(w), P diss = z(w) w ≥ 0, ∀w), 3) connection ports as control inputs u and their associated outputs y such as u y is the outgoing power P ext .

(a) Schematics of the complete circuit of the ondes No. 169 (as is from [START_REF] Courrier | Analyse de fonctionnement onde 169[END_REF]).

C k R k Cosc 6F5 1 : P Losc V bias V start 6F5 V out (b) C k R k C dem Rp V bias Lp 6C5 R 4 C 21 V out V in (c) C k R k Rp V bias Lp 6C5 V out V in (d) C k R k Rpw V bias 2A3 V out V in (e)
C dem Rp -V in -V out V in V out
Fig. 3: Equivalent small signal schematic of the demodulator. The nodes c, g and p denote the triode cathode, grid and plate respectively; r a is the triode plate resistor.

The variables are generally time-dependent. For each component, ẋ accounts for the incoming state flow. It is named "flow" in the following. 1 Jointly, ∇H(x) is named "effort" (w.r.t. the state) in the following, so that the power received by the energy-storing component ẋ∇H(x) is the product of an effort and a flow. If such a system is realizable [START_REF] Brockett | Finite dimensional linear systems[END_REF], [START_REF] Van Der Schaft | A realization procedure for systems of nonlinear higher-order differential equations[END_REF], the exchanges of generalized flows and efforts between the system components are coupled through a skew-symmetric interconnection matrix S = -S :

   ẋ w y    F (flows) = S    ∇H(x) z(w) u    E(efforts)
.

(

) 1
For convenience, to distinguish storage, dissipative, and connection ports exchanges, we introduce the following blockmatrix notation:

S = J x -K -G x K J w -G w G x G w J y . ( 2
)
The skew-symmetry of S guarantees the system power balance. Indeed,

∇H(x) ẋ dE dt + z(w) w P diss ≥0 + u y Pext = E F, (3)
and

E F = E SE= (E SE) = -E SE = 0, ( 4
) V I R i R L i L C i C
Fig. 4: Current-controlled parallel RLC circuit.

C q q = i C q 2 2C q C = v C L φ φ = v L φ 2 2L φ L = i L w z(w) R v R v R R = i R
so we have the following power balance:

dE dt = -P ext -P diss . (5)
Note that -P ext is the incoming power (since P ext is the outgoing power). The passivity of the system stems from this power balance and the positivity of P diss .

B. Example

Consider a linear parallel RLC circuit (Fig. 4). The capacitor C and the inductor L are storage components whose states are given by the variables q (charge) and φ (magnetic flux) respectively; the resistor R is a dissipative component described by Ohm's law. In order to represent it as a PHS, the system is chosen to be current controlled. Table I recaps the variables and the constitutive laws of the three components. Based on Kirchoff's laws, matrix S relating efforts (E) to flows (F) in Eq. ( 1) is found to be

    iC vL vR V =       0 -1 -1 -1 1 0 0 0 1 0 0 0 1 0 0 0     vC iL iR I .
Flows which are currents (resp. voltages) have associated efforts which are voltages (resp. currents). Therefore their products are powers and Eq. ( 3) is naturally retrieved. Note that in this simple example, the matrix S is sparse with constant coefficients, but the properties of the PHS formalism also hold for non-linear or coupled systems (state depending interconnection matrices, non-quadratic H, non-linear z).

IV. ONDES MARTENOT PHS MODELING

A. Ondes Martenot components

The ondes Martenot circuit contains resistors, capacitors, inductors, transformers, and triode vacuum tubes. Among those components, the ribbon-controlled capacitor and the triode need to be carefully modeled to satisfy a physical power balance and to admit a PHS formulation. 1) Ribbon-controlled capacitor: The variable capacitor C 15 (Fig. 2a) is controlled by a partly conductive ribbon: as the player moves the ribbon along, it hides or "activates" the teeth of the comb-shaped electrical conductor, located at different distances from the ribbon (Fig. 5). Previous work [START_REF] Najnudel | Simulation of the ondes Martenot ribbon-controlled oscillator using energy-balanced modeling of nonlinear time-varying electronic components[END_REF] showed that the LTI capacitors C 13 , C 14 but also the variable capacitor C 15 in parallel (Fig. 2a, stage 1) are equivalent to a single variable capacitor whose energy function is given by

H cap (q, d) = q 2 2C osc (d) , (6)
C osc (d) = 1 4π 2 F -A 1 2 d 12d 0 2 L osc , (7)
where d is the ribbon displacement relative to the position corresponding to the lowest note A 1 , d 0 is the ribbon displacement for one semitone, F is the fixed oscillator frequency and L osc is the inductance of the oscillator.

2) Triode vacuum tubes: Numerous models are available for triode vacuum tubes [START_REF] Leach | Spice models for vacuum-tube amplifiers[END_REF], [START_REF] Cardarilli | Improved large-signal model for vacuum triodes[END_REF]. For audio applications, the Norman Koren model [START_REF] Koren | Improved vacuum tube models for Spice simulations[END_REF] has proven to be accurate. Here we will rely on an enhanced Norman Koren model [START_REF] Cohen | Measures and parameter estimation of triodes, for the real-time simulation of a multi-stage guitar preamplifier[END_REF] which takes into account the triode grid current. Although this current is small compared to that of the plate, it must be considered to ensure the passivity of the component (since its power supply is external, the triode is indeed passive). As a result, the triode vacuum tube is a dissipative component with w = (v pc , v gc ) and z θ (w) = (i pc , i gc ) , where p, g and c respectively refer to the triode plate, grid and cathode. The plate current i pc and the grid current i gc are given by

i pc = 2[E 1 v pc , v gc ] Ex /K g if E 1 ≥ 0 0 else, (8)
i gc =      0 if v gc < V a v gc -V a R gk else, (9)
with

E 1 v pc , v gc = v pc K p ln      1 + exp     K p    1 µ + v gc + V ct K vb + v 2 pc             .
For each type of vacuum tubes used in the circuit of ondes Martenot No. 169, the set of constant parameters Fig. 6: Dissipated power for the triode 6C5. ] through a least squares minimization. The results of that minimization fulfill the passivity criterion, i.e. z θ (w) w ≥ 0 for all w (Fig. 6). They are given in Table II, section VI.

θ = (µ, E x , K g , K p , K vb , V ct , V a , R
3) Other components: Constitutive laws of R, L, C components are supposed to be linear, as shown in Table I. In the oscillators, the transformers are also supposed to be linear; their ratios P being unknown a priori, they are subjected to an estimation described in section V-A2.

B. PHS of isolated stages

For each stage, we denote R k the regulation resistor between the triode cathode and the ground, C k the bypass capacitor in parallel with R k , V bias the triode DC power supply, V in the input voltage and V out the observed output voltage. We denote L p the transformer magnetizing reactance, R p the transformer resistance modeling the core losses in the demodulator and the preamplifier and R pw the input impedance of the diffuseur. For each oscillator, we denote C osc the capacitor equivalent to the parallel capacitors and L osc the inductance equivalent to the coils in series (Figs. 2b-2e). Kirchhoff's laws yield the interconnection matrices for the oscillator (Fig. 7a), demodulator (Fig. 7b), preamplifier (Fig. 7c), and power amplifier (Fig. 7d).

C. Stages coupling

The connection of two PHS is still a PHS. The total state is the concatenation of the subsystems states, and the total Hamiltonian is the sum of the subsystems Hamiltonians. A conservative interconnection is then achieved by connecting each output port of one system to exactly one input port of the other. For two stages A and B connected through a transformer of ratio ρ, we have V inB = ρV outA and I outA = -ρI inB , that is, y outA

y inB = 0 1/ρ -1/ρ 0 u outA u inB .
The circuit A connected to the circuit B is then equivalent to a circuit C whose PHS is given by Eq. ( 1), with notations of Eq. ( 2): x C = (x A , x B ) and ∇H C , w C , z C are defined likewise,

J xC = J xA 0 0 J xB
and J wC , G xC , G wC , K C are defined likewise, u C = (u inA , u outA , u inB , u outB ) and y C is defined likewise, and

J yC =       J yA 0 0 1/ρ 0 0 -1/ρ
J yB 0 0 .

V. ESTIMATION OF THE INSTRUMENT PARAMETERS

Most parameters of the components are provided by the circuit specifications [START_REF] Courrier | Analyse de fonctionnement onde 169[END_REF]. The others have to be determined, in particular the oscillator transformer re-injection ratio P and bypass capacitor C k , as well as the characteristics of the transformers towards the amplification stages.

A. Oscillators (stages 1 and 2) 1) Voltage noise source: The oscillation arises from the voltage noise V start at the grid of the triode at its operating point. This broadband noise mainly results from thermal agitation [START_REF] Llewellyn | A study of noise in vacuum tubes and attached circuits[END_REF] and is modeled with a white noise source of peak voltage 1 mV2 .

2) Re-injection transformer ratio: The voltage noise is amplified by the triode, filtered by a LC circuit at the plate and sent back to the grid through a transformer of ratio P (Fig. 2b). To estimate a range of possible values for P , the oscillator circuit is simplified: L osc acts as the transformer primary and the cathode potential is considered constant. Moreover, the oscillator model is linearized around its operating point and it is assumed that the influence of the load on the second secondary is negligible. The oscillator can then be represented as a system of input V start and output V out . Under those assumptions, it can be shown [START_REF] Najnudel | Simulation of the ondes Martenot ribbon-controlled oscillator using energy-balanced modeling of nonlinear time-varying electronic components[END_REF], [START_REF] Rocard | Dynamique générale des vibrations[END_REF] that the system characteristic equation in Laplace domain is s 2 + s(1 -µP )/(r a C osc ) + ω 2 0 (µ being the triode amplification factor, r a the triode plate resistor and ω 2 0 = 1/(L osc C osc ) the oscillation pulsation), and that the condition on P for the system to oscillate is

1 µ ≤ P < 1 µ + 2ω 0 C osc r a µ . ( 10
)
The closer to 1/µ P is chosen, the more stationary the oscillation amplitude will be. Fig. 8 shows the range of P values for which an oscillation is possible.

3) Inductance of the LC circuit: For a parallel LC circuit of resonant frequency F , L osc = 1/ C osc (2πF ) 2 . With F set to 80 kHz and C osc = 544 pF, this yields L osc = 7 mH.

4) Cathode bypass capacitor: The bypass capacitor C k stabilizes the cathode potential, which otherwise would depend strongly on the plate current (through the regulation resistor R k ), itself depending on the input voltage. This stabilization ensures that the signal to be amplified v gc is the biased input voltage and not a modulation of it. In parallel with (a) (b)

∇H(x) z(w) u i Losc v Cosc v C k igc ipc i R k Vstart Iout V bias                         v Losc . 1 . . . . . . . ẋ i Cosc -1 . . -P 1 . . -1 . i C k . . . 1
∇H(x) z(w) u i Lp v C 21 v C k v C dem i R 4 i Rp igc ipc i R k Vin Iout V bias                                     v Lp . . . 1
∇H(x) z(w) u i Lp v C k v Rp igc ipc i R k Vin Iout V bias                         ẋ v Lp 1 . . . i C k . 1 1 -1 i Rp -1 . . . 1 . . -1 . vgc . -1 . . . . 1 . . w vpc . -1 -1 . . . . . 1 v R k . 1 . . . . . . . Iin . -1 . . y Vout 1 . . . I bias . . - 1 .
(c)

∇H(x) z(w) u v C k v Rpw igc ipc i R k Vin Iout V bias                     ẋ i C k . 1 1 -1 i Rpw . . . 1 . . -1 . vgc -1 . . . . 1 . . w vpc -1 -1 . . . . . 1 v R k 1 . . . . . . . Iin . -1 . . y Vout 1 . . . I bias . . -1 .
(d) Fig. 7: PHS interconnection matrices of the isolated stages: oscillators (Fig. 7a), demodulator (Fig. 7b), preamplifier (Fig. 7c), power amplifier (Fig. 7d). R k , it constitutes a current-current high-pass filter whose frequency is f c = 1/(2πR k C k ). For the oscillator amplification to be maximal at the frequency F , f c must verify F ≥ 10f c . For the given value R k = 7.5 kΩ, this yields C k ≥ 3 nF. In the following, C k is set to 0.22 µF, a typical capacitor value, for both oscillators.

B. Coupling transformer characteristics

Except for its ratio, the properties of the transformer between stage 3 (resp. 4) and 4 (resp. 5) are unknown and accessing a real circuit to perform measurements has proven to be delicate due to the instrument rarity and fragility. It is assumed that Maurice Martenot chose the transformers with audio applications in mind; the magnetizing reactance L p and the resistance modeling the core losses R p will be thus set as such in the following. The capacitor C dem in parallel with R p and L p (Fig. 2c) forms a current-current band-pass filter, whose high cutoff frequency is given by f ch = (r a + R p )/(2πr a R p C dem ), with C dem = 7 nF and r a = 10 kΩ. The instrument ambitus is C1-B7 corresponding to a fundamental frequency ranging from 32 Hz to 3951 Hz in A440 tuning. Therefore, choosing f ch = 7.9 kHz (the highest possible note second harmonic) yields R p ≈ 4 kΩ. Similarly, the low cutoff frequency is given by f cl = (r a R p )/(2π(r a + R p )L p ). Choosing f cl = 50 Hz yields L p ≈ 9 H. It is assumed that the transformer to the power amplifier has the same properties. The values of the other circuit parameters are provided by the condition report of the ondes No. 169 [START_REF] Courrier | Analyse de fonctionnement onde 169[END_REF].

VI. POWER-BALANCED SIMULATION This section is devoted to the power-balanced simulation of the ondes Martenot No. 169 complete circuit, composed of the 5 coupled stages (the power amplifier being loaded by a basic resistor equivalent to that of a diffuseur-around 1.5 kΩ).

A. Discrete-time system

The simulation is performed using the PyPHS library [START_REF] Falaize | Passive guaranteed simulation of analog audio circuits: A port-Hamiltonian approach[END_REF], [START_REF] Falaize | PyPHS: Passive modeling and simulation in python[END_REF]. It is based on the numerical method recalled below in the case of mono-variate storage components (H(x) = N n=1 H n (x n ) where N is the storage components number).

Define x[k] = x(k/F s ) where F s = 1/T s denotes the sampling rate. Moreover, define the state increment δx[k] = x[k + 1] -x[k] and the discrete gradient (see [START_REF] Falaize | Passive guaranteed simulation of analog audio circuits: A port-Hamiltonian approach[END_REF], [START_REF] Itoh | Hamiltonian-conserving discrete canonical equations based on variational difference quotients[END_REF])

∇H x[k], δx[k] , the n th component of which is [ ∇H x[k], δx[k] ] n =        H n x n [k] + δx n [k] -H n (x n [k]) δx n [k] if δx n [k] = 0 dHn dx n (x n [k]) otherwise. ( 11
)
The method consists of replacing ẋ with δx[k]/T s and ∇H(x) with ∇H x[k], δx[k] in Eq. [START_REF] Davies | The preservation of electronic musical instruments[END_REF]. As in Eqs. [START_REF] Leipp | Les ondes Martenot[END_REF][START_REF] Pakarinen | A review of digital techniques for modeling vacuum-tube guitar amplifiers[END_REF][START_REF] Sanctis | Virtual analog modeling in the wavedigital domain[END_REF], due to the skew-symmetry of S, the following discrete power balance is satisfied:

∇H x[k], δx[k] δx[k] T s δE[k]/Ts + z(w[k]) w[k] P diss [k] + u[k] y[k] Pext[k] = 0.
Moreover, due to the positivity of z(w[k]) w[k], the system is also still passive. Denoting ν = δx[k]F s , w[k] , the discretization of Eq. ( 1) yields the implicit equation

F (ν) := ν -S Ē(ν, x[k], u[k]) = 0, (12)
where

S = J x -K -G x K J w -G w and Ē(ν, x[k], u[k]) = ∇H x[k], δx[k] , z(w[k]), u[k] .
Denoting F the Jacobian of F , Eq. ( 12) is solved using the Newton-Raphson iteration

ν k+1 = ν k -F (ν k ) -1 F (ν k ), (13)
if F is invertible. Conditions for convergence of the method can be found in [START_REF] Muller | Power-Balanced Modelling Of Circuits As Skew Gradient Systems[END_REF], [START_REF] Deuflhard | Newton methods for nonlinear problems: affine invariance and adaptive algorithms[END_REF].

B. Numerical experiments

The circuit is simulated based on the concatenation of the stages PHS with component parameters values of Table II. The control law of the ribbon position d is built to generate a frequency sweep according to parameters in Table III and Eq. [START_REF] Ho | The modified nodal approach to network analysis[END_REF]. Fig. 9a-9c show the output voltages of stages 3, 4 and 5 during the simulation, and the effect of the successive harmonic distortions. Fig. 9d shows the spectrogram of the complete circuit output during the whole sweep. Fig. 9e shows the stored power, outgoing power and dissipated power of the complete circuit during the whole sweep. As shown in Fig. 9f, the total power, sum of these three contributions, amounts to less than 10 -13 W.

VII. AUDIO PLUGIN A. Circuit simplification

The voltage generated by the oscillators, even when connected to the other stages, is nearly sinusoidal: the harmonic distortion is about 0.03% for harmonic 2 (Fig. 10a-10b). For a real-time application, physically modeled oscillators may be thus replaced with a sinewave generator.

Additionally, although the simulated performances of the power amplifier are rather poor for a class A amplifier (almost 5% harmonic distortion for harmonic 2, Fig. 11), the stage contributions to the final sound are of much less importance compared to the demodulator and preamplifier. In practice, based on numerical experiments, this stage may also be removed for real-time application.

B. Oscillators frequency reduction

The original oscillator frequency F is around 80 kHz. This choice, alongside with heterodyning, can certainly be attributed to the cost and unreliability of high value capacitors in Maurice Martenot's time: an oscillator at an audible frequency, 3000 Hz for instance, would necessitate a capacitor of 4 µF (the inductance being equal for both oscillators). But in a virtual instrument context, the components may take

C. GUI, real-time results and performances

The dimension of the complete PHS is defined as n =dim E =dim F. For each sample and each element, the matrix update computes n multiplications and n -1 additions, and consequently F s ×n(2n-1) floating operations per second (flops). Before the circuit simplification, the PHS contains 32 component variables, 5 voltage supply ports, 2 inputs ports and 1 observation port. Using F s = 768 kHz and n = 40 requires 2.42 Gflops. After the circuit simplification and the frequency reduction, 17 component variables and 3 ports are removed leading to n = 20. For F s = 192 kHz, this requires 149 Mflops (gain of about 16). The simplified PHS is implemented in the JUCE framework using the PyPHS library to generate C++ code [START_REF] Falaize | PyPHS: Passive modeling and simulation in python[END_REF]. The player controls the pitch and sound intensity either with sliders or with commercial midi interfaces. Additionally, the user has the possibility to adjust the quantity of harmonics by changing the demodulator input gain (a feature not present on the original instrument). While running, the plugin consumes 85% CPU (i7-2720QM CPU at 2.20 GHz) and 0.9% RAM (4Go), which is still quite heavy but allows real-time usage.

VIII. CONCLUSION AND PERSPECTIVES

In this paper, a refined simulation of the ondes Martenot circuit has been proposed, allowing a numerical investigation. It is based on a power-balanced modeling adapted to LTI (capacitors, resistor, inductor) and non-LTI (vacuum tubes and the multi-physical time-varying capacitor mechanically-driven by a ribbon) components. Due to the nature of the ondes Martenot instrument (rare, fragile and expensive), setting up extensive measurements to evaluate the accuracy of the circuit modeling is a complex operation which is still ongoing. However, numerical experiments on the modeled circuit lead to the following observations:

1) the combination of the triode vacuum tube amplification properties and the LC-resonator produces a quasisinusoidal oscillation with stable amplitude; 2) the contribution of the power amplifier is minor in terms of coloration compared to that of the preamplifier; 3) the overall result does not depend on the carrier frequency value, provided that this frequency stays outside the audio domain. In consequence, the oscillators are replaced with sinewave generators, the power amplifier is replaced with a multiplying factor, and the sampling frequency is lowered by a factor of 4. These simplifications lead to the implementation of the circuit within a controllable audio plugin, functioning in real-time on a common computer.

Further work aims to model the diffuseurs and the control interfaces (ribbon and intensity key). Moreover, this work will rely on anti-aliasing methods [START_REF] Muller | Trajectory anti-aliasing on guaranteed-passive simulation of nonlinear physical systems[END_REF] to handle the non-linear heterodyne process at lower sampling rates.

Fig. 1 :

 1 Fig. 1: Different parts of the ondes Martenot : diffuseurs (Fig.1aon the right, from[START_REF] Laurendeau | luthier de l'électronique[END_REF]), circuit, ribbon, keyboard (Fig.1b), intensity key with the powder bag underneath circled in red (Fig.1c).

Fig. 2 :

 2 Fig. 2: Complete circuit (top) and isolated stages simplified for modeling (bottom): oscillators (stages 1 and 2, Fig.2b), demodulator (stage 3, Fig.2c), preamplifier (stage 4, Fig.2d) and power amplifier (stage 5, Fig.2e).

Fig. 5 :

 5 Fig. 5: Variable capacitor schematics where d denotes the ribbon position (up to 1.2 m for the highest note) and (d) the distance between the ribbon and the comb teeth (around 1 mm).

  gk ) is estimated by fitting the model to the datasheets [47, Fig. 1] [48, Fig. 1] [49, Fig. 1

1 µFig. 8 :

 18 Fig. 8: Location of s and regimes of oscillation for P ranging 1/µ and 1/µ + β.

Fig. 9 :Fig. 10 :

 910 Fig. 9: From top left to bottom right: simulated output voltages at stages 3 (9a), 4 (9b) and 5 with carrier residual (9c), spectrogram of the simulated output voltage at stage 5 (9d), power balance of the circuit during the simulation (9e-9f).

Fig. 11 :Fig. 12 :

 1112 Fig. 11: Spectral density of the simulated power amplifier output voltage for an input of 1 kHz, amplitude 20 V.

TABLE I :

 I State variables and constitutive laws for a linear parallel RLC circuit.

	x	dx dt	H(x)	∇H(x)

  Table II recaps all components parameters values used for simulation.

TABLE II :

 II Component parameters values used for simulation.

	Triode	µ	Ex	Kg	Kp	K vb
	6F5	98	1.6	2614	905	1.87
	6C5	20	1.5	2837	138	89
	2A3	4.3	1.5	1685	43	102
		V ct	Va	R gk
	6F5	0.5	0.33	1300
	6C5	0.8	0.33	1300
	2A3	-1.2	0.33	1300
	Stage	V bias	C k	R k
	osc	90 V	0.22 µF	7.5 kΩ
	demod	100 V	0.277 µF	1 kΩ
	preamp	180 V	0.277 µF	1 kΩ
	poweramp	230 V	10 µF	750 Ω
	Losc	Cosc	C dem	C 21	R 4	Lp
	7 mH	544 pF	7 nF	200 pF 1 MΩ	9 H
	Rp	Rpw	P	ρ	I out	V start
	4 kΩ	1.5 kΩ	1/µ + β/35	3	0 A	1 mV

TABLE III :

 III Simulation parameters.

	sampling frequency Fs	768 kHz
	fixed osc. frequency F	80 kHz
	playing frequency fm	55 Hz-3520 Hz
	duration	1 s

Note that, according to the component type, the "state flow" can be a current (i = q = ẋ for capacitors with charge q) or a voltage (v = φ = ẋ for induction coils with magnetic flux φ). Thus, the flow is not systematically understood as the "charge flow" that is the electric current.

The measured noise amplitude is around 1 µV RMS[START_REF] Pearson | Fluctuation noise in vacuum tubes[END_REF] for a vacuum tube. However, as the noise is sufficiently filtered out by the oscillator, its amplitude mainly affects the transient duration[START_REF] Edson | Noise in oscillators[END_REF].

IX. ACKNOWLEDGMENTS

This work was funded by the Collegium Musicae and is supported by the project ANR-16-CE92-0028, entitled Interconnected Infinite-Dimensional systems for Heterogeneous Media, INFIDHEM, financed by the French National Research Agency (ANR). Further information is available at https://websites.isae-supaero.fr/infidhem/the-project. The authors thank Thierry Maniguet and Stéphane Vaiedelich (Equipe Conservation Recherche, Musée de la Musique) for their invaluable input.


G Fubiani

Y Jiang

J P Boeuf

E × B electron drift current across the aperture of an ion source surrounded by a cusped magnetic field profile

In negative ion sources, a cusped magnetic field is generated by magnets placed around each aperture of the extraction grid in order to limit the co-extracted electron current. In spite of this suppression magnetic field, the co-extracted electron current is large, on the same order as the negative ion current extracted from the plasma. In this paper we study the mechanisms of electron extraction from the plasma through a cusped aperture in a simplified situation, in the absence of negative ions, with the help of a three-dimensional Particle-In-Cell Monte Carlo Collisions (PIC-MCC) model. The calculation results show that the electron current extracted from the plasma is small for an infinite slit aperture with suppression (cusped) magnetic field and significantly increases in the case of finite slit or circular grid apertures. We find that E × B electron drift plays an important role in the extraction of electrons through a finite slit grid aperture and that current driven micro instabilities are present in the aperture region. This work is relevant to negative ion sources and micro-ECR neutralizers designed for space propulsion.

I. INTRODUCTION

The question of electron extraction mechanisms from negative ion sources (and in particular of those used for fusion applications [1][2][3][4][5][6]) and in micro-ECR neutralizers for space propulsion [7,8] is currently an open problem. In fusion type plasma sources, a magnetic cusp field (also called suppression field) is used to deflect extracted electrons onto the first grid of an electrostatic accelerator in order to prevent the acceleration to high energies of these particles. Electrons are co-extracted with negative ions (about the same current densities for both species). The acceleration voltage is 1 MV in the case of ITER and electrons can in principle damage the accelerator as a result of a high parasitic power density deposition on the accelerator grids [9,10].

The magnetic field from the suppression magnets penetrates also inside the ion source plasma, a few centimetres in front of the plasma grid (PG). The latter has 1280 apertures in the case of ITER. We focus in this work specifically on the effect of the cusp magnetic field on the electron dynamics in that area (the cusp surrounds each apertures). The numerical model is hence simplified and we did not include negative ions. Their influence is briefly analysed in the discussion section. Furthermore, we do not simulate the whole ion source but rather a small area around the apertures. We also model exclusively slit apertures (the generalisation to circular geometry is straightforward). The plasma is generated artificially (a source term for the charged particles) on the left-hand-side (LHS) of the simulation domain. The first grid of the accelerator, which is commonly termed "extraction grid" (EG), is included in the calculation and is biased positively with respect to the PG. The bias voltage extracts an electron current from the plasma and induces the formation of a meniscus around each apertures. The meniscus separates the quasineutral plasma from the non-neutral region inside the accelerator. The cusp field amplitude decreases exponentially from the permanent magnet position, which is located in the model on the EG (same as in ITER and DEMO accelerator designs). Note that electrons are strongly magnetized while ions are somewhat magnetized but only over a short distance close to the aperture surface.

In the model, the cusp field is generated by two permanent magnet bars of reversed polarity parallel to the slit apertures. The magnetic vector field hence lies approximatively in the plane perpendicular to the magnets when the latter are much longer than the length of the aperture. We find that the electron flux diffusing across the aperture is driven by an E × B drift where E is the electric field in the vicinity of the plasma meniscus. The E × B drift occurs in the plane perpendicular to the cusp magnetic field lines. As a consequence, in an hypothetical configuration where a slit aperture is of infinite length, this E × B drift does not contribute to electron extraction. The flux of electrons through the aperture is much smaller in that case and is only due to collisional cross-field drift. We show that a 2D Particle-In-Cell (PIC) algorithm with Monte-Carlo-Collisions (MCC) is sufficient to capture the underlying physics of an infinite slit, and is consistent with results from a 3D PIC-MCC model. In the case of a finite slit, only a 3D model can describe properly electron extraction due to E × B drift through the aperture. Note that plasma electron extraction from a negative ion source has been previously studied analytically and with a ray-tracing approach in 1D and 2D by Cavenago et al. [START_REF] Cavenago | AIP Conference Proceedings[END_REF].

In the next section, we describe in details the simulation parameters. In sec. III, we analyse the mechanisms of electron extraction in the case of an infinite slit and for a finite slit aperture. In sec. IV, we describe some micro-instabilities occurring in the vicinity of the plasma meniscus in the case of a finite slit aperture. Lastly, sec. V and VI, provide a discussion and conclusions on the work presented in this paper.

II. OVERVIEW OF THE 2D AND 3D PIC-MCC MODELS

We developed 2D and 3D PIC-MCC electrostatic explicit models [START_REF] Fubiani | [END_REF]. The electric field is calculated on the mesh nodes by solving Poisson's equation iteratively via a multi-grid algorithm [START_REF] Briggs | A Multigrid Tutorial[END_REF]. The magnetic cusp field profile is derived from analytical considerations [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF],

B x = B 0 sin π (z -z 0 ) d exp - π (x 0 -x) d , (1)
B z = B 0 cos π (z -z 0 ) d exp - π (x 0 -x) d , (2)
where B 0 is the maximum amplitude of the field, d is the distance between the magnet bars,

x 0 is their axial coordinate and z 0 the center position between the two magnets, respectively.

The norm of the magnetic field vector, |B|, is decreasing exponentially from the location of the magnets. The particles trajectories are calculated by solving the Newton's equation of motion at the particle location (the electric and magnetic fields are interpolated),

m p dv p dt = q p (E + v p × B) , (3)
where p either stands for electrons or ions, E(x, y, z) is the electrostatic field, v p (t) the particle velocity and m p (q p ) its mass (charge), respectively. The force on the RHS is the Lorentz force. In a PIC-MCC algorithm, charged particles are actually macro-particles but with the same charge over mass ratio as real electrons and ions. Once the particle positions have been updated, the density is calculated on the mesh nodes. This is used as a source term for Poisson's equations. In addition, we implemented collisions between charged particles and the neutral background gas using the Monte-Carlo technique [START_REF] Fubiani | [END_REF][START_REF] Nanbu | [END_REF][16][17]. The PIC-MCC cycle described above is repeated iteratively. The collision processes are simplified in this problem. We considered only the elastic collisions between electrons and the positive ions or the neutrals (table I). The neutral background gas profile is assumed constant with a density of n n = 5 × 10 19 m -3 and the gas temperature is set at T n = 0.1 eV. The average electronneutral elastic collision frequency in the calculation is ν en ≃ 5 × 10 6 s -1 (reaction #1 of table I) while the electron-ion collisions frequency is of order of ν ei ≃ 1.2 × 10 5 s -1 (reaction #2), respectively. The average plasma density is n i = 5×10 15 m -3 and the plasma generation is modelled by mimicking ionisation through the reinjection of an electron-positive ion pair for each positive ion lost on the physical boundaries of the simulation domain (i.e, of Dirichlet type). Note that the electron Debye length is typically about eight times larger than in the extraction region of negative ion source prototypes developed for ITER and DEMO [3,18].

The particles are injected into the gray area shown schematically on the LHS of fig. 1 (dubbed "plasma source"). The electron distribution function is artificially maintained in that area as a Maxwellian at a temperature of T e = 2 eV via the periodic replacement of the electron velocity (we use an artificial "heating" frequency). The positive ions are injected assuming also a Maxwellian distribution at a temperature T i = 2 eV (although we do not enforce it in the plasma volume, only at the time when the ions are generated). The latter is the typical value for the average kinetic energy deduced from 3D PIC-MCC models of a fusiontype negative ion source [19]. We implemented a numerical resolution of 256 × 384 × 192 grid nodes associated with a simulation box of physical dimensions 3.2 × 4.8 × 2.4 cm 3

(length versus width and height). The number of particles per cell is 40 and the time step

ω p ∆t ≃ 0.2.
In addition, the size of the two slits is 4.3×0.8 cm 2 (length times width) and the amplitude of the cusp magnetic field on the plasma grid (PG) is 200G, which corresponds to B 0 = 875G, d = 1.2 cm, x 0 = 3.2 cm, z 0 = 0.6 cm in eqs. ( 1) and ( 2), respectively.

The cusped magnetic field profile encompasses each of the slit apertures, as shown on the RHS of fig. 1. The direction of the magnetic field vector is reversed between neighbouring apertures. In addition, the numerical model is parallelized in an hybrid manner using the OpenMP and MPI libraries [START_REF] Fubiani | [END_REF]. For the simulations reported in this paper, we typically Next, we compare the 2D calculation to a 3D PIC-MCC model where the simulation domain is identical to the one of fig. 1 except that the slit aperture is of infinite length along (Oy); its width being 8 mm in the z-direction. Any occurence of plasma instabilities in the plane perpendicular to the magnetic field lines, which could potentially increase the electron current transmitted across the aperture, will be modelled in this case. We found instead that the electron current collected on the extraction grid (EG) is identical in 2D and 3D, as shown in fig. 3. The latter plots the ratio (in percent) of the electron current impacting the EG divided by the total current flowing toward the PG inside the plasma volume. The level of numerical noise is logically higher in the 2D calculation due to the larger weight of the macroparticles. To further support the observation that the plasma properties are identical in 2D and 3D in the case of slit apertures of infinite length, we show in fig. 4 the plasma potential profile along (Ox) for a location on the PG in the middle of the interstice between slit apertures. We compare the results from the 2D calculation versus estimates from the 3D PIC-MCC algorithm for either the case with two apertures or using a reduced domain in 3D where we modelled only one cusp, i.e., half the domain of fig. 1 along (Oz).

B. 3D PIC-MCC model of a slit aperture of finite length

We now analyse the electron kinetics for the case of a slit aperture of finite length. The simulation domain is described in details in sec. II and fig. 1. Figure 5 shows together the electron flux profiles in (i) the plane of the magnetic cusp field lines (XZ plane at Y = 2.4 cm),

(ii) the plane where the electron drift occurs (XY plane at Z = 0.6 cm, crossing the middle of the slit aperture), (iii) the front of the PG and, lastly, (iv) the EG, respectively.

In (i), one can see the telltale sign of the electrons oscillating along the cusp field lines while in (ii), the electrons drift in the direction perpendicular to the magnetic field (E × B drift, see fig. 6 for an illustration using individual particle trajectories) and they can exit toward the accelerator only through one side of the slit aperture. Figure 6(a) shows one electron trajectory in the Z = 0.6 cm plane assuming that the Z component of its velocity is zero. When the electron reaches the edge of the plasma grid aperture, the electric field increases rapidly due to the large voltage of the extraction grid. This large non-uniformity of the electric field generates an electron drift across the magnetic field, i.e. the electron guiding center is no longer confined along an equipotential line as can be seen in fig. 6(a). This is a finite Larmor radius effect due to the electric field non uniformity (see for instance FF. Chen [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF]).

In the case of fig. 6 The E × B drift directs the electron toward the PG where it impacts its surface (green color)

and lastly, (iii) the electron drifts in the vicinity of the plasma meniscus and is extracted toward the EG (red color). Note that it may also drift some distance inside the accelerator before impacting the EG [9].

Comparing the case of a slit aperture of infinite (described in sec. III A) versus finite length, one observe than in the latter the electron current extracted toward the accelerator vessel is greatly enhanced by the E × B drift motion in front of the plasma meniscus. The ratio of the electron current impacting the EG divided by the total current flowing into the PG area (inside the plasma volume) is shown in fig. 7. We find that the current on the EG strongly varies versus time within a range for the transmission ratio between 50% and 100% (i.e., all the current). This is to be compared to the ∼ 22% in the case of a slit aperture of infinite length. Such a large extracted current is induced by a plasma potentiel significantly higher than the average electron kinetic energy in the model as shown in fig. 4 (the temperature is T e = 2 eV).

IV. UNSTABLE PLASMA BEHAVIOUR IN THE VICINITY OF THE PG

The large time variations of the current transmitted toward the EG in fig. 7 We can conclude that the current oscillations observed in the simulations are due to a drift instability induced by the charge separation associated with the electron loss to the extraction grid.

V. DISCUSSION

We modelled in this work a slit aperture of either finite or infinite length. We also performed simulations for cylindrical apertures and reached similar conclusions. In addition, for the case of negative ion sources for fusion applications, there is typically a large negative ion to electron ratio in the vicinity of the aperture. In Christ-Koch et al. for instance [START_REF] Christ-Koch | [END_REF] the measured ratio n -/n e is about 0.4 in hydrogen and as high as 1.5 in deuterium 2 cm from the PG, inside the plasma volume. In these sources, the negative ions are produced on the surface of the PG as a byproduct of the impact of neutral hydrogen or deuterium atoms on the grid. The latter is covered by a thin layer of caesium to lower the work function of the metal increasing, as a consequence, the negative ions production yield. The adjunction of negative ions in our model did not modify the conclusions that an E × B drift near the plasma meniscus drives the electron current across the aperture toward the EG. The plasma density in the extraction region of these ion sources is typically sixty times larger than in our model, i.e., n i ≃ 3 × 10 17 m -3 , and hence the electron-ion Coulomb collision and the electron-neutral collision frequencies are of the same order of magnitude. The E × B drift remains the dominant electron extraction mechanism in our simulations when we artificially increase ν ei such that ν ei ∼ ν en .

Lastly, in sec. III B, we found that a high fraction of the electron current flowing from the plasma volume toward the PG may be extracted (fig. 7). The extracted current can be significantly lowered in the model by biasing the PG with a voltage positive with respect to the ion source walls. In that case, the potential difference between the biased PG and the sheath potential is reduced and hence the electron current collected on the grid may be significantly increased. This technique is commonly employed in fusion-type ion sources used for Neutral Beam Injection (NBI) [18] .

VI. CONCLUSION

We analysed the transport of electrons across a grid with a slit aperture surrounded by a magnetic cups field profile. An electrode biased positively with respect to the grid is placed behind in order to produce an electric field which can extract negatively charged particles.

In the case of negative ion sources, the role of the cusp field is to act as a barrier, preventing electrons from being extracted. The slit was either of finite or infinite length in the model. In the former, one found identical results between 3D and 2D PIC-MCC calculations indicating that the transport across the aperture is solely collisional. The plane which includes the cusp field lines was simulated in the 2D PIC-MCC model, hence any transport due to the magnetised drift dynamics was not included. In the 3D case, we found that the electron current transmitted onto the bias electrode (EG) was greatly enhanced compared to the 2D simulations. This is due to an E × B drift motion for the electrons where E is the electric field of the plasma meniscus. The electrons escape parallel to one of the two edges of the slit aperture, being the only location where the drift velocity is directed toward the EG (the drift is generated by the B z and E y components in our configuration). Furthermore, we observed that the plasma parameters (potential, density, electron flux) are strongly fluctuating in front of the plasma meniscus resulting in a transmitted electron current ratio which varies within a factor of two over a time scale of 1 µs. The maximum corresponds to an extraction of most of the electrons flowing onto the PG from the plasma volume in the conditions that we modelled. These oscillations are microinsatbilities due to the charge separation induced by the loss of electrons to the extraction grid.

  used 48 OpenMP threads together with 4 MPI threads (one per node), i.e, 192 cores in total (on Intel Xeon CPU E5-2697 v4 @ 2.30GHz). The top, bottom, left and right boundary conditions are periodic and we are hence modelling the equivalent of an infinite number of slit apertures. Lastly, the space between the PG and the extraction grid (EG) is free of plasma (the gap is 3.5 mm and the width of the PG is 2 mm). The LHS of the simulation domain, PG, and EG are of Dirichlet type. The EG is biased at a voltage of 400V in the simulations while the other absorbing boundaries (PG and LHS) are set at 0V.III. DESCRIPTION OF THE ELECTRONS MOTION ACROSS THE APER-TURES A. 2D versus 3D PIC-MCC calculation of a slit aperture of infinite lengthIn this section, we compare the plasma characteristics derived from a 2D versus 3D PIC-MCC model of a slit aperture of infinite length. The simulation domain for the 2D calculation is shown in fig.2. We model the plane which includes the cusp field lines, i.e., the XZ plane of fig.1. The particle drifts are out of the simulation domain and drift-wave type plasma instabilities are hence not described. The electrons which are extracted toward the accelerator diffuse across the field lines solely through collisions.

  fig.6(b) and (c) are calculated in the 3D time averaged plasma potentiel. Electrons move both along the magnetic field lines, i.e, v z = 0, and drift toward the right edge of the aperture. B z decreases when the electron moves either upward in the direction of increasing Z or downward meanwhile the X component of the magnetic field, which is perpendicular to the EG, increases. The reason why the maximum electron current on the PG is seen on the two right corners of the aperture (see fig.5) is therefore due to the combined effect of

  are due to strong fluctuations in the plasma quantities in the vicinity of the plasma meniscus. Figures8 and 9plots the electron flux and plasma potential profiles versus several time shots which are also displayed in fig.7for reference. The time window corresponds to approximatively two periods of oscillations. One can see that the flux profile is strongly modified, which is correlated with stiff changes of the plasma potential. We observe that (i) the electron flux follows approximately the potentiel lines in the aperture mid-plane plane as discussed previously in sec. III B (see also fig.6), (ii) case A and I (D et G), which correspond to a maximum (minimum) transmission rate, share similar features, respectively. In D et G, the isopotential lines are twisted and we observe the appearance of potential wells (which attract positively charged ions). The well forms due to a local increase of the electron density over a short time scale. The area of the well is accompanied by a double layer, i.e., over a short distance (∼ λ De ), the space charge (n e /n i -1) switches from globally positive (excess of electrons) to negative. The direction of the electron flux is upward in fig.8and electrons experience a strong E × B drift in the vicinity of the aperture top edge in cases D et G (the 10V and 20V isopotential lines are very close in that area, hence the amplitude of the electric field is large, of order 20 kV/m for a magnetic field strength of 200G and a drift velocity v D ∼ 10 6 m/s, respectively). Such a situation is accompanied by a surge in transmission rate, as shown in fig. 7. In cases A and I, the electric field strength (and hence the drift velocity) is much lower in front of the aperture, which explains the drop in the electron current extracted toward the EG.

FIG. 1 . 9 FIG. 2 .

 192 FIG. 1. Schematic illustration of the simulation domain. On the left-hand-side (LHS) is the plasma volume, which is separated from a biased electrode (EG) by a grid (PG) with two slit apertures (light grey). The volume between the PG and the EG is free of plasma. The latter is generated numerically by the injection of electron-positive ion pairs in the gray area labelled "plasma source". We fixed the plasma density in the model and hence each positive ion lost on a physical boundary (Dirichlet) is compensated by the re-injection of an electron-ion pair. The magnetic cusp field profile on the PG (in yellow) is shown in the right-hand-side (RHS) of the figure (the slit aperture is highlighted in black). In practice, this field would be generated by permanent magnet bars of reversed polarities (blue and red colours, N stands for "north" and S for "south"). The boundary condition on the LHS together with the plasma grid (PG) and extraction grid (EG) are of Dirichlet type while the top/bottom and left/right boundaries are periodic in 3D. The simulation box is 3.2 × 4.8 × 2.4 cm 3 (x, y, z) corresponding to a numerical resolution of 256 × 384 × 192 grid nodes. The dimensions of the slit apertures are 4.3 cm along (Oy) and 0.8 cm in the z direction. The potential of the walls are 0V on the LHS plane x = 0 and PG (between x = 2.65 cm and x = 2.85 cm) together with V EG = 400V (at x = 3.2 cm), respectively.

FIG. 3 .FIG. 4 .

 34 FIG. 3. Electron current collected on the extraction grid (EG). 2D versus 3D PIC-MCC calculationfor a configuration with a slit aperture of infinite length.

FIG. 6 . 3 .

 63 FIG. 6. Examples of electron trajectories, (a) in the Z=0.6 cm plane only (see fig. 5), (b) and (c) trajectories calculated in the 3D time averaged potential and projected on (a), the XY plane, and (b), the XZ plane, illustrating three typical cases. Blue trajectory: the trajectory is trapped and the electron oscillates back and forth in the cusp magnetic field. Collisions are needed to un-trap the electron and eventually move it to a green or red type of trajectories. Green trajectory: the electron undergoes E × B drift that draws it to the plasma grid. Red trajectory: the electron undergoes E × B drift. Due to the increasing X component of the electric field, and decreasing Z component of the magnetic field in the vicinity of the plasma grid (PG), the electron is no longer trapped by the magnetic field and can escape to the extraction grid (EG). The maximum current of extracted electrons is located in the top Y region of the PG due to E × B drift (direction shown in (a)) and in the top and bottom Z corners of the PG because the Z component of the magnetic field decreases in these cusp regions. The initial conditions (index 0) of the four trajectories are as follows. Purple trajectory: initial positions x 0 = 1.5 cm, y 0 = 2.4 cm, initial kinetic energy E 0 = 2 eV, initial components of the velocities v x0 = eE 0 /m e and v y0 = -v x0 , respectively. Blue trajectory: x 0 = 2 cm, y 0 = 4 cm, z 0 = 0.6 cm, kinetic energy (in eV) E 0 = m e /2ev 2 0 = 3 eV, velocities v x0 = v y0 = v z0 = v 0 √ 3. Green trajectory: x 0 = 2 cm, y 0 = 0.3 cm, z 0 = 0.6 cm,E 0 = 2 eV, v x0 = -v y0 = v 0 √ 0.4, v z0 = v 0 √ 0.2. Red trajectory: x 0 = 2.3 cm, y 0 = 2 cm, z 0 = 0.7 cm, E 0 = 2 eV, v x0 = v y0 = v 0 √2, v z0 = 0. The isopotential contours (black solid lines) are displayed in (a) together with the amplitude of the potential (grey scale).

FIG. 7 .

 7 FIG. 7. Ratio of the electron current collected on the EG divided by the total current flowing toward the PG, I EG / (I EG + I P G ) (in percent) versus time. I EG (I P G ) is the current impacting the EG (PG), respectively.

FIG. 8 .FIG. 9 .

 89 FIG. 8. Electron flux profile in the plane Z = 0.6 cm of fig. 5 for different time shots labeled from A to I (which are also indicated in fig. 7 for reference). The isopotential lines for φ = 10V and 20V inside the plasma volume are also shown. Γ max = 2.3 × 10 21 m -2 • s -1 .

TABLE I .

 I Particle processes.

	#	Reaction	Cross section ref.
	1	e + H 2 → e + H 2 (elastic)	[22]
	2	e + H + 2 → e + H + 2

ACKNOWLEDGMENTS This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Support from CEA and from the French Research Federation for Fusion Magnetic Confinement is also acknowledged.

AVAILABILITY OF DATA

The data that support the findings of this study are available from the corresponding author upon reasonable request.

AUTHOR'S CONTRIBUTIONS

All authors contributed equally to this work.