Aurélien Drezet

A dual-Lagrangian description adapted to quantum optics in dispersive and dissipative dielectric media

Keywords: numbers: 42.50.Ct, 41.20.Jb, 73.20.Mf

We develop a dual description of quantum optics adapted to dielectric systems without magnetic property. Our formalism, which is shown to be equivalent to the standard one within some dipolar approximations discussed in the article, is applied to the description of polaritons in dielectric media. We show that the dual formalism leads to the Huttner-Barnett equations [B. Huttner, S. M. Barnett, Phys. Rev. A 46, 4306 (1992)] for QED in dielectric systems. More generally, we discuss the role of electromagnetic duality in the quantization procedure for optical systems and derive the structure of the dynamical laws in the various representations.

I. INTRODUCTION

The description of quantized electromagnetic fields propagating in a dielectric media has been the subject of intensive research since the early times of quantum optics. It was shown already in 1948 [1] that an Hamiltonian formalism adapted to non dispersive and dissipative model can be easily implemented by analogy with standard QED quantization procedures used in vacuum. With the development of quantum optics and integrated photonics several authors have extended the previous canonical quantization method to complex inhomogeneous systems with a spatially dependent dielectric permittivity ε(x) (see for example Refs. [2][3][4]). However, in order to consider the physical properties of realistic optical media the central issue is the inclusion of dispersion and dissipation in the quantum description [5]. Since these properties are connected by the well-known Kramers-Kronig relations [START_REF] Jackson | Classical Electrodynamics 3rd[END_REF][START_REF] Landau | Electrodynamics of Continuous Media[END_REF] related to causality it is clear that a correct description should not analyze dispersion and dissipation separately. Nevertheless, when dissipation is weak it is actually possible to give an Hamiltonian basis to the quantum formalism [START_REF] Constantinou | [END_REF][9][10][11] using relations obtained long ago by L. Brillouin [START_REF] Landau | Electrodynamics of Continuous Media[END_REF][START_REF] Brillouin | Wave propagation and group velocity[END_REF] for the energy density in dispersive media at frequency ω, at which absorption is negligible. These approaches involve 'quasimodal' expansions and were recently applied to the field of quantum plasmonics [START_REF] Tame | [END_REF] for the coupling of emitters to nano antennas [14][15][16] when dissipation is weak. Clearly these assumptions are however not generally valid, e.g., in plasmonics [START_REF] Raether | Surface Plasmons[END_REF][START_REF] Barnes | [END_REF][19], where optical resonances whith low quality factors are coupled to fluorescent emitters [20][21][22][23][24][25][26][27][28][29].

However, in the 1990's it already became clear that fundamental progresses could only occur by including more degrees of freedom in the quantum description, i.e., by modeling dissipation with thermal baths coupled to the photonic variables. The approach was used in quantum optics for the description of lossy beam splitters [30,31], for instance for modeling Casimir forces [32], and more recently was applied in quantum plasmonics by modeling decay channels as an information loss in the metal/dielectric environment [33][34][35]. Moreover, the most important progress came in 1992 when it was suggested by Huttner and Barnett to model any dispersive and absorbing linear bulk media satisfying Kramers-Kronig relations by using a modified Hopfield-Fano [36,37] model involving a continuous distribution of harmonic oscillator fields interacting with light [38][39][40][41][42][43][START_REF] Barnett | The quantum optics of dielectrics[END_REF]. This powerful approach was originally limited to infinite and homogeneous bulk media but continuous efforts have been done later on to extend its range of validity to the more interesting inhomogeneous medium case, which is needed for nanophotonics applications in general and for quantum plasmonics in particular [45][46][47][48]51]. Actually these formal developments of the original Huttner-Barnet model [39] are strongly motivated by the parallel development of the phenomenological Green function approach of Gruner and Welsch [52][53][54] based on the Langevin equation method [START_REF] Cohen-Tannoudji | Atom-Photon Interactions: Basic Processes and Applications[END_REF] and the dyadic Green function formalism [START_REF] Girard | [END_REF][START_REF] Novotny | Principles of Nano-Optics[END_REF]. In this strategy a noise current is added phenomenologically to Maxwell's equation in order to preserve unitarity of the full evolution and in particular the constancy of all conjugate canonical variables commutator with time [START_REF] Scheel | [END_REF]. The Green function strategy has been intensively used in the literature in the recent years [59][60][START_REF] Knöll | QED in dispersing and absorbing media[END_REF][START_REF] Matloob | [END_REF][63][64][65][66][67] and applied to several problems including dielectric or magnetic materials, and coupling of light with atoms in the regime of weak or strong coupling in presence of plasmonic nanoparticles [68][69][70][71][72][73][74][75][76]. Despite its success the Langevin noise method applied to macroscopic electrodynamics (unlike for atomic physics in vacuum [START_REF] Cohen-Tannoudji | Atom-Photon Interactions: Basic Processes and Applications[END_REF]) lacks a neat and clear quantum foundation that, like the Huttner-Barnett model [39], could be justified using an Hamiltonian description. In Refs. [45][46][47][49][50][51] general proofs were given for the formal equivalence between the Huttner Barnett model and the Langevin noise approach.

The present work contributes to this discussion by providing a different Lagrangian and Hamiltonian foundation to the approach considered in, e.g., Ref. [51]. More precisely, our aim is to bypass the usual canonical procedure based on the 'minimal coupling' Lagrangian in the Coulomb gauge [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF]. Generally speaking this is usually done in the literature by introducing the Power-Zienau unitary transformation [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF][START_REF] Power | [END_REF] leading to the multipolar representation of the electromagnetic field [START_REF] Craig | Molecular Quantum Electrodynamics[END_REF]. Here, we proceed differently by introducing an equivalent description of the electromagnetic system using a dual electric vector potential F(x, t) different from the usual magnetic vector potential A(x, t). With this formalism there is no unphysical separation between transverse quantized and longitudinal unquantized fields (indeed, this separation is no relativistically causal [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF]). The method of quantization we propose here is actually a direct way to introduce the multipolar Hamiltonian [START_REF] Power | [END_REF][START_REF] Craig | Molecular Quantum Electrodynamics[END_REF], where the fundamental electromagnetic quantities are now the displacement D and magnetic B fields. In turn this approach, adapted to neutral dielectric systems without magnetic properties, will shed some new light on the Huttner-Barnett model (which was based on the minimal coupling scheme) and therefore on the derivation of the Langevin Noise method [52][53][54]. We conclude the article by a discussion concerning electromagnetic duality between electric and magnetic physical quantities and compare different but equivalent representations of the dynamical laws.

II. THE GENERAL HUTTNER-BARNETT EFFECTIVE MODEL FOR A LINEAR DIELECTRIC MEDIUM

A. A Lagrangian formulation

We here follow the standard approach in quantum electrodynamics and start with writing a Lagrange-Hamilton principle adapted to our system. In other words, we choose the Lagrangian L(t) and the Lagrangian density L(x, t) such that the variational problem δ t2 t1 dtL(t) = δ t2 t1 dt d 3 xL(x, t) = 0 could be solved and lead to the Maxwell equations in a dissipative and dispersive inhomogeneous dielectric medium. However, in this work we will not consider the usual canonical 'minimal coupling' Lagrangian in the Coulomb gauge [START_REF] Cohen-Tannoudji | Atom-Photon Interactions: Basic Processes and Applications[END_REF] which was used historically by Hutner and Barnett [39] but instead work with the dual choice (in the Lorentz-Heaviside system of units):

L = B 2 -D 2 2 + F • ∇ × P - P 2 2 + L M . (1)
Here, the Lagrangian density depends on the electric vector potential F(x, t) ∈ R 3 and its derivatives defined as

B(x, t) = 1 c ∂ t F(x, t), D(x, t) = ∇ × F(x, t). (2)
B(x, t) and D(x, t) being the magnetic and displacement fields respectively. Additionally here we will work exclusively in the 'Coulomb' gauge ∇ • F(x, t) = 0. As detailed in the Appendix A and discussed in the section IV the Lagrangian is motivated by fundamental duality relations existing between electric and magnetic variables. This in turn explains why the magnetic field term B 2 /2 appears with a positive sign in Eq. 1 while D 2 /2 appears with a negative sign contrarily to the usual choice made in the minimal coupling representation. We point out that this duality symmetry was discussed long ago in a more classical context involving for example chirality or bi-anisotropic materials [START_REF] Stratton | Electromagnetic Theory[END_REF][START_REF] Serdyukov | Electromagnetics of Bi-anisotropic Materials (Gordon and Breach[END_REF][START_REF] Fernandez-Corbaton | [END_REF], here we applied this concept to the Lagrange-Hamilton language adapted to second quantization. Moreover, in this formalism, the material part L M reads

L M = +∞ 0 dω (∂ t X ω ) 2 -ω 2 X 2 ω 2 . ( 3
)
where X ω (x, t) defines the material oscillator vectorial fields which are coupled to the electromagnetic variables. These fields are labeled by the pulsation ω varying continuously from zero to +∞. We emphasize that in the total Lagrangian density we have also a term -P 2 2 which plays the role of an internal interaction and which is not included in L M to respect the conventions used in Ref. [39]. This choice for L leads to the same dynamical equations as those deduced from the usual formalism based on the minimal coupling Lagrangian [START_REF] Cohen-Tannoudji | Atom-Photon Interactions: Basic Processes and Applications[END_REF]. However, the dual formalism appears more convenient for the present study of neutral dielectrics system since it does not actually involve a non physical separation of the electromagnetic field between a transverse (quantized) part associated with photons and a longitudinal (un-quantized) part associated with non-causal Coulomb Fields (see later on this paper). The present description can actually be seen as an alternative way to access to the multipolar representation which is connected to the minimal coupling one through a unitary transformation as shown by Power and Zienau [START_REF] Cohen-Tannoudji | Atom-Photon Interactions: Basic Processes and Applications[END_REF][START_REF] Power | [END_REF][START_REF] Craig | Molecular Quantum Electrodynamics[END_REF]. Here the use of the F potential vector allows us to deduce similar results by introducing different canonical variables. From the previous definitions follow directly (i.e., using the Euler-Lagrange formalism) the two Maxwell equations:

∇ × B(x, t) = 1 c ∂ t D(x, t), ∇ • D(x, t) = 0. ( 4
)
The matter field X ω (x, t) ∈ R 3 is related to the polarization density by the formula

P(x, t) = +∞ 0 dω 2σ ω (x) π X ω (x, t) (5)
where the coupling function σ ω (x) ≥ 0 will represent the conductivity of the medium at the harmonic pulsation ω.

The fact that the conductivity is sufficient for describing the field will be shown to be consistent with Kramers-Kronig relations [39]. We emphasize that the part L M of the Lagrangian is not exactly the same as the one used in Hutner-Barnett's theory in which three contributions associated with the electromagnetic field, the dielectric and the oscillator bath were included. Here, we use the simpler model proposed by Philbin [51] where the electromagnetic field characterized in our approach by F is directly coupled to the oscillator bath variables X ω without including the additional and non necessary mechanical oscillator used in refs. [38][39][40][41][42][43][START_REF] Barnett | The quantum optics of dielectrics[END_REF].

From Eqs. 1,3 we deduce straightforwardly the Euler-Lagrange equations for the electromagnetic fields, i.e., the two missing Maxwell equations:

∇ × E(x, t) = - 1 c ∂ t B(x, t), ∇ • B(x, t) = 0 (6)
with the electric field

E(x, t) = D(x, t) -P(x, t). (7)
For the matter oscillator field X ω (x, t) we similarly obtain the second-order Euler-Lagrange equation

∂ 2 t X ω (x, t) + ω 2 X ω (x, t) = 2σ ω (x) π E(x, t), (8)
which contains a linear coupling to the local electric field (there is no magnetic coupling such as a Lorentz force). The term -P 2 2 in the Lagrangian density Eq. 1 is necessary for deriving this dynamical equation and obtaining a coupling proportional to E(x, t) and not to D(x, t) (inversely, in absence of coupling with D(x, t) the electric field in Eq. 8 is -P(x, t) which can be seen as an internal force acting on the oscillators X ω (x, t)).

Furthermore, in order to use the Hamilton formalism needed for quantization we define the canonical momenta associated with the different fields:

Π Xω = δL δ(∂ t X ω ) = ∂L ω ∂(∂ t X ω ) = ∂ t X ω Π F = ∂L ∂(∂ t F) = B c (9)
where δ[...] δ(∂tXω) is a functional derivative and where L = +∞ 0 dωL ω . This allows us to introduce the full Hamiltonian defined as

H = d 3 x[Π F ∂ t F + +∞ 0 dωΠ Xω ∂ t X ω )] -L.
We get:

H(t) = d 3 x[ B 2 + D 2 2 -F • ∇ × P + P 2 2 ] + H M = d 3 x[ B 2 + D 2 2 -D • P + P 2 2 ] + H M = d 3 x B 2 + E 2 2 + H M (10)
with the helpful condition:

d 3 xD • P = d 3 xF • ∇ × P,
and where the material contribution H M (t) reads:

H M (t) = d 3 x +∞ 0 dω (∂ t X ω ) 2 + ω 2 X 2 ω 2 . ( 11
)
For the purpose of canonical quantization the full Hamitonian can be written as a functional of the conjugate canonical variables and we deduce:

H(t) = d 3 x[ c 2 Π 2 F + ∇ × F 2 2 -F • ∇ × P + P 2 2 ] + d 3 x +∞ 0 dω (Π Xω ) 2 + ω 2 X 2 ω 2 (12)
At this stage it is useful to deduce Eq. 10 using a different approach. Indeed, the Hamiltonian used here is an integral form of the local energy-tensor conservation law associated with Noether's theorem and it can be preferable, for the sake of generality, to use such a local approach instead of a global one. Starting from Maxwell's equation obtained with Eqs. 1, 3 we get the Poynting conservation theorem

-∂ t ( B 2 + E 2 2 ) = ∇ • (cE × B) + J • E, (13)
which involves the local electric current J = ∂ t P associated with the polarization density P. Replacing P by its value in Eq. 5 and after direct integration leads to

J • E = ∂ t ( +∞ 0 dω (∂tXω ) 2 +ω 2 X 2 ω 2
) and therefore to the local conservation law

-∂ t u = ∇ • (cE × B), (14)
where u, the local energy density, is given by u =

B 2 +E 2 2 + +∞ 0 dω (∂tXω ) 2 +ω 2 X 2 ω 2
. Integration over an infinite volume leads to the Hamiltonian H(t) = d 3 xu(x, t) which becomes equivalent to the total conserved energy if the Poynting vector flow Σ∞ cE × B • dΣ over the infinite closed surface Σ ∞ surrounding our system vanishes sufficiently well.

B. From the Lagrangian to the polarisability of the linear dielectric medium

We should now summarize briefly the consequence of the Lagrangian choice Eq. 1 and show that it allows us to justify the usual Maxwell equations in a causal dielectric medium satisfying Kramers-Kronig relations. To do that we have to solve our dynamical coupled equations for matter and electromagnetic fields. In this section we start with the material field which is the easiest part. In order to solve the evolution Eq. 8 we introduce here other field variables:

Z (±) ω (x, t) = ∂ t X ω (x, t) ± iωX ω (x, t) = (Z (∓) ω (x, t)) * (15)
which obey to the following first order equations:

∂ t Z (±) ω (x, t) = ±iωZ (±) ω (x, t) + 2σ ω (x) π E(x, t). ( 16
)
These equations are easily solved using the method of the variation of constants and lead to

Z (±) ω (x, t) = Z (±) ω (x, t 0 )e ±iω(t-t0) + 2σ ω (x) π t-t0 0 dτ e ±iωτ E(x, t -τ ) = Z (±,0) ω (x, t) + 2σ ω (x) π t-t0 0 dτ e ±iωτ E(x, t -τ ) ( 17
)
with t 0 an initial time which eventually could be sent infinitely in the remote past, i.e., if t 0 → -∞.

We also introduced the 'free' solution Z

(±,0) ω (x, t) = Z (±)
ω (x, t 0 )e ±iω(t-t0) corresponding to the harmonic oscillation of the medium in the absence of coupling.

We can of course go back to the initial field variables by using the transformation X ω (x, t) =

Z (+) ω (x,t) 2iω
+cc. and

∂ t X ω (x, t) = Z (+) ω (x,t) 2
+ cc. This leads to

X ω (x, t) = X (0) ω (x, t) + 2σ ω (x) π t-t0 0 dτ sin ωτ ω E(x, t -τ ) ( 18
)
with the 'free' solution

X (0) ω (x, t) = Z (+,0) ω (x,t) 2iω + cc. (i.e., X (0) ω (x, t) = cos (ω(t -t 0 ))X ω (x, t 0 ) + sin (ω(t-t0)) ω ∂ t X ω (x, t 0 )
). The polarization density now reads

P(x, t) = P (0) (x, t) + t-t0 0 χ(x, τ )dτ E(x, t -τ ) (19)
with the free dipole density distribution:

P (0) (x, t) = +∞ 0 dω 2σ ω (x) π X (0) ω (x, t) (20)
and the linear susceptibility:

χ(x, τ ) = +∞ 0 dω 2σ ω (x) π sin ωτ ω Θ(τ ) (21)
For convenience we introduced the Heaviside unit step function Θ(τ ), defined as Θ(τ ) = 1 unless τ < 0 whereas Θ(τ ) = 0. Now, we get for the displacement field

D(x, t) = E(x, t) + P(x, t) = P (0) (x, t) +E(x, t) + t-t0 0 dτ χ(x, τ )E(x, t -τ ). ( 22
)
If as usually we consider the limit t 0 → -∞ we can write

D(x, t) = P (0) (x, t) + +∞ -∞ dτ ε(x, τ ) 2π E(x, t -τ ) (23)
where the dielectric permittivity of the polarizable medium is defined as

ε(x, τ ) 2π = δ(τ ) + χ(x, τ ). ( 24
)
In the following we will keep this useful definition even if we do not work in the limit t 0 → -∞. At this stage it is useful to introduce the Fourier transform of the fields defined as

A(ω) = +∞ -∞ dt 2π e iωt A(t) and A(t) = +∞ -∞ dωe -iωt A(ω).
If we temporally suppose that it makes mathematically sense to define a Fourier transform for D(x, t), E(x, t)and P(x, t) then in the t 0 → -∞ limit of Eq. 23 we get

D(x, ω) = P (0) (x, ω) + ε(x, ω) E(x, ω) ε(x, ω) = 1 + 2π χ(x, ω) (25)
Actually it is not at all obvious that such a Fourier transform can be defined unambiguously. Fluctuating and stationary quantum fields are not in general converging sufficiently fast in the future or past directions so that the Fourier transform is not in general a well defined mathematical object for real frequency ω. The problem can be more rigorously handled using Laplace's transform [45,46]. We will not consider this problem in this article keeping the complete analysis for a future work.

Here it is sufficient to observe that in optics the Fourier transform

ε(x, ω) = 1 + +∞ 0 dτ χ(x, τ )e iωτ (26)
is an analytical function in the upper part of the complex plane ω = ω ′ + iω ′′ , i.e., ω ′′ > 0, provided χ(x, τ ) is finite for any time τ ≥ 0. From this naturally follows the symmetry ε(x, -ω) * = ε(x, ω * ). From these analytical properties it is possible to derive the general Kramers-Kronig relations existing between the real part Re[ ε(x, ω)] ≡ ε ′ (x, ω) and the imaginary part Imag[ ε(x, ω)] ≡ ε ′′ (x, ω) of the permittivity. Remarkably Eq. 21 fully satisfies these conditions. Indeed, we can write the Fourier transform of Eq. 21 in the upper complex plane

χ(x, ω ′ + iω ′′ ) = +∞ 0 du π σ u (x)/π u 2 -(ω ′ + iω ′′ ) 2 = +∞ -∞ du 2π σ u (x) πu 1 u -ω ′ -iω ′′ ( 27
)
where we used the relation:

1 u 2 -(ω ′ + iω ′′ ) 2 = 1 2u [ 1 u -ω ′ -iω ′′ + 1 u + ω ′ + iω ′′ ] ( 28
)
as well as the symmetry (definition): σ u (x) = σ -u (x). We thus get in the limit

ω ′′ → 0 + ε(x, ω ′ + i0 + ) = 1 + +∞ -∞ du π σ u (x) u 1 u -ω ′ -i0 + = 1 + P +∞ -∞ du π σ u (x) u 1 u -ω ′ + i σ ω ′ ω ′ . ( 29
)
where P [...] denotes Cauchy's principal value. This is Kramers-Kronig relation if we use the identity

ε ′′ (x, ω ′ ) = σ ω ′ (x) ω ′
along the real axis. We have thus ε ′′ (x, ω) = -ε ′′ (x, -ω) and ε ′ (x, ω) = ε ′ (x, -ω) in agreement with the symmetry requirement (from Eq. 27 we deduce χ(x, -ω) * = χ(x, ω * )). Eq. 21 characterizing the generalization of Hutner-Barnett model [39] is thus a complete representation of a causal linear dielectric medium including both dispersion and dissipation.

III. THE FORMAL QUANTIZATION PROCEDURE

A. Quantizing the matter field equations

In order to obtain a quantized theory of the matter field it is useful to introduce the auxiliary fields f * ω (x, t), f ω (x, t) which correspond to rising-lowering operators:

2ω f ω (x, t) = Z (-) ω (x, t) -2iω , 2ω f * ω (x, t) = Z (+) ω (x, t) 2iω ( 30
)
with

X ω (x, t) = 2ω (f ω (x, t) + f * ω (x, t))
. From this we can obtain the following representation for the Hamiltonian H M :

H M = d 3 x +∞ 0 dω ωf * ω (x, t)f ω (x, t) (31)
Canonical quantization starts with the replacement of the vector fields f ω by operators acting on the Hilbert space associated with the quantum system under study. In particular, we have the replacement f * ω (x, t) → f † ω (x, t). The equal-time commutators between the conjugate canonical variables read:

[X ω (x, t), Π X ω ′ (x ′ , t)] = i δ(ω -ω ′ )δ 3 (x -x ′ )I (32) and [X ω (x, t), X ω ′ (x ′ , t)] = [Π Xω (x, t), Π X ω ′ (x ′ , t)] = 0(33) with I = x ⊗ x + ŷ ⊗ ŷ +ẑ ⊗ẑ the unit dyad. We used the definition [A(x), B(x ′ )] = µ,ν [A µ (x), B ν (x ′ )]x µ ⊗ xν .
This implies the commutation rules:

[f ω (x, t), f † ω ′ (x ′ , t)] = δ(ω -ω ′ )δ 3 (x -x ′ )I. ( 34
)
and

[f ω (x, t), f ω ′ (x ′ , t)] = [f † ω (x, t), f † ω ′ (x ′ , t)]
= 0 allowing a clear interpretation of f ω (x, t) and f † ω (x, t) as lowering and rising operators for the bosonic states associated with the matter oscillators. The quantized Hamiltonian operator is obtained by using the normal-ordered product H M →: H M : such as

: H M := d 3 x +∞ 0 dω ωf † ω (x, t)f ω (x, t) (35)
which allows the elimination of the unphysical vacuum infinite energy. At that stage it is useful to focus on the polarization density P (0) (x, t) which can be equivalently written as

P (0) (x, t) = +∞ 0 dω σ ω (x) πω [f (0) ω (x, t) +f †(0) ω (x, t) (36)
with by definition f -t0) . This mathematical expression allows for an unambiguous definition of the Fourier transform P (0) (x, ω) along the real axis which reads:

(0) ω (x, t) = f ω (x, t 0 )e -iω(t
P (0) (x, ω) = +∞ 0 dω ′ σ ω ′ (x) πω ′ [f (0) ω ′ (x, t 0 ) •e iω ′ t0 δ(ω -ω ′ ) + f †(0) ω ′ (x, t 0 )e -iω ′ t0 δ(ω + ω ′ )] (37)
The presence of Dirac distributions is key in the reasoning since it allows us to identify the frequency ω ′ in the integral with the pulsation ω of the Fourier transform. In other words we have

P (0) (x, ω) = σω(x) πω f (0) ω (x, t 0 )e iωt0
for ω > 0 while

P (0) (x, ω) = σ-ω (x) -πω f †(0)
-ω (x, t 0 )e iωt0 for ω < 0. Now, Maxwell's equations allow the definition of the current operator

J (0) (x, t) = ∂ t P (0) (x, t) which im- plies J (0) (x, ω) = -iω P (0) (x, ω). Furthermore we have [f (0) ω (x, t), f †(0) ω ′ (x ′ , t)] = δ(ω -ω ′ )δ 3 (x -x ′ )I, ( 38
)
and we thus deduce (for ω, ω ′ ≥ 0 or ω, ω ′ ≤ 0)

[ J (0) (x, ω), J †(0) (x ′ , ω ′ )] = ω σ ω (x) π δ(ω -ω ′ )δ 3 (x -x ′ )I, (39)
or equivalently [

J (0) (x, ω), J †(0) (x ′ , ω ′ )] = 0| J (0) (x, ω) ⊗ J †(0) (x ′ , ω ′ )|0
(note that Eq. 39 vanishes if ω and ω ′ have different signs). These formulas are in full agreement with the Gruner-Welsch [52][53][54] formalism which defines a quantum version of the fluctuation-dissipation theorem for the dielectric medium. We point out that in the Langevin noise approach the fundamental Hamiltonian is given by the fluctuating term

H (0) M = d 3 x +∞ 0 dω ωf †(0) ω (x, t)f (0) ω (x, t). ( 40
)
It was the aim of the original derivation by Huttner and Barnett [39] (see also ref. [47,51]) to demonstrate that the Hamiltonian H

M is for all practical needs sufficient for any QED calculations in a dielectric medium. The equivalence will not be studied in this article since it requires a specific study.

B. Formally quantizing the Maxwell equation

In order to quantize Maxwell's equations it is here sufficient to solve formally these equations by considering the polarization P(x, t) as an external source. We will give some details on the derivation here since the use of the F potential is not very common in the context of quantum optics. The strategy will be, like we did for the material field, to start from the Heisenberg picture in which priority is given to the field evolution equations but with 'c numbers' replaced by 'q numbers', i.e., operators, which are both time and space dependent. First, from Maxwell's equations we obtain the following second order differential equation

1 c 2 ∂ 2 t F(x, t) -∇ 2 F(x, t) -∇ × P(x, t) = 0 (41)
which can be formally solved by using a modal expansion of the potential into plane waves. For this we write

F(x, t) = α,j q α,j (t)ǫ α,j Φ α (x) (42)
with α a generic label for the wave vector k α , Φ α (x) = e ikα•x / √ V (here we consider as it is usually done the periodical 'Box' Born-von Karman expansion in the rectangular box of volume V ), j = 1 or 2, labels the two transverse polarization states with unit vectors ǫα,1 = k α × ẑ/|k α × ẑ|, and ǫα,2 = kα × ǫα,1 (conventions and more details are given in Appendix B). The method for solving Maxwell's equations is to transform the secondorder differential evolution Eq. 41 into a set of first-order equations in time. For this we use the variables

c 2 ω α c α,j (t) = d dt q α,j (t) -iω α q α,j (t) (43)
We thus obtain a modal expansion for the fields:

F(x, t) = α,j ic 2ω α c α,j (t)ǫ α,j Φ α (x) + cc. D(x, t) = α,j - ω α 2 c α,j (t) kα × ǫα,j Φ α (x) + cc. B(x, t) = α,j ω α 2 c α,j (t)ǫ α,j Φ α (x) + cc.. (44)
Quantization of those fields holds if we impose the commutation relations:

[c α,j (t), c † β,k (t)] = δ α,β δ j,k , [c α,j (t), c β,k (t)] = 0 and [c † α,j (t), c † β,k (t)
] = 0. We can easily deduce several useful commutation relations (see Appendix B) like for example

[B j (x, t), E k (x ′ , t)] = [B j (x, t), D k (x ′ , t)] = ic l ε j,k,l ∂ l δ 3 (x -x ′ ) (45)
This commutator plays an important role in the Langevin's equation [START_REF] Scheel | [END_REF] approach and it is here deduced directly from our canonical formalism. Furthermore, using this representation based on the F potential the Hamiltonian for the pure field

H F = d 3 x(B 2 + D 2 )/2 becomes: : H F (t) := α,j ω α c † α,j (t)c α,j (t) (46)
which has the usual form for free bosons. We conclude this section by commenting on the use of the electric potential F(x, t) instead of the more usual magnetic potential A(x, t). Standard canonical quantization of the electromagnetic field starts from the separation B(x, t) = ∇ × A(x, t) and E(x, t) = -1 c ∂ t A(x, t) -∇V (x, t) where V is the scalar potential. The usual standard Lagrangian density reads

L s = E 2 -B 2 2 + A • J c -ρV + L M ( 47
)
where L M is the same as in Eq. 2. This Lagrangian density allows us to derive the same equation of motion as done before. However, since the canonical momentum Π V associated to V is vanishing we cannot define commutators for those fields and the quantization procedure becomes tricky unless we use the Gupta Bleuer method [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF].

The usual solution to circumvent this difficulty is to work exclusively in the Coulomb gauge∇ • A(x, t) = 0 which allows a clear separation between physical and redundant electromagnetic variables (see Ref. [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF] for a clear analysis of this problem). The field is thus separated into a transverse contribution E ⊥ (x, t) = -1 c ∂ t A(x, t), B(x, t) = ∇×A(x, t), which can be nicely quantized, and into a longitudinal electric field E || (x, t) = -∇V (x, t) = -P || (x, t) which depends on material fields. Without rewriting here the complete analysis this formalism leads to the following second order equation

1 c 2 ∂ 2 t A(x, t) -∇ 2 A(x, t) - J ⊥ (x, t) c = 0 ( 48
)
with J(x, t) = ∂ t P(x, t). By using a modal expansion A(x, t) = α,j x α,j (t)ǫ α,j Φ α (x) and the transformation

-c 2 ω α a α,j (t) = d dt x α,j (t) -iω α x α,j (t) (49)
similar in spirit to Eqs. 44-50 we can rewrite the relevant quantized transverse fields as

A(x, t) = α,j -ic 2ω α a α,j (t)ǫ α,j Φ α (x) + cc. B(x, t) = α,j ω α 2 a α,j (t) kα × ǫα,j Φ α (x) + cc. E ⊥ (x, t) = α,j ω α 2 a α,j (t)ǫ α,j Φ α (x) + cc. ( 50
)
which must be compared with Eq. 44. The most important difference between the formalism using A and the one based on F is that the later use only local and causal electromagnetic properties such as D, B, and P while the former use a separation between transverse and longitudinal fields and currents which not are causal when taken separately [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF]. It is therefore an advantage of our method to eliminate such unphysical separation from the ground.

Quantization can be easily done by imposing the commutation rules [a α,j (t), a † β,k (t)] = δ α,β δ j,k and [a α,j (t), a β,k (t)] = [a † α,j (t), a † β,k (t)] = 0 from which we deduce the same field commutators as the one discussed previously and derived in the appendix (in particular Eq. 45) .

The canonical Lagrangian density

L c = E 2 ⊥ -B 2 2 + A • J ⊥ c -ρV /2 + L M ( 51
)
with ρ = -∇ • P = -∇ • P || (implying d 3 xE 2 || = d 3 xP 2 || = d 3 xρV
) allows us to introduce the canonical Hamiltonian:

H c (t) = d 3 x[ E 2 ⊥ + B 2 2 -A • J ⊥ c + E 2 || 2 ] + H c,M , (52)
where the canonical Hamiltonian for the matter field differs from H M , as given in Eq. 12, by the amount

H c,M -H M = A • J ⊥ c .
This results from a different canonical momentum for the matter field Π c,Xω = ∂ t X ω + 2σω (x) π A. The canonical Hamiltonian H c (t) can be also expressed as a functional of the canonical variables and we get

H(t) = d 3 x[ c 2 Π 2 A + ∇ × A 2 2 + P 2 || 2 ] + d 3 x +∞ 0 dω (Π c,Xω -2σω (x) π A) 2 + ω 2 X 2 ω 2 ( 53
)
which should be compared with Eq. 12. Despite these differences the terms proportional to A• J ⊥ c cancel out in Eq. 52 and we get the remarkable result

H c (t) = d 3 x B 2 + E 2 2 + H M (t) = H(t) (54)
with

d 3 x :B 2 +E 2 ⊥ : 2 = α,j ω α a † α,j (
t)a α,j (t). Eq. 54 shows that the two formalisms based on A or F should be equivalent (this is also confirmed by the fact that we obtain the same commutators for the electromagnetic fields in both formalisms). Rigorously the correspondence between these two languages can be done by equating the electromagnetic fields obtained with both methods. We thus obtain a relationship between the x α,j (t) and q α,j (t) variables. This is done easily by remarking that we have

D = E + P = E ⊥ + P ⊥ . As explained in Appendix C we get c α,1 (t) = -a α,2 (t) - 1 √ 2 ω α P α,2 (t) c α,2 (t) = a α,1 (t) + 1 √ 2 ω α P α,1 (t). ( 55
)
with the definition

P α,j (t) = d 3 x +∞ 0 dω 2σ ω (x) π X ω (x, t) •ǫ α,j Φ * α (x) (56)
We emphasize that the transformation between two formalisms based on F and A potentials can be handled differently using the so called Power-Zienau unitary transformation [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF]. In this article we did not use this approach (see however Appendix C) and we instead introduced the dual Lagrangian given in Eq. 1. We will however discuss further in the conclusion the relationship between the different formalisms.

Before leaving this section we emphasize that all the dynamical equations used previously could be equivalently obtained from the Heisenberg equation, which, for an operator A(t), reads as

i d dt A(t) = [A(t), H (t)
] where H(t) is the full Hamiltonian operator. This was checked for the equations used in this work. However, we point out that the Hamiltonian formalism has some intrinsic limitations since it relies on some convergence hypothesis which were briefly mentioned in the Section II. Indeed, since for fields we discuss the problem of radiating system the conservation of the energy in a fixed volume is in general not valid unless we accept some specific boundary conditions at spatial infinity (i.e., in general the Poynting vector flow at infinity does not vanish). At the opposite, we could take the quantized dynamical equations as fundamental postulates of the theory without relying on the Hamilton operator and checking the consistency of the formulas at hand (see also Ref. [83] for a wave-function analysis of the photon dynamics in vacuum). However, both approaches give of course similar results as far as the boundary conditions are taken into considerations.

IV. COMPARISON WITH THE MINIMAL COUPLING AND MULTIPOLAR REPRESENTATION

While the previous analysis was given in order to model macroscopic quantum electrodynamics in dielectric media, it is particularly important to discuss the microscopic physical origin of the model in order to evaluate the hypothesis, limitations and possible generalization of the approach. First, observe as a reminder that the minimal coupling scheme associated with the Lagrangian density Eq. 51 corresponds actually to the continuous limit of the non-relativistic Lagrange function [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF]:

L c = d 3 x E 2 ⊥ -B 2 2 + A • J ⊥ c - ρV 2 + n T n -U (57)
where n T n is the discrete sum of the kinetic energy terms T n = 1 2 m n ( dxn(t) dt ) 2 associated with the point-like particle with individual mass m n and spatial coordinates x n (t) (n = 1, ...N is an integer labeling the particles) [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF].

We also included an interaction potential U (x 1 , ..., x N ) which depends on the N particle coordinates. In this description the symmetrized electric current is

J = 1 2 n e n dx n (t) dt δ 3 (x -x n ) + 1 2 n e n δ 3 (x -x n ) dx n (t) dt (58)
(the symmetrization is necessary for satisfying the commutation relations) while the charge density is ρ = n e n δ 3 (xx n ) with e n the individual charge of each moving individual electrons. For neutral matter this total charge n e n is of course globally neutralized by the static charges of the atomic nuclei. Moreover, in the coulomb Gauge description the fundamental current is the transverse current J ⊥ where the transverse delta function δ ⊥ (xx n ), defined in Eq. B14, replaces δ 3 (xx n ). From Eq. 51 we deduce the direct generalization of the Hamiltonian H c given by Eq. 52

H c (t) = d 3 x[ E 2 ⊥ + B 2 2 -A • J ⊥ c + E 2 || 2 ] + n T n + U (59)
with

T n = (p n -enA(xn,t)/c) 2 2mn
is explicitly written using the particle canonical momentum p n = m n dxn(t) dt + e n A(x n , t)/c associated with x n . Moreover, the full evolution leads straightforwardly to Maxwell's equations [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF] and to the quantized Lorentz force dynamical equation:

m n d 2 x n (t) dt 2 = e n E(x n , t) - ∂ ∂x n U +e n 1 2c dx n (t) dt × B(x n , t) -e n 1 2c B(x n , t) × dx n (t) dt (60)
where the symmetrization is required from the ground. The usual approximation made to model a dielectric medium is to expand the Lorentz force as a power of the relative coordinate ξ n between the n th electron and the associated nuclei (we suppose only one electron per atom). The classical result [START_REF] Novotny | Principles of Nano-Optics[END_REF] implies the standard dipolar approximation leading to the force

M n d 2 X n (t) dt 2 ≃ e n i=3 i=1 ξ n,i ∂ ∂X n E i (X n , t) +e n ∂ ∂t [ξ n × B(X n , t) c ] ( 61
)
where X n is the center-of-mass coordinate of the electron-nucleus system with total mass M n . In many applications, e.g., with harmonic excitation, the time derivative in the second terms average to zero so that only the gradient force survives [START_REF] Novotny | Principles of Nano-Optics[END_REF]. This center-of-mass equation is not exhausting the dynamics of the dielectric system and for optical application the fundamental relation is the internal dynamics which is given by the equation

µ n d 2 ξ n (t) dt 2 ≃ e n E(X n , t) - ∂ ∂ξ n U n (ξ n ) +e n 1 2c dX n (t) dt × B(X n , t) -e n 1 2c B(X n , t) × dX n (t) dt (62)
where we have made the assumption U = n U n (ξ n ) and introduced the reduced mass

µ n = m (e) n m (n)
n /M n of the electron-nucleus pair. Clearly, Eq. 8 is a special case of this internal dynamics corresponding to an harmonic interaction potential U n and to the static condition dXn(t) dt = 0 removing the magnetic Lorentz force. We now point out that while the minimal coupling description is clearly sufficient for many purposes a different but rigorously equivalent way to describe the electromagnetic coupling in dielectric systems is to use the multipolar representation [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF][START_REF] Power | [END_REF][START_REF] Craig | Molecular Quantum Electrodynamics[END_REF]83]. It is obtained by adding a term -d dt [ d 3 x A•P c ] to the canonical Lagrangian given by Eq. 57. Through this canonical transformation and after introducing the density of electric and magnetic polarization P and M (which for neutral systems are connected to J and ρ by J = ∂ t P+c∇×M, and ρ = -∇×P) we obtain the Lagrange function:

L multi. = d 3 x[ D 2 -B 2 -P 2 2 + M • B] + n T n -U . (63)
With this definition we deduce the canonical momenta Π A,multi. = -D c , and p n,multi. = m n dxn(t) dt

+ ∂ ∂ ẋn [ d 3 xM • B]
(which differ from the usual minimal coupling values) and a different Hamiltonian

H multi. = d 3 x[ D 2 + B 2 + P 2 2 -M • B -P • D] + n T n + U + n ẋn • ∂ ∂ ẋn [ d 3 xM • B]. ( 64
)
These equations are sufficient to obtain directly the macroscopic Maxwell equations for any dielectric and magnetic media. In the general case the electric dipole density P is expressed in the classical problem by a line integral as [START_REF] Novotny | Principles of Nano-Optics[END_REF][START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF][START_REF] Craig | Molecular Quantum Electrodynamics[END_REF]] P(x, t) = n e n ξ n (t)

1 0 duδ 3 (x - uξ n -x (N ) n ) with x (N ) n
the coordinate of the n th nuclei of electric charge -e n and ξ n + x

(N ) n := x (e)
n is the coordinate of the electron of electric charge e n . We have

generally ρ = -∇ • P = e n [δ 3 (x -x (e) n ) -δ 3 (x -x (N ) n )] (65)
in order to preserve the total charge cancellation of the dielectric medium. We obtain similarly M(x, t) = n en c

1 0 duδ 3 (x -uξ n -x (N ) n )ξ n (t) × [u dξ n dt + dx (N )
n dt ] and thus the classical (non symmetric) current

J = e n [ dx (e) n dt δ 3 (x -x (e) n ) - dx (N ) n dt δ 3 (x -x (N ) n )].( 66
)
We point out that the rigorous extension of the line integral formulas from classical to quantum regime is not a trivial task since we should consider carefully the non commutation of conjugate variables, i.e., positions and velocities [86,87]. It is alternatively possible to use the second quantized formalism involving the wave function operators Ψ(x, t), Ψ † (x, t) associated with particles to remove ambiguities [88]. Hence, at the end we will obtain a symmetrized current like in Eq. 58. With these definitions it is not difficult to obtain

H multi. = d 3 x[ E 2 +B 2 2
]+ n T n + U allowing us to identify the multipolar Hamiltonian with the total energy. For the present purpose we are interested in the crude dipolar approximation

M(x, t) ≈ 0, P(x, t) ≈ n e n ξ n (t)δ 3 (x -X n ) with X n ≃ x (N ) n
the center-of-mass coordinate of the n th electron-nucleus pair (in order to remove the singularity of the delta function δ 3 (x) and the infinite self-interaction coming from the term P 2 /2 in the Lagrangian and Hamiltonian one can alternatively introduce a narrow function ∆(x) [START_REF] Jackson | Classical Electrodynamics 3rd[END_REF]peaked on the origin). Within this approximation the total force acting on the electron-nucleus system leads to the dynamical law:

M n d 2 X n (t) dt 2 ≃ e n i=3 i=1 ξ n,i ∂ ∂X n E i (X n , t) (67)
without the time derivative term in Eq. 61 (Eq. 61 can be justified by relaxing the approximation done on P and M [4, 84], i.e. by introducing a Rontgen current in the Hamiltonian [85]). These forces can be summed over the different particles and leads to the total force acting on the medium: F total = d 3 x i=3 i=1 P i (x, t)∇E i (x, t) (the inclusion of internal forces between particles n and m will not change the final result as far as these additional terms cancel out over the summation and integration). Importantly, in the multipolar framework the internal dynamical motion given by Eq. 62 is kept unchanged, meaning that all optical QED applications in a dielectric medium can be handled equivalently with the minimal coupling or multipolar scheme.

What is however key for the present article is to check the self-consistency of our dual representation discussed in previous sections and based on the electric potential F. To do this, observe that the obvious generalization of Lagrangian Eq. 1 reads:

L = d 3 x[ B 2 -D 2 2 + F • ∇ × P - P 2 2 ] + n T n -U n (68)
As explained in Appendix A this Lagrangian is actually obtained from the multipolar Lagrangian Eq. 63 by using some duality relations between electric and magnetic quantities. By using the electromagnetic dual-Lagrangian L m,multi. (see Eq. A5) with M=0 we have

L = L m,multi. + n T n -U n -d 3 x P 2
2 . The added terms only depend on the particle variables and are not affecting the Maxwell equations. The inclusion of the supplementary potential d 3 x P 2 2 is necessary in order justify the equation of motion Eqs. 62 and 67. As explained before the inclusion of internal forces leads to Eq. 62 and this is true in this representation as well. In this way the Huttner-Barnett Lagrangian L M (see Eq. 3) is only a particular case leading to Eq. 8. Still, in all cases we can easily obtain the macroscopic Maxwell equations as well as the dynamical motion laws for the particles, i.e., Eq. 67 without the time derivative term. Clearly this means that the approach outlined in this paper is equivalent to the minimal coupling or multipolar Lagrangian within the crude dipolar approximation, which is what we wanted to demonstrate (more on this and on the duality relation is given in Appendix A).

V. CONCLUSION

To conclude, we introduced a dual representation of the quantized electromagnetic field in dielectric media based on the transverse electric potential vector F instead of the usual magnetic potential vector A in the Coulomb gauge. The method, contrarily to the usual minimal coupling representation, involves only physical, i.e., local and causal quantities. We showed that our approach is well adapted to the analysis of neutral systems without magnetic property. The dual Lagrangian formalism is equivalent to both standard multipolar and minimal coupling representations within the crude dipolar approximation. This allowed us to re-derive the standard QED equations within the Huttner-Barnett model for quantized polaritons [39,47,51]. The equivalence with the usual commutation relations given in the Langevin noise approach [54,67] was discussed and the Hamiltonian formalism was justified from the ground within our dual approach. In future works we plan to apply the dual Lagrangian formalism to general QED questions in homogeneous and inhomogeneous dielectric media involving coupling with quantum emitters. We believe our work will motivate further studies concerning duality relations and quantization in dense optical media. which we summarize as ǫ-α,j = η j ǫα,j with η 1 = -1 and η 2 = +1. The field being real valued, i.e., F(x, t) = F * (x, t) we deduce q -α,j (t) = η j q * α,j (t) (B3)

In order to solve Eq. 39 it is useful to introduce the auxiliary variables:

Z (±)
α,j (t) = d dt q α,j (t) ± iω α q α,j (t) (B4)

i.e., q α,j (t) = Z with the commutators [c α,j (t), c † β,k (t)] = δ α,β δ j,k and [c α,j (t), c β,k (t)] = [c † α,j (t), c † β,k (t)] = 0. We thus have q α,j (t) = c (2 ω α ) η j c † -α,j (t) -c α,j (t) 2iω α d dt q α,j (t) = c (2 ω α ) η j c † -α,j (t) + c α,j (t) 2 (B11)

This means:

[q α,j (t), q † β,k (t)] = [q † β,j (t), qα,k (t)] = i c 2 δ α,β δ j,k (B12)

and [q α,j (t), q β,k (t)] = [q α,j (t), q † β,k (t)] = 0,

[ qα,j (t), qβ,k (t)] = [ qα,j (t), q † β,k (t)] = 0. From these we deduce [F(x, t), Π F (x ′ , t)] = [F(x, t), Π † F (x ′ , t)] = α,β,j,k

[q α,j (t), q † β,k (t)] c 2 ǫα,j ⊗ ǫβ,k Φ α (x)Φ * β (x ′ ) = i δ ⊥ (xx ′ ).(B13)

with the unit transverse dyadic distribution:

δ ⊥ (x -x ′ ) = α,j
ǫα,j ⊗ ǫα,j Φ * α (x ′ )Φ α (x). (B14)

We have also [F(x, t), F(x ′ , t)] = [Π F (x, t), Π F (x ′ , t)] = 0.

From these relations we deduce the commutation rules: However, X and F conjugate variables commute as well and we have [B(x, t), P(x ′ , t)] = [D(x, t), P(x ′ , t)] = [P(x, t), P(x ′ , t)] = 0. Therefore, from E(x, t) = D(x, t) -P(x, t) we deduce: [E(x, t), E(x ′ , t)] = 0 and [B j (x, t), E k (x ′ , t)] = [B j (x, t), D k (x ′ , t)] = ic l ε j,k,l ∂ l δ 3 (xx ′ ) (B17) Within this approach the electromagnetic fields can also be directly expressed as a function of the rising and lowering operators and we get Eq. 44.

Finally, In the configuration space and by analogy with Eq. B8, it is also possible to define the auxiliary 'photon' fields

Z (±) (x, t) = ∂ t F(x, t) ± ic -∇ 2 F(x, t) (B18)
which obey the following first order equations:

∂ t Z (±) (x, t) = ±ic -∇ 2 Z (±) (x, t) + c 2 ∇ × P(x, t) (B19)
The analogy with Eq. 48 allows us to introduce a photon field Ψ(x, t) such as

Z (-) (x, t) 2c = c 2 (-∇ 2 ) 1/4 Ψ(x, t) (B20)

  the reality requirement Eq. B3 we deduce:Z (±) -α,j (t) = η j Z ǫα,j Φ α (x) + cc. Φ α (x) + cc.. (B7)In this representation the field equations readŻ(±) α,j (t) = ±iω α Z (±) α,j (t) + S α,j (t) (B8)with the source termS α,j (t) = c 2 d 3 x∇ × P(x, t) • ǫα,j Φ * α (x) (B9)One strategy for quantizing the electromagnetic field is to introduce rizing c † α,j (t) and lowering c α,j (t) photon operators defined by Z

  [B(x, t), B(x ′ , t)] = [D(x, t), D(x ′ , t)] = 0 (B15) and [B j (x, t), D k (x ′ , t)] = ic l ε j,k,l ∂ l δ 3 (xx ′ ) (B16)

VI. ACKNOWLEDGMENTS

This work was supported by Agence Nationale de la Recherche (ANR), France, through the SINPHONIE (ANR-12-NANO-0019) and PLACORE (ANR-13-BS10-0007) grants. The author gratefully acknowledge stimulating discussions with G. Bachelier, S. Huant, and C. Genet.

Appendix A: Duality in dielectric and magnetic media

It is interesting to understand that the motivation for the Lagrangian given by Eqs. 1 and 68 is connected to a particular form of duality existing for Maxwell's equations. Consider indeed the most general set of Maxwell's equations for continuous media:

where the electric current and charge density are given by J e = ∂ t P+c∇×M, and ρ e = -∇•P. By introducing H = B -M and D = E + P Eqs. A1 transform as:

with the magnetic current and charge density given by: J m = ∂ t M-c∇×M, and ρ m = -∇•M. There is clearly a duality relation between Eqs. A1 and A2 and we go from the first to the second by the dual replacement:

where by definition H = 1 c ∂ t F + ∇V ′ which depends on the dual potentials F and V ′ . Now, Eqs. A1 can be derived from the standard electromagnetic canonical Lagrangian density

.e., Eq. 47 without the material part. This means that the dual set of equations Eqs. A2 can equivalently be derived from the 'standard' Lagrange function:

The multipolar analog of Eq. 63

and the Hamiltonian becomes

This duality is actually reminiscent of the early work made in the 19 th century when Coulomb and Biot proposed an interpretation of magnetism in term of magnetic charge (in analogy with electrostatic ) while Ampere proposed to give a electric current origin to magnetism [START_REF] Darrigol | From Ampère to Einstein[END_REF]. The equivalence is however not complete since there is no magnetic monopole in the standard electromagnetism approach. More precisely if we consider the total force acting on a medium characterized by P and M in the standard Maxwell representation Eqs. A1 we get from Noether's theorem:

while the dual representation Eqs. A2 leads to the total force

Clearly in general F e,total differs from F m,total unless the time derivatives cancel. This special case occurs in the magnetostatic limit but also for oscillating motions when time average removes the time derivative terms as in Eq. 67 (this is the case for fluctuating forces considered in nano-photonics [START_REF] Novotny | Principles of Nano-Optics[END_REF])

Appendix B: Plane waves modal expansion and quantization of the electromagnetic field

In order to quantize the electromagnetic field using a plane-wave modal expansion we first remind several mathematical properties of the Born von Karman expansion method. First, we have the normalization (obtained from

Second, we have the following symmetries for

The commutators read now

Appendix C: Transformation matrix between A and F representations

We have using Eq. B2:

α,j ẋα,j (t) c ǫα,j Φ α (x)

Similarly for the magnetic field

From Eqs. C1 and C2 we thus deduce the equivalence relations:

with P α,j (t) given by Eq. 56. Equivalently, we can write Eq. C3 in terms of rising and lowering operators and deduce Eq. 55. We can also rewrite Eq. [START_REF] Cohen-Tannoudji | Atom-Photon Interactions: Basic Processes and Applications[END_REF] This transformation can be used in different ways for defining equivalent representations of the electromagnetic field as discussed in Refs. [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF][START_REF] Ackerhalt | [END_REF].


Laurent Schuller
email: laurent.schuller@insa-lyon.fr

Ieee Student Member

Romain Delpoux
email: romain.delpoux@insa-lyon.fr

Jean-Yves Gauthier
email: jean-yves.gauthier@insa-lyon.fr

Xavier Brun
email: §xavier.brun@insa-lyon.fr

Online Estimation and Compensation of Back-Electromotive Forces for Synchronous Reluctance Machine

Keywords: Synchronous Reluctance machine, Back-Electromotive Forces, ReMa model, disturbances rejection

The synchronous reluctance machine is a promising technology in the electrical machine market. They are now considered in various applications such as transportation or industries. Despite the lack of magnet in such machines, a small residual flux may remains in the ferromagnetic material. This residual flux produces back-electromotive forces as the machine rotates. In this paper, two methods are proposed to determine the back-electromotive forces of the machine online. One is based on the measurements of the short-circuit currents, the other utilizes a linear disturbances observer. The estimations of the back-electromotive forces are then used in the control of the machine in order to compensate their undesirable effects such as current ripples. Both techniques require only the knowledge of the speed and classical electrical parameters of the machine. Experimental results are shown to prove the effectiveness of the proposed approaches.

I. INTRODUCTION

Synchronous Reluctance (SynREL) machines constitute a former technology when it comes to electrical actuators [START_REF] Kostko | Polyphase reaction synchronous motors[END_REF]. Because of phase switching and machine design issues, Syn-REL machines only became a serious competitor in the electrical machines market very recently. These past decades, technological advances in power electronic and computer science have made SynREL machines more effective and novel designs are continuously developed to improve them [START_REF] Bao | A novel concept of ribless synchronous reluctance motor for enhanced torque capability[END_REF]. SynREL machines are generally three-phase machines, they only differ from others by the particular shape of their rotor, while their stator are highly similar to those of Permanent Magnet Synchronous Machine (PMSM) and Induction Machines (IM) [START_REF] Donaghy-Spargo | Synchronous reluctance motor technology: opportunities, challenges and future direction[END_REF]. Although IM have a relatively low efficiency, they are currently the most popular technology in the industrial three-phase electrical machines. However, SynREL machines challenge them [START_REF] Boglietti | Induction and synchronous reluctance motors comparison[END_REF] both in terms of efficiency and weight [START_REF] Rassolkin | Efficiency map comparison of induction and synchronous reluctance motors[END_REF]. Compared to the PMSM, SynREL machine's performances are definitely inferior. Yet, the excellent performances of PMSM are mainly due to the rare-earth permanent magnet [START_REF] Rahman | History of interior permanent magnet motors [history[END_REF] used in such technology. In today's society, theses elements raise several issues in ecological and economical fields [START_REF] Jahns | Getting rare-earth magnets out of ev traction machines: A review of the many approaches being pursued to minimize or eliminate rare-earth magnets from future ev drivetrains[END_REF], [START_REF] Boldea | Automotive electric propulsion systems with reduced or no permanent magnets: An overview[END_REF]. The absence of rare-earth magnet in SynREL machines is a significant advantage. In addition, it allows safer operations of the SynREL machines, in particular if a short circuit appears in the stator winding [START_REF] Qin | Reluctance motor control for fault-tolerant capability[END_REF].

Despite the absence of magnet, small back-electromotive forces (EMF) are produced by the residual magnetism of the iron which constitutes the machine [START_REF] Wang | Investigation of self-excited synchronous reluctance generators[END_REF]. In generator application, these residual magnetic flux, that depend on the magnetic history of the machine, are discussed in [START_REF] Ibrahim | The loss of self-excitation capability in stand-alone synchronous reluctance generators[END_REF] [START_REF] Maroufian | Self-excitation criteria of the synchronous reluctance generator in stand-alone mode of operation[END_REF]. Considering the self start of an isolated SynREL generator, the back-EMF have a paramount importance. Indeed, at low current the back-EMF effect is predominant.

The study of the back-EMF was the subject of a previous work by the authors [START_REF] Schuller | Dynamical Model of Residual Magnetism for Synchronous Reluctance Machine Control[END_REF] where a dynamical model of the back-EMF was proposed. However due to the possible changes during the machine operation, it seems useful to propose online solutions to estimate the phenomenon which, to the knowledge of the authors, was never studied.

This article aims to develop two different methods that estimate the back-EMF in order to compensate their unwanted effects: while the first method is based on short-circuit currents measurements, the second one uses a linear observer so as to determine the disturbances caused by the back-EMF. Both methods require only commonly known machine's parameters. Another ambition of this paper is to confront the two proposed methods to the reference ReMa (Residual Magnetism) model presented in [START_REF] Schuller | Dynamical Model of Residual Magnetism for Synchronous Reluctance Machine Control[END_REF].

The paper is structured as follows: SynREL technology and dynamical models considering residual magnetism are presented in section II; in section III, the two new methods to determine back-EMF online are detailed and the compensation approach is presented; in section IV experimental results are used to compare both methods to a reference model for small currents, the reduction of the current ripples by compensating the back-EMF is also highlighted. Finally, the conclusion is presented in section V.

II. MODEL OF THE SYNREL MACHINE

The electrical equations of SynREL machine can be expressed by applying Kirchhoff's voltage law [START_REF] Lyshevski | Synchronous reluctance motors: nonlinear analysis and control[END_REF] as [START_REF] Kostko | Polyphase reaction synchronous motors[END_REF] :

dL(θ e )i abc dt = v abc -R s i abc -e abc , (1)
where i abc = i a i b i c is the vector of stator currents, v abc = v a v b v c is the vector of stator voltages, e abc = e a e b e c is the vector of back-EMF on each phase of the machine, θ e is the electrical position of the rotor, R s (scalar) is the resistance per phase and L(θ e ) is the inductance matrix defined as: 1

L(θ e ) =   L a M ab M ca M ab L b M bc M ca M bc L c   , (2)
with :

               L a = L 0 + L 2 cos(2θ e ), L b = L 0 + L 2 cos(2θ e + 2π 3 ), L c = L 0 + L 2 cos(2θ e -2π 3 ), M ab = M 0 + M 2 cos(2θ e -2π
3 ),

M bc = M 0 + M 2 cos(2θ e ), M ca = M 0 + M 2 cos(2θ e + 2π 3 ), (3)
here, L 0 , M 0 and L 2 , M 2 are respectively average values and amplitudes of phase inductances and mutual inductances. Leakage inductance have been neglected.

Remark 1. In this work, back-EMF denote voltages produced when no current flows through the stator windings. Sometimes this term is employed to denominate the voltages generated by the reluctance variations. Depending on currents in each phase, they are not strictly back-EMF from our point of view.

The back-EMF are generally neglected [START_REF] Lyshevski | Synchronous reluctance motors: nonlinear analysis and control[END_REF] because of the absence of magnet within the rotor, however, as for most electrical actuators, SynREL machines are mainly made of ferromagnetic materials. The hysteresis property of such materials leads to a small magnetic excitation after being subjected to an intense magnetic field. This results in back-EMF not being zero as modeled in [START_REF] Schuller | Dynamical Model of Residual Magnetism for Synchronous Reluctance Machine Control[END_REF]. It results in back-EMF composed of two frequencies in the conventional threephase frame. The amplitudes and the shapes of the back-EMF depend on the configuration during the magnetizing process, their expressions are given in ( 4)

               e a = -Φ rot ω e sin(θ e + δ 0 ) -3I stat ω e M 2 sin(2θ e -σ 0 ), e b = -Φ rot ω e sin(θ e + δ 0 -2π/3) -3I stat ω e M 2 sin(2θ e -σ 0 -2π/3), e c = -Φ rot ω e sin(θ e + δ 0 + 2π/3) -3I stat ω e M 2 sin(2θ e -σ 0 + 2π/3), (4)
where Φ rot and I stat represent the amplitudes of the two frequencies, δ 0 and σ 0 are the direction of the current and the position of the rotor during the magnetization process. Fig. 1 presents an example of simulation results. Unlike IM, the phase-to-phase voltages are not sinusoidal. Indeed, the presence of the second harmonic, due to the residual magnetism within the stator and the reluctance variations, on the phase-to-neutral voltages implies the same second harmonic on the phase-to-phase voltages. 1 Ignoring the homopolar behavior allows us to consider the hypothesis Fig. 2 shows a SynREL machine with two pole pairs (n pp = 2). The dark areas represent the iron while white parts are made of non-magnetic material such as air. In such machines, the variations of inductances are mainly responsible for the torque generation. The expression of the produced torque is deduced from a power balance (iron losses are neglected) and is expressed as follows:

M 0 = -1 2 L 0 and L 2 = M 2 [15]
Γ em = 1 2 n pp ∂L(θe) ∂θe i abc i abc + n pp e abc i abc ω e , (5)
where ω e is the electrical speed of the machine. The electrical position and speed are linked to their mechanical counterparts (θ m and Ω m respectively) by the number of pole pairs: θ e = n pp θ m and ω e = n pp Ω m . Equation [START_REF] Rassolkin | Efficiency map comparison of induction and synchronous reluctance motors[END_REF] shows that the back-EMF have an impact on the torque generated by the SynREL machine. In particular they can cause undesirable torque ripples.

Because the Park reference frame [START_REF] Park | Two-reaction theory of synchronous machines generalized method of analysis-part i[END_REF] is very often used to control electric machines, the equivalent dynamical model that takes into account the back-EMF in that frame is given in [START_REF] Rahman | History of interior permanent magnet motors [history[END_REF]. In the electrical reference frame, the direction of the rotor that involves the maximum of iron is called direction d whereas the so-called q direction is the one with large airgap (see Fig. 2). With these definitions, the d and q axes are unique and orthogonal in the electrical reference frame. The reference of the electrical position θ e is chosen such that θ e = 0 rad coincides with the alignment between the direction d and the phase A windings.

L dq di dq dt = v dq -R s i dq + ω e J L dq i dq -e dq , (6)
with J = 0 1 -1 0 . In [START_REF] Rahman | History of interior permanent magnet motors [history[END_REF], L dq = L d 0 0 L q is the inductance matrix in the Park frame, v dq , i dq and e dq are respectively the voltages, currents and back-EMF vectors in the Park reference frame. These two back-EMF are composed of a sinusoid p x2 with an offset p x0 such that e d = -p d0 -p d2 and e q = -p q0 -p q2 with:

               p d0 = ω e 3 2 Φ rot sin(δ 0 ), p d2 = 3 3 2 I stat ω e M 2 sin(θ e -σ 0 ), p q0 = -ω e 3 2 Φ rot cos(δ 0 ), p q2 = -3 3 2 I stat ω e M 2 cos(θ e -σ 0 ). (7)

III. ONLINE DETERMINATION OF THE BACK-EMF

Because the back-EMF may change during the operation of the SynREL machine, the ReMa model is quite sensitive. When the SynREL machine is operating, the back-EMF can not be directly measured. This section aims to develop two different methods for the online determination of the SynREL machine's back-EMF. The estimations of the two proposed methods can be used compensate the disturbances caused by the back-EMF.

A. Short-circuit currents measurement combined with a Goertzel algorithm

As mentioned previously, there is no problem to shortcircuit the stator of a SynREL machine. it is even an advantage of this technology compared to the PMSM. During a short-circuit (v abc = 0) the equation ( 1) boils down to:

e abc = - dL(θ e ) dt i sc abc -L(θ e ) di sc abc dt -R s i sc abc (8)
The inductance matrix L(θ e ) is considered known. Because the back-EMF are composed of only two frequencies, it can be demonstrated that the short-circuit currents contain these same frequencies (considering the remark in the previous footnote 1 ). They are hence modeled as follows:

i sc = A 0 cos(θ e + Φ 0 ) + A 2 cos(2θ e + Φ 2 ), (9)
and their derivatives are analytically deduced:

di sc dt = -ω e A 0 sin(θ e + Φ 0 ) -A 2 ω e sin(2θ e + Φ 2 ). ( 10
)
To find the parameters present in (9) a Goertzel algorithm is used.

The Goertzel algorithm has been introduced for the first time in 1958 [START_REF] Goertzel | An algorithm for the evaluation of finite trigonometric series[END_REF]. This algorithm is similar to the discrete Fourier transform applied only for some specific frequencies. Due to its low algorithmic complexity, this method is particularly useful when the number of frequencies is small. For this reason, it has been widely used in fault detection from high frequency injection [START_REF] Briz | Online diagnostics in inverter-fed induction machines using high-frequency signal injection[END_REF]. The operating details of the algorithm can be found in [START_REF] Najafi | A novel current mode controller for a static compensator utilizing goertzel algorithm to mitigate voltage sags[END_REF].

From that method, the short-circuit currents derivatives are analytically estimated, thus there is no noise issues. Finally the back-EMF are deduced using [START_REF] Boldea | Automotive electric propulsion systems with reduced or no permanent magnets: An overview[END_REF].

B. Disturbances observer in the Park reference frame

An other method to estimate the back-EMF of the SynREL machine online is presented here. The method is based on a sinusoidal disturbances observer as in [START_REF] Godbole | Performance analysis of generalized extended state observer in tackling sinusoidal disturbances[END_REF]. The observer is designed in the Park reference frame. In that frame, the back-EMF are only composed of a constant value and a sinusoidal signal for each axis as presented in [START_REF] Rahman | History of interior permanent magnet motors [history[END_REF]. The electrical speed ω e is slowly varying, thus considered constant, so the system is linear. Hence a linear Luenberger observer is designed from the system represented by the following matrix:

ẋ = Ax + Bu, y = Cx, (11)
with:

A =             -Rs L d ω e Lq L d 1 L d 0 0 0 1 L d 0 -ω e L d Lq -Rs Lq 0 0 1 Lq 0 0 1 Lq 0 0 0 1 0 0 0 0 0 0 -ω 2 e 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 -ω 2 e 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0             , (12)
B = 1 L d 0 0 0 0 0 0 0 0 1 Lq 0 0 0 0 0 0 , (13)
C = 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ,

the states of the system are:

x = i d i q p d2 ṗ d2 p q2 ṗ q2 p d0 p q0 , (15)
the inputs are:

u = v d v q , (16)
and the outputs are:

y = Cx = i d i q . ( 17
)
The observed states ( x) follow the dynamical equation:

˙ x = A x + Bu + L obs (y -C x). (18)
The gain observer matrix L obs is chosen so that the matrix (A -L obs C) is Hurwitz. This system can not be observed at zero speed. Moreover, some uncertainties, particularly on L d and L q may have an effect on the stability of the observer. If the speed varies, a gain scheduling approach may be considered.

C. Feedforward

In [START_REF] Scalcon | Decoupled vector control based on disturbance observer applied to the synchronous reluctance motor[END_REF] the observed disturbances are used to compensate the latter via the inputs of the system. Here, a similar approach is considered. The control diagram is represented on Fig. 3 and Fig. 4 for the two proposed methods.

Thanks to the switch "Sw" the perturbations caused by the back-EMF can either be compensated ("Sw" activated) by the back-EMF estimations or not ("Sw" inactivated).

IV. EXPERIMENTAL RESULTS AND DISCUSSIONS

In this section the experimental results are presented. The different methods are analyzed and compared to each other.

A. Description of the test bench

For the experimental tests, a BSR90LE154055FB5 SynREL machine from Bonfigioli is used. Its main characteristics are listed in Table I. A Leroy-Somer 95UMC300HAAAA PMSM with the same speed and power range is attached to the SynREL machine by the shafts. Both machines are fed through three-phase inverters composed of NTHL080N120SC1 SiC MOSFET. Position and speed measurements are provided by an incremental encoder with 5000 points per mechanical revolution. All acquisitions are extracted through a dSpace MicroLabBox Rapid Control Prototyping (RCP) System. The experimental setup is presented in Fig. 5.

B. Experimental back-EMF determination

1) Estimated back-EMF: the PMSM is driven so that the mechanical speed is kept constant. The SynREL machine is in a first time short-circuited and no feedforward compensation is considered ("Sw" inactivated). The short-circuit is accomplished by imposing the same constant duty-cycle on each phase of the inverter-fed machine. Despite the potential deadtime effect, it was experimentally verified that the behavior was the same as with a real short-circuit. During this, the Goertzel algorithm method is applied in order to determine the back-EMF from this first method. The identified parameters are stored in memory in order to be used hereafter. After that, the SynREL machine is controlled to regulate the current at zero. There is still no feedforward compensation ("Sw" inactivated). The Observed back-EMF are computed and another estimation of these latter is then available. On the figure Fig. 6 the two proposed methods results are compared to the ReMa model [START_REF] Schuller | Dynamical Model of Residual Magnetism for Synchronous Reluctance Machine Control[END_REF] on phase A for a constant mechanical speed of 72.2 rad/s. The two procedures provide satisfactory results both in term of shape and amplitude although the observer method is slightly better than the Goertzel algorithm based method. The similarity of the results make us more confident on the estimation of the back-EMF.

2) back-EMF compensation: the feedforward compensation is now considered. The estimated back-EMF are used to compensate the current ripples caused by the back-EMF. Again the PMSM is driven so that the mechanical speed is kept constant. The SynREL machine is controlled to keep its currents at zero. While the observer's estimated back-EMF are computed online, those of Goertzel's algorithm are based on the previous short-circuit process. The currents i d and i q are measured for different compensation methods as shown on Fig. 7 and Fig. 8.

The performances of the estimation methods are deduced from the efficiency of the disturbance compensation. Although visible for the direct current i d , the effect of the back-EMF compensation is particularly notable on i q . This is due to the fact that the disturbances are bigger on the q axis than in the direct one because of the lower value of the inductance L q compared to L d . On Fig. 7 the two proposed methods provide similar results. the methods implying a Goertzel algorithm consumes more computational resources. This algorithmic complexity may be improved deducing three from only one measured using the classical relations the currents.

The effects of the disturbances compensation is even better when the speed increases since the higher the speed, the bigger the disturbances. Furthermore, when the disturbances frequency increases, the natural disturbances rejection of the PI current controllers is less and less effective.

Fig. 8 shows the results for the back-EMF compensation methods when the machine rotates at a mechanical speed of 105 rad/s. On that figure the effect of the compensation is visible for both axes. The compensation from the observer method provides better results. However, the methods based on the Goertzel algorithm has no risk of instability and is easier to tune. Indeed, depending on the required frequency resolution, a simple relation relies sample frequency and the measure duration. While the choice of L obs matrix is no that trivial.

3) Perspectives: it is commonly accepted that the residual magnetism can change during operation. To confirm this, the back-EMF are observed for different current loads while the PMSM makes the SynREL machine rotate at a constant speed. No back-EMF compensation is considered here ("Sw" inactivated). Starting from a regulation of the i d and i q currents at zero. A current step is generated by imposing i d = i q = 1.5 A. After that the currents are once again regulated to zero.

The Fig. 9 shows the back-EMF (in the three-phase frame) on phase A observed by the second method for the described operations of the SynREL machine. From that figure, it can be Fig. 9. Observed back-EMF for different operations seen that the back-EMF change depending on the history of the machine. The ReMa model, that does not change online, is close to the solid line of Fig. 9 as proved on Fig. 6. However, after applying a current step on both d and q axis, the back-EMF change. When currents are flowing through the stator windings, the back-EMF are over estimated. Indeed, during such operation, the inductance L q and more likely L d can saturate. This generates errors in our model, these errors on the inductances are included into the observed perturbations. For such cases, the estimation of the back-EMF can be improved by taking into account the saturation of the inductances. More generally the both proposed methods are sensitive to inductances variations. This latter point is suggested as a further work to improve the present development.

V. CONCLUSION

Two methods to determinate the back-Electromotive forces of a SynREL machine online are presented in this paper. The two suggested strategies require only commonly known machine's parameters. The first proposed method is based on a short-circuit currents measurement. The second approach relies on a linear disturbances observer in the reference

The estimations of the two proposed methods can be used to compensate the disturbances caused by the back-EMF. Both methods are compared to a previously developed ReMa model of the SynREL machine. The experimental results show that both methods are effective to estimated the back-EMF of the SynREL machine. Experiments show good results for the compensations at low currents. Making the compensation more higher currents is a perspective that certainly requires precise inductance model.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Simulation of phase-to-phase back-EMF when both the rotor and stator are magnetized. ωe = 209 rad/s, M 2 = 0.058 H, Φrot = 0.0045 Wb, Istat = 0.0228 A, σ 0 = π/4 rad, δ 0 = -2π/5 rad

Fig. 3 .Fig. 4 .

 34 Fig. 3. Control diagram with the first proposed method

Fig. 5 .

 5 Fig. 5. Experimental setup

Fig. 6 .

 6 Fig. 6. Comparison of the two proposed methods and the ReMa model (Ωm = 72.2 rad/s)

Fig. 7 .Fig. 8 .

 78 Fig. 7. Current measurements for different compensations of the back-EMF (Ωm = 36.6 rad/s)

TABLE I MAIN

 I PARAMETERS OF THE SYNCHRONOUS RELUCTANCE MACHINE

	Parameters	Value	Units
	Base speed (mechanical) Ωm	157	rad/s
	Number of pole pairs npp	2	-
	Rated current In	4.5	A
	Power Pn	1500	W
	Nominal phase voltage Vn	230	V
	Resistance per phase Rs	2.6	Ω
	Unsaturated inductance amplitude L 2	0.078	H
	Unsaturated inductance average value L 0	0.144	H
	Unsaturated Mutual amplitude M 2	0.058	H
	Unsaturated Mutual average value M 0	-0.048	H
	Unsaturated d-axis inductance L d	0.289	H
	Unsaturated q-axis inductance Lq	0.095	H