Arthur Degasperi

Lucie Labied

Carole Farre

Emmanuel Moreau

Matteo Martini

Carole Chaix

Agnès Hagège
email: agnes.hagege@isa-lyon.fr

Probing the protein corona of gold/silica nanoparticles by Taylor Dispersion Analysis-ICP-MS

Keywords: Taylor dispersion analysis, Inductively coupled plasma mass spectrometry, Metalcontaining nanoparticles, Protein corona, Reversible interactions, Albumin

Despite the tremendous interest for nanoparticles (NPs) in the biomedical field, their transfer to the clinics is still hampered, in particular due to the lack of knowledge of their behaviour in a biological environment. Indeed, the protein corona formed as soon as NPs enter the bloodstream can drastically affect their properties. The use of Taylor dispersion analysis-ICP-MS as an efficient technique dedicated to metal-containing NPs was proposed to examine these NP-protein interactions and determine protein corona thicknesses in biological fluids. This method was applied on core-shell gold/silica NPs in the presence of proteins at high concentrations and serum. Protein corona around 4 nm were measured. Moreover, the versatility of the method allowed assessing the reversible/irreversible character of the interactions.

1.

Introduction

Metal-containing nanoparticles (NPs) are increasingly investigated in the biomedical field. This burgeoning technology enables to explore various alternative strategies, in the field of both diagnosis and therapy [START_REF] Huang | Inorganic nanoparticles in clinical trials and translations[END_REF][START_REF] Park | Advanced hybrid nanomaterials for biomedical applications[END_REF]. However, the development of effective therapeutic strategies is based on the understanding of the in vivo behavior of NPs. As they enter the bloodstream, NPs may be rapidly surrounded by a wide variety of biomolecules, in particular proteins, which adsorb on their surface and may completely cover them. This protein corona modifies the physicochemical characteristics of the pristine NPs and constitutes the first contact of the nanomaterial with the cells. It is therefore susceptible not only to modify the absorption pathways of NPs but also to generate an additional level of complexity [START_REF] De Crozals | Nanoparticles with multiple properties for biomedical applications: A strategic guide[END_REF][START_REF] Lesniak | Effects of the Presence or Absence of a Protein Corona on Silica Nanoparticle Uptake and Impact on Cells[END_REF][START_REF] Singh | In vivo protein corona on nanoparticles: does the control of all material parameters orient the biological behaviour?[END_REF]. Milani et al. demonstrated that proteins involved in this corona exhibit very distinct kinetics of dissociation. The so-called "hard corona" is assumed to be irreversible while the "soft corona" shows fast off-rates (in the range of a minute) [START_REF] Milani | Reversible versus Irreversible Binding of Transferrin to Polystyrene Nanoparticles: Soft and Hard Corona[END_REF]. Due to the difficulty to preserve the soft corona, its potential occurrence has been often overlooked. Apart from a few studies reporting that soft corona proteins may contribute to the stealth properties of nanoparticles [START_REF] Kari | In situ analysis of liposome hard and soft protein corona structure and composition in a single label-free workflow[END_REF] or affect the attachment efficiency of NPs to biological membranes [START_REF] Wang | Influence of Solution Chemistry and Soft Protein Coronas on the Interactions of Silver Nanoparticles with Model Biological Membranes[END_REF], the role of this soft corona remains elusive.

Much research has been performed to identify the proteins forming the corona by mass spectrometry [START_REF] Cedervall | Understanding the nanoparticle-protein corona using methods to quantify exchange rates and affinities of proteins for nanoparticles[END_REF][START_REF] Kurtz-Chalot | Impact of silica nanoparticle surface chemistry on protein corona formation and consequential interactions with biological cells[END_REF][START_REF] Fernández-Iglesias | Complementary mass spectrometric techniques for the quantification of the protein corona: a case study on gold nanoparticles and human serum proteins[END_REF]. However, the usual procedures imply the previous separation of the NPs from the biological fluids, leading then to the loss of proteins with fast dissociation kinetics.

The number of proteins involved in the hard corona of gold NPs was determined by laser induced breakdown spectroscopy, using the plasmonic features of gold NPs [START_REF] Dell'aglio | Sensing nanoparticle-protein corona using nanoparticle enhanced Laser Induced Breakdown Spectroscopy signal enhancement[END_REF]. An alternative approach used capillary-ICP-MS hyphenation to separate free NPs from Au-NP-protein conjugates. It demonstrated the replacement of initially adsorbed transferrin on the particle surface by albumin with time [START_REF] Matczuk | Speciation of metal-based nanomaterials in human serum characterized by capillary electrophoresis coupled to ICP-MS: a case study of gold nanoparticles[END_REF].

Moreover, the determination of the number of albumin involved in the conjugate was realized by quantifying Au, referring to the NP, and S referring to the protein [START_REF] Matczuk | Characterization of the protein corona of gold nanoparticles by an advanced treatment of CE-ICP-MS data[END_REF].

Characterization of protein adsorption can also be performed through NP size measurements. When NPs are suspended in biological media, protein adsorption on their surface results in an increased size of the NPs, due to protein adsorption. Therefore, corona thickness can be deduced from NP sizes measured in media either free or containing proteins. However, measuring a protein corona of a few nanometers is intricate. Only few techniques do not require sample treatment before analysis, so as not to disturb the equilibrium of the protein corona and its properties. The most popular technique for NPs size measurements is dynamic light scattering (DLS). However, the scattering intensity of NPs in DLS is often much smaller as compared to that of serum [START_REF] Fischer | Pitfalls and novel applications of particle sizing by dynamic light scattering[END_REF]. For that reason, DLS measurements are often conducted in diluted media. Due to the low contrast generated by proteins, transmission electron microscopy (TEM) requires the addition of a contrast agent (generally uranyl acetate), which precludes its use in situ [START_REF] Moya | Insights into Preformed Human Serum Albumin Corona on Iron Oxide Nanoparticles: Structure, Effect of Particle Size, Impact on MRI Efficiency, and Metabolization[END_REF]. Therefore, alternatives methods were developed. In a recent work, Galdino et al. studied the process of protein corona formation on silica NPs by core-loss energy-filtered TEM and cryo-TEM. They successfully measured corona thicknesses of 4 -9 nm. Small-angle X-ray scattering (SAXS) was also used to determine binding parameters in the presence of both BSA and serum [START_REF] Galdino | Inside the Protein Corona: From Binding Parameters to Unstained Hard and Soft Coronas Visualization[END_REF]. Although very sensitive, SAXS revealed insufficiently accurate measurements of the protein corona thickness. Nevertheless, it seems to be the only technique reported to distinguish between protein corona and aggregation. In the case of fluorescent NPs, fluorescence correlation spectroscopy (FCS) was used for determination of the protein corona thickness in media containing proteins with concentrations close to the ones encountered in serum. Corona thicknesses of a few nanometers were successfully measured on a wide variety of NPs [START_REF] Wang | The protein corona on nanoparticles as viewed from a nanoparticle-sizing perspective[END_REF][START_REF] Treuel | Impact of protein modification on the protein corona on nanoparticles and nanoparticle-cell interactions[END_REF][START_REF] Maffre | Characterization of protein adsorption onto FePt nanoparticles using dual-focus fluorescence correlation spectroscopy[END_REF].

We recently showed that Taylor dispersion analysis hyphenation to inductively coupled plasma mass spectrometry (TDA-ICP-MS) is particularly adapted for metal-containing NP size measurements [START_REF] Labied | Taylor dispersion analysis coupled to ICP-MS for ultrasmall nano-particle size measurement: from drug product to biological media studies[END_REF][START_REF] Labied | Biodegradation of metal-based ultra-small nanoparticles: a combined approach using TDA-ICP-MS and CE-ICP-MS[END_REF]. Furthermore, thanks to the selectivity of the technique to metals, TDA-ICP-MS allows size measurements to be performed in complex matrices, rich in proteins, whereas these measurements are often not possible by DLS. Herein, we devised a methodology based on the monitoring of the metal signal, which enables not only to measure the protein corona thickness onto NPs, but also to gain information about the reversibility of this adsorption. As a proof-of-principle, the method was applied using 60 nm core-shell gold/silica nanoparticles (Au@SiOx NPs) suspended in protein-rich model solutions and serum.

2.

Materials and methods

Chemicals and solutions

Triton X-100, n-hexanol, cyclohexane, HAuCl4

Synthesis of the Au@SiOx NPs

Au@SiOx NPs were prepared according to a method developed by Martini et al. [START_REF] Martini | How gold particles suppress concentration quenching of fluorophores encapsulated in silica beads[END_REF]. Quaternary W/O microemulsion was prepared by mixing 11.5 mL of Triton X-100 (surfactant), 11.5 mL of n-hexanol (co-surfactant), 48 mL of cyclohexane (oil). The water phase consists of mixture of 3 mL of 16.7 mmol.L -1

HAuCl4⋅3H2O, 3 ml of 32.8 mmol.L -1 MES and 1 mL of 412 mmol.L -1 NaBH4 in order to synthesize the gold nanoparticle. After 5 min, 0.050 mL of APTES was added to the microemulsion solution, followed by 0.2 mL of TEOS. After 30 min, the silica polymerization reaction was completed by adding 0.2 mL of NH4OH and then stirring for 24 h. In order to achieve the colloidal stability, a carboxylated protective shell was created by adding 0.02 mL of Si-20. The microemulsion was broken by the addition of ethanol to the microemulsion and followed by vortexing and centrifuging several times. Unreacted precursors were removed by ultrafiltration using 100 kDa polyethersulfone membranes. Particles were dispersed in aqueous solution (2 mg.mL -1 ) and stored at 4°C. Characterization of NPs is provided as Supplementary material (section 1).

Taylor dispersion analysis -ICP/MS analysis

A 7100 capillary electrophoresis system (Agilent) was used for TDA analyses.

Detection of gold was enabled by hyphenation with a Nexion 300X ICP-MS (Perkin Elmer). For that purpose, a T-connector was used as sheath-flow interface. The TDA capillary went through a Tconnector toward the tip of the nebulizer and the sheath liquid (HNO3 1 mmol.L -1 ) was introduced by self-aspiration to the third inlet of the interface. This interface was previously shown not to introduce any external dispersion that could enlarge peaks [START_REF] Labied | Taylor dispersion analysis coupled to ICP-MS for ultrasmall nano-particle size measurement: from drug product to biological media studies[END_REF].

Au signal was monitored at m/z 197, with a dwell time of 500 ms. Operating conditions used for the ICP-MS were: nebulizer gas flow rate, 1.09 L/min; plasma gas flow rate, 18 L/min; auxiliary gas flow rate, 1.2 L/min; radiofrequency power, 1600 W for the plasma. All other parameters were tuned to maximize the Au signal.

3D-CE Chemstation software version B.04.03 (Agilent) was used to control the CE system and Syngistix software version 2.3 was used to both control the ICP-MS and acquire the 197 Au signal.

Using sodium phosphate mobilization solutions, capillaries were conditioned 10 min at 1 bar. For protein mobilization media, capillaries were conditioned 5 min at 1 bar with the protein solution, then 10 min at 50 mbar.

One g.L -1 Au@SiOx NP suspensions were prepared by addition of an equivalent volume of the incubation medium (sodium phosphate buffer 10 mM, protein solutions, and serum). The Au@SiOx NP suspensions were kept at room temperature without stirring and aliquots were collected at different times for size measurements.

NPs were injected hydrodynamically under 40 mbar for 3 s in a capillary filled with the desired medium and a mobilization pressure of 50 mbar was applied.

Measurements of hydrodynamic radii in the incubation medium were performed in capillaries filled with the same medium. Assessment of the reversible character of the protein interactions were performed in capillaries filled with 10 mM phosphate buffer, pH 7.4.

Calculation of the nanoparticle hydrodynamic radius

The detected peak can be fitted by a sum of Gaussian distributions, according to the following equation:

𝑆(𝑡) = ∑ ! ! " ! √$% exp +- ('(' " ) # $" ! # - * +,- (1)
where i is the number of species, t is the experimental time, t0, the residence time, σi the mean temporal variance of the Gaussian distribution associated to the species i and Ai, a constant depending on the concentration of the species i.

Under experimental conditions described by Chamieh et al. [START_REF] Chamieh | Limits in Size of Taylor Dispersion Analysis: Representation of the Different Hydrodynamic Regimes and Application to the Size-Characterization of Cubosomes[END_REF], the molecular diffusion coefficient D can be obtained using the Stokes-Einstein equation:

𝐷 = . $ # ' " $/" # = 0 % 1 2%3. & (2)
where Rc is the capillary radius, kB is the Boltzmann constant, T is the temperature, η is the viscosity, and Rh is the hydrodynamic radius of the solute.

The diffusion coefficient, and therefore the hydrodynamic radius of a solute, are then based on the Gaussian fits of the peaks obtained in TDA-ICP-MS realized using the Origin software (version 8.5).

As a tailing deformation induced by retention may occur at the end of the peak, only the first half of the peak can be assimilated to a sum of Gaussians. In theory, a peak half contains all the necessary and sufficient information. However, the noise at the top of the peak can bias the determination of the residence time. Therefore, the time interval over which the adjustment must be made was carefully determined in order to limit both noise and retention effects. Preliminary experiments have shown that 15 to 30% of the second half of the peak must be taken into account.

DLS measurements

Colloidal solutions of Au@SiOx NPs at 2 mg.mL -1 in pure water were used. Measurements were performed using a Zetasizer Nano-S (633 nm He-Ne laser) from Malvern Instruments (USA).

Measurements were averaged over 10 individual measurements, performed at 25°C. The validity of the data was determined by examination of both the autocorrelation curves and the quality report.

Results and discussion

Hydrodynamic size of the nanoparticle

Hydrodynamic size measurement of the Au@SiOx NPs suspended in aqueous solution was performed using both TDA-ICP-MS and DLS (Figure 1). The hydrodynamic radius obtained with TDA-ICP-MS in 10 mmol.L -1 sodium phosphate buffer was found to be (30 ± 1) nm, over 12 repetitions (Fig. 1a). This result is in close agreement with the value obtained by DLS: (30 ± 2) nm (Fig. 1b). The goodness of the fit strongly supports the assumption of the main monodisperse character of the NP, as confirmed by TEM and DLS (see Supplementary material, section 1).

Measuring the protein corona thickness

TDA-ICP-MS was already demonstrated to successfully measure hydrodynamic radii of NPs in various complex media [START_REF] Labied | Biodegradation of metal-based ultra-small nanoparticles: a combined approach using TDA-ICP-MS and CE-ICP-MS[END_REF]. The selectivity of the ICP-MS detection, in our case 197 Au, allows the only detection of Au compounds and eliminates interference caused by other compounds such as free proteins present in the solution. Thanks to this features, this coupling allows to consider the measurement of the size of NPs in media with high protein contents. For that purpose, two model proteins have been chosen: BSA, which is one of the most studied protein and FET. Moreover, both proteins are important components of many cell culture media (i.e. media containing fetal calf serum). The isoelectric point of FET (pI=5.4) is higher than the one of BSA (pI = 4.7). Both proteins have close molecular weight.

However, in serum, FET concentration (0.4 -0.9 mg.mL -1 ) is rather low compared to that of albumin (35 -55 mg.mL -1 ), although it is one of the 15 most abundant proteins in serum.

The influence of the concentration of both BSA and FET was investigated. For that purpose, colloidal solutions of Au@SiOx NPs were incubated at room temperature for at least 24 h with increasing concentrations of both FET and BSA (0 to 30 and 50 mg.mL -1 , respectively) in sodium phosphate 10 mmol.L -1 , pH 7.4. Prior to incubation, the protein solutions were checked by TDA, using UV detection, in order to ensure the non-aggregation of proteins (see Supplementary material, section 2)

Aliquots were frequently collected (25 µL) for at least 6 h and analysed in a capillary, filled with the same medium than the incubation medium, to ensure the equilibrium between the NPs and the proteins.

Results are presented in Figures 2 and3 for BSA and FET, respectively. For each condition, Taylorgrams exhibit a peak, which can be fitted using a single Gaussian. Moreover, the peak areas were found to be constant over the timescale of the experiments, for all the protein concentrations, indicating that the colloidal stability of the NPs is maintained in such media.

The observed variations in average peak residence times as a function of protein concentration are due to the different viscosities of the studied media. As for DLS, the diffusion coefficient, and hence the hydrodynamic radius calculation, depends on viscosity. For that purpose, viscosities were measured in the different media (see Supplementary material, section 3) and used to determine hydrodynamic radii.

Whatever the protein concentration, a rapid increase of the hydrodynamic radius (in less than 20 min) was observed for both proteins. This increase was attributed to the formation of a protein corona around the particle. Values are presented in Table 1. The significance of differences between the data were determined based on P values < 0.001 by using t-tests (as described in Supplementary Material, section 4) Consequently, protein corona thicknesses can be obtained from the difference DRH, between the hydrodynamic radius obtained in a protein medium (either BSA or FET) and the one obtained in phosphate buffer, recapitulated in Table 1. For both proteins, an increase in hydrodynamic radius was observed.

The isoelectric point of the Au@SiOx NPs surface, obtained by measuring the zeta potential z as a function of pH (see Supplementary material, section 1), was about 4.3, confirming the main presence of carboxylate functions on surface. This value is higher than that conventionally obtained for silica NPs (pI = 3.5) [START_REF] Pacard | Elaboration of silica colloid/polymer hybrid support for oligonucleotide synthesis[END_REF]. This difference may be due to the additional presence of both carboxylated groups and ammonium ions originating from APTES on the NP surface [START_REF] Rouet | Colloidal chemistry with patchy silica nanoparticles[END_REF].

Owing to the pI of both proteins, the question of a heterogeneous spatial distribution of charges should be evoked to explain the protein adsorption onto the NPs, as both proteins and AuSiOx NPs are negatively charged at pH 7.4.

In the case of BSA, the values obtained for the corona thickness (DRH ~ 3 -7 nm) seem to be in agreement with the values reported in the literature. The formation of a BSA corona on silica NPs was studied by other groups and seems to be size-dependent. A corona thickness of 7 nm and 4 nm was reported by Galdino et al. for 40 nm and 80 nm -silica NPs respectively [START_REF] Galdino | Inside the Protein Corona: From Binding Parameters to Unstained Hard and Soft Coronas Visualization[END_REF]. Monopoli et al. observed a corona thickness of 4 nm for 200 nm silica NPs [START_REF] Monopoli | Physical-Chemical Aspects of Protein Corona: Relevance to in Vitro and in Vivo Biological Impacts of Nanoparticles[END_REF]. Smaller particles are thought to be more prone to aggregate in the presence of BSA. Orts-Gil et al. observed an hydrodynamic radius of ca. 100 nm for 15 nm-silica NPs in BSA solutions, resulting from the formation of NPs-BSA clusters [START_REF] Orts-Gil | On the role of surface composition and curvature on biointerface formation and colloidal stability of nanoparticles in a protein-rich model system[END_REF] and a small aggregation was also reported by Monopoli et al. for 50 nm -silica NPs [START_REF] Monopoli | Physical-Chemical Aspects of Protein Corona: Relevance to in Vitro and in Vivo Biological Impacts of Nanoparticles[END_REF].

A corona thickness of about 3-4 nm for the corona in the presence serum albumin was also observed for other NPs, such as 5 nm-FePt NPs [START_REF] Maffre | Characterization of protein adsorption onto FePt nanoparticles using dual-focus fluorescence correlation spectroscopy[END_REF], 5 nm-quantum dots [START_REF] Treuel | Impact of protein modification on the protein corona on nanoparticles and nanoparticle-cell interactions[END_REF] and even larger nanoparticles, i.e.

polystyrene nanoparticles between 20 and 100 nm [START_REF] Wang | Formation of a Monolayer Protein Corona around Polystyrene Nanoparticles and Implications for Nanoparticle Agglomeration[END_REF].

Given the hydrodynamic radius of BSA (4.3 ± 0.3 nm, see Supplementary material, section 2), this 4 nm-corona may reflect the formation of a BSA monolayer, also indicating the absence of multilayer adsorption. . However, the accuracy of TDA-ICP-MS does not allow us a more thorough interpretation of the recorded DRH until now.

On the contrary, FET adsorption seems to be concentration-dependent. For a FET concentration close to the one of serum (0.8 mg.mL -1 ), no significant adsorption was observed. Increasing concentration in the 5 -15 mg.mL -1 range yielded a corona thickness of ~3 -5 nm, supporting the assumption of a monolayer adsorption. Indeed, the hydrodynamic radius of FET, measured by TDA -UV (see Supplementary material, section 2) was found to be (4.6 ± 0.4) nm.

For the highest concentration, similar to the one of BSA in serum, an unexpected increase in RH was observed. As can been seen in Figure 3b, a gradual increase of the hydrodynamic radius seems to occur upon time, until a plateau is reached for values approximately twice those of the Au@SiOx NP radius after 80 min incubation.

In the literature, this behaviour was also encountered for 40 nm-silica particles in fetal calf serum, for which a progressive increase in hydrodynamic radius until a value of 37-46 nm was observed by DLS.

Based on microscopy data, this was attributed by Galdino et al. to the formation of small aggregates [START_REF] Galdino | Inside the Protein Corona: From Binding Parameters to Unstained Hard and Soft Coronas Visualization[END_REF]. As TDA experiments showed that FET aggregation did not occur in phosphate buffer using such a concentration (see Supplementary material, section 2), NPs aggregation mediated by proteins might be considered, even if the formation of protein multilayers cannot be excluded in the absence of additional proofs. However, it is worth noting that we did not observe protein aggregates by TDA-UV, in the 30 and 60 mg.mL -1 FET solutions (see Supplementary material, section 2).

The method was then applied to the measurement of Au@SiOx NPs protein corona formed in human serum. Results are presented in Figure 4. In the presence of serum, the hydrodynamic radius of the Au@SiOx NPs also increases. Although the standard deviation is more important, as for BSA 50 mg.mL -1 and FET 30 mg.mL -1 , the values obtained for the longest incubation times ((35 ± 2) nm, (36 ± 3) nm, and (39 ± 4) nm, for 120, 150, and 255 min incubation respectively) are significantly higher than the hydrodynamic radius of the NPs in the absence of proteins. As a consequence, the existence of a protein corona of ca. 5 nm-thickness can be postulated.

Assessing the high binding affinity protein corona

It is now commonly accepted that the protein corona consists in proteins with high binding affinities and therefore low exchange rate, known as «hard corona», but also in a highly dynamic «soft corona», consisting of further adsorbed proteins with low binding affinities and high exchange rates. In order to evaluate the affinity of the proteins forming the corona, Au@SiOx NPs were incubated in different media (phosphate buffer, protein solutions, diluted serum) for 2 h and were then analysed in a capillary filled with either the incubation medium or the phosphate buffer pH 7.4. Additionally, experiments were also performed in the phosphate buffer after incubation of Au@SiOx NPs in a BSA solution 50 mg.mL -1 for 4 h. Results are presented in Figure 5. The comparison of the observed radii indicates whether proteins were adsorbed as either a "hard + soft corona" (blue bars in Figure 5) or a hard corona (grey bars in Figure 5).

As can be seen in Figure 5, analyses performed in the phosphate buffer do not allow observing a hard corona composed of BSA on Au@SiOx NPs after 2 h incubation. Indeed, the radius obtained after analysis in the phosphate buffer of the sample incubated in a BSA-containing medium ((30.0 ± 0.9) nm) is very close to the one observed for the Au@SiOx NPs suspended in a protein-free incubation medium (phosphate buffer). However, this reversibility seems to depend on the incubation time (see Supplementary material, section 5). Increasing the incubation time up to 4 h leads to the recovery of the protein corona thickness observed in the BSA medium (4-5 nm), which may reflect the irreversible character of the protein corona. It can then be assumed that the corona formed by BSA on Au@SiOx NPs after a 2 h -incubation is reversible and dissociates when the medium is switched to phosphate buffer, pH 7.4, whereas the corona formed after a longer incubation time remains present at the NP surface when switching to a free protein medium. Mc Ritchie [START_REF] Ritchie | The adsorption of proteins at the solid/liquid interface[END_REF] reported a desorption at pHs far from the BSA pI. Similarly, although Rezwan et al. [START_REF] Rezwan | Lysozyme and bovine serum albumin adsorption on uncoated silica and AlOOH-coated silica particles: the influence of positively and negatively charged oxide surface coatings[END_REF] detected the formation of a corona on 97 nmsilica NPs in suspension in a BSA solution, and no BSA was found to be adsorbed after elimination of free proteins by centrifugation. A reversible corona was also observed by dilution experiments of carboxyl-coated quantum dots in the presence of human serum albumin [START_REF] Wang | The Nature of a Hard Protein Corona Forming on Quantum Dots Exposed to Human Blood Serum[END_REF].On the contrary, Galdino In our case, we assume that the protein layer is nearly totally released upon dilution in phosphate buffer after 2 h incubation and that the corona formed by BSA on Au@SiOx NPs is reversible and dissociates when the medium is switched to phosphate buffer, pH 7.4. However, increasing the incubation time leads to an increase of the hydrodynamic radius, which may reflect the irreversible character of the protein corona. A possible explanation could be that irreversible adsorption of BSA requires structural changes in the adsorbed BSA, which takes more time to occur, as suggested in several papers [START_REF] Norde | Adsorption, desorption and re-adsorption of proteins on solid surfaces[END_REF][START_REF] Giacomelli | The Adsorption-Desorption Cycle. Reversibility of the BSA-Silica System[END_REF][START_REF] Kondo | Structural Changes in Protein Molecules Adsorbed on Ultrafine Silica Particles[END_REF].

By dilution in a phosphate medium, a remaining monolayer of FET was found to be irreversibly adsorbed. For both 15 mg.mL -1 and 30 mg.mL -1 FET solutions, the corona thickness was found to be around 5 nm, which makes unlikely the hypothesis of an irreversible aggregation of the particles.

Indeed, the experiment previously carried out with a capillary filled with 30 mg.mL -1 FET showed an increment in the apparent radius of the NPs around 32 nm, after 80 min incubation in the FET solution.

As previously stated, this size, observed by TDA-ICP-MS performed in the 30 mg.mL -1 FET incubation medium, corresponds either to the reversible interactions of protein aggregates induced by the NPs or to the reversible formation of small NP clusters. By the dilution of NPs in phosphate buffer, the increment in radius was then reduced to a 5 nm thickness. Upon dilution in phosphate, a hard corona of 5 nm thickness was still tethered around the NPs. To conclude from these experiments, we proved that FET has a higher binding affinity than BSA for the Au@SiOx NP surface. FET was also noticed to exhibit a higher affinity than BSA for gold NPs [START_REF] Moustaoui | A protein corona study by scattering correlation spectroscopy: a comparative study between spherical and urchin-shaped gold nanoparticles[END_REF].

The Au@SiOx NPs were then incubated in 50% serum for two hours before TDA-ICP-MS analysis. As evidenced by the absence of RH variations in diluted serum experiments using both conditions (50% serum and phosphate buffer), an irreversible adsorption seemed to occur. This behaviour cannot be explained by BSA adsorption and therefore, this very abundant protein does not seem to dictate the fate of Au@SiOx NPs in serum. However, the presence of FET in serum cannot explain such a corona.

Indeed, as shown in Table 1, neither a hard nor a soft corona was observed for FET concentrations representative of concentrations found in serum (0.8 mg.mL -1 ).

Conclusion

In this work, TDA-ICP-MS sizing experiments was shown to provide insight into corona formation of metal-containing NPs from only 25 µL sample volumes. NPs size measurements were achieved in media containing high concentrations of proteins, similar to serum. Variations of hydrodynamic radii of about 15 % can significantly be recorded. The values of protein corona thickness reflect the adsorption of a protein monolayer for both proteins on the Au@SiOx NPs surface. This method was shown to successfully overcome the limitations of DLS in such concentrated media. However, like DLS, TDA-ICP-MS is not able to distinguish an increase in size due to the formation of a large corona from the one due to aggregation without additional experiments.

By a simple experimental design, the reversible/irreversible character of the protein adsorption can be assessed. Experiments performed in non-equilibrium conditions show that BSA dissociates rapidly after a 2h-incubation , unlike BSA after a longer incubation, FET, and serum constituents for which the monolayer remains unchanged in size. This seems to indicate that BSA cannot be used as a representative model of serum for these Au@SiOx NPs, as its binding affinity was revealed to be lower than that of some of the serum constituents.

Hence, TDA-ICP-MS can be a valuable tool to measure corona thickness of metallic NPs dispersed in solutions with high concentrations of proteins and can be applied not only to Au@SiOx NPs but also to other metal-containing particles or elemental-tagged nanoparticles. However, this method is at present limited to NPs under 70 nm (RH = 35 nm), due to the signal-to-noise ratio of the method, which limits accuracy when size increases. This size limitation allows anyway its use for the vast majority of metal-containing particles designed for nanomedicine's purposes.

Fig. 1 :

 1 Fig. 1: Comparison of size measurements for Au@SiOx NPs at 2 mg.mL -1 in water a) Representative Taylorgram

Fig. 2 : 1 Fig. 3 :

 213 Fig. 2: Influence of BSA concentration on the hydrodynamic radii of Au@SiOx NPs a) Taylorgrams of Au@SiOx NPs incubated in various BSA solutions. Experimental conditions: Au@SiOx NPs, either 2 mg.mL -1 suspensions in pure water or 1 mg.mL -1 suspensions in BSA solutions; injection 40 mbar 3 s; HPC capillary 75 µm i.d. x 70 cm total

Fig. 4 :

 4 Fig. 4: Hydrodynamic radii of Au@SiOx NPs as a function of time in different media. Experimental conditions:

Fig. 5 :

 5 Fig. 5: Mean hydrodynamic radius of Au@SiOx NPs incubated in various incubation media.Experimental conditions: Au@SiOx NPs, diluted at 1mg.mL -1 in various media; injection 40 mbar 3 s; HPC capillary 75 µm i.d. x 70 cm total length; analysis medium: phosphate buffer after 2h incubation (in light grey), after 4h incubation (in dark grey), and incubation medium (in blue); mobilization pressure 50 mbar

  et al. described both reversible and irreversible adsorption of BSA onto silica NPs upon incubation in a BSA solution at 60 µM (4.16 mg.mL -1 ). A thin layer of protein corona (~ 4nm) was still observed by energy-filtered transmission electron microscopy (EF-TEM) around the silica NPs after several washing steps[START_REF] Galdino | Inside the Protein Corona: From Binding Parameters to Unstained Hard and Soft Coronas Visualization[END_REF]. Norde et al. also reported an irreversible adsorption of BSA on silica upon dilution after a 20h-incubation[START_REF] Norde | Adsorption, desorption and re-adsorption of proteins on solid surfaces[END_REF].

Table 1 .

 1 Mean hydrodynamic radii for Au@SiOx NPs as a function of the incubation medium

		BSA		FET
	[BSA] mg.mL -1	Hydrodynamic radius (nm) [FET] (mg.mL -1 )	Hydrodynamic radius (nm)
	0 (phosphate buffer) 30.0 ± 0.7	0 (phosphate buffer) 30.0 ± 0.7
	5	32.8 ± 0.4	0.8	31 ± 1
	12.5	36.8 ± 0.7	5	33 ± 1
	25	35 ± 2	15	35 ± 1
	50	34.2 ± 0.8	30*	62 ± 7

*Hydrodynamic radius was calculated by averaging values at the plateau (time ³ 80 min)

Acknowledgement

This work was supported by the Institute of Chemistry Lyon (Coronano project) and by a French Government Grant managed by the French National Research Agency (ANR-18-CE17-0025-02).

CRediT authorship contribution statement

Arthur Degasperi: Investigation. Lucie Labied: Validation, Software. Carole Farre: Resources, Writingreview & editing. Emmanuel Moreau: Writing -Review & Editing. Carole Chaix: Funding acquisition, Writing -Review & Editing. Matteo Martini: Funding acquisition, Resources, Writing -Original Draft, Review & Editing. Agnès Hagège: Funding acquisition, Project administration, Conceptualization, Methodology, Writing -Original draft, review & editing, Supervision.

Declaration of competing interest

The authors have declared no conflict of interest.


Guillaume Delannoy
email: guillaume.delannoy@ifsttar.fr

Marielle Gueguen Minerbe
email: marielle.gueguen@ifsttar.fri.nourissam.nour@ifsttar.frs.marceausandrine.marceau@ifsttar.frd.diafidinarzed.diafi@ifsttar.frf.farcasfabienne.farcas@ifsttar.fr

Issam Nour

Sandrine Marceau

Philippe Glé
email: p.gléphilippe.gle@cerema.fr

Etienne Gourlay
email: etienne.gourlay@cerema.fr

Dinarzed Diafi

Sofiane Amziane
email: s.amzianesofiane.amziane@uca.fr

Fabienne Farcas

Marielle Guéguen-Minerbe

Impacts of alkalinity and environmental ageing on fungal growth in hemp concretes

Keywords: Hemp concrete, hemp shiv, micro-organisms, moulds, alkalinity

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

INTRODUCTION

Air pollution is an issue studied for decades. The indoor air quality is also more and more monitored and analyzed. Indeed, the way of life, air ventilation, materials used to build or to equip a building, or even cleaning products are the source of indoor pollution, partly due to VOC or to proliferation of microorganisms [START_REF] Koivula | Emissions from thermal insulations-part 2: evaluation of emissions from organic and inorganic insulations[END_REF]. To compare and assess the resistance of material to the developpement of microorganisms, various protocols are used, like visual assessment of mould growth or quantification of viable cells [START_REF] Stefanowski | A rapid screening method to determine susceptibility of bio-based construction and insulation products to mould growth[END_REF][START_REF] Bayon | Development of a laboratory test method to assess the resistance of bio-based insulation materials against moulds[END_REF]. However the results obtained with these protocols are not comparable and do not represent the conditions of use. The proliferation of microorganisms depends of a wide variety of parameters: the chemical composition of the carrier material, the presence of nutrients, the relative humidity in the environment, the water absorption capacity of the material, the amount of spores initially present in the material, the pH, the rugosity or the surface conditions [START_REF] Gradeci | Mould growth criteria and design avoidance approaches in wood-based materials -A systematic review[END_REF].

Biobased materials are increasingly used in buildings, as wood for structures, vegetal wools or biobased concretes for insulation. The higher presence of vegetal biodegradable materials raises questions about the indoor air quality and their durability. Indeed, cellulosic materials, which are a source of nutrients for microorganisms, can have a strong water absorption capacity and also contain a significant amount of spores [START_REF] Hoang | Resistance of green building materials to fungal growth[END_REF]. These viable cells can give rise to a colony in appropriate conditions, and consequently lead to fungal growth, associated with spores release in air and material degradation. Hemp concretes are an example of biobased materials used in buildings. These efficient insulation materials which ensure comfort in habitat are based on vegetal aggregates called hemp shiv, corresponding to the crushed woody part of the stem, covered and bonded by a coating of mineral binder. Therefore, hemp concretes contain vegetal aggregates, rather sensible to fungal growth as mentioned before, and a mineral part, less vulnerable to microorganisms thanks to its high pH value [START_REF] Weyman-Kaczmarkowa | The development of fungi as affected by pH and type of soil, in relation to the occurrence of bacteria and soil fungistatic activity[END_REF]. Fungal growth has been already observed in hemp concrete walls, often after incorrect implementation such as too important water content in formulation, insufficient drying conditions and time or wrong choices of render [START_REF]Isolants biosourcés : points de vigilance[END_REF][START_REF] Piot | Study of a hemp concrete wall exposed to outdoor climate: Effects of the coating[END_REF]. In laboratory, fungal growth has been noticed under environmental conditions considered as extreme: relative humidity higher than 90%, after inoculation of fungi and storage at high humidity and after climatic simulations in presence of vapour and liquid water [START_REF] Marceau | Durability of Bio-based concretes' in Bio-aggregates Based Building Materials[END_REF][START_REF] Arizzi | Predicting the long-term durability of hemp-lime renders in inland and coastal areas using Mediterranean, Tropical and Semi-arid climatic simulations[END_REF]. In other cases, no fungal growth was observed despite wetting and drying cycles with a maximum at 95% RH or successive inoculations of fungi and storage at 80% RH for months [START_REF] Marceau | Influence of accelerated aging on the properties of hemp concretes[END_REF][START_REF] Walker | Mechanical properties and durability of hemp-lime concretes[END_REF]. Besides water activity, the pH is also a significant parameter to consider. Marceau et al. showed that no fungal growth is observed after inoculation of fungi on fresh hemp concrete while a fungal development appears after inoculation on aged specimens [START_REF] Marceau | Influence of accelerated aging on the properties of hemp concretes[END_REF]. Carbonation of binders, particularly of portlandite, induces a decrease of pH from values higher than 12 to pH values lower than 10. Conditions become then more favorable to a fungal growth. The same hypothesis is proposed by Walker et al. [START_REF] Walker | Mechanical properties and durability of hemp-lime concretes[END_REF] to explain the absence of fungal growth despite successive inoculations of fungi.

In buildings, the stress of the environmental conditions is less important than that applied in laboratory experiments quoted above. However, it is important to study the risk of fungal growth on well-manufactured hemp concretes under environmental conditions closer to the real conditions, in order to effectively address the durability issue of hemp concretes. It is the aim of the present study, divided in three parts. Firstly, the quantification and differentiation of viable fungal cells from hemp shiv are performed. Secondly, for fresh hemp concretes, the impact of the alkalinity of the binder on the total amount of fungal cells is investigated. Finally, the fungal growth and the pH are monitored and followed in core and in surface of hemp concretes for 18 months under different ageing conditions.

MATERIALS AND METHODS

Materials

One type of commercial hemp shiv (HS) for insulation and one binder commonly used for hemp concretes are selected. The binder is a commercially formulated binder based on natural lime (70%) and hydraulic and pozzolanic binders (30%), referred as FL in the text.

2.2

Hemp concretes Manufacturing and curing of hemp concretes Hemp concretes for 'wall' application were manufactured as detailed in Delannoy et al. [START_REF] Delannoy | Influence of binder on the multiscale properties of hemp concretes[END_REF]. Hemp concrete is called HC. After manufacturing, 90 days of curing at 65% RH and 20°C ensure the hardening of the hemp concretes by both hydration and carbonation reactions of binders. The manufacturing and curing protocols are detailed in Delannoy et al. [START_REF] Delannoy | Influence of binder on the multiscale properties of hemp concretes[END_REF]. Specimens studied after curing time and before ageing are defined as the initial time A0 (ageing 0 month).

Protocols of ageing

After the manufacturing and the curing time, specimens of hemp concretes have been divided in two groups to be aged under different conditions for 18 months: a static ageing and a dynamic ageing.

The first set is stored in a conditioned room at 50% RH and 20°C. This set serves as reference and this ageing condition is called "Ref".

As the material is sensitive to humidity, an accelerated ageing protocol in climatic chamber has been chosen to study the durability of hemp concretes. The temperature was set up at 30°C, a good temperature for fungal growth. Humidification and drying cycles of 1 week are applied. During the 5 first days, a high relative humidity is applied (98% RH) follow-up by 2 days at 40% RH as a drying step. Transitions between the two RH occur in one hour. These humidification and drying cycles ageing protocol is referred as HD.

At different stages of ageing, specimens are characterized.

Quantification of viable fungal cells

The fungal cells present in 10 g of shiv and in 10 g of fresh hemp concrete (2 hours and 48 hours after the end of the manufacturing) are extracted with a sodium dodecyl sulfate (SDS) water solution (1% weight/volume) for 1 min with the stomacher 400 circulator (320 rpm). Successive 1/10 dilutions are then applied with the SDS solution and 100 µL of the liquid is pipetted onto PDA (Potatoes Dextrose Agar) growth medium in Petri dishes. Dishes are incubated for 2 days at 30°C, after which the colonies are counted, and the numbers of viable fungal cells are presented as colony forming unit per gram of hemp (cfu/ghemp). Two replicated are done. All the procedure takes place in a biological safety cabinet with sterile instruments.

2.4

Analysis of fungal growth To monitor a potential fungal growth on hemp concrete depending on time and on the ageing protocol, one 25 cm² square of adhesive tape is applied for 1 min on the surface or inside the specimens after having been split in two parts. In order to compare with hemp shiv, adhesive tape is also applied on hemp shiv in bulk. Then the adhesive tape is placed and vortexed in 1 mL of the SDS water solution. The same dilution and inoculation protocols presented in section 2.3 are used, with in this case 5 days of incubation. The results are presented in cfu/cm² of material.

2.5

pH analysis Phenolphthalein test is performed to follow the evolution of the pH. This indicator is used to determine if a mineral binder is carbonated or not. The colour indicator has its transition for a pH value around 10: below 10 the material is uncoloured and carbonated, higher than 10 a pink colour appears.

2.6

Microstructure observations and identification of fungi The microstructure of hemp concrete and of fungi developed on Petri dishes are observed on FEI Quanta 400 Scanning Electron Microscope (SEM) with low vacuum and in secondary electron mode (SE).

RESULTS AND DISCUSSION

3.1 Fungi present in hemp shiv and at fresh state in hemp concrete The total amount of viable cells determined with the experiment described in section 2.3 for hemp shiv is 37.8×10 6 cfu/ghemp (Table 1). 6 kinds of colonies are visually differentiable on growth medium. An example of Petri dish is presented Figure 1-a with SEM observations. Penicillium sp. (Figure 1-1) or Aspergillus sp. (Figure 1-4) seems to be present. Figure 1. Picture of a hemp shiv Petri dish and SEM pictures of fungi A high amount of viable cells is found in hemp shiv as described in the last paragraph. In alkali-conditions during the manufacturing of hemp concretes, these fungi can either: grow if they are alkali compatible, wait for more favourable conditions by producing spores or be killed [START_REF] Weyman-Kaczmarkowa | The development of fungi as affected by pH and type of soil, in relation to the occurrence of bacteria and soil fungistatic activity[END_REF]. The study of the impact of alkalinity of binders on the total amount of viable cells helps to determine the origin of fungi observed in buildings or in laboratory studies under high humidity: these fungi can be initially present in shiv and/or come from an external contamination if the initial fungi from the shiv are killed during the manufacturing.

Table 1. Viable fungal cells in hemp shiv and in sample of hemp concrete taken 2h and 48h after manufacturing

Hemp shiv HC 2h after manufacturing HC 48h after manufacturing cfu/ghemp 37.6×10 6 ± 2.3×10 6 4.3×10 6 ± 1.4×10 6 1.5×10 6 ± 0.7×10 6 The total amounts of viable cells for hemp concrete are presented Table 1 for samples taken 2 hours and 48 hours after the end of the manufacturing. A decrease of the cfu/ghemp about 1 log is observed for HC compared to hemp shiv. Even if the values for hemp concrete are getting closer to the limit between potential risk and important risk, fixed at 1×10 6 cfu/g [START_REF] Piot | Study of a hemp concrete wall exposed to outdoor climate: Effects of the coating[END_REF], the amount of viable fungal cells remains important. The conclusion is that the pH of the binder does not mainly kill fungi during the fresh state of the material. Then if a fungal growth appears later on hemp concrete with the appropriate environmental conditions, it could be coming from fungi initially present in hemp shiv as from external contamination.

3.2

Monitoring of the fungi development in and on hemp concretes during ageing As we saw that the alkalinity of binder does not kill fungi cells present in hemp shiv, fungal growth is possible if appropriate conditions are met, such as a good range of temperature, a sufficient water activity or a lower pH on the surface of the material. The amounts of viable fungal cells per cm 2 in core and in surface of specimens of hemp concrete with log scale are presented Figure 2 as function of time and of the type of ageing. First, whatever the type of ageing, at least a 3 log decrease of the amount of fungal cells is observed compared to the amount of cells for hemp shiv in bulk. In the food industry [START_REF] Sheridan | Monitoring CCPs in HACCP systems[END_REF], the amount of viable cells is considered bad for the hemp shiv, and most of the time fair, good or even excellent for hemp concrete. No significant differences can be concluded for all measurements in hemp concrete whatever the time and the type of ageing. Moreover, the fungal cells in hemp concrete do not seem active. Indeed, 5 days of incubation are necessary to allow the counting of colonies on Petri dishes for hemp concretes against 2 days for hemp shiv. Spores are sleeping and need more time to develop on the growth medium. Finally, only 2 different kinds of colonies are visually differentiable compared to the 6 for hemp shiv (Figure 1). The reasons of these results are different according to the ageing protocol applied.

Figure 2. Amount of viable cells on surface and in core of HC with time (at A0 and at 1, 3, 6, 12 and 18 months of ageing) under different ageing conditions (Reference in blue and HD cycles in red) compared to viable cells counted on hemp shiv in bulk (green)

Concerning the static ageing, the explanation to the lowest fungal content compared to hemp shiv is first that the relative humidity applied (50% RH) is lower to the 70% RH necessary to a fungal growth [START_REF] Gradeci | Mould growth criteria and design avoidance approaches in wood-based materials -A systematic review[END_REF], despite a pH lower than 10 in surface observed with the 1 More favourable conditions for fungal growth are brought by the HD ageing, where a high value of relative humidity is fixed during 70% of the time (mean RH per week = 81%). In a first step, this higher relative humidity applied than in Ref ageing impacts the carbonation of FL binder. Indeed, Figure 3.a for HC tested with phenolphthalein after 18 months in HD shows a carbonated depth higher than 2 cm. Even if more volume fraction of hemp concrete has a pH lower than 10, no fungal growth are observed, despite the high value of relative humidity applied. Then, it is the use of alternating periods of high and low relative humidity which could limit the fungal growth [START_REF] Johansson | The effect of cyclic moisture and temperature on mould growth on wood compared to steady state conditions[END_REF][START_REF] Kasanen | Fungal growth and survival in building materials under fluctuating moisture and temperature conditions[END_REF]. In fact, the two days at 40% RH seem to be enough to neutralize the development of hyphae which might start during the 5 days at 98% RH. Moreover when the RH is higher than 95%, a part of the mineral binder is dissolved in the condensate liquid water in hemp concrete, which can migrate into the shiv porosity, and precipitate again during the drying step of HD [START_REF] Delannoy | Influence of binder on the multiscale properties of hemp concretes[END_REF]. This mineralization inside the shiv, observed by SEM on cross-sections of aggregate of HC sample in Figure 3-b (binder in white, cell wall in light grey and pores in dark grey), hinders the accessibility of the lignocellulosic part of hemp, and interferes with the potential fungal growth.

CONCLUSIONS

Hemp concrete is an insulation material mainly made of vegetal aggregates, and consequently considered to be sensitive to fungal growth. In addition to being a source of nutrient thanks to the presence of lignocellulosic molecules, a high amount of viable fungal cells is measured on hemp shiv from the living period of the plant and the retting process.

One formulation of hemp concrete containing a binder based on lime is manufactured for this study. During the manufacturing, the alkalinity of binders does not significantly impact the number of viable fungal cells in the concretes.

After manufacturing, curing time and ageing in static or dynamic conditions, no fungal growth is observed on hemp concrete. Fungal cells remain inactive, despite more favourable conditions brought by the carbonation reaction, leading to a decrease of the pH lower than 10 in the surface. That is explained by a low relative humidity for the static ageing and by the drying step of the dynamic ageing coupled with a mineralization of the porosity of shiv.

In buildings with both outdoor and indoor environments, hemp concrete is submitted to important variations of relative humidity and temperature. Projecting the results obtained in this study to the case of hemp concrete implemented with a render in buildings, with values of relative humidity measured in real wall of hemp concrete (< 70% RH), no fungal growth should be observed without wrong manufacturing or accidental damages.

  on cross-sections of hemp concrete for A0 and after 18 months in Ref ageing (Figure3-a).

Figure 3 .

 3 Figure 3. (a) Pictures of HC tested with phenolphthalein at A0, 18 months after Ref and HD ageing -(b) SEM picture of hemp concrete cross-section


François Dubois

Olivier Lafitte

Analytic solutions and numerical method for a coupled thermo-neutronic problem

Keywords: ordinary differential equation, analytic functions, thermohydraulics, neutronics, Crank-Nicolson numerical scheme

come

Analytic solutions and numerical method for a coupled thermo-neutronic problem

Introduction

In this note, we recall briefly the construction of an analytic solution of the low Mach number thermohydraulic model obtained by Dellacherie et al. detailed in [START_REF] Dellacherie | A simple monodimensional model coupling an enthalpy transport equation and a neutron diffusion equation[END_REF]. Assuming that the neutronics problem, stated in a 1d setting as [START_REF] Abramowitz | Stegun -Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]

d dz (D(h(z)) dφ dz ) + Σ a (h(z))φ(z) = νΣ f (h(z)) k ef f φ(z), z ∈ [0, L], φ(0) = φ(z) = 0, φ ≥ 0
where φ is the neutron density, h is the enthalpy (a measure of the temperature), Σ a and νΣ f are respectively the absorption and fission cross section of the fissile material, D is the diffusion coefficient. This equation is coupled with an equation on the enthalpy: [START_REF] Coddington | Theory of ordinary differential equations -Chapter 8, International series in pure and applied mathematics[END_REF] h (z) = KνΣ f (h)φ(z) supplemented with the boundary conditions h(0) = h e , h(L) = h s . Recall that in [START_REF] Dellacherie | A simple monodimensional model coupling an enthalpy transport equation and a neutron diffusion equation[END_REF], we constructed an analytical solution of (1) and ( 2) for any set of continuous positive functions D, Σ a , νΣ f , as a generalisation of Dellacherie and Lafitte [START_REF] Dellacherie | Une solution explicite monodimensionnelle d'un modèle simplifié de couplage stationnaire thermohydraulique-neutronique[END_REF]. Indeed, introducing the two functions X and Y such that d dh (D(h)νΣ f (h) dX dh ) = Σa(h) νΣ f (h) , -d dh (D(h)νΣ f (h) dY dh ) = 1, X (h e ) = X (h s ) = Y(h e ) = Y(h s ) = 0, and defining Ψ(h, k ef f ) = 2( Y(h) k ef f -X (h)) we have shown that the problem of finding k ef f resumes at using the unique solution z → h(z, k ef f ), increasing on [0, L], such that h(0, k ef f ) = h e , of (h (z, k ef f )) 2 = Ψ(h(z, k ef f ), k ef f ) and at imposing h(L, k ef f ) = h s , which is the equation on k ef f to be solved. It has been proven in [START_REF] Dellacherie | Une solution explicite monodimensionnelle d'un modèle simplifié de couplage stationnaire thermohydraulique-neutronique[END_REF], [START_REF] Dellacherie | A simple monodimensional model coupling an enthalpy transport equation and a neutron diffusion equation[END_REF] that Lemma 1.1. The equation h(L, k ef f ) = h s has a unique solution k * ef f ∈ (0, (max [he,hs] Y X ) -1 ). It is a consequence of k ef f → h(L, k ef f ) decreasing from this interval to (+∞, 0). As, for industrial applications (namely the control of nuclear reactors) the neutronic model has been used first, numerical methods have been developed independently for the two models (neutronics and thermohydraulics). The classical methods used to solve this coupled problem [ref !] is to proceed iteratively (see Annex).

In this contribution, we shall, without losing the generality, simplify the problem (assuming that D and Σ a are constant functions), we shall choose h e = 0, h s = 1, denote by λ = 1 k ef f and consider the ODEs on [0, L] -(Dφ ) + Σ a φ = λΣ(h)φ, h = Kφ supplemented by φ(0) = φ(L) = 0, h(0) = h e , h(L) = h s , where Σ is a given continuous function. Further adimensionnalizing the problem, Σ a = 1, D = 1, L = 1, K = 1, h e = 0, h s = 1. Finally the system studied here is

(3) -d 2 ϕ dz 2 + ϕ(z) = λ Σ h(z) ϕ(z), dh dz = ϕ(z) , 0 < z < 1 , h(0) = 0, h(1) = 1, ϕ(0) = 0, ϕ(1) = 0, λ ≥ 0, ϕ(z) > 0 if 0 < z < 1.

Let ψ λ solution of ψ λ (h) = 2 -2λΣ(h) with ψ λ (0) = ψ λ (1) = 0.

We assume throughout this paper (4) Σ continuous, Σ(h) ≥ Σ * > 0, for 0 ≤ h ≤ 1.

Under this hypothesis, the result of Lemma 1.1 implies the following result Lemma 1.2. System (3) has a unique solution (λ * , h * , φ * ) where h * ∈ C 1 ([0, 1]), φ * ∈ C 2 ([0, 1]).

We intend to present a numerical method which solves the coupled problem without using the numerical methods traditionally used for solving each equation but rather concentrating on solving the equation h(L) = 1. This numerical method, as well as analytic and symbolic methods, are implemented when one knows only three values of Σ, and for simplicity again one assumes that one knows Σ(0) = σ 0 , Σ( 12 ) = σ 1 2 , Σ(1) = σ 1 . The analytic and symbolic methods are consequences of the finding of exact solutions using the incomplete Jacobi functions (which are its solutions) when Ψ is a polynomial of degree 3 or 4.

The following comparisons will be performed.

Analytical approach

We study in this section four representaions of the function Σ which lead to exact analytical solutions of (3) using the incomplete elliptic integrals. Indeed, we consider four cases: Σ constant (σ 0 = σ 1/2 = σ 1 ), Σ an affine polynomial (σ 1/2 = 1 2 (σ 0 + σ 1 )), Σ the interpolation polynomial of degree 2 defined by the three input data σ 0 , σ 1/2 and σ 1 and Σ continuous piecewise affine defined by the three previous input data.

Following e.g. Abramowitz and Stegun [START_REF] Abramowitz | Stegun -Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF], we define the incomplete elliptic integral of the first kind K(m, ϕ) by the relation

K(m, ϕ) = ϕ 0 dθ 1 -m sin 2 θ , 0 ≤ ϕ ≤ π 2 , m < 1.
The complete elliptic integral of the first kind K(m) is defined by

K(m) ≡ K(m, π 2
) . Lemma 2.1. (i) Given two positive reals a and b such that 0 < a < b, we have

(5) a -a dT (a 2 -T 2 ) (b 2 -T 2 ) = 2 b K(m) , m = a 2 b 2 . and (6) b a dT (T 2 -a 2 ) (b 2 -T 2 ) = 1 b K(m) , m = 1 - a 2 b 2 .
(ii) Given a positive real a and a non null real number b, we have

(7) a -a dT (a 2 -T 2 ) (b 2 + T 2 ) = 2 |b| K(m) , m = - a 2 b 2 .
Proof. First cut the first integral into two equal parts, between -a and 0, and between 0 and a. Secondly introduce the change of variable T = a sin θ with 0

≤ θ ≤ π 2 . Then √ a 2 -T 2 = a cos θ, dT = a cos θ dθ, √ b 2 -T 2 = b 1 -a 2 b 2 sin 2 θ and the first relation is established. The same calculus conducts to √ b 2 + T 2 = |b| 1 --a 2
b 2 sin 2 θ and the negative value for the parameter m in the third relation is clear. For the integral (6), we consider the change of variables [START_REF] Abramowitz | Stegun -Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]. We have on one hand T dT = -b 2 m sin θ cos θ dθ and on the other hand

T = b 1 -m sin 2 θ with T (0) = b and T ( π 2 ) = a. Then a = b √ 1 -m and m = 1 -a 2 b 2 ∈ (0,
T 2 -a 2 = b 2 m cos 2 θ. Then b a dT √ (T 2 -a 2 ) (b 2 -T 2 ) = 0 π/2 1 T (-b 2 m sin θ cos θ) dθ b √ m sin θ √ m cos θ = π/2 0 1 b dθ √ 1-m sin 2 θ = 1 b K(m).
In (3), we assume that Σ is known only through the three positive real numbers σ 0 , σ 1/2 and σ 1 which are respectively the values of Σ at 0, 1/2 and 1. In what follows, we establish that for five modelling of the function Σ from these values, it  Franc ¸ois Dubois and olivier Lafitte is possible to put in evidence an analytical approach to determine firstly the scalar parameter λ > 0 and secondly the functions z -→ ϕ(z) and z -→ h(z).

• Zero-th case: Σ is constant.

Then the system modelized by the previous set of equations is totally decoupled and an exact solution can be provided with elementary arguments. Proposition 2.1. If σ 0 = σ 1/2 = σ 1 ≡ µ > 0, the problem (3) admits a unique decoupled solution: we have

λ µ = 1 + π 2 , ϕ(z) = π 2 sin(π z) and h(z) = 1 2 (1 -cos(π z)).
Proof. From the relation -d 2 ϕ dz 2 = (λ µ-1) ϕ(z) and the conditions ϕ(0) = ϕ(1) = 0 with the constraint ϕ(z) > 0 if 0 < z < 1, we deduce that (λ µ -1) is the first eigenvalue of the Laplace equation on the interval (0, 1) wih Dirichlet boundary conditions. Then λ µ -1 = π 2 and ϕ(z) = k sin(π z) for some constant k. We integrate relative to z this relation and we get h

(z) = -k π cos(π z) + for 0 ≤ z ≤ 1. The condition h(0) = 0 imposes = k π and the datum h(1) = 1 implies 2 k π = 1.
If the function Σ is no-more constant, it has been proven in [START_REF] Dellacherie | Numerical results for the coupling of a simple neutronics diffusion model and a simple hydrodynamics low mach number model without coupling codes[END_REF][START_REF] Dellacherie | Une solution explicite monodimensionnelle d'un modèle simplifié de couplage stationnaire thermohydraulique-neutronique[END_REF] that the unknown of the problem can be obtained with the following process. First integrate twice the function Σ and obtain a convex negative function V such that

d 2 V dh 2 = Σ(h), V (0) = V (1) = 0. Second define ψ λ (h) ≡ h (h -1) -2 λ V (h).
Then the equation for the function z -→ h(z) can be written dh dz = ψ λ (h) and the condition h(1) = 1 gives a scalar equation for the unknown λ > 0:

1 0 dh √ ψ λ (h) = 1. A first difficulty is to compute the integral I λ ≡ 1 0 dh √ ψ λ (h)
. Then we can solve easily the equation I λ = 1 with a Newton-like algorithm. Once λ is determined, the explicitation of the functions z -→ ϕ(z) and z -→ h(z) is not difficult. Thus the method we propose is founded on an analytical determination of the integral I λ . We focus our attention to this question in the next sub-sections.

• First case: Σ ∈ P 1 .

In this case, the function Σ is a positive affine function on the interval (0, 1). We set σ 0 = µ (1α), σ 1/2 = µ and σ 1 = µ (1 + α). Then µ > 0 and |α| < 1 to satisfy the constraint of positivity. We introduce the notation ξ ≡ λ µ. Proposition 2.2. With the notations introduced previously, the functions V and ψ λ admit the algebraic expressions

V (h) = µ h (h -1) 1 2 -α 6 + α 3 h , ψ λ (h) = h (h -1) ξ (1 + α 3 ) -1 -2 α ξ 3 h . Then we have I λ = 2 ξ 1+ |α| 3 -1 K(m) with m = 2 |α| ξ 3 ξ+|α| ξ-3 .
• Second case: Σ ∈ P 2 .

In this case, the polynomial ψ λ is a polynomial of degree 4 with real coefficients and ψ λ is positive on the interval (0, 1). We have also ψ λ (0) = ψ λ (1) = 0. Recall that ψ λ (h) = h (h-1)-2 λ V (h) with V ≡ Σ the Lagrange interpolate polynomial such that Σ(0) = σ 0 , Σ( 12 ) = σ 1/2 and Σ(1) = σ 1 . All these coefficients are  supposed positive: σ 0 > 0, σ 1/2 > 0, σ 1 > 0. We introduce appropriate parameters µ, α and δ such that σ 0 = µ (1α), σ 1/2 = µ (1 + δ) and σ 1 = µ (1 + α). In this sub-section, we exclude the case of a linear interpolation, id est δ = 0, in coherence with the hypothesis that the degree of the polynomial ψ λ is exactly equal to 4.

Lemma 2.2. With the above notations and properties, the new parameters µ, α and δ are defined by µ = 1 2 (σ 0 + σ 1 ), α = 1 -σ0 µ and δ = -σ0-2 σ 1/2 +σ1 σ0+σ1

. They satisty the inequalities µ > 0, |α| < 1 and δ > -1. If we set γ ≡ 1 -1 3 |α| + 2 3 δ, we have the inequality γ > 0.

Proof. The two first inequalities for µ and α are clear because σ 0 > 0 and σ 1 > 0. From σ 1/2 > 0 and µ > 0 we deduce that δ > -1. Then γ > 1 -1 3 + 2 3 = 0. Lemma 2.3. With the above notations and properties, introduce the notation ξ ≡ λ µ; then ξ > 0. Introduce also the two other roots p and g of the polynomial ψ λ . We can write ψ λ (h) = a 0 h (h -1) (hp) (hg). Then a 0 = 2 3 ξ δ. Denote by σ ≡ p + g the sum and π ≡ p g the product and the two roots p and g. We have a 0 σ = 2 3 ξ (δ + α) and a 0 π = 1ξ + 1 3 α ξ -2 3 δ ξ. Proof. The inequality ξ > 0 is a consequence of the hypothesis λ > 0 and of the property µ > 0 established in Lemma (2.2). The other results are obtained with elementary algebra. We have used the Sagemath [START_REF] Stein | SAGE: System for Algebra and Geometry Experimentation[END_REF] package, an open source mathematical software system. Lemma 2.4. With the above notations and properties, if the integral I ≡

1 0 dh √ ψ λ (h)
is a positive real number, the roots p and g cannot be equal to 0 or 1.

Proof. If one of the two roots p or g is equal to 0 or 1, the polynomial ψ λ (h) has a double root and the function (0, 1)

h -→ 1 √ ψ λ (h)
is not integrable in the interval (0, 1). Proposition 2.3. We keep active the above notations.

(i) If the discrete positive family (σ 0 , σ 1/2 , σ 1 ) is the trace of a concave function, id est if σ 1/2σ 0 > σ 1σ 1/2 , we have δ > 0, a 0 > 0 and γ ξ > 1. The two roots of the function ψ λ are real with opposite signs. To fix the ideas, p < 0 < 1 < g. If we set

(8) ∆ ≡ σ 2 -4 π = 1 a 2 0 ξ (5 δ 2 + α 2 + 6 δ) ξ -6 δ , we have p = 1 2 (σ - √ ∆) and g = 1 2 (σ + √ ∆) with σ and π introduced in Lemma (2.3).
(ii) Conversely, if the function ψ λ are real zeros with opposite signs, then

σ 1/2 -σ 0 > σ 1 -σ 1/2 . Proof. (i) Recall that a 0 > 0 because δ > 0 in this case. Set f (h) = h 2 -σ h + π ≡ (h -p) (h -g). Then ψ λ (h) = -a 0 h (1 -h) f (h) and the function f must be negative on (0, 1) because ψ λ is positive on this interval. But a 0 f (0) = a 0 π = 1 -ξ 1 -1 3 α + 2 3 δ , a 0 f (1) = 1 -ξ 1 + 1 3 α + 2 3 δ . Then ξ inf 1 -1 3 α + 2 3 δ, 1 + 1 3 α + 2 3 δ > 1 and in other words γ ξ > 1.
The polynomial f is strictly negative on the interval (0, 1) and the inequality p < 0 < 1 < g is just a notation that distinguish p as the "petite" root and g as the "grande" root. 

(ii) Conversely, if the polynomial f is strictly negative on the the interval (0, 1), we must have a > 0 because the function ψ λ (h) is positive on this interval. Then δ > 0 and

σ 1/2 -σ 0 > σ 1 -σ 1/2 .
Proposition 2.4. We suppose as in Proposition (2.3) that the function ψ λ (h) has two real zeros p and g that satisfy p < 0 < 1 < g. Then the integral I ≡

1 0 dh √ ψ λ (h)
can be computed with the help of the following formulas:

A = -g (g-p) 1-p , B = (1-p) (g-p) g , a = A+1 A-1 , b = B+1 B-1 , m = 1 -a 2 b 2 I λ = 1 √ a0 1 -a b 1 √ g-p K(m) . Moreover, we have 0 < a < 1 < b and 0 < m < 1.
Proof. The idea is given in the book [START_REF] Abramowitz | Stegun -Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]. We introduce an homographic mappping

T (h) ≡ h-d h-c
and positive numbers a and b in order to transform the four roots of ψ λ (h) into the family of numbers {±a, ±b}. More precisely, we enforce the following conditions T (p) = b, T (0) = -b, T (1) = -a and T (g) = a, as presented on Figure 1. If h satisfies ψ λ h = 0, define T by the relation

T = h-d h-c . Then T ∈ {±a, ±b}. After three lines of elementary algebra, we have h-h = (h-c) T -T T -1 with T = h-d h-c . Then ψ λ (h) = a 0 (u -c) 4 (T 2 -a 2 ) (b 2 -T 2 ) (1-a 2 ) (b 2 -1)
. If we enforce the condition 0 < a < 1 < b, we still have ψ λ (h) > 0 for 0 < h < 1 and the transformed integral has to be computed on the interval (-b, -a):

I = √ (1-a 2 ) (b 2 -1) √ a0 (d-c) -a -b dT √ (T 2 -a 2 ) (b 2 -T 2 ) = √ (1-a 2 ) (b 2 -1) √ a0 (d-c) b a dT √ (T 2 -a 2 ) (b 2 -T 2 ) = √ (1-a 2 ) (b 2 -1) √ a0 b (d-c) K(m)
with m = 1 -a 2 b 2 due to the relation [START_REF] Dellacherie | Une solution explicite monodimensionnelle d'un modèle simplifié de couplage stationnaire thermohydraulique-neutronique[END_REF]. We have now to determine the four parameters a, b, c and d as a function of the data p and g. We have the constitutive relations

c -d = (b -1) (p -c) = (-b -1) (0 -c) = (-a -1)(1 -c) = (a -1) (g -c) .
From the third equality, we obtain c = -a+1 b-a < 0. We report this value inside the second and fourth equalities to obtain a system of two equations for the parameters a and b: p

(b-a)+a+1 = -(a+1) (b+1) b-1 and g (b-a)+a+1 = -(a+1) (b+1) a-1 . Taking the difference, we have -a+1 a-1 b+1 b-1 = g-p > 0. Then A ≡ a+1 a-1 is negative, B ≡ b+1 b-1
is positive and -A B = gp. We insert the transformed parameters A and B into the equation g (ba)

+ a + 1 = -(a+1) (b+1) a-1
and after some lines of algebra, we

obtain -A (1 -p) = g B. Since A < 0 and B > 0, it is clear that A = -g (g-p)
1-p and B = b-a > 0. We can achieve the evaluation of the integral:

I = √ (1-a 2 ) (b 2 -1) √ a0 b (d-c) K(m) = 1 √ a0 √ (1-a 2 ) (b 2 -1) (a+1) (b+1) b-a b K(m) = 1 √ a0 1 -a b 1 √ -A B K(m) = 1 √ a0 1 -a b 1 √ g-p K(m)
and the result is established.  ; case p < 0 < 1 < g.

If the sequence σ 0 , σ 1/2 , σ 1 is convex, id est if σ 1/2 -σ 0 < σ 1 -σ 1/2
, we have δ < 0, a 0 < 0 and the discriminant ∆ introduced in (8) can be negative or positive. We begin by the case ∆ < 0 and we have two conjugate roots for the polynomial ψ λ . To fix the ideas, we set p K(m) with m = 1 ∆ < 0, (ii) if σ < 1, the integral I can be computed with the following relations

= 1 2 (σ + i √ -∆) and g = 1 2 (σ -i √ -∆) with σ introduced in Lemma 2.3.
       ∆ = 1 16 (σ 2 -∆) (σ -2) 2 -∆ > 0 , c = -1 1-σ ∆ + 1 4 (σ 2 -∆) < 0 , a = 1 1-2 c ∈ (0, 1) , b = √ -∆ 2 c-σ , d = -a c > c , m = -a 2 b 2 < 0 , I = 2 √ (1-a 2 ) (b 2 +1) √ |a0| |b| (d-c) K(m) , (iii) if σ > 1,
the integral I is evaluated with the same relations that in the case (ii), except that we have now c > 0, a < 0 and d < c.

Proof. First recall that ψ λ (h) = |a 0 | h (h -1) (h -p) (h -g) = (-a 0 ) h (1 -h) (h -σ 2 ) 2 -∆ 4 . (i) If σ = 1, we make the change of variable: v = h -1 2 . Then I = 1 √ |a0| 1/2 -1/2 dv 1 4 -v 2 v 2 -∆ 4
and we can use the relation [START_REF] Dellacherie | A simple monodimensional model coupling an enthalpy transport equation and a neutron diffusion equation[END_REF]

with a = 1 2 and b = √ -∆ 2 . Then I = 1 √ |a0| 4 √ -∆ K(m) with m = -1 (-∆) = 1 ∆ < 0. (ii) If σ < 1,

we make an homographic change of variables T (h) = h-d

h-c as in Proposition 2.4. We enforce now the conditions

T (0) = -a, T (1) = a and T σ±i √ -∆ 2 = ∓i b  with real coefficients a, b, c and d. If h is a root of ψ λ , id est h ∈ {0, 1, p, g}, we have as previously h -h = (h -c) T -T T -1 with T ∈ {±a, i ± b}. Thus ψ λ (h) = |a 0 | (h -c) 4 (a 2 -T 2 ) (T 2 +b 2 ) (1-a 2 ) (1+b 2 )
and we must enforce a 2 < 1. We have also dT =

d-c (h-c) 2 dh. Then I = 1 √ |a0| a -a (h-c) 2 dT d-c √ (1-a 2 ) (1+b 2 ) (h-c) 2 √ (a 2 -T 2 ) (T 2 +b 2 ) = √ (1-a 2 ) (1+b 2 ) √ |a0| (d-c) a -a dT (a 2 -T 2 ) (T 2 +b 2 ) = 2 √ (1-a 2 ) (1+b 2 ) √ |a0| |b| (d-c) K(m)
with m = -a 2 b 2 due to the relation [START_REF] Dellacherie | A simple monodimensional model coupling an enthalpy transport equation and a neutron diffusion equation[END_REF]. The relation is established, but we have now to construct the homography h -→ T (h), taking into consideration the constraints d > c and a 2 < 1. The constitutive relations T (0) = -a, T (1) = a and

T σ±i √ -∆ 2 = ∓i b take the form c -d = (-a -1) (0 -c) = (a -1) (1 -c) = = (-i b -1) ( σ+i √ -∆ 2 -c) = (i b -1) ( σ-i √ -∆ 2 -c) . Then i b+1 i b-1 = i √ -∆+(2 c-σ) i √ -∆-(2 c-σ) and b = √ -∆ 2 c-σ . We have also a+1 a-1 = 1 c -1 then a = 1 1-2 c
. We can write now the equation relative to the parameter c: c

(a + 1) = (i b+1) (c-σ+i √ -∆ 2 ). But (i b+1) (c-σ+i √ -∆ 2 ) = (2 c-σ) 2 -∆ 2 (2 c-σ)
and a+1 = 2 1-c 1-2 c . An equation of degree 2 for the parameter c is emerging:

f (c) ≡ (1 -σ) c + 2 σ 2 -∆ 4 c -σ 2 -∆ 4 = 0 . We remark that the constant term σ 2 -∆ 4
is the square of the modulus of the two roots p and g. So we have two roots c -and c + for the equation f (c) = 0. We remark that f ( σ 2 ) = -1 4 ∆ (σ -1) < 0 in this case and σ 2 is between the two roots. We have also f ( 12 ) = 1 4 (1σ) > 0 and 1 2 is outside the two roots. We deduce the inequalities c -< 0, c -< σ 2 < c + < 1 2 , c + > 0. Before making a choice between the two numbers c ± , observe that we have cd = (a + 1) c and d = -a c. Moreover a = 1 1-2 c and the simplest choice is a > 0. Then the condition a 2 < 1 take the form a -1 < 0, but a -1 = 2 c 1-2 c and we must have c < 0. Then d > 0 and c < 0

< d ; the condition d -c > 0 is satisfied. The reduced discriminant ∆ is equal to σ 2 -∆ 4 2 + (1 -σ) σ 2 -∆ 4 = 1 16 (σ 2 -∆) (σ -2) 2 -∆ and is positive. Then c = c -= -1 1-σ ∆ + 1 4 (σ 2 -∆) . (iii) If σ > 1,
the beginning of the previous proof is unchanged. But now the second order equation for the parameter c takes the form

f (c) ≡ (σ -1) c -2 σ 2 -∆ 4 c + σ 2 -∆ 4 = 0 . We still have ∆ = 1 16 (σ 2 -∆) (σ -2)
2 -∆ > 0 and the two roots c ± have the same sign. We have now f ( σ 2 ) = -1 4 ∆ (1σ) < 0 and σ 2 is still between the two roots. Moreover, f ( 12 ) = 1 4 (σ -1) > 0 and 1 2 is outside the two roots. We deduce the inequalities 1 2 < c -< σ 2 < c + . Then whatever the choice between c ± , 2 c -1 > 0 and a < 0. The condition a 2 < 1 imposes now a + 1 > 0 and a + 1 = 2 1-c 1-2 c is positive if and only of 1c < 0. We observe also that

f (1) = -1 -σ 2 2 + ∆ 4 < 0, then c -< 1 < c + . The choice c + is mandatory and we have c = c + = 1 σ-1 ∆ + 1 4 (σ 2 -∆) .



An analytic and symbolic analysis of a coupled problem ; case 0 < 1 < p < g.

Proposition 2.6. With the notations recalled previously, in the case ∆ > 0 with two real roots with the same sign, the computation of the integral I ≡

1 0 dh √ ψ λ (h) is splitted into two sub-cases: (i) if σ > 2, the integral I can be evaluated according to        ∆ = 4 p g (p -1) (g -1) , b = 1 (g-p) (g+p-1) p (p -1) + g (g -1) + √ ∆ > 1 , a = b (p -g) + p + g -1 , c = a+1 2 a , d = a c < c , m = a 2 b 2 , I = 2 √ (1-a 2 ) (b 2 -1) √ |a0| b (c-d) K(m) ,
(ii) if σ < 0, the integral I is evaluated with the relations

       ∆ = 4 p g (p -1) (g -1) , b = 1 (g-p) (1-g-p) p (p -1) + g (g -1) + √ ∆ > 0 , a = b (p -g) + 1 -p -g , c = -1-a 2 a < 0 , d = -a c > 0 > c , m = a 2 b 2 , I = 2 √ (1-a 2 ) (b 2 -1) √ |a0| |b| (d-c) K(m) ,
that are very similar to the first case. In both cases, we have 0 < a < 1 < b.

Proof. (i) In the case σ > 2, we define an homographic function

T (h) = h-d h-c
represented in Figure 2 in order to enforce the inequalities 0 < a < 1 < b and 0 < d < 1 < p < c < g. The constitutive relations beween the roots take now the form T (0) = a, T (1) = -a, T (p) = -b and T (g) = b . We can write, with the same notations as previously,

h -h = (h -c) T -T
T -1 and after mutiplying the four analogous factors,

ψ λ (h) = |a 0 | (u -c) 4 (a 2 -T 2 ) (b 2 -T 2 ) (1-a 2 ) (b 2 -1) positive for T 2 < a 2 < b 2 . Then I = -a a (h-c) 2 √ |a0| (d-c) √ (1-a 2 ) (b 2 -1) (h-c) 2 √ (a 2 -T 2 ) (b 2 -T 2 ) dT = √ (1-a 2 ) (b 2 -1) √ |a0| |b| (c-d) a -a dT √ (a 2 -T 2 ) (b 2 -T 2 ) = 2 √ (1-a 2 ) (b 2 -1) √ |a0| b (c-d) K(m)
with m = a 2 b 2 due to the relation [START_REF] Dellacherie | Une solution explicite monodimensionnelle d'un modèle simplifié de couplage stationnaire thermohydraulique-neutronique[END_REF]. The explicitaion of the coefficients a to d is conducted as previously: We first have  (ii) In the case σ < 0, we follow the same method than previously. The homographic function T (h) = h-d h-c represented in Figure 3 and the parameters satisfy the conditions 0 < a < 1 < b and p < c < g < 0 < d < 1. We impose the following permutation between the roots:

0 < c -d = (a -1) (0 -c) = (-a -1) (1 -c) = (-b -1) (p -c) (b -1) (g -c
T (p) = b, T (g) = -b, T (0) = -a and T (1) = a . Then we have as previously

ψ λ (h) = |a 0 | (u -c) 4 (a 2 -T 2 ) (b 2 -T 2 ) (1-a 2 ) (b 2 -1) , positive for T 2 < a 2 < b 2 . Then we have I = a -a (h-c) 2 √ |a0| (d-c) √ (1-a 2 ) (b 2 -1) (h-c) 2 √ (a 2 -T 2 ) (b 2 -T 2 ) dT = 2 √ (1-a 2 ) (b 2 -1) √ |a0| |b| (d-c) K(m) with m =
≡ (g -p) (1 -g -p) b 2 -2 (p 2 + g 2 -p -g) b + (g -p) (1 -g -p) .
The reduced discriminant ∆ is still given by the expression ∆ = 4 p g (p -1) (g -1) > 0. The two roots are positive, we have f (1) = 4 g (1-g) > 0 and the only root greater that 1 is b = 1 (g-p) (1-p-g) p (p -1) + g (g -1) + √ ∆ . The proof of the proposition is completed.

• Third case: Σ is a continuous positive fonction, affine in each interval (0, 1 2 ) and ( 12 , 1). The function Σ is defined by its values Σ(0) = σ 0 , Σ( 1 2 ) = σ 1/2 and Σ(1) = σ 1 . We introduce new parameters, still denoted by α and β to represent the data:

σ 0 = σ 1/2 (1 -α) and σ 1 = σ 1/2 (1 + β) .
Then the inequalities α < 1 and β > -1 express the constraints σ 0 > 0 and σ 1 > 0. Moreover, σ 1/2 is positive. Lemma 2.5. If the continuous function Σ defined on (0, 1) by its values ; case p < g < 0 < 1.

Σ(0) = σ 0 , Σ( 1 2 ) = σ 1/2 , Σ(1) = σ 1 is affine in each interval (0,
d 2 ψ λ dh 2 = 1 -2 λ Σ and ψ λ (0) = ψ λ (1) = 0 admits the following expression ψ λ (h) = ψ 0 (h) if h ≤ 1 2 ψ 1 (h) if h ≥ 1 2 , with ψ 0 (h) = -2 3 α ξ h 3 + (α ξ -ξ + 1) h 2 + β-5 α 12 ξ + ξ -1 h ψ 1 (h) = 2 3 β ξ (1 -h) 3 + (-β ξ -ξ + 1) (1 -h) 2 + 5 β-α 12 ξ + ξ -1 (1 -h) , with ξ ≡ λ σ 1/2 .
We have in particular ψ λ ( 12 ) = β-α 24 ξ + 1 4 (ξ -1). Moreover, if ψ λ is positive on the interval (0, 1),

dψ λ dh (0) = β-5 α 12 ξ + ξ -1 > 0 , dψ λ dh (1) = -5 β-α 12 ξ + ξ -1 < 0. Remark 1.
We observe that the expression of ψ 1 is obtained from the expression of ψ 0 by making the transformations α ←→ (-β) and h ←→ (1h).

Proof. The function Σ admits the algebraic expression:

Σ(h) = (2 α h -α + 1) σ 1/2 if h ≤ 1 2 (2 β h -β + 1) σ 1/2 if h ≥ 1 2 .
We integrate two times, enforce the conditions ψ λ (0) = ψ λ (1) = 0 and impose the continuity of ψ λ and dψ λ dh at the specific value h = 1 2 . The result follows.

We have to compute the integral I ≡

1 0 dh √ ψ λ (h)
. We have the following calculus:

I = 1/2 0 dh √ ψ 0 (h) + 1 1/2 dh √ ψ 1 (h) = 1/2 0 dh √ ψ 0 (h) + 1/2 0 dh √ ψ 1 (1-h)
.

Due to the Remark 1, the determination of the second term relative to ψ 1 is very analogous to the term associated to ψ 0 . In the following, we will concentrate essentially to the evaluation of the first integral

I 0 ≡ 1/2 0 dh √ ψ 0 (h) . Proposition 2.7. If α = 0, the integral I 0 ≡ 1/2 0 dh √ ψ 0 (h)
admits several expressions parameterized by ξ ≡ λ σ 1/2 :

(i) if ξ = 1, then β > 0 and I 0 = 2 6 β ,  O h ψ 0 O h ψ 0 O h ψ 0 Figure 4
. Contribution ψ 0 to the function ψ λ . The case α > 0 is on the left and the case α < 0 is splitted into two figures: two complex roots in the middle and three real roots on the right.

1 h ψ 1 1 h ψ 1 1 h ψ 1
Figure 5. Contribution ψ 1 to the function ψ λ . The case β > 0 is splitted into two sub-cases: two complex roots on the left and three real roots in the middle; the case β < 0 is presented on the right.

(ii) if ξ > 1, I 0 = 2 √ ξ-1 arctan 1 ϕ0 with ϕ 0 = 1 + β ξ 6 (ξ-1) and (iii) if ξ < 1, I 0 = 1 c √ b0 log √ c 2 +2+c √ c 2 +2-c with b 0 = β 12 ξ + ξ -1 > 0 and c = 1-ξ b0 .
Proof. We have in this particular case ψ 0 (h) = -(ξ -1)

h 2 + β 12 ξ + ξ -1 h. (i) If ξ = 1, ψ 0 (h) = β 12 h and I 0 = 12 β 1/2 0 dh √ h = 24 β . (ii) If ξ > 1, we write ψ 0 (h) = (ξ -1) h 1 + β ξ 12 (ξ-1) -h . With a ≡ 1 + β ξ 12 (ξ-1) > 0, we make the classical change of variable h (a -u) = ϕ 2 h 2 . Then ϕ 2 = a h -1 and 2 ϕ dϕ = -a h 2 dh. If h = 1 2 , then ϕ = ϕ 0 > 0 with ϕ 2 0 = 2 a -1 = 1 + β ξ 6 (ξ-1) . In consequence, I 0 = ϕ0 +∞ -2 h 2 ϕ dϕ a √ ξ-1 1 ϕ h = 2 √ ξ-1 +∞ ϕ0 dϕ 1+ϕ 2 = 2 √ ξ-1 π 2 -arctan(ϕ 0 ) = 2 √ ξ-1 arctan 1 ϕ0 because ϕ 0 > 0. (iii) If ξ < 1, we can write ψ 0 (h) = b 0 (h + c 2 h 2 ) with b 0 = β 12 ξ + ξ -1 > 0 and c = 1-ξ b0 . We make now the change of variable h + c 2 h 2 = θ 2 h 2 . Then θ 2 = c 2 + 1 h > c 2 and 2 θ dθ = -1 h 2 dh. Then I 0 = 1 √ b0 √ c 2 +2 +∞ (-2 θ h 2 dθ) 1 θ h = 2 √ b0 +∞ √ c 2 +2 h dθ = 2 √ b0 +∞ √ c 2 +2 dθ θ 2 -c 2 = 1 c √ b0 log θ-c θ+c +∞ √ c 2 +2 = 1 c √ b0 log √ c 2 +2+c √ c 2 +2-c .
The reader will observe that if ξ tends to 1, each of the results proposed in (ii) and (iii) converge towards the expression proposed in (i).

Proposition 2.8. If α > 0, recall that ψ 0 (h) = -2 3 α ξ h 3 + (α ξ -ξ + 1) h 2 + β-5 α
12 ξ + ξ -1 h. It corresponds to the left part of Figure 4. The function ψ 0  admits three real roots and we set ψ 0 (h) ≡ a 0 (h-r -) h (r + -h) with a 0 = 2 3 α ξ > 0 and r -< 0 < r + . Then

I 0 ≡ 1/2 0 dh √ ψ 0 (h) = 2 √ a0 (-r-) K(ϕ 0 , m) with m = r+ r-< 0 and ϕ 0 = arcsin 1 √ 2 r+ .
Proof. The explicitation of the algebraic expressions of r -and r + has no interest and is not detailed here. To compute I 0 , we make the change of variables h = r + sin 2 θ. Then dh = 2 r + sin θ cos θ. Then

ψ 0 (h) = a 0 (r + sin 2 θ -r -) r + sin 2 θ r + cos 2 θ = a 0 r 2 + sin 2 θ cos 2 θ (-r -) 1 -r+ r-sin 2 θ = a 0 (-r -) r 2 + sin 2 θ cos 2 θ (1 -m sin 2 θ) with m = r+ r-< 0. Moreover, if h = 1 2 , then θ = ϕ 0 with sin 2 ϕ 0 = 1 2 r+ . Then I 0 = ϕ0 0 2 dθ √ a0 (-r-) 1 √ 1-m sin 2 θ
and the result is established.

Proposition 2.9. The case α < 0, corresponds to the middle and right pictures of Figure 4. Recall that ψ 0 (h) = -2 3 α ξ h 3 + (α ξξ + 1) h 2 + β-5 α 12 ξ + ξ -1 h. (i) If the polynomial ψ 0 (h) has two complex roots, we write it under the form ψ 0 (h) ≡ a 0 h (h 2 + µ h + ζ 4 ) with a 0 = 2 3 |α| ξ > 0 and we have the inequality

µ 2 -4 ζ 4 < 0; we introduce m = 1 2 -µ 4 ζ 2 that satisfies 0 < m < 1 and ϕ 0 = 2 arctan 1 ζ √ 2 . Then we have I 0 = 1 ζ √ a0 K(ϕ 0 , m). (ii) If the polynomial ψ 0 (h) has three real roots, we set ψ 0 (h) ≡ a 0 h (r -- h) (r + -h) with a 0 = 2
3 |α| ξ > 0 and 0 < r -< r + ; we consider m = r- r+ ∈ (0, 1) and ϕ 0 = arcsin

1 √ 2 r-
. The integral I 0 is computed with the following expression:

I 0 = 2 √ a0 r+ K(ϕ 0 , m).
Proof. This proof is directly inspired by the book [START_REF] Abramowitz | Stegun -Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]. Nevertheless, we give the details herein for a complete explanation of the final relations.

(i) We operate the change of variable h = ζ 2 t 2 with t = tan θ 2 . Then

ψ 0 (h) = a 0 ζ 2 t 2 ζ 4 (1 + t 4 ) + µ ζ 2 t 2 = a 0 ζ 2 t 2 ζ 4 (1 + t 4 ) + µ ζ 2 t 2 = a 0 ζ 6 t 2 (1 + t 2 ) 2 -2 t 2 + µ ζ 2 t 2 = a 0 ζ 6 t 2 (1 + t 2 ) 2 1 -1 2 -µ 4 ζ 2 4 t 2 (1+t 2 ) 2 . Because the discriminant µ 2 -4 ζ 4 < 0 is negative, we have |µ| 2 ζ 2 < 1 and -1 < -µ 2 ζ 2 < 1. In other terms, -1 2 < -µ 4 ζ 2 < 1 2 and 0 < m ≡ 1 2 -µ 4 ζ 2 < 1. We remark that sin θ = 2 t 1+t 2 and we get ψ 0 (h) = a 0 ζ 6 t 2 (1 + t 2 ) 2 (1 -m sin 2 θ). We have also dh = 2 ζ 2 t (1 + t 2 ) dθ 2 and dh = ζ 2 t (1 + t 2 ) dθ. With this change of variables, the upper bound is equal to ϕ 0 such that ζ 2 tan 2 ϕ0 2 = 1 2 , id est ϕ 0 = 2 arctan 1 ζ √ 2 .
We can now achieve the calculus of the integral:

I 0 = 1/2 0 dh √ ψ 0 (h) = ϕ0 0 ζ 2 t (1+t 2 ) dθ √ a0 ζ 3 t (1+t 2 ) √ 1-m sin 2 θ = 1 ζ √ a0 ϕ0 0 dθ √ 1-m sin 2 θ = 1 ζ √ a0 K(ϕ 0 , m).
(ii) If the polynomial ψ 0 (h) ≡ a 0 h (r -h) (r +h) has three real roots, we set h = r -sin 2 θ and ψ 0 (h) = a 0 r -sin 2 θ r -cos 2 θ (r +r -sin 2 θ) = a 0 r 2 -sin 2 θ cos 2 θ 1 -r- t+ sin 2 θ .  With 0 < m ≡ r- t+ < 1, we have ψ 0 (h) = a 0 r 2 -sin 2 θ cos 2 θ (1m sin 2 θ). The upper bound ϕ 0 of the integral is associated to the relation 1 2 = r -sin 2 ϕ 0 and ϕ 0 = arcsin

1 √ 2 r- . We have finally , I 0 = 1/2 0 dh √ ψ 0 (h) = ϕ0 0 1 √ a0 r+ 2 dθ √ 1-m sin 2 θ
and the result is established.

Remark 2. As emphazized previously, the comutation of the second term

I 1 ≡ 1 1/2 dh √ ψ 1 (h) = 1/2 0 dh √ ψ 1 (1-h)
of the global integral I = I 0 + I 1 is obtained from the change of parameters α ←→ (-β) in Propositions 2.7, 2.8 and 2.9. The different cases are illustrated on Figure 5.

Semi-analytical approximation

In this section, we revisit the results of Dellacherie et al. [START_REF] Dellacherie | Numerical results for the coupling of a simple neutronics diffusion model and a simple hydrodynamics low mach number model without coupling codes[END_REF]. In this paper, they solved the equation V = Σ, V (0) = V (1) = 0 using a finite element approximation in the two-dimensional space {V a,b , V a,b (h) = h(h -1) (a h + b)}, Σ being the interpolation polynomial of degree 2 associated with σ 0 , σ 1/2 and σ 1 . We consider again this specific case to be able to compare with the results of section 2 and we add the case where Σ is piecewise affine with the same approach.

• First semi-analytical case: Σ is a P2 polynomial approximated by an affine function More precisely, we consider the discrete vector space W generated by the functions h -→ h (1h) and h -→ h 2 (1h). Then we project the function

V (h) = σ 0 (1 -h) (1 -2 * h) + 4 σ 1/2 h (1 -h) + 2 σ 1 h * (2 * h -1)
onto this space in the sense of least squares, as described in [START_REF] Dellacherie | Une solution explicite monodimensionnelle d'un modèle simplifié de couplage stationnaire thermohydraulique-neutronique[END_REF]. We obtain after the resolution of a 2 × 2 linear system a new function V of degree 3 that correcponds to an affine function. We obtain Σ = 1 5 (3 σ 0 + 4 σ 1/2 -2 σ 1 ) (1h) + (-2 σ 0 + 4 σ 1/2 + 3 σ 1 ) h . Then the process follows analogously to the first case with an affine function.

• First semi-analytical case: Σ is a continuous positive fonction, affine in each interval (0, 1 2 ) and ( 1 2 , 1), approximated by an affine function This case is analogous to the previous one, except that the initial function Σ is piecewise affine:

Σ(h) = σ 0 (1 -2 h) + 2 σ 1/2 h if h ≤ 1 2 2 σ 1/2 (h -1) + σ 1 (2 h -1) if h ≥ 1 2 .
We obtain after some elementary calculus Σ = 1 16 (11 σ 0 + 10 σ 1/2 -5 σ 1 ) (1h) + (-5 σ 0 + 10 σ 1/2 + 11 σ 1 ) h . And the end of the process is analogous to the previous case.

Numerical method

In this section, we describe a general numerical method (valid for any continuous funtion Σ) which solves directly the coupled problem (3) and finds an equation for computing k * ef f introduced in lemma 1.1. Note that this numerical method, unlike the one described in this Annex, does not use any coupling of codes. It needs only to sole one equation with one single real unknown.



Recall that problem (3) is

-d 2 ϕ dz 2 + ϕ(z) = λ Σ h(z) ϕ(z), dh dz = ϕ(z) , 0 < z < 1 , h(0) = 0, h(1) = 1, ϕ(0) = 0, ϕ(1) = 0, λ ≥ 0, ϕ ≥ 0.
We introduce a nonregular meshing 0 = z 0 < z 1 < ... < z N -1 < z N = 1 of the interval [0, 1] and we set ∆z j+1/2 ≡ z j+1z j for j = 1, ... , N . We integrate the differential equation dh dz = ϕ(z) with the Crank-Nicolson scheme: h j+1h j ∆z j+1/2 = 1 2 ϕ j + ϕ j+1 for j = 0, ... , N -1 .

Then after two integrations, the first equation can be written as ϕ j = ψ λ (h j ) and the Crank-Nicolson scheme takes the form

h j+1 -h j 1 2 ψ λ (h j ) + ψ λ (h j+1 ) = ∆z j+1/2 (λ) .
Note that the choice of the Crank-Nicolson algorithm allows to recover a discretization of each sub-problem (nemely the idealized neutronic one and the simplified thermo-hydraulic one). We impose the values h j ≡ sin 2 π j 2 N in order to take into account the singularities and two boundary conditions of the problem at z = 0 and z = 1. The notation ∆z j+1/2 (λ) in the right hand side of the previous relation is justified by the fact that if the numbers h j are given, the left hand side is a simple function of the scalar λ. The number λ is a priori not known, but we have the natural relation 9) has a unique solution λ ∈ (λ * , +∞) when N -1 j=0

N -1 j=0 ∆z j+1/2 (λ) = 1 that takes the form (9) N -1 j=0 h j+1 -h j 1 2 ψ λ (h j ) + ψ λ (h j+1 ) = 1 . Lemma 4.1. Introduce V the unique solution of V = Σ, V (0) = V (1) = 0. We have observed previously that ψ λ (h) = h (h -1) -2 λ V (h). Define λ * = max [0,1] h(h-1) V (h) . (i) For λ > λ * , ψ λ (h j ) > 0 (ii) Equation (
hj+1-hj

1 2 √ ψ λ * (hj )+( √ ψ λ * (hj+1) ∈ R + > 1. Proof. (i) One notes that V < 0 on [0, 1] thanks to Σ ≥ Σ * . Hence ψ λ (h) = (-V (h))(λ -h(1-h) -V (h) ) > 0 on [0, 1]. (ii)
The limit of ψ λ (h i ) is +∞ for at least one value of i when λ → +∞, the function ψ λ (h i ) is increasing for any i, strictly increasing for i = 0, N . The inequality

N -1 j=0
hj+1-hj

1 2 √ ψ λ * (hj )+( √ ψ λ * (hj+1)
> 1 (where h j or h j+1 can be a point of maximum of h(h-1) V (h) in which case the sum is +∞) and λ -→

N -1 j=0
hj+1-hj

1 2 √ ψ λ * (hj )+ √ ψ λ * (hj+1)
strictly decreasing yields the result.

For a fixed discretization with N meshes, a Newton algorithm can be implemented without difficulty. With this procedure, we recover on one hand an approximated value λ N of the unknown λ and on the other hand the entire approximate solution of the problem h j ≈ h(x j ) and ϕ j ≈ ψ λ (h j ). Observe that at convergence of the Newton algorithm, the abscissas x j are a function of the solution λ N  Franc ¸ois Dubois and olivier Lafitte and the converged space mesh is a result of the problem. This coupled problem can be reduced to a single equation with only one real variable even after discretization!

Numerical results

They are presented in Figures 6 and7. They correspond to the decreasing data σ 0 = 8, σ 1/2 = 6 and σ 1 = 3. We obtain the following exact values for the scalar parameter: λ = 1.89036 in the decoupled case (case 0), λ = 1.99533 in the affine case (case 1), λ = 1.86593 in the parabolic case (case 2), λ = 1.89454 in the piecewise affine case (case 3) and λ = 1.85769 in the parabolic case approached by an affine functions (semi-analytical case 1) λ = 1.88614 in the piecewise affine case projected on affine functions (semi-analytical case 2). For each of these six cases, our numerical approach gives converging results at second order accuracy for the parameter λ.

If one wants to relate these calculations with the usual problemes solved in the neutronics world, by homogeneity, it is enough to consider the values σ 0 = 14.92744, σ 1/2 = 11.19558 and σ 1 = 5.59779. For this case, the following exact values for k * ef f are 0.98708 in the decoupled case (case 0), 0.93515 in the affine case (case 1), 1.00000 in the parabolic case (case 2), 0.98490 in the piecewise affine case (case 3), 1.00444 in the parabolic case approached by an affine functions (semi-analytical case 1) and 0.98928 in the piecewise affine case projected on affine functions (semi-analytical case 2). This investigation shows that some times, the exigence of accuracy of the operational calculations [START_REF] Tiberga | Results from a multi-physics numerical benchmark for codes dedicated to molten salt fast reactors[END_REF] could be lightened.

Conclusion

We considered in this paper a simplified idealized one-dimensional coupled model for neutronics and thermo-hydraulics. A numerical method allowing to find the unique solution of this coupled model (uniqueness of the multiplication factor k ef f 

We write L 2 ([0, 1]) = Rψ h ⊕ H and we consider the decomposition of S( h)P -1 0 S( h) on this decomposition of spaces. It writes, for any h ∈ C 0 ([0, 1]) A 11 ( h) A 12 ( h) A 21 ( h) A 22 ( h) .

Observe that A 12 (h) = 0, A 21 (h) = 0. One has L h -L h = S( h)P -1 0 S( h) -S(h)P -1 0 S(h) = S( h)P -1 0 (S( h) -S(h)) + (S( h) -S(h))P -1 0 S(h), which implies

||L h -L h || L(L 2 ,L 2 ) ≤ 2Σ 1 2 ∞ C 0 N |h -h| C 0 .
For h -h small enough in C 0 , the smallest eigenvalue of S( h)P -1 0 ( h), that is µ( h), satisfies the identity [START_REF] Stein | SAGE: System for Algebra and Geometry Experimentation[END_REF] A 11 ( h)µ( h) -A 12 ( h)(A 22 ( h)µ( h)I) -1 A 21 ( h) = 0.

Recall that, thanks to the inequality on L h -L h , one has

|A 11 ( h) -A 11 (h)| ≤ 2Σ 1 2 ∞ C 0 N |h -h| C 0 , ||A 12 ( h)|| L(L 2 ,L 2 ) ≤ 2Σ 1 2 ∞ C 0 N |h -h| C 0 , ||A 21 ( h)|| L(L 2 ,L 2 ) ≤ 2Σ 1 2
∞ C 0 N |h -h| C 0 , and that the smallest eigenvalue of A 22 ( h) is µ 2 ( h), hence A 22 ( h)-µI ≥ (μ 2 (h)-µI) (on L 2 ([0, 1])). Hence, from (10)

|A 11 ( h) -µ( h)| ≤ (2Σ 1 2 ∞ C 0 N |h -h| C 0 ) 2 ||(A 22 ( h) -µ( h)I) -1 || L(L 2 ,L 2 ) ≤ 4Σ∞(C0N ) 2
µ2( h)-µ( h) |h -h| 2 C 0 . There exists 0 > 0 such that for |h -h| C 0 ≤ 0 , µ 2 ( h)µ( h) ≥ µ2(h)-µ(h)

( 1 -

 1 p) (g-p) g . Then a = A+1 A-1 and b = B+1 B-1 . The first constitutive equation gives the value d = -c b. We observe that dc = (a+1) (b+1)

Figure 1 .

 1 Figure 1. Homographic transformation for the computation of the integral 1 0 dh √ ψ λ (h)

Proposition 2 . 5 .

 25 With the notations recalled previously, in the case ∆ < 0 with two conjugate roots, the computation of the integral I ≡1 0 dh √ ψ λ (h)is splitted into three sub-cases:(i) if σ = 1, we have I =

Figure 2 .

 2 Figure 2. Homographic transformation for the computation of the integral 1 0 dh √ ψ λ (h)

1 (

 1 ) and we get after some algebra conducted with the help of formal calculus, d = a c, c = a+1 2 a , a = b (pg) + p + g -1 and the parameter b is a root of the second degree polynomialf (b) ≡ (gp) (g + p -1) b 2 -2 (p 2 + g 2 pg) b + (gp) (g + p -1) whose leading coefficient (gp) (g + p -1) is positive. The reduced associated discriminant ∆ = 4 p g (p -1) (g -1)is positive and we have just to compare the roots with the value 1 because we want to enforce b > 1. But f (1) = 4 p (1p) < 0 and 1 is between the two roots. So we have b = g-p) (g+p-1) p (p -1) + g (g -1) + √ ∆ > 1 and the first part of the proof is completed.

a 2 b 2

 2 thanks to the relation (5). The explicitation of the coefficients a, b, c and d begin with the algebraic form 0 > cd = (-a -1) (0c) = (a -1) (1c) = (b -1) (pc) (-b -1) (gc) of the conditions T (p) = b, T (g) = -b, T (0) = -a and T (1) = a. We deduce that d = -a c, c = -1-a 2 a , a = b (pg) + 1pg and the coefficient b is solution of the equation of second degree f (b)

1 2 )Figure 3 .

 23 Figure 3. Homographic transformation for the computation of the integral 1 0 dh √ ψ λ (h)

Figure 6 .Figure 7 .

 67 Figure 6. function Σ on the left, V on the right.

2 ,0

 2 hence the required estimate for A 11 ( h)µ( h), which implies the estimate forµ( h)µ(h) = µ( h) -A 11 ( h) + A 11 ( h) -A 11 (h).An eigenvector associated to µ( h) writes ψ h + r, whereA 11 ( h)µ( h) A 12 ( h) A 21 ( h) A 22 ( h)µ( h)I ψ h r = 0, hence one gets r = -(A 22 ( h)µ( h)I) -1 A 21 ( h)ψ h .One has thusψ h = k(ψ h + r), where k -1 = 1 0 Σ( h)(ψ h (z) + r(z))dz. One has ψ hψ h = (k -1)ψ h + kr, and k -1 = k(1 -1 k ) = k( |h -h| C 0 , 

This contribution has been presented at LAMA Chambery and CEA Saclay in spring

and of the neutron flux profile) is based on the Crank-Nicolson scheme. Simultaneously, special function (incomplete elliptic integrals) are used for finding analytic solutions for five different methods of representation of the fission cross section known at three points. Both methods match perfectly. We observe important differences in the multiplication factor, even if the neutron flux is really similar. Future work concerns, for example, increasing the number of discretization points of cross sections.

Annex: Coupling algorithm

We derived in this paper an approach based on the use of an equation which takes into account the multiphysics of the problem, indeed constructing our method for solving numerically the problem using a totally coupled equation.

However, the method traditionally used for solving such a problem is based on the coupling of existing codes solving each equation.

The aim of this Annex is to describe the best result that can be obtained so far for such a coupling algorithm, and to give a sufficient condition for the convergence of this coupling algorithm. It is a follow-up of a first version (unpublished, on hal, see [START_REF] Dellacherie | Une solution explicite monodimensionnelle d'un modèle simplifié de couplage stationnaire thermohydraulique-neutronique[END_REF]) of this numerical coupling problem. Observe that the numerical method described in Section 4 converges with very weak hypotheses on Σ, while the coupling algorithm uses a very sharp condition on the coefficients and on the initial point chosen (see Lemma 6.3) The equation coming from the thermohydraulics model is here h (z) = (h sh e )φ(z), where 1 0 φ(z)dz = 1 (in the reduced model described in the paper h s -h e = 1 and one chooses to use this value from now on). This traditional numerical procedure, which couples codes rather than models reads as follows:

(1) Σ = Σ(0) ⇒ λ 0 , φ 0 , with

Observe that the unique solution (h * , φ * , λ * ) obtained in Lemma 1.2 for the system (3) satisfies that

When treating only the differential equation of the system (3), one has, for the equation for φ, λ, assuming that Σ(h) is a C 1 function bounded below by Σ 0 > 0 Lemma 6.1.

(1) The set of values of λ such that (-

Franc ¸ois Dubois and olivier Lafitte

(2) For h ∈ C 0 , the operator L h is self-adjoint compact on L 2 ([0, 1]), hence its spectrum is a countable sequence of eigenvalues µ j (h), decreasing to 0 when j → +∞.

One denotes by µ(h) the largest eigenvalue of L h . One denotes ψ h its associated unique eigenvector such that

We shall use the following regularity result, for (h, h) two C 1 functions satisfying the end conditions h(0) = h(0) = 0, h(1) = h(1) = 1. Assume M 0 > 0 be given and denote by

). Proposition 6.1. There exists 0 > 0 such that for all h, h in B such that |h -h| C 0 ≤ 0 :

• Control of the approximation of the difference of the largest eigenvalue for h and for h µ(h)µ( h) at order 2:

• control of the difference of the normalized eigenvectors

Let K be the constant

Lemma 6.2. The sequence h n , λ n described in the numerical procedure satisfies,

We deduce immediately Lemma 6.3. This algorithm converges under the sufficient condition

It is, evidently, a very strong condition, which is fulfilled when C 0 is small and when ||ψ n || L 2 is bounded uniformly in n for n large enough (meaning also that one is in a neighborhood of the solution).

Proof of Proposition 6.1. Consider the operator L h . The operator S(h

For C0

The inequality

follows. The proposition is proven. Note

In a similar way, (h n+1h n )(z) =

Under the hypothesis of Lemma 6.2 on h n+1h n , we deduce the estimate on | 1 λn+1 -1 λn |.