Ezzoubair Bendadesse

Charlotte Gervillié-Mouravieff

Cédric Leau

Kateryna Goloviznina

François Rabuel

Mathieu Salanne

Jean-Marie Tarascon
email: jean.marie-tarascon@college-de-france.fr

Ozlem Sel
email: ozlem.sel@college-de-france.fr

Spotting interface structuring during Na-insertion into the NaSICON Na3V2(PO4)3 by EQCM and operando fiber optic infrared spectroscopy

Keywords: Sodium ion battery, electrode/electrolyte interface, EQCM, Infrared Fiber Wave Evanescent Spectroscopy, Molecular Dynamics simulations 1

Reliable operation of Na-ion batteries requires a mastering of the interfacial processes governing intercalation kinetics. Ion desolvation, amongst others, is a crucial phenomenon as far as power applications are envisioned. Studying the intricacy of interfacial processes faces experimental challenges. Here we report the complementary use of EQCM and operando Infrared Fiber Wave Evanescent Spectroscopy (IR-FEWS) combined with MD simulations to investigate the electrode/electrolyte interface of polyanionic Na-insertion compounds Na3V2(PO4)3 (NVP) and Na3V2(PO4)2F3 (NVPF). Through screening of various electrolytes that differ in their solvation energy and viscosity, we demonstrate the existence of an unusual mass uptake/depletion signal at the NVP interfaces only, indicative of a dissimilarity in migration kinetics of solvent molecules between NVP and NVPF interfaces as supported by IR-FEWS. We rationalize our finding that sodium depletion is caused by the unusually rapid ion diffusion within NVP mediated by an intermediate metastable phase, which leads to opposite Na + concentration and free solvent molecule gradients at the interface that causes the unusual and anomalous feature observed on EQCM frequency profiles. This feature does not appear with NVPF showing lower ionic diffusion. Altogether, this fundamental study provides valuable fundamental insights on the interface structuring of insertion compounds depending upon their Na + mobility.

Introduction

Sodium ion batteries (SIBs) are an emerging alternative to the Li-ion technology. Their affordability coupled with a greater abundance of sodium containing precursors offer a nice trade-off in terms of cost vs. energy density. [1][2][3][4][5] Common to all alkali metal batteries, interfacial phenomena are of paramount importance. Notably, the processes governing the Electrode Electrolyte Interface (EEI) are of particular interest to the battery research community. Alike Li-ion batteries, an often unclear and elusive picture is painted when discussing interfacial phenomena in SIBs due to the limited number of non-invasive analytical techniques. [6] Moreover, amongst these phenomena is the charge transfer process leading to ion intercalation in the host structure, often evoked in discussions when it comes to high performances due to its important contribution to the kinetics of insertion. [7][8][9] Indeed, through past works, researchers have been able to distinguish the rate-limiting step [10][11][12] among the different processes constituting the charge transfer which are known to be: cation migration from the bulk of the electrolyte to the interface, desolvation, diffusion through possible interphase and finally through the bulk of the material to an electrochemically active site. However, to the best of our knowledge, no substantial correlation has been established linking the nature of the material with the processes leading to charge transfer at the interface. The solvation/desolvation process, known to be a key aspect of electrolyte engineering, was investigated using spectroscopy techniques consisting of Surface Enhanced Raman Spectroscopy (SERS) [13] or operando infrared spectroscopy. [14] More recently, our group reported the use of fiber optic evanescent wave sensors (FOEWS) for carrying operando infrared spectroscopy to unravel the dynamics of (de)solvation amongst other interface phenomena (Scheme S1). [15] Additionally, Electrochemical Quartz Crystal Microbalance (EQCM) based methods has become a tool of choice to quantify the solvation shell of alkali metal ions and its evolution through charge and discharge as well as to probe different salt/solvent combinations as it is illustrated in past works by our group and others. [16][17][18][19][20] Herein, we examine a NaSICON material Na3V2(PO4)3 (NVP) to find any evidences for a close connection between the high Na + mobility of the NaSICON structure and the desolvation process prior to sodium intercalation. To do so, we employ EQCM, operando IR fibre evanescent wave spectroscopy (IR-FEWS) and MD simulations as complementary tools to investigate the dynamics of the (de)solvation process in different electrolyte solvents. More specifically, we explore the importance of electrode material/electrolyte selection on the sodium ion desolvation process and on rate capabilities of the active material targeted.

Results and Discussion

The unusual frequency variation profile of Na3V2(PO4)3. The NVP spray-coated quartz resonators were cycled in three different electrolytes consisting of NaPF6 as a salt and PC, EC:DMC (1:1 mol) and DMC as a solvent, bearing in mind in light of our previous work that the dielectric constant is a varying parameter that affects the EQCM response of a Li + insertion compound. [16] Figure 1 a-c shows the frequency response of the NVP thin film in the three electrolytes during cyclic voltammetry between 2.8 and 3.6 V vs. Na/Na + (Figure S1a-c).

Capacities calculated from current vs time plots were found to be 70 mAh g -1 and 40 mAh g -1 for NVP and NVPF, respectively, at a corresponding C-rate of 3C (1 mV s -1 ). We can easily distinguish a cationic response in the form of a frequency increase/decrease during desodiation/sodiation, respectively. However, a peculiar feature of this response is the small frequency "bump" or peak that repeats upon cycling and that seems to gain in amplitude with a lower solvent polarity (from PC, EC:DMC to DMC). The peak appears exactly before the start of frequency decrease through a negative scan, which means that we have a punctual rapid mass loss followed by a mass increase due to sodiation (Figure 1 a-c, insets). We can speculate that we are spotting the Na ion-desolvation process, which involves the return of the solvent molecules from the coordination shells of the Na ions to the electrolyte bulk, based on the fraction of mass loss that changes inversely proportionate to the solvent polarity and viscosity.

The same trend was also recorded upon inserting lithium ions into the NVP structure (Figure S2a-c) and during galvanostatic charge discharge vs. Na/Na + (Figure S2d-f). Additionally, we also observe the reverse process of ion solvation at the end of the charge/oxidation with a corresponding frequency decrease (mass gain) due to solvent molecules diffusing back to the solvent-poor interface to resolvate the Na + cations (Figure S2d). The variation of the frequency gain (caused by desolvation mass loss) as a function of the solvent's polarity is shown in Figure S3. For both lithium and sodium insertion, a drop in solvent polarity is associated with an increase in the % frequency gain linked to desolvation. This ties nicely with the solvation energies calculated for Na + bound to the carbonate solvents of interest, where weaker interactions were found linking Na + and linear carbonates such as DMC than cyclic carbonates such as PC and EC. [21] With weaker solvation energies, DMC molecules detach more readily than its cyclic counterpart during sodium insertion into the NVP structure, and thanks to its lower viscosity [22] it diffuses back faster to the electrolyte bulk causing the mass loss anomaly upon reduction. This further supports our belief that the brief mass loss/frequency gain peak is the result of sodium ion desolvation. For the sake of comparison, a non-NASICON polyanionic compound in the form of Na3V2(PO4)2F3 (NVPF) was then evaluated alongside with NVP (Figure 1 d-f).

f.

As the redox process of NVPF (de)sodiation occurs at higher potentials, we limit the charge cutoff potential to 3.9V vs. Na/Na + to prevent any electrolyte degradation from having an impact on our frequency response, especially in the DMC containing electrolytes. The EQCM response was thus recorded for the first sodium removal/intake happening at ~3.6V vs. Na/Na + (~3.35 V vs. Na/Na + for NVP) (Figure S1 d-f). When cycled in the same electrolytes as NVP, the NVPF thin film exhibited the expected cationic response of (de)sodiation. None of the three CV-EQCM curves showed the mass loss peak attributed to desolvation (Figure 1 d-f). The drastic difference in the mass variation behavior of both materials leads us to consider the influence of the materials' properties on the desolvation process. The substitution of one (PO4) 3- with three F -, owing to the inductive effect, leads to an increase in the average potential of NVPF when compared to NVP and hence the energy density of the material. [23] However, it has been proven that Na + diffusion kinetics are faster in NVP than NVPF. For instance, Li et al.

have shown that the NVP phase dominated the insertion kinetics at higher rates in a NVP /NVPF blended cathode. [24,25] Operando IR fibre evanescent wave spectroscopy at the service of interface monitoring.

To gain further insights into the chemical origin of the punctual mass loss peak detected in EQCM analyses (Figure 1a-c), we ran operando IR fiber evanescent wave spectroscopy (IR-FEWS) measurement on the NVP system. Indeed, our group developed the possibility of applying such technique operando conditions using optical fibers as propagation media for IR wavelength. [15] In short, absorption spectroscopy through evanescent waves is performed at the surface of the Te2As3Se5 (TAS) fiber that is embedded in the active material in contact with the electrolyte. The penetration depth of the evanescent waves depends on the refractive index of the surrounding media (~1 µm for the electrolytes of this study) and exceeds that of EQCM acoustic wave (range from ~140-270 nm for the electrolytes of this study). More details about the IR-FEWS setup are in the Supplementary Information file (Scheme S1).

NVP was cycled in two different electrolytes (NaPF6 in PC and DMC 1%vol FEC) to highlight the difference in the dynamics of desolvation following the drastic change of solvent polarity and viscosity. Coordinated and free solvent bands were scrutinized with a special focus on the υOCO band situated between 1150 and 1400 cm -1 for both solvents as the solvation of Na + ions by DMC molecules resulted in a weakening of the 𝜈C=O vibrational band of DMC and a strengthening of the 𝜈O-C-O band (Figures 2a and3a). For clarity, only the 2 nd and 3 rd cycles were analyzed to avoid any electrolyte degradation related shifts of bands due to possible CEI formation. The evolution of 𝜈 OCO absorbance and its relative absorbance evolution after The υOCO band shift during cycling was monitored in DMC (Figure 2d) and PC (Figure 3d) containing electrolytes. We can clearly notice a decrease and an increase of the band absorbance of coordinated and free solvent molecules, respectively, following discharge (vice-versa during charge). This evolution depicts the corresponding desolvation and solvation processes occurring during Na + insertion and disinsertion. The ratio of coordinated and free band absorbances is plotted versus time to portray the dynamics of the separation of sodium ions from its solvation shell. A clear difference is spotted between both electrolytes as the DMC υOCO band ratio shows a rapid decay (increase) at the beginning of discharge (charge) which can be correlated to the fast depletion of free solvent molecules at the interface caused by their re-migration towards the bulk of the electrolyte upon sodium ion insertion. Meanwhile, a more linear evolution is observed in the PC electrolyte throughout both charge and discharge processes highlighting the slow diffusion of solvent molecules possibly due to higher solvent viscosity as it will be discussed later. At this stage, a legitimate question regards how such high sodium diffusion kinetics in NVP could account for the desolvation process observed by EQCM and IR-FEWS alike. While operando techniques are unable to replicate the exact real-life conditions of sodium-ion batteries (SIBs) operation due to equipment or technique limitations, the impact of NVP's intrinsic transport properties and the observed anomaly during the desolvation process should still be valid in actual battery conditions. This is because these processes occur at the nanoscale, suggesting that the fundamental effects would persist despite the imperfect replication of real-life conditions by operando techniques.

Estimating ionic the apparent diffusion coefficients. To test the aforementioned hypothesis, we estimated the apparent solid diffusion coefficient of Na + in both materials using two different techniques: Randles-Sevcik equation and Potentiostatic Intermittent Titration Technique (PITT) measurements (Figure 4) following the protocols described in Supplementary information. Both techniques were applied to thin-film configuration electrodes sprayed onto the QCM. When using Randles-Sevcik, we notice that the estimated coefficient for sodium insertion and extraction in both structures (NaSICON and non-NaSICON) are comparable, with a slightly higher value of DNa + during discharge in NVP and charge in NVPF.

Figure 4a

and 4b show that DNa + values, obtained by PITT analysis for both materials are at least one order of magnitude higher than DNa + estimated by Randles-Sevcik, this can be explained by its dynamic nature which takes into account several charge transfer steps enlisted in the Na + insertion process which can lead to an underestimation of the apparent diffusion coefficient. Differently, PITT analysis measures current relaxation during a potentiostatic step of a controlled duration under equilibrium conditions, thus, the estimated diffusion coefficient should reflect the kinetics of sodium diffusion in the crystal bulk. Furthermore, the sharp drop in DNa + observed during charge and discharge in the middle of the plateau for NVP is not an experimental artefact as it has been reproduced twice, but rather due to the stronger interactions between sodium ions and host structures during the two-phase transition reaction, as previously observed for LFP and NVP. [26,27] Nonetheless, whatever the technique used, it remains that the ion diffusion in NVP (~10 -9 cm 2 s -1 ) is higher than in NVPF (~10 -10 cm 2 s -1 ), comforting our hypothesis in the influence of fast sodium diffusion ability on the structuration of the EEI.

Correlating the EQCM frequency profile anomaly to superior rate capabilities of NVP.

As discussed earlier, NVP was proven to withstand high rates more efficiently than NVPF. [24,25] For a fair and direct comparison of high-rate properties between both materials, we performed power rate experiments over different rates (ranging from C/10 to 5C). We cycled both NVP and NVPF in the same electrolyte (EC:DMC 1%vol FEC) over the full potential windows for both materials (2-3.8V and 2-4.3V vs. Na/Na + for NVP and NVPF, respectively) (Figure 4c).

Spectacularly, NVP fully outpaces NVPF in terms of rate performances with namely a 65% of the initial capacity retention for NVP as compared to only 20% for NPVF at 5C. Such a drastic difference in the rate performance comes somewhat as a surprise bearing in mind the one order of magnitude difference, at most, in the DNa + between both materials. To further interrogate the underlying science beyond this difference, we limit the state of charge (SOC) to its half (1 exchanged Na + ) in order to closely match the EQCM testing conditions, so as to deduce meaningful correlations if any. Surprisingly, when Δx=1 (in Na3±xV2(PO4)3 and Na3±xV2(PO4)2F3), the NVP structure retains 100% of the initial capacity at very high rate (10C) while NVPF seems to display similar capacity fade to an unlimited SOC cutoff condition (Figure 4d and Figure S5c-d). This outstanding rate performance of NVP suggests that the first sodium ion extracted may suffers solely from minor kinetic limitations as opposed to the removal of the second Na + . Such an observation can find its meaning in light of recent work by Chotard et al. [28] that demonstrated, through operando synchrotron X-ray diffraction (SXRD), the existence of a metastable intermediate phase during the biphasic reaction of Na3V2(PO4)3 -Na1V2(PO4)3. Indeed, Na2V2(PO4)3 seems to appear in a more pronounced manner during nonequilibrium conditions (high C rates) and shows a reduction in symmetry from the usual Napoor/rich phases. This entails the splitting of Na(1) and Na(2) sites into several new positions which creates vacancies at the Na(2) sites. Additionally, due to a different connectivity to the Na(1) site that results from such site-splitting, a Na(2a) sodium site forms and may offer a faster diffusion path, explaining the cause of such high rate performances. Moreover, it was deduced that the intermediate phase considerably reduces the lattice mismatch and facilitates the Na3-Na1 phase transition. In contrast, the Na-driven phase transition diagram in NVPF is more complex and involves the emergence of many single phases at both SOC cutoffs (Δx=1 or 2).

Besides, the fact that NVPF does not have the suited NASICON structure for fast ionic diffusion, this multiple-phase mechanism during partial and full (de)sodiation could be an additional handicap to explain the inferior rate performances of NVPF when compared to NVP. [29] At this stage a legitimate question regards the importance of this metastable phase at Δx=1 Na + for NVP and its effect on observed EEI structuration. Interestingly, there have been reports mentioning the existence of a metastable phase occurring under non-equilibrium conditions [30] at Δx=0.5 Li + in the polyanionic LiFePO4 phase. We thus explore the rate capability of LiFePO4/Li half cells. Figure S6 demonstrates that the capacity recovered at a 5C discharge rate is only half that of the initial charge capacity. This suggests that the mere presence of a metastable phase is not enough to enhance the diffusion kinetics to such outstanding levels. Besides the intricacy of phase transitions, particle size is a legitimate concern in regards to rate capabilities. Although diffusion coefficients are an intrinsic property of each material, a reduced particle size can greatly enhance power rate performances by shortening the cation travelling distance as evidenced repeatedly for an array of materials. [31,32] In our case, NVPF particles were found to be smaller than NVP particles as evidenced by SEM images (Figure S7), hence suggesting a higher rate performance for NVPF; the opposite of what we observed. Lastly, it should be mentioned that all materials used in this study (NVP, LFP and NVPF) were carbon coated in-situ via solid state synthesis to counter the inherently poor electronic conduction of polyanionic materials. Altogether, these results seem to indicate that the unusual Na + uptake-removal process spotted by EQCM is inherent to NASICON-type NVP and its metastable phase. Work exploring this topic is ongoing in our group.

Impact of diffusion kinetics on rate limiting step.

To interrogate the effect of the enhanced diffusion kinetics on the rate-limiting step of each material, EQCM tests were conducted using CV at different sweeping rates (1 to 5 mV s -1 ) to monitor the frequency profile during different rates. Firstly, NVPF-loaded QCM resonator seems to display lower mass exchange amplitude (smaller Df) with higher rates (Figure 5a-b), meaning that the quantity of sodium intake is lower at higher rates which correlates nicely with the lower capacities obtained at higher Crates (Figure S5 c-d). Figure 5d show the frequency variation profile for NVP, surprisingly, the amplitude increases proportionally to the sweeping rate as opposed to NVPF. Higher Crates allow us to keep a larger portion of the initial capacity, but a capacity loss is still seen in NVP upon faster cycling (Figure S5a). We conclude that the desolvation process is the ratelimiting step during sodium (de)intercalation in NVP, as evidenced by the fact that the faster the sweeping rate is, the more mass appears to be transferred to the EEI (higher Df during cathodic scan, Figure 5d). This suggests that the desolvation process is hampered by the fast ion transfer kinetics to the interface. More specifically, faster sweeping rates cause the concentration gradients of bare Na + and uncoordinated solvent molecules to develop at the interface and to become more and more prominent, increasing polarization as seen when cycling at high C-rates (Figure S5a-b).

Confirmation of the desolvation by MD simulations.

In order to rationalize these experimental findings, we performed MD simulations of the 1 M NaPF6 in PC in the contact with the NVP electrode. In equilibrium condition, with no voltage applied, the closest layer is formed by PC molecules (Figure 6a) adopting two different orientation of carbonyl groups with respect to the electrode surface (the peaks at ∼3.1 Å and ∼5.9 Å). The second solvent layer includes not only the PC molecules (broad low-intense peak at ∼8.6 Å) but also the NaPF6 salt (peaks at 9.2 Å and 9.4 Å for Na + and P of PF6 -, respectively), as is illustrated by the simulation snapshot (Figure 6b). Such a structural arrangement in the equilibrium state creates a prerequisite for the formation of the concentration gradient during the interaction of the Na + ions into the NVP. The second effect traceable by MD simulations is the desolvation of sodium ions. In the bulk, a single Na + cation is surrounded by ∼0.95 anions and ∼4.75 solvent molecules, on average, in its first coordination shell (Figure S10b). On the contrary, while Each symbol corresponds to a separate snapshot of a trajectory during which a single Na + cation was pulled from the bulk to the NVP surface. approaching the NVP surface, its total coordination number continuously decreases from six at about 9-10 Å from the surface to one at about 2-3 Å (Figure 6c). The sodium ion is not able to preserve its solvation shell when diffusing through the dense solvent layer at the surface since an important energy barrier should be overcome. Therefore, the experimental findings are confirmed by the simulations.

Conclusion

We believe that thanks to the high sodium ion mobility in the bulk of NVP and the presence of the metastable intermediate phase upon sodium insertion, there is a depletion of bare (desolvated) Na + at close surface, creating a concentration gradient at the interface, which in turn induces an opposite concentration gradient of uncoordinated solvent molecules. This will drive the latter to diffuse back to the bulk of the electrolyte and thus be detected by the QCM resonator as a mass loss. As for the case of NVPF, the uncoordinated solvent molecules do not diffuse as far away from the interface as in NVP's case. The whole desolvation and diffusion process is fully confined within the penetration depth of the acoustic wave emanating from the QCM resonator and thus remains undetectable. The penetration depth, δn, is calculated to be 271, 198, 140 nm for PC, EC:DMC and DMC, respectively. [33] Such a scenario can be pictured via a qualitative model (Scheme 1) where opposing concentration gradients coexist at the interface are sketched together with the penetration depth variation superimposed. This explains how it is feasible to observe the mass loss from desolvating Na + cations before they are inserted into the NaSICON material. Moreover, the mass loss could potentially be attributed in part to the anions participating the charge shielding of the electrode's surface upon polarization and the resulting Electrical Double Layer (EDL).

But, this possibility is discarded as we have demonstrated in our recent work that anions populating the EDL were ejected during lithiation of LixMoO3 in a PC containing electrolyte but became less mobile when DMC was used as a solvent due to contact ion pair and aggregate formation. [16] Here, because of the pronounced concentration gradient at the interface and the lower viscosity of such a linear carbonate solvent, [22] free DMC molecules diffuse more quickly than PC molecules during cycling. This is in accordance with the IR-FEWS results that show a clear difference in solvent diffusion kinetics based on polarity and viscosity (Figure 2 and3), similar to the distinction between the two electrolytes' EQCM profiles, where we see a much more pronounced desolvation "peak" in DMC than in PC (Figure 1a-c andFigure S3). This is further confirmed by the MD simulations at equilibrium where a similar structuration of solvent molecules and cations is observed (Figure 6).

Based on these findings, one can strategize different approaches to deal with the shortcomings of each material by a careful selection of electrolyte/electrode combination. For instance, "softsolvating" electrolytes, as proposed by Wang et al [34] , could be considered to overcome the ratelimiting desolvation process in NVP's case and lower the concentration-driven polarization. On the other hand, NVPF would benefit from a more stable CEI to prevent capacity fade related to vanadium dissolution [35] issues by incorporating electrolytes that would favor the formation of an inorganic CEI such as Localized High Concentration Electrolytes (LHCE). [36,37] In conclusion, the desolvation process during the insertion of Na + in the NVP framework was captured thanks to EQCM and operando infrared spectroscopy measurements and confirmed by MD simulations. We infer that the origin behind the appearance of such punctual phenomenon is the high Na + diffusion kinetics properties of the NaSICON structured Na3V2(PO4)3 as the desolvation process creates opposing solvent and bare Na + concentration gradient, leading the solvent molecules to migrate in the opposite direction out of the probed boundary of the QCM resonator. As a comparison, NVPF, shown to have slower bulk diffusion kinetics, does not exhibit the mass loss peak attributed to desolvation from which we conclude that desolvation is confined within the penetration depth of the resonator's acoustic wave and becomes undetectable for classical EQCM measurements. Altogether, these findings provide a clearer picture of the sodium-ion interface structuring in polyanionic compounds that should enable to design better electrode/electrolyte combinations to enhance SIBs performances.

Associated Content

Experimental details of thin film elaboration (Na3V2(PO4)3, Na3V2(PO4)2F3 and LiFePO4) on quartz resonators, electrochemical and EQCM analyses, estimation of apparent solid diffusion coefficient of Na + by PITT and Randles-Sevcik and the operando IR-FEWS spectroscopy procedure, the MD simulation protocol and the force field details.

workstation coupled with a commercial SEIKO QCM922A microbalance which permitted the resonance frequency (f) along with the motional resistance (Rm) to be monitored during electrochemical cycling. Coated QCMs act as working electrode and metallic sodium acts as both reference and counter electrodes. The electrolyte concentration (1 mol L -1 ) was kept constant, using sodium hexafluorophosphate (NaPF6) in combination with 3 different solvents (propylene carbonate (PC), ethylene carbonate:dimethyl carbonate (EC:DMC) (1:1 mol) and dimethyl carbonate (DMC)) 1%vol addition of fluoroethylene carbonate (FEC). The cell was cycled using cycling voltammetry (CV) or galvanostatic cycling with potential limitation (GCPL) and the frequency response of the coated QCM was simultaneously monitored. The EQCM tests of NVP were also performed in 1M LiPF6 in DMC, EC:DMC and PC, for comparison purposes.

IR fibre evanescent wave spectroscopy (IR-FEWS): A 19 mm-diameter Swagelok cell's body was drilled from side to side, and the electrode slurry was applied using a stainless steel current collector (Figure S1a). The holes in the cell body were lined up with the positions of the stainless steel plunger and current collector so that the fiber could pass through the hole above the current collector. The fiber was once connected to the spectrometer to gather background information after being fixed on both sides with cured epoxy for 12 hours. The surface of the fiber and the current collector were then covered with a droplet of the electrode formulation.

Finally, the entire assembly was dried in a vacuum oven overnight at 80 °C. The cells were put together in a glove box using two glass fiber separators (GF/D, Whatman) and a disc of sodium metal.

The operando measurements were performed with a Fourier-transformed infrared spectrometer (T37 or Invenio S, Bruker) with an accessory connection system on the right side to focus the infrared beam on one entry extremity of the fibre. The optical signal at the output extremity of the fiber was recorded using a mercury-cadmium-tellurium detector with a spectral range of 12,000-600 cm -1 . Liquid nitrogen was used to chill the mercury-cadmium-tellurium detector throughout each experiment. The TAS fiber's edges were manually cut for each experiment.

The edge had to be flat, perpendicular and clean to have the beam going through the fibre with minimal loss. Then, using a bare fibre terminator and multimode connectors, the fibres were attached to the light source and the detector. The signal between the input and output of the fibre was controlled, with an amplitude around 10,000 for 150 µm TAS fibre. The IR-FEWS spectra were acquired with 16 scans and a resolution of 4 cm -1 for a total acquisition time of 4 seconds. The spectra were recorded every 30 seconds during all the electrochemical tests.

Part III. Origin of two-plateau behavior in the IR-FEWS measurements

The double plateau observed on the discharge curve of the IR-FEWS cycling data can find its explanation in a study conducted by our group where the effect of FEC as an additive is investigated on Na metal. [2] Indeed, it has been concluded that the sodium metal's reactivity towards FEC caused the apparition of a second plateau upon discharge. The second plateau was hypothesized to be caused by sodium plating on the protective SEI layer on sodium metal formed by FEC degradation. To put it more clearly, upon cathode sodiation (discharge) the "fresh" sodium metal plated on the SEI is stripped and once depleted, the bulk sodium metal is stripped instead which causes a sudden drop in cell voltage observed in our IR-FEWS cycling curve. To confirm that it is indeed a Na-metal artefact, we conducted a 3-electrode experiment to separate both working electrode (NVP) and counter electrode (Sodium metal) potentials. Half charged NVP was used as a reference with a stable voltage of 3.375 V vs Na/Na + corresponding to the V 3+ /V 4+ redox couple of NVP. Figure S4a shows that indeed upon NVP sodiation (discharge) the sodium metal potential jumps around 100 mV causing a distinct double plateau behavior on the cell voltage (EWE-ECE) (Figure S4b). The PITT procedure is based on a staircase potential signal applied to the system, as represented in Figure S4a. In the case of a planar electrode, Montella demonstrated that the current response of the system is dictated by the well-known Cottrell relationship for short times and an exponential decay of diffusion current with respect to time for long times, if a restricted diffusion condition and diffusion control (very fast insertion reaction kinetics) are assumed. [3] For short times, i.e. 𝑡 ≪ where 𝐹 is the Faraday's constant (96485 C.mol -1 ), 𝐴 is the specific surface are (0.196 cm²), 𝛥𝑐 is the concentration variation in the film (mol.cm -3 ), 𝐿 is the film thickness (cm) and 𝐷𝑎𝑝𝑝′ is the apparent diffusion coefficient in the solid phase (cm².s -1 ). Then, the apparent diffusion coefficient can be extracted from slope of the natural logarithm of the Equation S2 (Figure S7b), as it is independent of unknown variables.

Solvent polarity increase

ii. Randles-Sevcik Equation

Randles-Sevcik equation was employed to estimate the apparent diffusion coefficient of Na + into NVP and NVPF frameworks. Cyclic voltammetry measurements at multiple rates were performed (Figure S8), then, assuming a Nernstian redox reaction (i.e. a fast reversible electron transfer reaction), the absolute current peak value has a linear relationship (in oxidation and reduction) to the square root of the scan rate of the cyclic voltammetry, following this equation [4] :

𝐼 𝑝 = 0.4463 𝐹𝑆 √ 𝐹 𝑅𝑇 𝑆 𝑒 𝑆 √𝐷 𝑎𝑝𝑝 𝑥 𝐶 𝑚𝑎𝑡 𝜈 1/2 Equation S3
where 𝐹 is the Faraday's constant (96485 C.mol -1 ), 𝑅 is the gas constant (J.K -1 .mol -1 ), 𝑇 is the temperature (K), 𝑆𝑒 and 𝑆 are the specific surface and geometric area, respectively (cm²), 𝑥 is the maximal mole fraction of inserted cation that can be inserted in the host material, 𝐶𝑚𝑎𝑡 is the concentration of Na + cation in the host material (mol.cm -3 ) (0.0077 mol cm -3 for NVPF and 0.0069 mol cm -3 for NVP, respectively) [5] , 𝜈 is the scan rate (V.s -1 ), 𝐼𝑝 is the value of the current peak (A) and 𝐷𝑎𝑝𝑝 is the apparent diffusion coefficient of the inserted cation in the solid phase (cm².s -1 ).

Part VI. MD simulations

Molecular dynamics simulations were performed in the LAMMPS program. [6] Initial configurations and input files were generated using Packmol, [7] VESTA [8] and fftool. [9] The simulation box consisted of two electrodes, NVP and gold ones, separated by the NaPF6/PC solution. The NVP electrode was represented by rhombohedral R-3c space group [10] with the [001] plane exposed towards the solvent, following the recommendations from the literature. [11] The gold electrode was represented by cubic Fm-3m space group [12] with the [100] plane exposed to the solvent. The liquid phase contained 75 ion pairs of NaPF6 salt and 830 molecules of PC, matching the experimental composition. The total system size was approximately 15 000 atoms.

Figure 1 .

 1 Figure 1. CV and frequency profiles of Na3V2(PO4)3 (NVP) (a-c) and Na3V2(PO4)2F3 (NVPF) (d-f) in the three selected electrolyte: 1M NaPF6 in PC, EC:DMC and DMC with 1%vol addition of fluoroethylene carbonate (FEC). The oxidation (charge) and the reduction (discharge) processes are portrayed with the correspondent frequency change. Inset figures in (a-c) are the zoom of the current versus time plots during discharge, indicating a mass loss right before the sodiation process. The 20 th cycle is shown here for all experiments. All CVs were performed at 1 mV s -1 . Loading of the composite electrodes on the quartz resonators range from 50-80 µg cm -2 . Details of electrode preparation are given in Supporting Information file.

Figures

  Figures 2b-c and 3b-c. The origin behind the double-plateau feature on the discharge is discussed in the Supplementary Information file (FigureS4).

Figure 2 .

 2 Figure 2. Operando IR-FEWS measurement of fiber coated and embedded in NVP with NaPF6 DMC + 1vol % FEC electrolyte. (a) Absorbance IR-FEWS spectra acquired at t=0s. (b) Absorbance spectra in the νOCO absorbance region in function of the 2 nd charge (blue) 2 nd discharge (green), 3 rd charge (yellow) and 3 rd discharge (red). (c) Relative variations of the absorbance spectra in the νOCO absorbance region in function of the 2 nd charge (blue) 2 nd discharge (green), 3 rd charge (yellow) and 3 rd discharge (red). (d) Voltage evolution of the NVP/Na cells during 2 nd and 3 rd charge. Absorbance evolution of the 𝜈OCO bands of DMCfree (blue), of DMCcoor. (orange) and of the ratio of the coordinated and free bands (yellow).

Figure 3 .

 3 Figure 3. Operando IR-FEWS measurement of fiber coated and embedded in NVP with NaPF6 PC + 1 vol% FEC electrolyte. (a) Absorbance IR-FEWS spectra acquired at t= 0s. (b) Absorbance spectra in the νOCO absorbance region in function of the 2nd charge (blue) 2nd discharge (yellow), 3rd charge (yellow) and 3rd discharge (red). (c) Relative variations of the absorbance spectra in the νOCO absorbance region in function of the 2nd charge (blue) 2nd discharge (red), 3rd charge (yellow) and 3rd discharge (red). (d) Voltage evolution of the NVP/Na cells during 2nd and 3rd charge. Absorbance evolution of the 𝜈OCO bands of PCfree (blue), of PCcoor. (orange) and of the ratio of the coordinated and free bands (yellow).

Figure 4 .

 4 Figure 4. Apparent solid diffusion coefficients of Na + cations in NVPF (a) and NVP (b) estimated using PITT and Randles-Sevcik equation. Rate performances of NVP and NVPF with a 100% (c) and 50% (d) SOC cutoff. All experiments were conducted using a 1M NaPF6 in EC:DMC + 1%vol FEC electrolyte.

Figure 5 .

 5 Figure 5. EQCM results where CV and frequency variation profiles performed at different sweeping rates of potential on NVPF (a-b) and NVP (c-d) loaded QCM resonators in 1M NaPF6 in EC:DMC + 1%vol FEC. The arrows shows the scan direction.

Figure 6 .

 6 Figure 6. MD simulations of an equilibrium NVP//Na/NaPF6 PC system with no voltage applied. (a) Density profiles of the Na + cations, P atom of the PF6 -anions and C atom of the carbonate group of PC as a function of the distance from NVP electrode. The intensities of Na and P are multiplied by the molar ratio of PC/NaPF6 ≃ 11. Top Na layer of the NVP electrode is considered as z = 0 Å. (b) Simulation snapshot of the NVP//Na/NaPF6 PC system. (c) Coordination number (CN) of P atom of the PF6 -anions and O atom of the carbonyl group of PC around the Na + cation as the function of the distance of Na + to the NVP surface. The cutoff values for the first coordination shell of Na + are 4.55 Å for P and 3.20 Å for O atoms.

Scheme 1 .

 1 Scheme 1. Scheme representing the model of opposing concentration gradient accentuated by the fast intercalation kinetics of the host NVP framework during reduction. The penetration depth of the acoustic wave emanating from the resonator is superimposed.

Figure S2 .

 S2 Figure S2. CV-frequency profiles of NVP loaded resonators cycled in 1M LiPF6 in DMC, EC:DMC and PC (left to right, respectively (a-c) and GCPL profiles of NVP loaded resonators cycled in 1M NaPF6 in DMC, EC:DMC and PC (left to right, respectively) (d-f), experiments were conducted at 1mV.s -1 (a-c) and at 1C rate (d-f), 5 th cycle is shown for all experiments.

Figure S3 .

 S3 Figure S3. Evolution of the % frequency gain (relative to the full amplitude variation) caused by desolvation mass loss following the polarity of our selected solvents. i) in 1M LiPF6 in DMC, EC:DMC (1:1 mol) and PC and ii) in 1M NaPF6 in DMC, EC:DMC (1:1 mol) and PC.

Figure S4 .

 S4 Figure S4. Three electrodes cell cycling data, comprising of NVP:Csp powder (80:20 %wt) as a working electrode, sodium metal as a counter and half charged NVP as a reference electrodes. (a) Potential profile during discharge and charge of NVP and sodium metal and (b) cell voltage profile during the same discharge and charge cycle showcasing the FEC-sodium metal related anomaly.

Part IV.

  Figure S5. Charge/Discharge profile of NVP (a,b) and NVPF (c,d) cycled under two different SOC cutoffs: 100% (a,c) and 50% (b,d).

Figure S6 .

 S6 Figure S6. Charge/discharge profile of LiFePO4 (LFP) at limited SOC cutoff (50%). Charge was performed at slow rate C/5 followed by fast discharge (5C).

Figure S7 .

 S7 Figure S7. SEM images of NVPF (a) and NVP (b) powders. The particle size is indicated on each image.

  Figure S8. (a) Applied potential staircases and current response of the NVP system in a 1M NaPF6 in EC:DMC 1%vol electrolyte. Focus on (b) the natural logarithm of the current response for the 3.47 V vs. Na/Na + staircase.

Figure S9 .

 S9 Figure S9. CV profiles at different scan rates (0.2 to 1.5 mV s -1 ) of NVP (a) and NVPF (c) cycled in 1M NaPF6 in EC:DMC 1%vol electrolyte. Apparent Na + diffusion coefficient of NVP (b) and NVPF (d) calculated from Equation S3 during oxidation (top) and reduction (bottom).

Figure S10 .

 S10 Figure S10. Radial distribution functions (a) and coordination numbers (b) of P of PF6 -and O atoms of the carbonyl group of PC around the Na + ion.

Acknowledgements

EB, OS and JMT acknowledge the "Sorbonne Université -Programme doctoral Génie des Procédés" for the PhD scholarship. OS and JMT also acknowledge the funding from DIM-RESPORE. KG and MS acknowledge the funding from European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement 771294). Computational tasks were performed using resources from GENCI-IDRIS (grant 2022-A0120910463) and from the MeSU supercomputer of the Sorbonne University.

Conflict of Interest

The authors declare no conflict of interest. in weight %) and grinded using a mortar and pestle, n-methyl-2-pyrrolidone (NMP) was then added drop by drop until the desired consistency was achieved. The slurry was then casted onto an aluminum sheet using a doctor blade. The casted film was dried at 120°C in air overnight to allow solvent evaporation, 12.7mm electrodes were then punched and subsequently dried at 80°C in a vacuum oven overnight. The loading was kept at 2-3 mg cm -2 .

Supporting information file for

Spotting interface structuring during

Electrochemical and electrogravimetric characterization: A diluted slurry using NMP as a solvent (72:18:12 active material:Csp:PVDF wt%) was first prepared and sonicated to ensure homogenous particle dispersion. The platinum-patterned 9MHz QCM resonators (Bio-Logic)

were heated at 150°C during the spraying process to evaporate the solvent upon contact and were heated at 200°C during 30 minutes to melt the binder and ensure good film cohesion. The resonators modified with the Na3V2(PO4)3 (NVP) and Na3V2(PO4)2F3 (NVPF) thin films were mounted in an airtight EQCM cell developed in our previous work. [1] The cell is mounted in an Ar-filled glove-box and the measurements were conducted outside the glove-box. The EQCM with motional resistance monitoring (EQCM-R) was performed using a Biologic SP200 The loading of the electrodes was 50-80 µg cm -2 and 10 mg cm -2 for EQCM and IR-FEWS setups, respectively. The thickness estimated from the tap densities of both NVP and NVPF at similar particle size would be as follow: 400-600 nm for the EQCM thin film and ~60 µm for the drop-cast IR-FEWS electrode. The amount of active material utilized estimated from experimental and theoretical capacities would be 56% for NVP and 66% for NVPF. The Na + cation were represented using Aqvist force field, [13] while the parameters of PF6 -were taken from the CL&P model. [14] The intermolecular interaction parameters of the NVP and gold electrodes, given in Table S1, were considered from several sources. [13,[15][16][17] The bonded and non-bonded parameters of PC were taken from the OPLS-AA, [18] except for the partial charges, which were computed using the CHelpG procedure [19] on electron densities obtained at the MP2/cc-pVTZ level of theory at previously optimized geometries with the Gaussian16 program. [20] The PC model was validated by computing density and viscosity of bulk liquid NaPF6/PC system. The system had a total size of 34.92×30.24×160.0 Å, with the thickness of ∼20 Å and ∼12 Å for the NVP and Au electrodes, respectively. The simulation box was periodic in x and y direction, and a slab of 3×lz was added along the z direction. A cutoff of 12 Å was considered for non-bonded interactions.

Part II. EQCM analysis of NVP

The particle-particle particle-mash (PPPM) method was used to evaluate electrostatic energies (the accuracy of 1•10 -5 ), with the slab correction applied. Both electrodes were considered immobile. Bonds terminating with hydrogen atoms were constrained using the SHAKE algorithm. The time step was set to 1 fs. The Nosé-Hoover thermostat was used to keep the temperature of the electrolyte solution at 298 K. The system was equilibrated for 5 ns in the NVT ensemble.

The radial distribution functions around Na + ions and the density profiles along the z direction were obtained from the 50 ns trajectory, which followed the equilibration step. The coordination number plots were obtained from a separate production run, during which a single Na + ion was pulled from the solution towards the NVP surface. For this, the Na + atom was tethered by a harmonic spring of 125 kcal mol -1 Å -1 and moved from 10.0 Å to 2.0 Å (with the respect to the top Na layer of the NVP electrode) with a step of 0.25 Å. At each step, the system was equilibrated for 100 ps and the coordination number was recorded during following 200 ps. The coordination number cutoff was set to 4.55 Å for P of PF6 -and 3.20 Å for O atoms of the the carbonyl group of PC, being defined as the position of the first minimum of the corresponding radial distribution functions (Figure S9). The data analysis was performed using TRAVIS [21] and our self-written tools. The visualization of the simulation box was done in VMD. [22]


domain: math spi

text:
Reduced Order Modeling for Parameterized Electromagnetic Simulation Based on Tensor Decomposition

Xiao-Feng He, Liang Li, Stéphane Lanteri, Kun Li

Abstract-We present a data-driven surrogate modeling for parameterized electromagnetic simulation. This method extracts a set of reduced basis (RB) functions from full-order solutions through a two-step proper orthogonal decomposition (POD) method. A mapping from the time/parameter to the principal components of the projection coefficients, extracted by canonical polyadic decomposition (CPD) , is approximated by a cubic spline interpolation (CSI) approach. The reduced-order model (ROM) is trained in the offline phase, while the RB solution of a new time/parameter value is recovered fast during the online phase. We evaluate the performance of the proposed method with numerical tests for the scattering of a plane wave by a 2-D multi-layer dielectric disk and a 3-D multi-layer dielectric sphere.

Index Terms-electromagnetic simulation, model order reduction, proper orthogonal decomposition, canonical polyadic decomposition, cubic spline interpolation

I. INTRODUCTION

E LECTROMAGNETIC scattering problems are generally described as the parameterized time-domain Maxwell equations [START_REF] Jin | Theory and computation of electromagnetic fields[END_REF],where the parameters may represent geometric shapes, material properties, system configurations and so on. Solving the full-order model (FOM) repeatedly for varying parameter values by high-fidelity simulations is often required in many applications such as real-time system control, weather forecast and structure design, leading to great CPU time and memory demands. It is urgent to study the parametric surrogate modeling to accelerate the solution of the partial differential equations (PDEs) in an efficient and accurate way. Model order reduction (MOR) techniques have been proposed to solve this problem in the past decades. The main idea of MOR is to replace a high-dimensional FOM with a surrogate model, which has a lower dimensionality and acceptable accuracy loss [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF].

Some techniques such as the fast multipole method (FMM) [START_REF] Song | Multilevel fast multipole algorithm for electromagnetic scattering by large complex objects[END_REF], subdomain multilevel approach (SMA) [START_REF] Suter | A subdomain multilevel approach for the efficient mom analysis of large planar antennas[END_REF] and characteristic basis functions method (CBFM) [START_REF] Vvs Prakash | Characteristic basis function method: A new technique for efficient solution of method of moments matrix equations[END_REF] have been proposed to accelerate the solution of scattering problems. For example, CBFM has been developed for an efficient solution of Method of Moments matrix equations. Compared

The first two authors are with the School of Mathematical Sciences, University of Electronic Science and Technology of China (hexf@std.uestc.edu.cn, plum.liliang@gmail.com, plum liliang@uestc.edu.cn). Stéphane Lanteri is with Inria ATLANTIS project-team (Stephane.Lanteri@inria.fr). Kun Li is with the School of Mathematics, Southwestern University of Finance and Economics (rclikun@163.com, likun@swufe.edu.cn). The corresponding author is Liang Li.

Manuscript received March 25, 2023 to conventional RWG basis functions [START_REF] Rao | Electromagnetic scattering by surfaces of arbitrary shape[END_REF] bounded by the λ/20 domain discretization, this approach is based on high order basis functions over enlarged blocks, thus reducing the number of degrees of freedom (DoFs). A reduced basis (RB) method adopted in this work reduces the dimensionality of FOM in a completely different manner. In general, repeated simulations are usually not necessary because intrinsic similarities among the solutions of different parameters can be used to predict the solutions for new parameter values [START_REF] Guo | Data-driven reduced order modeling for time-dependent problems[END_REF].

The RB method with an offline-online procedure, is one of the most widely used techniques for the surrogate modeling of time-dependent parameterized problems. In the offline phase, a reduced space representing the main dynamics of the FOM is constructed from full-order snapshots of different time/parameter values. There are several methods to extract the reduced basis functions, such as the proper orthogonal decomposition (POD) [START_REF] Jan S Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF], [START_REF] Li | A reduced-order discontinuous galerkin method based on pod for electromagnetic simulation[END_REF], [START_REF] Liang | Proper orthogonal decomposition and its applications-part i: Theory[END_REF], greedy algorithms [START_REF] Lappano | A greedy sampling approach for the projection basis construction in parametric model order reduction for structural dynamics models[END_REF], [START_REF] Murata | Nonlinear mode decomposition with convolutional neural networks for fluid dynamics[END_REF] and an autoencoder [START_REF] Fresca | A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized pdes[END_REF], [START_REF] Francisco | Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems[END_REF]. As the most widely-used technique, POD is used in this paper to generate the RB functions with a singular-value decomposition (SVD) [START_REF] Guo | Data-driven reduced order modeling for time-dependent problems[END_REF]. As an approximation to the full-order solution, the RB solution, also named the reduced-order solution is recovered as a linear combination of the RB functions, where only the combination coefficients, also named projection coefficients will be calculated in the online phase [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF]. It is advised that the methods mentioned above can be combined to accelerate the solution of parameterized electromagnetic scattering problems. For example, we could build a data-driven ROM based on snapshots generated by CBFM offline, then recover the RB solution for a new parameter online. RB methods are divided into intrusive and non-intrusive (or data-driven) ones according to the methodology of calculating the projection coefficients. The former determines the projection coefficients by solving a low dimensional linear system, i.e., a projection of the full-order model onto the reduced space. And the Galerkin procedure is often used for the projection [START_REF] Couplet | Calibrated reduced-order podgalerkin system for fluid flow modelling[END_REF], [START_REF] Quarteroni | Certified reduced basis approximation for parametrized partial differential equations and applications[END_REF]. However, the main drawback is that intrusive ROMs require access to the high dimensional model, and provide limited computational gain especially for complex nonlinear problems [START_REF] Ballarin | Supremizer stabilization of pod-galerkin approximation of parametrized steady incompressible navier-stokes equations[END_REF].The non-intrusive ROMs approximate the projection coefficients over the parameter domain by using a database of reduced order information [START_REF] Casenave | A nonintrusive reduced basis method applied to aeroacoustic simulations[END_REF]. The high-fidelity solvers such as FEM, FDM, DGTD are only used to generate snapshots, leading to a full decoupling of the online stage and the high-fidelity scheme. Many excellent non-intrusive ROMs are proposed in the past. For instance, a feed-forward neural network (FNN) has been designed for parameterized unsteady flows [START_REF] Hesthaven | Non-intrusive reduced order modeling of nonlinear problems using neural networks[END_REF], [START_REF] Wang | Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem[END_REF], where a FNN is trained to approximate the projection coefficients. Interpolation and regression methods such as cubic spline interpolation (CSI) [START_REF] Bui-Thanh | Proper orthogonal decomposition extensions for parametric applications in compressible aerodynamics[END_REF], [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF], radial basis functions (RBF) [START_REF] Xiao | A parameterized non-intrusive reduced order model and error analysis for general time-dependent nonlinear partial differential equations and its applications[END_REF] and Gaussian process regression (GPR) [START_REF] Guo | Data-driven reduced order modeling for time-dependent problems[END_REF], [START_REF] Zhao | Non-intrusive reduced order modeling of parametric electromagnetic scattering problems through gaussian process regression[END_REF] have also been adopted to approximate the projection coefficients. By combining an FNN-based regression model with a physics-informed neural network (PINN) [START_REF] Lu | Deepxde: A deep learning library for solving differential equations[END_REF], [START_REF] Raissi | Physicsinformed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF], a hybrid surrogate model is devised in [START_REF] Chen | Physics-informed machine learning for reduced-order modeling of nonlinear problems[END_REF], which requires less training data.

Data-driven reduced-order modeling based on RB method has been applied to electromagnetic scattering simulations [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF], [START_REF] Zhao | Non-intrusive reduced order modeling of parametric electromagnetic scattering problems through gaussian process regression[END_REF] recently. These related works use the POD technique to extract a set of time-/parameter-independent orthogonal basis functions and compute a series of projection coefficient matrices, then perform SVD for each projection matrix to extract the dominating temporal and parametric modes. Interpolation or regression models are built to approximate the temporal and parametric modes respectively. These methods are accurate enough for time-dependent parametric problems, but require to build a large number of interpolation or regression models (determined by the number of basis functions and the number of truncated singular values of each projection coefficient matrix), which increases the test time online. Motivated by the above problem, a ROM composed of two-step POD, canonical polyadic decomposition (CPD), and cubic spline interpolation (CSI) is developed for parameterized electromagnetic simulation. We employ a canonical polyadic decomposition (CPD) [START_REF] Tamara | Tensor decompositions and applications[END_REF] to extract the principal components of the projection coefficients tensor, instead of repeated SVD for all projection coefficient matrices, which significantly reduces the number of interpolation or regression models, thus reducing the test time online with an acceptable loss of accuracy. The proposed method, named POD-CPD-CSI ROM is validated to be efficient and accurate for parameterized electromagnetic scattering problems.

The symbols and notations used in this paper are introduced as follows. Notations t and µ denote time and parameter respectively. S j u is the time trajectory matrix and S u is the snapshot matrix. N h is the geometric DoFs, N t and N p are the number of time and parameter samples respectively. u h (t, µ) and u r h (t, µ) are the full-order solution and reducedorder solution respectively. α u (t, µ) denotes the projection coefficient vector and V is the POD basis matrix. P u is the projection coefficient matrix and G u is the projection coefficient tensor.

II. TWO-STEP POD METHOD

Assume that

P tr h = {µ 1 , µ 2 , • • • , µ
Np } is a parameter sampling that we're interested in. For each µ j ∈ P tr h , we have a time trajectory matrix

S j u = u h (t1, µ j ) | u h (t2, µ j ) | • • • | u h (tN t , µ j ) ∈ R N h ×N t , (1)
with u ∈ {E, H}, which could be obtained either from highfidelity solutions or from experiments, for a time sampling

T tr h = {t 1 , t 2 , • • • , t Nt }.
Here, N h denotes the number of geometric DoFs. A snapshot matrix for all parameters in P tr h is then formed

S u = S 1 u | S 2 u | • • • | S Np u ∈ R N h ×Ns , (2)
with N s = N t • N p . To reduce the dimensionality of FOM, we seek a low-rank decomposition to the snapshot matrix and construct a reduced space

V u,rb = span {v u,1 , v u,2 , • • • , v u,du } of dimension d u min{N h , N s },
where {v u,i } du i=1 are time-and parameter-independent orthogonal RB functions [START_REF] Guo | Data-driven reduced order modeling for time-dependent problems[END_REF], [START_REF] Zhao | Non-intrusive reduced order modeling of parametric electromagnetic scattering problems through gaussian process regression[END_REF].

For the large snapshots matrix S u , i.e. when N s is large, it will be very computationally expensive to directly perform POD on S u . So, we adopt a two-step POD method [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF] shown in Algorithm 1 to reduce the demands on computational cost. Details of the computational complexity and error analysis of this method can be found in [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF]. The RB solution u r h (t, µ), an approximation to the full-order solution u h (t, µ), is recovered as

u r h (t, µ) = V u α u (t, µ) = du i=1 α u,i (t, µ)v u,i , (3)
where V u ∈ R N h ×du is the POD basis matrix generated by the two-step POD and vector

α u (t, µ) = [α u,1 (t, µ), • • • , α u,du (t, µ)]
T is called the coefficient vector [START_REF] Zhao | Non-intrusive reduced order modeling of parametric electromagnetic scattering problems through gaussian process regression[END_REF]. After calculating the POD basis functions, we need to build a mapping between the time/parameters and the projection coefficients.

Algorithm 1: Two-step POD approach

input : S 1 u , S 2 u , • • • , S
Np u and truncation parameters ρ t and ρ µ output:

POD basis matrix V u 1 Calculate matrices T j u =POD (S j u , ρ t ) for j = 1, 2, • • • , N p ; 2 Assemble matrix T u = [T 1 u | T 2 u | • • • | T Np u ]; 3 Compute matrix V u =POD (T u , ρ µ ).
(d) = d i=1 λ i / r i=1 λ i ; 7 Calculate the RB functions v i = 1 √ λi Au i for i = 1, 2, • • • , d, d r; 8 Assemble matrix V = [v 1 , v 2 , • • • , v d ]. 9 end

III. CP-DECOMPOSITION-BASED INTERPOLATION

Each vector in the snapshot matrix (2) can be approximated by a linear combination of the RB basis vectors (3) [START_REF] Zhao | Non-intrusive reduced order modeling of parametric electromagnetic scattering problems through gaussian process regression[END_REF].

In such a way, we create training data and rewrite all the projection coefficients with a matrix form

P u =      α u, (t1, µ 1 ) α u, (t1, µ 2 ) • • • α u, (t1, µ Np ) α u, (t2, µ 1 ) α u, (t2, µ 2 ) • • • α u, (t2, µ Np ) . . . . . . . . . . . . α u, (tN t , µ 1 ) α u, (tN t , µ 2 ) • • • α u, (tN t , µ Np )      = v T u, S 1 u ; v T u, S 2 u ; • • • ; v T u, S Np u T ∈ R N t ×Np , (4)
with = 1, 2, • • • , d u , u ∈ {E, H}. Instead of dealing with matrix P u directly [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF], [START_REF] Zhao | Non-intrusive reduced order modeling of parametric electromagnetic scattering problems through gaussian process regression[END_REF], we concatenate the coefficient matrices to form a coefficient tensor G u of dimension (N t , N p , d u ). The CP decomposition [START_REF] Tamara | Tensor decompositions and applications[END_REF] illustrated in Fig. 1 factorizes G u into a sum of component rank-one tensors

G u ≈ Ĝu = Ru r=1 ϕ u,r • ψ u,r • ζ u,r , (5)
where R u is a positive integer and determines the error of CP decomposition. The symbol • represents the vector outer product and ϕ u,r ∈ R Nt , ψ u,r ∈ R Np , and ζ u,r ∈ R du for r = 1, . . . , R u are vectors. Our goal is to compute a CP We adopt an alternating least squares (ALS) method [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an n-way generalization of "eckartyoung" decomposition[END_REF], [START_REF] Richard A Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF] restated in Algorithm 2 to find the CP decomposition in this work, where G u(n) is the mode-n unfolding of tensor G, and notation and denote the Khatri-Rao product and Hadamard product respectively. A † denotes the Moore-Penrose pseudoinverse of A. Readers are referred to [START_REF] Tamara | Tensor decompositions and applications[END_REF] for more details.

With the discrete data set, a cubic spline interpolation-based model u is trained to approximate the continuous modes as

t -→ ϕ u,r (t) with {(t i , (ϕ u,r ) i ), i = 1, 2, • • • , N t }, µ -→ ψ u,r (µ) with {(µ j , (ψ u,r ) j ), j = 1, 2, • • • , N p }, (7)
for r = 1, 2, • • • , R u . Elementwise, we have

(G u ) ijk ≈ ( Ĝu ) ijk ≈ Ru r=1 ϕ u,r (t i ) ψ u,r (µ j )(ζ u,r ) k , (8)
with

1 ≤ i ≤ N t , 1 ≤ j ≤ N p , 1 ≤ k ≤ d u .
For a new time/parameter value (t * , µ * ), the RB solution is recovered as

u r h (t * , µ * ) = V u α u (t * , µ * ) = V u Ru r=1 ϕ u,r (t * ) ψ u,r (µ * )ζ u,r . (9)
Algorithm 2: ALS method for CP decomposition of the projection coefficient tensor G u input : projection coefficient tensor G u ∈ R Nt×Np×du and a truncation parameter R u . output:

factor matrices Φ u = [ϕ u,1 , • • • , ϕ u,Ru ], Ψ u = [ψ u,1 , • • • , ψ u,Ru ] and Σ u = [ζ u,1 , • • • , ζ u,Ru ]. 1 Initialize Φ u ∈ R Nt×Ru , Ψ u ∈ R Np×Ru and Σ u ∈ R du×Ru ; 2 repeat 3 Φ u = G u(1) (Σ u Ψ u )(Σ T u Σ u Ψ T u Ψ u ) † 4 Ψ u = G u(2) (Φ u Σ u )(Φ T u Φ u Σ T u Σ u ) † 5 Σ u = G u(3) (Ψ u Φ u )(Ψ T u Ψ u Φ T u Φ u ) † 6
until fit cases to improve or maximum iterations exhausted;

The procedure of POD-CP ROM is decoupled into offline/online phases, which is shown in Algorithm3, Fig. 2 and Fig. 3. Compared to previous works [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF], [START_REF] Zhao | Non-intrusive reduced order modeling of parametric electromagnetic scattering problems through gaussian process regression[END_REF],fewer interpolation or regression models are built in the offline stage, thus reducing the test time to evaluate the RB solution online. Specifically, for each EM field component, 2d u du =1 Q u, interpolation-based models are built in [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF], where Q u, denotes the number of singular values truncated of coefficient matrix P u , while only 2R u interpolation-based models are built in POD-CP ROM.

The recovery error the POD-CPD-CSI ROM at (t * , µ * ) will be written as

ε = u h (t * , µ * ) -u r h (t * , µ * ) 2 = u h (t * , µ * ) -V u α u (t * , µ * ) 2 = ||u h (t * , µ * ) -V u V T u u h (t * , µ * ) + V u V T u u h (t * , µ * ) -V u α u (t * , µ * )|| 2 ε 2 POD + V u V T u u h (t * , µ * ) -α u (t * , µ * ) 2 = ε 2 POD + N i=1 (α u,i (t * , µ * ) -α u,i (t * , µ * )) 2 = ε 2 POD + ε 2 coe , (10)
where ε 2 POD is the truncation error of the two-step POD. And ε 2 coe denotes the error of predicting the projection coefficients, which includes the error of CP decomposition and the CSI error of ϕ u,r (t * ) and ψ u,r (µ * ).

IV. NUMERICAL RESULTS

In this section, we present four numerical tests to evaluate the performance of the proposed POD-CPD-CSI ROM. Moreover we compare the proposed ROM with the POD-ANN ROM [START_REF] Wang | Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem[END_REF] and POD-CSI ROM [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF] both in time-efficiency and errors. We consider non-magnetic materials (ν r = 1) in a vacuum for all numerical tests. The norm of a tensor G is defined as We determine the truncation parameter R u of the CP decomposition by the relative error indicator e u,CP = G u -Ĝu / G u . In the online phase, for a new time value t ∈ T te h and physical parameter value µ ∈ P te h , the relative projection error

G = Nt i Np j du k [(G) ijk ] 2 . ( 11
)
e u,Pro (t, µ) = u h (t, µ) -V u V T u u h (t, µ) 2 u h (t, µ) 2 , ( 12
)
provides a lower bound for the relative error of the POD-CPD-CSI ROM

e u,POD-CPD-CSI (t, µ) = u h (t, µ) -V u α u (t, µ) 2 u h (t, µ) 2 . (13
N te . ( 15
)
We use the function cp_als in the MATLAB library Tensor Toolbox [START_REF] Bader | Tensor Toolbox for MATLAB[END_REF] to decompose the coefficient tensor G u . The CSI-based approximation of the projection coefficients is built via the MATLAB functions spline and griddedInterpolant for single-layer cases and multilayer cases respectively. We refer to [START_REF] Sakib Hasan | Modeling of soi four-gate transistor (g4fet) using multidimensional spline interpolation method[END_REF], [START_REF] Wang | Multivariate spline functions and their applications[END_REF] for more details of the multivariate CSI method. The architecture of ANN used in the POD-ANN ROM is listed in Table VI. All simulations are run on a desktop with an Intel Core i9 10-core 2.8 GHz × 20 CPU and 64 GB RAM. A. Scattering of a plane wave by a dielectric disk First, we consider the scattering problem described by the 2D time-domain TM waves in the normalized form

     ν r ∂H ∂t + curl E = 0, ε r ∂E ∂t -curl H = 0, (16)
where the magnetic field H = (H x , H y ) T and the electric field E = E z . An incident plane wave takes the form H inc x (x, y, t) = 0, H inc y (x, y, t) = -cos(ωt -kx) and E inc z (x, y, t) = cos(ωt -kx) with a frequency f = 300 MHz [START_REF] Zhao | Non-intrusive reduced order modeling of parametric electromagnetic scattering problems through gaussian process regression[END_REF]. The computational domain is a square Ω = [-2.6λ 0 , 2.6λ 0 ] 2 and the radius of the inner disk is 0.6 m (background wavelength λ 0 = 1m). The range of the relative permittivity ε r is presented in TableI, i.e., the parameter domain we are interested in is

P = {µ = ε r ∈ [1, 5]}.
To determine the number of sampling parameter points in P, we create three training parameter sets uniformly with different sampling intervals ∆p = 1, 0.5 and 0.05, resulting in N p = 5, 9 and 81 points respectively. A DGTD solver [START_REF] Li | Podbased model order reduction with an adaptive snapshot selection for a discontinuous galerkin approximation of the time-domain maxwell's equations[END_REF] is used to generate a collection of high-fidelity solutions of the three training parameter sets with N t = 263 samples in the last period T tr h = {49.0024, 49.006, • • • , 49.966}. We evaluate the proposed method on a test parameter set P te h = {1.215, 2.215, 3.215, 4.215} and a test time set T te h = T tr h . Note that the truncation parameter R u is selected carefully such that the projection coefficient tensor G u lose little information in the CP decomposition. Fig. 4 shows the average relative errors ēu,Pro and ēu,POD-CPD-CSI on the training set P tr h × T tr h with a fixed truncation parameter ρ t = 0.1. The performance of the POD-CP ROM mainly depends on the size of training set as well as the truncation parameter ρ µ used. So the training parameter set P tr h = {1, 1.05, • • • , 5} is adopted in the following discussion.

The convergence curves of the average relative errors ēH,Pro and ēE,Pro are plotted in Fig. 5, from which we can see that the tolerance ρ t has little impact on the performance of the two-step POD method. By decreasing the tolerance ρ µ from 1×10 -1 to 1×10 -5 , there is a significant improvement of the two-step POD performance. Therefore, we set ρ t = 1 × 10 -1 and ρ µ = 1 × 10 -5 . The number of POD basis functions d u generated by the two-step POD method are listed in Table ??. The relative errors e u,CP for different truncation parameters R u adopted in the CP decomposition are plotted in Fig. 6. We take R Hx = R Hy = R Ez = 50 to minimize the information lost in CP decomposition while constructing as few interpolation models as possible. VII. It can be seen that the accuracy of POD-CSI and POD-CP-CSI are similar, and both are better than POD-ANN.

In addition, Table VIII shows the average computing time of a time trajectory matrix (1) between DGTD and three ROMs. The online test time of the ROM achieves 5520 times faster than that of the DGTD solver, which demonstrates the high efficiency of the proposed method. The non-intrusive ROMs features an offline-online framework, i.e., offline training is completely decoupled from online testing. Given a new parameter, the RB solution is rapidly recovered without an access to the high-dimensional FOM, which greatly improves the efficiency of the ROM. In addition, POD-CSI and POD-CPD-CSI take less training time offline compared with POD-ANN [START_REF] Wang | Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem[END_REF]. Finally, as shown in Table IV, POD-CPD-CSI generates less CSI models and takes less test time compared to POD-CSI.

B. Scattering of a plane wave by a multi-layer disk

Next, we consider the case of a multi-layer dielectric disk as shown in Fig. 11. A square Ω = [-3.2λ 0 , 3.2λ 0 ] 2 is the computational domain, and the radius and the range of material parameters of each layer are presented in Table I. The excitation is the same as the case in IV-A.

To determine the number of sampling parameter points in P = P 1 × P 2 × P 3 × P 4 with µ = (ε r,1 , ε r,2 , ε r,3 , ε r,4 ) and ε r,i ∈ P i , we create two training parameter sets uniformly with N ∆p = 2 and 3 in each layer, resulting in N p = 16 and 81 points in total. For each parameter point, N t = 253 snapshots are calculated in time at the last period ( t = 0.0038 m and T tr h = {49.0002, 49.0041, • • • , 49.9669}). We evaluate the proposed method on a test set P

C. Scattering of a plane wave by a dielectric sphere

We consider a plane wave that is polarized in the x direction, propagates in the +z direction, and is scattered by a dielectric sphere [START_REF] Zhu | Sphere scattering[END_REF]. The plane wave is expressed as

E inc = e -jkz e x , H inc = E inc x η e -jkz e y , (17)
where the magnitude of the electric part is 1 V/m with a frequency f = 300 MHz. The computational domain is a cube Ω = [-0.2λ 0 , 2.2λ 0 ] 3 with λ 0 = 1 m, and the range of material parameters is listed in The average computing time and average relative errors are presented in TableVIII and TableVII respectively. D. Scattering of a plane wave by a multi-layer sphere Finally, we consider a plane wave scattered by a 3-layer dielectric sphere [START_REF] Zhu | Sphere scattering[END_REF] as shown in Fig. 11. The computational domain is Ω = [-0.2λ 0 , 2.2λ 0 ] 3 , and the radius and range of material parameters of each layer are listed in Table I. The excitation is the same as the case in IV-C. In the offline phase, full-order solutions are calculated by the Mie series method on N h = 91125 nodes at N p = 125 uniformly distributed parameter points. A test parameter set P The computing time and average relative errors are presented in TableVIII and TableVII respectively.

V. CONCLUSION

A data-driven RB method based on POD, CP decomposition and CSI is proposed to accelerate the solution of parameterized time-domain Maxwell equations. RB solutions can be obtained at a very low cost from the online phase. We verify the performance of the proposed POD-CP ROM with scattering by a medium disk and a dielectric sphere. for scattering problems with complex geometric objects, more POD basis functions generated by the two-step POD might be required to capture the main dynamics of the full-order system.

We remark that although data-driven ROMs are effective and easy to implement when sufficient training data (experimental or computational data) is available, it may be still expensive to generate relatively large data in many application fields. In addition, data-driven ROMs fit observed data well, but generally have poor generalization performance for unseen data. A possible approach to address these issues is to take physical information into account. In fact, the governing partial differential equations (PDEs) behind the training data are clearly known, but often ignored by data-driven ROMs. Future research direction is to combine the data and physical information to improve the generalization capacity of ROMs when training data is limited.

4 5

 5 Function V=POD (A, ρ) Compute the eigenvalues {λ i } r i=1 and the eigenvectors {u i } r i=1 of matrix A T A with r = rank(A); 6 Determine the number of basis functions d based on d = arg min{π(d) : π(d) 1 -ρ} and π

Fig. 1 .

 1 Fig. 1. CP decomposition of a three-order tensor Gu.decomposition with R u components that best approximates G u , i.e., to find

Fig. 2 .Fig. 3 .

 23 Fig. 2. Offline satge of the proposed POD-CPD-CSI ROM.

Algorithm 3 : 2 3

 32 POD-CPD-CSI ROM 1 Function [V u , u ]=Offline(P, T ) Prepare a training parameter sampling P tr h = {µ 1 , µ 2 , • • • , µ Np } ⊂ P;Compute the high-fidelity solutions u h (t i , µ j ) in the time domain T h ⊂ T for each physical parameter µ j , j = 1, 2

  ) Furthermore, the average relative errors on a test set T te h ×P te h of size N te are calculated to measure the accuracy of the POD-CPD-CSI ROM, which are defined as ēu,POD-CPD-CSI = (t,µ)∈T te h ×P te h e u,POD-CPD-CSI (t, µ)

Fig. 4 .

 4 Fig. 4. The average relative errors of ēH,Pro and ēH,POD-CP (top), ēE,Pro and ēE,POD-CP (bottom) on the training set with different sampling intervals ∆p and truncation parameters ρµ.

Fig. 5 .Fig. 6 .

 56 Fig. 5. The average relative errors of ēH,Pro (top) and ēE,Pro (bottom) with different truncation parameters ρt and ρµ. For the parameter instance µ 3 = 3.215, comparison of the time evolution of the field H y and E z at a fixed point is displayed in Fig. 7. Moreover, the total field of |H x | calculated by DGTD and the absolute errors between POD-ANN, POD-CSI, POD-CPD-CSI and DGTD respectivily are shown in Fig. 8. The time evolution of the relative error e u,Pro and e u,POD-ANN , e u,POD-CSI and e u,POD-CPD-CSI is plotted in

Fig. 9 .

 9 Fig. 9. The average projection error ēu,Pro and three ROMs are listed in TableVII. It can be seen that the accuracy of POD-CSI and POD-CP-CSI are similar, and both are better than POD-ANN.In addition, TableVIIIshows the average computing time of a time trajectory matrix (1) between DGTD and three ROMs. The online test time of the ROM achieves 5520 times faster than that of the DGTD solver, which demonstrates the high efficiency of the proposed method. The non-intrusive ROMs features an offline-online framework, i.e., offline training is completely decoupled from online testing. Given a new parameter, the RB solution is rapidly recovered without an access to the high-dimensional FOM, which greatly improves the efficiency of the ROM. In addition, POD-CSI and POD-CPD-CSI take less training time offline compared with POD-ANN[START_REF] Wang | Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem[END_REF]. Finally, as shown in TableIV, POD-CPD-CSI generates

Fig. 7 .

 7 Fig. 7. Comparison of the time evolution of the Hy (top) and Ez (bottom) at a fixed point.

Fig. 8 .

 8 Fig. 8. The DGTD solution of |Hx| (a) at t = 49.966 m and the absolute errors between POD-ANN (b), POD-CSI (c), POD-CPD-CSI (d) and DGTD respectively.

Fig. 9 .

 9 Fig. 9. Comparison of the time evolution of the relative errors e u,Pro , e u,POD-ANN , e u,POD-CSI and e u,POD-CPD-CSI for the field H and E with the testing-parameter instances ε 3 = 3.215.

  te h ×T te h where T te h = T tr h and P te h = {µ 1 , µ 2 , µ 3 , µ 4 } with µ 1 = {(5.1, 3.3, 2.2, 1.3)}, µ 2 = {(5.3, 3.5, 2.2, 1.5)}, µ 3 = {(5.5, 3.4, 2.5, 1.7)} and µ 4 = {(5.6, 3.7, 2.4, 1.6)}.

Fig. 10 .

 10 Fig. 10. The average relative errors of ēH,Pro and ēH,POD-CPD-CSI (top), ēE,Pro and ēE,POD-CPd-csi (bottom) on the training set with different number of sampling points N ∆p and truncation parameters ρµ. The average relative errors ēu,Pro and ēu,POD-CPD-CSI on the training set P tr h × T tr h with a fixed truncation parameter ρ t = 0.1 are shown in Fig. 10, therefore we select N ∆p = 3 parameter points uniformly in each layer. The parameter setting in the POD-CPD-CSI ROM is listed in TableIII.The average computing time of a time trajectory matrix and average relative errors are presented in TableVIII and TableVII respectively.

Fig. 11 .

 11 Fig. 11. Configuration and scattering of a plane wave by a multi-layer disk (top) and multi-layer sphere (bottom).

  te h = {(2.1, 8.4, 4.2), (2.2, 8.3, 4.4), (2.3, 8.2, 4.1), (2.4, 8.1, 4.3)} and test time set T te h = T tr h = {0, 0.0101, • • • , 1} are used to evaluate the proposed method. The number of basis functions and truncation tolerances of two-step POD and CP decomposition are presented in TableIII.

TABLE VII AVERAGE

 VII RELATIVE ERRORS OF THREE ROMS.

	Case	Projection	POD-ANN	POD-CSI	POD-CPD-CSI
		ēH	ēE	ēH	ēE	ēH	ēE	ēH	ēE
	disk	1.11e-2 1.36e-2 3.24e-2	4.27e-2	1.16e-2 1.36e-2 1.19e-2 1.36e-2
	multi-layer disk	7.37e-3 8.23e-3 2.60e-2	3.03e-2	8.01e-3	8.26e-3	8.03e-3 8.86e-3
	sphere	6.64e-4 1.44e-3 3.97e-2 3.21e-2	6.73e-4	1.49e-3	6.89e-4 1.50e-3
	multi-layer sphere 7.34e-4 4.02e-3 3.10e-2	4.03e-2	7.50e-4	4.12e-3	7.58e-4 4.15e-3

TABLE VIII COMPARISON

 VIII OF THE CPU TIME OF FOM AND THREE ROMS. THE UNIT OF TIME IS SECOND.

	Case		FOM	POD-ANN	POD-CSI	POD-CPD-CSI
		DoFs	CPU time training	test	training	test	training	test
	disk	30264	258	32450	0.02	9	0.25	19	0.04
	multi-layer disk	37236	305	31784	0.03	7	0.19	16	0.05
	sphere	91125	343	43287	0.03	23	0.34	35	0.07
	multi-layer sphere 91125	3710	39749	0.03	27	0.39	38	0.06

  = 100 points in a period. We evaluate the proposed method on a test parameter set P te h = {2.215, 3.215, 4.215, 5.215} and test time set T te h = T tr h = {0, 0.0101, • • • , 1}. The number of basis functions and truncation tolerances of the two-step POD and CP decomposition are presented in TableIII.

Table I. During the offline phase, high-fidelity solutions are calculated by the Mie series method [9] on N h = 91125 nodes at N p = 81 uniformly distributed parameter points (P tr h = {2, 2.05, • • • , 5.95, 6}) with N t

ACKNOWLEDGEMENT

The last author is supported by NSFC, 12101511 and Fundamental Research Funds for the Central Universities (Grant No. JBK22YJ37).

Stéphane Lanteri is a senior research scientist at Inria center of Université Côte d'Azur. Sophia He holds a PhD in Engineering Sciences from University of Nice-Sophia Antipolis (Year of defense: 1991) and a HDR (habilitation to advise doctoral theses) from the University of Nice/Sophia Antipolis (Year of defense: 2003 -Title: High-performance numerical methods on unstructured meshes with applications to compressible fluid dynamics). From May 1992 to October 1993, he was a postdoctoral fellow at the Center for Aerospace Structure, University of Colorado at Boulder, under the supervision of Charbel Farhat. His current research interests are concerned with numerical modeling of physical problems related to computational electromagnetics and computational nanophotonics: design of high order discontinuous Galerkin type approximation methods on unstructured meshes and of parallel numerical algorithms for solving differential systems modeling electromagnetic wave interaction with complex media with a focus of nanoscale light-matter interaction. Since February 2020, he is the scientific head of the Atlantis project-team at the Inria Research Center of Université Côte d'Azur. He is also the scientific coordinator of the development of the DIOGENeS software suite, which is dedicated to computational nanophotonics. Stéphane Lanteri has authored or co-authored more that 80 publications in international journals in applied mathematics, scientific computing, computational physics and more recently, optics and photonics. Web: https://slanteri.github.io/ Kun Li received the Ph.D. degree in computational mathematics from the University of Electronic Science and Technology of China, Chengdu, China, in 2020. He is currently an assistant professor with the School of Mathematics, Southwestern University of Finance and Economics, Chengdu, China. His current research interests include the discontinuous Galerkin time-domain methods and model order reduction techniques for solving the electromagnetic, and local/nonlocal optical response problems.