Marianne Gaborieau

Robert Graf

Hans Wolfgang Spiess
email: spiess@mpip-mainz.mpg.de

Prof Wolfgang Hans

Spiess

Versatility of the dipolar filter selection: From 1H nuclear spin diffusion experiment to the measurement of nuclear Overhauser effect in homopolymer melts

Keywords: solid state NMR, dipolar filter, spin diffusion, nuclear Overhauser effect, NOE, local dynamics, homopolymer melt, poly(ethyl methacrylate), polymethacrylate

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

Introduction

Heterogeneous molecular mobility in polymeric systems plays an important role in various macroscopic application properties. Several solid-state NMR methods can be used to characterize these systems. [START_REF] Schmidt-Rohr | Multidimensional solid-state NMR and polymers[END_REF] For instance, the dipolar filter [START_REF] Egger | [END_REF] selects 1 H magnetization according to local dipolar dephasing, which corresponds to the site mobility. Combination with a conventional exchange experiment results in a so-called "spin diffusion" experiment. [3] The term "spin diffusion" describes the transport of spin polarization between spatially separated equivalent spins [4] and it usually occurs without material transport. The original idea of a spin diffusion experiment was introduced by Goldman and Shen. [5] Apart from local mobility, the selection in spin diffusion experiments can be done according to various other criteria, including T1ρ relaxation time or 1 H chemical shift. [6] [START_REF] Schmidt-Rohr | Multidimensional solid-state NMR and polymers[END_REF][8] The 1 H nuclear spin diffusion with dipolar filter has been applied already to numerous polymer samples, including block copolymers, [3] blends, [3] core-shell particles [9][10] and conetworks (polymer chains covalently bonded by blocks of another polymer) [11]. The different samples had in common a phase separation on the nanometer length scale in the material, leading to the formation of two phases composed of one homopolymer each. 1 H nuclear spin diffusion with dipolar filter allowed to quantify the size of these structures.

Heterogeneous molecular mobility was recently revealed in poly(n-alkyl methacrylate) melts. X-ray scattering studies [12][13] [14] as well as solid-state NMR studies [15] of the melt state of poly(n-alkyl methacrylates), and in particular poly(ethyl methacrylate), indicated a tendency of local nanophase separation between the polar, less flexible backbones and the non polar, more flexible side chains. The experimental results are consistent with an alternation of main chain layers and side chain layers a few monomeric units long as limiting structure, resulting in organized domains of the size of a few nanometers. This structure could be associated with a dynamic contrast between more mobile side chains rich areas and the less mobile main chain dominated parts of the sample.

Therefore, in the present work, poly(ethyl methacrylate) melts were investigated by a combination of dipolar filter with conventional exchange experiment. Poly(ethyl methacrylate) melts represent a new kind of sample for this technique. Previously investigated polymeric samples exhibited typically two glass transition temperatures, each corresponding to one phase, and the difference between the Tg values ranged from 80 K to 175 K. In all cases, the experimental results have been obtained at an intermediate temperature between the two Tgs. This implied a high difference of mobility between the two phases at the measurement temperature, i.e. a high mobility contrast and thus a significant difference in residual dipolar couplings. In those cases, the more mobile phase exhibited a line narrower than the one of the less mobile one by at least an order of magnitude. In the present work, the same technique was applied to homopolymers in the melt. All samples exhibited a single Tg and a low mobility contrast, as concluded from 1 H and 2D-WISE spectra. In homopolymer melts close to Tg, dynamic heterogeneities have been detected and their sizes have been determined by combining conventional spin diffusion with reduced 4D exchange experiments. [16] The question then arises, whether similar experiments can be performed at somewhat higher temperatures, applying the dipolar filter as a source of selection. Moreover, care has to be taken to distinguish magnetization exchange from cross-relaxation via the nuclear Overhauser effect (NOE), [START_REF] Ernst | Principles of nuclear magnetic resonance in one and two dimensions[END_REF] which occurs on a molecular level. Such experiments are described below emphasizing the difficulties of distinguishing these two cases from incomplete data sets and describing ways to avoid them.

Materials and Methods

Samples description

The synthesis of the poly(ethyl methacrylate) samples by conventional free-radical polymerization in solution was described elsewhere. [START_REF] Kulik | [END_REF] The samples investigated in the present work were, respectively, non isotopically labeled (PEMA), fully deuterated on the main chain (PEMADMC), and 20 % 13 C statically labeled on the carbonyl group (PEMA13C). Chemical structures are shown on Figure 1 and the results from characterization given in Table 1. The glass transition temperatures, Tg, were determined using differential scanning calorimetry (DSC) at 10 Kmin -1 and molar masses via size exclusion chromatography (SEC) in THF, calibrated with PMMA standards and the Mark-Houwink-Sakurada parameters (K=9.4410 -5 dLg -1 and =0.719 for PMMA, K=9.7010 -5 dLg -1 and =0.714 for PEMA) [19].

NMR experiments

Solid-state NMR measurements were carried out on a Bruker DSX spectrometer at Larmor frequencies of 300.13 MHz for 1 H and 75.47 MHz for 13 C using commercial 7.5 mm static and 4 mm magic angle spinning (MAS) double resonance probes from Bruker BioSpin GmbH. All spectra were externally referenced through adamantane ( 1 H chemical shift at 4.64 ppm, 13 C CH chemical shift at 38.5 ppm). For the measurements done under static conditions, the actual sample temperature was calibrated using lead nitrate [20] and a few melting points (due to the lack of proper standard to calibrate the 207 Pb ppm scale [21]). [22] 1 H spectra were recorded using single pulse excitation under static conditions in a 7 mm rotor, with 62.5 kHz RF nutation frequency corresponding to a 4 µs 90° pulse length, 16 transients and 5 s recycle delay. Temperatures ranged between Tg-45 K and Tg+115 K for PEMA, Tg-45 K and Tg+129 K for PEMADMC, Tg-46 K and Tg+128 K for PEMA13C.

A 2D-WISE spectrum [23] was recorded under static conditions on PEMA13C at Tg+81 K using a 7 mm rotor. A 5 µs 90° 1 H pulse was used, followed by a 180° pulse in the middle of indirect time to refocus the chemical shifts, in order to acquire the static 1 H line width in the indirect dimension. In order to avoid 1 H spin diffusion, polarization transfer from 1 H to 13 C was accomplished by Lee-Goldburg cross-polarization, with 500 µs contact time. During 13 C acquisition, 1 H continuous wave decoupling at 50 kHz RF nutation frequency has been applied. 160 increments were acquired in the indirect 1 H dimension and 256 transients were recorded in the direct 13 C dimension with 2 s relaxation delay between consecutive transients.

Experiments monitoring the decay of magnetization in the more mobile parts of the sample after a dipolar filter [START_REF] Egger | [END_REF] were carried out under static conditions (s. Figure 2 for pulse scheme). 62.5 kHz RF nutation frequency and 5 s recycle delay were used. The delay  between pulses in the dipolar filter was varied from 10 to 20 µs, and the number n of cycles in the dipolar filter from 1 to 12, depending on the sample and the temperature. The measurements were carried out at temperatures ranging from Tg+55 K to Tg+115 K on sample PEMA, from Tg+60 K to Tg+100 K on sample PEMADMC.

Lee-Goldburg CP MAS spectra were recorded at 3 kHz MAS spinning frequency using 4 mm rotors. In the 1 H channel the RF power level was set to a nutation frequency of 83 kHz (corresponding to a 90° pulse length of 3µs) for excitation as well as for 1 H continuous wave decoupling during acquisition. A relaxation delay of 3 s was chosen for these experiments.

The Lee-Goldburg irradiation [24] for the CP was adjusted on the 1 H nuclei by first calculating the corresponding irradiation offset frequency for the chosen RF power level and then finely adjusting the irradiation power by optimizing the resolution of 13 C multiplets under Lee-Goldburg decoupling conditions. The following series of experiments was conducted (s. Figure 3): first a simple LG-CP spectrum, second a LG-CP spectrum recorded immediately after a dipolar filter, third a LG-CP spectrum recorded after a dipolar filter with the same experimental settings and a subsequent mixing time. The experiments were conducted on PEMA at ca 390 K (ca Tg+45 K) with a CP contact time of 500 µs, a filter with 20 µs delay and 1 cycle, mixing times of 1 ms and 50 ms, and 1 536, 25 600, 32 768, 46 080 transients for the four spectra, respectively. The experiments were conducted on PEMADMC at ca 390 K (ca Tg+45 K) with a CP contact time of 500 µs, a filter with 10 µs delay and 1 cycle, a mixing time of 8 ms, and respectively 4 096, 8 192, 15 360 transients for the three spectra.

Results and Discussion

Dynamic contrast in PEMA

The dynamic contrast, the difference in mobility between the more mobile and the less mobile parts of a sample, can be characterized by 1 H static spectra and two dimensional wideline separation (2D-WISE) experiments. A few representative 1 H spectra are shown for sample PEMA in Figure 4. At none of these temperatures, a simple superposition of a broad and a narrow line is observed. Instead, the line becomes narrower in a visually homogeneous way with increasing temperature. This means that the whole sample becomes more mobile with increasing temperature, and exhibits no significant dynamic contrast. However, weak dynamic contrast within the sample might still be present and probed using the dipolar filter.

In that case, the static spectra would be a superposition of lines with similar line widths, so that it would be difficult to visually differentiate them. The same experimental results were obtained for samples PEMA13C and PEMADMC (data not shown [22]).

The 2D-WISE technique [23] was used in order to characterize more precisely the dynamic contrast in molten poly(ethyl methacrylate) with a higher structural resolution. In a 2D-WISE spectrum the different sites of the molecule are resolved according to their chemical shifts in the 13 C (direct) dimension and are correlated with the line width in the 1 H (indirect) dimension. This experimental procedure provides information on the mobility of the corresponding site: the narrower the line, the more mobile the chemical site. 1D 1 H spectra extracted from the 2D-WISE spectrum of sample PEMA13C at Tg+81 K are presented on Figure 5. The lines exhibited the following order for decreasing mobility: side chain CH3, then main chain CH3, then CH2 groups. It should be pointed out that apparent differences in mobility are obscured by possible 1 H nuclear spin diffusion during the CP contact time. [23] Therefore Lee-Goldburg CP, providing 1 H homonuclear decoupling during the CP contact time and thus suppressing spin diffusion, was used here. [24][25] However, an accurate adjustment of the Lee-Goldburg conditions is only possible under MAS, due to the presence of heteronuclear dipolar couplings under static conditions. The 1 H nuclear spin diffusion would thus not be completely suppressed but only weakened in the 2D-WISE recorded here under static conditions. It should be noted that main chain and side chain CH2 groups cannot be differentiated in this experiment because of the poor spectral resolution in the 13 C dimension.

Combination of the dipolar filter with a conventional exchange experiment applied to investigate poly(ethyl methacrylate) Poly(ethyl methacrylate) melts were investigated by a combination of the dipolar filter [START_REF] Egger | [END_REF] and a conventional exchange experiment (s. Figure 2 for pulse scheme), monitoring the decay of 1 H magnetization of the more mobile parts after the dipolar filter. This investigation of samples with low mobility contrast required several changes in the data processing and analysis, which are presented below, after a short description of the usual procedure for a spin diffusion experiment.

In the conventional 1 H nuclear spin diffusion experiment, several FIDs are recorded as a function of mixing time m and Fourier transformed. For increasing m, a broad line of increasing area appears below the narrow line present at the end of the dipolar filter, indicating that magnetization has diffused to a less mobile component. Furthermore, the total intensity (line area) decreases, due to longitudinal relaxation T1. In the usual case of a high mobility contrast, the decay of the intensity of the more mobile 1 H magnetization is monitored as the area of the narrow line by integration of the signal in the corresponding spectral window. [3] The monitored decay is due to two factors: magnetization transfer to less mobile parts and T1 relaxation. In order to distinguish magnetization transfer and T1 relaxation, the magnetization is alternatively stored along +z and -z during the mixing time m, and the corresponding transients are subtracted. [START_REF] Schmidt-Rohr | Multidimensional solid-state NMR and polymers[END_REF] To correct the data for T1 relaxation, the recorded intensity of the more mobile parts is divided by the intensity obtained for the same mixing time without application of the dipolar filter. [START_REF] Clauss | [END_REF][3] The so-corrected intensity is then normalized by the initial intensity (extrapolated at m=0 from the linear initial slope as a function of the square root of the mixing time). The final curves are plotted as a function of the square root of the mixing time; they exhibit a linear behavior for short mixing times and a plateau for the long ones. According to the Babinet's principle, the initial decay does not exhibit any difference between the case of less mobile domains in a more mobile matrix and the inverse case. [START_REF] Schmidt-Rohr | Multidimensional solid-state NMR and polymers[END_REF] The plateau value corresponds to the percentage of selected magnetization after application of the dipolar filter (it usually slightly depends on the filter parameters due to the presence of an interphase with mobility gradient rather than a sharp interface between the more and the less mobile parts [9]). From the intercept of the initial slope with the X-axis, * m  , the domain size dsize of the underlying structure (lamella thickness, cylinder or sphere diameter) can be calculated, according to equation 1, where ε is the dimensionality (1 for lamellae, 2 for cylinders and 3 for spheres), and Deff the effective 1 H nuclear spin diffusion

coefficient through coherent magnetization transfer. * 2 m eff size D d       = (1)
Deff is calculated from the 1 H nuclear spin diffusion coefficients of the more mobile and less mobile phases Dmob and Drig according to equation 2.

        + = rig mob eff D D D 1 1 2 1 1 (2)
Drig is usually taken equal to the value measured for polystyrene below its Tg (0.8 nm 2 ms -1 ), while Dmob can be determined using its correlation with the transversal relaxation time T2. [3] The T2 relaxation time of the more mobile phase has to be measured independently, either using the Carr-Purcell-Meiboom-Gill (CPMG) experiment [26][27], or by measuring the line width of the static 1 H spectrum assuming a Lorentzian line shape, in either case just after applying the dipolar filter.

Determining the decay of the intensity of the more mobile 1 H magnetization by integration of the signal in a narrow spectral window is not applicable to poly(ethyl methacrylate) melts, because the contrast in mobility (and thus the line width difference) is too low. Therefore, the maximum height of the total 1 H spectrum was monitored as a function of the mixing time, which is more sensitive to the narrower component. This estimate allowed us to monitor the time dependence of the magnetization transfer. The absolute values of the degree of exchange, e.g. the plateau value, should, however, be interpreted with care.

Concerning the correction for T1 relaxation, it should be pointed out that the way explained above is valid only if the T1 value is homogeneous within the sample, or for mixing times much shorter than the shorter T1 values. [START_REF] Schmidt-Rohr | Multidimensional solid-state NMR and polymers[END_REF][28] T1 relaxation was measured in the samples using the inversion recovery technique [29][30] (data not shown [22]). It exhibited a monoexponential decay, characteristic of a spatially homogeneous T1 relaxation, or of samples with a relaxation sink and extensive 1 H nuclear spin diffusion. Furthermore, all T1 times were longer than 800 ms. For all samples, the diffusion curves were processed (initial linear slope and beginning of the plateau) for m values smaller than 100 ms, thus much smaller that the determined T1 relaxation values. Therefore the corresponding parts of the 1 H nuclear spin diffusion curves could be corrected for T1 relaxation via division by the corresponding signal without dipolar filter, as detailed above. In poly(ethyl methacrylate) melts, the whole sample was much more mobile than polystyrene below its Tg, so that the less mobile phase could not be assumed to have a spin diffusion coefficient of 0.8 nm 2 s -1 . Experimentally, only a single value of spin diffusion coefficient could be determined. Since both phases have similar mobility and 1 H spin density, they have similar spin diffusion coefficients. Therefore the average of the spin diffusion coefficient is close to the diffusion coefficient of the more mobile phase. Furthermore, due to the way of averaging, the diffusion coefficient of the more mobile phase has a stronger weight on the average than the other one, since it is lower.

Finally, only the diffusion coefficient of the more mobile phase was determined, and assumed to be close to the one of the less mobile phase, thus similar to the average diffusion coefficient. The diffusion coefficient of the more mobile phase was determined via the T2 relaxation time. Both methods (CPMG experiment, line width) were used for all samples.

Only the value calculated from the line width measurement was considered for two reasons.

First, the line is broader at the temperature of the 1 H nuclear spin diffusion measurements than at higher temperatures (indicating an intermediate mobility). Second, the diffusion coefficients determined through CPMG were higher than the ones determined from the line width of the 1 H static spectrum (after a dipolar filter), and both methods overestimate the diffusion coefficient; therefore the lower value was assumed to be more accurate. There was no information on the geometry of the investigated system, since X-ray diffraction detects a smaller structure and the contrast is too low in homopolymers for electron microscopy. In order to minimize the uncertainty, the dimensionality  was assumed to be 2.

In sample PEMA at Tg+67 K, 1 H polarization transfer occurred after the dipolar filter, as indicated by the gradual broadening of the basis of the line with increasing mixing time. The evolution of the 1 H magnetization of mobile species with the mixing time after the dipolar filter is shown on Figure 6a. This plot exhibits the typical evolution of the normalized and corrected intensity as a function of m  in a 1 H nuclear spin diffusion experiment using the dipolar filter: a signal decay linear in the square root of the mixing time for small mixing times and a plateau for long mixing times. Processing the data as detailed above for the specific case of polymer melts yielded apparent domain sizes of 3 to 6 nm. Furthermore, the plateau values ranged from 53 to 70 %. Assuming that the dynamic heterogeneity detected by the dipolar filter is related to nanostructuring in PEMA at Tg+67 K, the quantified size would be in accordance with the typical length of 5 to 10 monomeric units determined in NMR [15] and X-ray scattering [14] studies. However, the actual selection of polarization due to the dipolar filter needed to be checked.

Actual selection done by the dipolar filter and actual subsequent magnetization transfer mechanism

In order to check the actual selection resulting from the application of the dipolar filter, the selected 1 H magnetization was transferred to 13 C nuclei and acquired in the 13 C channel, to gain chemical shift resolution. Obscuring spin diffusion during the magnetization transfer from 1 H to 13 C nuclei was reduced using Lee-Goldburg cross-polarization (LG-CP, s. above). It would have been best to carry out these experiments under the same conditions as the experiments using the dipolar filter described in the preceding paragraph, namely under static conditions and at Tg+67 K. However, due to the poor chemical shift resolution and the difficulties adjusting Lee-Goldburg irradiation on 1 H spins under static conditions, the investigations were conducted under MAS. Due to a limited temperature range of the available MAS probehead, the investigations were conducted at ca 390 K, i.e. ca Tg+45 K.

The pulse schemes used to investigate the selection done by the dipolar filter using Lee-Goldburg CP are presented in Figure 3 roughly corresponds to the plateau values of 0.5 to 0.7 observed on Figure 6). The 13 C LG-CP spectrum (c') was identical to the one shown in (a'), proving that the magnetization was back at equilibrium 50 ms after the application of the dipolar filter.

In the preceding paragraph, it was assumed for the data processing that the dipolar filter would select mobility contrast due to nanodomains present in PEMA. This would have led to the selection of some whole monomeric units. Therefore it would have resulted in the presence of all the chemical parts of the monomeric units in the 13 C LG-CP spectrum recorded after the dipolar filter and no mixing time. However, the LG-CP investigations conducted on PEMA at Tg+45 K clearly showed that the dipolar filter actually selected essentially the CH3 group of the side chain. This proved that the mobility contrast in PEMA at ca Tg+45 K observed using the dipolar filter did not result from structured and unstructured domains on the nanometer length scale. Instead, the dipolar filter selected the end group of the alkyl side chain only. Similar experiments conducted at ca 390 K (ca Tg+45 K) on sample PEMADMC proved also the selection of the CH3 group of the side chain only (data not shown [22]).

It was assumed above that the magnetization transfer would have occurred via coherent, energy conserving flip-flops, as it is the case in a typical 1 ( )

2 0 5 . 0 1 m z DC M M    -  = (3)
The quadratic behavior was experimentally observed for mixing times up to 400 s. In the present investigation of PEMA at Tg+67 K under static conditions, a coherent transfer of magnetization can be expected for even shorter mixing times, due to the lower molecular mobility of PEMA. However, the time of 400 s would have been between our second and third experimental point. Furthermore, the magnetization decay of interest in the present study was on the order of several ms (s. Figure 6a), thus of a factor 10 to 100 longer than the time range of coherent transfer. Indeed, plotting our experimental magnetization decays as a function of the square of the mixing time did not result in a linear behavior for very short mixing times (s. Figure 6b). Therefore, it was concluded that in PEMA at Tg+67 K the magnetization transfer after the dipolar filter occurred predominantly via incoherent zero-and double-quantum transitions. Thus the magnetization transfer data obtained in the present work should be processed considering incoherent zero-quantum and double-quantum transitions, i.e. a NOE mechanism. This transfer mechanism results in an exponential decay of the magnetization, as detailed in the following paragraph. Indeed, the linear behavior observed for the recorded magnetization decay after subtraction of the plateau value, normalization and plot on a logarithmic scale (s. Figure 6c) is in agreement with an exponential decay. It should be noted that the magnetization decay caused by ideal zero-quantum transitions in a dense spin ensemble turning into spin diffusion, as calculated from theoretical equations [START_REF] Schmidt-Rohr | Multidimensional solid-state NMR and polymers[END_REF], exhibits a quasi-linear behavior when plotted on a logarithmic scale as a function of the mixing time (s. Figure 7c). Therefore, the knowledge about the time dependence of the magnetization decay is not sufficient to discriminate the mechanism of the polarization transfer. However, in the case of PEMA at Tg+67 K, an incoherent cross-relaxation mechanism for the polarization transfer along alkyl side-chain appears to be more likely than the multi-step process, but some coherent contributions to the predominant cross-relaxation transfer in PEMA can not be excluded.

Equation governing the monitored magnetization decay via cross-relaxation

The magnetization transfer after application of the dipolar filter in PEMA at Tg+67 K occurred predominantly via incoherent zero-quantum and double-quantum transitions, i. all diagonal lines have magnetization. However, "having magnetization" only means that for this particular species, there is excess of spin orientation in one direction of the order of parts per millions. Therefore, it doesn't make a difference for the cross-relaxation at the level of the individual spins, if only one species A has initially magnetization or both species A and B.

(see [START_REF] Abragam | The Principles of Nuclear Magnetism[END_REF]). As a result, the mathematical equations for the time evolution of the diagonal lines intensity in a 2D NOE study [START_REF] Macura | [END_REF] directly apply to the present one-dimensional monitoring of the more mobile sites magnetization.

The most elaborate analytical expression describing the evolution of the diagonal line intensity aAA with mixing time m is given by Macura and Ernst [START_REF] Macura | [END_REF] for the case of a homonuclear system composed of a group of nA equivalent A nuclei and a group of nB equivalent B nuclei. In the investigated samples, the initial magnetization is mainly located at the CH3 end group of the alkyl side chain. Therefore a CH3-CH2 moiety was considered as source of cross-relaxation. It should be emphasized that this is strictly valid for sample PEMADMC, but does not correspond to the whole monomeric unit of sample PEMA; however, to our knowledge, no analytical equation is available for moieties larger than two groups of equivalent nuclei. Moreover, since the initial magnetization is mainly located at the CH3 end group, the CH3-CH2 contribution is expected to be predominant in the initial magnetization decay for sample PEMA. The investigated samples are macromolecules in bulk, therefore they are in the slow motion limit; moreover, they are homonuclear systems.

Macura and Ernst [START_REF] Macura | [END_REF] did not give explicitly the equation corresponding to that case if nA is different from nB; however it could easily be derived from equation 4 [START_REF] Macura | [END_REF], where M0 is the total equilibrium magnetization of the nA+nB nuclei, and the cross-relaxation rates Rij are defined below.

( ) ( ) ( ) ( )       -         - + +         - -  - + = m C C BB AA C BB AA m L B A A m AA R R R R R R R R n n M n a    exp 1 1 exp 2 0 with ( ) BA AB BB AA C R R R R R 4 2 + - = and
( )

C BB AA L R R R R 2 1 2 1 - + = (4)
The cross-relaxation rates Rij can be expressed by the transition probabilities W, resulting from AA, AB and BB interactions, and by the external relaxation rates R1A and R1B, which take into account possible interactions with further spins (s. equation 5 [START_REF] Macura | [END_REF]). In this equation, the subscripts 0, 1(A) and 2 relate respectively to zero-quantum transition, single-quantum transition (with flip on A nucleus), double-quantum transition; the superscripts ij relate to both nuclei involved in the corresponding transition.

( ) ( ) ( ) ( )

AB AB A AB A AB AB A AB B AA AA A AA W W n R R W W W n W W n R 0 2 1 2 1 0 2 1 2 1 2 - = + + + + +  - = (5)
The transitions probabilities W obey equation 6 [START_REF] Macura | [END_REF], with the spectral densities Jxy() and the constants qxy defined in equation 7 [START_REF] Macura | [END_REF]. In the latter, C xy is the correlation time of the isotropic motion which modulates the xy interaction, 0 is the permeability of space, x is the magnetogyric ratio of spin x, ħ is the reduced Planck's constant, and rxy is the internuclear xy distance. Please note that the equation given for the spectral density is obtained for a rigid molecule undergoing isotropic random motions. [START_REF] Likic | [END_REF] Moreover, it should pointed out that the investigated phenomenon is local due to its dependence on the inverse sixth power of the internuclear distance.

J q W J q W J q W J q W       + = - = = = 6
            = (7)
The investigated samples are macromolecules in bulk, therefore they are in the slow motion limit and XC >> 1. The investigated system is a homonuclear system, so that (X-Y)C  0, therefore W1 XY =W2 XY =0, and the zero-quantum transitions dominate crossrelaxation, with a transition probability of W0 AB =qABC AB . The elements of the crossrelaxation matrix can thus be reduced to RAA=nBW0 AB +R1A, RBB=nAW0 AB +R1B, RAB=-nAW0 AB and RBA=-nBW0 AB . Assuming equal external relaxation rates R1 for A and B nuclei, the quantities RC and RL can then be calculated as RC=(nB+nA)W0 AB and RL=R1, then (RAA-RBB)/RC=(nB-nA)/(nA+nB).

With this, the intensity of a diagonal signal is described by equation 8. ( )

( ) ( ) ( ) ( ) ( )        + -         + - - +         + - +  - + = m AB C AB B A B A B A B A B A m B A A m AA q n n n n n n n n n n R n n M n a    
( )       - +  = m AB C AB m AA q M a    5 exp 5 4 5 6 10 3 0 (9)
Thus, the processing of the data aAA(m) recorded after application of the dipolar filter allows to extract information on the molecular dynamics, via the determination of the correlation time of the involved molecular motion C AB , if the parameter qAB is known.

Several approaches can be used to estimate the coupling strength qAB, i.e. from equation 7 assuming the internuclear distance or considering the second moment for different geometries. As detailed in the Appendix, these values differ but are of the same order of magnitude. Considering the overall uncertainty of our approach, qAB= 9 kHz 2 for PEMA and 6.3 kHz 2 for PEMADMC as deduced from the second moment were used to extract correlation times from the relaxation rates. [37] Distinguishing different mechanisms of magnetization exchange in dense spin systems

The difficulty to distinguish the different mechanisms of magnetization exchange from the mixing time behavior in dense spin systems is demonstrated in Figure 7. Here the evolution of the 1 H magnetization is modelled from the equations given in. [START_REF] Schmidt-Rohr | Multidimensional solid-state NMR and polymers[END_REF] The plots demonstrate that magnetization exchange by "spin diffusion" also leads to quasi-exponential behavior at long exchange times. Thus, a definite conclusion how the magnetization exchange curves should be interpreted is not possible from such data alone. In polymers above Tg, however, considerable information about the relaxation times is usually available from other sources, such as dielectric relaxation, light scattering or other NMR methods. [START_REF] Schmidt-Rohr | Multidimensional solid-state NMR and polymers[END_REF] In the following we will therefore examine plots like Figure 7 for PEMA at different temperatures and relate them to literature data. [15] Quantification of local dynamics inside side chains of poly(ethyl methacrylate)

The NOE experiment with dipolar filter was carried out at temperatures ranging from Tg+55 K to Tg+115 K for PEMA, and from Tg+60 K to Tg+100 K for PEMADMC. For each temperature, various parameters were used for the dipolar filter.

At temperatures higher than Tg+85 K for sample PEMADMC and higher than Tg+60 K for sample PEMA, a mono-exponential decay was observed for the recorded magnetization.

The data were fitted as a single cross-relaxation process, which follows equation 9: an exponential decay followed by a plateau. In order to fit experimental data, programs were written using the Matlab ® (The MathWorks) software, which handled simultaneously a series of measurements done at one temperature using different parameters for the dipolar filter. The plateau value was first determined as the average of all experimental values measured for mixing times longer than a chosen threshold, and subtracted from all experimental values in order to obtain an exponential decay to zero. Then the experimental data were normalized so that the exponential decay begins at 1 for m=0, and fitted with a mono-exponential function (s. equation 10) to extract values for the product qABC AB .

( ) ( )

m AB C AB m q f    5 exp - = (10)
In the case of PEMA at Tg+67 K, the fits yielded as series of qABC AB values characterized by an average of 275.2 Hz, a range of 15.4 Hz and a standard deviation of 5.6 Hz for the mono-exponential decay. With qAB= 9 kHz 2 determined from second moment at a temperature well below the glass transition temperature (s. above), a value of 310 -5 s was obtained for the correlation time C AB . The shape of the curves for high temperatures (above Tg+50 K for PEMA and above Tg+85 K for PEMADMC) were identical to those obtained at Tg+67 K for PEMA, only the numerical values varied. The example of sample PEMA at Tg+115K is given on Figure 8: a linear behavior was observed for the logarithm of the recorded intensity as a function of the mixing time over 15 ms. The obtained correlation times C AB of the dynamic process occurring in the alkyl side chains and giving rise to the cross-relaxation are shown Figure 9 for samples PEMA and PEMADMC.

At lower temperatures, a more complex decay was observed for the recorded magnetization, as shown in Figure 10 for PEMA at Tg+55 K. This behavior could be interpreted as the superposition of two distinct NOE processes with different correlation times, resulting in a bi-exponential decay (s. Figure 10c). It could also be interpreted as a coherent spin diffusion process, as shown by the similarity of the curves shapes in Figures 7c and Figure 10c. The existing data did not allow us to decide between these two interpretations, and the magnetization transfer occurring in the sample was probably a superposition of both processes. The data were fitted using the Microcal TM Origin ® software and by a bi-exponential decay. The slow decay was fitted first as a linear decay of the logarithm of the magnetization versus mixing time. Then the fast decay was fitted as one component of a bi-exponential decay where the other component was set as the slow decay determined previously. Note that the slower magnetization decay corresponded to the faster molecular motion, and thus is designated as "fast" in the following. The extracted correlation times are plotted on Figure 9. Please note that the slow process was quantified, while the fast one could be estimated only. The slow process, corresponding to the fast decay, exhibited the same Arrhenius behavior as the single NOE process quantified at higher temperatures (s.

Figure 9). Therefore, the complex behavior at lower temperatures was interpreted as the superposition of two independent NOE processes with different correlation times.

Conclusion

The dipolar filter, combined with a conventional exchange experiment, was applied to homopolymer melts of different poly(ethyl methacrylate) samples. Recent x-ray and NMR studies of these melts indicated the presence of nano-structures in the melt state. However, the melts exhibited only low dynamic contrast originating from the difference in mobility between distinct parts of the monomeric units. In the experiments, the dipolar filter selected only the end of the alkyl side chain, which has been checked using Lee-Goldburg CP experiments. The subsequent 1 H polarization transfer was proved to occur via non coherent zero-quantum transition. The NOE processing of our experimental results allowed us to quantify correlation times of dynamics processes in the alkyl side chains. At the lower temperatures, two distinct processes were observed. The determined correlation times are in good agreement with values for the complex chain motion of this polymer known from literature. A sample related interpretation of the detected processes will be presented in a future article, together with results obtained for other poly(n-alkyl methacrylates) and poly(nalkyl acrylates) melts at various temperatures.

The dipolar filter is used routinely in polymeric samples to quantify domain sizes in phase separated samples. The present work points out processing problems that one should be aware of, if systems with poor phase separation or local structure on a sub nanometer scale are studied. On the other hand, the introduced NOE processing procedure opens a new experimental route to estimate correlation times of local dynamic processes in non isotopically labeled systems using solid state NMR.

internuclear distance rAB, which has to be estimated by e.g. molecular modeling. Second, it can be calculated as a value proportional to the second moment [START_REF] Slichter | Principles of Magnetic Resonance[END_REF] of the 1 H line recorded well below Tg under static conditions. The calculation of qAB was done as follows. In the homonuclear 1 H-1 H case, the qAB factor was calculated for an internuclear distance of 1 Å=10 - 10 m from the magnetogyric ratio of 1 H spins (26.7510 -7 rads -1 T -1 ) [START_REF] Harris | [END_REF] and the physical constants 0/4 and ħ (10 -7 NA -1 and 6.62610 -34 /2 Js -1 rad -1 respectively) [START_REF] Homann | Abbreviated list of quantities, units and symbols in physical chemistry[END_REF]. This yielded a value of 56.910 9 rad 2 s -2 i.e. qAB=1.44 s -2 for 1 H-1 H and rAB= 1 Å. The actual rAB distance between a 1 H nucleus of the CH3 and a 1 H nucleus of the CH2 group was calculated as a function of the involved dihedral angle, then the rHH -3 values (proportional to the dipolar coupling) were averaged for equi-distributed values of the dihedral angle between 0 and 360°, which yielded an average distance of rHH=2.67 Å. [22] The calculated qAB parameter was then qAB=3.98 kHz 2 for a CH3-CH2 moiety.

Independently, the qAB parameter was measured through the second moment of the 1 H line recorded below Tg under static conditions. The equation describing the qAB factor (s.

Starting from expressions of M2 in the MKSA unit system [42][43] [44], assuming that all nuclei are equivalent, [START_REF] Abragam | The Principles of Nuclear Magnetism[END_REF] taking the powder average over all orientations, [START_REF] Abragam | The Principles of Nuclear Magnetism[END_REF] and taking into account the value I=1/2 for 1 H nuclei, the equation 12 was obtained for M2, where rij is the internuclear distance between the observed spin and any other one. [22] 

        = j ij r M 6 2 4 2 0 2 1 4 20 9     (12)
The residual sum was simplified considering either an isolated spin pair [45][46] or a simple cubic lattice. [START_REF] Abragam | The Principles of Nuclear Magnetism[END_REF] The resulting expressions are shown on equation 13. [22] ( )

Tables and Figures

O O C CH 2 n CH 2 H 3 C CH 3 (a) O O C CD 2 n CH 2 H 3 C CD 3 (b) O O 13 C CH 2 n CH 2 H 3 C CH 3 (c)

  (left). In the first experiment (a), a simple LG-CP spectrum is recorded to obtain a reference spectrum. In the second experiment (b), a dipolar filter is applied directly followed by LG-CP, in order to determine the parts of the sample actually selected by the dipolar filter. The third experiment (c) corresponds to (b), where a mixing time m is introduced between the dipolar filter and the acquisition of the LG-CP spectrum, in order to observe the relaxation of the sample back to equilibrium. The corresponding spectra for sample PEMA at ca Tg+45 K are shown on Figure3(right). The carbonyl signal was too weak to be detected. The 13 C LG-CP spectrum (a') exhibited a chemical shift resolution high enough to resolve all the chemical sites of the monomeric unit of PEMA. Furthermore, it gave their reference intensities in a LG-CP spectrum. It is obvious from the 13 C LG-CP spectrum (b') that the dipolar filter actually selected essentially the CH3 group of the side chain of PEMA. Small amounts of other groups were also selected. (It should be noted that this selection of 3 out of 10 1 H nuclei of the monomeric unit of PEMA

  H nuclear spin diffusion experiment. The LG-CP experiments showed, however, that the dipolar filter selected essentially the end group of the alkyl side chain and that the equilibrating polarization transfer took place from the end group of the alkyl side chain along the alkyl side chain and further to the main chain. Such a magnetization transfer along an alkyl side chain can occur via either coherent or incoherent transfer. In the case of coherent transfer, the residual dipolar couplings would cause zero-quantum transitions, which in the limit of many flip-flop transitions would result in incoherent 1 H nuclear spin diffusion. Incoherent transfer, occurring via incoherent zero-quantum or double-quantum transitions as cross-relaxation caused by fluctuations of the dipolar coupling due to the chain motion, is often termed Nuclear Overhauser Effect (NOE).Since the physical mechanisms behind coherent or incoherent magnetization transfer are completely different, resulting in different time behaviors,[3][31] a quantitative analysis of the experimental data required a detailed knowledge about the dynamic processes in the samples.Fritzhanns et al.[32] investigated the magnetization transfer mechanism in multidimensional NOE experiments carried out on elastomer samples under MAS; they indicated coherent transfer in the case of static measurements, and incoherent transfer in the case of MAS measurements. Demco et al.[33] conducted a detailed study of SBR elastomers at Tg+70 K using NOE experiments under static conditions, and concluded a coherent transfer mechanism in the limit of short mixing times. This leads to a quadratic decay of the magnetization Mz with the mixing time m, described by equation 3, where M0 was the initial magnetization and <DC> the average residual dipolar coupling.

  e. a cross-relaxation or NOE mechanism. Different initial conditions and different transfer conditions are possible for NOE experiments, leading to different equations governing the magnetization transfer.[START_REF] Ernst | Principles of nuclear magnetic resonance in one and two dimensions[END_REF] In the present work, the line height was monitored, which is equivalent to the integral over a thin slice of the spectrum, or to integrate predominantly the mobile component of a superposition of a more mobile component and a less mobile component. Furthermore, in the present experiments, magnetization was present initially only at the more mobile sites, and was transferred in the course of our experiment to all sites depolarized by the dipolar filter. Thus it resembles a 2D NOE experiment where two components A and B would be resolved on the chemical shift scale, where the magnetization would have been selected at A sites initially, where the magnetization would then have been transferred between the sites A and B via cross-relaxation, but only the magnetization of component A would have been recorded. Finally, the present way of monitoring the magnetization decay is equivalent to the integration over time of a diagonal line (A) in a 2D NOE experiment. It should be noted that the present initial conditions (where only the monitored line had magnetization) differ from the classical 2D NOE initial conditions, where

  of a CH3-CH2 moiety, nA=3 and nB=2. Taking the experimental correction of the longitudinal relaxation compensating the exponential factor simplified further to equation 9.

  The obtained correlation time C AB of the local motion occurring in the alkyl side chains and giving rise to the cross-relaxation are shown Figure9for samples PEMA and PEMADMC. It was observed that in the temperature range studied the dependence of C AB on inverse temperature was nearly linear. Furthermore, the values determined for sample PEMADMC were in agreement with those determined for sample PEMA, validating the assumption of the predominance in the initial magnetization decay of the cross-relaxation from the CH3 to the CH2 group along the alkyl side chain. A linear regression of the C AB values yielded a prefactor of 3.10 -13 s and an activation energy of 27.5 kJmol -1 . The obtained correlation times are in good agreement with values from literature. The slow and the fast component correspond to chain isotropization and axial chain motion, respectively. Details, as well as complementary results obtained on other poly(n-alkyl methacrylates) and poly(n-alkyl acrylates) will be presented in future publications.

equation 7 ,

 7 in MKSA unit system) is very similar to the one describing the second moment M2 of the rigid lattice (s. equation 11,[START_REF] Van Vleck | [END_REF][START_REF] Abragam | The Principles of Nuclear Magnetism[END_REF] in CGS unit system, in the case of a powder of crystalline sample exhibiting a cubic lattice, where I is the spin quantum number, and d the lattice constant).

Figure 1 :

 1 Figure 1:

Figure 1 :

 1 Figure 1: Chemical structure of poly(ethyl methacrylate) samples; (a) PEMA, (b) PEMADMC, (c) PEMA13C. Selective isotopic enrichment is indicated by chemical symbols given in bold fonts.

Figure 2 Figure 2 :

 22 Figure 2:

Figure 3 :Figure 3 :

 33 Figure 3:

Figure 4 :

 4 Figure 4:

Figure 4 :

 4 Figure 4: Influence of the temperature on the shape of the 1 H spectrum of sample PEMA (recorded under static conditions); the dotted frame indicates the temperature at which the investigation using the dipolar filter was carried out first. Note that the very narrow line below Tg most probably originates from small molecules present in the sample and not from the sample itself.

Figure 6 :

 6 Figure 6:

Figure 7 :

 7 Figure 7:

Figure 8 :

 8 Figure 8:

1 ]Figure 9 :

 19 Figure 9:

Figure 10 Figure 10 :

 1010 Figure 10:

Table 1 :

 1

	Characterization of different poly(ethyl methacrylate) samples
		PEMA	PEMADMC	PEMA13C
	Tg [K]	342	345	338
	Mn [gmol -1 ]	115 000	77 700	55 40
	Mw [gmol -1 ]	156 200	107 600	122 200
	Mw / Mn	1.36	1.38	2.21

Table 1 :

 1 Glass transition temperatures from DSC and molar masses from SEC in THF of the investigated samples.

Acknowledgements

The authors gratefullly acknowledge the Arkema company (Total group), in particular Francois Beaume and Rene-Paul Eustache, for financial support. We thank Uta Pawelzik for samples synthesis (in particular the synthesis of selectively labeled monomers), Michael Wind and Mario Beiner for sample characterizations and for fruitful discussions, Manfred Hehn and Hans-Peter Raich for technical help with the NMR spectrometers, Manfred Wilhelm and Attila Domjan for discussions, Lothar Brombacher for providing the program allowing the simulation of the ideal spin diffusion curve.

Appendix

The qAB constant can be determined by two different approaches. First, it can be calculated according to equation 7, where all the terms are known factors, except the By simple comparison of equations 7 and 13, it was deduced that M2(spin pair)=4.5 qAB and M2(cubic lattice)=38.25 qAB. The factor of 8.5 found between the results of the two approaches originates mainly in the presence of 6 next neighbors (and some more remote ones) in the case of the cubic lattice versus a single one in the case of the spin pair. The present case was rigorously neither a spin pair nor a cubic lattice. The geometry of the investigated system was not sufficiently known to draw conclusions on the number and positions of the neighbors. However, several neighbors were present, so that the simple cubic lattice approximation was retained as more realistic. Thus the numerical value of qAB was obtained dividing the experimental value of M2 by a factor of 38.25.

The numerical integration of a 1 H spectrum recorded under static conditions at Tg-45 K (s. Figure 4 for sample PEMA, data not shown for sample PEMADMC [22]) yielded a second moment value of M2=342 kHz 2 for PEMA and M2=240 kHz 2 for PEMADMC; thus values of qAB=9 kHz 2 for PEMA and 6.3 kHz 2 for PEMADMC. Compared to the calculated qAB parameter for a CH3-CH2 moiety: qAB=4 kHz 2 , the values determined via the second moment were significantly different, but of the same order of magnitude. Considering the overall uncertainty of the present technique, the correlation times were calculated using the qAB parameter value obtained from the experimental second moment exclusively.


domain: phys spi

text:
Denis Gueyffier
email: denis.gueyffier@onera.fr

Sylvie Plot
email: sylvie.plot@onera.fr

Matthieu Soismier
email: matthieu.soismier@onera.fr

SoNICS: a New Generation CFD Software for Satisfying Industrial Users Needs

We briefly present a novel software called SoNICS (So Naturally Innovative CFD Solver) for high performance CFD simulations. SoNICS aims at replacing in the medium term, both at ONERA and at several industrial and academic partners, the elsA aerodynamics simulation software. This software has a functional scope close to elsA with some major additional features. SoNICS also leverages innovations which rely on recent parallel computing and software engineering paradigms and on automatic code generation and differentiation. Several drivers have emerged which have led ONERA to designing an entirely new architecture for the SoNICS software. While our partner Safran prepares for the launch of SoNICS at its design offices in the next few years, Safran also makes use of elsA for its current motor programs which generate new needs for the elsA software. Development and validation of both SoNICS and elsA are steered though a close cooperation, as well as a common governance between ONERA and Safran. Both solvers are mainly developed by the same teams within a project supported by DGAC, the French Civil Aviation Authority. For the development of both solvers, we have set up an agile organization targeted at satisfying industrial users' needs. This organization relies on three pillars: i-handling user' needs instead of delivering features, ii-having the end-user involved at all steps of the development process, iii-testing in the end users' environment. We finally briefly illustrate a few recent applications enabled with elsA thanks to this new organization and the first milestones achieved with the new SoNICS software.

INTRODUCTION

During the past three decades, ONERA has developed aerodynamic simulation software for the industry and it has put a lot of effort in aligning its CFD offer with industrial needs. During this time period the elsA software [START_REF] Cambier | The Onera elsA CFD software: input from research and feedback from industry[END_REF] (property of ONERA and Safran) for example has grown from an academic scale to an industrial scale compressible internal and external flow solver which is now currently in use at several ASD industries partners like Safran, Airbus and MBDA where hundreds of users at design offices make use of it on a daily basis. Hundreds of simulations with elsA run each day on a dozen of supercomputers and, more recently, on individual computers for smaller design calculations. elsA is nowadays the main aerodynamics simulation software in use at Safran for RANS simulations for turbomachinery design. Industrialization has been at the heart of ONERA CFD teams' work. These teams have gained experience over the years and learned how to ship software matching end-users needs in industrial software platforms. Solid production and industrialization teams with high quality standards have been put in place in order to integrate new features, deploy the software at industrial and academic partners and bring support to the users. These teams contain product managers who gather and steer industrial needs and application teams who test software and run uses cases with excellent level of technical expertise of the industrial field. elsA has also captured over time models originating from various research fields: numerical analysis [START_REF] Lerat | High order residual-based compact schemes for aerodynamics and aeroacoustics[END_REF] [3] [4] [START_REF] Bras | Development of Compressible Large-Eddy Simulations Combining High-Order Schemes and Wall Modeling[END_REF], turbulence and transition [START_REF] Aupoix | Turbulence and transition modeling[END_REF][7][8] [START_REF] Bégou | Database Approach for Laminar-Turbulent Transition Prediction: Navier-Stokes Compatible Reformulation[END_REF][10] [START_REF] Bouchard | Validation of correlations-based transition modeling strategies applied to the Spalart-Allmaras turbulence model for the computation of separation-induced transition[END_REF], CFD-CSM [START_REF] Dugeai | Aeroelastic capabilities of the elsA solver for rotating machines applications[END_REF][13], ... However, in the past few years, because of drivers ranging from competition with other CFD software providers addressing the same market to breakthrough innovations in the global CFD offer (e.g. LBM software) or in the hardware technology offer (e.g. GPU clusters), ONERA has questioned itself and has challenged at the same time both its technologies and its means to industrialize its CFD software. This has led ONERA on the one end to propose new ways to improve communication with industrial users to better respond to their needs and on the other end to challenge its technical offer to become faster at capturing breakthrough innovations in CFD software products. The effort goes into two directions: pulling needs from industrial users in order to capture in software products the features responding to these needs, and pushing in the same products scientific and technical innovations from its most innovative technical and scientific experts or academic partners. In the following, we describe why and how we started developing a new software with breakthrough innovations and building a new organization to respond to users' needs.

Safran has decided to partner with ONERA along this strategy and to make use of elsA for its current and near future aircraft motor designs, and then to switch its design offices in a few years to the next generation software SoNICS. Both elsA and SoNICS are property of ONERA and Safran. On the one hand, for its near ongoing motor programs, Safran needs additional features in elsA and on the other hand Safran needs to anticipate the changes, challenges and opportunities coming with SoNICS for its future programs.

ONERA thus follows two avenues simultaneously with the same development and application teams: continue development of the current elsA software and start building the future SoNICS CFD software. The main project supporting these developments is named SONICE. This collaborative project between Safran Aircraft Engines, Safran Helicopter Engines, Safran Tech and ONERA is supported by the French DGAC and is led by Safran Tech. Both companies, ONERA and Safran, bring resources and expertise into this project.

In the following, we present shortly the feature scope of elsA. We then present the drivers that led us to start building SoNICS. We then outline a few key points about its innovative architecture and the scope of new features and new capabilities of SoNICS compared with elsA.

FEATURE SCOPE OF THE CURRENT INDUSTRIAL ELSA SOLVER

The current software elsA has been in use to simulate internal and external flows for applications ranging from turbomachinery [START_REF] Plot | The High Level of Maturity of the elsA CFD Software for Turbomachinery Aerodynamics Applications[END_REF][15] [START_REF] Barrier | Linear Stability Analysis in Rotating Frames for Fan Blade Transition Prediction[END_REF] to aircraft [17][18], helicopters [START_REF] Lienard | RACER high-speed demonstrator: Rotor and rotorhead wake interactions with tail unit[END_REF][20] [START_REF] Boisard | Aerodynamic Analysis of Rotor/Propeller Wakes Interactions on High Speed Compound Helicopter[END_REF], missile and rockets for several industrial partners of the ASD domain. elsA offers flexibility in terms of geometry and meshes: structured, unstructured, hybrid [START_REF] Castillon | A hybrid structured/unstructured grid strategy for the CFD modeling of technological effects on complex turbomachinery applications[END_REF], chimera/overset [START_REF] Castillon | A Technological Effect Modeling on Complex Turbomachinery Applications With an Overset Grid Numerical Method[END_REF], body force [START_REF] Awes | Unsteady Body Force Methodology for Fan Operability Assessment under Clean and Distorted Inflow Conditions[END_REF] or Immersed Boundaries Method (IBM), see fig. 1.

Hybrid/unstructured Chimera Body force elsA contains efficient and robust finite volume methods with high order schemes like O(3) k-exact reconstruction [START_REF] Pont | Multiple-correction hybrid k-exact schemes for high-order compressible RANS-LES simulations on fully unstructured grids[END_REF], O(6) Lele [START_REF] Bras | Development of Compressible Large-Eddy Simulations Combining High-Order Schemes and Wall Modeling[END_REF] or O(5) MUSCL. It offers specific boundary conditions targeted at turbomachinery or propeller applications like mixing plane, chorochronic/phase-lagged [START_REF] Castillon | Multiple-Frequency Phase-Lagged Unsteady Simulations of Experimental Axial Compressor[END_REF], time spectral methods and reduced blade count approach, see fig. elsA contains turbulence models which capture complex phenomena with various levels of precision and cost, ranging from RANS to LES modeling [26][27]. It also offers numerous laminar-turbulence transition models, see fig. 3. In addition, elsA offers advanced CFD-CSM coupling capabilities for a large spectrum of aero-elastic simulations [START_REF] Placzek | Assessment of single passage simulations for the aeroelastic stability of a multistage centrifugal compressor[END_REF] and it handles shape optimization using linearized and adjoint methods [29][30]. elsA also interoperates with external HPC solvers and external users' scripts using python and the CGNS standard. Integration in industrial platforms also relies heavily on python and CGNS.

BUILDING AN INNOVATIVE CFD SOFTWARE

Motivations for building new software

For the past years, innovation has been mostly incremental within elsA although there have been noticeable improvements in numerical and physical modeling. However, there has been no major disruption within elsA regarding HPC and software architecture in the past years.

In the meantime, the ONERA CFD teams have acquired expertise in numerous fields. ONERA has been able to test very innovative software engineering and architecture concepts, or existing one which had not yet been used in industrial CFD software. The ONERA teams have developed prototypes in order to grasp all levels of parallelism on modern hardware architecture (distributed and shared memory, cache tiling [START_REF] Song | New tiling techniques to improve cache temporal locality[END_REF], loop nest optimization, vectorization, GPU), linear algebra (Krylov based methods and their preconditioning, communication optimal linear algebra [START_REF] Demmel | Communication-optimal parallel and sequential QR and LU factorizations[END_REF]), algorithmic differentiation [START_REF] Hascoet | The Tapenade automatic differentiation tool: principles, model, and specification[END_REF][34] [START_REF] Bradbury | JAX: Composable transformations of Python+NumPy programs[END_REF], data flow graph programming [START_REF] Wolf | Computers as components: principles of embedded computing system design[END_REF][37], software product line [START_REF] Damiani | A formal model for multi software product lines[END_REF] and automatic code generation and optimization [START_REF] Chelini | Progressive raising in multi-level IR[END_REF][40].

It has been quite complex to test some of the innovations mentioned above in an industrial scale software like elsA. Nevertheless, ONERA has been able to develop from elsA an HPC prototype, leveraging cache tiling, loop nesting optimizations and vectorization and has succeeded at obtaining a x7 performance speedup (with exact same machine, and exact same calculation) compared with elsA on implicit calculations on unstructured grids. However, the elsA developers have quickly realized that it would be extremely long and costly to manually extend these optimizations to the entire functional scope of elsA which contains millions of lines of code, thousands of features and highly complex cross compatibilities between these features. Moreover, the developers involved both in the production of industrial scale production software elsA and in the development of these innovations could realize that elsA with its current architecture would be unable to sustain the changes coming with the future hardware architecture of the exascale and post-exascale era.

Building SoNICS architecture

ONERA has decided to build an entirely new CFD software in order to overcome the limitations of elsA's architecture. SoNICS new architecture relies on several innovations, as summarized on fig. 4, which relies heavily on automation:

• Graph generation: each specific calculation with SoNICS depends on each user needs in terms of physical and numerical models, parameters, boundary conditions and also on the target machine (CPU/GPU/Vector engines/TPU clusters). For this specific calculation, a workflow based on an operator graph is generated. The graph will not be the same if the user wants to perform e.g. an Euler simulation on a structured grid on a CPU cluster, or if the user chooses to perform instead e.g. a k-w RANS calculation on an unstructured grid on an heterogeneous CPU/GPU cluster. Such a graph is created by an automatic graph generator. This generator automatically handles cross compatibilities between features and provides flexibility to integrate new numerical and physical models in the future. It eases coupling between solvers and provides the possibility to create external operators developed by users using high level python scripts. These external operators are automatically inserted into the operator graph.

• Code generation: a code generator is built to transform generic code developed in a high-level language into specific code with optimized loops and numerous HPC optimization for hitting all targets of every specific user context (mesh, modeling, target hardware). With this ingredient, we propagate HPC optimizations into the entire code and separate the HPC optimizations (into a so called HPC layer) from the code handling models, hiding the complications of such optimizations to developers building numerical or physical models.

• Automatic linearized and adjoint code generation: we make use of algorithmic differentiation throughout the entire SoNICS software to obtain the linearized code by differentiation of the direct graph and to get the adjoint code for optimization by differentiating the same graph in reverse order.

3.3

Feature scope of new SoNICS software

SoNICS targets an enlarged application spectrum compared with elsA, and aims at new capabilities which elsA could not offer because of its older architecture.

Enlarged application spectrum:

• SoNICS considers multiple fluid species, helping us simulate configurations with higher operating temperatures and pressure,

• SoNICS offers a shorter turnaround time thanks to dynamic mesh adaptation, initializing calculation on a coarse mesh and using metric-based adaptation to refine where necessary. The introduction of node centered finite volumes will enable adaptation with anisotropic mesh cells,

• SoNICS offers advanced Krylov based linear algebra, shape optimization, flow stability, flow control based on adjoint and embedded deep learning capabilities based on backpropagation (learning optimal parameters of numerical or physical models), all of these thanks to algorithmic differentiation. SoNICS natively handles automatic differentiation for linearized and adjoint code and for backpropagation, using forward and backward differentiation of the operator graph applied to the entire solver,

• SoNICS handles more easily cross-compatibilities between features thanks to its operator graph structure.

Shorter CPU time: with SoNICS, the generated source code is automatically optimized for each existing and future hardware. In addition to the design of an HPC layer leveraging all level of parallelism, SoNICS will strongly benefit from the impressive performance of current and future GPU hardware (e.g. NVIDIA's Hopper H100 and followers). This opens for example new possibilities for scale resolved LES and DES/ZDES simulations with SoNICS, and, together with optimization which we are applying to the entire distribution (solver plus external software tools), this will strongly reduce turnaround time at design offices for day-to-day RANS calculations in large design parameter spaces.

Precision and robustness: SoNICS's architecture eases integration of future numerical or physical models thanks to its operators graph structure and to the opportunity to fit external operators developed by users into the operator graph.

SoNICS also enables shorter development cycles because of the use of automation. The developer codes only once, and optimized source code will be generated for every type of hardware. Then algorithmic differentiation automatically delivers the linearized and adjoint code. This strongly reduces the number of lines of source code which need to be developed and maintained.

SoNICS will not only perform high fidelity calculations. It will also make use of low fidelity models like body force which we use for example for meridian initialization. It can also be coupled with low resolution models. Body force is one example of implementation of an external operator in python, which then fits in the operator graph and benefits subsequently from all properties (differentiation, HPC layer, cross compatibilities) of operators inside the operator graph.

BRIDGING THE GAP BETWEEN INDUSTRIAL NEEDS AND ONERA'S CFD OFFER

We have put recently a lot of efforts building a new, more agile organization in order to better satisfy Safran end users' needs. This has been possible because Safran brings sustained financial support into the organization but also directly participates in a lot of processes, in particular into industrialization.

This organization relies on three pillars presented in the following paragraphs.

Handling users' needs instead of delivering features

There exist high-level software platforms for the different industrial components (e.g. HP turbine, LP turbine, compressor, booster, nacelle, combustion chamber, …). Theses platforms have pre-tuned numerical and physical parameters and encapsulate the CFD solver and dedicated pre and post-processing tools tailored to each design need, see fig. 5. To feed the software development plan, instead of periodically centralizing a list of features to be developed in a given period of time, we have put in place an agile organization which aims at satisfying users' need. In fact, users are e.g. design engineers who are experts at their craft but less skilful in numerical analysis. In order to develop the right features in the right software layers, we organize teams around specific users' needs. Product owners collect users' needs, gather a team of developers for each software layer and select testers from the applied engineering/design offices teams. End users express their needs in the language of their craft. They ask for example "how can I perform an unsteady calculation for this multi-stage turbine with cooling holes". The product owner will decompose the need into features, and with the team of software developers and applied team testers, will decompose the resulting features into tasks to be performed in the various software layers, see fig. 6.

Having the end user involved at all steps

For each need, the end user is involved into all steps of the process. A use case, consisting of a test case and testing environment, is defined at the beginning of the process. Several meetings are organized periodically for each need which gather the end user, product owner, developers of the various software layer and testers.

Testing in the industrial end users' environment

ONERA does not ship only a solver and a few tools designed by its teams, but instead must deliver to the end users an entire distribution with solver and software tools designed by teams at different organizations (ONERA, Safran Aircraft Engines, Safran Helicopter Engines, Safran Aero Boosters, Safran Tech, other industries, academics and various sub-contractors) working at remote locations. One challenge is to ship a distribution of tools at each different layer (as shown above on fig. 5) which are coherent and operate together with perfect compatibility. The version of each tool is tagged at each release of the distribution (e.g. Distribution_v2.0 = elsA_v5.1.1 -Tool1_v1.0 -Tool2_v2.1 -Tool3_v0.5) and a set of tests must be passed in order to guaranty that each feature functions correctly.

To do so, we define use cases for each need which are industrial-scale test cases plus a well defined user environment, and simplify the test case in order to perform the development of each feature in the solver and software tools in ONERA's environment. These tests are stored in various databases both at ONERA and Safran. ONERA stores the simplified development test cases in a development non-regression basis. Each developer runs the entire non-regression basis at each source code commit. At each release of the distribution, each organization runs (usually performed by each tool integrator) its own non-regression basis and validation basis. These tests are performed on different hardware/operating systems.

The use of all tools in every layer of a common pivot python/CGNS standard eases a lot the development of interfaces between tools and the exchange of use cases between organizations. Sharing the same single standard also eases support for handling bugs and anomalies in the distribution.

In order to avoid issues when shipping distributions to Safran end users, ONERA adds to its validation basis larger size industrial test cases from Safran bases which are representative of industrial needs (e.g. multi-stage axial turbines, centrifugal compressor, axial compressors, …). However, Safran cannot share with ONERA most of its configurations because of confidentiality restrictions.

While ONERA ships the distribution of its tools to Safran three times a year, the planning of these releases is not synchronized with the acceptance tests of features and needs (illustrated above on fig. 7). In order to ship and test the distribution more often, we are planning to move the organization in the near future to Continuous Integration and Continuous Deployment (CI/CD). Opening access to common source repositories across organizations and providing access to shared computing resources for running acceptance cases are necessary steps which we are tackling right now.

The three principles above have guided the development of dashboards used by teams at Safran and ONERA to plan and follow the development of elsA and SoNICS so that these solvers satisfy users' needs.

RECENT APPLICATIONS WITH ELSA AND SONICS MATCHING USERS NEEDS

In the following paragraphs, we showcase a few recent new features developed in elsA thanks to the new agile organization. We then show the first milestones achieved with the newly born SoNICS.

A few recent features of elsA originating from Safran end users' needs

Because our partner Safran starts ambitious programs for new aircraft motor configurations, motor designers need many additional features in elsA in order to obtain better simulations representativity (geometry and mesh strategy, physical models like turbulence, transition, real gases, boundary conditions dedicated to turbomachinery, …), precision and CPU performance. We exhibit below a few recent features among numerous other features developed during the two past years.

Non-reflecting boundary conditions

Future motor configurations will be more compact with faster rotation speed. In order to simulate such configurations, it is difficult to avoid putting boundary conditions nearby blades, and artificial reflection of waves (shocks, rarefaction waves …) occur on these boundary conditions. Following and extending the theoretical work of Giles [START_REF] Giles | Nonreflecting boundary conditions for Euler equation calculations[END_REF], we have developed non-reflecting boundary conditions which strongly dampen spurious reflections on inlet, outlet and mixing plane boundary conditions, as shown on fig. This feature is currently validated both at ONERA and Safran using increasingly complex test cases (axial and centrifugal compressors) representative of design offices' needs.

Stage reduction boundary condition on unstructured meshes

Adapting the stage reduction boundary condition to unstructured meshes involves complex computational geometry algorithms, which have to be performed in distributed parallel environment. One needs to perform optimally mesh intersections on both sides of the join and distributed HPC communication strategies must also be optimized. This type of boundary conditions enables unsteady turbomachinery calculations. A large number of validation cases has been successfully passed before shipping this feature to the design offices who will shortly use it in production. Figure 9 illustrates validation on a Vega 3 test case.

CPU acceleration of implicit calculations with elsA

We illustrate a development where elsA has benefited from innovations of an early prototype of the SoNICS software. In this prototype, we have implemented an algorithm for ordering the unknowns along wavefront hyperplanes inside the mesh. These ordering and additional vectorization optimizations led to a x7 performance speedup of implicit LU-SGS calculations on unstructured meshes with the SoNICS prototype compared with elsA. We have transferred some of the ingredients of this algorithm (not all of them because of elsA design limitations) into the elsA solver. Calculations by Safran of a Tatef compressor (fig. 10) exhibit a x2.6 performance speedup while ONERA ZDES calculations of a shear layer exhibit a x3.1 performance speedup.

First milestones achieved with SoNICS

We have finalized in the past months many of the key elements of the SoNICS architecture. A few milestones have been successfully achieved. We have run test cases which concern HPC and numerical schemes:

• An operator graph and numerical operators have been generated both for CPU and GPU calculations.

• Two versions of the numerical schemes, one cell-centered, the other node-centered, have been set up. This second version allows for more robust calculations on anisotropic unstructured meshes which we target in our dynamic mesh adaptation strategy.

Figure 11 shows an Euler calculation with a Roe scheme on an M6 wing at low Mach number using a cellcentered discretization.

In addition, we have measured notable acceleration of the HPC performances of SoNICS compared with elsA:

• On CPU processors: we have measured a x14 acceleration between SoNICS and a 2020 version of elsA (on same hardware, same configuration, same calculation) thanks to the optimizations contained in the HPC layer.

• On GPUs: a x100 acceleration compared to the same calculation on CPUs. This is the expected result which corresponds grossly to the flops difference between the two hardware architectures. STO-MP-AVT-366

CONCLUSIONS

We have shown in this paper the drivers which have led ONERA to switch to the development of an innovative CFD software named SoNICS, while continuing development of the current elsA solver for the medium term needs of our partner Safran. While we progressively scale down elsA, SoNICS will evolve to an industrialscale software with a feature scope larger than elsA. SoNICS will also bring breakthrough capabilities which we believe will provide a competitive edge to Safran for its future motor programs. The robust organization around software development and industrialization will help aligning the software distribution to users' needs.

While we continue on our roadmap and increase momentum by switching more skilled developers and application engineers to SoNICS, progressively achieving major development milestones, we are confident that the differentiating capabilities brought by SoNICS will help Safran and other ONERA partners achieve successes in future industrial programs.

Figure 1 :

 1 Figure 1: various meshing/geometry strategies within elsA.

2 .Figure 2 :

 22 Figure 2: multi-chorochronic calculation with elsA.

Figure 3 :

 3 Figure 3: left: laminar/turbulent transition on nacelle using parabola model, right: comparison between SAS DRSM and ZDES model.

Figure 4 :

 4 Figure 4: generating specific code from generic code.

Figure 5 :

 5 Figure 5: software layers with various types of tools. Industrial end-users have access to the top-level tools with pre-tuned parameters and models. Developing a coherent distribution of tools involves coordination between end-users, developers and testers at the various levels.Coordination is handled by product owners.

Figure 6 :

 6 Figure 6: users' needs are split into features and tasks. Tasks are performed in various software layers.

  turbine with cooling holes Acceptance tests are performed on the end user machine within the software chain of the design office dedicated to the aircraft component. Planning is performed by splitting the need into features, themselves being split into tasks, see fig 7. Tasks are assigned to individuals or group of individuals. Because individual developers and testers work in various software layers and at different organizations, coordination and synchronization by the product owner is paramount.

Figure 7 :

 7 Figure 7: planning for team handling a specific need decomposed into features and tasks in the various software layers.

Figure 8 :

 8 Figure 8: left: older version of elsA exhibiting reflection of shocks on outlet boundary condition located near blade row, right: April 2022 version of elsA with dampened reflection on the outlet condition

Figure 9 :

 9 Figure 9: stage reduction join validated on Vega 3 case.

Figure 10 :

 10 Figure 10: left : waverfront ordering of the unknowns, left : Tatef compressor calculation showing that the current version of elsA is 2.6 times faster than older (before 2020) versions f l A hi fi i

Figure 11 :

 11 Figure 11: Euler calculation of M6 wing using SoNICS

a New Generation CFD Software for Satisfying Industrial Users Needs 10 -10 STO-MP-AVT-366

  [START_REF] Content | Application of the γ-Reθt laminar-turbulent transition model in Navier-Stokes computations[END_REF]. With these new boundary conditions, we obtain a strong dampening of the shock reflection together with almost negligible discrepancies on quantities of interest compared with reference calculations on a very long domain: flow rate < 0.05%, compression rate < 0.09% and efficiency < 0.47%.

SoNICS:

-2 STO-MP-AVT-366

-4 STO-MP-AVT-366

-6 STO-MP-AVT-366

STO-MP-AVT-366

-8 STO-MP-AVT-366

-14 STO-MP-AVT-366

-16 STO-MP-AVT-366

Acknowledgments: we thank the participants of the SONICE project and Safran and ONERA management for valuable input and discussions. Especially, we thank all those involved in the development of elsA and SoNICS, who worked on producing the results and many of the illustrations presented in this paper. We finally thank DGAC for its support within the framework of France Relance.