Intact NMR spectroscopy: slow high-resolution magic angle spinning chemical shift imaging

Covadonga Lucas-Torres, Alan Wong* High-Resolution Magic-Angle Spinning Chemical Shift Imaging (HR-MAS CSI) has recently been explored for nuclear magnetic resonance (NMR) metabolomics and showed a considerable promise in organism research. This is owing to its capability of offering a supplemental dimension -spatial metabolic distribution -in the profiling. However, HR-MAS CSI suffers from the large centrifugal stress exerted onto the sample, which inevitably hinders the metabolic assessments. Herein, a slow sample spinning strategy was implemented and evaluated. The results demonstrate a potential -with high quality data and feasibility -as a highly informative profiling approach to intact specimens.

Metabolomics is one of the significant OMICS fields which study the chemical profiles of the low-molecular weight metabolites (<1 kDa) those responses to the metabolism of interest. It plays a vital role of untangling many biocomplexities in life sciences. 1 This is, in fact, the consequence of the vast advancement of bioanalytical technologies; they play a significant role on the development of metabolomics. Among them, NMR spectroscopy has found tremendous success due to its simplicity and versatility of acquiring unbiased and rich metabolic information on diverse biospecimens. In particular, High-Resolution Magic-Angle Spinning (HR-MAS) is a valuable NMR tool for profiling the metabolic constituents inside the heterogeneous specimens such as whole cells and organisms, and tissues. It benefits from the capacity of offering high spectral quality data by rapidly rotating the sample at a 54.74° (magic-angle) to the NMR magnetic field. This suppresses the unwanted magnetic susceptibility line broadening ascribed within the sample. For this reason, HR-MAS NMR is a wellestablished technique and vastly applied. This includes intact NMR studies of a single mm-scaled organism, 2 and even a smallscaled animal. 3 However, despite its versatility, HR-MAS only offers a global metabolic view of the entire organism, masking the profiles of the specific anatomy.

Recently, an advanced development of HR-MAS methodology -with chemical shift imaging (CSI) -has been introduced 2b and explored 4 towards metabolomics. HR-MAS CSI offers an additional spatial dimension for the possibility of a homogeneous profiling (i.e. localized spectroscopy along the sample rotational axis) from a heterogeneous specimen such as an intact organism. The experiment is taking advantage of the sample spinning with HR-MAS for providing the essential high quality spectral data, and a gradient field for localizing the metabolic profiles along the rotational axis z' (magic-angle axis). Indeed, the capability of simultaneous profiling the different anatomies in a single HR-MAS CSI experiment should offer a valuable analytical approach for intact study of an organism. However, the induced centrifugal stress by MAS exerted onto the organism can inevitably disrupt the metabolic profiles (i.e. change in the metabolic activities or even in the morphology). A common strategy for subsiding the stress is by reducing the sample spinning. For example, spinning the sample at 500 Hz instead of 5000 Hz would reduce the force by two-order of magnitude. However, the issue is the presence of the dense spinning-sideband manifolds on the spectrum (denoted as ssb), which can complicate the analysis and, sometimes, impossible. Most efforts on slow spinning HR-MAS have been focusing on suppressing these ssb by (i) manipulating the rotor-phases in the pulse-experiments, 5 (ii) improving the water suppression efficiency, 6 (iii) utilizing the isotope 2 H detections, 7 (iv) optimizing the sample preparations, 8 and by (v) spectral manipulation. 9 This study intends to improve the HR-MAS CSI strategy for intact NMR spectroscopy by incorporating a slow sample spinning experiment together with the use of the integrated homospoil z'-gradient (for improving the feasibility). Both have yet to be explored. The experimental design is discussed and applied to a food phantom and an intact insect wasp (Vespula germanica). The detail descriptions of the experimental strategies and results will offer the basis for designing future intact studies with slow sample spinning HR-MAS CSI.

The integration of HR-MAS CSI pulse experiment with a ssb suppression is facilitated by the spin-echo nature in both experiments. It can be achieved by replacing the refocusing p pulse in CSI with a set of rotor-period composite p pulses (see Fig. S1). A major concern for metabolomics is the length of the experiment carrying out on susceptible specimens. This is especially of importance for a two-dimensional (2D) HR-MAS CSI with the chemical shift profile in F2 and the spatial profile in F1. As such, HR-MAS CSI can be a time demanding experiment. It depends on the number of t1 spatial encoding (denoted herein as t1CSI), and the detection sensitivity (i.e. metabolic abundance). Therefore, the choice of implementing the ssb suppression experiment must not dramatically lengthen the acquisition.

Total Suppression of Sidebands (TOSS) is a 1D ssb suppression experiment 10 and should be a beneficial facto to the integration; however, TOSS is highly susceptible to the imperfection of the p pulses and the instabilities of the spinning. Both can result in distorted spectral. 11 An extensive 243-steps phase cycling scheme are often necessary to overcome the effect of the pulse imperfections, 12 and has successfully applied to intact specimens. 6 The current preferred ssb suppression is a pseudo-2D Phase-Adjusted Spinning Sidebands (PASS), 13 with t1 associated to the signal evolution of the ssb (denoted herein as t1PASS, often with 8, 16 or 32 increments). Despite its 2D nature, PASS, in general, has better resistance to the pulse imperfection as compared to TOSS, 11 rendering a more effective approach with a less strenuous phase cycling scheme. 5 This includes an 11-steps cogwheel phase cycling scheme. 14 Fig. S2 shows the results of PASS with an 11-steps phase cycling.

This study incorporates the PASS experiment into the HR-MAS CSI (denoted herein as HR-PASS CSI) with an 11-steps cogwheel phase cycling scheme. HR-PASS CSI is essentially a pseudo-3D experiment by acquiring a series of individual 2D CSI spectra, each with an identical t1CSI spatial encoding but with different t1PASS steps. In this study, an 8-steps t1PASS is applied. As a result, the total acquisition time of HR-PASS CSI can be estimated with (11 ´ 8)•t1CSI. This is significantly shorter than that of an incorporation with TOSS, (243)•t1CSI.

In contrast to the previous HR-MAS CSI experiments, 2b,4 this study also explores the convenience of the experimental setup by exploiting the integrated homospoil gradient field Gz (along the magic angle z') of a conventional HR-MAS probe, without an external pulsed field gradient (PFG) system. In general, homospoil generates a weaker gradient field (and less durable) as compared to PFG; however, the field efficiency is compensated by the proximity to the sample. But one cautious note should be considered when applying the homospoil gradient. The duty cycle (the fraction of the pulse in one acquisition) should not exceed the specification limit. [START_REF]The upper limit for a Bruker HR-MAS probe is suggested to be of ~4 ms (δ) at 70 % (Gz ) of a 10 A amplifier[END_REF] As shown in Fig. S1b, rotor synchronization periods (tR) are carried out throughout the pulse experiment -the echo delays (TE), the 5p PASS pulses and the gradient pulse (δ). Consequently, this study limits the spinning frequency to 500 Hz and renders a fixed value of 2 ms for δ in the experiments. Subsequently, the field-of-view (FOV) in CSI is regulated by the gradient amplitude (GZ ) and the t1CSI encoding

(FOV = [t1CSI] × [k• 1H γ•δ•GZ ] -1 ,
where k is a scaling factor and γ is the 1 H frequency).

Food tissue is an excellent phantom to evaluate the capability of metabolic profiling with HR-PASS CSI owing to their uniform distributions of the metabolites in the tissue and, as necessary, it offers superior signal resolutions. The phantom consists of four distinct layers [zucchini (or courgette) -garlicapple -ground peanuts], each with a thickness of ~1 mm, placed transversely at the centre of the rotor. A single thin individual green leaf (of about 20 µm in thickness) is strategically placed between the layers of zucchini and garlic, and garlic and apple to evaluate the spatial separation by HR-PASS CSI.

As shown in Fig. 1a, both HR-MAS CSI and HR-PASS CSI exhibit a 2D spectrum with a spatial z' distribution (in mm) along F1, separating the individual sliced spectral profiles (in 1 H d ppm) in F2. Experimental details can be found in ESI. Both experiments were acquired with identical spatial encoding parameters, with a FOV of 8.3 mm attributed by a 64 t1CSI and 10 % of GZ . This results in a slice resolution DS (FOV / t1CSI) of 0.13 mm. As expected the acquisition time for HR-PASS CSI is 8 times longer (8 t1PASS) with a duration of two hours. However, based on the well-defined spatial profile (Fig. S3), the data acquisition could readily shorten by half (or more) with 32 t1CSI, which doubles the DS to 0.26 mm.

The HR-MAS CSI spectra (Figs. 1a and S4) reveal highly congested spectra with ssb signals masking some significant metabolites. For example, a large ssb appears at 6 ppm in the slice spectra associated with the garlic layer has overlapped with a set of multiplets of organosulfur allicin. In contrast, the results from HR-PASS CSI exhibit ssb-free spectra with isotropic metabolic signals. With one exception, a single ssb signal of the fatty acid methylene resonance is observed. This is ascribed from the large fatty acid content in the peanut layer. The stacked plot of the slice spectra in Fig. S5 shows distinct profiles among the different food layers. This is also supported by a principal component analysis. The superior peak resolution (with <0.02 ppm) allows for a rich metabolic identification for each layer (Tables S1-S2). It is noteworthy that a subsequent fast sample spinning (4000 Hz) HR-MAS CSI has been carried out (Fig. S5), and showed a dramatic change in the spectral profile for the individual layers, indicating a displacement of the metabolic fluids among the layers.

Exploiting the superior spectral quality data, a heatmap of the relative metabolic content (Rmet) along the z' is constructed (Fig. 1b). For each slice spectrum, Rmet can be estimated by the individual metabolic integral (Imet

): Rmet = [100% × Imet] / [NH-met × Ispectral sum],
where NH-met is the number of corresponding protons, and Ispectral sum is the total integral between 0.5-8.5 ppm excluding the water signal. This is followed by a colour coding scheme with red being the highest content and blue at zero. Subsequently, each colour box along z' corresponds to DS. The heatmap displays a few hidden variances that are not evident from the spectra. For examples, it reveals the different sugar distributions in the layers; allicin in the garlic layer; alanine, glutamine and fumarate in the zucchini layer; and the fatty acids associated within the thin layer of the green leaf (more discussion below). The z' projection of the heatmap (Fig. 1c) shows a greater extent of the metabolic distributions across the layers. Among the NMR observable metabolites, sugars comprise the highest -excluding fatty acids -metabolic content (ca. 1-2 %). The plot exhibits an apparent variation among the common natural sugars (sucrose, fructose and glucose) in the food tissues.

Moreover, it discloses that the apple layer consists of the highest total sugar level (i.e. sweetest among the layers) as compared to that of the garlic with the lowest. For fatty acids, the peanut layer, as excepted, consists of the highest content. However, a small -but visible -increased level of fatty acids can be located in between the food layers (zucchini | garlic and garlic | apple). Similarly, a subtle increased of choline is also found at the same z' location. These increased contents are associated with the green leaves between the layers. Since the acquired DS, 130 µm, is much greater than the leaf thickness ~20 µm, the observed metabolic signals (including the fatty acids and choline [START_REF] Augustijn | [END_REF] and others) are expected to overlap with the neighbouring profiles.

To illustrate the feasibility of HR-PASS CSI to a single intact organism, Figs. 2 andS6 show the resultant spectra of an insect wasp. The data were acquired with a DS of 0.26 mm and with a total acquisition time of 70 minutes. No visible body damage (or fluids) is found after the experiment. Experimental details can be found in ESI. The insect diameter (ca. 3.0 mm) is similar to the inner rotor diameter. This allows securing the insect with the longitudinal body axis aligned with the z', and permits for profiling the transverse plane of the distant anatomical sections of a wasp: tail, mid and head. The resultant HR-PASS CSI spectra were collected with a T2-contrast (TE = 100 ms). The benefits with a T2-contrast experiment on the wasp are shown in Fig. S6b. The long TE: (i) improves the spectral bassline by removing the broad signals of the macromolecules such as proteins; (ii) improves the water suppression and sequentially minimizes (or even eliminates) the residual ssb; and as importantly, (ii) enhances the metabolic signals by suppressing fatty acid signals (e.g. 0.86 ppm, 1.24 ppm, 1.59 ppm …) and ssb. It is noteworthy that if a T2-contrast experiment (instead of 8 ms) was applied to the food phantom above (Fig. 1a), it would probably subside the observed ssb of the methylene fatty acid.

About 20 metabolites are tentatively identified (Table S3), and a few resolved metabolic signals are exploited in the heatmap analysis (Figs. 2 andS7). The corresponding z' projection reveals the individual metabolic profile across the anatomical sections. The results show that the tail consists of a higher content of fatty acid as compared other sections. Interestingly, it also able to locate the venom region within the tail. Based on the heatmap analysis, the size of the venom region is about 1 mm along z'. It is known that venom contains three significant metabolites: 17 (1) amine-derivatives (such as histamine) involves in the inflammatory response by increasing the permeability of blood capillaries; (2) serotonin plays a role in irritation and contributes to pain; and (3) acetylocholine increases the perceived pain. Fig. S8 shows the presence of the acetylcholine and serotonin signals in the slice spectrum of the venom region. The unassigned peaks at 8.34 and 7.27 ppm could be related to the physiologically active amine derivatives in the venom.

In the mid-section, HR-PASS CSI reveals a high level of taurine and betaine (Fig 2). Interestingly, there is also a variance in the sugars. While, both glucose and fructose are significant in the head and tail sections, but only a hint of glucose is found in the mid-section. These observed differences could relate to the different in blood metabolisms 18 in the different anatomy. To evaluate the reproducibility of HR-PASS CSI experiment, an experiment with identical setup (spatial encoding and PASS evolutions) was performed on a second wasp individual (WASP #2). The heatmap analysis exhibits similar metabolic distributions across the anatomical sections (Fig. S7). Of course, more experiments and data are necessary to reaffirm the reproducible quantification of HR-PASS CSI profiling. Nevertheless, the results shown here offer a promising potential to metabolomic study of intact organisms.

The experimental strategy present here with HR-PASS CSI can be considered as an incremental -but significantimprovement to the previous HR-MAS CSI. It offers both spatial and chemical information on a conserved intact specimen. HR-PASS CSI could render a potent metabolic NMR profiling approach to organism research. This study provides the basis for designing intact NMR metabolomics. As importantly, it also offers the fundamentals for incorporating additional spinedited pulse experiments -such as diffusion or X-spin filtered -prior to HR-PASS CSI. Moreover, the PASS sequence could also be readily implemented in 1D slice localized spectroscopy [2b] for a direct profiling on a specific region of interest. This would permit a quick data acquisition and further enhances the conservation of the specimens, but the drawback is a prior knowledge of the spatial region of interest is necessary.

HR-PASS CSI could provide a new prospective for intact NMR-based metabolomics of research organism. The abilities of (i) unfolding the specific metabolic profiles from a heterogeneous specimen; and (ii) simultaneous profiling the different autonomic regions should provide a vital parallel investigation on the concurrent metabolic activities of interest, offering in-depth insights in life sciences.

All spectra data were acquired on a 500 MHz Bruker Avance II equipped with a conventional Bruker 4-mm HR-MAS probe. The detailed description of the experiment can be found in ESI. The study was supported by the French Agence Nationale de la Recherché (ANR-16-CE11-0023-01).

Conflicts of interest

Authors declare there are no conflicts to declare.

Notes and references

Experimental Details

Sample preparation

The samples were packing into a standard 4-mm ZrO 2 rotor with an internal spacer at the bottom. A top Kel-F insert was used to secure the samples. No buffer or D 2 O were added to avoid the displacement of metabolites during the sample spinning. Wasp: Three individual wasps was caught locally -within few minutes apart -and immediately snap-freeze and stored under -80 °C until the acquisition. The frozen wasp was placed entirely inside the rotor with the head at the bottom. The body dimensions of the three wasps were nearly identical with 3 mm in diameter. This allowed a tight fix inside the rotor preventing body movement during the spinning. Among the three wasps, one was sacrificed for setting up the experiments: from field shimming, testing the phase cycling schemes

NMR experiments

1 H experiments were performed on a Bruker Avance II spectrometer operating at 499. 16 MHz, with a three channel Bruker 4-mm HR-MAS probe integrated with a homospoil z'-gradient. The gradient was operated by a 40 A amplifier. The fast-and slow-sample spinning were carried out with MAS II at rate of 4000 and 500 Hz, respectively, with ± less than 2 Hz. NMR spectra were collected with 16 k time domain points, 20 ppm sweep width, 1 s low-power presaturation pulse as a recycle delay, and with an echo time (TE = 2 × N number of rotor period t R ) of 8 ms for the food phantom and 100 ms for the intact wasp. The number of accumulative scan was set to one-full phase cycling: 16 for HR-MAS CSI and 11 for HR-PASS CSI. In all CSI (including HR-PASS CSI) experiments, the gradient pulse length (δ) was set to a fixed value of 2 ms (n•t R ); subsequently, the desired field-of-view (FOV) was adjusted by the gradient steps t1 CSI together with the gradient amplitude G

Z (i.e. FOV = [t1 CSI ] × [k• 1H γ•δ•G Z ] -1 ,
k is a scaling factor which is calibrated experimentally with a rigid cylinder). All CSI experiments were performed with a FOV of 8.3 mm attributed by 64 t1 CSI and 10 % of G Z for the phantom; or 32 t1 CSI with 5 % for the wasp.

For the PASS sequence in HR-PASS CSI, 8-step t1 PASS was carried out. Consequently, a total of 8 consecutive 2D HR-PASS CSI spectra with identical CSI parameters (t1 CSI , δ and G Z ) were acquired, and each with a different set of interpulse delays for PASS. The delays can be found in Table 2 in For HR-PASS CSI, a summation of all 8 individual (i.e. 8-steps t1 PASS ) 2D FID was carried by a Topspin command 'adsu'.

• The resultant summed F2 FID (i.e. chemical shift profile) was processed with zero filling and with an exponential apodization (with a LB = 0.5 Hz). • The resultant summed F1 FID (i.e. spatial profile) was processed without 'zero' filling (TD = t1 CSI ) and with a Gaussian apodization (LB = -1 Hz; GB = 0.5). • A command 'xfb' in Topspin was applied to transform the summed 2D FID to the 2D spectrum with spatial profile in F1 and chemical shift in F2.

Each individual F2 slice projection spectrum was independently phased (except for those of wasp, which were phased with a magnitude mode), baseline correction 'abs', and referenced (to a doublet of alanine at 1.48 ppm). Peak assignments were assisted by Chenomx NMR suite 8.3, Human Metabolome Database Version 4.0 (www.hmbd.ca) and the corresponding references. Unfortunately, no supporting experiments (i.e. TOCSY, HSQC or J-RES) were performed to assist in the peak assignments.

Heatmap analysis

For each slice spectrum, R met can be estimated by the individual metabolic integral (I met ):

R met = [100% × I met ] / [N H-met × I spectral sum ],
where N H-met is the number of corresponding protons, and I spectral sum is the total integral between 0.5-8.7 ppm excluding the water signal. This is followed by a colour coding scheme with red being the highest content and blue at zero. The peak integration of the individual metabolic signal (listed in Table S4) was carried out manually using MestreNova v8.1. It should be noted that the deduced R met values should be considered as estimation because of the overlapping with the neighboring peaks.

Multivariate data analysis

The processed spectral data (by Topspin) were uploaded into MestreNova v8.1. The spectra (8.7 to -0.5 ppm, excluding the water signal) were aligned. Each spectrum was divided into 0.02-ppm buckets. The signal area under each bucket was integrated and normalized to a constant sum of all integrals. The resultant data in the matrix (sample × chemical shift bucket) correspond to observation and variable, respectively, of the model study. The matrix was exported to SIMCA P13 to proceed with the multivariate modeling. It was scaled using the Pareto scaling. Unsupervised PCA (principal component analysis) was performed. The models were validated by assessing the values of R 2 X, R 2 Y, and Q 2 . S1 andS2.

Fig. 1 1

 1 Fig. 1 1 H CSI spectra of a food phantom at 11.75 T with a sample spinning at 500 Hz and with a TE of 8 ms. Experimental detail can be found in ESI. (a) 2D contour plots of HR-MAS CSI and HR-PASS CSI with the z' spatial distribution (in mm) along F1 and the localized 1 H chemical shifts (in ppm) along F2. The F2 projection of HR-MAS CSI shows a highly congested spectrum with the ssb (highlight in red); however, HR-PASS CSI shows the contrary. (b) A heatmap showing the estimated relative metabolic content (Rmet) along z'. The description of Rmet is discussed in the text. Each colour box along z' corresponds to a DS of 0.13 mm with red represents the highest metabolic content and blue represents zero. (c) Heatmap z' projection of the individual metabolite showing the metabolic distributions across the different food layers. The shaded area in the sugar plot represents the sum Rglc + Rfru + Rsuc.

Fig. 2

 2 Fig. 2 (a) Representative F2 slice projection spectra of HR-PASS CSI (spinning at 500 Hz with a TE of 100 ms) for each anatomical section: green-tail; red-mid; blue-head. The projections are indicated by the dotted line in the 2D spectrum in Fig. S6a. Peak assignments are summarized in Table S3. (b) Heatmap showing the metabolic variations across the different anatomies. Each box along z' corresponds to a DS of 0 26 mm with red represents the highest metabolic content and blue represents zero. The z' projections are shown in Fig S7.

Food phantom :

 phantom The phantom (Fig 1a) with fresh local food layers was prepared by punching the individual pre-sliced layer (ca. 1mm) with the open-edge rotor.

  (Fig S2) to setting up the spatial encoding parameters and the different TE (Fig S6). Other two wasps (WASP #1 and #2) were used to evaluate the HR-PASS CSI experiment and its data reproducibility (Fig S7). Spectra data of WASP #1 were discussed in the text (Fig 2). It should be noted that no visible destructions (or fluids) on the body were found after the HR-PASS CSI experiment.

  [R.A. Wind et al., Prog. Nucl. Magn. Reson. Spectrosc., 2006, 49, 207.] Data acquisitions and processing were carried out with the Bruker software package Topspin 3.2 and MATLAB.

Fig S1 .

 S1 Fig S1. The NMR pulse sequence of (a) HR-MAS CSI and (b) HR-PASS CSI, along with the desired coherence transfer pathway which are determined by the phase cycling scheme: (a) j 1 = ({0°}*4, {180°}*4, {90°}*4, {270°}*4); j 2 = (0°, 90°, 180°, 270°); j R = (0°, 180°, 0°, 180°, 180°, 0°, 180°, 0°, 270°, 90°, 270°, 90°, 90°, 270°, 90°, 270°); and (b) 11-steps cogwheel scheme: j 1 = 90°; j 2 = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)  2p/11; j 3 = 0° and j R = (0, 6, 1, 7, 2, 8, 3, 9, 4, 10, 5)  2p/11. All interpulse delays are rotor synchronized (t R ), including the echo delay (TE), the gradient pulse and the ensemble composite 5p in PASS. The PASS interpulse delays can be found in Table 2 in [R.A. Wind et al., Prog. Nucl. Magn. Reson. Spectrosc., 2006, 49, 207.] The Bruker Topspin pulse sequences (HR-MAS CSI and HR-PASS CSI) are available upon request [alan.wong@cea.fr].

Fig S3 Phantom.

  Fig S3 Phantom. The F1 summed projection (i.e. 1 H density spatial profile). It shows the distant spatial profile for the different food layers along z'.

Fig S5 Phantom.

  Fig S5 Phantom. Comparison between slow (500 Hz, top) and fast (4000 Hz, bottom) spinning: (a) Stacked plots of the slice spectra; and its (b) PCA with R 2 X = 0.976, Q 2 = 0.94 for the slow spinning data, and with R 2 X = 0.993, Q 2 = 0.978 for the fast spinning data; and (c) F1 summed projection (i.e. spatial profile).

Fig

  Fig S6 Wasp. (a) 1 H HR-PASS CSI. A 2D contour plot showing the different anatomical sections. A significant venom metabolite, serotonine, has been identified in the tail section. (b) Spectral comparison of the different TE on a scarified intact wasp (i.e. Fig S2). It shows the effect of long TE (with 100 ms) for suppressing both the isotopic and ssb of H 2 O and fatty acids.

Fig S7 Wasp.

  Fig S7 Wasp. Heatmap comparison between wasp #1 (top) vs #2 (bottom). The data of wasp #1 are discussed in the text (i.e. Fig 2). The color lines in the z' projection corresponds to those labelled in the heatmap. Note that the venom metabolites (i.e serotonin) are located in the tail section.

Table S1 :

 S1 Tentative peak assignments from the individual HR-PASS CSI slice spectra of zucchini, garlic and apple

	Peak	[ppm] /	Metabolite	Zucchini	Garlic	Apple
	(multiplicity)
	0.94 (t)		Leucine (leu)	x	x	x
	0.96 (t)		Isoleucine (ile) x	x	-
	1.04 (d)		Valine (val)	x	x	-
	1.18 (t)		unknown	x	x	x
	1.33 (d)		Lactate (lac)	x	x	x
	1.48 (d)		Alanine (ala)	x	x	x
	1.91 (m)		Arginine (arg)	very low	x	-
	2.3 (t)		GABA	x	-	x
	2.14 (m), 2.45 (m)	Glutamine (gln) x	?	-
	2.38 (s)		Pyruvic	-	x	-
	2.37 (dd), 2.67 (dd)	Malic	?	-	?
	2.4 (s)		Succinate	?	-	-
	2.85 (s)		N(Me)3	-	x	-
	2.86 (dd), 2.96 (dd)	Asparagine	x	-	x
			(asn)
	3.36 (s)		methonal	x	x	x
	4.11 (m)		Fructose (fru)	x	x	x
	4.64 (d), 5.2 (d)	Glucose (glc)	x	x	x
	5.42 (d)		Sucrose	x	x	x
	5.54 (m)		Allicin	-	x	-
	6.53 (s)		Fumarate	x	x	-
	6.9 (d), 7.2 (d)	Tyrosine (tyr)	x	x	-
	7.3-7.4 (m)		Phenylalanine	x	x	-
	8.11 (s)		UMP	-	x	-
	Abbreviation: 'x' present; '?' questionable; '-' not present; (s) singlet; (d) doublet; (t) triplet;
	(m) multiple; GABA g-aminobutyric acid; UMP Uridine monophosphate
	Blue font: Peak used in the heatmap analysis

Table S2 :

 S2 Tentative peak assignments from the individual HR-PASS CSI slice spectra of peanut Peak used in the heatmap analysis. The shift range are summarized in Tab S4.

	Peak [ppm]	FA-peak
	0.88	-CH3
	1.28	-CH2-
	1.58	-CH2bCO
	2.02	-CH2-allyl
	2.25	-CH2aCO
	2.75	-CH2-bisallyl
	3.7	?
	3.83	?
	4.07	-CH2-glycerol
	4.28	-CH2-glycerol
	5.21	-CH-glycerol
	5.32	-CH = CH-
	Abbreviation: '?' questionable
	Blue font:

Table S3 :

 S3 Tentative peak assignments from the individual HR-PASS CSI slice spectra of head, mid and tail section (excluding the fatty acid peaks) Blue font: Peak used in the heatmap analysis. The shift range are summarized in TableS4.Reference:[1] Moreno M, Giralt, E. Three valuable peptides from bee and wasp venoms for therapeutic and biotechnological use: Melittin, apamin and mastoparan. Toxins (Basel). 2015, 7, 1126-1150; [2] Kapranas A, Snart CJP, Williams H, Hardy ICW, Barrett DA. Metabolomics of aging assessed in individual parasitoid wasps. Sci Rep. 2016, 6, 34848.

	Metabolite	peak (ppm)	Head	Mid	Tail
	Leucine (leu)	0.96		x	-	-
	Isoleucine (ile)	1.00		x	-	-
	Valine (val)	1.02		x	-	x
	Lactate (lac)	1.33		x	x	x
	Alanine (ala)	1.48		x	x	x
	Acetate (ace)	1.92		x	x	-
	Glutamate	2.03, 2.35		x	x	x
	Succinate (suc)	2.42		x	-	x
		3.08,	3.35,
	Serotonin*	6.86,	7.06,	-	-	x
		7.27, 7.37
	Choline (cho)	3.19		x	x	x
	Acetylcholine (ACh)	3.20, 3.73		-	-	x
	Phosphocholine (PC)	3.20		x	x	-
	Glycerophosphocholine (GPC)	3.22		x	x	x
	Taurine (tau)	3.26, 3.43		x	x	-
	β-glucose (β-glc)	4.65		x	x	x
	α-glucose (α-glc) / or trehalose	5.24		x	x	x
	Fructose	4.11, 4.01	4.04,	x	-	x
	Venom active animes (imidazole ?)*	7.27, 8.34		-	-	x
	Abbreviation: 'x' present; '?' questionable; '-' not present;
	* according to Chenomx Profiler assignment (Fig S9)

Electronic Supplementary Information

Peak assignment Table S1. Metabolites in the food phantom p12 Table S2. Fatty acids in the peanut layer p13 Table S3. Intact wasp (tail, mid, head) p14 Table S4. Chemical shift range used for the heatmap analysis p15


F Yamamoto

T Duc

O Gottschall

Z Mathon

Wang

S Chattopadhyay
email: s.chattopadhyay@hzdr.de

V Balédent

S K Panda

Sh Yamamoto

F Duc

T Herrmannsdörfer

M Uhlarz

T Gottschall

O Mathon

Z Wang

C Strohm

M Greenblatt

P Foury-Leylekian

J Wosnitza

4 f -spin driven ferroelectric-ferromagnetic multiferroicity in PrMn 2 O 5 under magnetic field

HAL is

I. INTRODUCTION

The research on functional materials, aiming at the next generation of smart devices, witnessed a massive upturn after the advent of materials showing multiple combined properties. Among such materials, an intriguing family of compounds known as magnetoelectric multiferroics (MEMF) with strongly coupled magnetism and ferroelectricity has attracted special attention [1][2][3] . For the MEMF family, reasons for being so highly sought after are twofold: (i) The microscopic origin of the magnetoelectric coupling is of fundamental interests that fascinates the condensed-matter community and (ii) from an applied perspective, the spintronics and data-storage technologies would greatly benefit from such functionalities.

However, to be well suited for device-applications, in addition to the room temperature functionality, MEMF materials are expected to fulfill two very important criteria: the coexistence of ferromagnetism with ferroelectricity and a strong coupling between the two order parameters. In reality, it is extremely difficult to satisfy both criteria simultaneously and only a few single-phase materials are known to show coupling between the ferroelectric and ferromagnetic component [4][5][6] . There are also a few theoretical works proposing strategies to couple ferroelectricity with weak ferromagnetism in a material 7,8 .

Multiferroic materials have been classified into two types. In type-I materials, ferroelectricity and magnetism have different origins. Whereas in type-II multiferroics, spin-driven ferroelectricity is caused by magnetic ordering itself, resulting an intrinsically strong magnetoelectric coupling 3,9,10 . However, the magnetism being frustrated in character, type-II multiferroics show the coexistence of complex antiferromagnetism and ferroelectricity. As multiferroics are relatively scarce, an important research activity focused on heterostructures with stacked ferromagnetic and ferroelectric layers 3 . However, as this eventually leads to type-I artificial multiferroics, the disadvantage of usually having small-coupling remains unsolved.

In this context, studies on RMn 2 O 5 (R = Bi, rare-earth) oxides are particularly interesting. The RMn 2 O 5 family is known for showing a series of magnetic transitions starting from R = Nd to Lu. The first magnetic transition around 40±5 K to an incommensurate antiferromagnetic (AF) state is followed by a second transition leading to a commensurate AF ordering of the Mn spins. An electric polarization (P) emerges either with the first (R = Tb, Sm, Nd) or the second magnetic transition along the b axis, asserting the type-II character of the multiferroicity [11][12][13][14][15][16][17][18][19][20] . The associated strong coupling has been reported for GdMn 2 O 5 and TbMn 2 O 5 where the electric polarization (P) can even be reversed by applying a modest magnetic field of a few Tesla 1,2,21,22 . With lowering the temperature, a third magnetic transition generally appears stabilizing another incommensurate AF order.

Neutron diffraction studies indicate that this spin-driven ferroelectricity is in general a consequence of a quasi-collinear ordering of the Mn spins (either Mn 3+ sublattice or both Mn 3+ and Mn 4+ sublattices) 11,12,14 . We note that, although the actual room-temperature crystal structure of RMn 2 O 5 is already polar with a Pm (monoclinic) space group 23 , the structural distortion with respect to the average non-polar Pbam space group 24,25 is so small that the induced polarization becomes extremely weak and cannot be measured directly in this regime. We will thus neglect this weak symmetry breaking in the following.

In RMn 2 O 5 , the long-standing debate on the microscopic origin of the spin-driven ferroelectricity has been resolved recently confirming Mn-Mn exchange-striction to be the responsible one 12 . This result portrays the dominant role of the 3d ions and their frustrated superexchange interactions in the emergence of the spin-driven ferroelectricity. Although known for a few other multiferroics, 4,6,26 in RMn 2 O 5 , the role of the 4 f ions has only been revealed recently to explain the magnetoelectric behaviors in some candidates 27,28 . Particularly in GdMn 2 O 5 and NdMn 2 O 5 , in addition to the dominant 3d-3d effect, a weak 3d-4 f interaction was proposed to explain the observed spin-driven ferroelectricity 15,29 .

Among all compositions, PrMn 2 O 5 appears to be an outstanding exception. In contrast to the other RMn 2 O 5 members with non-zero numbers of 4 f electrons, PrMn 2 O 5 does not show any spin-driven ferroelectricity. Powder neutron diffraction measurements show that Mn 3+ moments order at T N1 = 25 K, following a magnetic propagation vector q 1 = (0.5, 0, 0). With decreasing temperature, Mn 4+ ordering appears at T N2 = 18 K with q 2 = (0, 0, 0.5) 30 . Such distinct orderings of the Mn 3+ and Mn 4+ ions are indicative of a very weak exchange coupling between the two sublattices, explaining the absence of ferroelectricity. Pr 3+ sublattice does not fully order down to 1.5 K. Only a partial ordering of Pr 3+ below T N1 was reported 30 suggesting a coupling between Mn 3+ and Pr 3+ . A powder neutron study under high pressure (a few GPa) shows the emergence of a collinear magnetic phase favorable for spin-driven ferroelectricity 31 . However, electric polarization measurements were not possible to perform at this pressure regime to probe the possible onset of multiferroicity directly.

In this work, we report the emergence of a multiferroic phase under magnetic field in PrMn 2 O 5 . Using the combination of high-field electric polarization, x-ray magnetic circular dichroism (XMCD), and density functional theory (DFT) based calculations, we show that unlike the multiferroicity observed in other RMn 2 O 5 members, the spin-driven ferroelectricity in PrMn 2 O 5 originates from and couples to a ferromagnetic component. Moreover, the associated mechanism involved is no longer based on the 3d-3d coupling, rather the spin-driven ferroelectricity is solely a manifestation of 3d-4 f exchange interaction.

As mentioned, the presence of coupled ferroelectric and ferromagnetic components is very rare to find. In addition to that, PrMn 2 O 5 hosts 3d-4 f coupling and exchangestriction mechanism as well. As separate phenomena, these features have been observed. However, to have a material like PrMn 2 O 5 , which hosts simultaneous presence of all of these effects/mechanisms is rather unique to our knowledge.

Single crystals of PrMn 2 O 5 from the same batch as mentioned in reference 23 were used for this study. The crystals were grown using electrolysis method as described in references 32,33 . The as-grown crystals have a thin plate-like morphology with the plate-surface being perpendicular to the baxis.

II. RESULTS & DISCUSSION

A. High-field Electric Polarization & Magnetization

We carried out magnetic-field-dependent electricpolarization (P) measurements up to 45 T using a pyroelectric technique at the Dresden High Magnetic Field Laboratory (HLD-EMFL). Field pulses of ∼20 ms duration were applied along b to measure the spin-induced pyroelectric current (I) along the same direction. A schematic diagram of the measurement technique is shown in the inset of the top panel of Fig. 1. Since the sweep rate of the magnetic field (dH/dt) is large in pulsed-field, the dP/dt and hence the pyroelectric-current becomes detectable even for a small change in P. The field induced pyrocurrent (I) was recorded by measuring the voltage variation across a shunt resistor (R S ). This shunt resistor was connected in series with the measurement circuit by a digital oscilloscope (Yokogawa DL750). The oscilloscope was operated with a high sampling rate of 1 MSs -1 and a resolution of 16 bit. The top panel of Fig. 1 shows the field dependence of I measured at 1.5 K, 7 K, and 22 K. The electric polarization was then calculated by integrating the I(H) curves. The reproducibility of the data was verified carefully by repeating the measurements multiple times.

From the P(H) data shown in the bottom panel of Fig. 1, it is evident that at 1.5 K, an electric polarization emerges above ∼12 T along the b direction. P increases with H and attains a maximum around 27 T. With further increase in field, P decreases slowly and finally enters to a flat region above 35 T. Notably, the polarization amplitude in PrMn 2 O 5 (∼1 nC/cm 2 ) is of the same order as found in NdMn 2 O 5 , the adjacent member. The polarization gradually becomes weaker at higher temperatures and finally vanishes in the paramagnetic regime.

We also performed high-field magnetization measurements using pulsed magnetic fields at the HLD. Similar to the polarization measurements, the pulse duration was ∼20 ms. We repeated this measurement a few times to ensure the reproducibility of the data. Figure 2 shows M(H) data measured at 1.5 K up to 42 T with H applied parallel to the b direction. In contrast to the polarization data, no anomaly is seen in M(H) over the entire field range. The magnetization does not show any indication of saturation even at the highest field measured. In this system, the overall magnetization is dominated by the Mn moments. Therefore, the observed dissimilarities between the field dependence of the polarization and magnetization indicate that the magnetic-field-driven ferroelectricity does not seem to involve 3d ions directly. Rather, it indicates a possible 4 f ion involvement in this case.

B. High-Field X-Ray-Magnetic-Circular-Dichroism

In order to investigate the effect of a strong magnetic field on the rare-earth ordering, we performed x-ray-magneticcircular-dichroism (XMCD) measurements. The high-field XMCD measurements at the Pr-L 2 edge were performed at the energy dispersive x-ray absorption spectroscopy beamline ID24 at ESRF, Grenoble 34 . The pulsed-field magnet used for this purpose was connected to a 1.15 MJ portable power sup- ply developed at the LNCMI-Toulouse. The magnet produces a maximum field of ∼30 T with a rise time of ∼10 ms and a total duration of ∼23 ms (Fig. 3) in every 8 min [35][36][37] . Figure 4(a) shows a representative x-ray absorption spectrum (XAS) recorded at the Pr-L 2 edge at 2 K. XMCD spectra were obtained in transmission mode at 2 K as the difference of x-ray absorption spectra when changing the field direction (+b andb directions). The same measurement protocol was applied for both right-and left-handed circularly polarized x-ray beams.

To record the field dependence of the absorption spectrum, a multi-frame acquisition scheme with high-frame-rate detector FReLoN (Fast-Readout Low Noise) was used and a series of 50 full-energy spectra (i.e., acquisition windows of 1 ms) were recorded during each field pulse 38 . For our purpose, the crystal was mounted with the b axis parallel to the magnetic field and incident beam direction. The polished sample (thinned down to ∼15 µm) was sandwiched between two diamond windows. This assembly was mounted in a dynamic He-flow cryostat in which the sample was cooled down through forced convection.

Figure 4(b) shows XMCD spectra recorded at magnetic fields between ∼4.6 T and ∼30 T. The data presented here is a summation of 16 field-pulses to improve the signal-tonoise ratio. The main contribution at the L 2 edge comes from the 2p 1/2 → 5d 3/2 dipole transition according to the selection rules. As the magnetic field increases, the XMCD amplitude becomes more pronounced. The observed x-ray absorption spectra and dichroic spectra at the high-field regime are consistent with those reported for other Pr 3+ -based compounds (4 f 2 ) 39,40 .

The dichroic signal starts to appear about 10 T, reaches its maximum around 27 T, and then drops gradually with further increase in H. Although we cannot entirely exclude a contribution from the Mn ions, the Pr-L 2 edge XMCD signal observed in PrMn 2 O 5 is mainly caused by the Pr spins. The direct evidence can be seen from the peak shape of the XMCD signal. It does not change with field, only its amplitude gets modified unlike garnet oxides, 41,42 where the L edge XMCD line shape as well as the amplitude change with field due to contributions of comparable strength from both 4 f and 3d ions. In case of intermetallic compounds composed of a 4 f rare-earth and a 3d transition-metal, earlier study reported a sizable transition-metal contribution as well to the rare-earth L 2 edge 43 . This is because the 3d-4 f exchange interaction in such intermetallic compounds is mediated by the rare-earth 5d band which is directly accessed by the rare-earth L edge. Instead, in our system, the weaker super-exchange 3d-4 f interaction is mediated via the oxygen 2p state, which results in a marginal transition-metal contribution to the XMCD signals at the rare-earth L edges. Here, the mixed 4 f -5d state of Pr 3+ is the reason the transition toward 5d 3/2 gives access to the 4 f spin's ferromagnetic ordering. As a consequence, these results prove that the Pr ions contain a ferromagnetic component along the b axis under an external magnetic field applied along the same direction.

Notably, the field dependence of the integrated XMCD signal has a remarkable resemblance with the aforementioned polarization curve. In Fig. 5 we show the field dependence of both integrated XMCD signal and polarization at comparable temperatures (2 K for XMCD and 1.5 K for P). The error bars shown are standard deviations extracted from averaging procedure for 16-polarization-dependent spectra. The similarity between both quantities suggests that the two order parameters are coupled. The direct consequence is that PrMn 2 O 5 becomes a type-II multiferroic under magnetic field with coupled ferromagnetism and ferroelectricity. Above ∼27 T, the steeper decrease of the XMCD integrated intensity as compared to the P(H) behavior might come from the fact that the fluctuations of the signal at high magnetic field are larger, as reflected through the larger error bars when approaching 30 T. Therefore, the slope of the curve should be taken as an overall trend of the field dependence without considering it with high precision.

C. Density-Functional-Theory Calculations

In order to provide further credence to our experimental findings, the density-functional-theory (DFT) calculations were performed in the generalized-gradient-approximation + Hubbard U (GGA+U) approach by means of a full potential linearized muffin-tin (MT) orbital method (FP-LMTO) 44,45 as implemented in the RSPt code [START_REF] Wills | Electronic Structure and Physical Properties of Solids: The Uses of the LMTO Method[END_REF] . The Brillouin-zone (BZ) integration is carried out by using the thermal smearing method in a 8 × 7 × 10 k-mesh which corresponds to 560 k-points in the irreducible part of the BZ. For the charge density and potential angular decomposition inside the MT spheres, the value of maximum angular momentum was taken equal to l max = 8. To describe the electron-electron correlation within the GGA+U approach, we used U = 4 eV and J = 0.8 eV for the Mn-d states and U = 5 eV and J = 0.5 eV for the Prf states.

The calculated projected density of the states (PDOS) are displayed in Fig. 6 which shows that the system is insulating with a gap of 1.2 eV. As expected, the Prf states are highly localized. We observe that the majority spin channel of Prf states shows a valence band peak at around -1.3 eV binding energy and rest of the states appear above the Fermi level between 3 to 5 eV, while the minority spin channel is completely empty. Both the spin moment on the Pr-site and the projected density of states are consistent with the nominal 4 f 2 state of the Pr ion. On the other hand, the occupied Mn-d states are fairly delocalized within the binding energy range from -7 eV to the Fermi level. For both Mn 4+ and Mn 3+ ions, the majority channel is partially filled, while the minority channel is completely empty. The PDOS and the projected moments are consistent with the tetravalent (3d 3 ) and trivalent (3d 4 ) states of Mn 4+ and Mn 3+ , respectively.

The obtained spin and orbital angular momentum of the Pr 3+ ions are 1.81 µ B and 1.99 µ B , respectively and they are in opposite directions. The spin moment on Pr 3+ is large as expected from the nominal charge state ( f 2 ). However, the presence of a strong reverse orbital moment makes the net moment very small (0.18 µ B ). It provides the explanation for the observed low moments on Pr site in neutron experiments (∼0.5 µ B /Pr 3+ ) 30 . The moment of Mn 3+ and Mn 4+ are essentially spin moments with the estimated values 3.38 and 2.55 µ B , respectively.

Further, to establish a spin model to understand the field dependence of the electric polarization and its relation with the long range magnetic ordering, we estimated the interatomic magnetic exchange interactions from the converged GGA+U calculations using the formalism of Ref. [START_REF] Kvashnin | [END_REF] . We used the magnetic-force theorem 48,49 to extract the effective intersite exchange parameters (J i j ). In this method, we mapped the total converged energies of the magnetic system onto the following Heisenberg-type spin Hamiltonian

Ĥ = - i j J i j S i • S j . (1)
Here, the indices i and j span over the positions of the intrinsically magnetic ions, i.e. Pr 3+ , Mn 3+ , and Mn 4+ . The effective J i j is extracted in a linear-response manner via a Green's function technique. A detailed discussion of the implementation of the magnetic force theorem in RSPt is provided in Ref. 50 . All the calculations are carried out using the structural parameters given in Ref. 30 . The magnetic structure determination from neutron diffraction data shows that the magnetic moments of Pr 3+ are parallel to their nearestneighbour Mn 3+ spins and they make only a small angle with respect to the a axis. Therefore, in the first approximation, we can safely ignore the anisotropic terms in the spin-Hamiltonian. That is precisely the reason, we have only considered the isotropic Heisenberg exchange term in our model spin-Hamiltonian.

The estimated exchange interactions are listed in Table I. These results are fully compatible with the magnetic structure reported in 30 . The interactions J 3 to J 6 , which play a significant role in the present context, are presented in Fig. 7 along with the zero-field magnetic structure. The strongest interaction is J 5 , imposing a perfect antiparallel alignment of the Mn 3+ spins (represented in blue), along the direction of anisotropy. The Pr moments (in yellow) are only coupled to those Mn 3+ through a ferromagnetic J 6 . Owing to the perfect colinearity of the Mn 3+ and Pr 3+ spins, it is reasonable to assume that the Pr anisotropy energy is negligible compared to J 6 . Now, it is understandable why the Mn 4+ (in green) moments are not coupled to the Mn 3+ spins. Indeed, as can be seen in Table I, J 4 is almost two orders of magnitude smaller that J 5 , excluding this path to connect the two sub-lattices. The other possible connection between Mn 3+ dimers goes via J 3 which is nearly an order of magnitude smaller than J 5 . This interaction connects Mn 4+ to two antiparallel Mn 3+ spins. As a consequence, the two paths connecting Mn 3+ and Mn 4+ via J 3 exactly compensate each other. As for J 1 and J 2 , they couple two Mn 4+ moments along c. Symmetries of the magnetic structure are shown in Fig. 7(a), resulting the magnetic space group P a b2 1 a. This space group is the same as the one reported for other members of the series, such as GdMn 2 O 5 15 . From the experimental results 30 , we know that there are two distinct magnetic sublattices having separate ordering temperatures. The first one involves Mn 3+ and Pr 3+ spins and the second one is formed by the Mn 4+ spins. The fact that these two sublattices do not order at the same temperature, indicates that there is almost no coupling between the two. That is the reason why we considered the Pr 3+ -Mn 3+ interaction J 6 only TABLE I. The Mn-Mn and Mn-Pr exchange interactions (in meV) obtained from the converged GGA+U calculations using the magnetic-force theorem 51 . The negative values imply antiferromagnetic interaction.

J i Ions involved Distance (Å) Magnitude (meV) J 1 Mn 4+ -Mn 4+ 2.95 -0.110 J 2 Mn 4+ -Mn 4+ 2.74 0.182 J 3 Mn 3+ -Mn 4+ 3.39 -0.040 J 4 Mn 3+ -Mn 4+ 3.61 -0.004 J 5 Mn 3+ -Mn 3+ 3.61 -0.205 J 6 Pr 3+ -Mn 3+ 3.35 0.009

and not the Pr 3+ -Mn 4+ one. As for the J 6 interaction, every Pr 3+ ion is connected to two Mn 3+ ions: one at 3.35 Å distance and another, at 3.80 Å distance. Our calculations show that the Pr 3+ -Mn 3+ exchange corresponding to the distance of 3.80 Å is negligibly small due to the very large distance. In view of that, we have provided only the value of J 6 in the Table I which actually corresponds to the Pr 3+ -Mn 3+ distance of 3.35 Å. Interestingly, the Pr 3+ -Mn 3+ interaction J 6 is ferromagnetic as opposed to the antiferromagnetic Gd 3+ -Mn 3+ interaction in GdMn 2 O 5 15 . The primary reason of this could be attributed to the difference in the fillings of the f -orbitals of these compounds. According to the extended Kugel-Khomskii model 52,53 , the nature of the inter-atomic magnetic interaction primarily depends on three important parameters, namely (i) crystal-field splitting, (ii) effective hopping strengths between the relevant orbitals, participating in forming local magnetic moments, and (iii) the nominal fillings of those orbitals which determine whether the virtual hopping is allowed between them depending on their parallel or antiparallel alignment. The differences in the structural parameters of GdMn 2 O 5 and PrMn 2 O 5 will result slight differences in the first two parameter. However the most significant difference comes due to the third point: the nominal occupancy of Pr 3+ is f 2 , while for the Gd 3+ it is f 7 , which is exactly half-filled. It is well established that half-filled orbitals promote antiferromagnetic superexchange since virtual hopping is allowed only if they possess antiparallel alignments, making AFM ordering between Gd 3+ and Mn 3+ energetically favourable. Whereas for Pr 3+ , f orbitals are less than half-filled ( f 2 ) and thus virtual hopping between Pr 3+f and Mn 3+ -d are allowed for both parallel and antiparallel alignments and resulting exchange comes out to be ferromagnetic. A detailed calculation based on the Kugel-Khomskii model using calculated hopping and onsite energies is outside the scope of this work.

D. Origin of the Field-Induced Ferroelectricity

On the basis of the estimated exchange interactions, we propose a possible model to couple the experimentally observed ferromagnetic component of Pr 3+ to the induced ferroelectricity along the b direction. As mentioned, the essential components responsible for the magnetic order are J 5 connecting two Mn 3+ spins and J 6 , connecting Mn 3+ and Pr 3+ spins. The amplitudes of the arrows showing the magnetic moments are not to the scale for the sake of clarity. Mn 3+ , Mn 4+ , and Pr 3+ ions are represented in blue, green, and yellow, respectively.

Upon increasing the external magnetic field along the b direction, both Mn 3+ and Pr 3+ spins are expected to align along b. Due to the strong J 5 interaction, Mn 3+ will persist as antiparallel dimer. This dimer will start rotating towards b as soon as the magnetic field exceeds the anisotropy energy of Mn. For the Pr 3+ moments, there is a competition between two opposite interactions. On the one hand, the ferromagnetic interaction J 6 tries to align the Pr 3+ spins parallel to the neighboring Mn 3+ spins. As a consequence, Pr moments get antiferromagnetically aligned like the Mn 3+ pairs. Now comes the Zeeman interaction term induced by the external magnetic field. It will push the Pr 3+ spins to align along the b direction.

As the XMCD results show, this happens above ∼12 Tesla, when a ferromagnetic component from Pr moments appears along b. A compatible magnetic structure preserving the maximum of the zero-field symmetries is schematically proposed in Fig. 7(b). The magnetic symmetry operations (represented on the magnetic structure) correspond to the magnetic space group Pb 2 1 m (Pm c 2 1 #26.70 in conventional setting), a subgroup of P a b2 1 a. As a result, the two Mn 3+ -Pr 3+ pairs do not remain equivalent anymore. A slight displacement of the oxygen ions bridging the Mn 3+ and Pr 3+ is thus to be expected and would be different from one pair to the other. For the nearly ferromagnetically aligned Mn 3+ -Pr 3+ pair, this dis-placement will tend to maximize the non-frustrated J 6 interaction. On the contrary, the oxygen displacement for the nearly antiferromagnetic Mn 3+ -Pr 3+ pair will tend to minimize this J 6 exchange coupling. Displacements are thus different from one side to the other resulting an effective polarization along b, as observed experimentally. From a structural point of view, the average centrosymmetric space group Pbam at zero field is no longer compatible. The high-field space group is thus expected to be a non-centrosymmetric subgroup of Pbam and compatible with the Pm c 2 1 Shubnikov group. Among the two possible subgroups fulfilling the first constraint (Pba2 and Pmc2 1 ), only Pmc2 1 allows a polarization along the b direction. Fortunately, this group also satisfies the second requirement since it is a subgroup of Pmc2 1 1 . As a conclusion, it is most likely that the high-field space group is Pmc2 1 . This mechanism, leading to the coupling of a ferromagnetic Pr component and ferroelectricity, is thus mediated purely by a 3d-4 f super-exchange interaction. This is fundamentally different from the exchange-striction mechanism involved in the other members of the RMn 2 O 5 family where a dominant 3d-3d interaction leads to the coupling between antiferromagnetism and ferroelectricity.

According to our DFT calculation and neutron diffraction results reported earlier 30 , in PrMn 2 O 5 , the combined moment of Mn 3+ and Mn 4+ is much larger than the Pr 3+ . In the fielddependent-magnetization shown in Fig. 2, only the b components of these moments contribute. The weak nature of the ferroelectricity and the measured XMCD signal signifies that the field-induced ferromagnetic component of Pr 3+ spins along b, responsible for the breaking of centrosymmetry leading to the emergence of ferroelectricity, is rather a small fraction of its already tiny total moment of about 0.18 µ B . That is why in the bulk magnetization measurement, where Pr 3+ , Mn 3+ , and Mn 4+ contribute together, any possible characteristic contribution from Pr 3+ spins in the vicinity of the fieldinduced polarization, was not visible in presence of a much larger background caused by the Mn moments. In contrast, an element selective technique like XMCD was able to detect successfully even the small ferromagnetic component of Pr 3+ along the b direction. In this context, it is also important to clarify that Fig. 7(b), which depicts our model to explain the field-induced ferroelectricity, is an exaggeration of the real scenario. To explain the symmetry breaking mechanism with better clarity, we have chosen to show the emergence of the field-induced ferromagnetic component of the Pr 3+ spins in an amplified way compared to what one expects in reality for the present situation. This is only for the sake of better visualization of the proposed model.

III. CONCLUSION

In conclusion, the combination of electric polarization and XMCD measurements in high magnetic fields along with DFT calculations reveals the emergence of spin-driven ferroelectricity in PrMn 2 O 5 under magnetic field with strong magnetoelectric coupling. In contrast to other RMn 2 O 5 members (R = Nd to Lu), multiferroicity in PrMn 2 O 5 under magnetic field is characterized by coexisting ferroelectric and ferromagnetic components. The underlying mechanism for this spindriven ferroelectricity involves an exchange-striction mechanism solely originated from 3d-4 f coupling as opposed to the 3d-3d dominated multiferroicity in the other members. The observation of such a coupled ferroelectric-ferromagnetic state opens up new perspectives for technological applications. The present study evidences that there exists the possibility to stabilize a robust ferroelectric-ferromagnetic combination along with strong magnetoelectric coupling by manipulating the magnetic frustration using external parameters such as magnetic field and pressure.

ACKNOWLEDGMENT

We acknowledge the attribution of ESRF beam time under proposal HC-3783. We also acknowledge the support of the HLD at HZDR, member of the European Magnetic Field Laboratory (EMFL), and the DFG through SFB 1143 and the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matterct.qmat (EXC 2147, Project No. 390858490).

FIG. 1 .

 1 FIG. 1. Top panel shows magnetic field (H) dependence of the pyrocurrent (I) up to 45 T at different temperatures using H I b configuration. Inset shows schematic diagram of the pyroelectric technique used to record the field dependence of the pyrocurrent applying pulsed magnetic fields. Bottom panel depicts electric polarization (P) along the b axis as a function of magnetic field.

FIG. 2 .

 2 FIG. 2. Field dependence (H) of the magnetization (M) up to 42 T at 1.5 K with H b.

FIG. 3 .

 3 FIG.3. Time profile of the magnetic field pulse used for the x-ray magnetic circular dichroism measurements. The red symbols correspond to the average values of the magnetic field and error bars over a series of 10 pulses for each data-acquisition window of 1 ms.

FIG. 4 .

 4 FIG. 4. (a) X-ray absorption spectrum (XAS) and (b) x-ray magnetic circular dichroism (XMCD) spectra at the Pr L 2 -edge for PrMn 2 O 5 with magnetic field applied parallel to the b-axis at T = 2 K. Data sets for different magnetic fields are offset along the vertical axis for clarity.

FIG. 5 .

 5 FIG. 5. Field dependence of the XMCD integrated intensity over the entire Pr-L 2 edge range at 2 K compared with the field dependence of the electric polarization along the b axis at T = 1.5 K shown in Fig. 4. The dotted line is a guide to the eyes.

FIG. 6 .

 6 FIG.6. Theoretically calculated partial density of states for Prf and Mn-d in the magnetic ground state.

FIG. 7 .

 7 FIG. 7. (a) Schematic magnetic structure at zero field (magnetic space group P a b2 1 a) 30 . (b) Proposed schematic magnetic structure at high field in the ferroelectric phase (magnetic space group Pb 2 1 m ).The amplitudes of the arrows showing the magnetic moments are not to the scale for the sake of clarity. Mn 3+ , Mn 4+ , and Pr 3+ ions are represented in blue, green, and yellow, respectively.


M Le Ny

O Chadebec

G Cauffet

J M Dedulle

Y Bultel

A Three Dimensional Electrical Model of PEMFC Stack

Keywords: Cell Interaction, Fault Modeling, Internal Loop of Current, PEMFC Stack Modeling

Introduction

Proton exchange membrane fuel cells (PEMFC) are energy conversion devices that produce electricity and heat directly from hydrogen and oxygen. Nevertheless, in spite of undeniable progresses undergone both in scientific knowledge and technological development, PEMFC assemblies and parts still present certain drawbacks and improvements must be made such that the system is reliable for industrial applications. Reliability and durability appear to be the most important considerations to successfully achieve the commercialization of proton exchange membrane fuel cell (PEMFC). Fuel cells are prone to chemical, mechanical, or thermal degradation [1] and three degradation classes could be distinguished: "baseline degradation", as well as degradation due to transients and incident induced degradation. All these degradation classes induce fuel cell voltage/performance decrease and lifetime reduction. The baseline degradation is due to longterm material degradation. It is irreversible, unavoidable, and exists as long as the fuel cell is operating. The second and third degradation classes are linked to the variations of operating conditions or detrimental operating conditions. Such degradations are induced by temperature, water management, or under faulty control of operating parameters, such as reactant starvation which refers to an operation with an insufficient amount of gas in the active layer.

In the ideal situation, cells in a stack have to perform uniformly. Unfortunately, this cannot be observed experimentally. The reason for the non-uniform operation of the cells in the stack is still not well understood. Significant variations among the cells may accelerate long-term structural problems and may eventually lead to a total failure of the power generation unit. Possible causes are non-uniform fuel/air flow distribution to individual cells, non-uniform temperature, and current distribution within the stack, and material nonuniformities (compositions and microstructure).

Generally, the modes of membrane electrode assembly (MEA) degradation can be divided in three categories: cata--

[ * ] Corresponding author, mathieu.le-ny@g2elab.grenoble-inp.fr

Abstract

Mathematical modeling is an essential tool in the design of fuel cell systems, as it is important to understand the response of a stack under different conditions. This paper proposes a three dimensional electrical model of PEMFC stack. Such approach highlights the phenomena that occur at the macroscopic scale which are the electrical interactions taking place within the stack. At this scale, the microscopic phenomena are taken into account by averaging them over all the thickness of the membrane electrode assembly (MEA). This model consists of solving the charge transport equation in three dimensions. In order to avoid interfacial issues with the potential jump that occurs in the MEAs, a special source term is added in the transport equation. It makes the model robust, reliable, and helpful for understanding the effect of anomalies on the electric behavior of a PEMFC stack. After describing the equations and the numerical method used for the model, the paper shows some electrical phenomena such as cells interactions, non-equipotential bipolar plate, higher cell voltages and internal loop of current. lyst activity loss, carbon support porosity loss, and membrane conductivity loss. The first is caused primarily by loss of the platinum surface area through dissolution/precipitation platinum during voltage cycling. Additionally, catalyst activity can degrade through poisoning. Damage to the carbon support occurs through corrosion during high potential excursion that is generally caused either by hydrogen-air front on the anode during shutdown/start-up or localized fuel starvation. Finally, it is conceivable that membrane can undergo chemical degradation. All this degradation processes as well as MEA delamination would be manifested as an increase of the MEA resistance. As well as MEA degradation, membrane drying, or electrode flooding can also induce an increase of the cell resistance.

Flow-field channels obstruction by liquid water and fuel cell design anomaly can both induce uneven gas distribution or even local insufficient supply of reactants. Another cause of improper gas supply could be a low reactant flow rate aiming at maximizing reactant utilization in order to optimize fuel cell efficiency. In this situation, coming from oxygen fed in the fuel cell's cathodic compartment (respectively from hydrogen fed in the fuel cell's anodic compartment), oxygen (respectively hydrogen) molecules could diffuse due to the concentration gradient towards the anode (respectively cathode) where they accumulate. In addition to the problem associated with sudden load increase, start-ups and shut-down might be affected by the possible presence of air inside the anodic compartment, particularly after a prolonged shutdown. Finally, in addition to chemical degradation, mechanical degradation of membrane can occur due to water management. Cycling relative humidity causes the membrane to expand and contract by absorbing varying quantities of water. Both chemically and mechanically degraded membrane can induce localized pinhole formation. As a consequence, fuel or oxidant gas starvation and pinhole formation can induce a local decreasing of the electromotive force which can suddenly tend to zero.

PEMFC modeling is a powerful tool in improving the durability and the reliability. In many models, cells are often studied alone, regardless of the surrounding elements. But in reality, cells are integrated in large stacks and interact with their neighbor cells. Up to date, although lot of work has been devoted on single fuel cell modeling, less attention has been paid in the literature to the behavior of a whole fuel cell stack. In Ref. [2], Park et al. describes a 20-cell stack model in which each cell operate uniformly. In order to study the fluidic interactions, Karimi et al. propose a hydraulic network approach in Ref. [3], whereas Burt et al. developed a pseudo 2D model using Finite Volume Method in Ref. [4]. There is also a stack model using partial system approach [5]. In Refs. [START_REF] Gerard | Proceeding of FDFC Conference Nancy[END_REF][START_REF] Gerard | [END_REF], Gerard et al. present a multi-physic model with a complex 3D discretization of a stack. If some of these models highlight the fluidic and thermal interactions between cells, none of them takes into full account the electrical phenomena. Bipolar plates are not modeled or are considered equipotential. These approaches make impossible the study of many of the electrical interaction that can occur between the cells of a same stack. For diagnostic application and stack stabilization purpose, Hirschfeld et al. [8] developed a 3D electrical model of PEMFC stack using the Finite Volume Method and whose equations in the MEA are more detailed in Ref. [9]. This approach is really useful to see the effects of resistive faults inside a stack and to study the influence of aluminum plates that have stabilization purpose. Nevertheless they do not seem to take into account the potential jump that occurs in the MEAs. Even if this model is capable of predicting the effect of many resistive faults, getting the right potential distribution remains always impossible as well as the phenomena associated to inhomogeneous distribution of electromotive forces. Finally, in Refs. [10,11], Berg et al. developed a 1 + 1 dimensional electrical and fluidic model of PEMFC stack using their own numerical approach. Currently, it is the more advanced PEMFC stack model taking electrical coupling into account. But its numerical approach makes it difficult to reproduce.

The aim of this paper is to present an electrical model of PEMFC stack that provides the three dimensional steady state distribution of the potential and current density for various cells' configurations and operating conditions. In this model, the electrical phenomena at the stack's level are highlighted instead of the electrochemical and mass transport processes that take place at microscopic scale as it is done usually. The contributions of the anode, the cathode, and the membrane are not distinguished. All the physical phenomena that take place in these different layers are approximated and averaged over the MEA thickness. Fluidic, thermal transfer, and dynamical behavior are not yet taken into account but should be incorporated in our code in future work. The Finite Volume Method is used here because of some of its advantages. This technique allows some modularity with the definition of the problem (geometry, mesh, and physics) and involves only flux balance equations which are equivalent to the Kirchhoff's nodal rule. This makes the model relatively easy and quick to reproduce for someone writing his own code. Moreover, in order to take into account the potential jump that occurs in the MEAs, a special source term in the transport equation is proposed here. This term avoids interfacial issues, can be discretized, integrated, and approximated like any other term in the transport equation and can also be implemented in other numerical method such as the Finite Element Method. Finally the equations presented here can either be solved using a commercial package or using a homemade code. The model presented here is implemented within Matlab with a homemade library of functions that is able to solve problems involving transport equations with the Finite Volume Method.

The Section 2 of this paper presents the equations used for the model. Based on that approach, the electrical behavior of a PEMFC is simulated by considering classical stack geometry and physical parameters as described in Section 3. The model is also used to evaluate the system sensibility to design parameters. In Section 4, simulations are performed to ana-lyze the impact of MEA anomalies on the electrical characteristics. Phenomena such as cell interactions, non-equipotential bipolar plates, and internal loop of current are highlighted.

Formulation

Electrical charges are driven through fuel cells, thanks to complex electrochemical processes that involve many mass and charges transfers. The 3D modeling of these phenomena requires resolutions of many coupled partial differential equations. But here the objective is to model large fuel cells stacks with short time computation. As a consequence resolving coupled partial differential equations at the microscopic level in each MEA of the stack is not suitable. The choice was made to focus on the macroscopic behavior of the stack and to not distinguish the nature of the charge carriers (protons or electrons). An equivalent charge carrier is considered for the whole stack. Mass and heat transfers are not considered yet. This section describes how only one partial differential equation can be chosen in order to compute the 3D steady state macroscopic distribution of current and potential inside a stack. This way of modeling allows the study of the electrical behavior of large stacks with coarse meshes and very efficient time computation.

Transport Equation in Steady State Current Conduction

The electrical behavior of a stack is based on the equations of electromagnetism. Usually the steady state charge conservation equation and the local Ohm's law are used assuming that the electric field derives from a potential:

div j 0 ( 1 j r E ( 2
E Àgrad V (3
where j is the current density, E the electric field, V the potential, and r the electrical conductivity that can vary according the region of the stack. Combining these three equations leads to the well-known transport equation:

div r grad V 0 ( 4
In PEMFC single cell models, Eq. ( 4) is generally used to model the transport of some charges in some regions of a cell but it is not used for all the charge carriers in all regions. Here the objective is to consider only one charge carrier and to have only one equation for the whole stack. But if Eq. ( 4) is used alone, it cannot take into account energy sources inside a volume. With this kind of formulation, energy sources can only be introduced via the boundary conditions. It is a problem here because the energy comes from the fuel cell itself and not from the exterior. This energy source is related to the potential jump induced by the electrochemical phenomena occurring in the Membrane Electrode Assembly (MEA). In consequence, an additional term related to the electrochemical behavior of the MEA needs to be added, so Eq. ( 4) can be used as a stand alone equation for the whole stack.

Addition of an Electromotive Force Term in the Transport Equation

In the previous section, the electrical field was assumed to derive from a potential. But it is not always true and the electrical field can be written in a more general form:

E E m À grad V (5
where E m is the electromotive field. This term allows the modeling of a force that can drive charges through a medium.

In stacks this force arises because of the electrochemical phenomena occurring inside the MEA. In this paper, this term will be often referred as the electromotive force. The reason is that the electromotive field is the electromotive force per unit length. Using again Eqs. ( 1) and ( 2), the new charge transport Eq. ( 6) includes now an additional term. Thanks to this known source term, the electrical behavior of a stack can be easily modeled with a set of two distributions parameters: the conductivity and the electromotive field.

div r grad V À div r E m 0 ( 6
The electromotive field is null inside bipolar plates and end plates. Its value in the MEA is the open circuit voltage (OCV) divided by the thickness of the MEA in the linear case.

Addition of Materials With Non-Linear Behavior

For sake of simplicity, our approach was previously developed in the case of a linear behavior for the fuel cell. Nevertheless, it is well known that a PEMFC exhibits a non-linear electrochemical behavior. This response can be characterized by a MEA's polarization curve resulting from a model or from experimental measurement. In this paper, an empirical polarization curve of the form (7) was chosen:

DV j U mo À R so j À A ln j j eq 1 ! (7
where U mo is the OCV of one cell, R so the total resistance of the MEA and A and j eq are fitted parameters that represent the activation overpotential due to electrochemical reactions. The analytical function ( 7) is a source of information that models the behavior of MEAs. A Newton-Raphson algorithm was used to couple this equation with Eq. ( 6) in MEA regions.

The procedure is presented in the annex.

PEMFC Stack Model

In this section, equations presented previously are used to model a PEMFC stack. An overview of the model is given by describing the geometry and by giving an example of param-eter distribution in the situation of a healthy behavior. This can be done by assuming a uniform gas concentration distribution through the stack. This means that the oxygen and hydrogen gas concentration is constant along the gas channel. These operations conditions are achieved when the stochiometric coefficient is larger than the unity or the reactant utilization is very low.

Geometry, Mesh, and Physic

This model is built up of multiple single repeat units stacked on top of each other and sandwiched between two end plates along the z-axis. In Figure 1, each color represents a physical region with different parameters. In the examples shown in this paper, the stack consists of eight cells and the current collecting was considered to be performed on the edge of the PEMFC stack. It is important to remind that there is only one region for all the thickness of a MEA. The cathode, the anode and the membrane are not distinguished. All the physical phenomena that take place in these different layers are averaged over the MEA thickness and represented by the local polarization law. Finally, different kind of boundary conditions can be used at the terminal of the stack. In order to control the stack with voltage or with current or with a load, Dirichlet, Neumann or Robin conditions are respectively used.

In order to reduce the size of the mesh, surface elements were used for all the MEA regions. As a consequence, transverse currents in MEA are considered null. Currents can use transverse paths only in bipolar plates or in end plates. With the mesh shown in Figure 1 simulations are performed within a few seconds. Results of simulations are not sensitive to the mesh density along the z-axis and a coarser mesh can be used if simulations have to be performed on a stack with a large number of cells. The mesh density chosen in the x-and y-directions gives an error inferior to 2%.

PEMFC Stack Parameters

Table 1 gives the parameters values used for each region in the case of a healthy fuel cell stack. For sake of simplicity, a linear polarization law can be first considered to simulate the electrical behavior of the stack. For a MEA, we consider an OCV for single cell equal to 0.9 V and an effective MEA conductivity equal to 5 S m -1 (Table 2).

Using a linear behavior for MEAs is a good approximation if what is needed is the study of the mechanisms that explain the cells interactions. If accurate results are desired, a more complex polarization law for MEA can be used like the one given in Eq. (2.5) and whose parameters are given in Table 3. MEA electromotive force (U m0 ) 0.9 V MEA resistance per unit area (R s0 ) 21 0 -5 X m -2

The parameters given here correspond to a stack fed with pure oxygen and composed of graphite bipolar plate.

Electric Behavior of a Healthy PEMFC Stack

The next figures present the electrical behavior of healthy PEMFC stacks considering the standard parameters (Table 1). The Figure 2 presents the simulated stack voltage versus the current flowing through the stack for the both cases of linear and non-linear electrochemical behaviors. The Figures 3a-c presents the current density and cell voltage distribution through the PEMFC stack. Simulation is performed considering a current equal to 100 A flowing through the stack. The current goes from the lower potential to the higher potential and not the opposite as in passive materials. This is due to the electromotive force term added in the transport equation. As the gas phase distribution is not taken into account, homogeneous distributions are observed whatever the current collector position is. The distribution of the streamlines through a cross section (see Figure 3c) shows that the larger end plate conductivity allows the redistribution of the current. It is worth mentioning that the seal on the edge of the cells induce a curvature of the streamlines. For sake of simplicity, the simulated results are obtained for the case of a linear behavior. Note that similar results are recorded when considering a non-linear law.

This PEMFC stack model is able to investigate the effects of material properties and design parameter on the stack's operating conditions. For example if a material with a conductivity of 10 4 S m -1 is chosen for the end plate, Figure 4 shows what would be the effect of such choice on the electrical behavior of the stack. Similar results are obtained if the thickness of the end plates is decreased by two orders of magnitude. In this case, heterogeneous cell voltage and current density distributions are recorded. Due to cell-to-cell electrical coupling, each cell has a slightly different electrical behavior. The highest current densities and lowest cell voltages are shown close to the external current collector (Figures 4a,b). The cells in the middle of the stack demonstrate the smoothest behavior. As expected, the distribution of the streamlines (Figure 4c) shows that the lower end plate conductivity allows only partly the redistribution of the current. In consequence, the positioning of the current collector can impact the current density and potential distribution when considering lower end plate conductivity.

Electric Behavior of a Faulty PEMFC Stack

The simulations results described in this section are performed to predict the change in current density and potential through the stack taking into account anomalies into the PEMFC stack. The main capability of this model is the possi-

MEA electromotive force (U m0 )

1V MEA resistance per unit area (R s0 )

1.79 10 -5 X m 2 MEA non-linear parameter (A) 0.025 V MEA non-linear parameter (j eq ) 9 0Am -2 bility to induce inhomogeneous distribution of electrochemical parameters. This technique of modeling allows the computation of the effects of inhomogeneous distribution of conductivities and electromotive forces. Membrane dehydration, water flooding as well as materials degradation are assumed to generate a MEA resistance increase. Conversely, reactant crossover and starvation can induced a decay of the electromotive force. The Figure 5 then presents the position of the induced anomalous resistance and anomalous OCV. We consider an arbitrary MEA conductivity four times lower for the first case and a OCV close to zero in the second case on one corner of the cell number 6. The last case considers the transient distribution of OCV (Figure 5c) during fuel cell shutdown or start-up when ambient air may be introduced or removed in the anode gas channels, respectively. This section presents the electrical behavior of a faulty PEMFC stack by considering these two types of anomalous and their consequences on the potential, current and performances. First, the Figure 6 provides the voltage stack versus current when considering either the resistive anomalies (Figure 5a) or the anomalous electromotive force (Figure 5b). The global effect of these two kinds of anomaly is highlighted by comparison to the case of the healthy stack. As expect, the slope of voltage stack versus current curve increases when considering an anomalous resistance meanwhile the OCV of the stack shift toward lower value in the case of an anomalous electromotive force.

Anomalous Resistance

In order to illustrate the impact of an anomalous resistance on the electrical behavior (Figure 5a), Figure 7 shows the simulated cell voltage, current density and streamlines distribution through the PEMFC stack. The significant effect of the non-uniform resistance of the cell 6 is demonstrated. The both current density and cell voltage distributions (Figure 7a,b) are highly modified to force the current to pass through the healthy part of the cell 6. As it can be seen on Figure 7, this anomaly partially blocks the current flowing through the anomalous zone and induces low cell voltage. Conversely, this redistribution induces largest current densities close to the anomalous zone. Such behavior certainly causes hot spot in vicinity of the anomaly. The faulty cell 6 is severely affected but neighboring cells 5 and 7 are also impacted due to cell-to-cell electrical coupling. In consequence, healthy cells can suffer from inhomogeneous current due to a neighbor defective cell. Namely, fifth and seventh cells exhibit lower current densities and higher cell voltages close to the anomalous zone. Nevertheless, it is worth mentioning that the importance of this coupling between cells heavily depends on two parameters. When both conductivity and the thickness of the bipolar plate increase, the importance of the coupling decreases. In consequence the importance of cell interactions will heavily depend on the technology used for the PEMFC. The streamlines (Figure 7c) highlight another electrical phenomenon. This anomalous resistance partially blocks the current in the main direction and forces it to use transverse paths. These transverse currents create voltage drops inside a same bipolar plate which is not anymore equipotential. The figure shows that equipotential lines are bent out. These variations of potential on a same bipolar plate can reach about tens or hundreds millivolts depending on the current flowing through the stack and the bipolar plate material conductivity.

Finally, it is important to note that the neighbor's cells (cells 5 and 7) exhibit higher cell voltages close to the anomalous resistance (Figure 7a). These higher values have two explanations: the voltage drop created by the transverse currents in bipolar plates and a lower ohmic drop in AME regions with less current flowing through. As a consequence the healthy cells located close to anomalies can suffer from accelerating degradation processes such as carbon corrosion and platinum dissolution due to these induced higher cell voltages. This effect of anomalous resistance was already investigated in Refs. [12][13][14][15]. The results presented here are consistent with the ones shown in these papers.

In Figure 7a, the voltage of the sixth cell at the bottom is very low but still positive. Larger anomalous resistance can induce negative voltage, over a part of the cell or on the entire cell (Figure 8). This phenomenon occurs when the surface of the fault becomes too much important or if the conductivity of the anomalous zone becomes very low. In fact, a one dimensional model is enough to explain the phenomenon of negative voltage cell. The voltage of one cell is equal to the electromotive force of the cell minus the ohmic drop term. If resistance of the cell increases too much, the ohmic drop term becomes higher than the electromotive force term and the voltage of the cell becomes negative.

Anomalous Gas Phase Distribution

Figure 9 shows simulations performed with a local electromotive force (or OCV) assumed equal to zero at the bottom of the sixth cell (Figure 5b). The significant effect of a non-uniform gas distribution in the cell 6 is demonstrated. The both current density and cell voltage distributions are highly modified. The cell voltages exhibit lower values (Figure 9a) over the anomalous region meanwhile the current densities present negative values (Figure 9b). When few current is allowed to leave out the stack, the high electromotive force in the healthy region of cell 6 push the charge through area with lower electromotive force. The consequence here is this induced internal loop of current on the sixth cell (Figure 9c).

These simulations were performed using a MEA's linear response that allows currents in the opposite direction. So it has to be reminded that electrochemical non-linearity mainly hold back this phenomenon. Thus, the backward currents computed here are entirely over evaluated. But if the model is modified so that no reverse current is allowed, the results show that an electric field in the opposite direction remains in the anomalous region. This electric field produces a force that tries to move the charges backwards. If electrochemical processes allow some kind of charge transfer, a few current might flow in the opposite direction. The question about how much current will really flow backwards is beyond the scope of this paper. The point we want to emphasize here is that simple steady state electrical considerations can lead to configurations that create forces that try to push the charges in the opposite direction.

This behavior could be explained by the reverse current decay mechanism proposed by Reiser et al. [16] in the case of fuel starvation associated to improper gas supply. Because of the oxygen permeation through the membrane, oxygen reduction occurs at the anode meanwhile either oxygen evolution or carbon corrosion potentially occurs at the cathode. At the anode, electrons are then provided by the healthy part for oxygen reduction due to the transverse current whereas the healthy part of the cathode are partly fed by electrons coming from the faulty part of the cathode. Similar electrochemical behavior might be proposed in the case of oxidant starvation. This model predicts that the anodic regions filled with oxygen might create locally favorable reverse-current conditions in these exposed regions. The created reverse-current might cause local cathode degradation.

The model also shows that internal loops of current are amplified when considering a stack under reactive gases at the OCV (Figure 10). It also shows that if internal loops of currents exist (whatever are their origins), cells close to the anomaly have voltages higher than the OCV. The phenomenon is quite similar to the one encountered with anomalous resistance. As in this previous case, the transverse currents flowing in the bipolar plates create voltage drops that induce higher cell voltages on the neighboring cells. But here no current is flowing through the membrane of these healthy cells, so the local total voltage can become higher than the OCV. This means that these healthy cells can also suffer from accelerating carbon corrosion and platinum dissolution due to these higher cell voltages.

Unbalanced electromotive forces can also be encountered during shutdown or start-up. When a fuel cell is shutdown, ambient air may introduce in the anode gas channels where it replaces hydrogen. Conversely, during start-up, a transient fuel-air front is created along anode gas channels. In the both case, the anode catalyst layer is subjected to localized or even complete fuel starvation. When hydrogen pushes the air that was introduced to purge the stack (or respectively air pushes the hydrogen), there is a brief instant when a part of the MEA is subjected to an electromotive force while an other part is not. In consequence, this inhomogeneity creates an internal loop of current involving all the cells of the stack (Figure 11). Lastly it has to be reminded that other phenomena like the discharge of double layer capacity [17] can also contribute to reverse currents and mask this phenomenon.

Annex: Non-Linear electrochemical behavior of MEA

This annex describes how the differential Eq. (A.1) used to model the electrical behavior of a stack is coupled with a polarization law (Eq. A.2) that takes into account the non-linear response of the MEA. This approach is based on a Newton-Raphson algorithm that consists on linearizing an equation around some points until convergence is reached. The procedure given here can be used on any polarization law similar to Eq. (A.2) as far as derivation and convergence is possible. div r grad V À div r E m 0 ( A :1 DV j U mo À R so j À A ln j j eq 1 ! (A:2

Eq. (A.1) was obtained using the following general form of the local Ohm's law: j r E m À r grad V (A:3

The Ohm's law can also be written using global variables.

DV U m À R s j (A:4

If parameters are considered constants over the thickness t MEA of the MEA, the relations between local and global parameters can be written using Eq. (A.5).

r ≈ t MEA R s E mz ≈ U m t MEA

> < > :

(A:5 while Eq. (A.2) is non-linear, it is possible to state that it can be linearized around some point j n using a Taylor expansion of the first order.

DV j ≈ DV j n j À j n ∂DV ∂j j n (A:6

By identifying Eq. (A.4) with Eq. (A.6), the electromotive force and resistance per unit area around some point j n can be defined by Eq. (A. Using Eq. (A.2) and Eq. (A.7), it is possible to get analytical functions that give the variations of the resistance and electromotive force with the current density.

R s j n R so A 1 j n j ex U m j n U mo À A ln j n j ex 1 j n A 1 j n j ex 8 > > < > > :
(A:8

Eq. (A.1) is then solved iteratively using the parameters given by Eq. (A.8) and Eq. (A.5) at each iteration until convergence is reached.

Fig. 1

 1 Fig. 1 3D view of PEMFC stack's geometry and mesh with the localization of a cross section (green).

Fig. 2

 2 Fig. 2 Stack voltage versus current curves for the PEMFC stack: linear (᭺) and non-linear (*) electrochemical behaviors.

Fig. 3

 3 Fig. 3 Cell voltage in V (a), current density in A cm -2 (b) and streamlines on a cross section (c) for a 100A current flowing through the stack.

Fig. 4 Fig. 5

 45 Fig. 4 Cell voltage in V (a), current density in A cm -2 (b) and streamlines on a cross section (c) for a 100A current flowing through the stack: case of lower end plate conductivity (10 4 S m -1 ).

Fig. 6

 6 Fig. 6 Voltage stack versus current curve obtained in case of a healthy behavior (o), in the case of an anomalous resistance (+) and in the case of an anomalous OCV (ᮀ).

Fig. 7

 7 Fig. 7 Cell voltage in V (a), current density in A cm -2 (b) and streamlines on a cross section (c) for a 100A current flowing through the stack: case of anomalous resistance.

Fig. 8

 8 Fig. 8 Cell voltage distribution for a 150 A current flowing through the stack with a larger anomalous resistance.

Fig. 9

 9 Fig. 9 Cell voltage in V (a) current density in A cm -2 (b) and streamlines on a cross section (c) for a 100 A current flowing through the stack: case of anomalous OCV.

Fig. 10

 10 Fig. 10 Cell voltages (a) and streamlines on a cross section (b) for a 0 A current flowing through the stack in the case of an anomalous OCV.

Fig. 11

 11 Fig. 11 Streamlines on a cross section during start-up or shutdown operations (0 A of total current).

Table 1

 1 Stack parameters.

	Current collector conductivity	5 10 7 S m -1
	End plate conductivity	5 10 7 S m -1
	Bipolar plate conductivity	5 10 3 S m -1
	Current collector thickness	1 mm
	End plate thickness	2.5 mm
	Bipolar plate thickness	4 mm
	MEA thickness	0.1 mm
	Current collector area	0.02 × 0.02 m 2
	End plate area	0.14 × 0.14 m 2
	Bipolar plate surface area	0.14 × 0.14 m 2
	MEA surface area	0.1 × 0.1 m 2

Table 2

 2

Parameters used for a MEA for linear polarization law (healthy PEMFC stack).

Table 3

 3 Parameters used for a MEA for non-linear polarization law (healthy PEMFC stack).

Acknowledgment

The authors thank the Agence Nationale de la Recherche (ANR) for their financial support as well as the partners of the OMNISCIENTS project: PSA Peugeot Citroën, Helion, Commissariat à l'Energie Atomique et aux Energies Alternatives (CEA) and Euro Physical Acoustics (EPA).