František Němec
email: frantisek.nemec@gmail.com

O Santolík

M Parrot

Doppler Shifted Alpha Transmitter Signals in the Conjugate Hemisphere: DEMETER Spacecraft Observations and Raytracing Modeling

DEMETER measurements conjugate to Alpha transmitters show three distinct 10 ways of signal propagation.

11

• Two of them are identified as ducted and unducted propagation, the third is likely 12 only partially ducted.

13

• A raytracing analysis is employed to explain the observed Doppler shifts.

Introduction

Powerful ground-based very low frequency (VLF) transmitters operating at frequencies between about 10 and 30 kHz are used by military for long distance communications with submerged submarines and for navigation. They can propagate considerable distances in the Earth-ionosphere waveguide [START_REF] Budden | The wave-guide mode theory of wave propagation[END_REF]. Part of the signal power may eventually escape the waveguide and propagate to space, being somewhat attenuated when penetrating through the ionosphere, in particular during the daytime (Cohen et al., 2012;[START_REF] Graf | Analysis of experimentally validated trans-ionospheric attenuation estimates of VLF signals[END_REF]. Having made it to the space environment, the transmitter signals are detectable by low-altitude spacecraft both in the source hemisphere and in the magnetically conjugated region (Cohen & Inan, 2012;[START_REF] Zhang | The spatial distribution of three Alpha transmitter signals at the topside ionosphere[END_REF][START_REF] Zhao | Investigation of VLF transmitter signals in the ionosphere by ZH-1 observations and full-wave simulation[END_REF].

The signals are observable also at larger radial distances in the proximity of given magnetic field lines [START_REF] Foster | Observations of the impenetrable barrier, the plasmapause, and the VLF bubble during the 17 March 2015 storm[END_REF][START_REF] Ma | VLF waves from ground-based transmitters observed by the Van Allen Probes: Statistical model and effects on plasmaspheric electrons[END_REF][START_REF] Koronczay | VLF transmitters as tools for monitoring the plasmasphere[END_REF].

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Space Physics VLF transmitter signals may result in ionospheric heating in the vicinity of transmitters [START_REF] Bell | DEME-TER observations of ionospheric heating by powerful VLF transmitters[END_REF][START_REF] Němec | NWC transmitter effects on the nightside upper ionosphere observed by a low-altitude spacecraft[END_REF], as well as in formation of specific wave events due to wave-wave coupling and triggering [START_REF] Němec | Line radiation events induced by very low frequency transmitters observed by the DEMETER spacecraft[END_REF][START_REF] Gu | Alpha transmitter signal reflection and triggered emissions[END_REF].

They can also interact with energetic electrons in Van Allen radiation belts and eventually result in their precipitation (e.g., [START_REF] Karpman | Particle precipitation caused by a single whistler-mode wave injected into the magnetosphere[END_REF][START_REF] Abel | Electron scattering loss in the Earth's inner mag-netosphere: 1. dominant physical processes[END_REF][START_REF] Selesnick | Influence of a ground-based VLF radio transmitter on the inner electron radiation belt[END_REF], as was demonstrated by the analysis of precipitating electron fluxes measured at low altitudes [START_REF] Gamble | Radiation belt electron precipitation by man-made VLF transmissions[END_REF][START_REF] Sauvaud | Radiation belt electron precipitation due to VLF transmitters: Satellite observations[END_REF]. The efficiency of the transmitter induced pitch angle scattering is generally different for ducted and unducted propagation [START_REF] Rodger | Contrasting the efficiency of radiation belt losses caused by ducted and nonducted whistler-mode waves from ground-based transmitters[END_REF][START_REF] Ross | Effects of VLF transmitter waves on the inner belt and slot region[END_REF]. Nevertheless, it remains unclear how much energy actually propagates ducted or unducted. At least for low L-shell transmitters, a significant portion of the wave energy may propagate unducted [START_REF] Clilverd | Ground-based transmitter signals observed from space: Ducted or nonducted?[END_REF][START_REF] Zhang | Observed propagation route of VLF transmitter signals in the magnetosphere[END_REF]. However, at higher L-shells, the waves are possibly primarily ducted, as suggested by the lack of significant wave power above the electron half gyrofrequency for inter-hemispherically propagating waves [START_REF] Clilverd | Ground-based transmitter signals observed from space: Ducted or nonducted?[END_REF].

Ducted propagation should correspond primarily to low wave normal angles, while large wave normals are expected for unducted propagation in the magnetically conjugated hemisphere [START_REF] Bortnik | Temporal signatures of radiation belt electron precipitation induced by lightning-generated MR whistler waves: 1[END_REF]. The two propagation types could be thus in principle distinguished by evaluating the wave normal angles. However, this requires multicomponent magnetic field measurements, which are typically not available with a sufficient frequency resolution at frequencies well above 10 kHz. When the wave normal angles are high and close to the resonance cone, the refractive index is large and the signal is detected onboard a low altitude spacecraft with a noticeable Doppler shift [START_REF] Starks | Modeling of Doppler-shifted terrestrial VLF transmitter signals observed by DEMETER[END_REF]. This allows to distinguish between the ducted and unducted propagation paths, and eventually to estimate the wave normal angle variation with the spacecraft location.

We use electromagnetic wave data obtained by a low-altitude Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions (DEMETER) spacecraft to analyze Doppler shifts of Alpha navigation transmitter signals detected in the conjugate hemisphere. We distinguish between the ducted and unducted propagation, and we speculate about the respective propagation paths. The used data set is presented in section 2 and the observational results are presented in section 3. Wave trajectories calculated using a raytracing routine and their consistency with the observed Doppler shifts are eval-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Space Physics uated in section 4. The obtained results and the overall picture are discussed in section 5.

Section 6 contains a brief summary of the main findings.

Data Set

DEMETER was a French low-altitude spacecraft which operated between 2004 and 2010. It had a circular orbit with an altitude of about 660 km, which was nearly Sunsynchronous. The measurements were thus performed at local times either close to 10:30 hours ("daytime", north-to-south orbits) or 22:30 hours ("nighttime", south-to-north orbits).

Survey mode electromagnetic wave measurements were active nearly continuously at geomagnetic latitudes lower than about 65 degrees. In the VLF range (up to 20 kHz), these consisted of onboard calculated power spectra of one electric [START_REF] Berthelier | ICE, the electric field experiment on DEMETER[END_REF] and one magnetic [START_REF] Parrot | The magnetic field experiment IMSC and its data pro-Planet[END_REF] ,equal to about 11,905 Hz,12,649 Hz,and 14,881 Hz. A 0.4 s long pulse at a given frequency is always followed by a 0.8 s long gap, and then a pulse of a different frequency is radiated. The entire cycle thus takes 3.6 s, and then the sequence is repeated again (Cohen & Inan, 2012). south. The signals from Alpha navigation transmitters can be identified at frequencies of about 11.9, 12.6, and 14.9 kHz, both in the northern and southern hemispheres. Additionally, the signals in the southern hemisphere apparently split to form a V-like shape. The white bars at the top indicate the time intervals when the Burst mode was active and high resolution data were available. The white arrow indicates the magnetic latitude of the Khabarovsk transmitter.

Observations

north to south. Alpha navigation transmitter signals at frequencies of about 11.9, 12.6, and 14.9 kHz are identifiable. During the first half of the plotted time interval, the spacecraft was in the northern hemisphere, i.e., in the hemisphere of the transmitters. The signals are limited to the individual transmitted frequencies and they peak at the time when the spacecraft passes close to the transmitter. The white arrow in Figure 1 indicates the magnetic latitude of the Khabarovsk transmitter which is likely the source of the observed signal, as we show below.

The situation in the southern hemisphere, i.e., in the hemisphere geomagnetically conjugated, is more complicated. The signals at individual transmitted frequencies are still observable, albeit with a somewhat lower intensity. Additionally, an apparent gradual decrease of the signal frequencies followed by their subsequent relaxation back to the transmitted values, forming a V-like shape, is observed toward the end of the plotted time interval. It is noteworthy that the V-shapes at higher frequencies are observed slightly sooner (at lower geomagnetic latitudes) than the V-shapes at lower frequencies. The highest frequency V-shape is observed about 13 s sooner than the middle frequency V-shape, which is followed by the lowest frequency V-shape after about 6 s. The corresponding differences in geomagnetic latitude are lower than one degree. A more detailed view of the situation obtained using the Burst mode data is shown below in Figure 3. Note also that the intense vertical lines observable in the plot correspond to lightning generated sferics/whistlers, and they are not related to the topic of the present paper. There were two Burst mode intervals during the analyzed half-orbit, whose exact times are marked by the white bars at the top of the plot.

Figure 2a shows a geographic map of the situation during the time interval plotted in Figure 1. The projection of the satellite path on the ground is shown by the red -5-

Accepted Article

This article is protected by copyright. All rights reserved. curve. The red bars mark the locations where the Burst mode was active. The direction of the orbit (north to south) is marked by the blue arrow. The black crosses in the northern hemisphere show the locations of Alpha navigation transmitters (Novosibirsk and Khabarovsk). The black crosses in the southern hemisphere show their geomagnetic conjugates. Given the satellite path, it seems reasonable to assume that the observed signals come from the Khabarovsk transmitter. The geomagnetic longitudinal separation between the satellite and the Khabarovsk transmitter is lower than about 10 degrees during the entire half-orbit.

Accepted Article

This article is protected by copyright. All rights reserved. was located in the northern hemisphere. The discrete transmission pattern of the Alpha navigation transmitters can be identified at the three frequencies (about 11.9, 12.6, and 14.9 kHz). (b)

The same as (a), but for the second Burst mode time interval, when the spacecraft was located in the southern hemisphere. Transmitter signals with gradually decreasing and later increasing frequencies are observed in addition to the expected constant-frequency signals.

ing of 0.4 s long pulses alternating in frequency can be seen. The observed signal frequencies correspond well to the frequencies radiated by the transmitter. The occasional vertical lines over the plotted frequency range are lightning generated sferics/whistlers, not of importance for the present analysis.

Figure 3b shows a frequency-time spectrogram of power spectral density of electric field fluctuations corresponding to the second Burst mode time interval. This time the spacecraft was in the southern hemisphere close to the region geomagnetically conjugate to the transmitter location. Discrete transmission pattern at the three frequencies can be again identified. Additionally, there are bands of comparatively lower frequencies associated with each of the radiated frequencies. These appear first to decrease in frequency. After reaching a frequency minimum, they gradually increase in frequency again, forming V-like shapes. The frequency minima at higher frequencies occur at earlier times (lower geomagnetic latitudes). Just before reaching the frequency minima, each radiated transmitter frequency f 0 is observed to be apparently split into three distinct frequencies f obs .

The signals at frequencies lower than the transmitted signal frequency can be interpreted in terms of a Doppler shift of the original transmitted signal frequency:

f obs = f 0 1 - n • v c (1)
where n is a refractive index vector (directed along the wave vector), v is the DEME-TER velocity vector, and c is the speed of light. Considering that v ≈ 7.5 km/s and the maximum observed Doppler shifts are as large as about 2%, the corresponding re-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Space Physics fractive indices have to be very large, close to 1000. This is possible for large wave normal angles, close to the resonance cone.

The apparent splitting of each radiated frequency into three distinct frequency bands can be then interpreted as the radiated signal propagating to the spacecraft over three different propagation paths. Each of these eventually results in a different wave normal angle at the observation point, a different refractive index vector (both in direction and magnitude), and thus also a different Doppler shift. The propagation path resulting in principally zero Doppler shift can be interpreted as ducted propagation of the signal, with wave vectors oriented roughly along the ambient magnetic field. The values of refractive index corresponding to such a propagation are quite small and the angle between the wave vector and the satellite velocity vector is rather large. The expected Doppler shift is thus very small [START_REF] Starks | Modeling of Doppler-shifted terrestrial VLF transmitter signals observed by DEMETER[END_REF].

Identifying the two remaining propagation paths is less straightforward. There is a propagation path with the Doppler shift increasing toward larger geomagnetic latitudes.

This path ceases to exist at very large geomagnetic latitudes and it can be attributed to an unducted propagation. The unducted signals reach the conjugate hemisphere with quite oblique wave normals and exhibit a frequency-dependent latitudinal threshold beyond which they cannot reach [START_REF] Starks | Modeling of Doppler-shifted terrestrial VLF transmitter signals observed by DEMETER[END_REF]. The threshold latitude decreases with an increasing wave frequency [START_REF] Shklyar | On the origin of lower-an upper-frequency cutoffs on wedge-like spectrograms observed by DEMETER in the midlatitude ionosphere[END_REF]. This can explain why at higher frequencies the bottoms of the V-like shapes are observed closer to the magnetic equator.

The identification of the remaining propagation path is the most challenging. It allows the signals to propagate to large latitudes in the opposite hemisphere at not-too-large wave normals (small Doppler shifts). We suggest that it corresponds to only partially ducted propagation. This is supplemented by a raytracing analysis in section 4.

Nighttime (South-to-North) Orbits

Before progressing to the raytracing analysis, it is instructive to provide an additional experimental supporting argument for the Doppler shift explanation of the observed signals. During nighttime half-orbits the spacecraft moved in an opposite direction as compared to daytime half-orbits. Assuming that the propagation of the signals remains at least qualitatively similar, one might thus expect the Doppler shifts observed during the night to be opposite than during the day. This is indeed the case, as demonstrated

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Space Physics whistlers, the signals from Alpha navigation transmitters can be identified at frequencies of about 11.9, 12.6, and 14.9 kHz, both in the northern and southern hemispheres. Additionally, the signal at the frequency of about 14.9 kHz observed in the southern hemisphere apparently splits to form an inverted V-like shape. The white bars at the top indicate the time intervals when the Burst mode was active and high resolution data were available. i.e., in the hemisphere of the transmitters. In the first half of the plotted time interval, the observed frequencies closely match the transmitted frequencies only for the two lower frequency signals (11.9 and 12.6 kHz). A considerable increase of the signal frequency is apparent for the highest frequency (14.9 kHz) signal. The absence of the Doppler shifted, i.e., (partially) unducted, signals at lower frequencies can be explained by a magnetospheric reflection. The nighttime lower-hybrid frequency above the spacecraft is higher than the signal frequencies, resulting in their reflection above the satellite altitude [START_REF] Shklyar | Characteristic properties of Nu whistlers as inferred from observations and numerical modelling[END_REF][START_REF] Vavilov | Study of the lower hybrid resonance frequency over the regions of gathering earthquakes using DEMETER data[END_REF].

Accepted Article

This article is protected by copyright. All rights reserved. Note that another bright dot of statistically increased signal levels occurs above the NWC transmitter located at North West Cape, Australia (21.8 • S, 114.2 • E) and operating at 19.8 kHz, owing to spectral broadening [START_REF] Xia | Spectral broadening of NWC transmitter signals in the ionosphere[END_REF].

The example half-orbit contains two Burst mode intervals for which high-resolution frequency-time spectrograms can be obtained. These are depicted in Figure 6. Figure 6a corresponds to earlier times when the spacecraft was located in the southern hemisphere (i.e., conjugate to the transmitter). It shows the time coded signals at frequencies of about 11.9 and 12.6 kHz, which directly correspond to the transmitted signal. The signal at a frequency of about 14.9 kHz is partly obscured by strong natural (lightning related) background emissions. However, the signal apparently splits and extends to higher frequencies around the middle of the plotted time interval, forming an inverted V-like shape.

This is particularly clear in the left part of the plot, where the background emissions are

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Space Physics weaker and the transmitter signal is more pronounced. An indication of this transmitter signal shift to higher frequencies is identifiable also after 14:06:30 UT, corresponding to the second half of the inverted V-like shape. Note that the original signal at the frequency of 14.9 kHz does not disappear during the time of the inverted V-like shape.

It still remains identifiable, albeit being arguably weaker. Note also that no such deviation of the signal frequencies from the transmitted frequencies is observable for signals at 11.9 and 12.6 kHz, in agreement with the Survey mode plot from Figure 4. Figure 6b, shown for completeness, depicts the situation in the transmitter hemisphere. The signals at frequencies of about 11.9, 12.6, and 14.9 kHz are identifiable despite the strong natural background related to the lightning activity. They correspond directly to the frequencies transmitted by the Alpha navigation transmitters.

Raytracing Analysis

Having suggested the Doppler shift as an explanation of the frequency-time properties of the observed transmitter signals, we devote this section to a more quantitative analysis. Ducted signals, which propagate along the ambient magnetic field with not too large wave normals, have too low refractive indices and do not result in observable Doppler shifts. They can be thus identified as the signals in the conjugate hemisphere at frequencies corresponding to the transmitted frequencies. As such, their interpretation is straightforward and they are not further discussed in this section. Instead, we focus on the interpretation of the signals at frequencies apparently departing from the transmitted frequencies, which we interpret as unducted or only partially ducted signals.

We use a raytracing analysis which allows us to follow individual ray trajectories from the northern hemisphere to the conjugated hemisphere and to evaluate the respective wave normal angles and refractive indices upon reaching the DEMETER altitudes.

The used raytracing code is the same as the one used formerly by, e.g, [START_REF] Santolík | Propagation of equatorial noise to low altitudes: Decoupling from the magnetosonic mode[END_REF] and [START_REF] Martinez-Calderon | Spatial extent of quasiperiodic emissions simultaneously observed by Arase and Van Allen Probes on 29 november 2018[END_REF]. A simple dipole geomagnetic field model is used.

The density is modeled assuming a diffusive equilibrium with a reference altitude of 670 km, ion and electron temperatures 1000 K, and plasma number density at the reference altitude 18,000 cm -3 with 92% of oxygen ions and 8% of hydrogen ions. The propagation takes place well within the plasmasphere, so that no plasmapause density drop is considered in the density model. Note, however, that the exact choice of the density model does not significantly affect the raytracing results. The important point is that the max-

-11-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Space Physics imum lower hybrid frequencies along the propagation paths are lower than all the considered signal frequencies (11,905 Hz, 12649 Hz, and 14,881 Hz), corresponding to a typical daytime situation. This is not the case during the night, when the maximum lower hybrid frequencies along the propagation paths are typically higher than the lower two signal frequencies. They thus get magnetospherically reflected and do not reach the DEME-TER altitudes in the conjugate hemisphere [START_REF] Chum | Analysis of subprotonospheric whistlers observed by DEMETER: A case study[END_REF][START_REF] Shklyar | On the origin of lower-an upper-frequency cutoffs on wedge-like spectrograms observed by DEMETER in the midlatitude ionosphere[END_REF][START_REF] Vavilov | Study of the lower hybrid resonance frequency over the regions of gathering earthquakes using DEMETER data[END_REF]. This explains why during the nighttime half-orbits the (partially) unducted signals are observed only at the highest frequency .

Let us start the discussion with purely unducted signals, whose Doppler shifts observed by DEMETER were formerly analyzed by [START_REF] Starks | Modeling of Doppler-shifted terrestrial VLF transmitter signals observed by DEMETER[END_REF]. We assume that, before exiting to higher altitudes, the transmitter generated signals propagate some distance within the Earth-ionosphere waveguide. Referring to the maps of mean wave intensities in Figures 2b-d, the signals are considered to exit the waveguide within 8 degrees of the geomagnetic latitude of the transmitter (about 41.4 • ). Note that we limit the calculation to the magnetic meridian of the transmitter, but principally the same picture holds also for nearby geomagnetic longitudes. Considering that the refractive index within the Earth-ionosphere waveguide is close to 1, a wave exiting to the ionosphere bends according to the Snell's law toward the normal. It is thus reasonable to start the raytracing assuming vertical wave normals at altitudes of about 110 km [START_REF] Helliwell | Whistlers and related ionospheric phenomena[END_REF].

The calculated ray trajectories for three different starting geomagnetic latitudes

Accepted Article

This article is protected by copyright. All rights reserved. and 14,881 Hz, respectively. It can be seen that the waves starting at larger geomagnetic latitudes end at larger geomagnetic latitudes (i.e., the calculated ray trajectories do not intersect anywhere). However, there is a clear upper geomagnetic latitude limit for the ending points of the propagation paths [START_REF] Starks | Modeling of Doppler-shifted terrestrial VLF transmitter signals observed by DEMETER[END_REF][START_REF] Shklyar | On the origin of lower-an upper-frequency cutoffs on wedge-like spectrograms observed by DEMETER in the midlatitude ionosphere[END_REF]. Additionally, higher frequency waves end at lower geomagnetic latitudes than lower frequency waves started at the same geomagnetic latitude.

In order to understand the observed Doppler shifts, it is crucial to consider the refractive index value and wave vector direction upon reaching the DEMETER altitudes in the conjugate hemisphere. It is found that the wave vectors are typically rather close to the resonance cone, being oriented toward larger L-shells, i.e., approximately southward. This corresponds well to the observed Doppler shifts. During the daytime, when the spacecraft moves from north to south, the satellite velocity is oriented roughly along the wave vectors. The Doppler shift then results in observed frequencies being lower than the transmitted frequencies. During the nighttime, the spacecraft moves from south to north. The satellite velocity is thus oriented roughly opposite to the wave vectors, and the Doppler shift results in observed frequencies being larger than the transmitted frequencies. As the wave vectors are close to the resonance cone, the values of the refractive index are high and the Doppler shift significant.

Wave normal angles and refractive index values as a function of the geomagnetic latitude in the conjugate hemisphere are analyzed more in detail in Figures 8a and8b, respectively. These were obtained by tracing the waves started within a full range of considered geomagnetic latitudes in the northern hemisphere until they reach the DEME-TER altitudes in the southern hemisphere. The blue, green, and red color curves again correspond to individual transmitter signal frequencies. The horizontal dashed lines in Figure 8a show the respective resonance cone angles. Given that the waves propagate

-13-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Space Physics from north to south, the wave normal angles are larger than 90 degrees. However, they gradually decrease and approach the resonance cone as the final geomagnetic latitudes approach a given threshold value. Correspondingly, the refractive index (and thus also the Doppler shift) gradually increases and eventually diverges at the threshold final geomagnetic latitude. This threshold latitude is different for each wave frequency, being slightly closer to the geomagnetic equator for higher frequency waves. The unducted propagation can thus explain the left-hand parts of V-like shapes in Figure 1 and the righthand part of the inverted V-like shape in Figure 4.

The remaining problem is how to explain the right-hand parts of V-like shapes in Figure 1 and the left-hand part of the inverted V-like shape in Figure 4. These are formed by the signals which make it to larger geomagnetic latitudes than the unducted propagation threshold. Moreover, their refractive index decreases with increasing geomagnetic latitude. As both the purely ducted and purely unducted signals have been already discussed and identified in the frequency-time spectrograms, we need to look for a possible third way of the signal propagation. We suggest that the signal might be ducted only for a part of the propagation path. Significant Doppler shifts are observed upon reaching the conjugate hemisphere, i.e., the final wave normal angles are close to the resonance cone. The signal is thus clearly unducted when reaching the spacecraft. One could consider many different ducted/unducted propagation combination schemes. However, in order to keep things reasonably simple, we assume the signal to be ducted all the way to the geomagnetic equator, and unducted farther. The upper frequency limit on the wave ducting in the ducts of an enhanced density is 0.5f ce cos θ k [START_REF] Smith | Propagation characteristics of whistlers trapped in field-aligned columns of enhanced ionization[END_REF], where θ k is the wave normal angle and f ce the electron cyclotron frequency. Moreover, even the waves

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Space Physics Figure 10. The same as Figure 8, but for the propagation scheme from Figure 9.

with frequencies below the frequency limit or those ducted in density troughs are less likely to be ducted when the limit is lower. It is thus reasonable to assume that at least more oblique waves can become unducted close to the magnetic field minimum at the geomagnetic equator. We further assume a rather large range of wave normal angles leaving the duct [START_REF] Woodroffe | Whistler propagation in ionospheric density ducts: Simulations and DEMETER observations[END_REF]

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Space Physics

The time delays corresponding to the interhemispheric ducted propagation are on the order of a few tenths of second. When the wave normal angle upon arrival is not too close to the resonance cone, the respective time delays are very close to those obtained for the ducted propagation. This is the case for the unducted propagation at lower geomagnetic latitudes and for the partially ducted propagation at larger geomagnetic latitudes. However, the time delays increase considerably when the signals propagate close to the resonance cone, and they diverge when the final wave normal angle approaches the resonance. This explains the considerable time delays of individual transmitter pulses observed for the largest Doppler shifts in Figures 3b and6a.

Discussion

Electromagnetic wave measurements performed by the low-altitude DEMETER spacecraft allow us to observe signals from powerful ground-based VLF transmitters. Unfortunately, only a single electric and a single magnetic field components are measured, which does not allow us to perform a detailed wave analysis and to directly experimentally determine the wave normal angles. Moreover, as the magnetic field component of the analyzed signals is typically too weak to be detected by the used instrumentation, the observations are limited to a single electric field component.

Doppler shifts due to the spacecraft movement are typically rather small and difficult to identify, unless the refractive index becomes comparatively large. This is, however, the case of whistler waves propagating with wave normal angles close to the resonance cone. Then, the Doppler shifts are readily observable and their magnitude can be used to estimate the refractive index.

The observations performed in the hemisphere conjugate to the Alpha navigation transmitters revealed that during specific times/at some spacecraft locations, a given trans-

Accepted Article

This article is protected by copyright. All rights reserved. A rather surprising finding is the third propagation path necessary to explain the third observed signal frequency. We suggested that the corresponding signal might be ducted in the beginning and become eventually unducted close to the equatorial region.

Under such a propagation scheme, the raytracing results and the predicted Doppler shifts are in agreement with the observations. Most importantly, the Doppler shift steadily decreasing with increasing geomagnetic latitude is obtained. The wave has to be assumed to be initially ducted along a magnetic field line slightly northward from the transmitter location (i.e., at a larger L-shell). There may be two possible reasons why this should be the case. First, the ionospheric attenuation is lower at larger geomagnetic latitudes [START_REF] Helliwell | Whistlers and related ionospheric phenomena[END_REF], so that one might expect more VLF transmitter power to penetrate through in there. Second, and probably more important, the wave ducting depends crucially on the electron cyclotron frequency. It becomes less likely for low electron cyclotron frequencies [START_REF] Smith | Propagation characteristics of whistlers trapped in field-aligned columns of enhanced ionization[END_REF], suggesting that the waves tend to get unducted at larger Lshells close to the geomagnetic equator. This propagation scheme requires the waves to leave the duct with large range of wave normal angles. This may arguably correspond to the real situation. However, the large range of initial wave normal angles considered at a given location can also account for the waves leaving the duct over a range of latitudes/Lshells, which must be the case in reality.

Achieving an exact quantitative match between the calculated and observed Doppler shifts is difficult, as the wave propagation at large wave normal angles is very sensitive to even small variations of the conditions. However, the obtained qualitative agreement suggests that our interpretation of the three propagation paths is indeed correct. It is further supported by the fact that as the spacecraft moves in different directions during the daytime and nighttime half-orbits, the signs of the Doppler shifts are expected to be different, corresponding to the observations. Finally, it is worth to note again that in order to reach the DEMETER altitudes in the conjugate hemisphere the signal frequencies have to be higher than the maximum lower hybrid frequencies along the respective propagation paths. This immediately explains why the Doppler shifted signals are not observed at lower frequencies during the night, when the lower hybrid frequencies are higher [START_REF] Chum | Analysis of subprotonospheric whistlers observed by DEMETER: A case study[END_REF][START_REF] Shklyar | On the origin of lower-an upper-frequency cutoffs on wedge-like spectrograms observed by DEMETER in the midlatitude ionosphere[END_REF][START_REF] Vavilov | Study of the lower hybrid resonance frequency over the regions of gathering earthquakes using DEMETER data[END_REF].

  field component. The frequency resolution of the data was about 20 Hz, and the time resolution was about 2 s. Additionally, higher resolution data were available during sporadically active Burst mode intervals. These consisted of waveforms of one electric and one magnetic field component sampled at 40 kHz. However, some signals analyzed in the present study are mainly electrostatic and their magnetic field intensity is very low. They are thus generally undetectable in the magnetic field measurements, and only the electric field data are used. The Alpha navigation system consists of three transmitters: Krasnodar (45.4 • N, 38.2 • E), Novosibirsk (55.8 • N, 84.4 • E) and Khabarovsk (50.1 • N, 136.6 • E) which transmit specific discrete periodic patterns. Three different frequencies are used by each of the transmitters

3. 1

 1 Figure 1. Frequency-time spectrogram of power spectral density of electric field fluctuations measured on 1 March 2010 during a daytime half-orbit, when the spacecraft moved from north to

  Figure 2. (a) Geographic map showing the satellite orbit for the time interval from Figure 1 (red curve) and the locations where the Burst mode was active (red bars). The direction of the orbit is marked by the blue arrow. The black crosses in the northern hemisphere correspond to the Alpha transmitter locations, while the black crosses in the southern hemisphere correspond to their geomagnetic conjugates. (b) Average daytime power spectral density of electric field fluctuations at the frequency of 11,905 Hz is color coded as a function of the geographic location according to the color scale on the right-hand side. (c) The same as (b), but for the frequency of 12,649 Hz. (d) The same as (b), but for the frequency of 14,881 Hz.

Figures

  Figures 2b-d show geographical maps of average daytime power spectral densities of electric field fluctuations at respective transmitter frequencies. The entire DEMETER data set was used when constructing the plot. The central frequencies considered are marked in each of the figures. The frequency bandwidth of about 60 Hz was used (3 frequency bins of the instrument). The crosses mark again the locations of the Alpha navigation transmitters and their geomagnetic conjugates. The average wave intensities at all the three analyzed frequencies are significantly increased both above the transmitter locations and close to their geomagnetic conjugates. It is apparent that the waves detected during the example orbit from Figures 1 and 2a indeed come from the Khabarovsk transmitter.

  Figure 3. (a) Frequency-time spectrogram of power spectral density of electric field fluctuations corresponding to the first Burst mode time interval from Figure 1, when the spacecraft

Figure 4 .

 4 Figure 4. Frequency-time spectrogram of power spectral density of electric field fluctuations measured on 28 February 2010 during a nighttime half-orbit, when the spacecraft moved from south to north. Although background intensities are quite high due to lightning generated

Figure 5 .

 5 Figure 5. (a) Geographic map showing the satellite orbit for the time interval from Figure 4 (red curve) and the locations where the Burst mode was active (red bars). The direction of the orbit is marked by the blue arrow. The black crosses in the northern hemisphere correspond to the Alpha transmitter locations, while the black crosses in the southern hemisphere correspond to their geomagnetic conjugates. (b) Average power spectral density of electric field fluctuations at the frequency of 11,905 Hz is color coded as a function of the geographic location according to the color scale on the right-hand side. (c) The same as (b), but for the frequency of 12,649 Hz. (d) The same as (b), but for the frequency of 14,881 Hz.

  Figure 6. (a) Frequency-time spectrogram of power spectral density of electric field fluctuations corresponding to the first Burst mode time interval from Figure 4, when the spacecraftwas located in the southern hemisphere. The discrete transmission pattern of the Alpha navigation transmitters can be identified at the three frequencies (about 11.9, 12.6, and 14.9 kHz).Additionally, a gradually increasing signal frequency is observed in the upper left part of the figure. (b) The same as (a), but for the second Burst mode time interval, when the spacecraft was located in the northern hemisphere.

Figure 5

 5 Figure5uses the same format as Figure2to depict the situation during the example nighttime half-orbit and nighttime intensity maps. During nighttime half-orbits the spacecraft moved in a northward direction, as marked by the blue arrow in Figure5a. The longitudes of the example orbit are between the longitudes of the Novosibirsk and Khabarovsk transmitters marked by the black crosses. The transmitter locations and the average intensity maps in Figures5b-dindicate that in the beginning of the example time interval (i.e., in the southern hemisphere) the spacecraft likely detects predominantly the signal from the Khabarovsk transmitter. On the other hand, in the end of the example time interval (i.e., in the northern hemisphere) the spacecraft likely detects predominantly the signal from the Novosibirsk transmitter. Nevertheless, all the Alpha transmitters radiate at the same frequencies and comparable latitudes, only with a different time coding. The exact transmitter detected at a given time is thus not important for this study.

Figure 7 .

 7 Figure 7. Example raytracing trajectories corresponding to an unducted signal propagation. The signals were started at geomagnetic latitudes 33.4 • , 41.4 • , and 49.4 • , roughly corresponding to the geomagnetic latitude of the Khabarovsk transmitter and geomagnetic latitudes of ±8 • from it. The waves were started directly upwards (perpendicular to the ground surface) at an altitude of 110 km, as indicated by the short black arrows. The blue, green, and red curves were obtained for wave frequencies of 11,905 Hz, 12,649 Hz, and 14,881 Hz, respectively. The black dashed curves show the respective dipolar magnetic field lines corresponding to individual starting points.

(

  33.4 • , 41.4 • , and 49.4 • ) are shown in Figure 7. The black dashed curves show the dipolar magnetic field lines corresponding to individual starting points. The blue, green, and red color curves show the trajectories obtained for the wave frequencies of 11,905 Hz, 12,649 Hz,

  Figure 8. (a) Wave normal angles corresponding to unducted trajectories from Figure 7 at the times when the emissions reach DEMETER altitudes in the southern hemisphere are plotted as a function of the respective final geomagnetic latitudes. The blue, green, and red curves were obtained for wave frequencies of 11,905 Hz, 12,649 Hz, and 14,881 Hz, respectively. The horizontal dashed lines show the wave normal angles of respective resonance cones. (b) The same as (a), but for refractive index values.

Figure 9 .

 9 Figure 9. Example raytracing trajectories corresponding to only a partially ducted signal propagation. The waves were started at a geomagnetic latitude of 49.4 • , i.e., approximately 8 •northward from the Khabarovsk transmitter. They were assumed to propagate ducted, i.e., exactly along the geomagnetic field line all the way to the geomagnetic equator. They are then traced from the equator to the southern hemisphere, with wave normal angles of ±70 • , assuming the wave vectors lie in the meridional plane. The blue, green, and red curves were obtained for wave frequencies of 11,905 Hz, 12,649 Hz, and 14,881 Hz, respectively. The black dashed curves show the dipolar magnetic field lines.

  Figure 10, as they extend beyond the plotted wave normal angle range. However, such waves do not result in considerable Doppler shifts and are not of importance for the present analysis.

  mitter signal may split to up to three signals at different frequencies. These can be interpreted in terms of different propagation paths of the signal, resulting in different wave normal angles, refractive indices, and Doppler shifts at the observation point. The signal without an observable Doppler shift can be attributed to a ducted propagation, which occurs at low wave normal angles well away from the resonance cone. The signal with the Doppler shift steadily increasing with geomagnetic latitude until reaching a given latitudinal (frequency-dependent) threshold can be interpreted as an unducted signal(Starks

  manuscript submitted to JGR: Space Physics et al., 2009). The used raytracing analysis confirms this interpretation, achieving at least a qualitative agreement with observations.

Acknowledgments

We thank the engineers from CNES and scientific laboratories (CBK, IRAP, LPC2E, LPP, and SSD of ESTEC) who largely contributed to the success of the DEMETER mission.

DEMETER data are accessible from the https://sipad-cdpp.cnes.fr website. F. N. and O. S. acknowledge the support of GACR Grant 21-01813S.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Space Physics


Chloé Charbonné
email: charbonne.chloe@gmail.com

Maria-Lisa Dhuitte
email: maria-lisa.dhuitte@utbm.fr

Khadidja Bouziane
email: khadidja.bouziane@utbm.fr

Dominique Chamoret
email: dominique.chamoret@utbm.fr

D Candusso
email: denis.candusso@univ-eiffel.fr.

Yann Meyer
email: yann.meyer@univ-smb.fr

M.-L Dhuitte

DESIGN OF EXPERIMENTS ON THE EFFECTS OF LINEAR AND HYPERELASTIC CONSTITUTIVE MODELS AND GEOMETRIC PARAMETERS ON POLYMER ELECTROLYTE FUEL CELL MECHANICAL AND ELECTRICAL BEHAVIOUR

Keywords: Fuel cell, Mechanical contact pressure, Electrical contact resistance, Design of Experiments, Non-linear and Linear behaviour. Highlights

Design of experiments on the effects of linear and hyperelastic constitutive models and geometric parameters on polymer electrolyte fuel cell mechanical and electrical behaviour

Introduction

Proton Exchange Membrane Fuel Cells known for their high efficiency [START_REF] Luo | Overview of current development in electrical energy storage technologies and the application potential in power system operation[END_REF] -quick start at room temperature, and the variety of applications (stationary systems, mobile devices, transportation [START_REF] Kazim | Introduction of PEM fuel-cell vehicles in the transportation sector of the United Arab Emirates[END_REF] and space applications [START_REF] Giacoppo | 2 kW Modular PEM fuel cell stack for space applications: Development and test for operation under relevant conditions[END_REF]) is a promising substitute for combustion engines owing to its low level of output pollutants. Great efforts have been done toward the commercialization of PEMFC in energy market but a larger deployment involves enhancing their efficiency, lifetime [START_REF] Park | A review of the gas diffusion layer in proton exchange membrane fuel cells: Durability and degradation[END_REF][START_REF] Pei | Main factors affecting the lifetime of Proton Exchange Membrane fuel cells in vehicle applications: A review[END_REF][START_REF] Zhang | Enhancing fuel cell durability for fuel cell plug-in hybrid electric vehicles through strategic power management[END_REF], reducing their costs and improving electrical power production. The link between the energetic performance of a PEMFC and external-internal mechanical excitations [START_REF] Hou | AC impedance characteristics of a vehicle PEM fuel cell stack under strengthened road vibrating[END_REF] ,as well as environmental conditions is one of the significant factors influencing its industrialization [START_REF] Khetabi | Effects of mechanical compression on the performance of polymer electrolyte fuel cells and analysis through in-situ characterisation techniques -A review[END_REF].

The literature on these mechanical aspects and the selection criteria of cell components are relatively limited [START_REF] Saadat | Recent advances and future perspectives of carbon materials for fuel cell[END_REF]. To enhance these systems, it is thus necessary to have preliminary design models meeting the requirements of reliability and performance, reducing PEMFC cost and increasing lifespan. A preliminary design model is developed to allow the initial exploration of the design space with a limited computation cost. The early-stage design is important because the major design choices are made at this stage, in particular the product architecture. For example, this could enable to design new stack compression systems for assembly procedures and FC operation by reducing Ohmic losses in the electrical power generator. After this stage, the potentialities for product innovation or modification are quite limited [START_REF] Diampovesa | Optimal design for electromagnetic devices: A synthesis approach using intervals and constraint-based methods[END_REF]. Consequently, from the beginning of the design process, it is important to have a better understanding of the parameters that play a role in electrical power generation, such as the electrical Contact Resistance (CR) which influences the electronic losses.

Among the external excitations of a PEMFC, there are assembly procedures, operational and environmental aggressive conditions. Here, we will be interested in the mechanical clamping of stack [START_REF] Alizadeh | Development of contact pressure distribution of PEM fuel cell's MEA using novel clamping mechanism[END_REF] which is used in order to assemble different components of PEMFC. This mechanical stress influences the global FC performance [START_REF] Wang | Experimental study on clamping pressure distribution in PEM fuel cells[END_REF][START_REF] Hou | An investigation of characteristic parameter variations of the polarization curve of a proton exchange membrane fuel cell stack under strengthened road vibrating conditions[END_REF][START_REF] Zhou | Gas diffusion layer deformation and its effect on the transport characteristics and performance of proton exchange membrane fuel cell[END_REF]. To optimize the performance of a PEMFC, the mechanical excitations must be handled. The mechanical compression during mechanical clamping reduces Ohmic losses due to the resistance of the GDLs and rise of mass transport losses by reducing the porosity of the GDL.

Among the components of PEMFC, GDLs play a major role in the cell performances by allowing a homogeneous distribution of gas from the channels to the membrane. In fact, most papers treating PEMFC analysis and modelling use a linear GDL constitutive law. Articles focusing on GDLs sometimes consider a hyperelastic constitutive law [START_REF] García-Salaberri | Nonlinear orthotropic model of the inhomogeneous assembly compression of PEM fuel cell gas diffusion layers[END_REF]. However, to the knowledge of the authors, the comparison between these linear and hyperelastic behaviours is little explored. The experimental study described in [START_REF] Khetabi | Effects of mechanical compression on the performance of polymer electrolyte fuel cells and analysis through in-situ characterisation techniques -A review[END_REF] point out that "contact resistance was found to decrease non-linearly with compression: more than 75% of reduction was attained at 2.5 MPa". The representativeness of the preliminary models involves the study of nonlinearities, in particular for GDL [START_REF] Zhang | Study of the mechanical behavior of paper-type GDL in PEMFC based on microstructure morphology[END_REF]. In this paper, we are especially interested in the influence of the GDLs constitutive law (linear and hyperelastic constitutive laws) and the layer thicknesses on the contact pressure (CPRESS) on the performance of PEMFC. Geometry, especially all the layer thicknesses, plays an important role in the performance of PEMFC. The layer thicknesses influence the variation of the mechanical stresses and therefore affects performance. The PEMFC performance is associated with contact pressure CPress as in Eq. ( 1).

𝐶𝑃𝑟𝑒𝑠𝑠 = 𝑁 𝑎×𝑏 (1)

With

• N: Normal plan force [Newton],

• a×b: the real contact surface between interface [m²].

The performance in terms of power generating depends of protonic losses and Ohmic losses.

About 50% of electronic losses stem from Contact Resistance (CR) [START_REF] Netwall | Decreasing contact resistance in proton-exchange membrane fuel cells with metal bipolar plates[END_REF][START_REF] Ye | Bulk and contact resistances of gas diffusion layers in proton exchange membrane fuel cells[END_REF]. At each interface, there will be a loss of current intensity due to insufficient surface-to-surface contact. Studies have found that Mechanical CPress can strongly affect CR [START_REF] Alizadeh | Development of contact pressure distribution of PEM fuel cell's MEA using novel clamping mechanism[END_REF][START_REF] Wang | Experimental study on clamping pressure distribution in PEM fuel cells[END_REF] and a relation between CR and CPress was observed [START_REF] Hou | An investigation of characteristic parameter variations of the polarization curve of a proton exchange membrane fuel cell stack under strengthened road vibrating conditions[END_REF][START_REF] Zhou | Contact resistance prediction and structure optimization of bipolar plates[END_REF]. The most important electronic ohmic losses in FC are at the BPP -GDL interface with the electrical CR [START_REF] Hou | AC impedance characteristics of a vehicle PEM fuel cell stack under strengthened road vibrating[END_REF][START_REF] Alizadeh | Investigation of contact pressure distribution over the active area of PEM fuel cell stack[END_REF]. The aim is therefore to have the optimal contact pressure in order to reduce the CR as much as possible [START_REF] Atyabi | Effects of assembly pressure on PEM fuel cell performance by taking into accounts electrical and thermal contact resistances[END_REF][START_REF] Lai | A mechanical-electrical finite element method model for predicting contact resistance between bipolar plate and gas diffusion layer in PEM fuel cells[END_REF].

The paper goes further by using a Design of Experiments (DoE) with a 2D Finite Element (FE)

Model including contact friction to study the effect of geometry and constitutive law of materials, on the mechanical and electrical performance of a PEMFC stack. To validate and assess numerical results, we will make comparisons with experimental results from literature.

Our methodology is presented in the next section. The results of the parameters effects on the overall performance are shown in Section 3. The last section is a critical discussion of our study.

Methods

The main purpose of this study is to investigate the effects of geometry and GDL constitutive law on CPress. A 2D Finite Element (FE) model with contact friction has been developed to simulate the mechanical properties to obtain CPress and von Mises yield criterion. An effective approach is proposed to save CPU time with DoE methods, by using Taguchi approach. DoE provides also an optimal simplified model to maximize the pressures for each interface and minimize the constraints.

Finite Element Method

A FE model has been developed with the commercial software ABAQUS [START_REF] Abaqus | ABAQUS/Standard User's Manual[END_REF]. This software uses the Finite Element Method (FEM) to solve engineering and mechanical problems. The FE Modelling permits to build a representation of a problem in a model and to give the structure behaviour subjected to environmental conditions and to external / internal mechanical excitations.

Geometry and Model

A 2D axisymmetric plane strain FE model has been modelled on ABAQUS. This modeling is aimed to obtain a global mechanical behaviour and not a local one. In order to limit the computation (CPU) time, the model is simplified without BPP distribution channels. The elastomeric gaskets are not included in the model. This assumption is based on the fact that the force on the gasket used only represents 5% of the total clamping force [START_REF] Carral | A numerical analysis of PEMFC stack assembly through a 3D finite element model[END_REF]. In order to have The point of view adopted in this paper is the development of a preliminary model. The modeling assumptions result from this choice. The FE model is built with the following assumptions:

❖ Assembly pressure of 1 MPa applied uniformly on the collector plate.

❖ Boundary condition applied on BPP and EP is u1=0 (see Figure 1).

❖ Symmetry condition applied on MEA.

❖ Coefficient of adhesion friction is μ=0.3, Coulomb's law T= N.μ = N.tanφ [START_REF] Carral | A numerical analysis of PEMFC stack assembly through a 3D finite element model[END_REF].

The widths and the heights in direction x and z are constant as presented in Figure 1 and Table 1; the thicknesses in y direction are variants. Our first objective is to quantify the impact of thicknesses on CPress and so on CR. Our second purpose is to observe the consequence of the GDL constitutive law on the mechanical and electrical behaviour of a PEMFC. For this purpose, the model developed by Mishra [25] is applied to link mechanical contact pressure and electric contact resistance (and thus the PEMFC performance). Thus, it is possible to take into account analytically the geometry of the BPP and limit the loss of accuracy. A balance between computation time and accuracy was sought.

Material properties

The mechanical behaviour is included in the assembly model with homogeneous and isotropic materials. The constitutive laws applied derive from mechanical stress-strain tests. Gigos et al.

found that the GDL behaviour is linear for small deformations but non-linear in large deformations [START_REF] Gigos | Mechanical characterization and analytical modeling of gas diffusion layers under cyclic compression[END_REF]. Better still, the irreversible and non-linear nature of the hysteresis cycle depend on the number of loading-unloading cycles with a stabilization after 5 cycles. Note that the study focuses on the influence of GDL constitutive law on contact pressure with hyperelastic and linear constitutive law. The elastic material properties of EP, BPP, GDL, and MEA are listed in Table 2. The mechanical properties of the Nafion membranes (112,115 and 211) are given for 90 % RH at 40 °C [START_REF] Choi | Consideration of thermodynamic, transport, and mechanical properties in the design of polymer electrolyte membranes for higher temperature fuel cell operation[END_REF]. For simplicity, this hydration level is assumed to be constant throughout the compression process, although in real applications this may not be necessarily true. The GDL hyperelastic constitutive law was extracted from uniaxial experimental data with a Marlow deformation energy potential model [START_REF] Gigos | Mechanical characterization and analytical modeling of gas diffusion layers under cyclic compression[END_REF]. The GDL has a hyperelastic behaviour due to its microscopic structure. Therefore, using a linear constitutive law is an assumption to be carefully used. As shown in Figure 2, for a GDL deformation smaller than ε=0.14, the linear nominal stress will be lower (up to 56.12% max) than Marlow's stress (non-linear). And near to zero and for ε=0.14, the differences between these laws are slight.

Figure 2: Relation stress/strain of GDL -comparison between linear/hyperelastic laws

Meshing

In this FE modelling, the greatest difficulty is to keep a continuity in the mesh with a great variation in thickness. The MEA can be five hundred times smaller than BPP. The mesh has element shape quad dominated with CPE4RH element type a 4-node bilinear plane strain quadrilateral in linear interpolation [START_REF] Abaqus | ABAQUS/Standard User's Manual[END_REF]. The node and element numbers of the half-single fuel cell are given in Table 3. The FE number is investigated to have a good mesh quality on the entirety of the model. The aim is to have a sufficiently fine mesh so that the results are consistent and the CPU time short.

Hence the size mesh of the model is adapted in function of the variations in thickness as shown in Figure 3. To guarantee the quality of the results, we have chosen at least two elements in the thickness of each component. The compromise for the overall mesh size is between 0.01 and 2.5 depending on thickness size.

Parameterization of Abaqus

To set up our models, two steps are created to simulate an assembly procedure. One initial step at 0 "time" period and another one at 1 "time" period. In a nonlinear analysis, a step takes place over a finite period of "time", although this "time" has no physical meaning unless inertial effects or rate-dependent behaviour are important. The incrementation is automatic with a maximum number of 100 increments. The initial increment size is between 0.0001 and 0.01.

The maximum increment size is 1 and the minimum lies between 1E-30 and 1E-05. The other step parameters are the selection of the direct solver method and the full Newton solution technique. A linear extrapolation, in "time", of the previous incremental solution is used to begin with the non-linear equation solution for the current increment. The uniform pressure is 0 MPa at initial step and loaded of 1MPa at charge step. The boundary conditions are applied on BPP and EP. The BPP and EP have the boundary condition u1 blocked and a symmetry on the MEA with u2=0 as seen in Figure 1. The representations used in Abaqus for the interactions between components are: surface-to-surface contact (standard) type with master surface and slave surface. The contact interaction property is a tangential behaviour isotropic and friction coefficient about 0.3. To avoid interpenetration between each element, a non-penetration constraint is applied. The analysis is run in Abaqus/Standard with the Newton-Raphson method on 4 processors and 1 GPGPU acceleration. To ensure the convergence of the calculations, the initial increment size could be decreased.

Design of Experiments with Taguchi approach

Explanation for simplification

A DoE is a useful tool for designing a product or solving complex optimization systems. It allows saving time thanks to a series of tests organized in advance to understand and quantify multiple parameters in a minimum of tests. It avoids combining all the modalities of all the parameters. Using a DoE table, the number of tests is given in Equation [START_REF] Kazim | Introduction of PEM fuel-cell vehicles in the transportation sector of the United Arab Emirates[END_REF].

𝑚𝑜𝑑𝑒𝑙 𝑛𝑢𝑚𝑏𝑒𝑟 = 𝐿𝐸𝑉𝐸𝐿 𝑓𝑎𝑐𝑡𝑜𝑟 (2)
With:

❖ Factor: number of parameters to vary; factors are independent variables, ❖ Level: numbers of factor values.

The Taguchi approach permits to reduce model number thanks to standard tables which list the most common combinations of levels and factors [START_REF] Joneja | Design Principles[END_REF].

Modalities

The model has 4 factors: thicknesses of EP, GDL, BPP, and MEA. We chose to work with levels equal to 3. To highlight the influence of geometry with a classical DoE method, we would have had 3 4 run, so 81 different models to run. Moreover, as we also studied the influence of the constitutive law on the contact pressure, we would have two tables of DoE: one table with GDL hyperelastic constitutive law and the other one with GDL linear constitutive law.

Consequently, 162 models should have run.

For 4 factors with 3 levels, the table for the Taguchi approach is an experimental matrix L9 with 9 the model number. Tables 4 and5 respectively give the description of models and the details of the different factors and levels. The Taguchi method allows reducing the number of models from 162 to 18.

Table 4: Description of models in the Taguchi approach

Test Nr.

Levels of controlled factors Results

A

B C D 1 1 1 1 1 R1 2 1 2 2 2 R2 3 1 3 3 3 R3 4 2 1 2 3 R4 5 2 2 3 1 R5 6 2 3 1 2 R6 7 3 1 3 2 R7 8 3 2 1 3 R8 9 3 3 2 1 R9

Coupling FEM and DoE

As presented in Subsection 2, in a first step, the geometry and material characteristics are chosen

Results

The mechanical behaviours of the PEMFC components are described in the following section, which also includes a detailed comparison between our numerical results and those from the literature.

Geometry effects and most influential components

Geometry effects on the model convergence

Concerning the choice of the factor D-EP to level 3 (D3) in the selected Taguchi table, a first set of simulations with a value of 1.588 mm was achieved. Tests 3, 4, 8, and 9 could not converge. In the analysis of the inconclusive tests, we observed that there could not be a smaller End Plate (EP) than the BPP. On the one hand, with a BPP thickness below than 10% of the EP thickness, the numerical calculations did not run. We assume that is due to a significant penetration of the materials with each other. Thanks to a new value of 7.5 selected for D3, we were able to run Tests 3, 4, and 8. On the other hand, regarding Test 9, we concluded that a very thin MEA in front of the GDL and/or EP has a negative impact on the convergence of the calculation. Indeed, in Test 9, the MEA-GDL thickness ratio (7%), and especially the MEA-EP thickness ratio (0.08%), are smaller compared to the thickness ratios of the other tests. There are no results for the Tests 9 because their thickness ratios have remained disadvantageous for the calculation in Table 6. By analyzing the inconclusive tests, it is possible to establish preliminary design rules for the numerical model (and maybe the physical stack) based on thickness ratios: the MEA/GDL ratio must be greater than 9%, the MEA/EP ratio between 0.7% and 1.7%, the GDL/EP ratio less than 12% and the BPP/EP ratio less than 190%. A detailed study on the thickness ratios effects on the numerical model would surely be very interesting.

Geometric effects on the mechanical behaviour

First, let us note that we aim to minimize the von Mises yield criterion in order to limit the load on the components and thus avoid the definitive failure of the materials. On the other hand, the contact pressure at the interfaces needs to be maximized since the more surface will be in contact, the better the contact resistance and therefore the FC performance.

Table 6 shows the maximum von Mises yield criterion and the maximum contact pressures at the interfaces of our models with a hyperelastic constitutive law for the GDL. The remarkable values come from Tests 2 and 8. In Test 2, the smallest von Mises yield criterion is between the EP and the BPP: 2.29 MPa. In Test 8, the highest stress pressure is between the BPP and the GDL: 115.2 MPa.

By only changing the thicknesses, there is a variation of 97.9% between extreme values of von

Mises yield criterion at the EP-BPP interface and 95.1% for the contact pressure at BPP-GDL interface. This observation confirms that thickness plays a key-role in the performance of the FC.

Most influential constitutive components

In Table 7 summarizing the factor effects, the EP factor (called D in the DoE) has the highest numerical values. So the EP thickness has a greater influence on stress and contact pressure.

The EP thickness controls the behaviour of stress and pressure in the other components. A thin EP considerably dramatically increases the CPress and especially the stress as shown in Table 7. It is certain that the EP is an important element to obtain a maximum contact pressure. The smaller thickness of the EP (7.5 mm) considerably increases the contact pressure and especially the stress. For EP thickness of 30 mm, there is an effect of 54.29 MPa on the contact pressure and 56.08 MPa on von Mises yield criterion.

To find the best EP thickness, it is therefore necessary to choose an EP thickness that can increase the contact pressure without overstraining the EP and the BPP. It will certainly be necessary to look for the best trade-off in the interval [7.5mm; 15mm].

Effect of the GDL constitutive law

The comparison between the linear and hyperelastic constitutive laws is presented in Table 8.

For the linear constitutive law, there are high reductions or little rises of the stresses and pressures, comparatively to the tests done with hyperelastic constitutive law. This is due to the different evolutions of the constitutive law in linear and hyperelastic modes. Surprisingly, the change in the GDL constitutive law is not directly related to its thickness. It appears that the influence of thickness ratios between them are predominant for the stress. It seems that the contact pressure is influenced by the thickness of the EP. When it is small or equal to 15, there is a decrease at the GDL-MEA and EP-BPP interfaces, and an increase at the BPP-GDL interface, in comparison with the results of the hyperelastic behaviour law. When the thickness of the EP is 7.5mm, another distribution can be observed.

The GDL constitutive law does not influence the characteristics of the best thicknesses which optimize the FC performance. On the contrary, the constitutive law affects the extremum of the DoE according to the thickness ratio. This confirms that the GDL constitutive law plays a role in the FC performance. Nevertheless, it seems to have less impact than the variation in thickness.

Behaviour of compression and stress at interfaces

In order to understand the performance, the contact between the interfaces is studied with the von Mises yield criterion. The behaviour of CPress and stress changes according to the thickness ratios and the GDL constitutive law for a given interface. For the EP-BPP interface, the simulations led to inhomogeneous but symmetrical behaviours along the surface, as shown in Figures 5 and6. There is a pressure drop on the edge and in the middle of the interface. A pressure drop means a bad contact between the two interfaces, leading therefore to a debonding. This bad contact can cause an increase in the electric contact resistance inside the fuel cell and therefore induce a loss of performance. As presented in Figure 5, the pressures vary between 5 and 22 MPa as a function of the thickness ratios.

Overall, for the tests with low pressures under 13.5 MPa, there is a correlation between pressure and stress as presented in Figures 5 and6. When the pressures are low, so are stresses. This is the case for Tests 1, 6, and 7. Others, such as Tests 2, 3 and 6, have medium pressures, and the stresses extend from the maximum pressure area to the centre of the interface. This can sometimes be noticed for a surface contact of 40 mm, as it is the case for Test 2. This suggests that there is slippage and shear between the EP and BPP interfaces. But in the middle of the edge, the mechanical pressures are close to 0. Overall, these pressures range from 0 to 118 MPa. It is two times more than the maximum stress of the BPP-GDL interface. The behaviours in pressure and stress at this interface are similar, except at the ends.

Further analysis of the tests shows that for each test the highest value of pressure and stress is at the edges of interface, as depicted in Figures 7, 8, 9 and 10. In the range [0 ; 5 mm], the contact pressure is greater than the stress. After 5 mm, the values are near to 0. The contact behaves in a similar way for the tests: with an overpressure at the edges, then with a depression close to the edge, and a new overpressure until approaching 0 after 5 mm as in Figure 7. The behaviour of the stresses is much more fluctuating between the different tests, with the same overstress and under-stress pattern as the contact pressures, except for Test 7

given in Figure 9. For example, in Figure 9 and for Test 7, there is no more under stressing close to the edge and for Test 9, there is considerable overstressing close to 4 mm. There is therefore a real interest in choosing thickness ratios that maximize contact pressures and minimize stress, in order to reduce electronic losses as much as possible and to avoid any mechanical rupture in the PEMFC. Therefore, we have a pressure concentration at the ends and then a bonding and a debonding in the middle of the structure. This creates a high variation in resistance throughout the structure.

A closer examination at the ends shows that, for this interface, the behaviour of the stresses is identical to the pressures. When the pressure is at its highest, the stresses are at their highest.

Even when it is minimal, so are the stresses. This observation is similar for each test.

Best set of parameters

According to our effect analysis table, the best EP thickness seems to be in the [7.5mm; 15mm] interval. We chose the thickness of 15 mm (D2) with the least negative influence on the contact stress and pressure as presented in Table 5. The best compromise for the MEA is the thickness of 0.05 mm (A1) in Table 4. The contribution of the thickness of the MEA to the performance of the FC is quite significant after the EP. As far as GDL is concerned, it is certain that the best choice is that of B2 and therefore the thickness 0.28 mm which will reduce the stress by maximizing the contact pressure. Finally, for the BPP, the best parameter according to Table 7 is C1 with a thickness equal to 3 mm.

Our most unfavourable case is with the linear behaviour because the lowest value of von Mises yield criterion increases and drastically reduces the greater value of CPress, compared to the tests with the hyperelastic constitutive law. That is why for the model of best sets of parameters, the constitutive law for GDL is hyperelastic.

An optimal is performed with 182 189 nodes and 161 725 elements. The maximum values of von Mises yield criterion do not exceed the elastic limit that each component can withstand as given in Tables 9 and10.

The von yield criterion in contact the GDL-MEA and EP-BPP it allows to have a high contact pressure.

Effect of the pressure distribution on the electrical contact resistance

The work of Bouziane et al [START_REF] Bouziane | Impact of cyclic mechanical compression on the electrical contact resistance between the gas diffusion layer and the bipolar plate of a polymer electrolyte membrane fuel cell[END_REF] makes possible to compute the electrical contact resistance between GDL and BPP. The TLM method is used to measure the contact resistance under 5 compression cycles over several mechanical pressure stages ranging from 0 to 8 MPa. The most interesting values to compute the contact resistance is obtained for the loading phase applied during the first cycle. We focused on the data linked to the GDL ref.

Toray H90 with a 0.28 mm thickness. From Figure 5, it is possible to obtain the distribution of contact resistances shown in Figure 15 by applying the approach described in [START_REF] Mishra | Measurement and prediction of electrical contact resistance between gas diffusion layers and bipolar plate for applications to pem fuel cells[END_REF][START_REF] Zhang | Estimation of contact resistance in proton exchange membrane fuel cells[END_REF]. The electrical contact resistance of the BPP-GDL interface is a non-linear function with respect to the thickness of the BPP. The model can also be used to compute the electrical resistance of the BPP-GDL interface: 3.49 mΩ for Test 2, 3.09 mΩ for Test 5 and 19.9 mΩ for Test 8.

Comparison with literature

To verify the validity of our model, we compare the pressure distributions on GDL and MEA from experiments and numerical models.

Experimentally, there are two major in-situ techniques used to measure contact pressure distributions [START_REF] Khetabi | Effects of mechanical compression on the performance of polymer electrolyte fuel cells and analysis through in-situ characterisation techniques -A review[END_REF]. The first technique uses pressure sensitive films composed of two polyester based sheets which one of them is coated with a layer of micro-encapsulated colour forming material and the other one is coated with a layer of colour developing material [START_REF] Lin | Experimental study of clamping effects on the performances of a single proton exchange membrane fuel cell and a 10-cell stack[END_REF]. Pressure is approximated depending on the density of red spots. The second tool for evaluating the pressure distribution is the use of piezoelectric sensors, as in the work of Dey et [START_REF] Dey | Study of PEM fuel cell end plate design by structural analysis based on contact pressure[END_REF].

Most the works are not usable since the clamping system is with nuts and bolts. The distributed pressure is therefore not uniform. As pointed out in the work of Wen et al. [START_REF] Lin | Experimental study of clamping effects on the performances of a single proton exchange membrane fuel cell and a 10-cell stack[END_REF],

positions of nuts and their number have considerable influence on the pressing process. This is partly due to the buckling of the endplate. [START_REF] Alizadeh | Investigation of contact pressure distribution over the active area of PEM fuel cell stack[END_REF].

The experimental work of Alizadeh

In their work, Bates et al [START_REF] Bates | Simulation and experimental analysis of the clamping pressure distribution in a PEM fuel cell[END_REF] have obtained the same behaviour of the von Mises constraints, although the clamping system does not match our model. In their model, the clamping system consists of nuts and bolts. They have performed two simulations: a single-cell and a 16-cell stack were run under various clamping pressures. In both models, a linear behaviour law is applied to the GDL and to all other elements. Only the single-cell model is used and compared to our test with the closest thicknesses. The same stress distribution is seen on the EP/BPP and GDL/MEA interfaces as shown in Figure 17, but the presence of the BPP channels considerably changes the mechanical behaviour of the BPP/GDL interface. The load corresponding to our model corresponds to a uniform load distribution. The numerical results of the pressure distribution on the BPP/GDL interface are in line with the results obtained in this paper although there are channel print marks.

Conclusions

The article presents a study of the mechanical behaviour of a PEMFC with a focus on the contact pressure distribution and electrical contact resistance at different layer interfaces in the stack.

A numerical development was initiated in order to quantify the influences of thickness of cell components and constitutive law of the GDL on the PEMFC performances expressed in terms of electrical resistances at the early-stage design of a PEMFC. To reach this aim, a 2D finite element model including contact friction was developed and combined with a DoE. Concerning the GDL, a hyperelastic experimental constitutive law and a linear constitutive law were implemented. An optimal set of geometric and mechanical parameters were determined. During the computation process, it appeared that the ratios between the thicknesses of the cell components influence the computation stability. After reviewing the simulation results, we noted that the electrical and mechanical responses are more impacted by the layer thicknesses than by the GDL constitutive law.

  an efficient and effective model, a single fuel cell was simulated. Due to symmetry of a cell, we had only modelled half of the single fuel cell. These assumptions allow reduced computation times. The model components are EP (End Plate), BPP, GDL, and MEA as shown in Figure. 1.

Figure 1 :

 1 Figure 1: Model of a single fuel cell

Figure 3 :

 3 Figure 3: Mesh used for computations

  in the literature dedicated to PEMFCs[29,30,31,32]. In a second step, the mechanical properties given in subsections 2.1 and 2.2 are simulated with the software ABAQUS[START_REF] Abaqus | ABAQUS/Standard User's Manual[END_REF]. Finally, as explained in Subsection 2.2, the DoE factors are extracted following the flowchart in Figure4. This methodology implied to simulate a configuration of mechanical properties individually. 9 tests have been modelized. For every model, CPress and the von Mises yield criterion are extracted for each step of increment. All the curves presented in this paper are related with the loading step at 1 s.

Figure 4 :

 4 Figure 4: Methodology for coupling FEM and DoE

Figure 5 :Figure 6 :

 56 Figure 5: CPress at the EP-BPP interface with hyperelastic GDL

Figure 7 : 8 Figure 8 :Figure 9 : 8 Figure 10 :

 7889810 Figure 7: CPress at the BPP-GDL interface with hyperelastic GDL

Figure 11 :Figure 12 : 8 Figure 13 :Figure 14 :

 111281314 Figure 11: CPress at the GDL-MEA interface with hyperelastic GDL

Figure 15 :

 15 Figure 15: Distribution of the electrical contact resistance at the GLD-BPP interface for the GDL ref. Toray H90

  [START_REF] Alizadeh | Investigation of contact pressure distribution over the active area of PEM fuel cell stack[END_REF] makes possible to validate the stress distribution for the MEA and the pressure distribution on the GDL. The loading conditions are identical to our model. They used pneumatic clamping system which can compress the cell at 7 bars. The contact pressure distribution over the active area of the PEMFC is obtained using super lowpressure Fujifilm used with a range of [0.5; 2.5] MPa inserted between two GDLs as depicted in Figure16. Their work confirms the behaviour obtained on the MEA and the GDL for both the numerical model and the experimental test. However, the values of pressures and stresses are not similar. This may possibly be explained by the hyperelastic constitutive law used to model the GDL behaviour.

Figure 16 :

 16 Figure 16: Contact pressure distributions over the MEA of PEM Single cell: (a) Stainless steel end plate with 30 mm thickness, (b) Stainless steel end plate with 50 mm thickness [20].

Figure 17 :

 17 Figure 17: Von Mises stress distribution (single cell, active area = 100 cm 2 ): a) deformed GDL stress plot, (b) deformed MEA stress plot [37].

Table 1 :

 1 Description of the model geometry

		Width [mm] Height [mm]
	EP	350	350
	BPP	295	295
	GDL	200	200
	MEA	200	200

Table 2 :

 2 Mechanical properties of the PEMFC components

	Cell component	Material	Young modulus [MPa] Poisson's ratio
	EP	Aluminium alloy	70 000	0.33
	BPP	Graphite	5 100	0.25
	GDL	TGP-H-90	60	0.25

Table 3 :

 3 Mesh description

	Test	Node numbers	Element numbers
	1	142 294	121 976
	2	181 601	161 284
	3	11 602	8 735

Table 5 :

 5 Description of the levels (thicknesses) in the DoE

	Factors	Level 1	Level 2	Level 3
		Nafion-112	Nafion-115	NR 211
	A -MEA
		0.05	0.127	0.0254
		TGP-H-060	TGP-H-090	TGP-H-120
	B -GDL
		0.19	0.28	0.37
			Bipolar Graphite Plate
	C -BPP
		3	6.35	12.7
	D -EP	30	15	7.5

Table 6 :

 6 Maximum pressure and von Mises yield criterion for hyperelastic GDL

	Hyperelastic GDL

Table 7 :

 7 Total effects of the DoE study with hyperelastic GDL

			Total effect on von Mises	Total effect on contact pressure
			yield criterion [MPa]	[MPa]
		A1	-28,74	-21,96
	MEA	A2	16,63	-3,96
		A3	18,18	38,88
		B1	14,93	-9,42
	GDL	B2	-4,38	16,12
		B3	-15,82	-10,04
		C1	-2,64	18,12
	BPP	C2	37,16	11,43
		C3	-22,13	-25,74
		D1	-39,51	-40,27
	EP	D2	-29,74	-27,45
		D3	56,08	54,29

Table 8 :

 8 Difference for maximal pressure and maximal von Mises yield criterion between the

			hyperelastic and linear laws
			Difference between hyperelastic and linear laws
	Results	von Mises yield criterion [MPa]	Contact pressure [MPa]
		Contact	Contact BPP-	Contact EP-	Contact GDL-	Contact BPP-	Contact EP-
		GDL-MEA	GDL	BPP	MEA	GDL	BPP
	R1	2.17	-4.79	-1.74	2.32	-2.94	2.39
	R2	2.53	-3.67	-5.85	3.50	-1.26	2.57
	R3	3.67	-2.14	0.65	5.74	0.21	-0.415
	R4	5.39	10.83	-0.84	6.35	23.17	0.46
	R5	9.31	-4.89	-2.18	2.76	-4.16	4.09
	R6	3.19	-6.92	3.23	3.97	-4.65	3.61
	R7	4.21	4.02	0.08	2.05	-9.06	0.13
	R8	10.73	21.15	4.6	15.33	52.6	4.35

Table 9 :

 9 Max von Mises yield criterion and CPress

		Optimal results
	Contact type	Max. von Mises yield criterion [MPa]	Max. contact pressure [MPa]
	Contact GDL-MEA	15	16
	Contact BPP-GDL	6	16
	Contact EP-BPP	29	13

Table 10 :

 10 Elastic limits of different constitutive materials

	Component	Material	Elastic limits (MPa)
	End Plate (D)	Aluminium alloy	220 [29]
	BPP (C)	Graphite	50 [30]
	GDL (B)	Carbon paper (TGP-H-90)	40 [31]
	MEA (A)	Nafion	23-32 [32]

  increases the contact resistance, thus decreasing the contact pressure. Adopting thin components increases the contact resistance. Therefore, the most favourable case seems to have thick and equivalent thicknesses between the different components.

				35,00
						Test 2 (6.35
	electrical contact Specific	resistance	[mΩ.cm²]	25,00 30,00		mm thick BPP) Test 5 (12.7 mm thick BPP) Test (3mm thick
				0,00	50,00	100,00	150,00	200,00
					Contact surface [mm]

Acknowledgements:

The "Région Bourgogne-Franche-Comté" is gratefully acknowledged for its support through the ELICOP Project (Ref. 2015C-4944 and 2015-4948). The French "Ministère de l'Enseignement supérieur, de la Recherche et de l'Innovation" is gratefully acknowledged for the funding of Khadidja Bouziane's PhD Thesis (Contrat doctoral Univ. Paris Sud n°2017-225).


Dr Ulrich Baisch

Prof Rinaldo Poli

Dr Liming Tang

Ms Julia L Conway

Prof M Kevin

Smith

Keywords: Chromium, homolytic bond cleavage, ligand stereoelectronic effects, controlled radical polymerization, OMRP, poly(vinyl acetate)

Compounds 3-5: Compounds 3, 4 and 5 were prepared using the same procedure as for compound 2, above, but using the DppNHC(Me)CHC(Me)NAr ligands, where Ar = C6H5, p-C6H4OMe, or p-C6H4CF3, respectively. Controlled polymerizations. A general procedure was as follows. Under strict exclusion of air and moisture, the reaction components (e.g. metal complex and V-70 for the OMRP runs) were placed in a Schlenk tube equipped with a stirring bar. Freshly distilled monomer and the solvent (when this was used) were added at 0 °C, followed by three freeze-pump-thaw cycles. The tube was then placed in a preheated oil bath at the desired temperature. Samples were periodically withdrawn by glass syringe after quenching the solution with an icebath. The conversion was calculated by weight difference, after removing all residual monomer by evaporation under vacuum to constant weight. The polymer residue was directly analyzed by size exclusion chromatography after dissolution in THF.

Chromium(III)-alkyl bonds

DOI: 10.1002/anie.200((will be filled in by the editorial staff))

Homolytic bond strengths and formation rates in half-sandwich chromium-alkyl complexes: relevance for controlled radical polymerization** Yohan Champouret, Ulrich Baisch, Rinaldo Poli,* Liming Tang, Julia L. Conway, Kevin M. Smith,* † This last decade, controlled/living radical polymerization (CRP) processes have seen a considerable surge of interest due, in part, to their relevance to access to a variety of well-defined polymer structures (e.g. predetermined molecular mass, narrow molecular weight distribution). [START_REF] Matyjaszewski | Macromolecular Engineering: Precise Synthesis, Materials Properties, Applications[END_REF] We have been interested in one-electron reactivity of transition metal complexes and its relevance in CRP. [START_REF] Poli | [END_REF] One way in which transition metal complexes can be used to control radical polymerization is through a reversible deactivation consisting of the formation of an organometallic dormant species, a metal-capped polymer chain, by trapping the growing radical chain with formation of a metal-carbon bond, see Figure 1. We refer to this particular control mechanism as "organometallic radical polymerization" (OMRP). [START_REF] Poli | [END_REF] One of the outstanding challenges in this area is the possibility to control the polymerization of less reactive monomers (e.g. vinyl chloride, vinylidene dichloride, vinyl acetate, etc.), for which activation is made difficult by the relative strong bonds established with the common radical traps.

Reasonable control for the radical propagation of poly(vinyl acetate) (PVAc) has been achieved on the basis of another control mechanist (degenerative transfer, DT, based on the use of xanthates or dithiocarbamates). [3,4] Results obtained by atom transfer radical polymerization (ATRP) have not been nearly as good, [5][6][7] while good control was recently achieved in the presence of Co(acac)2 (Mw/Mn as low as 1.1). [8,9] Recent studies carried out partly in one of our laboratories have shown that this process occurs either by DT or by OMRP depending on the presence of additional ligands such as pyridine or water. [10,11] <Figure 1>

In search for other transition metal complexes capable of controlling the polymerization of less reactive monomers, we have considered the use of half-sandwich -diketiminato system of Cr II , or the type shown in Scheme 1. Previous studies have shown that stable Cr III complexes with methyl ligands can be prepared by oxidation of 1 with silver triflate, followed by alkylation with methyl lithium, [12] similar to the preparation of related chromium diketiminato organometallic complexes. [13][14][15] However, attempts to synthesize half-sandwich Cr III complexes with larger alkyl ligands led to unexpected products presumably resulting from homolysis of the Cr III -alkyl bond and subsequent hydrogen atom abstraction reactions of the solvent. [16] Thus, this metal system shows promise for applications in OMRP. In this report, we will show by a combination of experimental and computational methods that the Cr III -PVAc bond strength is dramatically influenced by the steric effect of the -diketiminato aryl substituents and can be tuned to a suitable range for the CRP of the vinyl acetate monomer (VAc). Early reports of the use of Cr II acetate and benzoyl peroxide for controlled radical polymerization were discussed on the basis of what we now call an OMRP mechanism, but the polymerization activity was low for vinyl acetate. [17][18][19] In addition to the already known [12] compound CpCr II [Ar 1 NC-(CH3)CHC(CH3)NAr 2 ] with Ar 1 = Ar 2 = iPr2C6H3 (Dipp), 1, we have now prepared the analogue with Ar 1 = Ar 2 = 2,6-Me2C6H3 (Xyl), 2, as well as asymmetric analogues with Ar 1 = Dipp and Ar 2 = Ph (3), p-C6H4OMe (4), and p-C6H4CF3, (5). The structures of complexes 2, 4 and 5 have been confirmed using single crystal Xray diffraction. Synthetic and structural details are provided in the supporting information.

<Scheme 1>

Initial studies showed that compound 2 traps the growing polystyrene (PS) radical chain inefficiently. Indeed, the V-70initiated polymerization in the presence of 2 [bulk; styrene/V-70/2 = 250:0.8:1; V-70 = 2,2'-azobis(4-methoxy-2,4-dimethylvaleronitrile)] gave PDI = 2.5-3.7, and a final polymer with Mn = 7.9•10 4 g/mol (Mn(th.) = 1.7•10 4 g/mol) at 65 % conversion (additional data are presented in the Supporting Information). Since the growing PVAc radical is expected to establish a stronger bond that the growing PS radical, [20,21] we anticipated a more efficient radical trapping for the growing PVAc radical chain.

In agreement with the above expectation, a preliminary screening of VAc polymerization initiated by V-70 in the presence of 2 at 50°C (VAc/2/V-70 = 500/1/0.8) yielded essentially no polymer in 4 h (ca. 6 half-lives of V-70). Subsequent warming to 90°C yielded an initial polymerization (11% conversion after an additional 4 h), and then essentially no further conversion in the following 60 h (Mn = 1.12•10 4 g/mol, Mw/Mn = 1.68 at 12% conversion; Mn(th.) = 5.16•10 3 ), see Table 1, showing a weak propensity for the Cp[XylNC(CH3)CHC(CH3)NXyl]Cr III -PVAc bond to break homolytically.

<Table 1>

With the idea in mind that the Cr III -PVAc BDE could be weakened by a greater ligand steric bulk, we turned our attention to compound 1. Indeed, when an experiment was carried out with complex 1 at 50°C in the presence of V-70 (VAc/1/V-70 = 500/1/0.8), the conversion reached 70 % after only 20 h, but control of the polymerization was poor (Mn = 5.67•10 4 g/mol; Mw/Mn = 2.5; Mn(th.) = 3.01•10 4 g/mol), suggesting that the reversible trapping equilibrium of Figure 1 is not sufficiently displaced toward the dormant state under these conditions. Upon running the same process at lower temperature (VAc/1/V-70 = 500/1/0.8 at 30°C), the polymerization was slower, see Table 1 (kapp = 3.3•10 -6 s -1 ). Although the level of control is not yet ideal (at 70% conversion: Mn = 6.73•10 4 g/mol, Mw/Mn = 1.80, Mn(th.) = 3.01•10 4 g/mol), the continuous growth of Mn demonstrates the occurrence of an OMRP process. When the same experiment was run in the presence of a larger excess of V-70 (VAc/V-70/1 = 500:1.5:1), the polymerization was faster an uncontrolled. This result clearly shows that, contrary to the Co(porphyrin) [22] and Co(acac)2 [10] systems, the present Cr II complex is not capable of mediating an associative DT process, consistent with the absence of vacant sites on the Cr atom in the Cr III -capped dormant chain to promote an associative radical exchange.

The effect of the ligand steric encumbrance (Xyl vs. Dipp) and of the alkyl nature (PS vs. PVAc) on the Cr III -C BDE was probed by DFT calculations on system CpCr III [ArNC(CH3)CHC(CH3)NAr](R), where Ar = Ph or Xyl and R = CH2Ph, CHMePh (a model of the PS growing chain) and CHMeOOCMe (a model of the growing PVAc chain). All calculations were carried out by full QM, thus we avoided the larger Dipp-substituted -diketiminato system, the steric effects being conveniently probed by the comparison of Ph and Xyl systems. Amongst different rotamers, the lowest energy ones are those where the Ph or OOCMe group is oriented endo relative to the Cp ring, as shown in Figure 2. The optimized geometries of the separate CpCr II [ArNC(CH3)CHC(CH3)NAr] and R, as well as fuller details on the DFT calculations, are given as supporting information.

<Figure 2>

The calculated BDE's are reported in Table 2. They show, in agreement with the experimental evidence, a stronger Cr III -R bond for the PVAc model than for the PS model, whereas the Cr III -CH2Ph bond has an intermediate strength. Introducing o-Me substituents on the Ph rings of the -diketiminato ligand considerably weakens the Cr III -R bonds (by 7.5 kcal mol -1 for the secondary CH2Ph group and by 9.8 and 8.7 kcal mol -1 for the tertiary PS and PVAc models).

<Table 2>

The effect of steric bulk on the Cr III -C BDE is further revealed by trends in the bond lengths. Addition of the o-CH3 groups to the -diketiminato Ph substituents lengthens the bond for all R. Addition of the -CH3 group to R (going from benzyl to 1phenylethyl) also lengthens the bond. Conversely, replacement of Ph with OOCMe (going from the PS to the PVAc model) strengthens the bond, this being an electronic rather than a steric effect related to the smaller delocalization of the radical spin density. These results are in perfect harmony with the inefficient PS trapping and with an irreversible PVAc trapping by the Xyl system. Extrapolation of the calculated BDEs to the Dipp system rationalizes the reversible trapping observed for the radical polymerization of VAc, showing signs of controlled growth. On the basis of these results, we are now developing new CpCr II [Ar 1 NC-(CH3)CHC(CH3)NAr 2 ] systems with intermediate bulk.

Preliminary experiments in the presence of the isosteric complexes 3-5 indicated that the steric pressure of one Dipp and one Ph substituent does not sufficiently labilize the Cr III -PVAc bond, but revealed an unexpected electronic effect in the radical trapping rate, see Figure 3. Use of compound 3 yielded almost no polymer at 90°C (4% conversion after 30 h), although the Mn and Mw/Mn were low (5.6•10 3 and 1.15, respectively) in agreement with a controlled process. Compound 4 gave 12% conversion within the initial 5 h at 60°C, then no further increase after 40 h at 80°C. However, further heating to 100°C yielded an increased conversion (40% after an additional 65 h). Finally, in the presence of complex 5 an initial burst of polymerization (41% conversion) occurred within the first 2 h at 50°C, but then conversion reached a plateau at 50% for the next 52 h and subsequent heating to 80°C took the conversion up to 71% after an additional 15 h. It seems that the para substituent on the group (X in Scheme 1) has an important effect on the alkyl radical trapping ability of the Cr II complex (barrier to bond formation in Figure 1). While the complexes containing the Xyl, Dipp and Ph substituents (1, 2 and 3) are able to trap the growing PVAc radical chain rather effectively, those containing the p-C6H4OMe and especially the p-C6H4CF3 substituents (4 and 5) are much less efficient. The origin of this electronic effect is currently unknown; it is the subject of continuing investigations.

<Figure 3>

In conclusion, we have shown for the first time that a Cr II /Cr III system can provide a platform for controlled radical polymerization of less reactive monomers. The Cr III -C bond strength can be tuned by accurate choice of the steric bulk on the -diketiminato ligand. However, electronic factors affecting the rate of radical trapping must also be considered for the optimization of the OMRP process.

Experimental Section

Ligands: The symmetric β-diketiminato ligand, XylNHC(Me)-CHC(Me)NXyl (Xyl = 2,6-Me2C6H3) was prepared according to the literature procedure. [23] The mixed N-aryl β-diketiminato ligands, DppNHC(Me)CHC(Me)Ar (Ar = C6H5, p-C6H4OMe, or p-C6H4CF3), were prepare by reacting the appropriate aniline with DppNHC(Me)CHC(Me)O, according to the literature procedure previously used to prepare the corresponding Ar = o-C6H4OMe derivative. [24] Compound 2: In a glovebox, XylNHC(Me)CHC(Me)NXyl (2.1834 g, 7.124 mmol) dissolved in 20 mL of THF was reacted with BuLi (3.6 mL of 2.0 M solution in pentane, 7.2 mmol) at -30 ºC. In a separate Schlenk flask, NaCp (3.6 mL of 2.0 M solution in THF, 7.2 mmol) was added to CrCl2(tmeda) [25] (1.7034 g, 7.124 mmol) suspended in 25 mL of THF at 25 ºC. After 30 min, the yellow solution of the deprotonated ligand was added dropwise to the reaction mixture in the Schlenk flask. After the reaction mixture had stirred overnight at 25 ºC, the solvent was removed in vacuo, and the residue was extracted with pentane and filtered through Celite. The solvent was again removed in vacuo, and the residue was extracted into a minimum of pentane, filtered, and cooled to -30 ºC overnight to yield 1.2185 g of black crystals (41.7% yield).

Legends

Scheme 1. Compounds used in this study. Behavior of complexes 3 (circles; VAc/V-70/3 = 250:0.9:1), 4 (triangles; VAc/V-70/4 = 500:0.9:1) and 5 (squares; VAc/V-70/5 = 250:0.9:1) as OMRP traps for PVAc. Further details are given as Supporting Information.

Tables

Table 1. Radical Polymerization of vinyl acetate initiated by V-70 with complexes 2 and 1.

Complex t (h) Conv (%)

Mn SEC [a] Mn(th) [b] Mw/Mn [a] Determined by size exclusion chromatography (SEC). [b] Mn(th) = ([M]0/[Cr]0)*Mmono*conv. [c] VAc/2/V-70 = 500/1/0.8; 50°C for 4h, then 90°C.

[c] VAc/1/V-70/ = 500/1/0.8; 30°C. The steric properties of the aryl substituents in chromium -ketiminato complexes can be tuned in order to achieve reversible radical trapping of a growing poly(vinyl acetate) radical chain.

Entry for the Table of Contents

Chromium(III)-alkyl bonds

Graphical material

Figure 1 .

 1 Figure 1. Free energy profile of controlled radical polymerizations by OMRP.

Figure 2 .

 2 Figure 2. DFT-optimized structures of CpCr[ArNC(CH3)CHC(CH3)NAr](R) compounds.

Figure 3 .

 3 Figure 3.Behavior of complexes 3 (circles; VAc/V-70/3 = 250:0.9:1), 4 (triangles; VAc/V-70/4 = 500:0.9:1) and 5 (squares; VAc/V-70/5 = 250:0.9:1) as OMRP traps for PVAc. Further details are given as Supporting Information.

  Baisch, Rinaldo Poli,* Liming Tang, Julia L. Conway, Kevin M. Smith,*__________ Page -Page Homolytic bond strengths and formation rates in half-sandwich chromium-alkyl complexes: relevance for controlled radical polymerization

  Scheme 1

Figure 3

 3 Figure 2

Table 2 .

 2 Calculated Cr III -C Bond Dissociation Enthalpies (BDE) a in kcal mol -1 and bond lengths in Å for compounds CpCr III [ArNC(CH3)-CHC(CH3)NAr](R). 31G**; ZPVE and PV corrections on the basis of the ideal gas model.

		Ar = Ph	Ar = Xyl
	R	BDE	Cr-C	BDE	Cr-C
	CH2Ph	20.8	2.136	13.3	2.146
	CH(Me)Ph	11.8	2.173	2.0	2.197
	CH(Me)OOCMe	28.4	2.109	19.7	2.124
	a B3LYP/6-

RP thanks the ANR (contract No. NT05-2_42140) and IUF for financial support and the CICT (Project CALMIP) for granting free computational time, and KMS thanks NSERC of Canada for financial support and Dr. Brian O. Patrick (UBC) for performing the single crystal X-ray structural determinations and refinements.

Received: ((will be filled in by the editorial staff))