Nadejda Bouldi

Christian Brouder

Gauge invariance and relativistic effects in photon absorption and scattering by matter

Keywords: numbers: 78.70.Ck X-ray scattering, 78.70.Dm X-ray absorption spectra, 11.15.Bt Perturbation theory, applied to gauge field theories, 31.30.jx Nonrelativistic limits of Dirac-Fock calculations

There is an incompatibility between gauge invariance and the semi-classical time-dependent perturbation theory commonly used to calculate light absorption and scattering cross-sections. There is an additional incompatibility between perturbation theory and the description of the electron dynamics by a semi-relativistic Hamiltonian.

In this paper, the gauge-dependence problem of exact perturbation theory is described, the proposed solutions are reviewed and it is concluded that none of them seems fully satisfactory. The problem is finally solved by using the fully relativistic absorption and scattering cross-sections given by quantum electrodynamics. Then, a new many-body Foldy-Wouthuysen transformation is presented to obtain correct semi-relativistic transition operators. This transformation considerably simplifies the calculation of relativistic corrections. In the process, a new light-matter interaction term emerges, called the spin-position interaction, that contributes significantly to the magnetic x-ray circular dichroism of transition metals.

We compare our result with the ones obtained by using several semi-relativistic time-dependent Hamiltonians. In the case of absorption, the final formula agrees with the result obtained from one of them. However, the correct scattering cross-section is not given by any of the semi-relativistic Hamiltonians.

I. INTRODUCTION

There is a well-known conflict between time-dependent perturbation theory and gauge invariance (see section II for a non exhaustive list of references). Indeed, | φ n |ψ(t) | 2 gives the probability to find the system described by the state |ψ(t) in the eigenstate |φ n at time t, where |ψ(t) is a solution of the time-dependent Schrödinger equation for the Hamiltonian H(t) = H 0 + H 1 (t), while |φ n is an eigenstate of the time-independent Hamiltonian H 0 . A time-dependent gauge transformation of H(t) will be assigned to H 1 (t) but not to H 0 , which must remain independent of time. The transition probability is then modified because the state |ψ(t) is gauge-transformed and not the state |φ n (see section II for a more detailed argument).

As a (not so well-known) consequence, the absorption and scattering cross-sections derived by semi-classical arguments are not gauge invariant. Since a basic principle of quantum physics states that an observable has to be gauge invariant to be physically meaningful, we meet a serious difficulty: "Until this problem is understood, therefore, it seems that no calculation can be trusted at all." 1 There is a further conflict between time-dependent perturbation theory and semi-relativistic physics. 2 The semi-relativistic approximation of |φ n is obtained by applying to it the time-independent Foldy-Wouthuysen transformation: |φ FW n = U H0 |φ n . The semi-relativistic approximation of |ψ(t) is derived from the timedependent Foldy-Woutuysen transformation |ψ FW (t) = U H(t) |ψ(t) .

Since U H(t) = U H0 , the transformed transition probability | φ FW n |ψ FW (t) | 2 is not equal to | φ n |ψ(t) | 2 , even if the Foldy-Wouthuysen transformations U H0 and U H(t) are known to all orders.

In this paper we discuss and solve these two conflicts. In a nutshell, the gauge problem is solved by deriving relativistic absorption and scattering cross-sections from quantum electrodynamics instead of the usual semi-classical argument where the incident light wave is described by a time-dependent potential. The semirelativistic problem is solved by applying a new manybody Foldy-Wouthuysen transformation to the relativistic cross-sections instead of describing the dynamics of the system with a semi-relativistic Hamiltonian. The final result is a semi-relativistic absorption and scattering cross-section involving a new term that couples the spin and the position operators. In a companion paper, we show that this new term contributes significantly to the x-ray absorption of magnetic materials. 3 We now describe the outline of this paper. Section 2 discusses the gauge transformation of transition probabilities and reviews the solutions to the gauge-dependence problem proposed in the literature. Since none of them was widely accepted, we turn to the quantum electrodynamics framework in section 3, where we derive the relativistic electric dipole, quadrupole and magnetic dipole relativistic transition operators. In section 4, we derive a new many-body Foldy-Wouthuysen transformation that we apply to the transition operators. They are used to obtain semi-relativistic absorption and scattering cross-sections in sections 5 and 6, where a new spin-position term is derived. In section 7, the conflict between time-dependent perturbation theory and semi-relativistic methods is described and illustrated by the calculation of spin-position term using four different semi-relativistic Hamiltonians commonly used in the literature. The conclusion presents possible extensions of the present work.

II. GAUGE INVARIANCE

The gauge invariance of the absorption and scattering cross-sections of light is a long-standing problem. It started in 1952 when Willis Lamb calculated the spectrum of Hydrogen in two gauges and obtained different results. 4 This gave rise to a long series of papers up to this day. 1,2, In 1987, the same Lamb (then Nobel prize winner) still considered this as "one of the outstanding problems of modern quantum optics." 37 We quickly describe the meaning of gauge invariance and then consider its failure in semi-classical perturbation theory.

A. The principle of gauge invariance

The two homogeneous Maxwell equations ∇× E + Ḃ = 0 and ∇ • B = 0, where the dot denotes time derivative, imply the local existence of a vector potential A and a scalar potential Φ such that B = ∇ × A and E = -∇Φ -Ȧ. We denote A = (Φ, A). The same E and B are obtained from the potentials A ′ = (Φ-Λ, A+∇Λ), that we also denote A ′ = A -∂Λ, where Λ is any smooth function of space and time. In classical electromagnetism, gauge invariance means that the physics described by A and A ′ is the same.

In quantum mechanics, consider a non-relativistic Hamiltonian

H A = (p -eA) 2 2m + eΦ,
or a relativistic (Dirac) Hamiltonian

H A = cα • (p -eA) + mc 2 β + eΦ,
where α = (α x , α y , α z ) and β are the Dirac matrices. Both Hamiltonians are of the form

H A = f (p-eA)+eΦ,
where f is some function. For such Hamiltonians it can be checked that

M † Λ (i ∂ t -H A ′ )M Λ = i ∂ t -H A , where M Λ = e ieΛ/ . As a consequence, if ψ is a solution of the time-dependent Schrödinger equation i ψ = H A ψ, then ψ ′ = M Λ ψ is a solution of i ψ′ = H A ′ ψ ′ .
In quantum mechanics, a gauge transformation consists in both a change of the potentials and a change in the phase of the wavefunctions. An observable O A depending on the electromagnetic potential A is said to be

gauge invariant if M † Λ O A ′ M Λ = O A for every function Λ(t, r
). An observable must be gauge invariant to be considered a true physical quantity. [START_REF] Cohen-Tannoudji | Photons and Atoms. Introduction to Quantum Electrodynamics[END_REF] The principle of gauge invariance has become a cornerstone of particle physics. Since general relativity can also be considered as a gauge theory [START_REF] Blagojević | Gauge Theories of Gravitation[END_REF] , it may be safely said that gauge invariance was the guiding principle of most of the fundamental physics of the twentieth century. Therefore, we need to check that the cross-section formulas are gauge invariant to ensure their true physical nature.

Note that the time-dependent Dirac or Schrödinger equations are always gauge invariant but the timeindependent ones are not because H A is not gauge invariant due to the the scalar potential Φ. Indeed, under a gauge transformation Φ becomes Φ -Λ and the term Λ cannot be compensated for in the absence of a time derivative.

B. Gauge dependence of transition probabilities

In time-dependent perturbation theory, a system is assumed to be in the ground state |φ g of a timeindependent Hamiltonian H a0 . Then, at time t 0 , an electromagnetic wave represented by the time-dependent potential a is added to the system (with total potential A = a 0 + a), which is represented at time t by the state |ψ(t) . A good way to take both the initial state and the dynamics into account is to use the evolution operator U A (t, t 0 ), which is the solution of

i ∂ t U A (t, t 0 ) = H A (t)U A (t, t 0 ) with the boundary con- dition U A (t 0 , t 0 ) = 1. Thus, |ψ(t) = U A (t, t 0 )|φ g .
The probability of a transition to the eigenstate |φ n of H a0 at time t is

P ng (t) = | φ n |ψ(t) | 2 = | φ n |U A (t, t 0 )|φ g | 2 . (1)
Since we want to ensure that the breakdown of gauge invariance is not due to an approximation, we work with exact (i.e. to all orders) perturbation theory. If we carry out a gauge transformation of the perturbation from a to a ′ = a -∂Λ, then the evolution operator becomes 33

U A ′ (t, t 0 ) = M Λ (t)U A (t, t 0 )M † Λ (t 0 ), (2)
where

A ′ = a 0 + a ′ . Therefore φ n |ψ ′ (t) = φ n |U A ′ (t, t 0 )|φ g is drdr ′ φ * n (r)e ieΛ(r,t) U A (r, t; r ′ , t 0 )e -ieΛ(r ′ ,t0) φ g (r ′ ),
which is generally different from φ n |ψ(t) since Λ(r, t) is an arbitrary function (take for example Λ(r,

t) = r • k(t),
where k(t) is an arbitrary function of time). Therefore, in general, P ′ ng (t) = P ng (t) and the transition probabilities calculated in the two gauges are different. Moreover, since the transition rate entering cross-sections is the derivative of the transition probability with respect to time, [START_REF] Bohm | Quantum Mechanics: Foundations and Applications[END_REF] the arbitrariness of the transition rate is increased by the fact that an arbitrary function Λ(r, t) enters the integrand. Indeed, several papers evaluate the discrepancy between the probability calculated with two different gauges, and they find that it is generally not small. 4,39,51 By properly choosing Λ, the discrepancy can even be made arbitrary large. 49 The absence of gauge invariance is due to the fact that the operator is transformed but not the eigenstates of a 0 . This is called a hybrid transformation in the literature. 13

C. Proposed solutions

The lack of gauge invariance of transition probabilities is an alarming problem to which several solutions have been proposed. Since no clear consensus appears to have emerged, 46 we present a critical review of these solutions.

The first one, called the consistent procedure, was proposed by Forney and coll. and Epstein. 9,11 It is based on the observation that, if instead of gauge-transforming a we transform the potential a 0 of the initial Hamiltonian to get H a ′ 0 , where a ′ 0 = a 0 -∂Λ, then the evolution operator becomes again

U A ′ (t, t 0 ) (because a ′ 0 + a = a 0 + a -∂Λ = A ′ ) but the eigenstates |φ g and |φ n are also transformed into time-dependent states |φ ′ g = M Λ (t 0 )|φ g and |φ ′ n = M Λ (t)|φ n .
Therefore, the transition probability is now conserved. In other words, gauge invariance is lost if we subtract ∂Λ from the perturbation but not if we subtract it from the unperturbed Hamiltonian H 0 = H a0 .

Since the reference states |φ n become time-dependent, we leave the standard framework of time-dependent perturbation theory where the initial Hamiltonian H 0 does not depend on time. Moreover, it is not physically clear why the gauge transformation should be applied to H 0 and not to the perturbation.

In spite of these difficulties, many authors proposed to use the consistent procedure. However, as noticed by Yang, 2 this does not really solve the problem because, if we start the calculation with the initial potential H a0 and the perturbation a ′ , the transition probability is P ′ ng (t). If we then use the consistent procedure to come back to the perturbation a, then we still find P ′ ng (t) and we do not recover the result P ng (t). In other words, the transition probability is now gauge invariant (in the sense that a change of gauge does not modify the result) but it is gauge-dependent (in the sense that the result depends on the gauge we use in the perturbation to start the calculation). This gauge dependence would be a serious problem because we would have to select the "true" physical gauge for the perturbation.

A second solution appeared in a series of papers starting in 1976, 2,7, 15,16,21,23,24,29,33,35,40 , where Yang and collaborators proposed to define a gauge invariant transition probability. His idea is to start from the gauge-invariant (but time-dependent) initial Hamiltonian

H 0 (t) = (p -ea 0 -ea(t)) 2 2m + eV, (3)
where V describes the electron-electron and electronnuclear interactions so that H = H 0 + eφ: the perturbation is only the scalar potential φ. Then, the Hamiltonian H 0 (t) is diagonalized at every time t: H 0 (t)|φ n (t) = E n (t)|φ n (t) and the transitions are calculated between the time-dependent states |φ n (t) . The corresponding transition probabilities are indeed gauge invariant. This solution has been used up to this day, 38,50,53 although it was also strongly criticized. 10,17,18,[25][26][27][30][31][32][START_REF] Feuchtwang | [END_REF]58 The main arguments against Yang's interpretation are: (i) the quantity E n (t) is not physical because you cannot measure an energy at a given time with arbitrary precision;

(ii) the time-dependent states |φ n (t) can be neither prepared nor detected; (iii) the term V in Eq. (3) should be removed from H 0 (t) because it is a scalar potential and, as such, not gauge invariant. But if V is removed, then H 0 (t) is so far from the true Hamiltonian that perturbation theory is no longer valid. Following Goldman, 59 Feuchtwang, Kazes and coll. proposed the following alternative solution. 25,26,31,32,60 They started from the well-known fact that the equations of motion of a Lagrangian are not modified by the addition of the total time derivative of a function. [START_REF] Cohen-Tannoudji | Photons and Atoms. Introduction to Quantum Electrodynamics[END_REF] Thus, two Lagrangians that differ by a total time derivative are equivalent. [START_REF] Weinberg | The Quantum Theory of Fields. I Foundations[END_REF] Then, they remark that the addition of a total time derivative e Λ to the Lagrangian induces a gauge transformation A → A -∂Λ of the Hamiltonian. 47,[START_REF] Cohen-Tannoudji | Photons and Atoms. Introduction to Quantum Electrodynamics[END_REF]60 Finally, they use such a total derivative to compensate for the electric potential that is the cause of the gauge variance of the Hamiltonian. However, it is difficult to distinguish this procedure from picking up a specific gauge, namely the Weyl or temporal gauge where the scalar potential vanishes. We can conclude this short review by stating that no solution was found fully satisfactory.

To determine when gauge invariance can be achieved at the first order of perturbation theory, we consider a Dirac Hamiltonian in two gauges A and A ′ = A-∂Λ and we calculate the difference

ψ|H A -H A ′ |ψ ′ = e ψ|cα • ∇Λ + Λ|ψ ′ .
The advantage of the Dirac Hamiltonian is that the difference H A -H A ′ does not depend on A, but a similar calculation can be carried out in the non-relativistic case. 58 Then, we notice that cα • ∇Λ = (i/ )[H D , Λ] for any Dirac Hamiltonian H D . Thus, if |ψ and |ψ ′ are eigenstates of H D with energy E and E ′ , we obtain

ψ|H A -H A ′ |ψ ′ = e ψ| Λ|ψ ′ + ie E -E ′ ψ|Λ|ψ ′ . ( 4
)
If we consider the absorption cross-section of a photon of energy ω, then energy conservation implies that 40,58 In other words, by restricting the gauge transformations to those satisfying Λ = -iωΛ, the absorption cross-section, calculated up to first order in perturbation theory, is gauge invariant. However, in the resonant scattering cross-section, energy conservation does not apply to the transition involving intermediate states, and the cross-section is not gauge invariant even for those gauges. 8,40,51 Equation (4) shows that the matrix elements are also gauge invariant for a time-independent gauge transformation and energy conserving processes (i.e. E ′ = E). However, the gauge invariance principle is not supposed to restrict to gauges satisfying specific constraints such as Λ = -iωΛ or Λ = 0.

E ′ = E + ω. Thus, if Λ satisfies Λ = -iωΛ, then φ|H A -H A ′ |φ ′ = 0.
This rapid overview shows that, in the published semiclassical approaches where the photon is represented by an external potential, the transition probabilities are not gauge-invariant and no proposed solution has reached general acceptance. Therefore, we turn now to a framework where both electrons and photons are quantized: quantum electrodynamics (QED).

D. Quantum electrodynamics

In QED the incident light is no longer described by an external electromagnetic field but by a photon, i.e. a state in a bosonic Fock space. Therefore, a scattering experiment is now described by the transition from an initial state involving both the electronic system in its ground state and the incident photon, to a final state involving both the electronic system in its (possibly) excited state and the scattered photon. Thus, the energy of the initial and final states is the same and, in the Schrödinger picture, the gauge transformation is expressed in terms of time-independent operators instead of a time-dependent function Λ. [START_REF] Nakanishi | Covariant Operator Formalism of Gauge Theories and Quantum Gravity[END_REF] Equation (4) suggests that transition probabilities, which are now described through the so-called S-matrix, could be gauge invariant. This is indeed the case, although a review of the literature on the gauge invariance of QED might look ambiguous because the kind of gauge transformation considered in different works can vary. In standard textbooks, "the S-matrix is gauge invariant by construction" 63 because only the so-called ξ-term is modified. In the most general gauge transformation, the space of states change from one gauge to the other. [START_REF] Steinmann | Perturbative Quantum Electrodynamics and Axiomatic Field Theory[END_REF] For example, in the Coulomb gauge, only the transverse degrees of freedom are quantized and the photon states form a Hilbert space built by acting on the vacuum with creation operators of left and right polarized photons, while in the Lorenz gauge four degrees of freedom are quantized and the states (built by acting on the vacuum with creation operators of the left, right, longitudinal and scalar photons) can have a negative norm. In the Lorenz gauge, the Lorenz condition cannot be satisfied as an operator equation, [START_REF] Strocchi | An Introduction to Non-Perturbative Foundations of Quantum Field Theory[END_REF] it becomes a subsidiary condition used to determine a subspace of physical states with positive norm.

In other words, the state spaces of the Coulomb and Lorenz gauges have a quite different nature and the relation between them is delicate. Haller managed to show that the usual gauges are equivalent by devising a common framework containing all of them. 41 Note also that the gauge-invariance can only be expected for the renormalized S-matrix. [START_REF] Bialynicki-Birula | [END_REF][67][68] The gauge invariance under a general infinitesimal gauge transformation is well established within the Becchi-Rouet-Stora-Tyutin (BRST) approach: matrix elements of gauge-invariant operators between physical states are independent of the choice of the gauge-fixing functional if and only if the physical states |α satisfy Q|α = 0, where Q is the BRST charge. [START_REF] Weinberg | The Quantum Theory of Fields II[END_REF]70 The case of finite BRST transformations is in progress. 71,72 To summarize the discussion, the gauge invariance of the renormalized S-matrix is established for infinitesimal gauge transformations and for a reasonably large classe of gauges. 13,41,68,[73][74][75][76][77][78] In other words, it is proved at the physicist level of rigour.

The most studied gauges are the Lorenz and Coulomb gauges. Renormalization is perfectly established for the Lorenz gauge, but in most practical calculations the subsidiary condition (Gauss' law) is not enforced. 79 Although it was proved that the S-matrix elements are often the same with and without the subsidiary condition, 10,13,26,80,81 this fails when the Hamiltonian is suddenly changed, 82 as in the sudden creation of a core hole in photoemission or x-ray absorption 83,84 . In that case, Gauss' law has to be imposed in the Lorentz gauge and the Coulomb gauge result is recovered. 82 We choose to use quantum electrodynamics in the Coulomb gauge because it is the most accurate gauge for low-energy many-body calculations. [START_REF] Cohen-Tannoudji | Photons and Atoms. Introduction to Quantum Electrodynamics[END_REF][START_REF] Grant | Relativistic Quantum Theory of Atoms and Molecules[END_REF]

III. RELATIVISTIC MATRIX ELEMENTS

Since we now have a gauge-invariant framework, we can calculate the relativistic matrix elements that will be used in x-ray scattering and absorption cross-sections.

A. The Hamiltonian

The quantum field Hamiltonian describing the interaction of light with matter in the Coulomb gauge is: 13,[START_REF] Cohen-Tannoudji | Photons and Atoms. Introduction to Quantum Electrodynamics[END_REF][START_REF] Greiner | Field Quantization[END_REF][START_REF] Bialynicki-Birula | Quantum Electrodynamics and Quantum Optics[END_REF]

H = H e + H γ + H eγ ,
where

H e = drψ † (r) cα • (-i ∇ -ea) + βmc 2 + eφ ψ(r) + drdr ′ ρ(r)ρ(r ′ ) 8πǫ 0 |r -r ′ | ,
where φ is a time-independent scalar external potential (for instance the nuclear potential), a is a timeindependent vector potential (describing an external magnetic field) and ψ are fermion field operators. Normal ordering is implicit in H e . It is the QED form of the Dirac Hamiltonian in the Coulomb gauge. The manybody version of this Hamiltonian is

H N = N n=1 cα n • (-i ∇ n -ea(r n )) + β n mc 2 + eφ(r n ) + m =n e 2 8πǫ 0 1 |r m -r n | ,
where α n and β n act on the nth Dirac electron. It can be given a well-defined mathematical meaning if the electronic system is described with respect to the Dirac sea, 88 although the physical validity of the Dirac sea is sometimes disputed. 89 The photon Hamiltonian is

H γ = ǫ 0 2 dr|E ⊥ | 2 + c 2 |B| 2 = k,l ω k,l a † k,l a k,l ,
where l stands for the polarization of a mode (there are two independent directions for a given wavevector k) and

H eγ = -ec drψ † (r)α • A(r)ψ(r),
describes the photon-matter interaction in the Coulomb gauge. According to Bialynicki-Birula, the Hamiltonian H also describes the dynamics of gauge-invariant states in any gauge. [START_REF] Bialynicki-Birula | Quantum Electrodynamics and Quantum Optics[END_REF] The many-body version of this interaction Hamiltonian is

H I = -ec N n=1 α n • A(r n ).

B. S-matrix elements

Since we saw that the S-matrix is gauge invariant, we calculate its matrix-elements. We recall that

S = lim ǫ→0 T (e -i ∞ -∞ Hǫ(t)dt ), (5)
where H ǫ (t) = e -ǫ|t| e iH0t H eγ e -iH0t . The adiabatic switching factor e -ǫ|t| enables us to describe physical processes as matrix elements of S between eigenstates of H 0 = H e + H γ . The limit ǫ → 0 can be shown to exist up to technical assumptions. 90 Note that H 0 is not quadratic because of the Coulomb interaction term in H e . The eigenstates of H e are correlated multi-electronic wavefunctions. As a consequence, we are not in the textbook framework, the time-dependence of H ǫ (t) cannot be calculated explicitly and the Feynman diagram technique is no longer available to describe electrons.

We can bypass this problem with the so-called "noncovariant" approach, [START_REF] Heitler | The Quantum Theory of Radiation[END_REF] The sign of the damping term γ has been the object of some controversy. [93][94][95][96] Let us stress again that, since H e is not quadratic, we essentially work in the Schrödinger picture, where the operators are independent of time, instead of the standard interaction picture which is used in most textbooks. Both approaches are equivalent. 97 A modern version of the Schrödinger picture of QFT is given by Hatfield. [START_REF] Hatfield | Quantum Field Theory of Point Particles and Strings[END_REF] Our purpose is now to calculate the matrix elements m|H eγ |n , where H eγ is independent of time. The second quantized expression for the photon field in the Schrödinger picture is: 99

A(r) = k,l 2ǫ 0 V ω k ǫ k,l a k,l e ik•r + ǫ ⋆ k,l a † k,l e -ik•r .
Note that we do not assume the polarization vectors ǫ k,l to be real.

We denote |n = a † k,l |0 |Ψ n an eigenstate of H 0 where one photon is present in mode k, l and the electrons are in state |Ψ n with energy E n . The energy of |n is e n = ω k,l + E n . The interaction Hamiltonian H eγ is linear in A which is linear in photon creation and annihilation operators so that only one-photon transitions are possible. The state |n can make transitions towards

|a = |0 |Ψ m by absorption and |e = a † k,l a † k ′ ,l ′ |0 |Ψ m by emission. From now on, we denote ω = ω k,l , ǫ = ǫ k,l , ω ′ = ω k ′ ,l ′ and ǫ ′ = ǫ k ′ ,l ′ .
The corresponding matrix elements are:

a|H eγ |n = -ec 2ǫ 0 V ω ǫ • Ψ m | ψ † αψe ik•r |Ψ n , and
e|H eγ |n = -ec 2ǫ 0 V ω ′ ǫ ′⋆ • Ψ m | ψ † αψe -ik ′ •r |Ψ n ,
where ψ † αψe ±ik•r = ψ † (r)αψ(r)e ±ik•r dr.

C. Electric dipole and multipole transitions

To carry out a multipole expansion of the previous matrix elements, we shall continue working with quantum fields instead of the usual many-body expressions. In that framework, the expressions are simpler because there is no electron index and we can use the following well-known trick. [START_REF] Strocchi | An Introduction to Non-Perturbative Foundations of Quantum Field Theory[END_REF]100 Let F = ψ † (r)f (r)ψ(r)dr, where f is some function of r. To calculate the commutator of F with some Hamiltonian H 0 , we go to the interaction picture and define F I (t) = e iH0t/ F e -iH0t/ . Then, the timederivative ḞI of F I is given by -i ḞI (t) = [H 0 , F I (t)]. Now, we notice that F is related to the density operator ρ(r) = ψ † (r)ψ(r) by F = ρ(r)f (r)dr. Thus,

-i ḞI (t) = -i ρ(r, t)f (r)dr = [H 0 , F I (t)]
. If H 0 conserves the electric charge, the continuity equation e ρ(r) = -∇ • j holds, where j is the electric current operator. By taking t = 0 to recover the operators in the Schrödinger picture, we obtain

[H 0 , F ] = i e ∇ • j(r)f (r)dr = - i e j(r) • ∇f (r)dr = -i c ψ † (r)αψ(r) • ∇f (r)dr. ( 7
)
To find the electric dipole transition term we apply Eq. ( 7) with f (r) = ǫ • r and

H 0 = H e : [H e , ψ † (r)ǫ • rψ(r)dr] = -i c ψ † (r)αψ(r) • ǫdr,
and we obtain in the dipole approximation e ik•r ≃ 1

a|H eγ |n = e(E m -E n ) i 2ǫ 0 V ω Ψ m | ψ † ǫ • rψ|Ψ n .
To deal with electric quadrupole and magnetic dipole transitions, we expand to the first order: e ik•r ≃ 1 + ik•r. We apply Eq. (7) with

f (r) = ǫ • rk • r and H 0 = H e : [H e , ψ † ǫ • rk • rψ] = -i cψ † αψ • (ǫk • r + kǫ • r),
where we removed the integral sign for notational convenience. Thus,

ψ † ǫ • αk • rψ = i c [H e , ψ † ǫ • rk • rψ] -ψ † ǫ • rk • αψ. If we add ψ † ǫ • αk • rψ to both terms we obtain 2ψ † ǫ • αk • rψ = i c [H e , ψ † ǫ • rk • rψ] -ψ † (ǫ × k) • (r × α)ψ.
Finally, up to electric quadrupole transitions

a|H eγ |n = e∆E i 2ǫ 0 V ω Ψ m | ψ † T ψ|Ψ n , (8)
where ∆E = E m -E n and

T = ǫ • r + i 2 ǫ • rk • r - c 2∆E (ǫ × k) • (r × α). (9)
The first term of T is the usual electric-dipole operator, the second one is the electric-quadrupole operator and the third one will turn out to be the magnetic-dipole operator (see section IV D). Similarly,

e|H eγ |n = e∆E i 2ǫ 0 V ω ′ Ψ m | ψ † T ′ ψ|Ψ n , (10)
where

T ′ = ǫ ′⋆ • r - i 2 ǫ ′⋆ • rk ′ • r + c 2∆E (ǫ ′⋆ × k ′ ) • (r × α).

IV. SEMI-RELATIVISTIC REPRESENTATION

In the previous sections, we have shown that gauge invariance is ensured by describing the interaction of light and matter with quantum electrodynamics, where photons are quantized and electrons are described by fourcomponent Dirac spinor quantum fields.

However, in most solid-state calculations, we do not use Dirac spinors but two-component (Pauli) wavefunctions. Moreover, semi-relativistic expressions are often physically clearer. Therefore, we need to link the two representations by using a generalization of the Foldy-Wouthuysen transformation.

In this section, we first describe the Foldy-Wouthuysen transformation and its new many-body extension. Then, we use this framework to calculate the relativistic corrections to the dipole and quadrupole transitions. The calculations are considerably simpler than the usual approach, where the relativistic corrections are derived from a semi-relativistic Hamiltonian.

A. The Foldy-Wouthuysen transformation

The idea of the Foldy-Wouthuysen transformation is the following. If H D is a time-independent relativistic Hamiltonian, it has the form

H D = H 0 + H 11 H 12 H 21 H 22 ,
where H 0 = mc 2 β and each H ij is a 2x2 matrix. We write H D as the sum of even and odd parts This transformation does not satisfy Eriksen's condition U = βU † β discussed in the Appendix. This is because the product U = e iS (2) e iS (1) does not satisfy this equation even if e iS (1) and e iS (2) do. Silenko recently derived the correction that must be applied to go from Foldy-Wouthuysen to Eriksen transformations, 103 and he showed that the correction is at an order beyond the one we consider in this paper.

H D = H 0 + E + O, where E = H 11 0 0 H 22 , O = 0 H 12 H 21

B. Many-body Foldy-Wouthuysen transformation

To generalize the Foldy-Wouthuysen approach to the many-body Dirac Hamiltonian we face the following problem. The generalization of H 0 is imposed by the many-body Dirac Hamiltonian:

H 0 N = N n=1 β n mc 2 ,
where β n is the matrix β acting on the nth electron (i.e.

β n = 1 ⊗(n-1) ⊗ β ⊗ 1 ⊗(N -n)
). This definition is valid because H 0 N commutes with the projector P N onto the space of antisymmetric N -body states.

We show in the Appendix that a Foldy-Wouthuysen transformation can be defined whenever we have a selfadjoint operator η (with η 2 = 1) to define parity. In the one-body case, β 2 = 1 and η = β defines parity. But in the many-body case the operator N n=1 β n suggested by H 0 N cannot be used for that purpose because its square is not proportional to the identity (it contains products

β n β m ). It turns out that η = β 1 ⊗ • • • ⊗ β n
is the natural many-body generalization of β. Indeed, η † = η and η 2 = 1. Moreover, η commutes with P N , which allows us to work with tensor products instead of antisymmetric tensor products.

In the literature, the Foldy-Wouthuysen transformation was studied for two-body Hamiltonians, 104-106 , but the results were rather complicated and not easy to extend to the many-body case.

The even and odd parts of H N are then H 0 N + E and O, respectively:

E = e N n=1 φ 0 (r n ) + e m =n V (r m -r n ), O = N n=1 cα n • π n = N n=1 O n ,
where V (r) = e 8πǫ0|r| is the Coulomb potential and

π n = -i ∇ n -ea 0 (r n ).
At first order in c -1 , the Foldy-Wouthuysen operator is U = e iS (1) where

S (1) = - i 2mc 2 n β n O n .
Indeed, it can be checked that i[S (1) , H 0 N ] = -O removes the odd term of H D . At this order U = U 1 ⊗ • • •⊗ U N is a tensor power of one-body Foldy-Wouthuysen operators, as proposed by Moshinksy and Nikitin. 107 However, this tensor-power form does not hold at higher orders. Indeed, we show now that at the next order, the many-body Foldy-Wouthuysen Hamiltonian is the sum of one-body and two-body contributions. The usual formal Foldy-Wouthuysen transformation U = e iS (1) e iS (2) can be carried out almost unchanged and we find, with m as expansion parameter, at order m -2 :

H FW = H 0 N + E + 1 2mc 2 N n=1 β n O 2 n - 1 8m 2 c 4 N n=1 [O n , [O n , eϕ n + eV ]] + 1 8m 2 c 4 p =n β p β n [O p , [O n , V ]].
This Hamiltonian obeys ηH FW η = H FW which makes it a Foldy-Wouthuysen Hamiltonian.

It rewrites

H FW = N n=1 H n FW + H MB FW . (11)
where each H n FW is the usual one-body Foldy-Wouthuysen Hamiltonian:

H n FW = β n mc 2 + eφ 0 (r n ) + p =n eV (r n -r p ) + 1 2m β n π 2 n -e Σ n • b 0 (r n ) - 2 e 8m 2 c 2 ∇ • E n + e 8m 2 c 2 Σ n • (π n × E n -E n × π n ) where b 0 (r n ) = ∇ × a 0 (r n ) and E n = -∇φ n (r p ) - p =n ∇V (r n -r p ).
The mass-velocity term βn 8m 3 c 2 (p n • p n ) 2 would be obtained by expanding to higher order.

The new two-body term H MB FW arises because V (r mr n ) = V (r mn ) is a two body operator:

H n,p FW = e 8m 2 c 2 N p =n ∆V (r np ) -Σ n • (π n × ∇V (r np ) -∇V (r np ) × π n ) + 2 β n β p (α n • ∇ n )(α p • ∇ p )V (r np ) .
By using: 108

∂ j ∂ k V (r) = e 2 8πǫ 0 -δ jk 4π 3 δ(r) -δ jk 1 r 3 + 3r j r k r 5 ,
the derivatives in the last term can be rewritten

(α n • ∇ n )(α p • ∇ p )V (r np ) = jk α j n α k p ∂ j ∂ k V (r np ) = e 2 8πǫ 0 - 4π 3 α n • α p δ(r np ) - α n • α p |r np | 3 + 3 α n • r np α p • r np |r np | 5 .
This expression looks superficially like some contributions to the Breit interaction as presented by Bethe and Salpeter. 109 However, they are different since the Breit interaction is due to the exchange of a photon and not to a semi-relativistic effect. Note that the last two terms are singular. It is known that the expansion of the Foldy-Wouthuysen transformation as a power serie in 1/c 2 becomes more and more singular because of the presence of the Coulomb potential. 110 At order m -2 , the transformation writes

U = 1 + 1 2mc 2 n β n O n - 1 8m 2 c 4 n β n O n 2 + 1 4m 2 c 4 n β n m β m O m , E
and it obeys U = ηU † η. We also checked that U 2 is odd in H D after paying attention to the discontinuity at zero discussed in the Appendix. Thus, the positive (negative) energy eigenstate of H D are transformed into even (odd) states by the action of U .

e iS T e -iS = T + i[S, T ] + ∞ n=2 i n L n (T ) n! ,
where L(T ) = [S, T ] and

L n (T ) = L(L n-1 (T )). If U = U 1 ⊗ • • • ⊗ U N ,
where U i = e iSi , we can calculate the action of U on each variable independently. Removing temporarily the constant -i/2mc 2 , we take the one-body operator S = βO and compute

L(ǫ • r) = c[βα • (p -ea 0 ), ǫ • r] = c ij βα i ǫ j [p i , r j ] = -i c ij βα i ǫ j δ ij = -i cβα • ǫ,
and

L 2 (ǫ • r) = -i c 2 [βα • (p -ea 0 ), βα • ǫ] = -i c 2 ij (p i -ea 0i )ǫ j [βα i , βα j ] = i c 2 ij (p i -ea 0i )ǫ j [α i , α j ],
where we used βα i = -α i β and β 2 = 1. We compute

[α i , α j ] = 2i k ǫ ijk σ k 0 0 σ k = 2i k ǫ ijk Σ k ,
which defines Σ k the components of Σ. Therefore,

L 2 (ǫ • r) = -2 c 2 (p -ea 0 ) • (ǫ × Σ).
So that, for each particle, and up to O(m -2 ),

U n ǫ • r n U † n = ǫ • r n -i 2mc β n α n • ǫ -4m 2 c 2 π n • (ǫ × Σ n ).
The many-body version is obtained by summing the right-hand side over n.

In the matrix elements D = φ|U ǫ • RU † |ψ , recall that |ψ = η|ψ and |φ = η|φ because |Ψ and |Φ are positive energy states, as shown in the Appendix. Therefore, φ|U ǫ • RU † |ψ = φ|ηU ǫ • RU † η|ψ and all the terms that are odd in U ǫ • RU † are eliminated by the matrix elements. This eliminates the term proportional to β n α n and we are left with

D = N n=1 φ|ǫ • r n -4m 2 c 2 π n • (ǫ × Σ n )|ψ .

D. Semi-relativistic multipole transitions

From Eq. ( 9), we write the multipole transitions

M = i 2 M 1 - c 2∆E M 2 ,
where

M 1 = n φ|U ǫ • r n k • r n U † |ψ , M 2 = n φ|U (ǫ × k) • (r n × α n )U † |ψ ,
correspond to the electric quadrupole and magnetic dipole transitions, respectively. Since multipole transitions are smaller than dipole ones, it is enough to use the first two terms of the Baker-Campbell-Hausdorf formula.

The term [S n , ǫ • r n k n • r n ] is odd and disappears in the matrix element. Thus, at the order we consider,

M 1 = n φ|ǫ • r n k • r n |ψ . Let T 2 = (ǫ × k) • (r × α). We write [βO, T 2 ] = c[βα • p, T 2 ] -ec[βα • a 0 , T 2 ] = cβ({α • p, T 2 } -e{α • a 0 , T 2 }).
The anticommutators are

{α • p, T 2 } = ijkl ǫ jkl (ǫ × k) j (α i α l p i r k + r k p i α l α i ) = 2(ǫ × k) • ( Σ + L),
and

{α • a 0 , T 2 } = ijkl ǫ ikl (ǫ × k) i a j r k {α j , α l } = 2(ǫ × k) • (r × a 0 ).
Note that Σ = gS with g = 2 (because the spin operator is S = Σ/2). Thus, we recover the fact that the Dirac equation gives a gyromagnetic factor g = 2 to the electron. Moreover, L + Σ = L + 2S is the total magnetic moment of the electron. Finally, since r n × α n is odd,

M 2 = n β n mc φ|(ǫ × k) • ( Σ n + Λ n )|ψ ,
where

Λ n = L n -er n × a 0 (r n ) is
M = n φ| i 2 ǫ • r n k • r n - β n 2m∆E (ǫ × k) • ( Σ n + Λ n )|ψ . V. ABSORPTION CROSS-SECTION
The absorption cross section is calculated by assuming that initially the system of electrons is in state |I that can be transformed into Foldy-Wouthuysen eigenstate |i , with energy E i , and that a photon k, ǫ is present. In the final state there is no photon and the system is in state |F (|f after transformation).

The transition probability per unit time from state m to state n is related to the T-matrix elements by: 114

w = 2 δ mn Im m|T |m + 2π δ(e n -e m )| n|T |m | 2 . ( 12
)
and must be divided by c/V (rate at which the photon crosses a unit of surface) to obtain the cross section.

Since we consider real transitions (i.e. m = n), only the second term is present. From ( 8) and using the result of transformation derived in the previous section:

σ = 4π 2 α 0 ω f | f |T FW |i | 2 δ(E f -E i -ω),
where T FW is:

T FW = n ǫ • r n + i 2 ǫ • r n k • r n -4m 2 c 2 π n • (ǫ × Σ n ) - β n 2mω (ǫ × k) • ( Σ n + Λ n ),
with α 0 the fine structure constant and

∆E = E f -E i = ω.
It corresponds to the usual formula for the cross section 115 with two more terms: the third one and the last one.

The third term was already found by Christos Gougoussis in his PhD thesis, 116 but his final result was not in agreement with ours because of his use of the commutation relation, as described in section VII F. We rewrite it by using π = (m/i )[r, H FW 0 ] + O(c -2 ), where H FW 0 is the Foldy-Wouthuysen Hamiltonian, to get:

-4m 2 c 2 f |π • (ǫ × Σ)|i = i 4mc 2 (E i -E f ) f |r • (ǫ × Σ)|i = i ω 4mc 2 f |(ǫ × r) • Σ|i .
We call spin-position operator the operator (ǫ×r)•Σ. Its evaluation at the K-edge of materials will be presented in a companion paper. 3 The amplitude of the last term depends on the choice of the space origin in the Coulomb gauge for a 0 . It does not make the cross section gauge dependent because the states are changed accordingly when choosing the origin of the gauge. If the origin of the gauge is chosen at the atom position, fields larger than 10 6 T are required for this term to be significant. Such fields are way beyond laboratory accessible values.

VI. SCATTERING CROSS-SECTION

The scattering cross section is calculated by assuming that initially the system of electrons is in state |I with a photon k i , ǫ i and that in the final state the system is in state |F with a scattered photon k f , ǫ f . We do not consider the special case when k i , ǫ i = k f , ǫ f . Eqs. ( 8), ( 10) and ( 12) yield:

w = 2π F δ(E f + ω f -E i -ω i ) L e 2 c 2 2ǫ 0 V 1 √ ω i ω f F |e -k f ψ † α • ǫ ⋆ f ψ|L L|e ki ψ † α • ǫ i ψ|I E i -E l + ω i + iγ + F |e ki ψ † α • ǫ i ψ|L L|e -k f ψ † α • ǫ ⋆ f ψ|I E i -E l -ω f 2 ,
where γ > 0 and

e k ψ † α • ǫψ = 3 j=1 e ik•r ψ † (r)α j ψ(r)ǫ j dr.
The scattering cross-section is related to w by: 117

d 2 σ dΩdω f = V 2 (2π) 3 ω 2 f 1 c 4 w.
Since the electric charge is related to the classical electron radius r e by e 2 = 4πǫ 0 r e mc 2 , we obtain the relativistic Kramers-Heisenberg scattering cross-section:

d 2 σ dΩdω f = (r e mc 2 ) 2 ω f ω i F δ(E f + ω f -E i -ω i ) L F |e -k f ψ † α • ǫ ⋆ f ψ|L L|e ki ψ † α • ǫ i ψ|I E i -E l + ω i + iγ + F |e ki ψ † α • ǫ i ψ|L L|e -k f ψ † α • ǫ ⋆ f ψ|I E i -E l -ω f 2 .
In this expression, the sum over |L involves a complete set of states, with positive and negative energies. Since E i is usually the positive energy of the ground state including the electron rest energy, we have

E i = mc 2 +E ′ i > 0, where E ′ i is the usual (negative) ground state energy. If |L is a positive energy state, we have E l = mc 2 +E ′ l with E ′ l > E ′ i and the first term is resonant at ω i = E ′ l -E ′ i . If |L is a negative energy state, then E l = -mc 2 -E ′ l and E i -E l -ω f = 2mc 2 + E ′ i -E ′ l -ω f cannot be resonant in standard experimental conditions.
We show that the resonant scattering term has a semirelativistic expansion close to, but different from, the standard one. 118 If we are interested in the resonant part of the scattering cross section, then E l > 0 and

d 2 σ dΩdω f = ( r e m 2 ) 2 ω f ω i f δ(E f + ω f -E i -ω i ) L > (E l -E i )(E f -E l ) f |T ′f l FW (ǫ f )|l l|T li FW (ǫ i )|i E i -E l + ω i + iγ 2 . with T ij FW (ǫ i ) = n ǫ i • r n + i 2 ǫ i • r n k i • r n -4m 2 c 2 π n • (ǫ i × Σ n ) - β n 2m∆E ij (ǫ i × k i ) • ( Σ n + Λ n ),
and

T ′ij FW (ǫ f ) = n ǫ ⋆ f • r n - i 2 ǫ ⋆ f • r n k f • r n -4m 2 c 2 π n • (ǫ ⋆ f × Σ n ) + β n 2m∆E ij (ǫ ⋆ f × k f ) • ( Σ n + Λ n ), where ∆E ij = E i -E j .
As in the absorption case, the spin-position term in the transition operator is not present in the usual formula. 118

VII. OTHER METHODS

In this section, we compare our semi-relativistic transition matrix elements with the ones obtained by using time-dependent perturbation theory where the timeevolution is described by several time-dependent semirelativistic Hamiltonians: the one proposed by Blume, the "gauge-invariant" Foldy-Wouthuysen one, the textbook Foldy-Wouthuysen one and the effective Hamiltonian derived in non-relativistic QED (NRQED). Before making this comparison, we first explain why using a time-dependent semi-relativistic Hamiltonian in a perturbation calculation can lead to incorrect results.

A. Foldy-Wouthuysen subtelties

In this section, we assume that the exact timedependent Foldy-Wouthuysen operator U is known. Thus, the following difficulties are not related to the use of an approximation, but to the interplay of the Foldy-Wouthuysen method with perturbation theory.

The first subtelty was noticed by Nieto: 59,119 If |Ψ is a solution of the time-dependent Dirac equation (i ∂ t -H)|Ψ = 0, then the Foldy-Wouthuysen transformation turns it into |ψ = U |Ψ , where U is a unitary timedependent operator. The time-dependent Dirac equation for |Ψ implies that |ψ is a solution of the timedependent Schrödinger equation (i ∂ t -H ′ )|ψ = 0, where H ′ = U HU -1 +i (∂ t U )U -1 is the time-dependent Foldy-Wouthuysen Hamiltonian. In the following, an uppercase Greek letter (|Φ or |Ψ ) refers to a solution of the Dirac equation and the corresponding lowercase letter (|φ or |ψ ) to its Foldy-Wouthuysen transformation.

As a consequence, a matrix element Φ|H|Ψ is not equal to φ|H ′ |ψ , but to φ|H ′ -i (∂ t U )U -1 |ψ . In other words, H ′ has to be used to calculate the states |φ and |ψ but not to calculate the matrix elements of the Hamiltonian.

The second subtelty was observed by Yang. 2 and concerns the most straightforward way to use the Foldy-Wouthuysen Hamiltonian H ′ (t), where the time dependence is now explicit, to compute transition probabilities. This Hamiltonian is split into a time-independent part H ′ 0 and a time-dependent one

H ′ 1 (t), so that H ′ (t) = H ′ 0 + H ′ 1 (t). The scalar product φ ′ n |ψ(t) , where |φ ′ n is an eigenstate of H ′ 0 , cannot be equal to the relativis- tic scalar product Φ n |Ψ(t) . Indeed |ψ(t) = U (H 0 + H 1 (t))|Ψ(t) but |φ ′ n = U (H 0 + H 1 (t))|Φ n because |φ ′ n
and |Φ n are independent of time whereas U (H 0 + H 1 (t)) depends on time. Since only the QED relativistic matrix elements where found to be gauge invariant, φ ′ n |ψ(t) is generally not physically meaningful.

The two problems combine if first-order perturbation theory is naively applied with Foldy-Wouthuysen eigenstates and Hamiltonian. The Foldy-Wouthuysen interaction Hamiltonian

H ′ 1 (t) = H ′ (t) -H ′ 0 = U (H 0 )(H(t) - H 0 )U † (H 0 ). As a consequence, φ ′ n |H ′ 1 (t)|φ ′ g is not equal to Φ n |H 1 (t)|Φ g .
To illustrate the variety of results that can be obtained by using first-order perturbation theory with semirelativistic Hamiltonians, we now examine four Hamiltonians used in practice. To help comparing these Hamiltonians, we express them in a common one-particle framework.

B. The Blume Hamiltonian

Blume discussed the interaction of light with magnetic matter by starting from the Hamiltonian: 118,120

H B = π 2 2m + eV - e 2m σ • B - e 4m 2 c 2 σ • (E × π), (13)
where π = p -eA. This Hamiltonian is the sum of four terms: (i) the kinetic energy of the electron, (ii) an external potential, (iii) the Zeeman interaction between the electron and a magnetic field and (iv) the spin-orbit interaction (because, for a spherical V and a static A, σ

• (E × p) = -1 r dV dr σ • (r × p) = -1 r dV dr σ • L).
There are several differences between our notation and Blume's: he considers a many-body Hamiltonian (involving sums over electrons) and writes ij V (r ij ) for our eV , he adds the Hamiltonian H γ of the free photons, he uses A/c, ∇ × A/c and s where we use A, B and σ/2, finally, his Zeeman term is wrong by a factor of 2 in his first two papers on the subject, 118,120 but this was corrected in the third one. 121 In this third paper, Blume also replaces E by -Ȧ. This is not compatible with his quantized description of the photon field. Indeed, the time-derivative Ȧ is present in the Lagrangian but, after the Legendre transformation leading to the Hamiltonian, Ȧ is replaced by its canonical momentum -E. Note that Blume does not sketch any derivation of his Hamiltonian.

C. Foldy-Wouthuysen Hamiltonian

We consider now the so-called "gauge-invariant" Foldy-Wouthuysen Hamiltonian for positive-energy states up to order 1/(mc) 2 : 122

H FW = H B + mc 2 - e 2 8m 2 c 2 ∇ • E - ie 2 8m 2 c 2 σ • (∇ × E).
The difference between the Foldy-Wouthuysen and the Blume Hamiltonians consists of three terms: the rest energy mc 2 of positive-energy eigenstates, the Darwin term proportional to ∇ • E and a last term, proportional to σ • (∇ × E) and called the curl-term, that we discuss presently. A basic difference between H B and H FW must first be stressed: the former is a QED expression where the quantum fields A, B and E are independent of time because they are written in the Schrödinger representation, while the latter was derived under the assumption that A and V are external time-dependent potentials. In particular, the curl-term disappears if the external field A is independent of time. [START_REF] Strange | Relativistic Quantum Mechanics[END_REF] In the semi-classical treatment of light-matter interaction, the photons are represented by an external time-dependent potential and this term is present. These Hamiltonians can be written H(A, Φ), where the total vector potential A and scalar potential Φ are a sum A = a 0 + a, Φ = φ 0 + φ, of static external potentials a 0 and φ 0 (representing the static internal and external fields) perturbed by dynamical potentials a and φ representing the incident electromagnetic wave. We write the interaction Hamiltonian as H I = H(A, Φ) -H(a 0 , φ 0 ). The two Hamiltonians H B and H FW lead to two different interactions:

H B I = h 1 + h 2 + h 3 + h 4 + h 5 + h 6 and
H FW I = H B I + h 7 ,
where

h 1 = e 2 2m a 2 , h 2 = - e m a • π 0 , h 3 = - e 2m σ • (∇ × a), h 4 = e 2 4m 2 c 2 σ • (e × a)
,

h 5 = - e 4m 2 c 2 σ • (e × π 0 ), h 6 = e 2 4m 2 c 2 σ • (e 0 × a), h 7 = - ie 2 8m 2 c 2 σ • (∇ × e),
with π 0 = p -ea 0 . The curl-term in H FW is the origin of the presence of h 7 in H FW

I

. It originates from the term i (∂ t U )U -1 in the time-dependent Foldy-Wouthuysen Hamiltonian. The Darwin term gives no contribution to the interaction because ∇ • e is zero for the electromagnetic wave. The terms h 5 and h 6 were omitted by Blume, who considered them to be small. 120 We shall see that h 5 is the source of a spin-position term which is not negligible in x-ray magnetic circular dichroism (XMCD) spectra. 3

D. Textbook Foldy-Wouthuysen Hamiltonian

Standard textbooks often derive a Foldy-Wouthuysen Hamiltonian H TFW which is the same as H FW , except for the fact that π is replaced by p in the spinorbit term. 123,124 A mass-velocity term -(p • p) 2 /8m 3 c 2 is often added 123 but its contribution to the radiationmatter interaction is zero. The difference with the Foldy-Wouthuysen Hamiltonians is a term in σ • E × A. This results in the absence of h 4 and h 6 in the perturbation Hamiltonian, which changes the transition probabilities.

E. NRQED

To deal with QED calculations involving bound states, Caswell and Lepage proposed an alternative approach to relativistic effects, called non-relativistic QED (NRQED), which turned out to be highly successful. 125 They wrote the most general gauge-invariant nonrelativistic Lagrangian terms and fitted the coefficients of these terms to known QED processes. 126 The corresponding NRQED Hamiltonian is the same as H FW up to order c -2 , but its interpretation is different. 126 Indeed, NRQED is a quantum field theory, and the fields are independent of time in the Schrödinger representation. However, the curl-term is present in time-independent NRQED although it is generated by a time-dependence in H FW . In particular, the curl-term must not be removed from the Hamiltonian to calculate matrix elements of the Hamiltonian operator, in contrast to the example of section VII A.

Besides these four different Hamiltonians, we consider an additional source of discrepancies between authors: the commutators.

F. Commutators

To derive the multipole expansion of the matrix element of H I , it is useful to replace π by a commutator with H 0 = H(a 0 , φ 0 ). The derivations that start from Blume's interaction Hamiltonian usually use the relation: 127,128

p = mi [H 0 , r]. (14)
However, if one considers the static Hamiltonian given by Blume (13), its commutator with r is:

[H B 0 , r] = - i m π 0 + e 4m 2 c 2 (i )(σ × e 0 ),
which is different from Eq. ( 14) because p is replaced by π 0 = p -ea 0 and because of the term proportional to c -2 . The commutator of r with H TFW 0 and H FW 0 are the same. In H FW I and H B I , when π 0 in h 2 is rewritten as a function of the commutator, the extra relativistic term leads to the cancellation of h 6 , which is important in XMCD. On the other hand, it leads to a contribution

e 2 4m 2 c 2 σ • [∇v 0 × a] in H TFW I .
If the mass-velocity term -(p • p) 2 is present in H 0 , the additional contribution to the commutator, i (p•p)p 2m 3 c 2 is small compared to i m p if the order of magnitude of the kinetic energy of the core state satisfies E k << mc 2 .

For all the Hamiltonians presented here, using the relation [p, v 0 ] = i ∇v 0 , the electric field in matter writes at zeroth order in c -2 as a function of the commutator of π 0 with H 0 :

e 0 = -∇v 0 = -i e [π 0 , H 0 ].
In the case of absorption, the commutator transforms into a factor ∆E = -ω in the cross section so that h 5 and h 6 lead to the same contribution to the matrix element:

-ie ω 4m 2 c 2 σ • (a × π 0 ),
which corresponds to the spin-position interaction. Explicit calculations showed that this contribution can appear two times, one time or cancel completely, according to which Hamiltonian and which commutator was used.

Starting from H FW , the same absorption cross section as in our new approach can be derived. However, in the case of scattering, even with H FW , there is a factor ∆E/ ω which is not correct. The same kind of discrepancy was already observed in the literature. 8,40 VIII. CONCLUSION

This paper was written because of the gaugedependence of transition probabilities in the semiclassical approach and because we observed, after other authors, 127,129 that different semi-relativistic Hamiltonians lead to different cross-sections.

Our solution makes essential use of quantum electrodynamics as the correct gauge-invariant framework to discuss the interaction of light with matter. It is wellknown that the semi-classical and QED absorption crosssections are identical in the Coulomb gauge. 130 This is compatible with our discussion because, to go from the Coulomb gauge to another gauge, the semi-classical approach only involves the operator M Λ , while QED involves a redefinition of the space of states, including in an essential way non-physical polarizations and even ghost states in the BRST approach. This redefinition is able to maintain gauge invariance where the semi-classical M Λ fails to do so.

In the present paper, the stationary states of the electronic system was taken to be eigenstates of H e . The interaction Hamiltonian H eγ can modify these states through various QED effects, for example the Breit interaction discussed by Bethe and Salpeter. 109 We expect these contributions to be small in x-ray spectroscopy.

The explicit calculation of the spin-position contribution at the K-edge of Fe, Co and Ni will presented in a forthcoming publication. 3 It was known since Heisenberg in 1928, 131 that the Thomson cross-section which is due to the A 2 term in the non-relativistic approach, can be derived from the relativistic framework by using a sum over negative-energy states. [START_REF] Strange | Relativistic Quantum Mechanics[END_REF]132 We intend to provide a more accurate discussion of the contribution of negative-energy states to the scattering cross-section by using our many-body Foldy-Wouthuysen approach.

To derive a many-body Foldy-Wouthuysen transformation, we first notice that, in the one-body case, β endows the space of spinors with the structure of a Krein space, where β is then called a fundamental symmetry. 133 For quite a different purpose, 134 we investigated the tensor product of such spaces and showed that the fundamentaly symmetry of the N th tensor power is essentially η = β ⊗N . The abstract Krein-space framework leads us naturally to the following theorem:

Assume that H D and η are self-adjoint operators and η 2 = 1. Then, there is a unitary operator U such that

U = ηU † η and ηU H D U † η = U H D U † . Moreover, if |ψ D is an eigenstate of H D with positive (resp. nega- tive) eigenvalue, then |ψ = U |ψ D satisfies |ψ = η|ψ (resp. |ψ = -η|ψ ).
The condition U = ηU † η does not appear in Foldy and Wouthuysen works. It was added by Eriksen. 103,106,135 It means that U is self-adjoint for the Krein-space structure.

Let us start with general considerations involving a self-adjoint operator η such that η 2 = 1. It can be used to define projectors

B ± = (1 ± η)/2. It is clear that B + + B -= 1, B 2 ± = B ± , B † ± = B ± and B + B -= B -B + = 0. A vector |ψ is said to be even (odd) if η|ψ = |ψ (η|ψ = -|ψ ).
Then, any vector |ψ can be written as the sum of its even part B + |ψ and its odd part B -|ψ . An operator H is said to be even (odd) if it transforms an even state into an even (odd) state and an odd state into an odd (even) state. An operator H is even (odd) if and only if ηHη = H (ηHη = -H). Thus, the theorem states that U H D U is an even operator. Any operator H can be written as the sum of its even part B + HB + +B -HB -and its odd part B + HB -+B -HB + .

Our proof of the theorem is essentially a generalized and rigorous version of Eriksen's proof. 135 We use the fact that H D is self-adjoint to define λ = signH D by functional calculus. The operator λ is called the flat band Hamiltonian in topological insulator theory. 136 In physical terms, let |ψ D be an eigenstate of H D for the energy E, then λ|ψ D = |ψ D if E ≥ 0 and λ|ψ D = -|ψ D if E < 0. Since η and λ are self-adjoint and η 2 = λ 2 = 1, they are bounded and ηλ is unitary: ηλ(ηλ) † = ηλλη = η 2 = 1 and (ηλ) † ηλ = 1. By the spectral theorem for unitary operators, 137 there is a unique family of orthogonal projections P t such that ηλ = π -π e it P t dt.

In the finite dimensional case we could write this 138 We now show that this U satisfies the intertwining relation ηU = U λ. Indeed, the relation U2 = ηλ implies U = U † ηλ. By multiplying from the left with η and using ηU † η = U we find ηU = U λ. This important relation implies that H = ηHη and that |ψ = U |ψ D is even if |ψ D is a positive energy state and odd if |ψ D is a negative energy state.

The first property is easy to show:

ηHη = ηU H D U † η = U λH D λU † = U H D λ 2 U † = H,
because λ commutes with H D since it is a function of H D .

To show the second property, let Γ ± = (1 ± λ)/2, so that Γ + projects onto the space of positive energy and Γ -of negative energy, and recall that B ± = (1 ± η)/2. For a one-body system, B ± projects onto the large/small components. Then, U Γ ± = U/2 ± U λ/2 = U/2 ± ηU/2 = B ± U , which can be used to show that the so that the small components of |ψ are zero for a positive-energy state. This is not true for many-body systems. For example, if we neglect antisymmetrization for notational convenience, a two-body state can be obtained as the tensor product of one-body wavefunctions:

|ψ D = φ 1 ψ 1 ⊗ φ 2 ψ 2 .
Then, the even part of |ψ D is

φ 1 0 ⊗ φ 2 0 + 0 ψ 1 ⊗ 0 ψ 2 ,
while its odd part is

φ 1 0 ⊗ 0 ψ 2 + 0 ψ 1 ⊗ φ 2 0 .
The characterization of U as the square root of ηλ is not easy to handel. We give now a much simpler characterization:

Let U be a unitary operator continuously defined (outside zero) in terms of H D such that: (i) U = ηU † η;

(ii) ηU H D U † η = U H D U † ; (iii) U 2 (-H D ) = -U 2 (H D ).
Then U † ηU = ±sign(H D ).

To prove this, define Z = U † ηU . Clearly, Z † = Z. Moreover, Z is defined in terms H D since U is. However, for Z to be a function of H D in the sense of functional calculus, Z needs to commute with H D : 141 if we multiply condition (ii) from the right by ηU we find ηU H D = U H D U † ηU = U H D Z. Hence,

ZH D = U † ηU H D = U † U H D Z = H D Z.
Thus, there is a real function f (t) and a family of orthogonal projections P t corresponding to the eigenstates of H D such that

Z = ∞ -∞ f (t)dP t .
Moreover, Z 2 = 1 because Z 2 = U † ηU U † ηU = U † η 2 U = U † U = 1. Therefore, f 2 (t) = 1 for every t. Finally, observe that Z = η 2 U † ηU = ηU 2 , and condition (iii) implies that Z is an odd function of H D : f (-t) = -f (t). To conclude that f (t) = ±signt, we need to add the condition of continuity on f outside zero. Indeed, functional calculus is valid for measurable functions and we could build a non-continuous odd function f such that f 2 = 1 outside the origin. In practice this does not take place because U is smoothly defined in terms of H D , except at zero. No odd continuous function can satisfy f 2 = 1 over R. It has to be discontinuous at zero. Since it is crucial that f 2 = 1 everywhere, we can choose either sign0 = 1 or sign0 = -1. Both solutions are valid.

  0 , satisfy βEβ = E and βOβ = -O. Note that H 0 is also even. If |ψ D is a solution of the Dirac equation H D |ψ D = E|ψ D , where H D is the Dirac Hamiltonian, then the upper two components of |ψ D are called the large components and the lower two the small components. The Dirac equation couples the large and small components of |ψ D through the odd terms of H D . Foldy and Wouthuysen 101 looked for a unitary operator U that decouples the large and small components of |ψ = U |ψ D . In other words, H = U H D U † has only even components: H = βHβ. The method proposed by Foldy and Wouthuysen consists in successive transformations of the form U = e iS . 101,102 .

C

  . Semi-relativistic dipole transitions Matrix elements such as D = Φ| ψ † ǫ • rψ|Ψ are now evaluated by expressing the positive energy Dirac wavefunctions |Φ and |Ψ in terms of the Foldy-Wouthuysen ones |φ and |ψ : |Φ = U † |φ and |Ψ = U † |ψ . Since U is written as a many-body operator, we translate the quantum field expression for D the many-body formula D = Φ|ǫ • R|Ψ , where R = N n=1 r n . 111,112 We calculate D = φ|U ǫ • RU † |ψ , where U = e iS by using the Baker-Campbell-Hausdorff formula

ee

  itn |φ n φ n |. Thus, λη = (ηλ) † = π -π e -it P t dt = π -π e it P -t dt,and, by unicity of P t , ηλ = η(ηλ) † η implies P t = ηP -t η.We can now define a unitary square root U of ηλ by functional calculus: -it/2 ηP t ηdt = π -π e -it/2 P -t dt = U.

  Foldy-Wouthuysen wavefunctions |ψ = U |ψ D corresponding to positive energy have only even components. Indeed, let |ψ D be an eigenstate of H D corresponding to a positive energy. By definition of λ we have Γ + |ψ D = |ψ D and Γ -|ψ D = 0. Thus, U Γ + |ψ D = U |ψ D = |ψ and U Γ + = B + U implies |ψ = B + U |ψ D = B + |ψ . Thus η|ψ = ηB + |ψ = B + |ψ = |ψ and |ψ is even. Similarly 0 = B -|ψ , so that the odd part of |ψ is zero.For a one-body system, even components and large components are identical. Indeed a Dirac one-body wavefunction can be written|ψ D = φ ψ ,If η = β, then the even part and the odd parts of |ψ D are, respectively,

  using matrix elements of H ǫ (t) between eigenstates of H 0 . Then, cross-sections are expressed in terms of the S-matrix and T-matrix elements related by: m|S|n = δ mn -2iπδ(e m -e n ) m|T |n . |p and |n are eigenstates of H 0 with energy e m , e p and e n , respectively. The term iγ was added as a heuristic way to avoid divergence at resonance (i.e. when the states |n and |p are degenerate). More sophisticated methods exist to deal with such degeneracies[START_REF] Brouder | [END_REF] but they would bring us too far. From the physical point of view, γ describes the life-time of the state |p , which can decay by radiative or non-radiative relaxation.

	Up to second order,
	m|T |n = m|H eγ |n +

p m|H eγ |p p|H eγ |n e p -e n + iγ ,

(6)

where |m ,

Appendix: General Foldy-Wouthuysen transformation

D. H. Kobe and A. L. Smirl, Am. J. Phys.46, 624 (1978).

K.-H. Yang, J. Phys. A: Math. Gen.15, 437 (1982).

N. Bouldi, N. J. Vollmers, C. G. Delpy-Laplanche, Y. Joly, A. Juhin, P. Sainctavit, C. Brouder, M. Calandra,

Acknowledgments

We are very grateful to Uwe Gerstmann, Matteo Calandra and Nora Jenny Vollmers for encouraging us to work on the problem of the relativistic effects in x-ray absorption spectroscopy. We thank Alexander Silenko for this help concerning the Foldy-Wouthuysen transformations. Discussions with Amélie Juhin, Sergio Di Matteo, Yves Joly and Philippe Sainctavit are gratefully acknowledged. We are very grateful to Maria Esteban for her guidance through the mathematical literature on the many-body Dirac equation.

This work was supported by French state funds managed by the ANR within the Investissements d'Avenir programme under Reference No. ANR-11-IDEX-0004-02, and more specifically within the framework of the Cluster of Excellence MATISSE led by Sorbonne Universités.


Yannick Morel

Vincent Lebastard

Frédéric Boyer

Neural-based Underwater Spherical Target Localization through Electrolocation

Navigation of cluttered underwater environments remains to this day a challenging task in mobile robotics. Applying an electric field to a mobile robot's direct environment and measuring perturbations of this field, one is able to detect the presence of obstacles in close proximity of the system. In addition, one is also able to infer a range of information relative to the detected objects, such as their position or electrical characteristics. Extracting such information from available measures typically requires a model (analytical, numerical or heuristic) descriptive of the relationship from geometry of the scene to measures performed (typically referred to as forward model), or of the inverse relationship (inverse model). In the following, we directly extract one such model from experimental data, and capture a forward model using a neural formalism. Then, using an iterative procedure, we are able to estimate the position of a detected object and assess the degree of confidence one can place on this estimate. Merit of the approach is illustrated using experimental data for a spherical obstacle.

I. INTRODUCTION

Autonomous or remotely controlled robotic operations in unstructured, cluttered environments have long been recognized as constituting a significant challenge ( [START_REF] Thorpe | Field robots[END_REF]). The situation is to a large extent due to difficulties in using acoustic methods, such as SOund Navigation And Ranging (SONAR), in such a setting, as the clutter leads to severe multipath issues (see, among many others, [START_REF] Knight | Digital signal processing for SONAR[END_REF] for a discussion of the manner in which multipath affects SONAR). Operational problems are only compounded in the case of murky, muddy, or turbid water, in which conditions vision is of little assistance in perceiving the environment. As a result, a wide range of applicative situations, in which robotics could provide invaluable support, have remained beyond the capabilities of existing technology, such as for instance in active environmental monitoring for rivers, ponds, and lakes, exploration of underwater wrecks, or search and rescue missions in capsizing ships.

To overcome such limitations, roboticists have over the years investigated novel sensing modalities. Of particular interest in the aforementioned setting is the electric sense ( [3]- [START_REF] Der Emde | Active electrolocation in Gnathonemus petersii: Behaviour, sensory performance, and receptor systems[END_REF]). This active sensing modality, exploited by a number of fish species, is based on the manipulation of electric fields. Specifically, the fish generates one such field and continuously measures deformations of this field using electroreceptor located along its body. The presence of objects (inert or biotic), whose electric properties differ from those of the surrounding fluid, will perturb the field. Based on the information collected, the fish is typically able to infer detailed representations of its immediate surroundings, including the presence, location, and general geometry of foreign objects (obstacles), or of living organisms (prey). This sensing modality is characterized by a fairly short range (for the fish, typically a one-body length detection bubble around its body), however it has the advantages of being omnidirectional, well-suited to cluttered environments, and able to operate in turbid water. For these reasons, a number of research groups have investigated the application of a similar sensing principle to robotics (see for instance [START_REF] Mciver | Towards a biorobotic electrosensory system[END_REF], [START_REF] Boyer | Model for a sensor inspired by electric fish[END_REF]).

The typical functionality pursued has been that of underwater obstacle detection and localization (referred to hereafter as electro-location). Similarities between the considered problem and that investigated in Electrical Impedance Tomography (EIT, [START_REF] Cheney | Electrical impedance tomography[END_REF]) have motivated the use of EIT methods (see [START_REF] Snyder | Underwater object tracking using electrical impedance tomography[END_REF]). However, the number of electrodes (and therefore of available measures) on electro-location sensors is fairly limited. As a result, the inverse problem (of estimating the impedance map in the sensor's surroundings as a function of electric measures performed) is generally ill-posed ( [START_REF] Calderón | On an inverse boundary value problem[END_REF]). To reduce the dimension of the solution being estimated, it has in some instances been assumed that, in the sensor's direct environment, a significant portion of the volume is occupied by the surrounding fluid, whereas only a discrete number of obstacles are present. Therefore, the problem can be simplified by estimating the position of these obstacles, rather than attempting to construct a detailed impedance map of the environment. Such approaches have been very successful, with in particular techniques based on probabilistic estimation methods ( [START_REF] Solberg | Active electrolocation for underwater target localization[END_REF]), or stochastic approaches ( [START_REF] Lebastard | Environment reconstruction and navigation with electric sense based on a kalman filter[END_REF]).

In the following, we pursue a similar approach, in which we attempt to reconstruct the position of a given obstacle, present in the sensor's environment. The methods builds upon the technique used in [START_REF] Morel | Estimation of relative position and coordination of mobile underwater robotic platforms through electric sensing[END_REF], [START_REF] Morel | Actionperception trade-offs for anguilliform swimming robotic platforms with an electric sense[END_REF] to estimate relative positions in a group of robotic swimmers using the electric sense. In a first step, a forward map, from relative position of an obstacle with respect to the sensor to the electric measures performed by the sensor, is constructed. This forward map can assume a variety of forms, such as an analytical model (as that in [START_REF] Boyer | Model for a sensor inspired by electric fish[END_REF] or in [START_REF] Boyer | Exploration of objects by an underwater robot with electric sense[END_REF]), or a numerical model (e.g. based on Boundary Element Methods, BEM, [START_REF] Banerjee | The Boundary Element Methods in Engineering[END_REF], [START_REF] Wrobel | The Boundary Element Method[END_REF], such as that in [START_REF] Porez | Multi-physics model of an electric fish-like robot: Numerical aspects and application to obstacle avoidance[END_REF]). In the following, we use experiments to collect data descriptive of this forward map. Then, we train a polynomial neural network ( [START_REF] Shin | Approximation of multivariate functions using ridge polynomial networks[END_REF], [START_REF] Chak | Orthogonal polynomials neural network for function approximation and system modeling[END_REF]) to learn this forward map. In a second step, when attempting to localize an obstacle, series of candidate (or prospective) obstacle positions are presented as inputs to the trained neural net. The computed neural outputs (which correspond to predicted electric measures for each considered candidate positions) are then compared with the actual electric measures. The comparison allows to discern likelier positions of the obstacle.

The main contribution of the paper lies in presenting a novel electro-location technique for obstacle localization, technique which had to this day only been applied to the relative position estimation problems for two mobile systems equipped with the electric sense (or, equivalently, to two electric sensors). One of the main strengths of the approach lies in its expected generalizability. In particular, the problem of developing analytic models describing the effect of the presence of simple (in a geometrical sense, e.g. spherical) obstacles on electric measures has been addressed (see [START_REF] Boyer | Model for a sensor inspired by electric fish[END_REF], [START_REF] Rasnow | The effect of simple objects on the electric field of Apteronotus[END_REF]). However, that same problem remains open for more complex object geometries, be it for localized obstacles of arbitrary shape, or continuous, non-homogeneously curved walls. The corresponding analytical derivations are likely to be rapidly become untractable, and computing numerical BEM models of arbitrarily detailed scenes can prove prohibitively time-consuming. Yet, the electric sense is, to a large extent, a rather myopic sense. Specifically, if the presence of an obstacle can be discerned from available measures at a given range, distance to the object must be reduced to close-to-zero to distinguish geometrical features. For example, presence of a cube or of a sphere produces electric measures which are essentially identical at a reasonable range (see [START_REF] Lebastard | Localization of small objects with electric sense based on Kalman filter[END_REF]). In other words, from far enough, most objects look alike from an electric point of view. Accordingly, one could expect to be able to successfully address a large number of localization problems (i.e. localization for a large class of different objects) by exploring similarities between the electric measures corresponding to a number of generic shapes. This library of measures could be constructed using analytical or numerical models. However, following the approach used hereafter and collecting data experimentally instead, one is able to circumvent issues stemming from complexity of derivation or computation, and achieve greater flexibility in terms of the types of shapes considered. Establishing one such library and using it to estimate geometry of a complex scene is however beyond the scope of the present paper. Instead, in the following, we show that the proposed approach can be used to localize a spherical obstacle using experimental data.

In Section II, we provide a rapid description of the experimental setup and compare collected measures to values obtained from an analytical model. Then, in Section III, we use a neural approach to capture the map going from relative position of the obstacle to electric measures. Inverting this forward map, we are able to localize a sphere as a function of performed electric measures, as discussed in Section IV. Results obtained using experimental data are provided. Section V concludes this paper.

II. EXPERIMENTAL SETUP AND COLLECTED MEASURES

In later sections, we address the problem of localizing an object of interest, underwater, through manipulation of electric fields. The setting in which we have developed and Fig. 1. Experimental setup, top down perspective. The system is shown with a pair of electric probes, visible as cylindrical objects in white and grey. The probe in the foreground is mounted on a three DoF actuation system, allowing to control its position and orientation in the horizontal plane. tested the approach is described hereafter. At later stages of development, the technology is intended to equip mobile underwater vehicles. However, at this early developmental stage, the principles are tested on a dedicated experimental setup, shown in Figure 1.

A. Experimental Setup

The physical principle behind the object localization procedure is the following. A set of electrodes are used to generate an electric field and measure relevant electrical quantities. Specifically, in the setup used in the following, these electrodes are integrated within a slender cylinder (of length 21.6cm, diameter 2.5cm, probes seen in Figure 1 and Figure 2, see [START_REF] Boyer | Model for a sensor inspired by electric fish[END_REF] for greater details). The electrodes are made of steel, the body of the cylinder itself is plastic and non-conductive. This cylinder is hereafter referred to as the probe or the electric probe. These probes are used within a 1m×1m×1m glass tank, filled with fresh water. They are affixed at the extremity of a perch, rigidly connected to a three Degree of Freedom (DoF) actuation system, allowing movements in a horizontal plane (two translations, one rotation). A top-down view of the setup is provided in Figure 1. In the following, we will consider localization of a single object, a metallic sphere of radius 1cm.

The probe used features four pairs of electrodes. Two of these pairs form the probe's hemispheric extremities (see Figure 1 or Figure 2), the other two pairs are located on the fore part of the probe, their centers at 6.8cm and 3cm from the probe center, respectively. Each pair follows a strict leftright symmetry (see [START_REF] Boyer | Model for a sensor inspired by electric fish[END_REF], [START_REF] Lebastard | Localization of small objects with electric sense based on kalman filter[END_REF]). A signal generator is used to apply a sinusoidal signal of frequency 22.5KHz to the aftmost hemispherical pair of electrodes. The (peak-to-peak) value of the current flowing through the three fore-most pairs, Fig. 2. Schematics of the experimental setup, the origin of the xy-frame of reference is placed at the center of the tank, a conductive (metallic) sphere of radius 1cm is located at (xo, yo) = (0m, 0.25m), the position (x, y) of the probe may span the area in a darker shade of blue.

which are connected to a common ground, is measured. For a more detailed description of the setup, the interested reader is referred to [START_REF] Lebastard | Localization of small objects with electric sense based on kalman filter[END_REF], [START_REF] Boyer | Model for a sensor inspired by electric fish[END_REF].

Position of the obstacle and of the probe in the horizontal plane is described using an xy-frame of reference, the origin of which is located at the center of the setup (see Figure 2). To collect data to construct a forward map (from relative position to collected measures), we dispose the obstacle and probe as follows. The spherical obstacle remains fixed, its center located at (x o , y o ) = (0m, 0.25m). The probe is free to move about the horizontal plane. However, for security purposes (in particular to prevent collisions with walls or the spherical obstacle), we will in general constrain its movement to the darker blue area in Figure 2; that is, to (x, y) ∈ [-0.4m, 0.4m]×[-0.4m, 0.23m]. More specifically, to collect data the probe was made to travel along series of straight lines parallel to the x-axis in Figure 2, from y = 0m to y = 0.225m with 5mm steps (further measures in the direction of y decreasing are of limited use, as the effect of the obstacle's presence is not seen in the measures). For each such y-position, the probe is made to travel from x = -0.37m to x = 0.37m, with a speed of travel and a sample frequency such that we obtain a sample every 1.2mm traveled in the x direction. To isolate the contribution of the obstacle and eliminate influence of the walls, we repeated the same procedure without obstacle, and subtracted that reference value from the measures collected with the obstacle present. An example of the obtained measures is shown in Figure 3.

B. Analytical Model

For completeness and to offer additional insights into the manner in which presence of the obstacle impacts performed measures, we briefly present the analytical model of the currents measured, whose derivation is presented is detail in [START_REF] Boyer | Model for a sensor inspired by electric fish[END_REF] (alternately, see [START_REF] Boyer | Exploration of objects by an underwater robot with electric sense[END_REF] for fewer details and a higher level of abstraction). It is convenient to consider each pair of electrodes along the probe's body separately, and for any such pair to estimate (or measure) the total currents flowing through it, and the difference between the current flowing into the left electrode and that flowing into the right one. The former (the sum) produces signals of greater magnitude, which facilitates detection of a perturbation of the electric field. The latter (the lateral difference) provides more detailed information regarding the detected object's relative position with respect to the probe, in particular allowing to straightforwardly assess whether it is located port or starboard. For ease of presentation, we group the four current sums within i lo ∈ R 4 , where lo stands for longitudinal and in which current sums are included from fore-most to aftmost pair. Similarly, we group the lateral current differences within i la ∈ R 4 ( la stands for lateral) in the same front-toback order. In practice, we are interested in the perturbation due to the presence of an obstacle. Accordingly, we typically subtract from i lo the values obtained in the absence of any obstacle (which serve as reference). That operation is not necessary for i la , due to the sensor's left-right symmetry. The value of i lo (due to the presence of a spherical obstacle) can be computed as follows ( [START_REF] Boyer | Model for a sensor inspired by electric fish[END_REF]),

i lo = C 0 KC 0 u, (1)
where C 0 ∈ R 4×4 is a conductance matrix whose expression depends on the probe's geometry, u ∈ R 4 is the vector of electrical potentials applied to the four pairs of electrode (in our case u = 0 0 0 u 0 T , where u 0 represents the signal applied to the aft-most electrode pair), and K ∈ R 4×4 is a matrix whose entries are of the form

k ij = - 1 4πγ r T i χa 3 r j r i 3 r j 3 , i = 1, . . . , 4, j = 1, . . . , 4, (2)
where k ij is the i th line-j th column entry of K, γ ∈ R is the conductivity of the surrounding fluid, a ∈ R is the spherical obstacle's radius, χ ∈ R its relative conductivity, and r i ∈ R 3 , i = 1, . . . , 4, is the vector going from the obstacle's position to the center of the i th pair of electrodes.

Similarly, the i th entry of i la can be computed as follows,

i lai = S ⊥ γχa 3 Āi r i 5 2(n T i r i )r i + (n T i r i⊥ )r i⊥ E 0 , i = 1, . . . , 4, (3)
where S ⊥ ∈ R is a scaling factor, Āi ∈ R characterizes the surface of electrode pair i, n i ∈ R 3 is a vector characterizing the direction of the surface of the i th electrode pair, r i⊥ ∈ R 3 is a vector of same magnitude as and perpendicular to r i , and E 0 ∈ R 3 is the electric field produced by the sensor, evaluated at the position of the obstacle, which can be Comparison between model and measures

x [mm] i1(x, y) [A]
Fig. 3. Comparison between experimental measures (represented using a linearly interpolated surface) and model predictions (marked with red dots).

computed as

E 0 = 1 4πγ 4 i=1 C 0 u i r i 3 r i , (4)
where u i is the i th entry of vector u. We used the model given by ( 1)-( 4) to compute the currents expected in a configuration identical to the one for which we collected experimental data. Experimental measures (shown as a linearly interpolated surface) are compared to model values (represented with red dots) in Figure 3. Represented is the first component of i lo , with the position of the obstacle kept constant at (x o , y o ) = (0m, 0.25m), and values of current shown for the sensor traveling over (x, y) ∈ [-0.3m, 0.3m] × [-0.20m, 0.225m]. The current is represented on the z-axis as a function of the xy-position of the probe.We note that the exact position of the peaks is slightly different between model and experiment, and the background current value (away from the peaks) is not as uniform on the experimental data as it is in the case of the analytical model. The latter is certainly to a large extent due to imperfect compensation of the walls' influence on experimental measures. Nevertheless, we note that general trends are a good match, which tends to confirm that measures are representative of the phenomenon investigated, as described by the above analytical model.

III. FORWARD NEURAL MODEL

The experimental data collected, as discussed in the previous section, is descriptive of the forward map f (•) going from relative position of the sensor with respect to that of the obstacle, to the current measures obtained for that relative position,

f : D x × D y → D i , x, y → i, (5)
where D x , D y ⊂ R, D i ⊂ R 6 , x, y, represent the position of the probe in the frame of reference shown in Figure 2 (position which, assuming the obstacle remains fixed at (x o , y o ), is representative of the sensor's relative position with respect to the obstacle), and i represents the vector of useful current measures obtained for a sensor position of (x, y). Note that the analytical model's formulation leads, for the type of probe considered (featuring four pairs of electrodes) to eight different current measures available. However, measures obtained on the fourth pair are typically a linear combination of that obtained from the other electrode and are thus non-informative. In the measure vector i(x, y), we therefore only include the first three values in i lo (x, y) and i la (x, y) (corresponding to values obtained from the foremost three pairs of electrodes).

Following the approach proposed in [START_REF] Morel | Estimation of relative position and coordination of mobile underwater robotic platforms through electric sensing[END_REF], [START_REF] Morel | Actionperception trade-offs for anguilliform swimming robotic platforms with an electric sense[END_REF], we use a neural formalism to capture the input-output map (5) as described by the experimental data collected. Specifically, from the form of model ( 1)-( 4), we note that so long as the obstacle's position differs from that of the center of any electrode pair (which in practice is always verified), then f (x, y) as given in ( 5) is a continuous function of its arguments. Hence, from Weierstrass' approximation theorem ( [START_REF] Jeffreys | Methods of Mathematical Physics[END_REF]), we can conclude that for any ε > 0 there exists a polynomial function f (x, y) on

D x × D y such that f (x, y) -f (x, y) < ε, (x, y) ∈ D x × D y . (6)
Hereafter, we use the following structure for f (x, y) (which corresponds to that of a Single Hidden Layer Neural Network, SHL-NN, [START_REF] Spooner | Stable Adaptive Control and Estimation for Nonlinear Systems[END_REF]),

f (x, y) ϕ(x, y)W, (x, y) ∈ D x × D y , (7)
where we choose ϕ(x, y) ∈ R p as a row vector of polynomial functions of x, y, and W ∈ R p×6 is a matrix of neural weights to be selected. The learning (of the forward map ( 5)) task is accomplished by selecting W in such a manner that (7) becomes descriptive of ( 5) and of the experimental data collected. A wide range of methods can be used to achieve that goal. Specifically, numbering the experimental measures discussed in the previous section from

1 to n, let X x 1 . . . x n T ∈ R n , Y y 1 . . . y n T ∈ R n , I(X, Y ) i T 1 = f T (x 1 , y 1 ) . . . i T n = f T (x n , y n ) T ∈ R n×6 , and define Φ(X, Y ) ϕ T (x 1 , y 1 ) . . . ϕ T (x n , y n ) T ∈ R n×p .
Then consider

W * Φ † (X, Y )I(X, Y ), (8)
where • † denotes the usual pseudo-inverse. The particular choice of weights W = W * minimizes in the Least Mean Square (LMS) sense the reconstruction error f (x i , y i )f (x i , y i ), i = 1, . . . , n. That choice of weights is used in the following. Note that, alternately, one could use Bayesian inference in the parameter space to estimate the maximum (a posteriori) likelihood W matrix, as a function of available measures (as in [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]). Such an approach presents the advantage of allowing to explicitly account for measurement noise when estimating W . However, in the case considered here, comparisons show that the accuracy of the obtained forward map is not noticeably improved over the LMS solution.

IV. INVERSION OF THE FORWARD MODEL: LOCALIZATION

Inverting the forward model identified in the previous section can be accomplished by considering the distance between different sample points in the measurement space D i . Specifically, for a given measure i(t), performed at time instant t 0, consider the reconstruction error

e(i(t), η c ) i(t) -f (η c ) , t 0, (9)
where η c x c y c T ∈ D x × D y represents a candidate position of the sensor. Then, consider the relative position estimate chosen as

η(t) argmin ηc∈Dx×Dy e(i(t), η c ), t 0, (10)
where η x ŷ T , η x y T ∈ D x × D y , are the estimated and actual positions, respectively. It can be shown (see [START_REF] Morel | Estimation of relative position and coordination of mobile underwater robotic platforms through electric sensing[END_REF]) that, assuming the inverse map is sufficiently smooth, the more accurate our neural model (i.e. the smaller ε in ( 6)), and the more closely f (η) matches the actual measure i, the closer our estimate η is to the actual relative position η.

In practice, the approach only works if the sensor is within detection range of the object; that is, if it is close enough that the perturbation to the electric field due to the object's presence is important enough to be registered (which depends on measurement noise levels, magnitude of the electric field produced by the sensor, and geometry of the object). When outside of this detection range, the object's presence has no influence on measures. The approach however, generalizes elegantly to such cases. In particular, to exploit the above insights to perform localization, we use grids of candidate positions. Computing the reconstruction error (9) over the entire grid, it becomes straightforward to assess likelier areas. Specifically, in the case that the information measured is sufficiently informative (i.e. if close enough), the approach yields results that explicitly discriminate the most likely candidate position (for which the error is minimal), from unlikely candidates (with high error). In the case that the measured information is not sufficient, the ambiguity is made apparent by the results. Typically, large swaths appear as being good possible candidates. In that manner, considering how specifically the method is discriminating a likelier candidate position, it becomes possible to not only estimate the obstacle's position, but also assess the degree of confidence we can have in the estimate produced. The results can be represented in a manner that is visually intuitive. In particular, considering the contour plot of the reconstruction error over an area of interest, assessing likely and unlikely positions can be done by distinguishing areas of low reconstruction errors from areas of higher reconstruction error. One such contour plot is provided in Figure 4.

In the case considered, a metallic sphere of radius 1cm is present in position (x o , y o ) = (0m, 0.25m). The probe is made to travel on a straight line, parallel to the x-axis. The probe's actual trajectory is shown in green (with initial and final positions marked with an × and a •, respectively). As the probe moves along the green trajectory, samples of i lo and i la are measured. Using the forward map discussed in Section III and the above inversion procedure, we estimate the position of the probe with respect to that of the obstacle. These estimates are marked with ×'s in Figure 4. They were obtained using an iterative procedure. Specifically, for each sample we began with a candidate grid spanning (x, y) ∈ [-0.3m, 0.3m] × [0m, 0.22m] with steps of 10mm in the x-direction and 5mm in y. Then, we used a second grid, centered on the most likely candidate of the previous grid (chosen according to (10)), with a grid width of 60mm in x and y, and 1mm steps. A third gird is then used, centered on the estimate obtained from the second grid, with widths of 10mm in x and y, and 0.1mm steps. We note that the obtained estimates are reasonably faithful, as they remain close to the actual trajectory. The contour levels shown in Figure 4 are those of the reconstruction error [START_REF] Snyder | Underwater object tracking using electrical impedance tomography[END_REF] for the final position shown (marked with a • ). The lower the error (dark blue), the more similar the reconstructed measures are to the available set of measures. The greater the error (bright red), the more dissimilar. Accordingly, positions appearing in dark blue are likely candidate positions, while positions appearing in other colors (from green, to yellow and red) are unlikely. This information can be exploited in different manners. In particular, the candidate of least error represents in some sense the most likely candidate position (marked with ×'s). However, in cases that there are vast areas of low error (such as is the case in Figure 4, with a swath of low values at about x = 100mm, from y = 0mm to 200mm, approximately), the specific position of this minimum is not necessarily very informative, and it is useful to consider the gradient of the error (or slope of the corresponding error surface, if representing e as a surface over x-y). Both aspects can be combined using a Luenberger-form estimator ( [START_REF] Luenberger | Observer for multivariable systems[END_REF]). In black, we show the result obtained following one such approach. The model used is limited to kinematics and is descriptive of the fact that the probe is traveling in a direction parallel to the x-axis, at an a priori known speed (of 5mm/s). Within this model, we included system noise on speed and heading (additive Gaussian noises, with means 0deg and 0.1mm/s, standard deviations 20deg and 0.75mm/s) to reflect movement uncertainty. In addition, we exploited the information inferred from electric measures to design corrective terms. A first term is chosen proportional to the difference between current xy-position of the Luenberger observer and that of the position of least reconstruction error. A second term (akin to a derivative term in a Proportional Derivative controller) is chosen proportional to the gradient of e(x, y) evaluated at the xy-position of the observer. While the probe travels along the line y = 175mm with an initial x-position at x = -296mm, the initial condition of the observer is chosen as (x, y) = (-270mm, 40mm), to illustrate convergence of the estimate produced to the actual trajectory of the probe. As seen in Figure 4, the information is sufficiently accurate for the observed trajectory (in black) to converge to a neighborhood of the actual one (in green).

Finally, note that, as previously discussed, the closer the probe is to the obstacle, the greater the amplitude of the measures and the easier the estimation of the position. To quantify the impact of range, we made the probe travel in straight lines parallel to the x-axis (from x = -300mm to 300mm), and increasing lateral distance progressively, starting at y = 220mm (lateral distance of 30mm from the obstacle), up to y = 150mm (lateral distance of 100mm). The mean positioning error for each trajectory is represented with a • in Figure 5, and the standard deviation is represented with a bracket (in red). Note that, as seen for instance in Figure 2, the probe is dissymmetric, with three pairs of electrode on the fore part, and only one located aft. The idea behind that specific design choice is that, if the probe figuratively represents a mobile vehicle, it is more interested in discerning what lies ahead (in the probable direction of movement) than what is behind. As a result, as the probe travels from x = -300mm to x = 300mm, information collected (and therefore localization), is significantly better when the obstacle (situated in x = 0mm) lies in front of the probe. To verify this, we show in blue in Figure 5 the same error means and standard deviations, but only accounting for travel from x = -300mm to x = 0mm. As expected, the localization errors are noticeably smaller.

V. CONCLUSION

The work presented concerns the localization of spherical objects using electrolocation; that is, through the process of applying an electric field to the sensor's direct environment, measuring perturbations of this field due to the presence of foreign objects, and inferring information relative to these objects based on measured information. The approach uses experimental data to capture, using a SHL-NN, the forward map from obstacle relative position to the corresponding expected electric measures. Then, when detecting the presence of an obstacle and attempting to reconstruct its position, the neural forward map is inverted using an iterative procedure. A series of grids of candidate positions are presented to the neural map, which outputs the corresponding expected current measures. Comparing measures associated to different candidate positions to the actual measure performed allows to discern likelier candidate positions. We show, using experimental data, that the approach allows to estimate the obstacle position when within range, with a mean positioning error that is function of the distance between sensor and obstacle (the closer the sensor, the better the estimate). The main contribution of the paper lies in presenting a new solution to this electrolocation problem. In addition, this solution presents a number of interesting features when compared to existing alternate techniques. In particular, the approach is readily generalizable to any particular shape of obstacle, provided the corresponding data can be collected (to build the forward map). Conversely, typical methods rely on either analytical or numerical models, which, for complex geometries may become either untractable or require prohibitive computational efforts. In addition, the approach provides a clear map of space around the sensor, which allows to intuitively discern (a) what is the obstacle's likeliest position (if any), (b) what degree of confidence one may place on the estimate, and (c) if the estimate is uncertain, what areas of space are likelier to be occupied and which ones are likelier to be free. In the future, the approach will be extended to consider scenes composed of several objects. Specifically, the approach presented here will be supplemented by a neural classifier ( [START_REF] Haykin | Neural Networks: A Comprehensive Foundation[END_REF]) and rely on a library of captured neural models corresponding to different generic objects' forward maps (such as the one used here, which corresponds to spheres of a given radius). When presented with a new scene, in a first stage, the approach will use the neural classifier to classify the type of scene perceived (in terms of number and type of objects). Then, the localization approach proposed here will exploit the available library of forward maps to compose an estimate of the scene under consideration.

Fig. 4 .

 4 Fig. 4. Level curves of the reconstruction error e(x, y), indicating likelier (blue, low error) and less likely (in red, high error) positions of the probe; actual probe trajectory in green, instantaneous estimates marked with cyan ×'s, Luenberger observer trajectory in black.

Fig. 5 .

 5 Fig.5. Localization error for increasing lateral distances of the obstacle, means marked with dots, standard deviations represented with brackets; in red, values obtained with the probe traveling from x = -300mm to 300mm, in blue, with the probe traveling from x = -300mm to 0mm.