Julien Vialle

Arnaud Badel

Pascal Tixador

2 T insulated REBCO magnet used as 1 MJ SMES: protection strategies and preliminary assembly tests

Keywords: SMES, Superconducting magnet, Insulated REBCO magnet, High Magnetic field, Protection system

In the framework of the BOSSE project, a 12 T insulated REBCO solenoid is being manufactured. This magnet will be used as a SMES with the objective to reach 1 MJ and a specific energy of 20 kJ/kg for the winding. To manufacture this solenoid and reach such performances, 21 insulated REBCO Double Pancakes (DPs) will be assembled and cooled in a liquid helium bath at 4.2 K. To validate the solenoid's electromagnetic design, a prototype DP tested in self-field up to its critical current (972 A), thanks to a sensitive protection system, was presented in a previous work. Here we present the background magnetic field test of this prototype DP which allowed to validate the coil's mechanical design. In order to validate the electromagnetic performances of the 21 DPs, each of them is tested in liquid helium up to its rated current. The performances of the 17 DPs already tested and validated will also be presented. Finally, we will present the tests and results of 2 preliminary assemblies of 3 and 5 DPs. These assemblies were tested with two different protection philosophies. The results of these tests show that an individual protection of each DP is to be preferred in order to obtain a better sensitivity on the protection signal.

I. INTRODUCTION

12 T REBCO solenoid with a 188.5 mm wide inner bore tapes was designed in the framework of the BOSSE Project, using 12 mm-wide single insulated SuperOx tape. This 814 mm high coil can be used as a compact 1 MJ -2 MW inductive energy storage, with an 850 A rated current. A first objective is to surpass the record for mass energy density, held by a NbTi coil (13 kJ/kg) [START_REF] Yamamoto | A thin superconducting solenoid magnet for particle astrophysics[END_REF] by reaching 14 kJ/kg (mass of conductor only). The final goal is to, push the operating current further to reach 20 kJ/kg in the winding, which requires 890 A. Tab. I summarizes the different important quantities according to the objective in specific energy density to be reached.

This solenoid is composed of 21 Double Pancakes (DPs) with a soldered inner joint. Its geometry is optimized at the coil ends where the spacing between the pancakes is progressively increased in order to reduce the transverse field distribution and thus improve the critical current. The important rated values of the SMES for the most ambitious operation point, 20 kJ/kg, are listed in Tab. II. More specifications on the solenoid design and DP manufacture can be found in [START_REF] Vialle | Preliminary Tests of Pancakes From a 12 T REBCO Insulated Solenoid Magnet[END_REF].

Insulated REBCO windings are known to be difficult to protect in case of local thermal runaway. In order to discharge the magnet before destruction, it is necessary to detect at an early stage the beginning of the dissipation in the winding [START_REF] Badel | Detection and Protection Against Quench/Local Thermal Runaway for a 30 T Cryogen-Free Magnet[END_REF][START_REF] Ciceron | Test in Strong Background Field of a Modular Element of a REBCO 1 MJ High Energy Density SMES[END_REF]. This remains a delicate operation because the inductive voltage when ramping up the magnet is often several orders of magnitude higher than the voltage related to the onset of a transition. The electromagnetic noise itself is often higher in amplitude than the signal of interest. A detection concept using pickup coils to compensate the inductive component of the signal and some of the noise during current ramps, as well as numerical low pass filtering was introduced [START_REF] Ciceron | Test in Strong Background Field of a Modular Element of a REBCO 1 MJ High Energy Density SMES[END_REF]. The sensitivity of the quench detection system is in the 100 µV range, which made it possible to test a first full-scale prototype of DP at 4.2 K up to its limiting critical current (972 A with 0.18 A/s) without any degradation [START_REF] Vialle | Preliminary Tests of Pancakes From a 12 T REBCO Insulated Solenoid Magnet[END_REF]. The test of this DP prototype under background field is presented in part II, validating the mechanical design. This project started in 2015 and all conductor lengths were procured at that time. Considering the uncertainty on tape performances from this period in low temperature/high field and the low critical current margin of the design (14 % see [START_REF] Vialle | Preliminary Tests of Pancakes From a 12 T REBCO Insulated Solenoid Magnet[END_REF]), it was decided to test all DPs in liquid helium to the rated current before assembly. The results of the evaluation campaign for the 21 DPs are presented part III. Finally, part IV presents partial assembly tests. Our approach to operate safely the full-size magnet in spite of its high operating current density is to detect thermal runaway in each DPs independently or grouped by two or three, using the same detection concept used for the prototype pancake. Two partial assemblies of 3 and 5 DPs were tested in order to configure and validate the simultaneous protection of a growing number of DPs. During these tests, the transient voltage is carefully monitored to detect dissipation. The various phenomena contributing to this transient voltage have been studied with the help of detailed electromagnetic modelling results [START_REF] Vialle | Transient voltages and energy balance in REBCO insulated magnet: experimental and numerical studies[END_REF].

A

II. MECHANICAL DESIGN VALIDATION OF THE BOSSE MAGNET

According to the Virial theorem, the specific energy is proportional to the stress on the conductor [START_REF] Moon | The virial theorem and scaling laws for superconductingmagnet systems[END_REF]. In order to reach ultimately 20 kJ/kg the BOSSE solenoid must operate at high levels of mechanical stress and current density. Fig. 2 shows the stress distribution in the winding at the rated current. Fig. 2 shows the mechanical stress calculated with the JBR (independent turns) and Wilson formulas where the latter assumes "dependent" turns [START_REF] Arp | Stresses in superconducting solenoids[END_REF][START_REF] Wilson | Superconducting magnets[END_REF]. For these calculations, the rated values of the SMES were used. That is to say a current density of 525 A/mm², a magnetic field of 11.55 T at the inner turn and a magnetic field of -0.77 T at the outer turn. The mechanical stress calculated with Wilson is smaller at the inner turn and larger at the outer turn than with the JBR formula. This can be explained by the fact that the inner turns are subject to much larger Laplace forces and are therefore supported on the outer turns. In a way we can say that the mechanical stress is homogenized over the width of the winding in this magnetic field configuration. The radial stress is therefore negative and indicates a winding in compression.

From the conductor's point of view, it is mainly made of Hastelloy C 276, a very rigid material that has excellent mechanical properties with a high Young's modulus (YM) of 220 GPa. However, the copper stabilizer, which has a lower YM (124 GPa), weakens the mechanical properties of the conductor. This leads the conductor to an average YM of 172 GPa.

Calculations of the winding mechanical stress as well as the conductor YM allow us to estimate the conductor elongation. For the SMES to achieve the first target of 1 MJ, the conductor elongation on the central windings of the solenoid will reach 0.345%, and 0.38% for 1.1 MJ, corresponding to a specific energy density of 20 kJ/kg. This is very close to the maximum mechanical stresses allowed for the conductor since irreversible degradation may occur above 0.4 % in conductors from this period (2015) [START_REF] Barth | Electro-mechanical properties of REBCO coated conductors from various industrial manufacturers at 77 K, self-field and 4.2 K, 19 T[END_REF].

For this reason, before launching the production of the 21 DPs, it was necessary to validate the design and manufacturing process by testing electro-mechanical performances. Therefore, after the test in self-field operating current limit (mentioned in introduction), the prototype DP already was also tested in background magnetic field under larger mechanical stress. As the pancake is large (240 mm outer diameter), the 376 mm diameter resistive magnet (12 MW/30 kA/10 T) from the LNCMI in Grenoble was used. Additional pick-up coils were added to this test (ramp rate: 2 then 1 A/s Fig. 3) in order to compensate the important electromagnetic noise (hundreds of mV) generated by the resistive background magnet and to keep a high sensitivity on the protection signal [START_REF] Ciceron | Test in Strong Background Field of a Modular Element of a REBCO 1 MJ High Energy Density SMES[END_REF]. A photo of the prototype DP with its pick-up coils and the axisymmetric magnetic field map showing the positioning of each coil is shown in Fig. 1. The prototype DP reached a current of 623 A under 6 T, for a total magnetic field at the center of the DP of 7.6 T, and 9.7 T on the conductor (Fig. 3). We stopped the experiment at this value, as the stress on the outer diameter could become higher than in the final assembly, and the outer contacts are not designed for it [START_REF] Ciceron | Test in Strong Background Field of a Modular Element of a REBCO 1 MJ High Energy Density SMES[END_REF]. The stress could be concentrated on the outer contact ends (Fig. 4). Still, under 623 A the stress reached 350 MPa on the inner turns of the DP, corresponding to an elongation of 0.3% of the conductor (calculated data). This represent already 83% of the maximum hoop stress at the rated current for the SMES solenoid. The main cause of concern for the series-produced DPs is the inner joint. The inner resistance of the prototype was very low (< 20 nΩ), but the resistance values obtained in series-produced DPs was much higher. This is surprising as a soldered-joints validation campaign was organized with the industrial partner prior to the winding phase, to define the procedure and tooling to be used. During this campaign, joint resistance of 4 nΩ were routinely achieved, in good agreement with the prototype and literature [START_REF] Fleiter | In-Field Electrical Resistance at 4.2 K of REBCO Splices[END_REF]. It is suspected that a tension spike when starting the winding could be the origin, or an excessive time spent soldering the joint due to the more awkward position of the operator when soldering on an actual pancake. The average resistance of the 17 validated DPs reaches 367 nΩ, spreading from a minimum of 19 nΩ to a maximum of 700 nΩ. Several of them had higher values during their first tests and were rewound for this reason, and 4 are still waiting to be rewound.

IV. PROTECTION OF DPS IN ASSEMBLIES

In principle, two voltage taps across the superconducting winding and a single pickup coil for compensation are sufficient to ensure the protection of any DP assembly, just as it is done for a single DP. Nevertheless, the length of the winding between two voltage taps is important. Indeed, the amplitude of the residual voltage noise increases with the winding length. On the contrary, dissipative state in REBCO tapes will appear very locally, so that the amplitude of the dissipative voltage that needs to be detected does not depend on the coil length.

In order to increase sensitivity, it was decided for this project to equip each DP with voltage taps and an independent pick up coil. This makes it possible to protect the DPs either independently or by sub-groups. Individual protection of DPs also has the important advantage of accurately locating a faulty DP among the assembly in case of an early protection discharge event, making it possible to replace it or repair it.

For the final assembly the protection system presented in [START_REF] Vialle | Preliminary Tests of Pancakes From a 12 T REBCO Insulated Solenoid Magnet[END_REF] will therefore be multiplied by up to 21 (Fig. 5). The compensated voltage signals of each DP or sub groups of DPs are analyzed simultaneously in real time in order to detect the occurrence of a dissipative zone. The BOSSE solenoid will be discharged under 5 kV (± 2.5 kV) to get to its rated power; the acquisition electronics is disconnected when the discharge is triggered to protect it from overvoltage, making it possible to operate without high voltage differential units that are costly and generate noise [START_REF] Badel | Detection and Protection Against Quench/Local Thermal Runaway for a 30 T Cryogen-Free Magnet[END_REF].

The signals from the DPs and subassemblies tested in selffield, which are presented in this paper, were monitored using an oscilloscope controlling the protection switch through a 0 -5 V TTL output. The compensated voltage signals were manually adjusted using a voltage divider bridge. This method Stress concentration is easily manageable for testing a single DP or "small" assemblies. For future, larger assemblies, the signals will have to be treated automatically, taking into account the transient components of the signals as well as the voltage peaks due to winding movements. Modeling efforts have been made in this direction [START_REF] Vialle | Transient voltages and energy balance in REBCO insulated magnet: experimental and numerical studies[END_REF][START_REF] Rozier Blandine | Calculation of the local current density in hightemperature superconducting insulated rare earth-barium-copper oxide coils using a volume integral formulation and its contribution to coil protection[END_REF].

V. TEST OF A 3 DPS ASSEMBLY

A. Assembly characteristics and test result

The main step left before the final magnet is to succeed in protecting several DPs in an assembly. Two partial assemblies were tested. The first one used 3 DPs. This assembly was made using DPs from the BOSSE solenoid heads, so they are separated by a larger gap than the central DPs. Moreover, these coil head DPs were wound with the tape lengths having the best average Ic at 77 K (>500 A) to keep sufficient current margins in transverse field. The 2 DPs at the ends of this assembly have a spacing of 30 mm between pancakes and the DP at the center a spacing of 15 mm. Between these three DPs a 15 mm spacing is chosen to homogenize the field distribution as well as possible. As in the final assembly, each DP is connected to the adjacent one by their outer contacts, which are pressed together using three M4 screws. Indium is used to reduce the contact resistance. An axisymmetric magnetic field map of the assembly of these 3 DPs and their respective pick-up coils is shown in Fig. 6 as well as a picture of a DP pancake coil where the inner and outer contacts of the DP can be distinguished. Fig. 6. The magnetic field map of the 3 PDs assembly and its pick-up coils as well as the picture of a DP pancake winding and its outer crescent contact.

The pickup coils are inserted in the center of the solenoid. This saves space and minimizes the size of the cryostat, at no cost for SMES application, as the center bore does not have to be accessible for energy storage. They are however designed to be removable if the full magnet was to be used for other purpose in the future. They have 1200 turns of Φ 0.36 mm insulated copper wire with an inner radius of 56 mm, an outer radius of 63 mm and a height of 28 mm. The number of turns of the pickup coils is much larger than the number of turns of the DPs, in order to enable full inductive voltage cancellation even though the coupling is limited due to the large thickness of the DPs mandrels.

For this first assembly, we evaluate the possibility to protect the three DPs together. Their three pick-up coils were connected in series in order to perform a global compensation of the assembly and analyze the sensitivity of the protection signal. In order to minimize the risks of damage in such preliminary work, it was decided to limit the operating current to 745 A, corresponding to the first project target of 14 kJ/kg. We exceeded the first target of 745 A and reached 763 A (see summary Table IV) which corresponds to a stored energy of 52.3 kJ and a specific energy density of 6.35 kJ/kg. A first current ramp at 2 A/s is performed up to 600 A where a short plateau of about 100 s is maintained (to evaluate steady state dissipation), then the current ramp is increased with a ramp of 1 A/s up to 763 A where another plateau of about 160 s was achieved before the magnet discharge. During both current plateaus, the voltage returned close to zero, proving the safe operation of the assembly. Fig. 7 shows the individual voltages of the 3 DPs as well as current cycle of this test. The voltages of the DPs at the ends of the assembly are expected to be superimposed as they are identical and symmetrical. The slight difference observed in their voltage is due to a small variation in turn number: 360 for the upper one and 368 for the lower one. The voltage of the central DP is logically higher: its inductance is higher due to its smaller gap between pancakes (15 mm instead of 30 mm). It also has higher coupling with both others. Fig. 8 shows the compensated voltage signal of the assembly used to protect the winding (in blue). The triggering threshold outer contact of the protection is set at 1.5 mV. This triggering value is fixed based on simulations of the thermal runaway phenomenon developed in our group [START_REF] Badel | Detection and Protection Against Quench/Local Thermal Runaway for a 30 T Cryogen-Free Magnet[END_REF][START_REF] Badel | Modeling of 'quench' or the occurrence and propagation of dissipative zones in REBCO high temperature superconducting coils[END_REF]. Of course, the lower the transition detection threshold, the longer and more comfortable the time for the protection to be triggered. The voltage drift observed on the protection signal from the beginning of the ramp on the whole voltage profile is not due to a resistive component. It is due to the variation of the current density distribution in the width of the superconducting layer of the conductor as the current or the magnetic field varies, inducing a change of the superconducting coil inductance as well as AC losses. As such, it rapidly decays to zero during a current plateau (stabilization of the current density distribution). This phenomenon, first observed in [START_REF] Ciceron | Test in Strong Background Field of a Modular Element of a REBCO 1 MJ High Energy Density SMES[END_REF], has been modelled and studied in detail [START_REF] Vialle | Transient voltages and energy balance in REBCO insulated magnet: experimental and numerical studies[END_REF][START_REF] Rozier Blandine | Calculation of the local current density in hightemperature superconducting insulated rare earth-barium-copper oxide coils using a volume integral formulation and its contribution to coil protection[END_REF]. This drift forces us to perform a dynamic signal compensation. This explains the discontinuities in the protection signal curve Fig. 8. Each of these discontinuities correspond to a readjustment of the compensation in order not to exceed the triggering threshold due to this signal drift.

The remaining voltage of about 250 µV at the end of the plateau at 763 A is due to the inner and outer connection of the windings which are resistive. This gives us a global resistance of the assembly of 327 nΩ at this operating current. The outer pressed contacts connecting the DPs together are very good and paradoxically better than the inner soldered contacts of the DPs. The latter are still satisfying, with an estimated resistance of about 60 nΩ in average.

The main result of this assembly test is that it was possible to reach sufficient sensitivity (well below 200 µV) in the thermal runaway detection system when observing a group of three DPs. This way the number of signals to monitor in a larger assembly can be reduced. Of course, the drawback of this solution is that identifying the defective DP is impossible if the protection triggers during a current transient. Moreover, an excellent energy efficiency, 98.6 %, was obtained over the full chargedischarge cycle of this assembly (the discharge time constant was set to 6 s, slower than the 1 s target in the final configuration [START_REF] Vialle | Preliminary Tests of Pancakes From a 12 T REBCO Insulated Solenoid Magnet[END_REF]). The losses are AC losses and connection losses.

B. Inner contact damage and discussion

During a discharge of the 3 DPs assembly test, one of the inner contacts was degraded. Nevertheless, we reached several times this operating current of 763 A showing that the superconducting winding was not damaged. As mentioned in part III, the inner contacts were a source of concern from the beginning of the series DP manufacturing, with joint resistance measured well above what was expected from the prototype. Still, this degradation occurred at the discharge and therefore neither at the maximum applied current nor stress. The reason may come from the inner contact structure (Fig. 9). It was originally designed with large copper thickness to guarantee thermal stability, acting as a radiator in the helium bath. Eddy current heating during fast discharge was expected but the estimated temperature rise was acceptable [START_REF] Ciceron | High energy density Superconducting Magnetic Energy Storage with second generation high temperature superconductors[END_REF]. Superconducting strips were added on the copper to reduce the joint resistance, but this may have been counterproductive. The improvement in terms of joint resistance was negligible but during transients, the induced current circulating in the copper contacts pass through these strips. When the transport current is decreased, this induced current is added to the transport current, with two potentially detrimental effects.

First, the total current may exceed the critical current of the strips. Second, this high current in the strips may cause abnormally high mechanical stress and the delamination of the conductor. A simpler, thinner inner contact without REBCO tape strips will be considered in the future.

VI. TEST OF A 5 DPS ASSEMBLY

A. Assembly characteristics and test result

The main goal of the 5 DPs assembly was to test the concept of individual protection of each DP. The deterioration of the inner contact on one of the DPs was judged acceptable to perform tests with reduced operating current, and a good test case for the protection system. As for the previous assembly the DPs corresponding to the heads of the complete solenoid were used. Starting from the bottom, the stack included DPs with inner gaps of 30, 15, 10, 15, and 30 mm. The spacing between the DPs starting from the bottom was the following: 15, 11, 11 and 15 mm. The position of the DPs and their respective pick-up coils can be found on the magnetic field map Fig. 10, with a picture of the assembly.

Triggering Threshold

Baxis

REBCO tape strips

We did not set any current target for this experiment, but rather tested the coil up to the current for which the dissipation appeared unstable (inner contact heating). A current ramp of 8 A/s was performed up to 592 A where a plateau of about 300 s was achieved (Fig. 11). You will find below a summary table of this test and Fig. 10 In order not to degrade any other inner contact of DPs a very slow full discharge of 320 s was performed. The display of the compensated voltage signals of each DP is shown in Fig. 12. We can see that the signals of 4 DPs return to 0 V once the plateau has started except for one corresponding to the degraded DP of the previous assembly. The voltage is quite high and approaches 4 mV. Yet it remains stable, which means that it is due to the inner contact resistance of the DP. If a thermal runaway was occurring it would have shown another dynamic. However, the dissipation of this contact was so high (2.37 W) that we considered it too dangerous to go higher in current.

The sensitivity of the individual DP compensation signals is very good: the noise envelops after compensation and filtering (numerical low pass with 15 Hz cut-off) is about 45 µV.

B. Discussion

Here, in this assembly of 5 DPs, an individual protection of each DP has been successfully realized. We were able to perform several tests up to a current density of 365 A/mm² and this despite a degraded internal contact, which proves the robustness of the winding. The major benefit of this solution (in addition to the ability to clearly locate a defective DP) is the sensitivity level achieved below 50 µV. This advantage is all the more important as we are currently obliged to perform slow discharges. It is therefore of the utmost importance to have the finest possible transition detection sensitivity. This is why this protection solution will be adopted for the final solenoid protection.

VII. CONCLUSION

In this work, the electromechanical design of the BOSSE SMES was validated by testing a full scale Double Pancake under background field. All series-produced DPs were tested at 4.2 K up to their rated current for quality control. This test campaign study highlighted the difficulty to obtain low resistance consistently in an industrial setting with the chosen inner contact design. Further work will be conducted on that point. An important step was the development and validation of a detection/protection system to protect an assembly of several DPs. It highlighted the benefit of monitoring DPs in small groups or individually using multiple pick-up coils to obtain high detection sensitivity. The high sensitivity we obtained practically during the 5 DPs assembly testin the 50 uV range-makes us confident in the possibility to test the full assembly of the BOSSE solenoid reliably to its operation limits.

Fig. 2 .

 2 Fig. 2. Mechanical stress distribution in the winding of a DP at the center of the SMES for a stored energy of 1MJ. Hoop stress and radial stress are calculated with the Wilson formula.

Fig. 3 .Fig. 1 .

 31 Fig. 3. The current cycle and the compensated voltage signal during the test under magnetic field background of 6 T.

Fig. 4 .

 4 Fig. 4. Outer crescent contact structure

Fig. 5 .

 5 Fig. 5. Electrical circuit diagram of the BOSSE solenoid protection system.

Fig. 7 .

 7 Fig. 7. The current cycle and the individual voltage across each DP.

Fig. 8 .

 8 Fig. 8. The current cycle and the compensated voltage signal of the 3 DPs assembly

Fig. 9 .

 9 Fig. 9. Inner contact structure

5 Fig. 10 .

 510 Fig. 10. The axisymmetric magnetic field map of the 5 DPs and their pickup coils as well as a photo of the assembly.

Fig. 11 .

 11 Fig. 11. The current cycle and the voltage measured across the 5DPs assembly.

Fig. 12 .

 12 Fig. 12. The current cycle and the compensated voltage signals of the 5 DPs.

TABLE I :

 I OBJECTIVES TO ACHIEVE AT 4.2 K

	Energy	0.775 MJ	1 MJ	1.105 MJ
		14 kJ/kg	18.1 KJ/KG	20 KJ/KG
	I (A)	745	850	890
	J (MA/m²) B z max (T) B r max (T) Max. hoop stress (MPa)	460 10.14 3.7 306	523 11.53 4.2 396	550 12.12 4.4 437

TABLE II

 II

	: RATED SMES OPERATION VALUES FOR
	20 KJ/KG AT 4.2 K
	Energy	1100 kJ
	Height	814 mm
	Inner diameter	188.5 mm
	Outer diameter Bz max Br max	240 mm 12 T 4.4 T
	I operating	890 A
	Je	550 MA/m²
	Max. hoop stress (full cond.)	421 MPa
	Specific energy	20 kJ/kg
	Pancake number	42

TABLE III VALUES

 III REACHED DURING THE VALIDATION TESTS OF THE CENTRAL DPS (4 MM GAP)

	Average number of turns per pancake	180
	L (mH)	42.64
	I (A) Bz max (T) Br max (T)	863 5.28 4.17
	Max. hoop stress (MPa)	263


X Montiel

D Gusakova

M Daumens

A Buzdin

Proximity effect in atomic-scaled hybrid superconductor/ferromagnet structures: crucial role of electron spectra

Keywords: 74.45.+c -Proximity effects. Andreev effect, SN and SNS junctions 74.78.Fk -Multilayers, superlattices, heterostructures

We study the influence of the configuration of the majority and minority spin subbands of electron spectra on the properties of atomic-scaled superconductor-ferromagnet S-F-S and F-S-F hybrid structures. At low temperatures, the S/F/S junction is either a 0 or junction depending on the energy shift between S and F materials and the anisotropy of the Fermi surfaces. We found that the spin switch effect in F/S/F system can be reversed if the minority spin electron spectra in F metal is of the hole-like type.

1. Introduction-Phenomena arising in the superconductor-ferromagnet (S/F) heterostructures attract a growing interest due to their potential application in spintronic and quantum computation devices [1][2][3]. For example, the so-called spin-switch effect (SSE) (also called spin-valve effect) occurs when the critical temperature of the superconductor in the F/S/F structure depends on the mutual orientation of the ferromagnetic layers magnetizations. Moreover the S/F/S π-junctions with the changing superconducting order parameter phase shift are good candidates for the quantum computer elements, q-bit [4].

Critical current oscillations as a function of the exchange field and thickness of the ferromagnetic spacer in a S/F/S Josephson junctions have been predicted [5] and associated with the 0-π phase transition. Experimentally the evidence of such 0-π phase transition has been obtained first for the nanoscaled Josephson junctions S/F/S with a weak ferromagnetic interlayers in Cu x Ni 1-x [6] and PdNi alloys [7]. The 0-π phase transition was also predicted in layered compounds within a simple model of alternating F and S atomic layers [8][9]. This prediction is relevant for ruthenocuprates RuSr 2 GdCu 2 O 8 which are natural S/F layered compounds [10].

The F/S/F trilayer also exhibits interesting spin dependent phenomena. It was predicted [11][12][13] and observed [14][15][16][17][18][19] that the critical temperature is higher for antiparallel (AP) magnetization configuration than for the parallel (P) one. Besides this so called normal SSE, several recent experiments report on an inverse SSE whereby the parallel configuration is more favorable for the superconductivity than the antiparallel one [20][21][22][23][24][25].

The situation with inverse SSE investigated in [20] is somewhere controversial because in the similar Py/Nb/Py system in [21] the normal SSE was observed. The possible explanation of the inverse SSE observed in [20] has been recentlty proposed in [26] and it is related with a p-1 stray magnetic field generated in the AP configuration. The similar arguments may explain inverse SSE observed in [22,23]. However this reasoning do not work for the inverse SSE revealed in [24] where the magnetizations were perpendicular to the layers. therefoe due to the demagnetization the influence of orbital effect was excluded.

In this letter, we investigate the 0π phase transition and SSE for more realistic band structure. In particular, the π-phase might be suppressed when the ferromagnet bands are sufficiently energy-shifted with respect to the superconductor ones. The anisotropy of the quasi particle spectra is also of primary importance for the 0π phase transition. This could explain why the π-phase was not observed in Ruthenocuprates [27]. Furthermore, if the majority and minority spin subbands have opposite electron/hole character, we predict the inverse SSE providing a possible explanation of results [24].

2. Model -We start with an exactly solvable model [8] of alternating superconducting and ferromagnetic atomic metallic layers. The electron's motion is described in the F layers by the spin-dependent energy spectrum ξ σ (k) and by spin independent energy spectrum ξ s (k) in the S layers. Three basic parameters characterize the system : t is the transfer energy between the F and S layers, λ is the Cooper pairing constant which is assumed to be nonzero in S layers only, and h is the constant exchange field in the F layers only. It is supposed that the coupling between the layers is realized via the transfer integral t, which is relatively small (t ≪ T c ), so the superconductivity can coexist with ferromagnetism in the adjacent layers. The Hamiltonian of the system can be written as H = H 0 + H BCS + H t with:

H 0 = n,σ,k ξ n σ (k) ψ + σ,n (k) ψ σ,n (k) , (1)
H BCS = n,k ∆ * n ψ + ↓,n (k) ψ + ↑,n (-k) + ∆ n ψ ↑,n (k) ψ ↓,n (-k) , (2)
H t = t n,σ,k ψ + σ,(n+1) (k) ψ σ,n (k) + ψ + σ,n (k) ψ σ,(n+1) (k) , (3)
where ψ + σ,n (k) is the creation operator of an electron with spin σ and momentum k in the nth layer. The BCS pairing in the S layers is treated in H BCS in a mean-field approximation [START_REF] Abrikosov | Methods of quantum fiel theory in statistical physics[END_REF]. The superconducting order parameter ∆ n is non zero only in the S layers. Note that the electrons spectra in (3) are calculated from the Fermi energy. As usual, we introduce the normal and anomalous Green functions [START_REF] Abrikosov | Methods of quantum fiel theory in statistical physics[END_REF] k) which satisfy the system of equations:

G n,m σ,σ ′ = -T τ ψ σ,n (k) ψ + σ ′ ,m (k) and p-2 F n,m σ,σ ′ = T τ ψ + σ,n (k) ψ + σ ′ ,m (-
(iω -ξ n σ ) G n,m σ,σ ′ -tG n-1,m σ,σ ′ -tG n+1,m σ,σ ′ + ∆ * n F n,m -σ,σ ′ = δ nm , iω + ξ n -σ F n,m -σ,σ ′ + t F n-1,m -σ,σ ′ + t F n+1,m σ,σ ′ + ∆ n G n,m σ,σ ′ = 0, (4)
where ω = (2l + 1) πT c are the fermion Matsubara's frequencies, ξ n σ = ξ n σ (k) and n and m the layers indices. The superconducting order parameter in the nth layer satisfies the standard self-consistency equation

∆ * n = |λ| T ω k F n,m -σ,σ ′ . (5)
The anomalous Green function F n,m -σ,σ ′ of the system (5) can be be expressed by means of normal Green functions

F n,m -σ,σ ′ = p ∆ * n-l G ω;n-p,m -σ,σ ′ G -ω;p,m 0,σ,-σ ′ ,
where G 0 is the Green function in the absence of superconductivity pairing, i.e. for ∆ n = 0.

3. 0π phase transition in S/F/S system -In this section we examine the phase difference between the order parameters in the adjacent superconducting layers of the S/F/S system (see figure 1 (a)). Here n = 0, 2 and ξ s,0 = ξ s,2 = ξ s for the S layers, and n = 1 with ξ σ,1 = ξ σ for the F layer. Due to the symmetry reason the order parameters of both S layers may differ from each other in the phase prefactor e iϕ and |∆ 0 | = |∆ 2 | (in fact only the situation with ∆ 0 = ±∆ 2 is possible).

Solving the system of equation ( 4) in the case of the trilayer , we find the anomalous Green functions of the S layers. Using them, the self-consistency equation ( 5) can be written as :

∆ * 0 = -λT ω,k G ω;0,0 ↑↑ ∆ * 0 G -ω;0,0 0↓↓ + G ω;2,0 ↑↑ ∆ * 2 G -ω;2,0 0↓↓ . (6)
where

G 0,0 0↑↑ = (a * ↑ -t 2 ) (iω-ξ s )[(a * ↑ -2t 2 )] , G 2,0 0,↑↑ = t 2 (iω-ξ s )[(a * ↑ -2t 2 )]
, with a σ = (iω + ξ σ ) (iω + ξ s ). The transition between 0 and π states implies a change of the relative sign between ∆ 0 and ∆ 2 . Near T c , taking into account that in 0-phase, ∆ 0 = +∆ 2 and in the π-phase, ∆ 0 = -∆ 2 , the combination of the self consistency equations (6) written in 0-phase and π-phase provides one with the following expression ln

T 0 c T π c = 2T c0 ω,k G ω;2,0 0↑↑ G -ω;2,0 0↓↓ , (7)
where T 0 c is the superconducting critical temperature when the system is in the 0-phase, T π c is the superconducting critical temperature when the system is in the π-phase and T c0 is the bare mean-field critical temperature of the single superconductivity layer. The logarithm can be simplified as ln

T 0 c /T π c ≃ T 0 c -T π c /T 0 c = ∆T /T 0 c because the critical temperature variation is small T 0 c -T π c /T 0 c ≪ 1 .
The sign of the right hand side in (7) determines whether it is 0 or π-phase which is realized. Indeed, the transition occurs in the state with higher critical temperature. We are interested in the situation when the right hand side in (7) is negative, i.e. T π c > T 0 c and we deal with the π-phase difference of order parameters between the two superconducting layers.

p-3

First, we concentrate on the case with isotropic electron's dispersion when Fermi surfaces are circular. The calculation of the critical temperature, in the case of isotropic Fermi surface and symmetrical ferromagnet energy band splitting ξ ↑(↓) = ξ s ∓ h has been performed in [9]. Here we address the more general case with ξ ↑ = ξ s + E ↑ and ξ ↓ = ξ s + E ↓ where E ↑(↓) is the energy difference of the electron energy spectrum with spin up (down) compared to the superconductor energy spectrum ξ s . The calculation of the sum over momentum k in (7) may be transformed into an integration over ξ s , i.e. k -→ N (0) dξ s dθ, where N (0) is the density of state at the Fermi energy E F . Performing the integration over energy ξ s we confronted with several possible situations according to the values and signs of E ↑ and E ↓ that may be classified as follows : 1)

E ↑ ≫ T c0 , but E ↓ is of the order of T c0 , 2) E ↑ ≫ T c0 , E ↓ ≫ T c0 and 3) E ↑ ≫ T c0 , E ↓ ≪ T c0 .
In the case 1) when E ↑ ≫ T c0 and E ↓ is of the order of T c0 , the electrons up energy band is strongly shifted in comparison to the superconductor one and electrons down have energy band close to the superconductor one. The critical temperature difference in the limit t ≪ T c ,with ω > 0, becomes :

∆T T 0 c = 7πt 4 ζ(3) 8E ↑ E ↓ π 3 T 2 c0 - t 4 8E ↑ E ↓ 2 πTc0 Im Ψ ′ 1 2 -iE ↓ 4πTc0 . ( 8
)
where Ψ ′ is the first derivative digamma function. From ( 8) we deduce that ∆T = T 0 c -T π c > 0 and the S/F/S system is in the 0-phase if E ↑ and E ↓ have both the same sign. However, ∆T < 0 and the S/F/S system is in the π-phase if E ↑ and E ↓ have opposite sign e.g. E ↑ positive and E ↓ negative. Thus, the existence of the π-phase depends on the energy band shift between the ferromagnetic spectra and the superconductors spectra.

If we switch the role of E ↓ and E ↑ in ( 7) i.e. E ↓ ≫ T c0 and E ↑ is of the order of T c0 , the conclusions are similar. When E ↑ and E ↓ have both the same sign, the S/F/S system is in the 0-phase. When E ↑ and E ↓ have opposite sign, the S/F/S system is in the π-phase.

In the case 2), when E ↑ ≫ T c0 and E ↓ ≫ T c0 , the electrons up and down have both a strong energy band shift compared to the superconducting energy spectrum. The critical temperature difference becomes ∆T /T 0 c = 7πt

4 ζ (3) / 8E ↑ E ↓ π 3 T 2 c0
. The π-phase only appears if E ↑ and E ↓ have an opposite sign.

In the case 3) when E ↑ ≫ T c0 and E ↓ ≪ T c0 , the electrons with spin up have a strong energy band shift compared to the superconducting energy spectrum and the electrons with spin down a very small one. The critical temperature difference becomes ∆T /T 0 c = 31E ↓ t 4 ζ (5) / 64π 4 T 4 c0 E ↑ . Again the π-phase appears if E ↑ and E ↓ have an opposite sign. Here, we note that ∆T depends linearly on E ↓ . As a consequence, when E ↓ goes to zero (superconductor and electrons down ferromagnet energy band coincides) then T 0 c = T π c in our approximation and to determine the type of the ground state it is needed to calculate ∆T at higher order over t approximation.

As a conclusion, the π-phase does not exist if the energy band edges for both spin orientation of the ferromagnet are higher or lower than that of superconductor one. Consequently, the energy shift between the ferromagnet electron band and the superconductor electron band strongly influences the π-phase existence.

In real compounds, the results depends on the details of the exact form of the Fermi surfaces of superconductor and ferromagnet. Hence, the anisotropy may have an effect on the 0π phase transition.

To illustrate the anisotropic case, we choose elliptic Fermi surface of F layer. The majority and minority energy spectra ξ ↑ (θ) and ξ ↓ (θ) depend on the polar angle θ in the (k x , k y ) plane (see figure 2).

In the case of strong energy band splitting between ferromagnet and superconductor electron bands, the spectrum of majority (minority) ferromagnet electrons can be written as ξ ↑(↓) (θ) = ξ s +E ↑(↓) (θ) where ξ s is the superconducting spectrum considered as isotropic one. In spite of the dependence on θ, the situation depicted in the figures 2 a) and 2 b) p-4 are similar to the isotropic situation and we can estimate the sign of ∆T without doing the integration over θ. Thereby, the conclusions are the same as in the isotropic case.

From ( 8), we deduced that ∆T > 0 and the S/F/S trilayer is in the 0-phase if E ↑ (θ) and E ↓ (θ) have both the same sign for all values of θ e.g. E ↑ (θ) > 0 and E ↓ (θ) > 0 which corresponding sketch of Fermi surfaces is presented on the figure 2a). Thus, ∆T < 0 and the S/F/S trilayer is in the π-phase if E ↑ (θ) and E ↓ (θ) have opposite sign for all values of θ e.g. E ↑ (θ) > 0 and E ↓ (θ) < 0 which corresponding sketch of Fermi surfaces is presented on the figure 2b). Consequently, the anisotropy of electrons band do not influence the 0π transition phase when F layers and S layers Fermi surfaces do not intersect.

Let us consider now the case when the F layers and S layers Fermi surfaces intersect (see figure 2 c)). In this case, the intersections delimit the regions where ∆T is positive or negative. In order to estimate the resulting sign of ∆T without integrate over θ, we divide the Fermi surface in three regions. And then, according to the total surface of "positive" and "negative" region, one may conclude which state is realized.

In the regions (1) and (3), both E ↑ (θ) and E ↓ (θ) have the same sign (positive in (1) and negative in (3)). Consequently, ∆T is positive and the S/F/S junctions exhibit a 0-phase behavior in these regions. In the region (2), E ↑ (θ) and E ↓ (θ) have opposite sign (E ↑ (θ) > 0 and E ↓ (θ) < 0) so, ∆T is negative and the S/F/S junctions exhibit a π-phase behavior in these regions.

The sign of the contributions of the each region to ∆T is depicted in the figure 2 (d). If the overall size of the regions where π-phase is realized is greater than the overall size of the region where the 0-phase exists, we can conclude that the phase difference between the two superconducting layers is π. In the contrary case, the phase difference between the two superconducting layers is 0. These results implies that the anisotropy of the Fermi surfaces with intersection and energy band shift has an important effect on the 0-π phase transition.

We may generalize the S/F/S hybrid system by considering an arbitrary number of F layers between the two S layers. In [START_REF] Andersen | [END_REF], the authors found that the 0π phase transition of superconductor-antiferromagnet-superconductor (S-AF-S) junctions manifests a dependence on the number of magnetic atomic layer. The AF interlayer is composed by an even or odd number of monoatomic F layer where each F adjacent layer have an opposite magnetization. At low temperatures the junction S-AF-S is either a 0 or π junction depending on wether the AF interlayer consists of an even or odd number of atomic layers. In this paper, we show that in fact the 0π phase transition depends only on the number of ferromagnetic layers between the two S layers but not on the relative orientation of the magnetization between adjacent F layers. We generalize our model by adding N -2 F layers between the two S p-5 layers with n = 0 and n = N . For simplicity, we consider that the ferromagnet spectrum can be written as ξ ↑ = E ↑ and ξ ↓ = E ↓ where E ↑ and E ↓ have opposite sign and E ↑(↓) ≫ T c0 . Thus in the limit t ≪ E ↑(↓) we can express the normal Green function between the layer 0 and the N th layer like :

G ω 0↑↑ (0, N ) = t (iω-ξ s ) 2 N -1 n=1 t E ↑ n .
The critical temperature difference following the formular similar to (7) is :

∆T T 0 c = 7t 2 ζ(3) 8π 2 T 2 c0 N -1 n=1 t 2 E ↑ n E ↓ n . ( 9
)
As it follows from ( 9) the insertion of each supplementary F layers leads to additional phase shift of π. As a consequence, the resulting phase difference (0 or π) is determined by the number of ferromagnetic layers between the superconductors. If the number of ferromagnetic layers is even the total phase difference ϕ t induced by the ferromagnetic layers is ϕ t = 2nπ that correspond to a 0-phase shift. For the odd numbers of ferromagnetic layers the total phase difference is ϕ t = (2n + 1) π that means the π-phase shift. Hence, the magnetization of adjacent ferromagnetic spacers have no influence on the fact that the system is in 0 or -phase in contrast with the result in [START_REF] Andersen | [END_REF]. 4. Spin-switch effect in F/S/F system -In this section, we study the possible inversion of the SSE in a superconductor layer sandwiched by two ferromagnetic adjacent layers in the F/S/F structure (see figure 1b)). We label n = 0, 2 so ξ σ,0 and ξ σ,2 are the electron spectra for the F layers and n = 1 so ξ s,1 = ξ s for the S layer. The BCS Hamiltonian H BCS and the kinetic energy Hamiltonian H 0 are described by (3). In [24], the authors observed the inversion of the SSE (they observed that the critical temperature in the parallel state is higher than in the antiparallel one, T cP > T cAP ) in [Co/P t] /N b/ [Co/P t] structure and proposed an explanation for the inversion of the SSE. In the parallel (P) state the spin-polarized carriers can migrate from one F layer across the intermediate S layer to the other F layer whereas in the antiparallel (AP) state this effect is reduced and the spin polarized carriers is reflected back into the S layer. This means that the transparency of the S-F interfaces depends on the electron spin orientation. In our model, that implies that the transfer energy between two layers may depend on the electron spin orientation so in (3), instead of t one, we introduce two different energy transfer t ↑ and t ↓ with t ↑ = t ↓ . Near T c , in the isotropic case, the critical temperature equation is [START_REF] Abrikosov | Methods of quantum fiel theory in statistical physics[END_REF] :

ln TcP TcAP = 2T cP ω>0 Re F ↓↑,P -F ↓↑,AP dξ s , (10)
where T cP (AP ) is the superconductor critical temperature in the P (AP) case. Using the Green function formalism and the equations (4) modified with the new tunneling Hamiltonian H P (AP ) t

, we calculate the anomalous Green function in the S layer in the P state F + ↓↑,P (1, 1) and in the AP state F + ↓↑,AP (1,1). Expanding this function in power of t ↑ and t ↓ , and keeping the terms proportional to t 4 , we obtain the difference of the anomalous Green function F + ↓↑,P -F + ↓↑,AP :

F ↓↑,P -F ↓↑,AP = 1 ω 2 +ξ 2 2 (iω+ξ) 2 t 2 ↑ iω+ξ ↑ - t 2 ↓ iω+ξ ↓ 2 -1 (ω 2 +ξ 2 ) t 4 ↑ (ω 2 +ξ 2 ↑ ) + t 4 ↓ (ω 2 +ξ 2 ↓ ) + 2t 2 ↑ t 2 ↓ (ω 2 +ξ 2 )(iω-ξ ↑ )(iω+ξ ↓ ) .
We choose for simplicity, the energy spectra for the P case like ξ ↑ ≡ ξ ↑,0 = ξ ↑,2 and ξ ↓ ≡ ξ ↓,0 = ξ ↓,2 , and for the AP case, ξ ↑ ≡ ξ ↑,0 = ξ ↓,2 and ξ ↓ ≡ ξ ↓,0 = ξ ↑,2 . The sign p-6 of the right hand part in ( 10) determines whether we deal with normal or inverse SSE. In particular, if ln (T cP /T cAP ) is positive then T cP > T cAP and so the SSE is reversed.

To verify the hypothesis of [24] , we consider first the spectra ξ ↑ = ξ ↓ = ξ s , see figure 3 b) and we find :

∆T TcP = -31 64 ζ(5) π 4 T 4 cp t 2 ↓ -t 2 ↑ 2 (11)
where we have introduced for this section ∆T = T cP -T cAP and ζ (5) = 1.03. We see that the right hand side of equation ( 11) is always negative. It is the normal SSE. We choose another spectra ξ ↑ = ξ s and ξ ↓ = ∆E with ∆E ≫ T c , then the critical temperature difference becomes :

∆T TcP = -1 64 [31ζ(5)t 4 ↑ ∆E 2 +112t 4 ↓ ζ(3)π 2 T 2 cP ] ∆E 2 π 4 T 4 cp . (12)
The right hand side of ( 12) is always negative so T cP < T cAP . Consequently, the superconductor critical temperature in the P case is always lower than in the AP case, T cP < T cAP . It is the normal SSE. Thus, the hypothesis on the difference between t ↑ and t ↓ cannot explain the inversion of the SSE. Nevertheless, in strong ferromagnetic, the inversion of the electron spectrum provides another possibility to explain this unusual behavior.

Indeed, in the case of strong ferromagnet, the band splitting can be so important that one of the electron band (here the spin up band) becomes a hole-like (see figure 3 c)). Here, the exchange field triggers the hole-like energy band for spin up electrons but do not appear explicitly in the calculation. That's allows us to emphasize the important influence of inverse spectra on T c . In the very special case where the electron spectra can be written as ξ ↑ = -ξ s and ξ ↓ = ξ s , we find :

∆T TcP = 31 128 ζ(5) π 4 T 4 cp 3t 4 ↑ -2t 4 ↓ -4t 2 ↓ t 2 ↑ . (13)
If we put t ↑(↓) = t 0 ± ∆t then ∆T is positive (T cP > T cAP ) and the spin-switch effect is reversed if ∆t 0.1t 0 (see figure 4(a)). We notice that the spin-switch effect is reversed if t ↓ = 0 and normal if t ↑ = t ↓ . Hence, we verify that the condition t ↑ = t ↓ cannot explain alone the inversion of SSE and have to be associated with the inversion of majority electron energy band. The minority electron band can be shifted considerably compared to the superconductor spectrum (see figure 3 d)). In this case, the spectra are ξ ↑ = -ξ s and ξ ↓ = ∆E ≫ T c . This p-7 situation gives :

∆T TcP = 93 128 ζ(5) π 4 T 4 cp t 4 ↑ - 7ζ(3) 4(∆E) 2 π 2 T 2 cP t 4 ↓ ( 14
)
Therefore ∆T is positive (T cP > T cAP ) for high ∆E and any relation between t ↑ , t ↓ and the opposite situation realizes if ∆E is small (see figure 4(b)). In the situation of high ∆E, the inversion of SSE is related to the inversion of spectrum of majority electrons. We see that strong energy shifts ∆E between F and S layers favored the inversion of the SSE.

We perform numerical calculation of ∆T in the case of spectra like ξ ↑ = -α.ξ s + ∆E and ξ ↓ = ∆E ≫ T c where α is a positive parameter which represents the difference of effective mass between majority and minority electron and ∆E is energy shift between electrons up and the superconductor spectrum. For simplicity, we consider ∆E ≃ T c and t ↑ = t ↓ = t. For important or small difference of effective mass, there is no inversion of SSE (see figure 4(c)). The inverse SSE is favored when α ≈ 1. Hence, the inversion of the SSE is more efficient if the effective mass are close.

The dependence on ∆E (see figure 4(d)) shows that the inverse SSE disappears if ∆E 2T c . Thus, if the energy shift on the reversed electron energy band is too important then the inverse SSE disappears. This condition is quite restrictive for the appearance of the inversion of the SSE.

5. Conclusions-We have demonstrated an important influence of the spectrum anisotropy of F layers and energy shift between F and S energy spectra on the F/S heterostructures properties.

For the S/F/S sandwiches, we have shown that the π-phase cannot appear if the energy shift between S and F spectra is too important. Moreover, on the example of elliptic Fermi surfaces, we show that the anisotropy of F layer spectra may provide rather restrictive conditions of existence of the π-phase. We have also considered a more general S/F/S structure by including several F layers between the two S layers. The 0-π phase transition is influenced by the even or odd number of F layer independently of the relative magnetization between adjacent F layer.

For the F/S/F trilayers, we propose a possible explanation of the previously observed p-8 inverse SSE. In strong ferromagnet, the high energy band splitting may imply an inversion of one of the electron energy band. This type of electron spectrum in F layer could explain the anomalous SSE. We also study the influence of the effective mass difference and the energy shift on the inverse spectra.

Fig. 1 :

 1 Fig. 1: Geometry of the three-layer system. Here ψ 0 , ψ 1 and ψ 2 are the electron annihilation operators in the corresponding layer.

Fig. 2 :

 2 Fig. 2: 2-D Fermi surfaces for superconductor ξ s (solid line) and for majority ξ ↑ and minority ξ ↓ bands of ferromagnet (dashed and dotted lines correspondingly). In (a) is depicted the situation where ξ s < ξ ↓ < ξ ↑ . In (b) is depicted the situation where ξ ↓ < ξ s < ξ ↑ .In (c) are represented the different areas where the integration are done. In area (1) ξ s < ξ ↑ < ξ ↓ , in area (2) ξ ↑ < ξ s < ξ ↓ and in area (3) ξ ↑ < ξ ↓ < ξ s . (d) represents the sign of ln `T 0 c /T π c ´in the different areas.In the shaded areas, the integral is negative and in the other one, the integral is positive. In the frontier, the integral equals zero.

Fig. 3 :

 3 Fig. 3: (a) Spectrum of a superconducting metal (b) Spectrum of a normal metal where ξ ↑ = ξ ↓ = ξ s . (c) Spectrum of a strong ferromagnet with reversed spectrum where ξ ↑ = -ξ s and ξ ↓ = ξ s . (d) Spectrum of a strong ferromagnet with a reversed spectrum where ξ ↑ = -ξ s and ξ ↓ = ∆E ≫ kBTc.

Fig. 4 :

 4 Fig. 4: (a) ∆T /Tc as a function of ∆t/t0 is presented in the units of (t0/Tc) 4 in the situation ξ ↑ = -ξ s and ξ ↓ = ξ s . (b) ∆T /Tc as a function of ∆E/Tc presented in the units of (t/Tc) 4 with t ↑ = t ↓ = t, in the situation ξ ↑ = -ξ s and ξ ↓ = ∆E ≫ Tc. (c) ∆T /Tcas a function of α presented in the units of (t/Tc) 4 with t ↑ = t ↓ = t, ∆E = 10Tc and g ∆E = Tc in the situation ξ ↑ = -α.ξ s + g ∆E and ξ ↓ = ∆E ≫ Tc. (d) ∆T /Tc as a function of g ∆E/Tc presented in the units of (t/Tc) 4 with ∆E = 10Tc, t ↑ = t ↓ = t and α = 1 in the situation ξ ↑ = -α.ξ s + g ∆E and ξ ↓ = ∆E ≫ Tc.

Acknowledgements: The authors thank J. Cayssol for useful discussions and helpful comments. This work was supported by the French ANR Grant No. ANR-07-NANO-011:ELEC-EPR.


domain: spi math info phys

text:
J Floquet

Frédéric Champagnat
email: frederic.champagnat@onera.fr

Samir Beneddine

Physics-driven PTV interpolation through learned regularization

Keywords: PTV processing, Neural Networks, Unrolled optimization

This communication deals with spatial interpolation of instantaneous Particle Tracking Velocimetry (PTV) data. Such a task can be handled using generic signal processing tools like B-splines, radial basis functions (RBF) or Kriging.

We deal with non time-resolved PTV as gathered by dual frame particle image velocimetry (PIV) systems followed by dual frame PTV algorithms. In such a dual frame context, Physics-based interpolation techniques typically rely on enforcement of divergence free constraints on the interpolated field. The aim of this communication is to explore the wider field o f P hysics-based P TV i nterpolators b ased o n l earning a d eep n eural n etwork. T hese interpolators are learned on a database of flow e xamples t hat m ay be g enerated by h igh resolution P IV or D NS. We p ropose an instance of such a Physics-based PTV interpolator based on the machine learning concept of unrolled optimization. In this communication, this concept is extended in order to deal with random measurement location, characteristic of PTV data. We also show how to use this concept to perform pressure estimation from PTV data, along with velocity interpolation. The concept is illustrated on a 2D laminar flow, DNS data is used to train and benchmark the proposed method against RBF interpolation. The proposed method is shown to be much more robust than RBF interpolation against low density seeding and noise.

Introduction

This communication deals with spatial interpolation of instantaneous Particle Tracking Velocimetry (PTV) data. Such a task can be handled using generic signal processing tools like B-splines [START_REF] Gesemann | From noisy particle tracks to velocity, acceleration and pressure fields using b-splines and penalties[END_REF] or spatial statistic tools such as Kriging and radial basis functions (RBFs) [START_REF] Azijli | Solenoidal filtering of volumetric velocity measurements using gaussian process regression[END_REF], [START_REF] Vennell | A divergence-free spatial interpolator for large sparse velocity data sets[END_REF]. In the context of Data Assimilation, Physicsbased interpolation can be handled for time-resolved PTV along an horizon of ten times steps [START_REF] Schneiders | Dense velocity reconstruction from tomographic ptv with material derivatives[END_REF] [START_REF] Gesemann | From noisy particle tracks to velocity, acceleration and pressure fields using b-splines and penalties[END_REF] [START_REF] Yang | Kernelized lagrangian particle tracking[END_REF]. Here we deal with non time-resolved PTV as gathered by dual frame particle image velocimetry (PIV) systems followed by dual frame PTV algorithms [START_REF] Fuchs | Double-frame 3D-PTV using a tomographic predictor[END_REF][START_REF] Cornic | Double-frame tomographic ptv at high seeding densities[END_REF][START_REF] Mikheev | Enhanced particle-tracking velocimetry (EPTV) with a combined two-component pair-matching algorithm[END_REF][START_REF] Takehara | A kalman tracker for superresolution piv[END_REF]. In this dual frame context, Physics-based interpolation techniques merely rely on enforcement of divergence free constraints on the interpolated field [START_REF] Vennell | A divergence-free spatial interpolator for large sparse velocity data sets[END_REF] [START_REF] Azijli | Solenoidal filtering of volumetric velocity measurements using gaussian process regression[END_REF].

The aim of this communication is to explore the wider field of Physics-based PTV interpolators based on learning a deep neural network. These interpolators are based on a database of flow examples that may be generated by high resolution PIV or DNS. We propose an instance of such a Physics-based PTV interpolator based on the machine learning concept of unrolled optimization. In this communication, this concept is extended in order to deal with random measurements characteristic of PTV data. The concept is illustrated on a 2D laminar flow, DNS data is used to train and benchmark the proposed method against RBF interpolation. The proposed method is shown to be much more robust than RBF interpolation against low density seeding and noise.

PTV interpolation is by nature an ill-posed problem. Whereas RBF interpolation or Kriging can be interpreted as instances of regularization using hand-crafted penalty functions (e.g. Tikhonov penalties), we propose to learn this penalty function based on flow samples. Therefore, our method intends to provide a Physics-based interpolation without resorting to Data Assimilation. The proposed method is an instance of unrolled optimization methodology [START_REF] Gilton | Neumann networks for linear inverse problems in imaging[END_REF] and extends the concept to the non trivial case of random location of measurements, which is typical of PTV data.

It is well known that pressure can be inferred from velocity fields for compressible flows [START_REF] Van Oudheusden | PIV-based pressure measurement[END_REF]. Since our DNS yields the pressure fields in addition to velocities, we also use the pressure field to train our networks in order to show that the present approach can also infer pressure from PTV velocities, without explicitly solving the associated Poisson equation.

Proposed method

Unrolled optimization is a general framework to learn a regularizer for solving an inverse problem [START_REF] Gilton | Neumann networks for linear inverse problems in imaging[END_REF]. The following content is specialized to the recovery of the velocity field u from a partial measurement y. Let us introduce the measurement operator h such that

y = h(u) + n.
(1)

Depending on the experimental setup, h may be a filtering (PIV) or a sparse sampling (PTV) operator, n is a measurement noise. Operator h and noise statistics can be accurately calibrated at the beginning of the experimental phase, and are therefore considered as known (in a statistical sense) hereafter. The present study is restricted to linear operators h. This assumption is satisfied for classical PIV/PTV setup. A first approach for inversion of a measurement y is to define the best estimate u as the outcome of the following optimization problem

u = arg min u 1 2 ∥h(u) -y∥ 2 + r(u), ( 2
)
where r(u) is a penalty function required to stabilize the optimization problem, and is typically hand-crafted. The computation of u is classically given by an iterative minimization algorithm. In the simple case of a steepest gradient descent, it reads

u (k+1) = u (k) -η h * h( u (k) ) -y -η∇r( u (k) ), (3)
where η is a predetermined descending step, (•) * denotes the adjoint operator (which reduces for a linear operator, in the simplest case, to the transpose). Therefore, solving the optimization problem actually requires to specify ∇r rather than r. The present work proposes to define ∇r through supervised learning techniques.

Unlike gradient descent schemes that rely on a stopping rule based on a threshold on some residual, unrolled techniques specify a fixed number of iterations set to B. The overall Gradient Descent Network (GDN) then consists of B identical and sequential blocs. B is typically of the order of ten. Each block implements one iteration of (3). Last but not least, ∇r is implemented as a convolutional neural network ∇r(•) = R(θ, •) whose weights θ are tuned during a supervised learning phase: the training database is made of pairs (u, y) and the weights θ are calibrated through standard deep learning optimization techniques so that the output of the GDN minimizes a given error with respect to the database. The architecture of the overall GDN is depicted on Fig. 1. These ideas may be generalized to any descent-like optimization algorithm, variations about this template may be found in [START_REF] Adler | Learned primal-dual reconstruction[END_REF]; [START_REF] Aggarwal | MoDL: Model-based deep learning architecture for inverse problems[END_REF]; [START_REF] Chen | Trainable nonlinear reaction diffusion: A flexible framework for fast and effective image restoration[END_REF]; [START_REF] Gilton | Neumann networks for linear inverse problems in imaging[END_REF].

(•) -η h * h(•) -y -η∇r(•) = R(θ, •) u (0) + + (•) -η h * h(•) -y -η∇r(•) = R(θ, •) u (1) + + u (2) ... u = u (B)
Figure 1. Architecture of a GDN corresponding to direct unrolled optimization. The upper part of each block is specified by the measurement part of the criterion, the only parameter to adjust is η (which is a trainable parameter in the unrolled algorithm). The lower part of each block is made-up of a convolutional neural network that will learn the gradient of the penalty function.

Including the pressure in the reconstruction

Optical measurement techniques such as PIV/PTV do not yield the pressure fields, which has to be estimated from the velocity field. Small adjustments from the previously-introduced unrolled approach may provide a data-based estimate of the full flow state q = (u, p) t . The pressure p(x, y) may easily be added to the reconstruction procedure by replacing u by the augmented vector s = (u, p) t , and defining an augmented operator H such that H(s) = (f (u), 0) t . Then the new unrolled approach is based on the iterative relation

s (k+1) = s (k) -η H * y -H(s (k) ) -η∇r(s (k) ). (4)
Note that the operator H does not provide any physical or measurement information about the pressure. Therefore, the pressure estimation is entirely left for the term ∇r(•) to be accounted for: the neural network will have to learn from data how to model the pressure field. However, if pointwise pressure measurement were available, then they could be easily taken into account in H by changing its definition. In the following, we focus on the particular case where no pressure measurements are available and assess the performances of the unrolled approach for flow reconstruction from PTV-only measurements. [START_REF] Gilton | Neumann networks for inverse problems in imaging[END_REF] proposed a variant of the previous approach that is represented in figure 2.

Neumann Network variant for unrolled optimization

The general structure mainly differs from that of a classical unrolled approach by the presence of skip connections, represented in red. Mathematically, this new architecture is based on a Neumann series expansion of the operators appearing in the gradient descent iterations, details may be found in [START_REF] Gilton | Neumann networks for inverse problems in imaging[END_REF]. [START_REF] Gilton | Neumann networks for inverse problems in imaging[END_REF] showed that for image processing tasks, these so-called Neumann Network (NN) architectures yield better results compared to classical unrolled methods.

(•) -η H * • H(•) -η∇r(•) = R(θ, •) s (0) + + + + (•) -η H * • H(•) -η∇r(•) = R(θ, •) s (1) + + + + s (2) ... ... + + s (B) s = B i=0 s (i)

Modelling PTV measurement operators

A PTV measurement is a sparse and random measurement of the flow, sparse because it is attached to a particle that has been seeded in the flow, random, because there is no way to predict the position of that particle at a given instant.

One possibility for modelling such measurement is through the introduction of a mask. Presently, the flow is represented by a three-channel image. Selecting some location on the flow can be represented as multiplication with a predefined mask, as sketched in the following figure 3.

This is an approximation of the actual situation because there is no reason that all particle positions would actually match pixel positions. A more realistic situation could be modelled as well, but it would require more development without providing more information about the performances of the present approach. The goal of this work is to study the concept of sparse measurement with unrolled approaches, and the proposed model seems sufficient for this goal.

Since the operator H has sparse support due to local measurement, the linear part of the network produces corrections only at pixels of the mask. This produces small spikes in the output U and V channels, these spikes are located at pixels of the mask. This could be corrected in principle by using a large number of blocs B. Instead we relied on median filters that are very efficient in suppressing impulse noise [START_REF] Huang | A fast two-dimensional median filtering algorithm[END_REF]. In practice we found more efficient to keep a low number of blocks, but to add a final 3 × 3 median filter stage to the U and V channels output of the network.

Note that locations of flow measurements change from one PTV snapshot to another. Therefore, contrary to the classical unrolled approach from [START_REF] Gilton | Neumann networks for linear inverse problems in imaging[END_REF], the measurement operator H changes for each sample from the training database. This requires an original training strategy developed hereafter.

Learning strategy

The neural network has two kinds of inputs: the mask and the velocity measurements. It has also two kinds of outputs: the interpolated velocity field and the associated pressure field. The classical unrolled strategy deals with fixed H operator during training and testing. As mentioned earlier, this strategy is useless for PTV applications where measurement locations are random. Thus we consider random masks that change during the learning process in order to span the variety of measurement locations.

The training Data Base (DB) is made-up of T snapshots (U, V, P ). Each snapshot is associated to a random mask with a prescribed number of non zero values (this number amounts to the Nppp parameter that characterises the seeding density in PTV). The learning procedure proved to be more efficient when the training quantities have similar range. In order to achieve this we remove a single constant value to the pressure field before training. Therefore, pressure is presently recovered up to a constant value. Such a constant could be recovered using a side neural network, a task we leave here to future work.

An epoch of the learning process consists in feeding the network with data generated by the sample and the mask then comparing the network output with the original sample considering the three components (U, V, P ) of the flow. The loss function to minimize is the Mean Squared Error defined as follows:

MSE(X, Y ) = 1 C • M • N C,M,N c,m,n=1 (X[c, m, n] -Y [c, m, n]) 2 , ( 5
)
with C representing the number of output channels (in our case C = 3 for the three flow variable (U, V, P )); M, N denote the horizontal and vertical width of the network output.

Experimental setting

Test case

The database used for training, validation and testing has been generated using the finite-volume compressible solver elsA [START_REF] Cambier | The onera elsa cfd software: input from research and feedback from industry[END_REF]. The physical configuration corresponds to a bidimensional laminar wake behind a bluff body. The computational domain as well as the mesh and boundary conditions may be seen in figure 4. The mesh is made of approximately 50000 points, which are distributed in order to properly resolve the regions with strong velocity gradient.

The computation is carried out using non-dimensional quantities: the free-stream velocity and static pressure is set to 1, the Mach number is set to M = 0.15 such that compressibility effects (which scale with M 2 ) are negligible and the flow may be treated as incompressible. The only parameter that is varied is the Reynolds number (28 values considered between 200 and 1050) through the value of the viscosity. The spatial scheme is a Jameson scheme and the temporal integration is performed using the Gear algorithm.

Snapshots are extracted at a given sampling frequency, that is set to have at least 80 snapshots per period of vortex shedding (which depends on the considered Reynolds number). For each Reynolds number considered, the computed snapshots are interpolated on a regular grid spanning the domain (x, y) ∈ [0, 8] × [-2, 2] and made of 70 × 41 points. The interpolation is performed through a 3 rd order bi-variate spline using the Python package sciPy. An example of resulting fluctuation field may be seen in figure 5. The final database contains a total of 4200 tuples of 3 images associated with (U, V, P ).

Network architecture, learning and hyperperameter tuning

We consider hereafter either the Gradient Descent Network of figure 1 or the Neumann Network of figure 2.

The quantities to be learned are the regularizer R(θ, •), and the step size η. Following Gilton et al.

(2020) R(θ, •) is a 7-layer convolutional neural network with an auto-encoder structure. the activation functions are ELU, the filters are square with sizes (5, 3, 2, 3, 3, 5, 3), and strides (1, 2, 1, 2, 1, 1, 1).

The number of channels in each layer is (3,64,256,256,256,64,3).

All the presented results where performed with 32 × 32 images on 3 channels (U, V, P ) as inputs and outputs, plus 3 mask channels as input. The P input channel is always set to zero, since no pressure data is available.

The database (DB) used for training a model is performed on 32 × 32 images cropped randomly from the eslA DB, see figure 5.

Several databases of masks (MDB) were also generated with 4200 tuples of three 32 × 32 binary images coupled with DB. Each tuple have two identical binary matrix generated randomly and a third images with zero for all elements. Each MDB is characterised by an average percent of non-zero values, which amounts to Nppp in a PTV experiment.

The DB and MDB are split in 80% for training and 20% for validation, so 3360 couples of tuples for one and 840 for the second. This step enabled to give a unique mask for each tuple of DB and avoids overfitting for a particular mask.

The network has been trained with an ADAM optimizer with initial step equal to 5.010 -3 and 0.99 Figure 6 shows a typical learning curve in terms of MSE and Peak SNR (PSNR) defined as

P SN R(X, Y ) = -10 log MSE(X, Y ).
On figure 6, the MSE decreases fast on the first epochs, then this decrease slows down and is better seen on the PSNR logarithmic metric. After about 200 epochs, the validation loss starts to increase, yielding an early stopping (vertical red dashed line in figure 6).

Reference RBF method

A RBF method has been selected as a competitor to the NN or GDN. Interpolation is then performed independently on U and V channels. The best performing kernel for the DB is the inverse kernel:

K ϵ (r) = 1 √ r 2 + ϵ 2 , ( 6
)
where ϵ is a regularizing parameter tuned as the average distance between interpolating nodes. The method was implemented using the RBF interpolation package of SciPy.

Experimental results

We assess hereafter the proposed methods (GDN, NN) and RBF reference on the proposed test case. These comparisons are performed using the error RMS on the velocity magnitude as main performance metric [START_REF] Sciacchitano | Main results of the first Lagrangian particle tracking challenge[END_REF]. Our first experiment assess GDN vs. NN. GDN is shown to perform best on a variety of seeding density and noise levels. Then we follow on by comparing GDN and RBF both on velocity and pressure estimation. The generalization capability of GDN is explored through decoupling of learning DB and testing DB.

Comparison of GDN and NN unrolling strategies

We trained the GDN and NN architectures for four N ppp values: 0.02, 0.05, 0.1, 0.2. During the test, a random 32 × 32 crop is performed in each of the basic eslA DB, then a random mask is generated in order to produce pseudo PTV measurements. A Gaussian random noise of standard deviation σ is added (since the dynamic of flow is about 1 in these experiments, σ = 0.01 amounts to 1% noise). Random crop and random masks are performed on-the-flight during the test, therefore generating the same pattern of measurement on the same sample as during training is very unlikely.

Adapting from [START_REF] Sciacchitano | Main results of the first Lagrangian particle tracking challenge[END_REF], the main performance metric is the error RMS on the velocity magnitude. Indeed, in the present case, this is almost equivalent to random error in Sciacchitano et al. ( 2021) since the bias is negligible.

The evolution of the RMS metric is shown on figure 7 with respect to seeding density and for two noise configurations (σ = 0.0 and σ = 0.05): GDN always performs consistently better than NN. We observed that the Neumann series associated to sparse measurement operators converges very slowly, compared to gradient iteration counterparts. Note also that the RMS increases for NN with no noise for the larger Nppp. The behaviour of NN could be improved using preconditionning as suggested by [START_REF] Gilton | Neumann networks for linear inverse problems in imaging[END_REF]. We leave this conjecture to further work and concentrate on GDN in the following.

Comparison of GDN with respect to RBF reference

The GDN obtained in section 4.1 for four N ppp values: 0.02, 0.05, 0.1, 0.2 are tested and compared to the RBF approach, the RMS of the velocity magnitude with respect to seeding density is represented on figure 8.

As expected all quantities decrease as Nppp grows. It can be seen that the GDN approach is much more robust to addition of noise to PTV measurements. Adding 5% noise implies a slight increase in RMS for the GDN method, whereas it significantly degrades the RMS performance of RBF, in particular for the largest density, where the RMS increases by 50%. Moreover, the RMS performance of GDN interpolation is about 70% lower than its RBF competitor at the lowest density, whatever the noise level.

Qualitative results

This section shows a sample of the output of the GDN approach compared to the RBF method, with N ppp = 0.02 and no added noise. The first column of figure 9 provides the ground truth for (U, V, P ), and black circles mark the location of velocity data used by the GDN and RBF to recover the full field. Note that there is no circle on the pressure ground truth since there is no pressure measurement. The second row of figure 9 provides the output of the GDN, and the third row provides the RBF result. Since there was no input in terms of pressure, no output can be delivered by the RBF approach.

Figure 9 shows that the overall structure of (U, V, P ) estimated by the GDN is much more faithful to the original field than the RBF interpolation. Note that the V field is particularly ill-recovered by the RBF method. One drawback of RBF at low density is clearly apparent here: due to random measurement location there is no measurement of the V field around its lowest values. The driving principle of RBFs is to spread measurements around their location consistently with respect to each other measurement, this principle fails here because there is no V measurement in the central part of the flow where two local minima occur. It is likely that in the present case, the prior learned by the Neumann Network can compensate the lack of data on V by using the correlation between U and V channels. Note that indeed inter-channel correlation enabled the main features of the P field to be recovered despite the lack of pressure data. The performance of pressure reconstruction is assessed in more detail in the next section.

Pressure from PTV measurements

Here we reuse the unrolled architectures learned in section 4.1 and focus on assessment of pressure estimation.

Figure 10 shows the RMS error for the pressure estimation with respect to the seeding density for the NN and GDN methods. The behaviour is consistent with the velocity estimation case of figure 7, indeed the RMS error decreases with Nppp, except for the NN with no noise. In the latter case, RMS increases slightly, consistently with the behaviour observed in figure 7. We believe this behaviour is linked to the slow convergence of the Neumann series associated to the measurement operator. GDN always perform better than NN also for pressure RMS. Noise impact is most noticeable for the lowest Nppp, where the RMS of GDN increases by 25%. During test, a random crop is performed in the selected flow image and a random mask is generated on the flight. The RMS of the velocity magnitude with respect to Re groups is shown in figure 11.

Generalisation capability of unrolled optimization

First, we observe that the RMS of GDN is about twice as large as that of figure 8, but note that the training DB here is about half as small as the one used in figure 8. One can observe a slight increase of RMS with respect to Re for all methods. The RMS of GDN stays about 40% below that of RBF method, whatever the noise level. So GDN extrapolates fairly well in the Reynolds number range between 625 to 1050.

Conclusions

We proposed an instance of Physics-based PTV interpolator based on the machine learning concept of unrolled optimization. We extended the concept in order to deal with random measurements characteristic of PTV data. We also show that such neural networks can be trained to infer pressure from velocity only measurements. The proposed method is illustrated on a 2D laminar flow, DNS data is used to train and benchmark the proposed method against RBF interpolation. The proposed method is shown to be much more robust than RBF interpolation against low seeding density and noise. Two variants of unrolled optimization have been tested, the basic Gradient Descent Network and the Neumann Network variant. The latter proved to be less accurate for recovering velocity magnitude and pressure. The Gradient Descent Network variant has also passed a preliminary test of generalization capability, where testing and learning database are made of simulations with non overlapping Reynolds number range.

This work can be completed and extended in many directions. First, the present approach recovers pressure up to a constant additive factor, that is fixed in the present study. This is because the different (U, V, P ) channels of the network should have similar dynamic range. The remaining unknown constant could be estimated separately and added to the actual output of the NN or GDN. Second, the present implementation performs interpolation on 32 × 32 windows, the size of which can be augmented in relation to the memory available on training hardware. However, this augmentation will not reach Mega pixel images. Following (Gilton et al., 2019a), we propose to split such large images into overlapping patches (of size to be defined) each patch is then processed in parallel by the network and the output patches are then fused in order to yield a full image

Figure 2 .

 2 Figure 2. Architecture of a Neumann Network. At the input of the NN s (0) = ηH * y.

Figure 3 .

 3 Figure 3. PTV measurement model: the sparse velocity field on the left may be considered as the product of a full velocity field and a zero-one mask.

Figure 4 .

 4 Figure 4. Mesh used for the computations. Boundary conditions are shown in color: uniform inflow conditions are in blue, adiabatic no-slip walls in red, and non-reflective farfield conditions in green.

Figure 5 .

 5 Figure 5. Snapshot of V component of the flow at Re = 1000, from the eslA database. An instance of a 32 × 32 pix.random crop is depicted with a white dashed line square.

Figure 6 .

 6 Figure 6. MSE and PSNR of training and validation in function number of epoch until the early stopping (red line).The number of non-zero values per image in the DBM is equal to 2% (N ppp = 0.02).

Figure 7 .

 7 Figure 7. RMS of the velocity magnitude with respect to seeding density. For NN and GDN methods and two noise levels σ = 0.0 and σ = 0.05. GDN allways performs better than NN.

Figure 8 .

 8 Figure 8. RMS of the velocity magnitude with respect to seeding density. For GDN and RBF methods and two noise levels σ = 0.0 and σ = 0.05.

Figure 9 .

 9 Figure 9. Sample of (U, V, P ) fields recovered by Gradient Descent Network and Radial Basis function methods with noise level σ = 0.00. Black circles mark the location of velocity measurements. The central part of the V component of the flow is erroneous for RBF due to missing measurements in that area that contains to large basins.

Figure 10 .

 10 Figure 10. RMS of pressure estimation with respect to seeding density, for NN and GDN for two noise levels σ = 0.00, 0.05.

20th LISBON Laser Symposium 2022 Figure 11 .

 202211 Figure 11. RMS of the velocity magnitude with respect to Re. For GDN and RBF methods and two noise levels σ = 0.0 and σ = 0.05. Seeding density is set to Nppp = 0.02.

  So far, training and test has been performed on DB covering the DNS flows for all the available Reynolds numbers. We now explore the generalisation capability: training and test DB will be based on non overlapping Reynolds numbers range. The training DB will be made on all simulations with Re number between 200 and 575. During training the seeding density is fixed to Nppp= 0.02. All other training parameters are kept as in previous experiments. the testing DB is made of all simulations with Reynolds number between 625 to 1050. This testing DB is split in four groups labelled 650, 750, 850 and 1000 gathering neighbouring Reynolds number simulations.

result. Third, although the training database fulfils the divergence free constraint, this is not strictly the case for output velocities. The current approach can be reformulated in order to impose the divergence free constraint on output velocities.