Anne J E Jimenez-Lopez

Nadine Leconte

Fabienne Lambrouin

F Garnier-Lambrouin

A Bouchoux

F Rousseau

G Gésan-Guiziou
email: genevieve.gesan-guiziou@rennes.inra.fr

Ionic strength dependence of skimmed milk microfiltration: Relations between filtration performance, deposit layer characteristics and colloidal properties of casein micelles

Keywords: Microfiltration, skimmed milk, ionic strength, fouling, colloidal interactions

published or not. The documents may come

Ionic strength dependence of skimmed milk microfiltration: Relations between filtration performance, deposit layer characteristics and colloidal properties of casein micelles

Introduction

Skimmed milk is a complex mixture classically described as a suspension of casein micelles (association of αs1, αs2, β, and k-caseins with an average radius of 100 nm, 26 g L -1 ) in an aqueous phase containing soluble proteins (mostly b-lactoglobulin and a-lactalbumin, with an average radius of ~5 nm, 7 g L -1 ), minerals (mostly calcium and phosphate, 8 g L -1 ), and lactose (50 g L -1 ). Casein micelles are negatively charged colloids in balance with the aqueous phase. Physico-chemical modifications of milk, such as changes in ionic strength which is classically performed in the dairy industry, induce direct changes in casein micelles characteristics in terms of average diameter, charge, hydration, mineral and protein composition [START_REF] Gaucheron | The minerals of milk[END_REF]. This particular behavior makes it difficult to control the various operations used in dairy processing.

Crossflow microfiltration (0.1 µm) of skimmed milk is largely used in dairy industry to separate casein micelles from soluble proteins in order to produce a retentate enriched in casein micelles and a permeate enriched in soluble proteins, minerals and lactose [START_REF] Brans | Membrane fractionation of milk: state of the art and challenges[END_REF][START_REF] Maubois | New applications of membrane technology in the dairy-industry[END_REF]. The overall performance of the process, in terms of permeability and selectivity, is governed by the accumulation of the casein micelles at the membrane surface [START_REF] Gésan-Guiziou | Critical stability conditions in crossflow microfiltration of skimmed milk: transition to irreversible deposition[END_REF][START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF][START_REF] Berre | Skimmilk crossflow microfiltration performance versus permeation flux to wall shear stress ratio[END_REF]. As in colloid filtration, the matter accumulated in the concentration polarization layer turns into an irreversible deposit (= a gel) under given operating conditions commonly named "critical hydrodynamic conditions" (CHC) [START_REF] Gésan-Guiziou | Critical stability conditions in crossflow microfiltration of skimmed milk: transition to irreversible deposition[END_REF][START_REF] Bacchin | Modelling of filtration: from the polarised layer to deposit formation and compaction[END_REF][START_REF] Bacchin | A unifying model for concentration polarization, gel-layer formation and particle deposition in cross-flow membrane filtration of colloidal suspensions[END_REF][START_REF] Pignon | Structure and rheological behavior of casein micelle suspensions during ultrafiltration process[END_REF][START_REF] Chen | Particle deposition during membrane filtration of colloids: Transition between concentration polarization and cake formation[END_REF]. In terms of performance, the CHC correspond to the transition between steady runs, i.e., high transmission of soluble proteins and low fouling resistance, and divergent runs, i.e., sharp decrease in transmission of soluble protein with high and irreversible fouling resistance [START_REF] Gésan-Guiziou | Critical stability conditions in crossflow microfiltration of skimmed milk: transition to irreversible deposition[END_REF][START_REF] Metsamuuronen | Critical flux in ultrafiltration of myoglobin and baker's yeast[END_REF]. The CHC can be defined through the critical ratio (J/τw)crit, corresponding to the ratio between convective mass transport towards the membrane, related to the permeation flux J, and the so-called "back-transport" from the membrane to the bulk, related to the wall shear stress τw [START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF][START_REF] Gésan-Guiziou | Wall shear stress: effective parameter for the characterisation of the cross-flow transport in turbulent regime during skimmed milk microfiltration[END_REF].

To control the performance of skimmed milk microfiltration, it is necessary to determine the mechanisms behind the CHC and so to understand how the casein micelles accumulate at the membrane and how they interact with each other until forming a gel. To explore these questions, one approach consists in modifying the physico-chemical characteristics of the casein micelles and investigating the changes in the filtration performance. To this respect, interesting results have already been obtained by changing the composition of the milk and performing filtrations [START_REF] Kühnl | Impact of colloidal interactions on the flux in cross-flow microfiltration of milk at different pH values: A surface energy approach[END_REF][START_REF] Berre | Microfiltration (0,1 µm) of milk: effect of protein size and charge[END_REF][START_REF] Rabiller-Baudry | Limiting flux in skimmed milk ultrafiltration: impact of electrostatic repulsion due to casein micelles[END_REF][START_REF] Youravong | The effects of physico-chemical changes on critical flux of skimmed microfiltration[END_REF]. For instance, the addition of citrate to milk increases dramatically the total and irreversible fouling resistance in filtration. This is presumably because citrate induces a decrease in the casein micelles diameter and consequently leads to a deposit of lower porosity and membrane pore blocking [START_REF] Berre | Microfiltration (0,1 µm) of milk: effect of protein size and charge[END_REF][START_REF] Youravong | The effects of physico-chemical changes on critical flux of skimmed microfiltration[END_REF]. Decreasing the milk pH down to 5.5-5.6 or adding high concentrations of NaCl up to 0.4-0.7 M is also known to decrease the "limiting" and "critical" flux during skimmed milk filtration [START_REF] Rabiller-Baudry | Limiting flux in skimmed milk ultrafiltration: impact of electrostatic repulsion due to casein micelles[END_REF][START_REF] Youravong | The effects of physico-chemical changes on critical flux of skimmed microfiltration[END_REF] and to increase the fouling resistance remaining after filtration [START_REF] Kühnl | Impact of colloidal interactions on the flux in cross-flow microfiltration of milk at different pH values: A surface energy approach[END_REF][START_REF] Berre | Microfiltration (0,1 µm) of milk: effect of protein size and charge[END_REF][START_REF] Youravong | The effects of physico-chemical changes on critical flux of skimmed microfiltration[END_REF]. These results are mainly explained through the decrease in the average surface charge of the casein micelles (measured through the casein micelles zeta potential), leading to lower electrostatic repulsion forces between them. More recently, Kühnl et al., [START_REF] Kühnl | Impact of colloidal interactions on the flux in cross-flow microfiltration of milk at different pH values: A surface energy approach[END_REF] propose an explanation based on a shift in the balance between hydrophilic repulsion and hydrophobic attraction between the micelles as the pH is lowered.

In our opinion, the studies mentioned previously still fail in giving clear and unambiguous relationships between the physical properties of the casein micelles and the performance of milk filtration. A first important reason is that the characteristics of the deposit layer (cohesion, dynamic behavior during filtration, transmission properties) are poorly investigated in these works. Such information is crucial if one wants to accurately describe the mechanisms involved during filtration. A second reason is that the nature of the dispersions used in the quoted studies as well as the control of the performed physico-chemical modifications, are often questionable. Indeed, the solutions used are either strongly heat-treated milks (Ultra High Temperature milks [START_REF] Rabiller-Baudry | Limiting flux in skimmed milk ultrafiltration: impact of electrostatic repulsion due to casein micelles[END_REF]) or milk powder suspended in water [START_REF] Youravong | The effects of physico-chemical changes on critical flux of skimmed microfiltration[END_REF], i.e., solutions in which the casein micelles are not in their native solvent and/or are not strictly in their "native state". Moreover, in the studies in which the effects of NaCl or citrate are investigated [START_REF] Berre | Microfiltration (0,1 µm) of milk: effect of protein size and charge[END_REF][START_REF] Youravong | The effects of physico-chemical changes on critical flux of skimmed microfiltration[END_REF], the variation of pH resulting from the addition of salts is generally not corrected. For example, the addition of 40 mM of citrate increases the pH from 6.6 to 7.1 which modifies the milk particle size distribution [START_REF] Youravong | The effects of physico-chemical changes on critical flux of skimmed microfiltration[END_REF]. Also the addition of respectively 700 mM and 920 mM of NaCl to milk leads to a decrease in pH of 0.3 and 0.35 respectively and then results in partial solubilization of colloidal phosphate [START_REF] Rabiller-Baudry | Limiting flux in skimmed milk ultrafiltration: impact of electrostatic repulsion due to casein micelles[END_REF][START_REF] Youravong | The effects of physico-chemical changes on critical flux of skimmed microfiltration[END_REF]. Then both effects of pH and salt addition on the characteristics of casein micelles cannot be dissociated, which makes uncertain some of the reported interpretations. This paper focuses on the effect of ionic strength (through NaCl addition) on the performance of skimmed milk microfiltration. In that aim, fresh skimmed milk is used and ionic strength modifications are done with systematic pH adjustment. The properties of the casein micelles in terms of charge, average diameter and calcium content are measured as a function of the ionic strength. Filtrations of unmodified and modified milks are then performed and the properties of the casein micelles deposits are determined (setting up, structuration, transmission properties and fouling resistances). Some additional elements of understanding on the properties of the gel layer are also obtained through a careful treatment of the transmembrane pressures get during the CHC. The results are discussed with respect to the generic colloidal behavior of the casein micelle in concentrated regime as investigated through osmotic stress and rheological measurements [START_REF] Bouchoux | Rheology and phase behavior of dense casein micelle dispersions[END_REF][START_REF] Bouchoux | Casein Micelle Dispersions under Osmotic Stress[END_REF]. Hence this work tends to connect and explain the relations between the physicochemical properties of the casein micelles and the microfiltration performance through a new and original approach.

Experimental

Fluids

Dairy fluids

The present work is performed with heat-treated (68 ± 1 °C, 25 s) skimmed milk obtained from raw milk provided by Entremont Alliance (Montauban de Bretagne, 35, France). After fresh milk is collected (one day before filtration), skimming (50 °C, 315 L h -1 , Elecrem type 125, Châtillon, 92, France) and thermal treatment (Actini, Evian les bains, 74, France) are performed at UMR1253. To prevent micro-organism growth, 0.2 g L -1 of sodium azide is added to the milk.

The ionic strength is increased by direct addition of NaCl (0 -0.3 M, Fisher Scientific) at room temperature, under strong stirring. During the procedure, the pH variation is simultaneously corrected through the addition of small amounts of 2 M NaOH solution (the induced milk dilution is negligible with a dilution factor < 6 10 -3 ). The final ionic strength, after addition of 0, 0.085, 0.175 and 0.300 M of NaCl, is estimated to be 0.072, 0.157, 0.247 and 0.372 M in equivalent NaCl concentration (the estimation of ionic strength in equivalent NaCl has been checked through conductivity measurements [START_REF] Jimenez-Lopez | Structuration du dépôt de matière à la membrane lors de la microfiltration tangentielle (0,1 µm) de lait écrémé[END_REF]). All these physico-chemical modifications are performed the day before the microfiltration in order to reach the equilibrium state between the micellar and aqueous phases. The modified milk is then stirred one night at room temperature and the following day, the milk is warmed at 48°C (classic temperature used in the dairy industry for milk filtration) and then stirred for 1 hour before filling the microfiltration loop.

In all this work, unmodified milk and modified milks have comparable protein composition with a concentration in casein micelles of 26 ± 1 g L -1 and a concentration in total soluble proteins of 6.4 ± 0.3 g L -1 .

Water and cleaning solutions

The water used to warm up the installation to 48 °C before the experiment and to rinse the filtration rig after filtration is filtered on serial cartridges of 1.0 and 0.2 µm. The rinsing phase is performed in crossflow mode according to the method previously described by Jimenez-Lopez et al. [START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF].

Cleaning is performed at 50°C according to a procedure previously described [START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF]. Two detergent solutions are used: Ultrasil 25 F 1 % (v/v) (alkaline solution, Henkel-Ecolab SNC, Issy les Moulineaux, 92, France) and nitric acid (acid solution, HNO3; 58% purity) at 1% (v/v) (Langlois Chimie, St Jacques de la Lande, 35, France). Permeation flux is systematically recovered (± 5 %) after cleaning.

Microfiltration set-up 2.2.1. Membrane

The membrane used is a 0.1 µm tubular ceramic Kerasep TM membrane with 7 channels (NOVASEP PROCESS, 01, France) and an active layer of ZrO2 + TiO2: membrane area = 0.0744 m 2 , inner channel diameter = 4.13 ± 0.02 mm, length = 0.856 m. The hydraulic resistance of the clean membrane is 0.14 ± 0.02 ´ 10 12 m -1 , as determined through Darcy's law.

Filtration rig

The experimental rig, previously described [START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF], is fitted out with a uniform transmembrane pressure system (circulation of the permeate co-current to the retentate) in order to avoid any transmembrane pressure (TMP) difference along the hydraulic path [START_REF] Sandblöm | Filtering process[END_REF]. This allows performing experiments at high τw and low and homogeneous TMP, simultaneously.

Real time data are monitored and stored by a computer through a data acquisition card.

The control software is specially developed with LabVIEW 7.1 (National Instruments, USA) for this filtration rig.

Operating procedure

The conditions of the deposit set-up and its characteristics are determined following a methodology previously described [START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF][START_REF] Gésan-Guiziou | Wall shear stress: effective parameter for the characterisation of the cross-flow transport in turbulent regime during skimmed milk microfiltration[END_REF]: at constant permeation flux (J) and constant average retentate pressure (PAR), the cross-flow velocity is decreased by steps and increased again by the same steps [START_REF] Berre | Skimmilk crossflow microfiltration performance versus permeation flux to wall shear stress ratio[END_REF], thus performing successive filtrations at controlled wall shear stresses (tw). In these conditions, the TMP becomes a variable, directly related to the fouling phenomenon.

The experiments are performed at 48 °C. After concentration of the casein micelles at a volume reduction ratio (VRR) of 2, the parametric study (= run) is performed. The tw is decreased by steps of 30 minutes until reaching the critical hydrodynamic conditions (CHC, tw = twcrit) at which a sharp increase in TMP is observed [START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF]. After ~1 hour at τwcrit, τw is then increased step-by-step in a similar way and using the same values of tw that at decreasing tw.

From 120 to 80 Pa, τw is decreased by increments of 20 Pa at 6 Pa min -1 . From τw = 80 Pa to τw = τwcrit, in order to be more accurate in the determination of τwcrit, τw is decreased by increments of 5 Pa with a slope of 1 Pa min -1 .

The results reported are those obtained for one run with one given membrane. The experiments were repeated and similar results were obtained for the same given conditions (not shown). From the results obtained, the differences in TMP and protein transmission before and after the CHC are evaluated (what we call the "hysteresis" phenomenon in the rest of this work).

Physico-chemical analyses 2.4.1. Protein contents

The Kjeldahl method is used to determine the concentrations of casein micelles and total soluble proteins in milks, retentates and permeates after adapted preparations of samples.

Casein micelles concentration is calculated after determination of total nitrogen (TN) on raw samples and non-casein nitrogen (NCN) after acid precipitation of caseins to their isoelectric point at pH 4.1-4.6 (casein concentration = TN-NCN). Total soluble proteins are obtained by difference between non-casein nitrogen (NCN) and non-protein nitrogen (NPN) measured after the total protein precipitation with a trichloroacetic acid solution (15 % w/v). The serum proteins concentration obtained with the Kjeldahl method is compared to a chromatographic method of determination of the 5 major soluble milk proteins (RP-HPLC), adapted from Resmini et al. [START_REF] Resmini | Determinazione delle sieroproteine solubili del latte per HPLC (cromatographia liquida ad alta prestazione) in fase inversa[END_REF]. The RP-HPLC is performed using a PLRP-S column (gel of polystyrene divinylbenzene, 300 Å, 8 µm, 150×4.6 mm, Polymer Laboratories, Touzart et Matignon, 22, France). The sum of the concentrations of the 5 proteins obtained through HPLC and Kjeldahl gives similar results with a slight discrepancy (5 -8 %).

Calcium content

The determination of total and soluble calcium (Catotal, Cas) is performed using an atomic absorption spectrophotometer (AA300, Varian France, Les Ulis). The total calcium corresponds to the calcium contained in the casein micelles and in the aqueous phase, whereas the soluble calcium corresponds to the calcium contained in the aqueous phase only.

Soluble calcium is measured in the permeate obtained after ultrafiltration of the milks at 48 °C with a laboratory pilot (Filtron GP 141, Pall Gelman Laboratoire, St Germain en Laye, 92, France) using a 5 kg mol -1 polyethersulfone plate membrane (Pall Gelman Laboratoire, St Germain en Laye, 92, France).

Casein micelles characteristics

The average diameter and the electrophoretic mobility (global charge) of the casein micelles are determined in milk at 48 °C prior to filtration. The measurements are performed using the Zetasizer 3000 HS (Malvern, UK) after dilution of the casein micelles fractions in their corresponding aqueous phase obtained through ultrafiltration (5 kg mol -1 ) with the laboratory pilot described previously.

The turbidity of the microfiltration permeates is measured using a Turbidimeter Hach (2100AN IS, Hach, Namur, Belgium). This turbidity is attributed to a slight leak of casein micelles though the membrane, due to both distributions of membrane pore and micelle sizes [START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF]. In order to compare the turbidity data, and because there is a slight evolution in the membrane permeability with time, the turbidities are normalized with respect to that of the initial permeate turbidity obtained at τw = 120 Pa.

Calculations

Wall shear stress and hydraulic resistances

The wall shear stress, τw, is calculated from the inner diameter of the membrane channels as proposed by Le Berre and Daufin [START_REF] Berre | Skimmilk crossflow microfiltration performance versus permeation flux to wall shear stress ratio[END_REF].

The hydraulic resistances considered in these work are:

-Rm, the resistance of the cleaned membrane corrected by the irreversible fouling resistance which is assumed to be mainly attributed to protein adsorption phenomena that take place during the concentration phase. Rm is estimated from the water flux obtained at the end of an experiment performed in conditions where no deposit was formed. We find Rm = 0.51 10 12 m -1 and we assume this value is the same for all our experiments.

-Rp and Rf the polarization layer and the total fouling resistances, respectively. They are calculated from the Darcy's law at τw = 120 Pa, at the beginning and at the end of the parametric study, respectively. Thus Rp corresponds to the polarization resistance in given hydrodynamic conditions in which there is no deposit at the membrane surface. Conversely Rf corresponds to the resistance of the polarization layer (Rp) plus the resistance of the deposit (Rd) built during the run.

-Rd is the resistance of the deposit remaining at the membrane surface at the end of the parametric study. It is linked to Rp and Rd through:

Rd = Rf -Rp (1)
-Rid is the irreversible part of the deposit as determined through Darcy's law [START_REF] Berre | Skimmilk crossflow microfiltration performance versus permeation flux to wall shear stress ratio[END_REF] after a water flux at τw = 120 Pa at the end of the experiment.

-Rrd corresponds to the part of the deposit which is eliminated by water rinse after the experiment. It is calculated through:

Rrd = Rd -Rid (2)

Soluble protein transmission (Tr)

The transmission of soluble proteins, Tr, is calculated according to:

(3)

with CP and CR (g L -1 ), the protein concentrations in the permeate and in the retentate, respectively. The calculation of CP is corrected by the dilution factor of the permeate loop, assuming a perfectly stirred reactor. CP and CR are smoothed using the 5 points adjacent averaging method implemented in Origin ® 7.5 (OriginLab Corporation, USA). Tri and Trf are the average soluble protein transmissions for τw=120 Pa at the start and the end of the parametric study, respectively (Table 2). At the CHC (τw = τwcrit), sharp TMP oscillations are observed [START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF], either during the whole period of this "critical" step or during a fraction of this step only. Those variations in TMP are first analyzed according to the two following parameters:

Characterization of TMP evolution under critical conditions

-the overall increase in TMP during the oscillation phase, which is characterized by the slope of the TMP signal versus time after linearization (SlTMP in Pa s -1 ) (Fig 1 and4).

-the ratio , which gives the duration of the TMP oscillation phase in comparison with the total duration of the critical step (Dtcrit) (Fig 1).

We also characterize the TMP oscillations in terms of total number, frequency, and amplitude. Only the TMP oscillations with amplitudes that exceeded 0.5 ´ 10 5 Pa are considered in our calculations (all made using Origin ® 7.5).

-the number of oscillations (Nosc) is calculated from the number of maxima and minima (N) in TMP values (Fig 1), according to:

Nosc = (4)
-the TMP oscillations frequency (Fr, in s -1 ) is calculated according to:

Fr = (5)
with Δtosc (s) the duration of the TMP oscillation phase (Fig 3).

-the amplitude of the TMP oscillations (Amposc in Pa, Fig 1) is calculated according to:

Amposc = (6)
Assuming that the TMP oscillations are symmetrical over and under the plot of the linear

Results

Casein micelles characteristics versus ionic strength (I) at 48°C

The increase in ionic strength leads to changes in the characteristics of the casein micelles that are well documented at room temperature: at 20-25°C, the addition of NaCl to milk leads to the decrease in the global charge of the casein micelle, and to the partial release of the colloidal calcium into the aqueous phase [START_REF] Gaucheron | The minerals of milk[END_REF][START_REF] Dalgleish | Measurement of Electrophoretic Mobilities and Zeta-Potentials of Particles from Milk Using Laser Doppler Electrophoresis[END_REF][START_REF] Huppertz | Effect of NaCl on some physico-chemical properties of concentrated bovine milk[END_REF][START_REF] Graët | Effects of pH and ionic-strength on distribution of mineral salts in milk[END_REF][START_REF] Mekmene | A model for predicting salt equilibria in milk and mineral-enriched milks[END_REF]. However, to our knowledge, the modifications of the casein micelles characteristics in milk are still poorly known in our experimental conditions, i.e., 48°C and NaCl addition with simultaneous pH correction. It is then necessary to characterize the casein micelle modifications in such conditions to make sure that no unexpected phenomenon is induced. The results are presented in Table 1. With the increase in ionic strength (I) from 0.072 to 0.372 M, the casein micelle electrophoretic mobility (which is related to its global electric charge) linearly decreases from 31±1 to 15±1 10 -7 m 2 s -1 V -1 , while the concentration in soluble calcium (Cas) and the ratio Cas/Catotal linearly increase. Regardless of I, the average diameter of the casein micelle is constant at 181 ± 6 nm (Table 1). These observed trends and values are pretty similar to those obtained at room temperature [START_REF] Youravong | The effects of physico-chemical changes on critical flux of skimmed microfiltration[END_REF][START_REF] Dalgleish | Measurement of Electrophoretic Mobilities and Zeta-Potentials of Particles from Milk Using Laser Doppler Electrophoresis[END_REF][START_REF] Graët | Effects of pH and ionic-strength on distribution of mineral salts in milk[END_REF], indicating that no dramatic changes are induced by a 48°C temperature.

Microfiltration performance 4.2.1. Unmodified milk

Performance of the filtration of unmodified skimmed milk is presented in Fig 2 andTable 2. Before the CHC, in the range τw =120-80 Pa, the TMP is stable at 0.12 ± 0.01 ´ 10 5 Pa and the soluble protein transmission (Tr) is roughly constant at ~65 ± 4%. In the same time, the permeate turbidity increases, indicating a higher transmission of casein micelles through the membrane. When the CHC are reached (τwcrit = 55 Pa, Table 2), the TMP quickly increases and starts to oscillate (Fig 2). This increase in TMP lasts all along the CHC phase (Δtosc/Δtcrit = 1, Table 3), with an average slope of 4.25 Pa s -1 (Table 3 andFig [START_REF] Gésan-Guiziou | Critical stability conditions in crossflow microfiltration of skimmed milk: transition to irreversible deposition[END_REF]. We count ~56 of these oscillations, with an average frequency of 0.036 s -1 and an average amplitude of 1.6 ´ 10 -8 Pa (Table 3). Still in the CHC, the soluble protein transmission decreases from 65± 4 % down to 25 ± 3 % and the normalized turbidity also sharply decreases from 2.1 to 0.8, indicating a rapid decrease in casein micelles transmission.

After the CHC, the TMP oscillations stop immediately with the increase in τw. However, this increase does not fully restore the initial performance of the filtration. Indeed, TMPf is increased by 66 % (TMPf = 0.20 10 5 Pa, Table 2) compared to the average transmembrane pressure for τw = 120 Pa at the beginning of the parametric study (TMPi, Table 2). The transmission of the serum proteins is also affected since it decreases from Tri = 66±1% to Trf = 53±1% at the end of the parametric study (Table 2). This "hysteresis" in TMP and protein transmission indicates that the deposit formed at the membrane surface during the CHC partially remains after the CHC, and this until the end of the parametric study. At this point, the resistance of the casein micelles deposit Rd is about 2 times the resistance of the membrane (Table 2) and Rid corresponds to 60 % of the deposit resistance.

After the CHC, the permeate turbidity shows a quite surprising behavior with a strong increase at τw = 80-100 Pa and then a global decrease until reaching its initial value (Fig 2).

As it is already shown [START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF], the casein micelles are preferentially accumulated at the membrane surface. Hence, when the τw is increased again after the CHC, the sharp and quick decrease in TMP can lead to a "decompression" of the accumulated matter but the casein micelles accumulated do not returned instantaneously into the bulk. This phenomenon can result in a higher casein micelle concentration at the membrane surface and then in the increase in casein micelle transmission after the CHC for τw = 80-100 Pa.

Modified milks Impact of ionic strength on the overall filtration performance

The filtration performance obtained with modified milks is presented in Tables 2 and3 and in Fig 3 for the three conditions of I investigated. Regardless of the ionic strength, the TMP is always low (TMP ≈ 0.1-0.2 ´ 10 5 Pa) and stable before the CHC (120-80 Pa). Also the resistance of the polarization layer (Rp) is pretty much the same at the beginning of each run (Table 2) and this irrespective of I. These results confirm that it is reasonable to compare the deposit resistance between the experiments performed at different I.

It is noticeable that for the higher ionic strength (I = 0.372 M), the TMP shows a quite significant increase (up to 0.5 ´ 10 5 Pa) at τw = 80 Pa (Fig 3 c). This increase is still small compared to what is usually observed at the CHC so that we decided to further decrease the τw. But this small increase in TMP suggests that the CHC are more probably located in between τw = 80 Pa and τw = 72 Pa (Table 2, Fig 3 c).

With the increase in I, and for a constant convective mass transport (J constant) and a similar polarization layer before CHC (Similar TMP), (J/τw)crit decreases from 0.9 to 0.65 (Fig 5). As this critical ratio characterizes the formation of the irreversible deposit of casein micelles [START_REF] Gésan-Guiziou | Critical stability conditions in crossflow microfiltration of skimmed milk: transition to irreversible deposition[END_REF][START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF], this indicates that the deposit is formed at higher wall shear stress (τw) in high ionic strength conditions. Interestingly, the decrease in (J/τw)crit is linearly correlated to I and then to the electrophoretic mobility of the casein micelles and to the release in soluble calcium (Table 1).

During the CHC, regardless of I, a sharp and quick increase in TMP is observed, as already described with unmodified skimmed milk. However the TMP behavior at the critical step is strongly different according to the level of ionic strength. For I = 0.072 M (unmodified milk) the TMP oscillations are observed all along the critical step (Dtosc/ Dtcrit = 1), but for I = 0.157 and 0.247 M (Figs 3a andb) the TMP oscillations are observed only during the first part of the critical step (named TMP oscillation phase, Dtosc/ Dtcrit = 0.76 and 0.55 for I = 0.157 and 0.247, respectively); after the oscillation phase, the TMP reachs the value of 3.8 10 5 Pa, which corresponds to the maximum transmembrane pressure of the system (average retentate pressure set at 4.0 10 5 Pa and permeate compartment at atmospheric pressure). The TMP oscillations traduce a cyclic phenomenon of set-up/removal of the fouling layers accumulated at the membrane surface explained by a disengagement of the matter which returns into the bulk. A similar phenomenon has previously been observed with spikes of permeation flux during blood filtration [START_REF] Mcdonogh | Criteria for Fouling Layer Disengagement During Filtration of Feeds Containing A Wide-Range of Solutes[END_REF]. For the highest I (0.372 M), no oscillations are observed (Δtosc/Δtcrit = 0); TMP increases regularly during the first part of the critical step and then finally stabilizes at 3.8 10 5 Pa (Table 3

, Fig 3c).

With the increase in I, the number of TMP oscillations decreases from Nosc = 56 to 0, the oscillation frequency changes from Fr = 0.036 to 0 s -1 and the amplitude of oscillations, Amposc is multiplied by 10 from I= 0.157 to 0.247 M (Table 3). Moreover, the slope corresponding to the increase in TMP, SlTMP, calculated before TMP stabilization, changes from 4 to 22 Pa s -1 as I increases from 0.072 to 0.372 M (Table 3, Fig 4). Then the higher the ionic strength, the sharper the increase in TMP (SlTMP), the shorter the TMP oscillations phase, the lower the frequency of pressure fluctuations (Fr) and the larger the amplitude of pressure oscillations (Amposc) at the critical step.

After the CHC, the TMP decreases as τw increases. Interestingly, the higher the I, the more pronounced the pressure hysteresis (Fig 3), indicating a more cohesive deposit. Indeed, the TMP hysteresis (TMPf-TMPi, Table 2) increases from 0.08 to 1.70 10 5 Pa for I increases from 0.072 M to 0.372 M (Table 2). This TMP hysteresis is directly correlated to the increase in the total resistance of the deposit that remains at the membrane surface at the end of the parametric study (Rd = 0.98 10 12 m -1 for unmodified milk to 19.38 10 12 m -1 for milk at I= 0.372 M, Table 2). The absolute values of both the irreversible (Rid) and the reversible (Rrd) parts of this total resistance increase with ionic strength. However, the increase in Rrd is much more pronounced as Rid and Rrd is multiplied by 4.9 and 43.2 respectively for I = 0.072 to 0.372 M. As a result, Rrd represents 40 % of the total resistance for the unmodified milk, versus 84-90 % for I = 0.157-0.372 M (Table 2).

Impact of the increase in I on protein transmission

Before and during the CHC, the normalized turbidities and soluble protein transmissions show similar evolutions at all the ionic strengths investigated. Before the CHC, the casein micelles transmission regularly increases with a maximum at the beginning of the critical step and then sharply decreases during the critical step. Before the CHC, the soluble protein transmission is constant around 65-70 % (Figs 3) and then sharply decreases when the CHC are reached. The higher I, the lower Tr at the critical step: the average Tr is 38 ± 11 % for the unmodified milk , 33 ± 3 % for the milk at 0.247 M and 25 ± 5 % for the milk at 0.372 M.

After the CHC, the turbidity evolutions are different between the unmodified milk and the modified ones. Indeed for the unmodified milk the turbidity increases with the re-increase in τw, while for the modified milks, the turbidity regularly decreases down to 0.1 NTU (casein micelles concentration in permeate < 0.6 10 -3 g kg -1 ), indicating an irreversible alteration of casein micelle transmission, as previously observed by Le Berre and Daufin (1998) [START_REF] Berre | Microfiltration (0,1 µm) of milk: effect of protein size and charge[END_REF].

Beyond CHC, Tr increases except for the milk at the highest I, for what Tr still continues to decrease (final Tr = 17 %). Regardless of I, the re-increase in τw does not restore the filtration performance: the higher the ionic strength, the lower the final soluble protein transmission (Trf), (Table 2).

Discussion

As shown previously, the increase in ionic strength through NaCl addition drastically modifies the performance of milk microfiltration when the casein micelles are sufficiently concentrated at the membrane surface, i.e. near and at the CHC: lower critical ratio (J/τw)crit, sharper increase in TMP and higher alteration of protein transmission at high ionic strengths.

In this discussion we show that those changes in behavior can be related to differences in the properties of the casein micelles deposit because of the different ionic strength conditions. We also explain our results based on previous observations done on the behavior of casein micelles in concentrated regime.

Deposit formation at lower casein micelle concentrations

A first important result in our experiments is that the CHC are located at a lower critical ratio (J/τw)crit when ionic strength is increased . This indicates that the casein micelle deposit is formed at lower casein concentrations at the membrane when adding salt. Indeed and as it has been already shown by several authors [START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF][START_REF] Berre | Skimmilk crossflow microfiltration performance versus permeation flux to wall shear stress ratio[END_REF], the casein micelles accumulated turn into a gel when their concentration exceeds a critical value at the membrane surface. Moreover, it seems reasonable to assume that the polarization phenomenon is not affected by a change in ionic strength: first the TMP values are pretty similar before the CHC, between the experiments performed at different I and in addition the casein micelle diameter being the same for each milk studied, hence its diffusivity is constant irrespective of the milk ionic strength. Hence, the fact that the CHC are located at higher τw (i.e., lower (J/τw)crit) for higher ionic strengths indicate that the deposit is formed in weaker conditions of polarization and consequently at a lower "critical" concentration. This result is fully consistent with the observations of Bouchoux et al. [START_REF] Bouchoux | Casein Micelle Dispersions under Osmotic Stress[END_REF] who prepared casein micelle gels through osmotic compression and who found that an increase in ionic strength leads to a shift in the gel transition to lower osmotic pressures and volume fractions.

The formation of a gel layer at a lower concentration in casein micelle can be explained by a modification of the shift in the balance of interactions between the micelles. With the increase in ionic strength, the liquid-gel transition is linearly correlated to the decrease in the global charge of the casein micelles (represented by electrophoretic mobility), and the increase in soluble calcium. The gel transition takes place when the distance between the micelles is small enough to induce changes in interactions: the attractive forces between the micelles become stronger than repulsive ones, thus coagulation of the accumulated micelles takes place and the deposit is formed [START_REF] Jimenez-Lopez | Role of milk constituents on critical conditions and deposit structure in skimmilk microfiltration (0.1µm)[END_REF][START_REF] Krstics | Static turbulence promoter in cross-flow microfiltration of skim milk[END_REF]. With the increase in ionic strength, there is a decrease in the overall repulsions between casein micelles as:

-Steric repulsions become weaker because of the collapse of the κ-caseins brush, located at the surface of the casein micelles [START_REF] Tuinier | Stability of casein micelles in milk[END_REF] -Repulsive electrostatic interactions are affected through the screening of the surface charge of the casein micelle (represented by the electrophoretic mobility, Table 1). Note that some authors consider that, electrostatic interactions must be neglected at the natural ionic strength of milk [START_REF] Kühnl | Impact of colloidal interactions on the flux in cross-flow microfiltration of milk at different pH values: A surface energy approach[END_REF]. Nevertheless, Tessier et al. [START_REF] Tessier | Correlation of diafiltration sieving behavior of lysozyme-BSA mixtures with osmotic second virial crosscoefficients[END_REF] clearly show that electrostatic interactions still exist between proteins at 0.3 M of ionic strength.

As the overall repulsions decrease, the attractive interactions can take place and make easier the liquid-gel transition. Indeed, the hydrophobic forces are probably enforced at high ionic strength [START_REF] Bonnet-Gonnet | Osmotic pressure of latex dispersion[END_REF]. Moreover, as the soluble calcium concentration increases with I, the establishment of calcium bonds between the micelles [START_REF] Müller-Buschbaum | Effect of calcium concentration on the structure of casein micelles in thin films[END_REF] through the carboxylic acids or through sialic acids worn by the κ-caseins [START_REF] Mellema | Molecular mechanism of the renneting process of casein micelles in skim milk, examined by viscosity and lightscattering experiments and simulated by model SCF calculations[END_REF] can also take place.

A more cohesive deposit

Once the deposit is formed at the CHC, its characteristics clearly depend on the ionic strength. More particularly, the deposit appears to be more cohesive at higher ionic strengths. This increase in cohesiveness is visible through the TMP evolution (average value and oscillations) and the deposit irreversibility (Rid). At the CHC, the oscillations in TMP can be seen as the result of a cyclic phenomenon of matter accumulation and removal at and from the membrane surface. With the increase in I, these TMP oscillations are larger in amplitude but decrease in number, frequency and duration; thus indicating that the deposit is "harder" to remove from the membrane and so more cohesive. Moreover, this increase in the cohesiveness is also shown through the delay observed (after the CHC) between the increase in τw and the decrease in TMP which is much more pronounced for the higher I. Finally, this assumption is strengthened through the more pronounced TMP hysteresis and the less reversible deposit, as seen through the value of Rid (Table 2). In accordance with these results, Bouchoux et al. [START_REF] Bouchoux | Casein Micelle Dispersions under Osmotic Stress[END_REF] point out the fact that gels of casein micelles made through osmotic compression at high ionic strengths are more difficult to re-disperse, and consequently more cohesive. Once again, this phenomenon can be explained through a shift in the balance of attraction and repulsion forces between the micelles. Indeed, the hydrophobic groups that may be hidden by the κ-casein brush at milk ionic strength may become more accessible when this brush is partially collapsed at higher ionic strengths. So the micelles can stick to each other more efficiently as they come into direct contact. In addition, the casein network that constitutes the deposit may be strengthened by calcium bridges between the micelles and/or between the micelles and the membrane which can participate to the increase in the cohesiveness and the irreversibility.

A Deposit with modified selectivity

When the deposit is formed, the increase in ionic strength induces an irreversible alteration of the proteins transmission: the casein micelle transfer is immediately and irreversibly affected as the casein micelle transmission is close to zero (permeate turbidity < 0.1 NTU) at high ionic strengths while the transmission of soluble proteins decreases more gradually from 53 % to 17 % at I = 0.072 to 0.372 M, respectively. This change in the selectivity of the casein deposit can be attributed to several phenomena. First, a deposit more cohesive is probably less porous, i. e., made of a "tighter" casein network. This would result in some stronger size exclusion effects that partly explain the alteration of the protein transfer, which is much more pronounced for the larger proteins (average casein micelle diameter ≈ 100 nm VS soluble protein average diameter < 10 nm).

Second, the deposit is made of a quantity of matter accumulated that is more important at higher ionic strength (increase in Rrd with I, Table 2). This high density of negative charges possibly generates strong repulsive electrostatic interactions between the casein deposit and the proteins in solution in the bulk. Some authors already showed that significant electrostatic interactions still exist between proteins even for I=0.372 M [START_REF] Tessier | Correlation of diafiltration sieving behavior of lysozyme-BSA mixtures with osmotic second virial crosscoefficients[END_REF]. Hence, the casein micelle deposit could constitute a material with a high negative charge density at the membrane surface which could repel other negatively-charged proteins (casein micelles, major soluble proteins such as α-lactalbumin and β-lactaglobulin) and consequently alter their transfer. To investigate these two hypotheses, size exclusion phenomena and/or electrostatic repulsions, further experiments aiming at the transmission of individual soluble proteins as a function of their diameter and charge properties will be carried out.

Conclusions

The microfiltration 0.1 µm of skimmed milk at various ionic strengths allows a better understanding in the structuration of the casein micelle deposit which rules the performance (permeation flux, transmembrane pressure and protein transmission) of the operation. The experiments performed in this study, allow to put in light the role of the interactions between casein micelles in the set-up and in the properties of the deposit. With the increase in ionic strength, the global charge of the casein micelles is reduced and the concentration of soluble calcium is increased. These physico-chemical modifications lead to a decrease in repulsive interactions between micelles by the screening of the casein micelle charge and by the collapse of the κ-casein. The hydrophobic interactions between micelles are then favored and are possibly strengthened by calcium bonds between micelles. As a consequence, the micelles can come into direct contact more easily and can form stronger intermolecular links. Hence, with the increase in ionic strength, the deposit is set-up for a lower casein micelle concentration at the membrane surface and it is harder to remove by the increase in its cohesiveness. Finally, the higher the ionic strength, the more pronounced is the decrease in protein transmission. This can be related to both an increase in steric repulsion mechanism and an increase in electrostatic repulsion. The deposit appears less porous possibly because of the reduction of the interstitial spaces inside the casein network. Additionally, some repulsive electrostatic interactions between the deposit and the proteins can still exist even at high ionic strength; hence the negative charged proteins in suspension in the bulk are possibly rejected when they approach the deposit. All these results are very consistent with previous studies of the casein micelle behavior in concentrated regime (osmotic compression).

  regression TMP = f (t) (Fig 4), Aplot corresponds to the filled area between the two plots (Fig 1) as:

12 m - 1 .

 1 TMPi = average TMP for τw = 120 Pa before CHC TMPf = average TMP for τw = 120 Pa after CHC Tri = average Tr for τw = 120 Pa before CHC Trf = average Tr for τw = 120 Pa after CHC I (M) t wcrit (Pa) TMP i 10 5 (Pa) TMP f 10 5 (Pa) R f 10 12 (m -1 ) R p 10 12 (m -1 ) R d 10 12 (m -1 ) R id 10 12 (m -1 ) R rd 10

Figure 1 :Fig 2 :Figure 3 :

 123 Figure 1: Main parameters for the characterization of TMP oscillations at the critical hydrodynamics conditions (CHC). SlTMP (Pa s -1 ): slope of the TMP signal versus time after linearization corresponding to the overall increase in TMP during the oscillation phase Dtosc (-): normalized duration of the TMP oscillation phase Dtcrit (-): normalized total duration of the critical step Nosc (-): number of oscillations Fr (s -1 ): TMP oscillations frequency as Fr =

Figure 4 :

 4 Figure 4: Evolution of the transmembrane pressure (TMP) at the critical step (τwcrit) versus time (see calculation 3.3) during microfiltration of skimmed milk at different ionic strengths (I = 0.072-0.372 M).Operating conditions: T = 48°C, J = 50 L h -1 m -2 , pH = 6.55 Time scale corresponded to the normalized time with reference (t= 1) to the beginning of the oscillation phase.

Figure 5 :R 2

 52 Figure 5: Impact of the ionic strength on the critical ratio (J/τw)crit during microfiltration of skimmed milk.conditions: T = 48°C, J 50 L h -1 m -2 , pH = 6.55

Table 1 :

 1 Physico-chemical characteristics of the casein micelles and concentration of total (Catotal) and soluble calcium (Cas) at 48 °C in thermized skimmed milk presenting various ionic strengths (I) I: Ionic strength *Catotal concentration slightly changed because of differences in milks (season, cow species, food of cows…) during the duration of the experimental period.

	Thermized skimmed milks	Casein micelles characteristics		Calcium content
		Average diameter nm	Electrophoretic mobility µ (10 -7 m 2 s -1 V -1 )	*Ca total g L -1	Ca s g L -1	Ca s /Ca total
	reference I = 0.072 M	174±13	31.1±1.2	1.19±0.01	0.23±0.01 0.19±0.02
	I = 0.157 M	187±14	27.3±1.2	1.18±0.01	0.24±0.01 0.20±0.02
	I = 0.247 M	181±1	22.9±1.2	1.18±0.02	0.25±0.01 0.21±0.02
	I = 0.372 M	184±11	15.5±1.2	1.21±0.01	0.28±0.01 0.23±0.02

Table 2 :

 2 Critical wall shear stress (τwcrit), average of the transmembrane pressure (TMP), hydraulic resistances (total fouling (Rf), polarized layer (Rp), deposit (Rd), irreversible (Rid) and reversible (Rrd) part of the deposit) and total soluble protein transmission (Tr) during microfiltration of skimmed milk at various ionic strengths. The initial membrane resistance is Rm = 0.51 ± 0.02 ´ 10

Acknowledgements

The authors wish to thank L. Fromont and C. Geneste for technical support and precious advices.

Nomenclature

Aplot :

Filled area under two plots (Pa s) Amposc :

Amplitude of oscillation at CHC (Pa) CHC :

Critical Hydrodynamic Conditions Catotal :

Total calcium (g L -1 ) Cas :

Soluble calcium (g L -1 ) CP :

Soluble protein concentration in the permeate (g L -1 ) CR :

Soluble protein concentration in the retentate (g L -1 ) Fr :

Oscillations frequency (s -1 ) I :

Ionic strength (M) J :

Permeation flux (L h -1 m -2 ) NCN :

Non casein nitrogen NPN:

Non protein nitrogen N :

Number of maxima and minima of TMP during the oscillation period Nosc :

Number of TMP oscillations PAR :

Average retentate pressure (Pa) Rd :

Deposit resistance (m -1 ) Rid :

Irreversible deposit resistance (m -1 ) Rrp :

Reversible deposit resistance (m -1 ) Rf :

Total fouling resistance (m -1 ) Rm :

Membrane resistance (m -1 ) Rp :

Polarization layer resistance (m -1 ) SlTMP :

Slope of the linear regression of the increase in the TMP oscillation phase (Pa s -1 ) TMP :

Transmembrane pressure (Pa) TMPi :

Transmembrane pressure at the beginning of the parametric study for τw = 120 Pa TMPf :

Transmembrane pressure at the end of the parametric study for τw = 120 Pa Tr :

Soluble protein transmission (%) Tri :

Soluble protein transmission (%) at the beginning of the parametric study for τw = 120 Pa Trf :

Soluble protein transmission (%) at the end of the parametric study for τw = 120 Pa VRR :

Volume reduction ratio Greek symbols: τw : Wall shear stress (Pa) τwcrit :

Critical wall shear stress (Pa) Dtosc :

Oscillation phase duration (s) Dtcrit :

Critical step duration (s)


> La fusion membranaire est un mécanisme ubiquitaire, qui est à la base de nombreux phénomènes physiologiques fondamentaux tels que la fécondation d'un ovule par un spermatozoïde, le développement musculaire, l'infection virale, ou encore l'ensemble du trafic vésiculaire au sein des cellules eucaryotes [1,2]. Tous ces processus impliquent que deux compartiments cellulaires, initialement séparés et chacun délimité par une membrane biologique, entrent en contact puis mettent en commun leurs composants lipidiques et cytosoliques. Les membranes biologiques étant des structures naturellement très stables, les événements de fusion ne peuvent se produire spontanément et nécessitent l'intervention de protéines spécialisées, qui vont aider à surmonter les barrières énergétiques de chaque étape, conduisant à la fusion membranaire [3]. Ces étapes comprennent : (1) le rapprochement des membranes ; (2) la rupture de leur structure en bicouche (avec passage éventuel par un intermédiaire d'hémifusion où les feuillets externes des bicouches ont fusionné tandis que les feuillets internes demeurent séparés), et (3) la création d'un pore de fusion permettant le passage des composés solubles.

Les SNARE orchestrent la fusion membranaire lors du trafic vésiculaire

Comment se forment les complexes SNAREpins ?

Pour répondre à cette question, il est nécessaire de pouvoir mesurer les interactions, à l'échelle nanométrique, entre SNARE incorporées dans des membranes lipidiques. Nous avons utilisé l'appareil de forces entre surfaces (SFA), qui permet de mesurer simultanément la distance et l'énergie d'interaction entre deux bicouches lipidiques fonctionnalisées [7]. Les surfaces d'étude du SFA sont deux feuillets de mica, placés dans une configuration de type cylindres croisés, dont une face est recouverte d'argent et l'autre de la bicouche lipidique. Un système de réflexions multiples de la lumière blanche sur les couches d'argent permet d'obtenir optiquement la distance, tandis que la déflexion d'un ressort plan relié à l'un des deux micas donne accès à la force entre les deux surfaces (Figure 1) ; cette force est proportionnelle à l'énergie d'interaction entre deux surfaces planes équivalentes. Au cours d'une expérience SFA, les surfaces sont tout d'abord rapprochées l'une de l'autre jusqu'au contact afin d'identifier la nature des interactions agissant à différentes distances, puis elles sont séparées afin de mesurer leur énergie d'adhésion. Les profils d'interaction SFA entre v-SNARE et t-SNARE insérées dans des bicouches lipidiques [8] ont révélé que les SNAREpins commencent à s'assembler lorsque les membranes sont à environ 8 nm l'une de l'autre (Figure 2). In vivo, les protéines d'arrimage (Rab et facteurs d'attachement) devront donc positionner la vésicule à cette distance de la membrane plasmique cible afin que les SNARE puissent agir. L'énergie emmagasinée dans les SNAREpins -énergie qui pourra être utilisée pour la fusion -augmente au fur et à mesure que les membranes sont rapprochées l'une de l'autre et atteint son maximum lorsqu'elles se trouvent à 4 nm de distance. Dans cette configuration, les SNAREpins sont partiellement déstructurés au niveau de leur partie carboxy-terminale, et constituent un complexe très stable, irréversible à l'échelle de temps des phénomènes biologiques.

Ces complexes permettront donc de maintenir un contact très étroit entre les membranes durant l'organisation moléculaire de la fusion membranaire, et seront très probablement les substrats de certaines protéines de régulation, telles que Synaptotagmine et Complexine [9]. Nos mesures SFA montrent que l'énergie de ces SNAREpins partiellement assemblés (35 k B T) est très proche de l'énergie nécessaire à la formation des intermédiaires d'hémifusion (40-50 k B T), et que leur partie carboxy-terminale interagit physiquement avec la structure membranaire. Cette interaction permettra de transmettre l'énergie libérée durant l'assemblage des SNAREpins à la bicouche lipidique, et ainsi de favoriser la formation des structures intermédiaires de type non bicouche conduisant à la fusion.

Figure 1 .

 1 Figure 1. Schéma de principe de l'appareil de forces entre surfaces (SFA).Le SFA mesure la force entre deux surfaces de mica recouvertes des bicouches à étudier en fonction de la distance qui les sépare. La distance est obtenue par un système d'interférométrie de type Fabry-Pérot, avec une résolution de 0,1 nm. La force est déduite de la déformation d'un ressort à flexion connecté à l'une des deux surfaces, avec une précision de 1 μN.

Figure 2 .

 2 Figure 2. Dynamique et énergétique d'assemblage desSNAREpins. D'après nos mesures SFA, les protéines SNARE, initialement faiblement structurées, commencent à s'assembler sous la forme de complexes inter-membranaires SNAREpins lorsque les bicouches sont distantes de 8 nm. L'énergie de ces SNAREpins augmente lorsque la distance inter-membranaire diminue, et culmine à 35 k B T, une valeur très proche de l'énergie permettant de faire hémifusionner les membranes.

Distance Force Lumière blanche Déflexion du ressort Franges d'interférence Bicouches déposées sur du mica argenté

  La plupart des études biochimiques et biophysiques réalisées à ce jour ont porté sur les complexes « post-fusion » cis-SNARE [6]. Les complexes SNAREpins, en raison de leur caractère intrinsèquement dynamique, n'ont jamais pu être ni isolés, ni caractérisés. Quelle énergie est déployée durant l'assemblage des SNAREpins et comment cette énergie estelle transmise aux membranes afin de les rapprocher et les remodeler restent des questions actuellement sans réponse.

Article disponible sur le site http://www.medecinesciences.org ou http://dx.doi.org/10.1051/medsci/2008242142

M/S n° 2, vol.24, février 2008


Damien F G Minenna
email: damien.minenna@univ-amu.fr

Yves Elskens

Frédéric André

Alexandre Poyé

Jérôme Puech

Fabrice Doveil

DIMOHA: a time-domain algorithm for traveling-wave tube simulations

Keywords: numbers: 45.20.Jj (Lagrangian and hamiltonian mechanics), 52.40.Mj (Particle beam interaction in plasmas), 84.40.Fe (Microwave tubes) Traveling-wave tubes (TWTs), DIMOHA, simulation, AM/AM, AM/PM, nonlinear signals, time domain, wave-particle interaction, helix, folded waveguide, gratings, slow-wave structure, Hamiltonian, N-body dynamics, communication systems, vacuum electronics

To simulate traveling-wave tubes (TWTs) in time domain and more generally the wave-particle interaction in vacuum devices, we developed the DIscrete MOdel with HAmiltonian approach (dimoha) as an alternative to current particle-in-cell (PIC) and frequency approaches. Indeed, it is based on a longitudinal N -body Hamiltonian approach satisfying Maxwell's equations. Advantages of dimoha comprise: (i) it allows arbitrary waveform (not just field envelope), including continuous waveform (CW), multiple carriers or digital modulations (shift keying); (ii) the algorithm is much faster than PIC codes thanks to a field discretization allowing a drastic degree-of-freedom reduction, along with a robust symplectic integrator; (iii) it supports any periodic slow-wave structure design such as helix or folded waveguides; (iv) it reproduces harmonic generation, reflection, oscillation and distortion phenomena; (v) it handles nonlinear dynamics, including intermodulations, trapping and chaos. dimoha accuracy is assessed by comparing it against measurements from a commercial Ku-band tapered helix TWT and against simulations from a sub-THz folded waveguide TWT with a staggered double-grating slow-wave structure. The algorithm is also tested for multiple-carriers simulations with success.

I. INTRODUCTION

Nowadays, simulating traveling-wave tubes [START_REF] Minenna | The Traveling-Wave Tube in the History of Telecommunication[END_REF] (TWTs) is still a challenge due to the large number of parameters involved and because, in nonlinear regime, the large power of high frequency TWTs and their broad band spectrum generate critical instabilities. The physical process at the heart of the device is the wave-particle interaction, whereby electrons transfer their momentum to amplify telecommunication signals, as, e.g., an inverse Landau damping: synchronization [START_REF] Doveil | Experimental observation of nonlinear synchronization due to a single wave[END_REF] occurs when the phase velocity of the wave (determined by the tube geometry) is close to the particles speed.

Currently, to model the propagating wave interacting with the beam, two kinds of numerical approaches are used: time domain and frequency (steady-state) domain modeling. The first kind allows, in principle, to investigate the majority of TWT problems such as reflections, oscillations, harmonic and intermodulation generations, multi-frequency signals. But its main drawbacks are the computation size and duration. This is especially true when studying oscillations and intermodulations that require strong resolution in time to be accurate. Generally, time domain codes are based on three-dimensional (3D) meshes (thousands to millions of points) based on Maxwell's equations such as particle-in-cell (PIC) codes relying on kinetic (e.g. Vlasov) equations (like cst [START_REF] Safi | A Traveling-Wave Tube Simulation Approach With CST Particle Studio[END_REF][START_REF]CST Particle Studio[END_REF], karat [START_REF] Tarakanov | User's Manual for Code KARAT[END_REF] or mafia [START_REF] Weiland | MAFIA Version 4[END_REF]). To reduce the computational costs, one can use specialized frequency (harmonic or steady-state) models (like mvtrad [START_REF] Waller | Modélisation numérique de l'interaction et diagnostic expérimental du faisceau d'électrons dans un tube à ondes progressives spatial[END_REF], christine [START_REF] Antonsen | CHRISTINE: A multifrequency parametric simulation code for traveling wave tube amplifiers[END_REF] or bwis [START_REF] Li | Theory and Design of Microwave-Tube Simulator Suite[END_REF]) which rely on reduction models such as equivalent circuits envelope models.

We propose a third option, using an N -body (a.k.a. many body) description [START_REF] Minenna | Electronwave momentum exchange and time domain simulations applied to traveling wave tube[END_REF][START_REF] Minenna | Proc. IEEE Int. Vacuum Electron. Conf. (IVEC)[END_REF][START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF] to design a specialized theory in time domain and large-signal regime as an alternative to both PIC and frequency approaches. By nature, this description is deemed slower and computationally far more expensive than PIC since all particles should be considered. For periodic structures, we overcome this slowness thanks to a drastic degrees-of-freedom (dof) reduction [10-13] (Kuznetsov's discrete model). The reduction is extreme (∼ 100 000 dofs for space TWTs) compared to finite difference techniques in PIC codes that often involve several millions of dofs to reach the same accuracy. Moreover, our approach is one-dimensional (1D), while generally time domain codes require simulating the whole slow-wave structure (SWS) of the tube in 3D. Therefore, the model consistency and conservation properties are more easily achieved. This model reduction is combined with the N -body hamiltonian approach [START_REF] Minenna | Electronwave momentum exchange and time domain simulations applied to traveling wave tube[END_REF][START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF] leading to a better control on conservation properties, including Poincaré-Cartan invariants (dynamics geometry)

Pre-processing

Cold-field parameters

Beam management

Discrete model

Initial positions and velocities

Precomputed space charge

Electrodynamics and Rowe's model

Beam inputs

Signal inputs

Time integrator

Wave propagation

Beam propagation I. [START_REF] Arnold | Mathematical Methods of Classical Mechanics[END_REF]. From this, we built a symplectic [START_REF] Hairer | Geometric numerical integration[END_REF][START_REF] Saitoh | The Symplectic Finite Difference Time Domain Method[END_REF] algorithm dimoha (DIscrete MOdel with HAmiltonian approach), enabling us to increase the numerical time step without incurring too much error on results. To assess its validity, we compare this algorithm with measurements from a commercial 140W Ku-band tapered helix TWT and with simulations from a sub-THz staggered TWT. This paper is dedicated to presenting our algorithm, its capacities and some results obtained with it. However, to keep the paper short and focus on results, we defer the overview of the underlying theory to a forthcoming report.

N-body model

Wave coupling terms

Shape fields

Boundaries

Symplectic integrator

Tube splitting parallelisation

II. THE ALGORITHM

A. Algorithm description

The current version of dimoha is developed jointly between Aix-Marseille Université-CNRS, Thales AVS/MIS and the Centre National d' Études Spatiales in France. It is written in modern Fortran and parallelised with the Message Passing Intergace (MPI) librairies. dimoha simulates, in 1D, the interaction between electrons and waves in the slow-wave structure (SWS) of industrial TWTs.

In time domain, a comparison with mvtrad [7] (a 2.5D frequency-domain (field) particle-in-cell (beam) code specialized for helix TWT) has already shown [START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF] perfect agreement at power saturation, in nonlinear regime (trapping), in the middle of the transmission band of a basic tube without taper nor attenuators, using a preliminary version of our algorithm. The novelty, here, is that we adapted the model to take into account non-periodic structure (taper) and attenuators that are major compo- nents of industrial TWTs. The main vocation of dimoha is to provide a complete tool that can be used for industrial design and research activities as well as for device assessement by telecom operators.

As shown on Fig. 1, dimoha is composed of three parts (see Sections III-IV). The first part (pre-processing) computes the cold-field parameters from the tube geometry. The second part (time integrator) is the heart of the algorithm that simulates the wave-electron interaction in the slow-wave structure up to a given time t max . The last part (post-processing) performs final diagnostics after the interaction. The parametric inputs are given by Table I. In continuous wave regime, only the frequency and power of the injected wave are required.

B. The discrete model

Our algorithm dimoha is based on the Kuznetsov discrete model [START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF][START_REF] André | Hamiltonian description of self-consistent wave-particle dynamics in a periodic structure[END_REF][START_REF] Bernardi | Efficient Time-Domain Simulations of a Helix Traveling-Wave Tube[END_REF][START_REF] Bernardi | Control of the Reflections at the Terminations of a Slow Wave Structure in the Nonstationary Discrete Theory of Excitation of a Periodic Waveguide[END_REF][START_REF] Ryskin | Nonstationary nonlinear discrete model of a coupledcavity traveling-wave-tube amplifier[END_REF][START_REF] Terentyuk | Discrete model of a folded-waveguide traveling-wave tube[END_REF][START_REF] Minenna | Recent discrete model for smallsignal analysis of traveling-wave tubes[END_REF][START_REF] Kuznetsov | On one form of excitation equations of a periodic waveguide[END_REF][START_REF] Ryskin | Nonstationary nonlinear modeling of an electron beam interaction with a coupled cavity structure. I. Theory[END_REF] which uses an exact field decomposition for periodic structures for large-signal regimes. The SWS is directed along the z-axis with pitch d and cell index n ∈ Z. We assume the fundamental mode of propagation is dominant, so we neglect other modes. The scalar potential φ sc of the beam space charge satisfies Poisson's equation. Within the discrete model, the RF electric (circuit part) and magnetic fields in the whole structure read

E(r, t) = n∈Z V n (t)E -n (r) , (1)
H(r, t) = n∈Z iI n (t)H -n (r) , (2)
that satisfy Maxwell's equations. V n , I n are the celldiscretized temporal field amplitudes, and E n , H n are the "shape" functions of the fields depending only on the structure geometry (see Eq. ( 5)). The imaginary number (i) and the sign in -n are technical conventions [START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF][START_REF] André | Hamiltonian description of self-consistent wave-particle dynamics in a periodic structure[END_REF] allowing terms V n , I n , E n , iH n to be real. To compute time variables V n , I n , we need an interaction model like the N -body description [START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF] (used below) or a fluid model [START_REF] Minenna | Recent discrete model for smallsignal analysis of traveling-wave tubes[END_REF]. Since ∇ × A = µ 0 H, we also express the vector potential A = n iI n A -n with the same I n as for the magnetic field [START_REF] Doveil | Experimental observation of nonlinear synchronization due to a single wave[END_REF].

Equations ( 1)-( 2) highlight the interest of the discrete model for dof reduction. Indeed, the time-independent functions E n , H n are given by the SWS geometry. For a space TWT comprising n max = 200 cells, there are only 200 variables V n , I n , so that the field involves only 200 dofs.

Beside the time domain approach, one can use a frequency domain version of Kuznetsov's discrete model [START_REF] Terentyuk | Discrete model of a folded-waveguide traveling-wave tube[END_REF][START_REF] Minenna | Recent discrete model for smallsignal analysis of traveling-wave tubes[END_REF][START_REF] Ryskin | Nonstationary nonlinear modeling of an electron beam interaction with a coupled cavity structure. I. Theory[END_REF][START_REF] Ryskin | [END_REF] in small signal regime. Its validity was assessed [START_REF] Minenna | Recent discrete model for smallsignal analysis of traveling-wave tubes[END_REF] by comparing it with Pierce's four-waves theory [START_REF] Pierce | Traveling Wave Tubes[END_REF].

III. INITIALISATION OF COLD-FIELD PARAMETERS

The first part of dimoha computes cold-field (timeindependent) parameters from tube and beam inputs in Table I. Those cold-field parameters, depending on the SWS, can be obtained with a numerical solver or measured.

A. Coupling terms and attenuators

Using the decomposition (1)-( 2), we rewrite sourceless Maxwell equations (equivalently Hamilton equations in our formulation [START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF]) as

dV n dt = - n Ω n-n I n (t) , (3)
dI n dt = n Ω n-n V n (t) . (4)
Those equations represent the wave propagation along the SWS behaving as a harmonic oscillators chain. Coupling terms Ω n-n between cells n and n are obtained from the Fourier series of the dispersion relation ω(β) = N ph n=-N ph Ω n e inβd , with β the wavenumber (propagation constant). We take the range of coupling N ph , estimated in [START_REF] Bernardi | Efficient Time-Domain Simulations of a Helix Traveling-Wave Tube[END_REF], at 5 for helix TWTs and 1 for folded waveguides. To take into account any pitch tapering, coupling terms are recomputed at each different cell. Note that, to compute accurately Ω n , the discrete model requires a broader dispersion relation (e.g. all harmonics) than just the data for the transmission band. This is an advantage since the wave propagation takes into account the whole dispersion relation and not a mere single point at the frequency of interest as for envelope models.

To consider attenuators [START_REF] Bernardi | Efficient Time-Domain Simulations of a Helix Traveling-Wave Tube[END_REF][START_REF] Bernardi | Control of the Reflections at the Terminations of a Slow Wave Structure in the Nonstationary Discrete Theory of Excitation of a Periodic Waveguide[END_REF], we add -α n V n to the right side of Eq. ( 3). dimoha can process spatially distributed attenuations in tubes. We can either provide in input the coefficients in Np/m or dB/wavelength from the tube permittivity or provide the total transmission loss S 12 from measurements and estimate the average coefficients. For tubes with severs and local attenuators, the α n 's are adjusted consequently.

B. Tube's meshes

The simulation uses two meshes. The position of each cell, spaced by pitch d, is given by z cell . To reproduce accurately the space charge field and the coupling between the field and particles, an adapted mesh z field oversamples the basic mesh by N os points per cell. For TWTs, N os does not need to be large (e.g. N os = 20 in our simulations).

Without beam, Eqs ( 3)-( 4) can be seen as the equations of a harmonic oscillators chain, so an abrupt stop of the SWS at the edge of the tube would lead to a total reflection of the signal. To avoid this, we lengthen the tube and use another coefficient α n that slowly increases when moving away from the physical part (perfectly matched layer (PML) method).

C. Eigenfields

Before computing time-dependent amplitudes, we also need the "shape" basis fields E n , H n , A n . In one dimension, the system is projected along the longitudinal z-axis and we can assume H z,n (z) = 0. For each cell, the on-axis shape fields E z,n for the electric field and A z,n for the vector potential are obtained as E z,n (z) = (2π) -1 π -π E z,β (z) e -inβd d(βd), from their discrete Fourier transforms (or inverse Gel'fand βtransform) E z,β and A z,β = -iE z,β /ω [START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF]. From the electromagnetic power, those eigenfields are computed as [START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF][START_REF] Ryskin | Nonstationary nonlinear modeling of an electron beam interaction with a coupled cavity structure. I. Theory[END_REF]

E z,β (z) = e -iβz β 1 d ω(β)Z c (β)v g (β) , (5)
where ω the angular frequency, v g = ∂ω ∂β the group velocity and Z c the impedance. To take into account any pitch tapering, eigenfields are recomputed at each different cell. For TWTs, shape fields E z,n , A z,n are similar to cardinal sine functions centered at the n th cell position [START_REF] Bernardi | Efficient Time-Domain Simulations of a Helix Traveling-Wave Tube[END_REF]. In particular, E z,n (z), A z,n (z) are computed for each points of the mesh z field along the whole structure. Those functions are needed to express the coupling with particles in the time integrator.

D. Beam and space charge

Since RF shape functions A z,n (responsible for the coupling with particles) are continuous and smooth, we can ease the computation by aggregating particles as macroelectrons. In our simulations, each macro-electron (indexed by k) is characterized by its position q k (t) and velocity v k and contributes a source -Ne k iev k A z,n (q k ) to [START_REF] Safi | A Traveling-Wave Tube Simulation Approach With CST Particle Studio[END_REF]. Their electric charge is e = -I 0 δ/v 0 , with δ their initial spacing and I 0 the cathode current. Even if the relativistic correction is small, we build the model and algorithm with relativistic variables. The initial beam velocity is v 0 = c(1 -(1 -e e V 0 /(m e c 2 )) -2 ) 1/2 , with V 0 the cathode potential, e e the electron charge, m e the electron mass, assuming the initial beam to be a line of equally spaced particles. The mass of macro-electrons is m = m e e/e e . We can also modify the velocity of each macro-electron to generate pulsed emissions or distribution functions like a gaussian beam.

As the discrete model decomposes only the circuit part of fields, it misses the space charge field. For N e macroelectrons, the axial space-charge field in one dimension E sc (q k ) = -Ne k =k ∂ z φ sc (q k -q k ) follows Rowe's model [START_REF] Bernardi | Efficient Time-Domain Simulations of a Helix Traveling-Wave Tube[END_REF][START_REF] Rowe | Nonlinear Electron-Wave Interaction Phenomena[END_REF]

E sc (q k ) = Ne k =k K exp - Q k -Q k b/2 sign Q k -Q k b/2 , (6)
with K = 2e/(π 0 b 2 ) < 0, for a cylindrical beam with radius b (effects of the helix radius are neglected [START_REF] Rowe | Nonlinear Electron-Wave Interaction Phenomena[END_REF]), where Q k is the mesh point nearest to q k (t) on the mesh z field . The space charge field E sc is pre-computed before the interaction (the charges and masses of particles are constants). After each displacement of macroelectrons, the space charge force is applied on them with this nearest-meshpoint approximation, allowing to drastically reduce the computational time.

IV. TIME INTEGRATOR AND DIAGNOSTICS

The second and main part of dimoha is the numerical time integrator. The evolution of our system rests on a total hamiltonian [START_REF] Jackson | Classical electrodynamics[END_REF] initially composed of the electromagnetic term V ( 0 |E| 2 + µ 0 |H| 2 )/2 dV , the particles term k γ k mc 2 (for macro-electrons), with γ k the Lorentz factor, and the space charge term Ne k =k φ sc (q k -q k ) /2. This hamiltonian is re-expressed [START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF] with the discrete model and projected in 1D along the longitudinal axis. Associated Hamilton's equations

dV n dt = - n Ω n-n I n (t) - Ne k=1 ie qk A z,n (q k ) , (7)
dI n dt = n Ω n-n V n (t) , (8)
v k = 1 γ k m (p k -eA z (q k )) , (9) ṗk
= 1 γ k m (p k -eA z ) ∂A z ∂z (q k ) + eE sc (q k ) , (10)
are the evolution equations of our sytem and are solved with our algorithm. The dynamical variables become the amplitudes V n (t), I n (t) of fields and the 1D positions q k (t) and velocities v k (t) of macro-electrons. This dynamics is consistent with Maxwell equations and Lorentz force.

Our numerical integrator solves the evolution system from our N -body self-consistent hamiltonian model [START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF] and provides V n (t), I n (t), q k (t) and v k (t). This evolution is solved from time t = 0 to t max with a given step size ∆t. Typically, t max is set to reach the stabilisation of the signal inside the tube. In continuous waveform (CW) regime, the wave is excited (at prime cell n = 1) at the frequency F by setting V n=1 (t) = U cos(2πF t), with U depending on the input power. To simulate multicarriers, we simply set V n=1 (t) = χ U χ cos(2πF χ t), , mvtrad (dashed lines) and measurements (dotted red line) from a 140W Ku-band tapered helix TWT [START_REF] Minenna | Proc. IEEE Int. Vacuum Electron. Conf. (IVEC)[END_REF]. Output power at saturation versus frequency. The cathode potential V0 of mvtrad is shifted by a small V h (between 0 to -100 V) to adjust either the power output or the phase shift but not both at the same time. The phase shift of dimoha is identical to measurements (see Fig. 2).

with the wanted frequencies and powers. A major advantage of dimoha is that we can inject any wave of interest, e.g. telecom signals with phase shift keying. From our N -body self-consistent hamiltonian model [START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF], we built a symplectic integrator solving (7)- [START_REF] Minenna | Electronwave momentum exchange and time domain simulations applied to traveling wave tube[END_REF]. A symplectic transformation [START_REF] Arnold | Mathematical Methods of Classical Mechanics[END_REF], usually associated with Hamilton's principle, preserves areas of the phase space. For numerical simulations, it allows [START_REF] Hairer | Geometric numerical integration[END_REF] the increase of step size ∆t without inducing too large errors. The algorithm is parallelised to employ multi-processors. This parallelisation is done by dividing spatially the tube between the allocated threads, associating each to a spatial domain with corresponding z thread cell and z thread field . During the time evolution, those threads compute amplitudes V n (t), I n (t) for n ∈ z thread cell and particle motions (coupling + space charge effects) for Q k ∈ z thread field . The last part of dimoha computes the output power, the frequency spectrum of the output field and the distribution functions (velocities, positions) of particles. The electromagnetic power (used in the next section) for a given frequency F is [START_REF] Pierce | Traveling Wave Tubes[END_REF]

P z = | Êz | 2 /(2β 2 Z c ),
with Êz the time-Fourier transform of Eq. (1) over a period 1/F .

V. NUMERICAL RESULTS WITH A HELIX TWT

In this section, we perform several comparisons between dimoha and measurements from a commercial 140W Ku-band pitch tapered helix TWT.

A. Power and phase shift

We first compare dimoha with the output power and phase shift from a commercial 140W Ku-band tapered helix TWT [START_REF] Minenna | Proc. IEEE Int. Vacuum Electron. Conf. (IVEC)[END_REF]. Figure 2 displays the RF output power (AM/AM) and the nonlinear phase shift (AM/PM) at various frequencies. The phase shift in this context is the difference of phase between the input power P in of interest and the input power at -20 dBm under the saturation P in,sat . The horizontal axis reads 10 log 10 (P in /P in,sat ) in dB. Small variations between dimoha and measurements may result from our incomplete management of tube defects (like hot voltage standing-wave ratio VSWR, defect in the beam confinement, and 3D asymmetry), from the beam model or from the violation of periodicity.

With initial spacing δ = 10 -6 m, the model involves approximately 150 000 macro-electrons, leading to a total of 150 000 + 200 dofs. We set t max = 5 ns with step size ∆t = 334 fs. The algorithm takes less than 10 min to reach steady state on a small desktop computer (8 threads with one Intel R Xeon R E5640 at 2.66GHz) and 2 min (32 threads with one Intel R Xeon R Gold 6142 at 2.6GHz) on a cluster [28], with δ and ∆t small enough to ensure the convergence of the power to ±0.1 W. The parallelisation ensures a good scaling up to 8 processors. This is a significant improvement over the 12 hours needed [START_REF] Safi | A Traveling-Wave Tube Simulation Approach With CST Particle Studio[END_REF] with cst for a similar tube.

B. Comparison against a frequency code

Figure 3 displays the saturated output power over the complete tunable range of frequencies for the same Ku-band TWT and shows the comparison with mvtrad (an industrial standard for TWT design). As shown, the cathode potential V 0 of mvtrad is shifted (between 0 to -100 V) to adjust either the power output or the phase shift but not both at the same time. dimoha is almost overlapping with measurements, showing a better accuracy than mvtrad. Ongoing investigations with different TWTs support our claim that dimoha is more accurate than mvtrad. Again, we surmise that small variations of dimoha versus measurements stem from the management of tube defects.

C. Multiple carriers dimoha's prime advantage is that it considers the complete fields in time (not just wave envelopes) as well as the electron dynamics (positions and velocities). This enables our algorithm to work in nonlinear regime (trapping, chaos) and to inject any wanted type of wave, like multiple carriers and telecommunication signals. This is achieved by simply modifying V n=1 (t) at the input cell.

We consider a simple two-carrier regime, setting V n=1 (t) = U cos(2πF 1 t) + U cos(2πF 2 t), with F 1 = 11.65 GHz and F 2 = 11.75 GHz, and with the adequate U ensuring P in to be -10 dBm. Figure 4 displays the frequency spectrum at the same Ku-band TWT output. The power at each frequency is computed from the time-Fourier transform of Eq. ( 1) over t max = 100 ns (minus 5 ns to reach steady state). All intermodulation products are reproduced near the two original carrier waves and their second harmonics. For second harmonics (lower plot), we do not have an accurate calibration of the bench and this TWT is not adapted in this band. The computation time is about 5 h for 100 ns on a small desktop computer with the step size ∆t = 83 fs. Further investigations are in progress for multiple carrier and digital modulations.

VI. FOLDED WAVEGUIDES

Since the theory is valid for any periodic structure, the algorithm is not limited to helix SWS but should handle all slow-wave structures. We propose an example. Ref. [START_REF] Karetnikova | Gain Analysis of a 0.2-THz Traveling-Wave Tube With Sheet Electron Beam and Staggered Grating Slow Wave Struc-ture[END_REF] provides inputs to model their sub-THz staggered double-grating TWT and includes results, used in Fig. 5, from 3D electromagnetic PIC codes cst [START_REF]CST Particle Studio[END_REF] and karat [START_REF] Tarakanov | User's Manual for Code KARAT[END_REF], and their 1D frequency model [START_REF] Karetnikova | Gain Analysis of a 0.2-THz Traveling-Wave Tube With Sheet Electron Beam and Staggered Grating Slow Wave Struc-ture[END_REF]. Figure 5 displays the comparison between those approaches and dimoha in small signal regime. A truncation of coupling terms Ω n to the nearest neighbour is enough.

The agreement with cst in the middle of the band (200-240 GHz) is great for a computation time about 23 min on a small desktop computer with the step size ∆t = 18 fs. The band edge discrepancies may be due to various reasons. To compute the coupling terms Ω n , the discrete model requires the dispersion relation be- 1D code (from [START_REF] Karetnikova | Gain Analysis of a 0.2-THz Traveling-Wave Tube With Sheet Electron Beam and Staggered Grating Slow Wave Struc-ture[END_REF]) CST (from [START_REF] Karetnikova | Gain Analysis of a 0.2-THz Traveling-Wave Tube With Sheet Electron Beam and Staggered Grating Slow Wave Struc-ture[END_REF]) CST (from [START_REF] Karetnikova | Gain Analysis of a 0.2-THz Traveling-Wave Tube With Sheet Electron Beam and Staggered Grating Slow Wave Struc-ture[END_REF]) KARAT (from [START_REF] Karetnikova | Gain Analysis of a 0.2-THz Traveling-Wave Tube With Sheet Electron Beam and Staggered Grating Slow Wave Struc-ture[END_REF]) DIMOHA Figure 5. Comparison of dimoha (squares), and simulations from Ref. [START_REF] Karetnikova | Gain Analysis of a 0.2-THz Traveling-Wave Tube With Sheet Electron Beam and Staggered Grating Slow Wave Struc-ture[END_REF] using cst (circles and diamonds), karat (triangles) and their 1-D code (dashed blue lines), from a sub-THz folded waveguide TWT with a staggered double-grating SWS [START_REF] Karetnikova | Gain Analysis of a 0.2-THz Traveling-Wave Tube With Sheet Electron Beam and Staggered Grating Slow Wave Struc-ture[END_REF] in small signal regime. The agreement with cst in the middle of the band (200-240 GHz) is remarkable.

yond the transmission band to be more accurate. Therefore, the data reading from the dispersion relation from Ref. [START_REF] Karetnikova | Gain Analysis of a 0.2-THz Traveling-Wave Tube With Sheet Electron Beam and Staggered Grating Slow Wave Struc-ture[END_REF] (limited from 177 to 260 GHz) is sensitive to any deviation, especially near the band edge. Moreover, space charge models are different (Ref. [START_REF] Karetnikova | Gain Analysis of a 0.2-THz Traveling-Wave Tube With Sheet Electron Beam and Staggered Grating Slow Wave Struc-ture[END_REF] considers a rectangular beam and we simulate a cylindrical beam). However, experimental folded waveguides [START_REF] Hu | Design and Experimental Study of a Widebandwidth W-Band Folded Waveguide Continuous-Wave TWT[END_REF] exhibit no peaks near edges. The comparison on Fig. 5 is not meant to explain those discrepancies nor to provide a full analysis of folded waveguides. It rather shows that dimoha is highly flexible [31] and it can handle various geometries in a fraction of the time required for PIC codes. In addition, ongoing analyses are showing an excellent agreement between dimoha and cst for a W-band (92-95 GHz) folded waveguide TWT [START_REF] André | Fabrication of W-band TWT for 5G small cells backhaul[END_REF] with a sever.

VII. CONCLUSIONS AND PERSPECTIVES

dimoha is a promising tool, for both industrial and research activities. It also provides a new approach to analyse telecommunications signals for operators. Discrepancies exist between measurements and dimoha. A more complete management of defects and a better space charge model could reduce this. However, our algorithm is more accurate than the frequency domain code mvtrad and is greatly faster than classical time domain PIC codes.

Our algorithm bears five main advantages. It can simulate complex signals like multiple carriers. It is fast and robust. It supports both helix and folded waveguide designs and it should be able to reproduce other kinds of SWS. It reproduces harmonic generation, distortion phenomena and gain. Finally, it can simulate nonlinear dynamics [START_REF] Minenna | Electromagnetic power and momentum in N-body Hamiltonian approach to wave-particle dynamics in a pe-riodic structure[END_REF].

We are investigating the cases of multiple carriers and digital modulations. A more accurate version for folded waveguides and a 2.5D or 3D version of dimoha are under consideration.

Figure 1 .

 1 Figure 1. Process workflow of dimoha. The algorithm is divided in three parts: (i) the pre-processing to compute cold-field parameters; (ii) the time integrator; (iii) the postprocessing (not shown) that performs diagnostics and computes the RF power. Inputs and outputs are listed in TableI.

Inputs

  Tube : Pitch d, dispersion relation ω(β), (per cell) impedance Zc(β) and attenuations αn Beam : Current I 0 and potential V 0 (cathode) beam diameter b Signal: Injected 1D electric field in time at the tube inlet Simulation parameters tmax Duration of the interaction ∆t Time step size δ Initial spacing between macro-electrons Nos Number of mesh points per cell N ph Range of coupling between cells Outputs 1D electric field in time inside the tube 1D distribution function of electrons (positions, velocities)

Figure 2 .

 2 Figure 2. Comparison of dimoha (solid lines) and measurements from a Ku-band tapered helix TWT (dashed lines). a) Output power (AM/AM) depending on frequency and input power. b) Nonlinear phase shift (AM/PM) depending on frequency and input power. Our attenuator model is set for the frequency F = 11.7 GHz.

Figure 3 .

 3 Figure3. Comparison of dimoha (solid black line), mvtrad (dashed lines) and measurements (dotted red line) from a 140W Ku-band tapered helix TWT[START_REF] Minenna | Proc. IEEE Int. Vacuum Electron. Conf. (IVEC)[END_REF]. Output power at saturation versus frequency. The cathode potential V0 of mvtrad is shifted by a small V h (between 0 to -100 V) to adjust either the power output or the phase shift but not both at the same time. The phase shift of dimoha is identical to measurements (see Fig.2).

Figure 4 .

 4 Figure 4. Comparison of dimoha (solid black lines) and measurements (dotted red lines) from a 140W Ku-band tapered helix TWT [11]. Output spectrum in two-carriers regime of dimoha (after tmax = 100 ns) and the TWT. Input: F1 = 11.65 GHz and F2 = 11.75 GHz at Pin = -10 dBm. Upper plot: fundamental signal. Lower plot: second harmonic.

Table I .

 I Inputs and outputs of dimoha

ACKNOWLEDGMENTS

The authors thank Lénaïc Couëdel, Guillaume Fuhr, Magali Muraglia, Valentin Pigeon and Geoffroy Soubercaze-Pun for fruitful discussions, Élise Duverdier and Telma Pereira for their help and anonymous reviewers for constructive comments. This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements d'Avenir" supervised by the Agence Nationale de la Recherche.


Abdelaziz Mouhou
email: abdelaziz.mouhou@gmail.com

Abdelmajid Badri

Abdelhakim Ballouk

Genetic Algorithm based Model Predictive Control of Fractional Order Systems

Keywords: Model Predictive Control, Fractional Order Systems, Genetic Algorithm

Fractional calculus is a generalization of integration and derivative operators to no integer order. In the Last decades fractional systems receives a great attention in the research community, the reason lies in the fact that many real systems can be described accurately using dynamic models of fractional order. In our work a design of genetic algorithm based fractional order model predictive controller is simulated and applied to a given fractional order model plant.

INTRODUCTION

Fractional calculus (FC) was found for the first time more than 300 years, it is a field of theoretical mathematics that deals with non-integer order integrals and derivatives. Since it was applied in different fields of engineering and applied physics, such as viscoelasticity, heat diffusion, dielectric polarization, electromagnetic waves, Schrodinger equation, colored noise [START_REF] Journal | SURVEY PAPER A review and evaluation of numerical tools for fractional calculus and fractional order controls Zhuo Li[END_REF]. In the last decades the control of fractional order systems is the same thing as integer order ones, where researchers use an integer approximation of FO systems to synthesis the adequate controller. However with the evolution of technologies and to increase the performances of controlled FO systems we notice the appearance of FO controllers such as fractional PID, fractional order model predictive control (FOMPC).

Predictive control was introduced by J. Richalet in 1978 and 1979 by D.W. Clarke and Gawthrup. Generalized predictive control (GPC) has been proposed by D.W. Clarke in 1987, since it was applied in many industrials processes, the aim of this strategy is to use a model in order to predict the future behavior of the plant [START_REF] Mouhou | Genetic Algorithms Optimization of Tuning Parameters of Generalized Predictive Control[END_REF].

In this paper we present the application of GPC to control a given fractional non commensurate order system. The remainder of the paper is organized as follows the formulation of adaptive generalized predictive control is presented in Section 2, the optimization of GPC parameters using genetic algorithms is detailed in Section 3, simulation results by applying GA based GPC to control given fractional order system are presented in Section 4, the conclusion is reported in the last section.

FRACTIONAL ORDER SYSTEMS

Grunwald-Letnikov fractional derivative operator

There are several fractional derivatives operators that is used in the literature, the most popular one used for discrete-time systems is the Grunwald-Letnikov deriva-tive operator.It is introduced by A. K. Grunwald 1867 and A.V. Letnikov 1868. It allows the generalization of non-integer order derivative [START_REF] Guermah | Discrete-Time Fractional-Order Systems : Modeling and Fractional-Order Stability Issues Systems : Modeling and Stability Issues[END_REF].

0 0 0 0 ( ) lim ( 1) ( ) t t GL j h t t j h D f t h f t jh j                        (1)
Where the so called generalized Newton symbol is given by:

1 0 ( 1)...( 1) 1, 2,3,.... ! for j j for j j j                      (2)

Oustaloup approximation of fractional order systems

For a given fractional order differentiator salpha the oustaloup arrpoximation for a given franquency boundries omega and omegab is defined as follows:

1 ' () 1 N ui iN h i s Gs s              (3)
Where

11 2 2 2 2 2 1 2 1 '; i N i N NN uu i b i b bb                                , bu
 are frequency boundaries. N is the order of Oustaloup approximation.

MODEL PREDICTIVE CONTROL

The aim of GPC is to create an anticipatory effect in order to provide a stable corrected system. The CARIMA (Controlled Autoregressive Integrated Moving Average) model is considered: q  related to disturbances. The GPC criterion to be minimized is represented by a quadratic cost function considering the tracking error and the control effort on a receding horizon as follows:

1 1 1 1 () ( ) ( ) ( ) ( ) ( ) () d k A q y k q B q u k C q q         (4)

  k are respectively the input and the output of the system.  random centered signal uncorrelated with the input.

1

 1

GA optimization of synthesis parameters

The choice of the fitness function is very important, it allows the convergence of the algorithm to the optimal solution.

  it is clear from the previous figures that the controlled system presents oscillations. To increase the controller performances we optimize the synthesis parameters using GA MATLAB toolbox. From the results illustrated in Fig4 ad Fig5 it is clear that the reduced approximated model of FO system presents no overshoot and the settling time and rise time are minimized.

Use MPC with the FO model Calculate the model output Get the elements of fitness function Begin Initialize Pop (0) While k<max_gen do Evaluate Pop(k) using fitness function Select New_Pop from Pop(k) Applying genetic operators (cross & mutation) to

CONCLUSION

This paper gives the application of optimized GPC using GA to control fractional order system approximated with Oustaloop method. However the Oustaloup approximation is a high order filter we used the SPA and BT algorithms to reduce the order of the obtained model. The application of GA based GPC of a noncommensurate fractional order model show the effectiveness of this method.