domain: phys chim

text:
Amelia Domínguez-Celorrio

Carlos Garc-Fernandez

Sabela Quiroga

Peter Koval

Veroniq Langlais

Diego Peña

Daniel Sánchez-Portal

David Serrate

Jorge Lobo-Checa

On-surface synthesis of Mn-phthalocyanines with optically active ligands †

The synthesis of novel organic prototypes combining different functionalities is key to achieve operational elements for applications in organic electronics. Here we set the stage towards individually addressable magneto-optical transducers by the on-surface synthesis of optically active manganese-phthalocyanine derivatives (MnPc) obtained directly on a metallic substrate. We created these 2D nanostructures under ultra-high vacuum conditions with atomic precision starting from a simple phthalonitrile precursor with reversible photo-induced reactivity in solution. These precursors maintain their integrity after powder sublimation and coordinate with the Mn ions into tetrameric complexes and then transform into MnPcs on Ag(111) after a cyclotetramerization reaction. Using scanning tunnelling microscopy and spectroscopy together with DFT calculations, we identify the isomeric configuration of two bi-stable structures and show that it is possible to switch them reversibly by mechanical manipulation. Moreover, the robust magnetic moment brought by the central Mn ion provides a feasible pathway towards magneto-optical transducer fabrication. This work should trigger further research confirming such magneto-optical effects in MnPcs both on surfaces and in liquid environments.

I. Introduction

Organic electronic devices are nowadays an industrial reality due to their high efficiency, low-cost and synthetic design tunability. 1-3 Among them, phthalocyanines (Pcs) and porphyrins stand out as prototypical examples of organic semiconductor building units due to their electric, magnetic and optical adjustability by metallation of their central cores. 4 These metal-organic complexes are used as pigment catalysts and photo-sensitive semiconductors besides being promising candidates for molecular spintronic applications. 5,6 When organic complexes are adsorbed on surfaces, important interfacial interactions emerge that influence their electronic and magnetic properties. [5][6][7] In the case of metallated Pcs and porphyrins with first-row transition-metal (TM) ions, significant conformational and spin state modifications have been reported. [7][8][9][10] These are strongly dependent on the local adsorption sites of metal centers and the arrangement of their external ligands, which alter the relative position of the metal center with respect to the π-conjugated macrocyclic plane. Similarly, selective adsorption of external and axial ligands on the metal centers modify these metal-organic complex properties. [11][12][13][14] Such local variability limits the scalability of single-molecule device production since building blocks require dynamical control and operational reversibility while maintaining their overall chemical and physical environments.

A way to overcome these local conformation limitations could be inserting active molecular switches into molecular complexes so that morphological control can be exerted using external stimuli, in the form of photons, 15,16 focused electric fields 17,18 or electric currents. 19,20 Light illumination is particularly simple to technically control, but demands that the spin centers of molecular complexes are strongly coupled to switchable ligands whose conformations can be changed by photon adsorption. The synthesis and deposition on the sur- † Electronic supplementary information (ESI) available. See DOI: https://doi.org/ 10.1039/d2nr00721e a Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain. E-mail: serrate@unizar.es, jorge.lobo@csic.es faces of such magnetic and optically active molecules by direct thermal evaporation is still pending, as it is compromised by the structural fragility of these metal-organic complexes. Such molecular integrity limitation has prevented its fabrication and implementation into magneto-optical transducers.

In this work, we achieved an adequate strategy to create metal-organic complexes by on-surface synthesis (OSS) featuring optically active moieties. In particular, we generated active Pcs with a central magnetic atom by means of cyclotetramerization of a functionalized phthalonitrile precursor on the Ag (111) surface. We chose Mn atoms as metallic centers based on their expected magnetic response that is evidenced in their reported Kondo signals on Ag substrates. 21,22 We find that the produced manganese phthalocyanines (MnPcs) have the central magnetic ions coupled to four diarylethylene (DAE) moieties. These MnPcs on Ag(111) show at each branch two possible conformations that we image and study electronically by scanning tunneling microscopy (STM) and spectroscopy (STS) techniques. Density functional theory (DFT) calculations allow us to identify the relevant configurations of the open and closed DAE isomers on the Ag(111) surfaces, and confirm the robust magnetic moment localized at the Mn ion. Notably, we induce reversible changes between the open and closed conformations by lateral tip-molecule manipulation. Although magneto-optical experiments were not attempted here, this OSS process ensures that each MnPc we formed contains four switchable branches.

II. Results and discussion

Optical activity of precursor 1

As a molecular precursor we designed and synthesized 4,5-bis (benzo[b]thiophen-3-yl)phthalonitrile that we name hereafter compound 1a. This compound was obtained in two steps by solution chemistry from commercially available compounds (see the ESI † for details). Compound 1a includes two vicinal nitrile (CN) groups to induce the on-surface cyclotetramerization reaction leading to phthalocyanine. It presents two vicinal benzothiophene groups forming a diarylethene moiety that should be involved in a reversible electrocyclization reaction leading to compound 1b (closed isomer), i.e. these moieties are introduced as an active switchable structure. 23 When a UV light of 365 nm wavelength irradiates a solution of compound 1a in CH 2 Cl 2 , the optical activity is photochromatically evident: the initial colorless solution corresponding to the open isomer 1a changes into blue as these molecules transform into the closed isomer 1b (Fig. 1a). The photochemically induced electrocyclic reaction involves six π-electrons, and therefore, according to the Woodward-Hoffmann rules the cyclization should be conrotatory, giving product trans-1b. The isomerization is reversible since the solution changes back to colorless once the UV irradiation is stopped (see ESI Fig. S1 †).

Dehydrogenation of the closed isomer 1b leads to the formation of conjugated dicyanonaphthalene 2 (Fig. 1b), which does not vary its color in solution when irradiated with the same UV light. Its unique isomeric configuration and planar geometry will be used as a control to ensure the presence of the optically active ligands of 1.

Integrity and isomeric configurations of 1 and 2 upon Mn codeposition on Ag(111) Isomeric molecules exhibiting active reversible switching have already been reported on metal surfaces responding to different stimuli, e.g. photons, 15,16 electrons 19,20 or electric fields. 17,18 These delicate ligand groups can withstand UHV evaporation whenever the molecules are small, since this permits sublimation at low enough temperatures that preserve their integrity. By contrast, larger molecules containing heavier metal centers tend to break apart at the active (and weakest) end groups.

The pathway that could lead to the generation of magnetooptical transducers requires to fulfill two conditions: first, deposition of optically active precursors that preserve their moiety integrity and, second, their combination with magnetic atoms. To achieve this, we evaporated small, optically active precursor molecules under UHV conditions so they self-assemble with Mn atoms on Ag(111). This surface is then used to catalyze a more complex structure by OSS. In order to prevent the Mn from clustering, we simultaneously deposited Mn atoms and precursor 1a on a Ag(111) surface kept at 50 °C (see the Methods section in the ESI †). As shown in Fig. 2a, squaredlattice islands are formed by the self-assembly of tetramer units that consist of Mn atoms (imaged as dark depressions) coordinated to four molecules by their cyano groups. 24,25 These metal-organic tetramers are quite stable since they conserve their integrity after lateral atomic manipulation with an STM tip (see the inset of Fig. 2a). The islands of Mn + four 1a ligands are stabilized by the formation of hydrogen bonds with adjacent molecules (see Fig. S2 †). In the coordinated islands, STM topography reveals two types of asymmetric ligands with different external lobes, also found for other diaryl derivatives on Ag(111). 26 For convenience, we will refer to the brighter lobes as A-configuration and the darker lobes as B-configuration (cf. red and green arrows in the inset of Fig. 2a). These lobes are randomly distributed and for an island with 452 molecules of 1 we find an A/B ratio of 0.72, as if the interactions between neighboring tetramer units were independent of their external conformations. This randomness suggests that both configurations are structurally very similar, as expected from the two isomers 1a and 1b sketched in Fig. 1a.

An identical process of co-depositing Mn + compound 2 on Ag(111) leads to very different structures. The similarity with the previous metal-organic complex is restricted to the formation of coordinated units of Mn atoms and four molecules by their cyano groups 24,25 (cf. Fig. 2b). However, these complexes are monodispersed and do not self-assemble into islands. Moreover, the ligands do not present visible differences in the external groups. We assign the absence of island aggregation to the lack of the two additional hydrogen atoms at their external DAE moieties compared to 1a and to the enhanced structural conjugation of the 2 precursors. All these findings confirm the molecular integrity of both 1 and 2 on the Ag(111) surface and the coordination with Mn into tetramer units that sets the stage for the MnPc formation by OSS.

MnPc formation by cyclotetramerization of metal-organic complexes

Post-annealing of these Mn complexes leads to the cyclotetramerization process of the ligands around the Mn atom 24,25,27 (see the Methods section in the ESI †). The cross shape observed in STM (cf. inset Fig. 2c) follows the MnPc structure sketched in Fig. 1c with four branches associated with the precursor molecules. Interestingly, the central Mn is now imaged as a bright protrusion, implying a significant electronic reconfiguration after the annealing process. The synthesized switchable MnPcs (3) still arrange into squared selfassembled islands (Fig. 2c), whereas the conjugated MnPcs (4) continue to be spread out on the surface (Fig. 2d). This difference should relate once again to the presence or absence of hydrogen atoms at the bisbenzothiophene moieties.

Focusing on 3, from the high resolution images of individual molecules extracted from the island by atomic manipulation (Fig. 3a-d), we still identify two apparent configurations of the peripheral DAE groups. We keep naming A for the brighter lobes (e.g. top branch in Fig. 2c inset) and B for the darker ones. These lobes must correspond to different benzothiophene arrangements since these two configurations are randomly present in the MnPc branches. These MnPcs (3) are far from being planar, with A lobes being further corrugated than B lobes. Indeed, recording constant-height and bondresolved images with a functionalised CO-tip (see Fig. S3 †) is only feasible by vertically adjusting the tip height for each kind of lobe. Interestingly, the external branches in 3 do not show mirror symmetry with respect to their adjacent pyrrole ring (see Fig. 3). The superimposed chemical structure on a switchable MnPc with two A and two B branch configurations duplicated in Fig. 3c andd includes the branch configurations with the lowest energies (see DFT calculations in Table 1). The open-rot configuration fails when fitting the underlying topo- By contrast, the cyclotetramerized molecules 4 in Fig. 2d † present four similar and symmetric branches with a rather homogenous intensity distribution. The formation of a cyclodehydrogenated MnPc with its expected planar structure is confirmed by bond resolved STM images (see Fig. S3 †), corroborating the cyclotetramerization reaction process of these MnPcs on the surface. Moreover, the differences in the STM images of 3 confirms the integrity of DAE moieties in the latter, based on the A/B branch lobe detection that maintains a bimodal distribution.

Electronic structure by STS of the switchable MnPc

The identified A/B conformers in the four branches of 3 must have associated with different electronic states, even if they are isomers. To confirm this, we acquired constant height differential tunneling conductance (dI/dV) spectra at selected positions on several individual molecules extracted from a compact island (see Fig. 4a). For the B lobe (darker branch), we observed two occupied states at -1.3 V and -0.3 V, an increase of the conductance around 0.41 V and two clearly resolved peaks at 1.6 V and 2.3 V. This is in contrast to the rather featureless A branch spectrum in this energy range, which presents its dominant peak at 3.1 V (see ESI Fig. S5 †). Such electronic differences between branches are evident in the topo- graphic images and the dI/dV map at -1.3 V (Fig. 4b ande). The HOMO-LUMO gap is significantly smaller for the B branch (0.7 eV) than for the A branch, as can be expected for the closed isomer. Interestingly, the Mn spectrum is independent of the configuration of the branches and presents a sharp peak right below the Fermi level with a ∼50 mV linewidth and a broad resonance at 2.3 V.

DFT calculations of isomeric MnPc (3) adsorbed on Ag(111)

At this point we must correlate the A and B electronic structures of the DAE moieties of 3 with their atomic arrangement. To do so, we performed DFT calculations of fully relaxed isomeric MnPc on a 9 × 10 × 3 Ag(111) supercell setting an equal atomic arrangement to the four branches (see the Methods section in the ESI †). Table 1 shows the relative energies of the three most stable configurations, where the minimum corresponds to the open arrangement (one sulfur atom pointing towards the surface and the other away from it). The second most favorable configuration is the open-rot that features a benzothiophene group rotated around its σ-bond by about 195°so that both sulphur atoms point to the surface. The third most stable configuration is the trans-closed, where the external hydrogen atoms are vertically opposed with respect to the MnPc central plane (one pointing towards the surface and the second away). Another identified possible configuration is the cis-closed isomer (see Fig. S8 †), where both hydrogen atoms point away from the surface. Note that this structure confers a distinct mirror symmetry to the branch, which was experimentally absent.

These DFT atomic arrangements must agree with the STM/ STS datasets. We can start by discarding the symmetric cisclosed arrangement shown in Fig. S8 † since none of the branches of 3 exhibits mirror symmetry with respect to the pyrrole ring (cf. Fig. 3a-d). In a similar way, the more compressed structural model of the open-rot isomer (superimposed to the MnPc in Fig. 3c) does not match the topographic STM footprint neither for A nor B lobes. Not only is energetically less stable than the open configuration, but it also shows almost identical LDOS for the molecular orbitals near the Fermi level (see ESI Fig. S9 †). Thus, it is justified to rule out the presence of the open-rot isomer and restrict our external atomic arrangements to the trans-closed and open configurations. Their sideview comparison in Table 1 shows a larger apparent height and a more pronounced asymmetry for the latter, which agrees with the corrugation found in the constant-current images in Fig. 3. Consequently, we associate the open arrangement with the A configuration and the transclosed with the B lobes.

To further increase the confidence in the assignment of the atomic arrangements of A/B lobes we calculated the local density of states (LDOS) in the energy window (-1, -0.5) eV corresponding to the highest fully occupied molecular orbital of 3 with four identical branches. The resulting LDOS for the open and trans-closed MnPcs is shown in Fig. 4c andd. The open branched MnPc exhibits a higher electronic density mainly localized at the pyrrole rings, whereas the trans-closed shows more weight towards the external branches. Moreover, our calculated PDOS (see ESI Fig. S9 †) support a shift to lower energies in the HOMO position for the trans-closed structure due to a higher conjugation of its branches with the pyrrole ring center. Then, based on this feature of the HOMO, we correlated the calculated molecular orbital of Fig. 4d for the transclosed geometry with the experimental dI/dV signal observed at -1.3 V of Fig. 4e. In contrast to constant current images, the darker B branches (number 4 in Fig. 4b) display a dominant character in the charge density distribution in the dI/dV maps of occupied states. Indeed, the trans-closed intensity maximizes closer to the sulfur atoms for the HOMO (cf. Fig. 4e). In contrast, the A branches (numbers 1, 2 and 3 in Fig. 4b) featuring brighter signals in the constant current, become practically irrelevant and only reflect the higher corrugation inherent to the open configuration in the dI/dV map. The LDOS of this occupied molecular orbital qualitatively agree with the features described for the calculated HOMO shown in Fig. 4c.

Isomer branch switching by mechanical manipulation

After identifying these two A/B configurations, we must demonstrate the feasibility of switching between the open and Table 1 DFT calculations of the isomeric MnPc structures of 3 adsorbed on Ag(111). After force relaxation, three atomic arrangements are identified (based on identical branch configurations) when adsorbed on Ag(111): open, open-rot and trans-closed. The minimum total energy structure corresponds to the open configuration and is used as a reference for the other cases. The last two columns indicate the maximum vertical separation omitting the hydrogen atoms (Δh max ) and Mn ion distance (Δz Mn ) relative to the open configuration and with respect to the last Ag slab of the surface. The insignificant differences in the Mn height result in practically identical total magnetic moments. The last two bottom rows show the chemical models of the corresponding branch configurations and the DFT calculated side views on a Ag(111) substrate (with the back branch removed for visualization purposes). A fourth conformation, the cis-closed configuration (shown in Fig. S8 † We also tested the use of bias pulses and slow ramps in the current or electric field to drive the molecular switches. 17,18,28,29 However, we discarded these methods since no changes were observed below a 3V threshold, whereas biasing the molecules above this value systematically led to non-reversible structural changes of the branches, in the form of dehydrogenation with the appearance of 4-type branches (see ESI Fig. S6 †).

Magnetic fingerprint of isomeric MnPcs

The second condition that we previously mentioned for the pathway to generate magneto-optical transducers is that our isomeric MnPc carries a well-defined magnetic moment. Based on previous work, 21,22 such properties can be validated through the detection of a Kondo effect on the Mn central atom. In MnPcs, the fingerprint of the Kondo effect has been identified by other groups as a resonance at the Fermi level with a linewidth of 13 meV convoluted with the Mn d xz/yz orbitals (characteristic of the D 4h squared planar symmetry group) appearing at -50 mV. 9,21,22 Therefore, using high-resolution dI/dV spectra in proximity of the Fermi energy we should verify the existence of the Kondo feature. The experimental data are shown in the inset of Fig. 4a. Indeed, our Mn STS reproduces the expected deconvolution analysis (see ESI Fig. S7d †), which can also be considered a further evidence of the correct MnPc formation by OSS.

DFT theory sheds light on the origin of this Kondo resonance that we locate at the pyrrole ring center, in agreement with previous work. 21,22 The calculated density of states pro- jected onto the Mn d orbitals indicate that the d x 2 -y 2 orbital is emptied, whereas the other four d orbitals are partially filled (see ESI Fig. S10 †). This electronic arrangement results in a total magnetic moment slightly less than 2μ B where the main contribution is associated with the Mn center but with an opposite magnetic contribution from the molecular framework. Previously, the Kondo resonance was associated with the many-body screening of the d z 2 orbital by the conduction electrons of the substrate, 21 which in our case is substantiated by the nearly half-filling observed for this orbital. Also, the degenerated d xz/yz orbitals exhibit a strong contribution just below the Fermi level, accounting for the tail to the left of the Kondo resonance.

A critical parameter for the Kondo state is the screening strength, which affects the Kondo temperature and resonance line shape when, for instance, altering the distance of the central Mn to the metallic substrate. Despite the possibility of switching the branch configurations between open and transclosed, we experimentally find for our switchable MnPcs that the low energy Mn spectra do not change (cf. ESI Fig. S7 †). This can be understood if the Mn vertical distance turns out to be insensitive to the branch configuration, as deduced by theory (cf. Δz Mn in Table 1). Furthermore, the identical planar environment of the Mn atoms in both configurations leads to nearly the same d orbital decomposition (ESI Fig. S10 †). The use of bulkier groups at the hydrogen positions that are still optically active could lead to significant Mn height changes when switching between different open/closed isomers, thereby bearing the prospect of fulfilling the magnetic transducer functionality.

III. Conclusions

We have demonstrated that the branches of individual MnPc synthesized on Ag(111) can be mechanically switched between open and trans-closed isomeric forms. Preparation of the molecular complex has been possible by the OSS of a smaller precursor (1), which exhibits optical activity in solution, but most importantly maintains its integrity upon UHV sublimation. Indeed, the assignment of the two distinct configurations (open and closed atomic arrangements) was possible by combining STM images, STS spectroscopy and DFT calculations. Notably, we demonstrate the prospect of active functionality of these DAE branches by the configuration switching using tip manipulation processes. We also find that the Mn ions of 3 exhibits a significant magnetic fingerprint, so an entirely new playground to fabricate individually addressable magnetooptical transducers at the atomic scale is envisioned and requires further experimental validation.

The flexibility of the OSS process to produce complex structures, where the optical activity is brought by diarylethylene (DAE) moieties and the magnetism by metal centers exchangeable by metallation, opens vast opportunities. For instance, we foresee the implementation of different readout signals changing the magnetic state of the ion, the magneto crystalline anisotropy, 30,31 the energy resolved spin polarization of the local DOS, 32,33 or the Kondo screening. 21,22 Also, we envision the use of bulkier ligands to alter the adsorption geometry and height of the pyrrole ring around the metal center, thereby correlating the Kondo resonance to the different combinations of isomer ligands. All these possibilities expand the well-stablished spin-crossover photo-magnetism approach at the individual molecular level.

Fig. 1

 1 Fig. 1 Molecular precursors and MnPc complexes of this work. (a) Scheme of the reversible reaction between open 1a and closed 1b isomers. (b) Molecular structure of the fully conjugated and non-switchable precursor dicyanonaphthalene 2 used for verification. (c) Molecular structure of the MnPcs synthesized by cyclotetramerization reactions. Complex 3 is generated from 1 and presents four switchable DAE moieties, whereas 4 is a rigid MnPc produced from 2. Green shading marks the identical phthalocyanine rings containing the central Mn, while blue denotes the bisbenzothiophene moieties.

Fig. 2

 2 Fig.2Representative STM images of the metal-organic structures. (a) Survey and closeups of precursors 1a (a) and 2 (b) co-evaporated with Mn atoms on Ag(111). Each Mn atom is surrounded by four ligands forming stable metal-organic coordination units. In contrast to (b), the Mn + 1 units cluster into islands (see the molecular model in Fig.S2 †) where each complex with four ligands can be manipulated and extracted without damage (cf. inset). The DAE moieties present two types of lobes (named A and B) which are marked by red and green arrows in the inset. Post-annealing the above structures to ∼300 °C and ∼350 °C induces cyclotetramerization of the ligands forming the MnPcs 3 (c) and 4 (d), respectively. As evidenced in the insets, the Pc core in 3 displays branches with two distinct appearances, whereas 4 is fully conjugated and structurally flat in comparison (see bond-resolved STM images in Fig.S3 †). STM images: (a) 40 × 40 nm 2 , 0.6 V, 60 pA, inset (a) 3.5 × 3.5 nm 2 , 0.6 V, 60 pA; (b) 40 × 40 nm 2 , -1 V, 100 pA, inset (b) 3.5 × 3.5 nm 2 , -1 V, 100 pA; (c) 50 × 50 nm 2 , 10 mV, 50 pA, inset (c) 2.5 × 2.5 nm 2 , 0.4 V, 100 pA; and (d) 50 × 50 nm 2 , 10 mV, 50 pA, inset (d) 2.5 × 2.5 nm 2 , -0.4 V, 20 pA.

Fig. 3

 3 Fig. 3 High resolution imaging and DFT structural models of 3. (a)-(c) STM images with a functionalized tip (2.5 × 2.5 nm 2 , -50 mV, 20 pA) of three different switchable MnPcs (3) isolated by lateral atomic manipulation (manipulation set point 5 mV and 15-40 nA). (d) is duplicated from (c) to compare two different superimposed molecular models from the three DFT branch configurations (see Table 1 and Fig. S4 †). The positions where the model does not fit their underlying topographic STM footprints are indicated by red arrows in (c). This corresponds to the open-rot, which does neither fit A nor B lobes (cf. Fig. S4 †). The correct molecular structure is overlaid in (d). The perpendicular dotted lines in (a) and (b) indicate the symmetrical axis expected for square-planar molecules.

Fig. 4

 4 Fig. 4 Experimental and calculated electronic structures of 3. (a) Position dependent dI/dV spectra (V mod = 10 mV) recorded at representative points on several isomeric MnPcs: Mn is shown in violet, the A lobe in red, the B branch in green and the substrate reference in gray. (b) Constant current STM image (V = -50 mV, I = 20 pA, 3.5 × 3.5 nm 2 ) with color dots indicating the reference positions of the STS in (a). (c) and (d) LDOS calculated isosurfaces of 3 for the highest fully occupied molecular orbital in the energy window (-1, -0.5) eV with the four branches in open (c) and trans-closed (d) atomic arrangements. The electronic density becomes more delocalized over the external branches on the trans-closed, indicative of conjugation at the external branches in comparison to the open configuration. (e) dI/dV map at -1.3 V (V mod = 10 mV) corresponding to (b), showing the spatial distribution of the occupied states.

Fig. 5

 5 Fig. 5 Mechanical switching of the DAE moieties by tip manipulation. STM images recorded before and after the controlled lateral drag of 3 over the Ag surface (-100 mV, 50 pA, 8.75 × 8.75 nm 2 ). The arrows point to the switched branches. In the case of (a) and (b) one DAE branch changes from the open to trans-closed configuration. Alternatively, the sequence in (c) shows a reversible switch.

d Donostia International Physics Center, Paseo Manuel de Lardizabal 4, E-20018 San Sebastian, Spain e Centro de Investigación en Química Biolóxica e Materiais Moleculares (CIQUS), Departamento de Química Orgánica, Universidade de Santiago de Compostela, Santiago de Compostela 15782, Spain f Simune Atomistics S.L., Avenida Tolosa 76, 20018 San Sebastian, Spain g

This journal is © The Royal Society of Chemistry 2022 Nanoscale, 2022, 14, 8069-8077 | 8073

Acknowledgements

We sincerely thank Marten Piantek and Carlos Casas for their contribution to this work. Dr Piantek is also acknowledged for providing ideas for this work and for in-depth discussions. We gratefully acknowledge financial support from the Spanish Ministries of Economy, Industry and Competitiveness (MINECO, grants MAT2016-78293-C6-4-R and MAT2016-78293-C6-6) and of Science and Innovation (MICINN, grant no. PID2019-107338RB-C62, PID2019-107338RB-C64, PID2019-107338RB-C66/AEI/10.13039/501100011033 and RED2018-102833-T), the regional Governments of Aragon (E13-20R and E12-20R), the Xunta de Galicia (Centro de Investigación de Galicia acreditation 2019-2022, ED431G 2019/03), the Dept. of Education of the Basque Government and UPV/EHU (grant no. IT1246-19), the Spanish Research Agency (AEI), the European Regional Development Fund (ERDF) under the program Interreg V-A España-Francia-Andorra (contract no. EFA194/16 TNSI) and the European Union through Horizon 2020 (FET-Open project SPRING grant. no. 863098). Finally, we acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

Author contributions

Conflicts of interest

The authors declare no competing financial interests.


Zebo Zhou

Zeliang Zhang
email: zeliangzzl@163.com

Shanhui Mo

Jin Wu

Hassen Fourati
email: hassen.fourati@gipsa-lab.grenoble-inp.fr

Online Calibrated, Energy-aware and Heading Corrected Pedestrian Navigation with Foot-Mounted MARG Sensors

Keywords: Pedestrian Navigation, Energy Management, Kalman Filter, Heading Correction, Zero-Velocity-Update Pedestrian Navigation, Energy Management, Kalman Filter, Heading Correction, Zero-Velocity-Update

The objective of this paper is to propose an energy-aware and ambulatory gait corrected pedestrian dead reckoning (PDR) approach using foot-mounted magnetic, angular rate and gravity (MARG) sensors. Compared with existing algorithms of PDR, the proposed method aims to solve three main problems for real pedestrian applications. First, to avoid limitations of off-line calibration for personal step length parameters, we utilize the zero-velocity-update (ZUPT) aided pedestrian MARG performance to continuously compute one's pose information. Meanwhile, it accumulates the moving distance for further estimation of one's step length during the initialization process. Secondly, due to different pedestrian gaits implicating the heading deviation angle between one's moving direction and heading, there are non-negligible impacts on pedestrian dead reckoning accuracy. The linear Kalman filter is used to recursively estimate the deviated heading during aforementioned initialization process. The third problem is related to the energy consumption. Following three aspects of adaptive energy saving work are devised: (i) energy-ware strategy for gyroscopes measurements acquisition is adopted to guarantee lower energy consumption. (ii) Modes switching mechanism of navigation computation is applied to the initialization and dead reckoning processes. (iii) De-sampling after initialization process has been invoked. Finally, real-world experiments are carried out to evaluate the performances of developed PDR system. The results show the efficiency of the suggested approach. The personal moving deviation angle obtained from the developed pedestrian navigation system can be potentially used for monitoring the patients' walking rehabilitation training and provide essential data for doctors to make constructive suggestions.

Introduction

Motivation

Nowadays, with the increasing demands of real-time personal location based service, location-aware consumer electronics have become more and more indispensable and popular in daily life [START_REF] Bancroft | Estimating MEMS Gyroscope g-Sensitivity Errors in Foot Mounted Navigation[END_REF][START_REF] Shi | A Robust Pedestrian Dead Reckoning System Using Low-Cost Magnetic and Inertial Sensors[END_REF]. Pedestrian navigation technology has been intensively studied and widely applied in numerous military and civilian applications, in particular, pedestrian localization, indoor navigation, motion capturing and body tracking [START_REF] Jin | A Robust Dead-Reckoning Pedestrian Tracking System with Low Cost Sensors[END_REF][START_REF] Tian | A Low-Cost INS and UWB Fusion Pedestrian Tracking System[END_REF][START_REF] Wang | A Real-Time Pedestrian Dead Reckoning System With FM-Aided Motion Mode Recognition[END_REF][START_REF] Tong | A Double-Step Unscented Kalman Filter and HMM-Based Zero-Velocity Update for Pedestrian Dead Reckoning Using MEMS Sensors[END_REF]. Global navigation satellite system (GNSS) provides users with satisfactory positioning accuracy in outdoor environments. However, it cannot work in GNSS challenging environments, e.g. canyons, tunnels, indoors, signal jamming or interruptions [START_REF] Fan | Improved Pedestrian Dead Reckoning Based on a Robust Adaptive Kalman Filter for Indoor Inertial Location System[END_REF]. Wireless signals, e.g. Wi-Fi, Bluetooth, and ZigBee generated by local area networks have been commonly utilized as potentially feasible supplementaries through fingerprinting and trilateration for augmenting navigation solutions under those GNSS denied environments [START_REF] Tian | A Low-Cost INS and UWB Fusion Pedestrian Tracking System[END_REF][START_REF] Gu | Trajectory Estimation and Crowdsourced Radio Map Establishment from Foot-Mounted IMUs, Wi-Fi Fingerprints, and GPS Positions[END_REF][START_REF] Steinhoff | Dead Reckoning from the Pocket -An Experimental Study[END_REF][START_REF] Liang | An Automatic Site Survey Approach for Indoor Localization Using a Smartphone[END_REF]. However, the radio signals are prone to suffer from multi-path effects, e.g. signal attenuation, diffraction and refraction, thus significantly degrading the signal quality [START_REF] Tong | A Double-Step Unscented Kalman Filter and HMM-Based Zero-Velocity Update for Pedestrian Dead Reckoning Using MEMS Sensors[END_REF][START_REF] Nilsson | Cooperative Localization by Dual Foot-Mounted Inertial Sensors and Inter-Agent Ranging[END_REF]. Moreover, state-of-the-art ground based radio augmentation systems require a complex sensor-network infrastructure to be deployed in advance.The network system usually contains numerous beacons for good accuracy and coverage purpose such that it is neither a lowcost nor an instant-deployment solution for pedestrian navigation applications [START_REF] Piccinni | An Improved Technique based on Zadoff-Chu Sequences for Distance Measurements[END_REF][START_REF] Piccinni | A Novel Distance Measurement Technique for Indoor Positioning Systems based on Zadoff-Chu Sequences[END_REF]. Therefore, the beacon-free and self-contained navigation technology is much more preferable without any pre-installed infrastructure.

Related Work

In the past decades, with the development of micro-electro-mechanical (MEMS) technologies, MEMS inertial sensors have been extensively used for pedestrian navigation [START_REF] Bancroft | Estimating MEMS Gyroscope g-Sensitivity Errors in Foot Mounted Navigation[END_REF][START_REF] Fan | Improved Pedestrian Dead Reckoning Based on a Robust Adaptive Kalman Filter for Indoor Inertial Location System[END_REF][START_REF] Brahms | Stride Length Determination during Overground Running Using a Single Foot-Mounted Inertial Measurement Unit[END_REF][START_REF] Valenti | Keeping a good attitude: A quaternion-based orientation filter for imus and margs[END_REF][START_REF] Meng | Self-contained pedestrian tracking during normal walking using an inertial/magnetic sensor module[END_REF].These sensors consisting of a tri-axial gyroscope and a tri-axial accelerometer are usually attached to a certain part of one's body for motion monitoring or biomedical rehabilitation purposes [START_REF] Prateek | Modeling, Detecting, and Tracking Freezing of Gait in Parkinson Disease Using Inertial Sensors[END_REF][START_REF] Gómez-Espinosa | Foot-Mounted Inertial Measurement Units-Based Device for Ankle Rehabilitation[END_REF][START_REF] Trojaniello | Foot Clearance Estimation during Overground Walking and Vertical Obstacle Passing using Shank-Mounted MIMUs in Healthy and Pathological Subjects[END_REF][START_REF] Ma | Cyclepro: A Robust Framework for Domain-Agnostic Gait Cycle Detection[END_REF][START_REF] Turcot | New accelerometric method to discriminate between asymptomatic subjects and patients with medial knee osteoarthritis during 3-D gait[END_REF][START_REF] Rampp | Inertial Sensor Based Stride Parameter Calculation from Gait Sequences in Geriatric Patients[END_REF].

A great deal of endeavor has been paid to improving the performances of pedestrian navigation. Existing representative methods of pedestrian navigation with inertial sensors can be mainly categorized into two sorts. One is strapdown-inertial-mechanization based method. Due to the inherent drift of gyroscopes and accelerometers, the navigation error of strapdown inertial system rapidly grows, inevitably leading to an important error [START_REF] Foxlin | Pedestrian Tracking with Shoe-Mounted Inertial Sensors[END_REF]. To solve this problem, a considerable research on zero-velocity-update (ZUPT) method has been conducted by introducing zero velocity as a pseudo-measurement during one's stance phase of walking. It is proven to be sufficient to compensate for the inertial sensor drifts and reset the accumulation errors [START_REF] Jiménez | Indoor Pedestrian Navigation Using an INS/EKF Framework for Yaw Drift Reduction and a Foot-Mounted IMU[END_REF][START_REF] Zhang | Adaptive Zero Velocity Update Based on Velocity Classification for Pedestrian Tracking[END_REF][START_REF] Wahlstrom | Zero-Velocity Detection-A Bayesian Approach to Adaptive Thresholding[END_REF]. To maintain good heading accuracy and limit the vertical axis drift of gyroscope, magnetometers are increasingly common to be assembled with inertial sensors on massive low-cost devices [START_REF] Wang | A Real-Time Pedestrian Dead Reckoning System With FM-Aided Motion Mode Recognition[END_REF][START_REF] Nilsson | Foot-Mounted Inertial Navigation Made Easy[END_REF][START_REF] Fourati | Heterogeneous Data Fusion Algorithm for Pedestrian Navigation via Foot-Mounted Inertial Measurement Unit and Complementary Filter[END_REF]. It is referred to as MARG sensor that contains a triad of magnetic, angular rate and gravity components. Many researchers are concerned with magnetic disturbances and inertial drifts estimation for MARG sensor module applications [START_REF] Wu | Real-time Magnetic disturbance Determination for Micro Air Vehicles via Gravity and Global Navigation Satellite System Measurements[END_REF][START_REF] Jirawimut | A Method for Dead Reckoning Parameter Correction in Pedestrian Navigation System[END_REF][START_REF] Roetenberg | Compensation of Magnetic Disturbances Improves Inertial and Magnetic Sensing of Human Body Segment Orientation[END_REF][START_REF] Wu | MARG Attitude Estimation Using Gradient-Descent Linear Kalman Filter[END_REF]. To optimally fuse data from this module, various linear and nonlinear Kalman filters (KFs) e.g. conventional KF, extended KF, unscented KF and particle KF are intensively developed to estimate the one's position and orientation [START_REF] Tian | A Low-Cost INS and UWB Fusion Pedestrian Tracking System[END_REF][START_REF] Nilsson | Foot-Mounted Inertial Navigation Made Easy[END_REF][START_REF] Xu | An Indoor Pedestrian Localization Algorithm Based on Multi-Sensor Information Fusion[END_REF][START_REF] Wu | A Pedestrian Dead-Reckoning System for Walking and Marking Time Mixed Movement Using an SHSs Scheme and a Foot-Mounted IMU[END_REF][START_REF] Xia | Indoor Localization on Smartphones Using Built-In Sensors and Map Constraints[END_REF]. The other sort of pedestrian navigation method is dead reckoning based method. It recursively integrates one's step length and orientation during each detected step. Large amount of previous studies have been performed focusing on step detection [START_REF] Guo | Indoor Positioning based on Foot-Mounted IMU[END_REF][START_REF] Jimenez | A Comparison of Pedestrian Dead-Reckoning Algorithms Using a Low-Cost MEMS IMU[END_REF] and step length model estimation [START_REF] Tian | A Multi-Mode Dead Reckoning System for Pedestrian Tracking Using Smartphones[END_REF][START_REF] Lan | Using Smart-Phones and Floor Plans for Indoor Location Tracking-Withdrawn[END_REF][START_REF] Ho | Step-Detection and Adaptive Step-Length Estimation for Pedestrian Dead-Reckoning at Various Walking Speeds Using a Smartphone[END_REF][START_REF] Zhang | A Handheld Inertial Pedestrian Navigation System With Accurate Step Modes and Device Poses Recognition[END_REF]. To further enhance the navigation accuracy and reliability, [START_REF] Tian | A Multi-Mode Dead Reckoning System for Pedestrian Tracking Using Smartphones[END_REF] developed a multi-mode pedestrian navigation algorithm by taking multi-placement of smartphones into account. More recently, multi-mode behaviors including walking, running, swing and standing still are considered as well in recent work [START_REF] Zhou | Behaviors Classification based Distance Measuring System for Pedestrians via a Foot Mounted Inertial Sensor[END_REF]. However, the existing dead reckoning methods mentioned above need to calibrate a group of parameters with sufficient personal off-line datasets prior to use. This is a main impediment for its being off-the-shelf product. Moreover, though various types of misalignment angles caused by inconsistency of installed sensor frame and defined body frame have been previously discussed and fixed [START_REF] Shi | A Robust Pedestrian Dead Reckoning System Using Low-Cost Magnetic and Inertial Sensors[END_REF][START_REF] Wang | Error Analysis of ZUPT-Aided Pedestrian Inertial Navigation[END_REF][START_REF] Manos | Gravity-Based Methods for Heading Computation in Pedestrian Dead Reckoning[END_REF], to the best of our knowledge, there are quite few related works discussing the heading deviation angle compensation related to personal gait in foot-mounted pedestrian navigation yet. It is worthy of mentioning that such heading deviation angle varies from persons and their gaits. Energy management is another critical issue for the wearable devices. For pedestrian navigation applications, gyroscope measurements are proven to be more reliable than accelerometers in either orientation determination or step detection [START_REF] Zhang | Novel Drift Reduction Methods in Foot-Mounted PDR System[END_REF]. However, in practice, it is universally acknowledged that the gyroscopes consume much more energy compared with accelerometers or magnetometers [START_REF] Makni | Energy-Aware Adaptive Attitude Estimation under External Acceleration for Pedestrian Navigation[END_REF]. For instance, according to the official user manuals of STMicroelectronics, for its typical product of MEMS gyroscope LSGD20 and accelerometer LIS3DH, the former consumes about 6.1mA while the latter is merely about 11 µA [START_REF] Makni | Energy-Aware Adaptive Attitude Estimation under External Acceleration for Pedestrian Navigation[END_REF]. On the other hand, sampling rate is also of great importance for MEMS inertial measurement unit (IMU) performances.

The authors in [START_REF] Wang | Error Analysis of ZUPT-Aided Pedestrian Inertial Navigation[END_REF] present an interesting results that a 2.5 times increase of velocity uncertainty by reducing the sampling rate from 1000 Hz to 100 Hz during ZUPT aided pedestrian inertial navigation. This strongly motivates us to develop an efficient and applicable energy saving strategy of MARG sensor module for pedestrian navigation without evidently affecting its performance.

Main Contributions

Up to now, much attention has been paid to pedestrian navigation particularly concentrating on step detection, step length modeling and sensor bias compensation but rarely with regard to efficient energy management, online pedestrian parameter calibration and gait related heading compensation. In this paper, the pedestrian navigation problem is solved under the following contributions (1) A joint ZUPT/dead reckoning based pedestrian navigation method. It is free of off-line calibration for personal PDR parameters. The ZUPT-aided Pedestrian MARG Navigation is utilized to continuously compute one's pose information meanwhile it accumulates the moving distance during the initialization process for further estimating one's step length in the way of dead reckoning.

(2) A gait related heading deviation angle compensation. Due to different pedestrian gaits implicating the heading deviation angle between one's moving direction and heading, there are non-negligible impacts on PDR accuracy. To solve such problem, a linear KF will be proposed to recursively estimate the deviated heading for adapting one's different behavior modes.

(3) An adaptive and efficient energy management strategy. Three aspects of adaptive energy saving work are devised: (a) Energy-aware strategy for gyroscope measurements acquisition is adopted to ensure the lower energy consumption. (b) Modes switching mechanism of navigation computation is applied to the initialization and dead reckoning processes. (c) Adaptive de-sampling during pedestrian dead reckoning process.

The proposed method has the following advantages:

(1) It is based on the energy-aware strategy thus consumes less energy.

(2) It identifies step length parameters for different users without any off-line tests.

(3) Most importantly, it adaptively corrects one's moving direction according to personal gait. Inversely, it can be potentially used for monitoring the patients' walking rehabilitation training.

Arrangement of Contents

This paper is organized as follows. In Section II, the principle and structure of proposed method is given and explained. In Section III, a gait related heading deviation angle is analyzed and compensated through a linear KF. The energy management strategy is devised in Section IV. In Section V, real experimental tests and comparisons are carried out. Finally, concluding remarks in Section VI end this paper.

Principle and Structure of Proposed Pedestrian Navigation Method

The overall structure of proposed pedestrian navigation method is presented in Fig. 1. It is clearly shown that two major parts are involved in pedestrian navigation computation: on-line initialization and gait based dead reckoning.

During online initialization process, the ZUPT-aided pedestrian MARG navigation is performed to continuously compute one's pose information meanwhile it accumulates the moving distance to achieve on-line calibration for personal step length parameters. Besides, due to different pedestrian gaits implicating the heading deviation angle (see δγ in Fig. 2) between one's moving direction and foot's heading, there are non-negligible impacts on pedestrian dead reckoning accuracy. Actually, it is not compulsory for sensor placement. If the sensor module is placed as shown in Fig. 3, the misalignment angle between forward axes of shoe and sensor can be still roughly measured by protractor beforehand.

To recursively estimate the heading deviation, a linear KF is constructed. After completing initialization, the pedestrian position will be updated with dead reckoning based method. Such structure can greatly saves the computation resources than conventional ZUPT-aided inertial mechanization methods.

ZUPT-Aided Pedestrian MARG Navigation

ZUPT-aided MARG navigation not only continuously estimates ones' position, velocity and attitude (PVA) information [START_REF] Nilsson | Foot-Mounted Inertial Navigation Made Easy[END_REF][START_REF] Nilsson | Foot-mounted INS for Everybody -An Open-source Embedded Implementation[END_REF] , but also plays a role in initialization part for calibrating personal parameters (e.g. step length and heading deviation). At epoch k, one's PVA can be easily computed with inertial measurements through

         p n k = p n k-1 + v n k ∆t + (C n bk • f k + g) (∆t) 2 /2 v n k = v n k-1 + (C n bk • f k + g) ∆t C n bk = C n bk-1 (I 3 + [ω k ×] ∆t) (1)
where subscripts k and k-1 represent epoch indexes; superscript n and subscript b represent navigation frame (denoted as n-frame, North-East-Down) and body The accelerometer and gyroscope models are established respectively as

y G = ω + b G + w G (2)
y A = f + b A + w A (3)
where subscripts G and A denote the gyroscope and accelerometer respectively;

ḃA = w rA (5)
where w rG and w rA are the zero mean Guassian white noises with variances Q wrG and Q wrA respectively. To limit the inertial accumulative errors, the ZUPT-aided MARG KF is constructed with the following 15-dimensional error state

x I = φ δp δv b G b A (6)
where φ, δp and δv denote misalignment angles, position error and velocity error vectors, respectively; superscript I and II are introduced to distinguish different Kalman filters. Then the state model can be established as

x I k = Φ I k-1,k x I k-1 + Γ I k w I k (7)
where the state transition matrix Φ I k-1,k and the noise transition matrix Γ I k are given by

Φ I k-1,k =            I 3 0 3×3 0 3×3 -∆tC n b 0 3×3 0 3×3 I 3 ∆tI 3 0 3×3 0 3×3 ∆t [f n ×] 0 3×3 I 3 0 3×3 ∆tC n b 0 3×3 0 3×3 0 3×3 I 3 0 3×3 0 3×3 0 3×3 0 3×3 0 3×3 I 3            Γ I k =            -∆tC n b 0 3×3 0 3×3 0 3×3 0 3×3 0 3×3 0 3×3 0 3×3 0 3×3 ∆tC n b 0 3×3 0 3×3 0 3×3 0 3×3 ∆tI 3 0 3×3 0 3×3 0 3×3 0 3×3 ∆tI 3           
The process noise w I k can be written as

w I k = w G w A w rG w rA
which obeys zero mean Guassian distribution, with variance

Q I k = diag Q w G,k , Q w A,k , Q w rG,k , Q w rA,k
where diag represents diagonal matrix. The predicted state and its covariance matrix are computed by,

x I k|k-1 = Φ I k-1,k x I k-1|k-1 (8)
P I k|k-1 = Φ I k-1,k P I k-1|k-1 Φ I k-1,k + Γ I Q I k Γ I (9)
Considering ZUPT and magnetic sensor output, accordingly the observation model is formulated as

z I k = H I k x I k + ε I k ( 10
)
where the observation vector z k and the corresponding design matrix H k are given by

z I k = ψ M,k -ψ G,k -v k H I k =   D 0 1×3 0 1×3 0 1×3 0 1×3 0 3×3 0 3×3 I 3 0 3×3 0 3×3  
where ψ G,k denotes heading angle derived from gyroscope while ψ M,k is computed from magnetometer outputs; D is the connection matrix between misalignment angles and Euler error angles

D = [ sin ψ tan θ -cos ψ tan θ 1 ] (see Appendix A). The observation model noise ε I k ∼ N 0, R I k is ε I k = ε ψ,k ε v,k in which ε ψ,k ∼ N (0, R ψ,k ) denotes the heading correction model noise; ε v,k ∼ N (0, R v,k )
denotes the velocity correction model noise, which is highly determined by the static threshold setting during the zero velocity detection. The observation variance can be written as

R k = diag{R ψ,k , R v,k }
Then the KF solution and its covariance are estimated by

x I k|k = x I k|k-1 + K I k z I k -H I k x I k|k-1 (11)
P I k|k = I 15 -K I k H I k P I k|k-1 (12)
with the gain matrix K

K I k = P I k|k-1 H I k H I k P I k|k-1 H I k + R I k -1 (13)
Finally by correcting the PVA with estimated δp, δv, φ and resetting error state, the ZUPT-aided MARG navigation is completed.

Gait Based Dead Reckoning

The principle of gait based dead reckoning is generally depicted as:

  p x,k p y,k   =   p x,k-1 p y,k-1   + d k   sin γ k cos γ k   (14)
where [p x,k , p y,k ] denotes the position after the k-th step ; d k and γ k represents the k-th step length and moving orientation respectively. For step length estimation, there are numerous models with various parameters. The most popular one is a bi-parametric linear model [START_REF] Sun | Activity Classification and Dead Reckoning for Pedestrian Navigation with Wearable Sensors[END_REF], in which the step length is assumed to be related to the step frequency. Apart from the step frequency, other parameters are introduced into the refined step length model as well, for instance, acceleration variance [START_REF] Zhang | Adaptive Zero Velocity Update Based on Velocity Classification for Pedestrian Tracking[END_REF], acceleration boundary [START_REF] Ho | Step-Detection and Adaptive Step-Length Estimation for Pedestrian Dead-Reckoning at Various Walking Speeds Using a Smartphone[END_REF], height [START_REF] Tian | A Multi-Mode Dead Reckoning System for Pedestrian Tracking Using Smartphones[END_REF] and leg length [START_REF] Lan | Using Smart-Phones and Floor Plans for Indoor Location Tracking-Withdrawn[END_REF]. Finally, the following step length model [START_REF] Zhou | Behaviors Classification based Distance Measuring System for Pedestrians via a Foot Mounted Inertial Sensor[END_REF] is chosen

     d W,k = 2 j=0 K j × f -j/2 s k + K 3 × max ãk (walking) d R,k = K 4 × 3 ãk (running) (15)
where subscripts 'W' and 'R' represent the walking and running behaviors, respectively. To remove the influences of high-frequency noise from acceleration, the Hamming-window based, linear-phase finite impulse response (FIR) filter is introduced to generate the filtered acceleration ãk . ãk denotes the mean value of ãk ; f s k denotes step frequency during the k-th step period; K j for j = 0, 1,

are coefficients that need to be calibrated beforehand. Note that the calibration can be done online during the initialization process with the given distance estimated by ZUPT-KF.

The heading angle during one motion period can be obtained based on various MARG sensors data fusion strategies [START_REF] Valenti | Keeping a good attitude: A quaternion-based orientation filter for imus and margs[END_REF][START_REF] Tian | A Multi-Mode Dead Reckoning System for Pedestrian Tracking Using Smartphones[END_REF][START_REF] Stirling | Evaluation of a New Method of Heading Estimation for Pedestrian Dead Reckoning Using Shoe Mounted Sensors[END_REF]. However, the performances in real applications will be degraded due to the existence of gait related heading deviation which differs from individuals. For this reason, the heading should to be cautiously aligned to the moving direction before using ( 14) by

γ i = ψ i + δγ i (16)
where δγ i is the heading deviation angle between moving direction and heading; γ i and ψ i represent one's actual moving direction and heading angle, respectively. The compensation method for δγ i will be detailed in Section III.

Working Mode Switching

It needs to be clarified that (i) initialization and dead reckoning are two basic working modes during pedestrian navigation.

(ii) when the changing of motion behavior is detected, the initialization mode will be triggered for performing the ZUPT-aided MARG KF and calibrating the personal parameters.

(iii) while once the initialization is successfully completed, the dead reckoning 165 mode will be carried out with the calibrated parameters in (ii).

The working mode switching diagram is exhibited in Fig. 3. Therefore, the Based on that, the real-time state machine for multi-mode behaviors recognition is designed and developed to make decision on working mode. More details can be found in reference [START_REF] Zhou | Behaviors Classification based Distance Measuring System for Pedestrians via a Foot Mounted Inertial Sensor[END_REF].

Gait Related Heading Deviation Angle Compensation

Recalling the heading deviation in [START_REF] Meng | Self-contained pedestrian tracking during normal walking using an inertial/magnetic sensor module[END_REF], it has a strong relation with personal gait. Inversely, the deviation also reflects one's gait and potentially can be applied to walking rehabilitation training. As a matter of fact, personal gait has a very good repeatability thus such deviation can be estimated within aforementioned initialization period. Assuming the following state model through KF of δ γ as

δ γ i = δ γ i-1 + w 1,i (17)
Letting x II i = δ γ i δγ i , the state model can be reformulated as

x II i = Φ II i-1,i x II i-1 + w II i ( 18
)
where the state transition matrix Φ II i-1,i is

Φ II i-1,i =   1 0 ∆t 1  
The two dimensional process noise can be written as

w II i = w 1,i w 2,i
where w 1,i ∼ N (0, Q w1 ) describes the reliability of model [START_REF] Prateek | Modeling, Detecting, and Tracking Freezing of Gait in Parkinson Disease Using Inertial Sensors[END_REF];

w 2,i ∼ N (0, Q w2 )
is caused by discretization of first order deviation δ γ = ∂ (δγ)/∂t. The variance of process noise can be written as

Q II = diag {Q w1 , Q w2 }
Accordingly, the observation model is established as

z II i = H II i x II i + ε II i ( 19
)
where the design matrix is H II i = [ 0 1 ] and the observation vector is constructed as z II i = γ i -ψ i . ψ i denotes the heading angle obtained from ZUPT-KF, which is taken at the instant of hitting the ground. γ i can be approximately calculated with the position increment or the velocity vector by

γ p,i = arctan ∆p x,i ∆p y,i (20)
γ v,i = arctan v max x,i v max y,i (21)
where [ ∆p x,i ∆p y,i ] are the two-dimensional position increment vector dur-175 ing the i-th step, [ v max x,i v max y,i ] are the maximal value of velocity during the i-th step. ε II i is the zero mean Guassian noise, and its variance R II can be determined after obtaining the variances of heading angle and the position (or velocity). Then the deviation angle can be recursively estimated by implementing the proposed KF.

The Adaptive Energy Management Strategy

Three aspects of adaptive energy management (mentioned in Fig. 3) are devised including energy-aware for gyroscope, mode switching mechanism and adaptive system de-sampling.

Energy-Aware Strategy for Gyroscopes Measurements Acquisition

As is pointed out in [START_REF] Makni | Energy-Aware Adaptive Attitude Estimation under External Acceleration for Pedestrian Navigation[END_REF], gyroscopes consume much more energy compared with that of accelerometers or magnetometers. To ensure the lower energy consumption, an energy-aware strategy for gyroscope measurements acquisition is given in Table 1. Behavior recognition method is utilized to identify one's motion behavior and non-static behavior will wake up the gyroscope from its sleeping status.

For these non-static behaviors, the gait based dead reckoning (Gait-DR) mode does not employ any gyroscope outputs while ZUPT-KF utilizes gyroscope measurements. Furthermore, one thing should be noticed that motion behavior changing may lead to corruptions of one's step length model. To avoid such problem, without significantly increasing the computations, a periodical step length model examination procedure is introduced with a long time-interval.

Mode Switching Mechanism of Navigation Computation

Due to different computation complexity for Gait-DR and ZUPT-KF, modes switching mechanism of navigation computation mode is designed to minimize energy consumption (see Table 2). Implement the following verification condition:

IF LZUPT-KF < L0
⇒ keep computation mode 'ZUPT-KF'.

ELSE ⇒ switch the computation mode 'ZUPT-KF' to 'Gait-DR'.

Gait-DR mode:

Implement the following verification condition: At the initialization stage, 'ZUPT-KF' calibrates the step length model and estimates the heading deviation angle till the accumulated moving distance L ZUPT-KF derived from ZUPT-KF does not exceed the predetermined threshold L 0 . In another word, if L ZUPT-KF exceeds L 0 , then the initialization is complete and the computation mode 'ZUPT-KF' will be switched to 'Gait-DR'.

IF
On the other hand, there are two conditions to trigger the switching 'Gait-DR' to 'ZUPT-KF': i) a new motion behavior is detected; ii) for a certain behavior, the motion frequency changes may lead to the original calibrated step length model inaccurate. For example, for running behavior, jogging and running fast are different, then one set of step length model parameters may not be incompatible with these two scenarios. Therefore, even one's motion behavior is kept, we still need to carefully monitor the possible corruption of step length model.

The details of step length model examination is presented in Table 3.

σ d = n i=1 (dKF,i-dSL,i) 2 n-1
dKF,i denotes the distance derived from ZUPT-KF at i-th step;

dSL,i denotes the distance with calibrated step length model of i-th step

Step length model examination of 'Gait-DR':

do the examination on dKF and dSL with a regular time interval T0:

by the following testing condition:

IF |dKF -dSL| ≤ 3σd
⇒ keep the mode 'Gait-DR'.

ELSE ⇒ switch the mode 'Gait-DR' to 'Initialization'

Adaptive De-sampling for Gait-DR Mode

To further reduce system power consumption, an adaptive de-sampling strategy is proposed when one is with 'Gait-PDR' mode. As is widely acknowledged that, the strapdown-inertial-mechanization based method, e.g. ZUPT-KF requires a high sampling rate of inertial measurements. This is because it essentially calculates the real-time PVA information by continuous integral computation, where a lower sampling rate cannot adequately capture one's motion thus inevitably degrading the PVA accuracy. In contrast, the dead reckoning based methods calculate the position only through heading and step length of each step which only requires the sampling rate satisfying Nyquist theorem and obviously does not need a high sampling rate. Inspired by this point of view, system de-sampling is feasible to perform the pedestrian navigation during the Gait-DR working period.

Experimental Results and Analyses

In this section, pedestrian experiments are carried out to verify the proposed method and fully evaluate its performances. The MARG module HI229, includes accelerometers, gyroscopes and magnetometers [50]. The hardware device is mounted on one of user's foot. The module is communicated with a memory card through TTL-USB module (Fig. 5). The computation and evaluation of proposed method is post-processed under Matlab r2015a on a Dell Inspiron14 laptop with a CPU of 4-core i5-4210U. during each detected step, such as [START_REF] Jimenez | A Comparison of Pedestrian Dead-Reckoning Algorithms Using a Low-Cost MEMS IMU[END_REF][START_REF] Tian | A Multi-Mode Dead Reckoning System for Pedestrian Tracking Using Smartphones[END_REF][START_REF] Lan | Using Smart-Phones and Floor Plans for Indoor Location Tracking-Withdrawn[END_REF][START_REF] Ho | Step-Detection and Adaptive Step-Length Estimation for Pedestrian Dead-Reckoning at Various Walking Speeds Using a Smartphone[END_REF] (3) The proposed online calibrated and heading deviation corrected Gait-PDR method. To verify the accuracy of proposed method, the user is arranged to walk along the corridor of report hall twice with an arbitrary pace. For the first circular moving, ZUPT-KF is used for online calibrating step length parameter and estimating heading deviation angle which is estimated according to Section III. Then during the second circular trajactory, the working mode is switched to Gait-PDR. Fig. 7 clearly shows the heading deviation between moving direction and heading angle. γ p and γ v in Fig. 7 are derived from Eqs. ( 16) and ( 17)

respectively. γ p shows a better result. As a consequence, ( 16) is chosen in KF II to estimate heading deviation angle and the result is given in Fig. 8.

Besides, four testers' heading deviation angle estimation results are computed and presented in Table 4. Seen from this table, the gait related deviation angles are indeed varied with different persons. The largest deviation angle among testers even reaches 30 degrees. Consequently, its trajectory estimation error is much larger than others with smaller deviation angles. For online calibration of step length model, we can point out that: to gurantee an reliable calibration process, the step frequency is used as a key indicator to discriminate different walking patterns (usually two or three patterns, i.e. walking slowly and walking fast, are adequate for describing one's walking modes). This is clearly verified through previous experiments (see Fig. 9). The threshold of step frequency variation can be empirically set as ±0.2 according to our previous tests. Similarly, calibration in running mode also can be done in the way of calibration implemented in walking mode.

Next, we start to evaluate and compare the presented schemes, i.e. ZUPT- KF, PDR and Gait-PDR with the data of second cycle trajectory (shown in Fig. 10). The reference trajectory is obtained with a distance-measuring tape.

The default sampling rate is set as 200Hz. With this high sampling rate, both ZUPT-KF and Gait-PDR methods perform well, while PDR method is inferior to ZUPT-KF and Gait-PDR suffering from the heading deviation angle. Comparisons of de-sampling cases are shown in Fig. 11. Result shows that Gait-PDR method maintains a good performance even at a lower sampling rate. In contrast, ZUPT-KF method gradually works worse even diverges as the sampling rate decreases as a result of sparse inertial integral computations. We also More detailed errors evaluation in walking and running experiments are listed in Table 5. It needs to be clarified that the average time consumption criterion of one step is computed based on the number of total steps and their time consumptions. The relative error of distance is defined as

σ d = dtol -d ref d ref × 100% (22)
where d ref is the reference distance; dtol is the total distance derived from ZUPT-KF or Gait-PDR methods. The position error of starting point is obtained according to the closure trajectory. Besides, the position accuracy and maximum heading deviation are also given for comparison purpose. Since the 'groundtruth' position of each step cannot be acquired, the position uncertainty is evaluated by feature points consisting of eight corner points #Pk (k = 1, 2, . . . , 8)

and the end point #P9. The positions of these points (as shown in Fig. 10) are utilized as references to compute the root mean square error (RMSE) of horizontal position through where N = 9 denotes the number of referenced (subscript 'ref') feature points.

RMSE (p) = 1 N N k=1 p ref,x,k -px,k 2 + p ref,y,k -py,k 2 (23)
These results statistics are fully consistent with the advantages and disadvantages of listed methods. Furthermore, for the aspect of energy saving, Table 5 respectively). Table 6 presents the gait deviation influences imposing on position accuracies of four testers under a given testing scenario. Once again, the results show that the deviation angle has non-negligible impacts on pedestrian dead reckoning and our deviation angle compensation method is effective in improving the position accuracy. The total distance error is given in Table 7.

It intuitively shows the influence of de-sampling on two methods. Thus conventional ZUPT-KF is not suitable for the platform with low sampling rate.

One the contrary, Gait-PDR is very computationally efficient compared with ZUPT-KF. Meanwhile it still maintains a good accuracy even in the scenarios of low sampling rate. It is therefore very promising for its calibration convenience (online), computation efficiency and high accuracy in future wearable consumer electronics devices.

To verify the mode switching mechanism and gyroscopes energy-ware strategy, the users are allowed to walk or run in their own paces to collect a new dataset. Behaviors identification and working mode switching results are shown in Fig. 15. It shows gyroscope is only employed in initialization process and also used to examine step length model in Gait-PDR mode. It is switched off for energy saving purpose in the rest time.

There are three key factors related to power consumption, i.e. communication, computational complexity and execution time. Usually power consumption of system can hardly be accurately quantized with the first factor. It is universally known that communication cost is highly related to sampling rate of hardware system. Therefore, qualitatively speaking, our method with a lower sampling rate inevitably consumes much less power than ZUPT-KF method.

The second factor is also a common indicator to reflect the system power consumption. To clearly show the computational complexity difference between two methods, an example is given as follows:

Given matrix A with dimension of m × n and B with dimension of n × p, the computation complexity of A × B is defined as

T (A m×n × B n×p ) = O(m × n × p) (24)
Table 8 presents complexity comparison between ZUPT-KF and PDR during one step walking period where zero velocity lasts about 0.3 seconds. Actually, after successfully initializing, our proposed method reduces to PDR method.

It is obviously superior to ZUPT-KF method in aspect of computation complexity. Next we focus on the computation of initialization which contains two main parts: step length model calibration and heading deviation correction.

The whole computation complexity for one time of initialization in our method consumes O(88n), where n denotes the counting number of steps during the Intuitively, the last factor has the most straightforward connection with power consumption. The average time consumption of one step is presented in Table 9. Owing to the mode switching mechanism and de-sampling strategy, the consumed time of proposed method in one step is reduced by 80.75%. Therefore, considering overall facts above, we cautiously come to a conclusion that our proposed method outperforms the ZUPT-KF method for power consumption.

To further verify the proposed online calibrated step length model, the representative step length models with offline calibration are compared with our step length model with online calibrated parameters [START_REF] Levi | Dead Reckoning Navigational System Using Accelerometer to Measure Foot Impacts[END_REF][START_REF] Ladetto | On Foot Navigation: Continuous Step Calibration Using both Complementary Recursive Prediction and Adaptive Kalman Filtering[END_REF][START_REF] Woodman | Pedestrian Localisation for Indoor Environments[END_REF][START_REF] Krach | Cascaded Estimation Architecture for Integration of Foot-Mounted Inertial Sensors[END_REF]. The estimated step length results are presented in Fig. 16. It is obvious that the proposed method is much more accurate and reliable than others especially when undergoing transitions of different behaviors. Proposed method 20

Conclusion

In this paper, an online calibrated, energy-aware and ambulatory gait corrected pedestrian dead reckoning (Gait-PDR) method is proposed for pedestrian navigation via the foot-mounted magnetic, angular rate and gravity (MARG) for low sampling rate scenarios. Overall, it is accurate, energy-saving and easyto-implement thus is very promising for pedestrian navigation applications. In addition, the personal deviation angle obtained during the pedestrian navigation process would be an essential indicator for evaluating the efficiency of patients' walking rehabilitation training. In the future, more complex pedestrian behavior modes e.g. going upstairs, downstairs and climbing will be extended in theory.

Appendix A

Due to the existence of sensor error, there is a small difference between 'true' attitude matrix C n b and the estimated one C n b (n and n represent the 'true' navigation frame and calculated navigation frame). According to the chain rule of coordinate frame,

C n b = C n n C n b ( 25
)
where C n n is defined as misalignment attitude matrix and it is very close to identity matrix. Therefore it can be written as:

C n n = I 3 + [φ×] =      1 -φ z φ y φ z 1 -φ x -φ y φ x 1      (26
s ψ -c ψ 0           φ x φ y φ z      (29)
Thereby the heading errors is written as

δψ = Dφ (30)
where D = sinψtanθ -cosψtanθ 1 is only related to heading and pitch angles.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: The diagram of proposed pedestrian navigation method

y

  denotes the output of sensor; w G and w A are the zero mean Guassian white noises with variances Q wG and Q wA respectively; b G and b A represent the gyroscope and accelerometer biases respectively. The models of b G and b A are formulated as ḃG = w rG (4)

  145

Figure 4 :

 4 Figure 4: Diagram of working mode switching

  170

  180

13

 13

  new motion behavior is detected ⇒ switch the computation mode to 'ZUPT-KF'. ELSE IF the step length model cannot adapt to the current motion ⇒ switch the computation mode 'Gait-DR' to 'ZUPT-KF' ELSE ⇒ keep 'Gait-DR' computation mode.

Figure 5 :

 5 Figure 5: HI229 and TTL-USB module

Figure 6 :

 6 Figure 6: 2D trajectory of pedestrian experiment

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Comparison of moving direction and heading (walking)

Figure 10 :

 10 Figure 10: Trajectory comparisons of different methods (walking with the sampling rate of 200Hz)

Figure 11 :

 11 Figure 11: Results of de-sampling cases (walking)

Figure 12 :

 12 Figure 12: Gait related heading deviation angle estimation results (running)

Figure 13 :

 13 Figure 13: Trajectory comparisons of different methods (running with the sampling rate of 200Hz)

Figure 14 :

 14 Figure 14: Results of de-sampling cases (running)

Figure 15 :

 15 Figure 15: Illustrations of behaviors identification, working mode switching and gyroscope mode switching.

Figure 16 :

 16 Figure 16: Step length estimation results in the scenario of behaviors transition

c

  ) in which φ is referred to as misalignment angle. Herein, the rotation matrix C b n can be expressed a set of Euler anglesβ c ψ -s β s θ s ψ s ψ c β + c ψ s β s θ -s β c θ -c θ s ψ c θ c ψ s θ s β c ψ + c β s θ s ψ s ψ s β -c β s θ c ψ c θ cβ angles are with the rotation order of '3-1-2', i.e. yaw-pitch-roll ( ψ -θ -β ); c and s denote cosine and sine functions respectively. Denoting the calculated Euler angles as ψ θ β , and then the Euler angle errors are ψ s θ -c ψ s θ c θ -c ψ c θ -s ψ c θ 0

Table 1 :

 1 Energy-aware strategy for gyroscope measurements acquisition

	Behavior recognition:

IF

recognition result is non-static (e.g. "walking" or "running")

⇒

turns to Working mode identification ELSE ⇒ gyroscope 'off' Working mode identification: IF working mode is "ZUPT-KF" ⇒ gyroscope 'on' ELSE ⇒ turns to Step length model examination Step length model examination: IF ∆TG ≥ T0 (∆TG: the time interval from time instant of last step length model examination; T0: the time duration threshold of step length model examination) ⇒ gyroscope 'on' till a whole step is completed ELSE ⇒ gyroscope 'off'

Table 2 :

 2 Mode switching mechanism of navigation computation

ZUPT-KF mode:

Table 3 :

 3 Details of step length model examination

	Priori precision statistics of step length model in 'Initialization' part:

Table 4 :

 4 Gait related deviation angles of four testers

			Age	Height Weight Deviation angle (deg)
	Testers	Gender
			(years)	(cm)	(kg)	Range	Average
	Tester A	Male	24	175	65	27.5 ∼ 30.3	29.2
	Tester B	Male	25	180	72	5.1∼7.1	5.8
	Tester C Female	24	158	44	15.2∼17.7	15.6
	Tester D Female	23	160	51	7.3∼9.5	8.7

Table 5 :

 5 Statistics results for pedestrian navigation accuracy and time consumption

	Working	Sampling Average time consumption	Distance	Misclosure RMSE (p) Max orientation
	Behavior
	mode	(Hz)	for one step (×10 -3 s)	relative error	(m)	(m)	error (deg)
		200	100	0.57%	0.10	0.46	0.9
		100	48	1.10%	0.45	0.59	1.2
	ZUPT-KF	75	36	0.80%	0.51	0.89	3.0
		70	33	0.57%	5.31	3.12	3.1
		50	24	5.30%	11.60	8.95	7.1
	Walking
		200	23	0.57%	1.13	1.50	1.1
		100	12	0.59%	1.51	1.28	1.3
	Gait-PDR	75	6	0.59%	1.11	1.42	1.6
		50	4	0.57%	1.56	1.21	1.6
		20	2	0.58%	1.61	1.39	2.0
		200	78	1.46%	0.51	0.65	1.4
		100	42	3.06%	4.37	2.67	9.9
	ZUPT-KF
		80	35	3.18%	10.35	5.65	10.1
	Running	200	16	2.01%	1.16	1.73	1.4
	Gait-PDR	100	7	2.64%	2.21	1.78	1.6
		80	6	3.46%	2.00	1.87	1.6

initialization. Supposing n as 40, the computation complexity of one time of initialization is O(3520). It should be pointed out that the times of initialization is limited and usually does not exceed 3 ∼ 5 considering one's typical motion be-340 havior modes. Therefore, likewise, our proposed method has a power-consuming advantage on computation complexity.

Table 6 :

 6 Gait deviation influences imposing on position accuracies of four testers with sampling

	rate of 100Hz (unit: m)
	Testers	RMSE(p): PDR RMSE(p): Gait-PDR
	Tester A	20.34	1.34
	Tester B	4.10	1.28
	Tester C	10.26	1.15
	Tester D	7.28	1.30

Table 7 :

 7 Total distance error (walking)

	Working mode Sampling (Hz) Total distance error (m)
	200	1.00
	ZUPT-KF
	50	9.29
	200	1.00
	Gait-PDR
	50	1.01

Table 8 :

 8 Complexity comparison between ZUPT-KF and PDR during walking period of 1.2

	seconds
	Method	Sampling rate (Hz) Computation cycles	Complexity
	ZUPT-KF	200	60	O(1.0377 × 10 6 )
	PDR	50	1	O(6)

Table 9 :

 9 Average time consumption for one step

	Method	Time Consumption (×10 -3 s)
	ZUPT-KF (200Hz)	104

Acknowledgements

This work was supported by the National Natural Science Foundation of China (42074038) and the Department of Science and Technology of Sichuan Province (2021YJ0104).


domain: spi phys

text:
Valentin Ritzenthaler
email: valentin.ritzenthaler@onera.fr

Pierre Cantin

Xavier Ferrières

Sébastien Pernet

Guillaume Puigt

Spectral Difference method On Structured-Grids for Maxwell's Equations in Time Domain

Keywords: High-order method, Spectral Difference, Maxwell's equations, Time domain

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The Spectral Difference (SD) method is a discontinuous high-order method based on local polynomial interpolations. This method was introduced in [START_REF] Liu | Spectral Difference method for unstructured grids I: Basic formulation[END_REF] as a scheme that is conservative, high-order, geometrically flexible, computationally efficient and simply formulated. It solves the strong form of the equation, as in Finite Difference. This method has been widely explored in Computational Fluid Dynamic (CFD) as an alternative to the Discontinuous Galerkin method (DG) [START_REF] Yu | On the accuracy and efficiency of discontinuous Galerkin, Spectral Difference and correction procedure via reconstruction methods[END_REF]. The present article starts from the formulation in [START_REF] Liu | Spectral Difference method for unstructured grids I: Basic formulation[END_REF] for Maxwell's equations but accounts for specific location of the degrees of freedom associated with stable SD schemes [START_REF] Van Den Abeele | On the stability and accuracy of the Spectral Difference method[END_REF]. This paper is organized as follows. Section 1 explains how the SD method works for a 1D conservation law and how it generalizes to Maxwell's equations. Section 2 shows a testcase, on a non-Cartesian grid, with well-known analytical solution.

Spectral Difference method

In 1D, the main idea of the SD method, for a conservative system (∂ t u = -∂ x F (u)), is to compute the conserved variables u as a polynomial of degree p (p ∈ N * ) and to compute the flux as a polynomial of degree p + 1 so that solution and flux divergence are both polynomials of degree p. Solution and flux polynomials are defined by Lagrange interpolation using two sets of points, p + 1 solution points and p + 2 flux points (Figure 1). The stability of the method only depends on the location of the latter [START_REF] Van Den Abeele | On the stability and accuracy of the Spectral Difference method[END_REF]. In 1D, two of the flux points must be taken as the segment end points, the p remaining are taken inside the element. The formulation is shown linearly stable for any degree p [START_REF] Jameson | A proof of the stability of the Spectral Difference method for all orders of accuracy[END_REF]. The SD algorithm consists of four steps:

x k x k+1 • • •
1. The solution is extrapolated to the flux points. 2. For the two boundary points, the flux is the solution of a Riemann problem. On the internal points, the flux is a linear combination of the unknowns. 3. The flux polynomial is interpolated from the flux points. 4. The flux is differentiated at the solution points in order to update the solutions.

For 2D/3D configurations, problems are considered using tensorized 1D formulation. An efficient implementation can easily be obtained using matrix/vector products.

Suggested members of the Scientific Committee:

Bruno Després, Julien Diaz.

The SD formulation is applied to Maxwell's equations written in conservative form; they consist of two coupled linear advection equations (in 1D, u = (E, H) and F (u) = (H, E)). From there, everything described above applies: each component of the fields E and H is approximated by a polynomial of degree p, and the interpolations are taken in every directions using a tensorial rule.

Numerical results

As a primary validation, the method is tested by considering a propagating mode (here the (1, 1, 0) mode) inside the unit cube of R 3 . For computation, p = 3. The mesh is a structured non-Cartesian grid (Figure 2), with 32,768 degrees of freedom and 512 cells. Figure 3 shows the good agreement between the SD numerical solution and the analytical solution. Other configurations will be presented during the conference.

Conclusion

In this article, a new way (the SD method) to discretize Maxwell's equations is introduced and validated for a cavity mode by using a non-Cartesian grid. For the presentation further details will be given on the formulation of the SD method, in particular stability results and comparisons with both FDTD [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF] and DG schemes [START_REF] Cohen | A spatial high-order hexahedral discontinuous Galerkin method to solve Maxwell's equations in time domain[END_REF].

Figure 1 :

 1 Figure 1: Example of point positions over a segment, for p = 2: solution points (•) and flux points ( ).

Figure 2 :

 2 Figure 2: Structured non-Cartesian mesh used for the simulation. 32,768 degrees of freedom, 512 cells.

Figure 3 :

 3 Figure 3: Comparison of the numerical solution given by the SD method for p = 3 (red curve), and the analytical solution (blue dots) for the mode (1, 1, 0).