Dr B Mehdaoui

Dr A Meffre

J Carrey

Dr S Lachaize

Dr L M Lacroix

Dr M Gougeon

Prof B Chaudret

Prof M Respaud

Dr J Carrey Corresponding Author

Optimal size of nanoparticles for magnetic hyperthermia: A combined theoretical and experimental study

Keywords: hyperthermia, magnetic nanoparticles, iron, hysteresis loop calculation 2

consequence of theories deriving from the Stoner-Wohlfarth model; they also appear clearly in numerical simulations. These results open the path to a more accurate description, prediction and analysis of magnetic hyperthermia.

Progresses in the prediction and optimization of the heating of magnetic nanoparticles in an alternative magnetic field are highly desirable for their application in magnetic hyperthermia.

Here a model system consisting of metallic iron nanoparticles with a size ranging from 5.5 to 28 nm is extensively studied. Different regimes as a function of the nanoparticles size are evidenced: single-domain superparamagnetic, single-domain ferromagnetic and multi-domain.

Ferromagnetic single-domain nanoparticles are the best candidates and display the highest specific losses reported in the literature so far (11.2±1 mJ g -1 ). Measurements are analysed using state-of-the-art analytical formula and numerical simulations of hysteresis loops.

Several features expected theoretically are observed for the first time experimentally: i) the correlation between the nanoparticle diameter and their coercive field ii) the correlation between the amplitude of the coercive field and the losses iii) the variation of the optimal size with the amplitude the magnetic field. None of these features are predicted by the linear response theory -generally used to interpret hyperthermia experiments-but are a natural

Introduction

Magnetic hyperthermia (MH) is a promising therapeutic method based on the use of magnetic nanoparticles (MNPs). MH is a two step procedure: MNPs are first injected into the tumour; the patient is then immersed in an alternating magnetic field with a frequency f and an amplitude µ0Hmax appropriately chosen. Due to the excitation of the MNPs, the temperature of the tumour rises. This temperature increase improves the efficiency of chemotherapy (T~42-45°C) or can even directly kill the tumor cells by necrosis (T>50°C) [1,2,3,4,5,6].

When applying an alternative magnetic field to MNPs, the area of their magnetic hysteresis loop A corresponds to a dissipated energy, which is used in MH. The power generated by the MNPs is evaluated by this specific area A (specific losses) or, more commonly used, by their specific absorption rate (SAR), the two parameters being linked by the equation SAR = Af. Exploiting this effect in oncology requires biocompatible MNPs with large heating power, thus one of the challenges in this topic is to maximize SAR values. The SAR being an increasing function of f and µ0Hmax, one could think of increasing these parameters. However, human body cannot be exposed to alternating magnetic field of large µ0Hmaxf product, due to medical limitation [7]. Therefore, optimizing the MNPs requires optimizing their SAR at finite and specific values of the applied frequency and magnetic field.

Typical values used in medical treatments so far are 100 kHz and 20 mT [8].

For experimentalists, there are three main ways to optimize the SAR. The first one is to increase the saturation magnetization MS of MNPs since the maximum SAR achievable is directly proportional to MS. The second way is to optimize the anisotropy of the MNPs since the importance of this parameter has been recently shown [9]. These first two points help in determining the optimal materials. Finally, for a given material, the next major parameter is the MNP size. An experimental study of this parameter requires synthesizing monodisperse MNPs with a varying mean size and studying the size-dependence of SAR. This has been realised in a few articles centred on the properties of iron oxide MNPs [10,11,12,13,14]. 3

In the previous studies, when a quantitative analysis was tempted, [11,12,13,14] the linear response theory (LRT) was always used, which is a major problem. Indeed, this theory has a restricted domain of validity [9], but the studies mentioned above have at some point used LRT out of this domain, which casts doubt on the results obtained (on this issue, see [ 15 ]). For a rigorous analysis of experimental data, theories derived from the Stoner-Wohlfarth model in their domain of validity should be used in complement to LRT [9,16].

Out of the domain of validity of these two models, only numerical simulations can confidently calculate the SAR of single-domain MNPs [9,17]. However, such simulations have never been applied so far for a direct comparison with experiments.

The aim of the present article is precisely to compare experiments performed on a model system consisting of Fe MNPs, with a tuneable mean size ranging from 5.5 to 28 nm, with state-of-the-art methods of analysis consisting of numerical simulations and adequate theories. Iron is a material with a weak anisotropy and a high MS; its magnetic properties are less sensitive to unwanted size and surface effects such as dead layers than its oxide counterparts due to its metallic character. The synthesized Fe MNPs studied here are highly crystalline and display a MS value close to the bulk one in the whole range of size studied. Thus, these MNPs constitute a perfect system for a fine study on the influence of the MNP size on MH.

Here we show that the magnetic-field dependence of the SAR and its amplitude strongly depends on the MNP size. In particular the specific losses are maximized for MNPs in a narrow range of diameters, and reach unprecedented values. Behaviour typical of singledomain MNPs in the superparamagnetic regime, in the ferromagnetic regime and of multidomain MNPs are observed. Moreover, several features in the size-dependence of MH are reported here for the first time and can be well understood using theories derived from the Stoner-Wohlfarth model and numerical simulations. An original method for the quantitative analysis of MH is proposed and applied to our experimental results, permitting their fine analysis. Thus, this article confirms the limitation of LRT theory to explain experimental results on systems of interest for MH. It evidences that the use of suitable theories opens the path to more accurate description, prediction and analysis of MH. 4

Synthesis and characterisations

Our group has been investigating the elaboration of pure metallic iron(0) nanoparticles using an organometallic approach for about 10 years. [18,19,20] This approach is based on the decomposition of an iron source in a reductive medium in the presence of surfactants. By varying the reactants (iron source, surfactants) and the experimental conditions (concentration, temperature or pressure), we covered a large panel of size and shape for iron(0) nanoparticles. In this article, the iron sources were either the iron dimer {Fe[N(SiMe3)2]2}2 or the ultra-small Fe(0) NPs whose preparation has already been reported elsewhere. [20] The surfactant system can be one of the two following combinations: palmitic acid (PA) / hexadecylamine (HDA) or hexadecylamonium chloride (HDA⋅HCl) / HDA. In all the cases, the reactants were reacted in mesitylen, under a reductive atmosphere of dihydrogen at 150°C for 2 days. Depending on the conditions, we forced different growth mechanisms to take place in order to change the morphology of the final NPs. Thus, the spherical NPs of samples 1 and 2 were formed thanks to an isotropic growth regulated by the PA/HDA ratio.

Increasing the quantity of carboxylic acid over the amine led to a size increase because the growth is then favoured over the nucleation. [19] The same trend was observed for samples 3 to 7 that were prepared with the ammonium chloride HDA⋅HCl at various concentrations. [21] Samples 3 and 5 were prepared from the {Fe[N(SiMe3)2]2}2 iron source whereas samples 4, 6 and 7 were formed from ultra-small Fe(0) NPs. Changing the iron source nature helped us to produce different samples by playing with the different reduction kinetics. Tuning the experimental conditions finally allowed us to produce 7 samples of different shapes and mean sizes displaying a magnetization close to the bulk value. One main difference between the spherical and the cubic NPs is their aggregation in solution. Due to the shape and the larger mean size, the samples with cubic NPs are not colloids: the NPs self-organize into large superlattices. The latter are strongly stabilized by Van der Waals interactions between the alkyl chains of the surfactants at the surface of the NPs, and by magnetic interaction between the ferromagnetic NPs. These two contributions prevent any good dispersion of the NPs in compatible solvents.

Transmission electron microscopy (TEM) characterisations of the samples are shown in Figure 1. Magnetic parameters extracted from SQUID measurements (see Supplementary Information) are summarized in Table 1. For all samples, except sample 4, the saturation magnetization per unit mass σS is close to the bulk Fe magnetization (212 A m 2 kg -1 ).

Magnetic hyperthermia measurements

MH measurements at 54 kHz were performed on a specially designed frequencyadjustable electromagnet [22]. This setup has been recently improved and allows us to work with a magnetic field up to µ0Hmax = 60 mT for frequencies f in the range 5-100 kHz.

Measurements at 274 kHz were performed using the coil of an induction oven, with a magnetic field up to µ0Hmax = 73 mT, similarly to experiments described in Ref. [23] (see Experimental Section). Preliminary results on the MH properties of samples 1 and 4 were published in Refs [24] and [23], respectively.

The SAR dependence as a function of magnetic field at frequencies f = 54 kHz and f = 274 kHz for samples 1-7 are shown in Figures 2a and2b respectively. For sample 1 (d = 5.6 nm) the magnetic field dependence of the SAR follows a power law with an exponent of 2.8 and 3.1 at frequencies f = 54 kHz and 274 kHz respectively. Their behaviour has been analysed in details in a previous article and is typical of MNPs in the superparamagnetic regime [24].

In samples 3-7 (d = 8.9-27.5 nm), the SAR displays at 274 kHz a saturation at large magnetic field (see Figure 2b), which is typical of the ferromagnetic regime [16]. At 54 kHz, the saturation is only visible on sample 3 (d = 8.9 nm) because of the smaller maximum magnetic field available on this setup. On these samples the SAR shows a small amplitude at low magnetic field and a sharp increase above the coercive field of the MNPs, which is another feature of the ferromagnetic regime. The evolution of the coercive field with size is more clearly evidenced in Figure 2c and2d, where the normalized SAR values are shown for these samples. In samples 3-6, the coercive field increases with the size of the MNPs. It can be observed that the amplitude of the SAR at large magnetic field is directly related to the amplitude of the coercive field : the larger the coercive field, the larger the SAR. In sample 7 (d = 27.5 nm), both the coercive field and the SAR displays an abrupt drop (see Figure 2a-d).

Another interesting way to plot the same experimental data is to display the size dependence of SAR for various applied magnetic fields. Such a plot is shown for the two frequencies in Figures 2e and2f. It can be observed that, for each magnetic field, an optimal 6 diameter dopt maximizes the amplitude of SAR. Interestingly, dopt is not constant but increases with the amplitude of the applied magnetic field.

Finally, the degree of optimization α achieved in a given MNPs system for an application in MH can be extracted from the SAR value using [9,23]:

f H SAR S σ αµ max 0 4 = ( 1
)
where σS is taken here as the bulk Fe value. α is a dimensionless parameter which characterizes the relative area of the hysteresis loops with respect to the ideal square (α = 1).

In Figures 2g and 2h the evolution of α as a function of the magnetic field and MNP diameter is shown and its value at saturation is summarized in Table I.

For a possible application of such nanoparticles in medical applications, the µ0Hmaxf product is limited to be safe for the patient. We have performed hyperthermia measurements on the various samples at 35 kHz and 57 mT so as to keep the same µ0Hmaxf product as the one currently used at the Charité Hospital [8]. The evolution of SAR, A and α with size are shown as Supporting Information. Sample 6 (d = 19.7 nm) displays the largest SAR = 290 W/g, corresponding to A = 8.3 mJ/g and α = 0.17.

A few comments on the data shown in Figure 2 can be made. First, the SAR values measured on the optimized samples are very large. At 54 kHz, the maximum SAR = 427±42 W g -1 occurs for sample 6 (d = 19.7 nm), while at f = 274 kHz, the maximum SAR = 3066±306 W g -1 occurs for sample 5 (d = 13.7 nm). In the latter case, this corresponds to specific losses A = 11.2±1 mJ g -1 , which are the highest reported in the literature so far.

Interestingly, large values of α are obtained only for values of the magnetic field greater than 50 mT. In a previous article, we had theoretically shown that iron MNPs -due to the value of their anisotropy field around 50 mT-should in principle be optimized nanoparticles only when magnetic fields of this order of magnitude are applied [9]. This is experimentally confirmed here. The largest α = 0.19 is observed in sample 5 (d = 13.7nm) under µ0Hmax = 43 mT and f = 274 kHz. Though being the highest among our experiments, this value is still far 7 from the maximum value achievable in a system with randomly oriented anisotropy axis α = 0.39 [9]. Thus our system could be optimized further, as will be discussed in Section 3.2.2.

Analysis

Qualitative analysis

Magnetic field dependence of SAR

In this part, MH properties are calculated numerically without attempt to fit the experimental data. The results are then compared to experimental results for a qualitative analysis. The hysteresis loop area A is calculated for an assembly of spherical, uniaxial MNPs with their anisotropy axis randomly distributed in space, as detailed in Ref. [9]. The magnetic parameters of bulk iron have been used (MS = 1700 kA/m, the effective anisotropy Keff = 4.8×10 4 J m -3 ). The other parameters are T = 300 K, f = 100 kHz, and the frequency factor of the Néel-Brown relaxation time τ0 = 5×10 -11 s As an illustration, the calculated hysteresis loops for d = 16 nm are shown in Figure 3a as a function of the maximum field applied.

In Figure 3b the evolution of A with magnetic field for MNPs with a diameter d ranging from 4 nm to 24 nm is shown. In Figure 3c normalized data illustrating more clearly the evolution of the curve shape for the different diameters are shown. It can be seen that the shape of the A(µ0Hmax) evolves significantly with the diameter. For the smallest MNPs (d = 4 nm), A shows a square dependence with the magnetic field in agreement with LRT for superparamagnetic MNPs [9,25]. For sizes up to d =10 nm, A follows a power law, whose exponent vary with MNPs size [9]. For d ≥ 14 nm, the shape of the curves is characteristic of the ferromagnetic regime: the losses are very weak at low magnetic field followed by an increase around the coercive field and finally a progressive saturation [16,26].

The comparison between numerical calculations (Figure 3) and experiments (Figure 2)

shows a qualitative agreement on three main features: i) the difference of behaviour between the superparamagnetic regime and the ferromagnetic one, the former displaying small SAR values and no saturation at large magnetic field, ii) the correlation between the amplitude of the coercive field and the amplitude of SAR for MNPs in the ferromagnetic regime and iii) 8

the increase of SAR with the MNP diameter at high magnetic field. The main point of discrepancy comes from the sample 7 (d = 27.5 nm), which displays an abrupt drop of coercive field and SAR. This drop in coercive field can be explained by the transition toward a multi-domain state of the MNP magnetization. Indeed, a transition from single-domain to vortex configuration has been calculated by micromagnetic simulations for Fe particles of roughly d = 26 nm [27]. Moreover, magnetic holography measurements performed on a 30 nm nanocube, similar to the ones composing this sample, have evidenced a vortex state of the magnetic configuration [28].

Optimal size

To study the optimal diameter dopt or the optimal volume Vopt for MH, the evolution of A as a function of the diameter at a fixed µ0Hmax value must be evaluated. In order to do so, hysteresis loops are calculated numerically for various sizes, the parameter values being chosen as previously. Figure 4a shows an example for µ0Hmax = 20 mT -which is the magnetic field used during treatments [8]. For very small nanoparticles (d < dopt), the hysteresis loop is almost reversible and specific losses A almost null, MNPs are superparamagnetic. When the volume increases, the hysteresis loop progressively opens and displays a typical ferromagnetic shape. The area of the hysteresis loop reaches a maximum for d = dopt = 17.4 nm, plotted in thick line. For d > dopt, the coercive field is too large compared to the applied magnetic field, so the magnetization of the MNPs cannot be switched, leading to minor loops only. Thus, the losses decrease with increasing size, and vanish at very large size.

The evolution of the area A as a function of the MNP size is summarized in Figure 4b for various magnetic field amplitudes. The presence of an optimal diameter is clearly visible, permitting to plot the evolution of dopt as a function of the magnetic field (Figure 4c). At very low magnetic field amplitude, the optimal diameter dopt =11.7 nm is independent of the magnetic field amplitude. With strengthened µ0Hmax, dopt increases. The optimal size, derived here above from numerical calculations, can also be predicted analytically using two models.

At low field, for µ0Hmax < 2.7 mT, LRT is valid, so the optimum volume is given by the following equation [9]:

) ln( 0 τ πf K T k V B opt = (2) 9
The optimal volume is thus constant at low field, as represented in Figure 4c by the horizontal solid line.

In the ferromagnetic regime, LRT is no more valid, therefore theories derived from the Stoner-Wohlfarth model should be used instead. In this framework, one can have access to the optimal volume through the optimization of the coercive field µ0HC which reads [9]:

), 1 ( 48 . 0 8 . 0 0 0 κ - = K C H µ H µ ( 3
)
where

S K M K H 2 =
is the anisotropy field and . 4 ln

0 max 0         = τ κ Vf M H µ T k KV T k S B B (4)
Note that Equation ( 3) is only valid when κ < 0.7. The maximum value for A is obtained for an optimal coercive field, the value of which is slightly below the applied magnetic field [6] :

( ) . 04 . 0 81 . 0 max 0 0 H µ H µ C ± = (5)
Using Equations ( 3) and ( 5), one can then deduce the optimal volume and thus, the optimal diameter. The result of this resolution is shown in Figure 4c as a dashed line in the domain of validity of Equation ( 3). This figure illustrates that the evolution of the optimal size at large magnetic field is a natural consequence of theories derived from the Stoner-Wohlfarth model.

Here again, comparing Figure 4b with Figure 2e-f shows that this feature is observed experimentally, since the optimal diameter increases with the magnetic field. This evolution comes with an increase of SAR. This is the first time that this feature is observed experimentally.

Quantitative analysis

Method and hypothesis

The experimental and theoretical results presented above are in good qualitative agreement, insofar as all the main features predicted theoretically are observed experimentally. However, there are some quantitative disagreements on the values of the losses, on the optimal diameters and on the coercive field values. These differences cannot be explained only by the frequency used. For instance, a coercive field of ca. 10 mT is 10 determined numerically for 14nm MNPs at 100 kHz (cf Figure 3b) while values of ca. 40 and 35 mT are respectively obtained at 54 and 274 kHz (cf Figure 2c-d). This discrepancy might result from several effects: different values of τ0 and/or Keff, influence of size distribution and/or of magnetic interactions. In order to determine the origin of this discrepancy, a quantitative analysis of the experimental results is required.

Due to the fact that the calculation time of a typical graph in Figures 3 or 4 is that it can be calculated analytically for a sample composed of single-domain MNPs in the ferromagnetic regime. ii)

Thus our analysis is restricted to the samples 3-6 (d = 8.9-19.7 nm), given that these samples match these requirements. iii) We assume that these four samples have the same Keff value, which it is not unreasonable since surface effects are limited in this size range, and have also the same τ0 value. iv) The evolution of CHyp H µ 0 with diameter is fitted analytically to determine τ0 and Keff. v) These two parameters are then injected into numerical simulations to obtain the complete and exact SAR values.

The calculation of

CHyp H µ 0 as a function of the experimental parameters requires solving numerically:

                          - = 8 . 0 0 0 0 0 4 ln 1 463 . 0 τ Vf M H µ T k V K T k H µ H µ S CHyp B eff B K CHyp (6)
Details on the origin and domain of validity of this equation can be found in Supporting Information.

Results

The quantitative analysis of the data has been performed for measurements at 55 and at 274 kHz and lead to compatible results. We will focus here on the results obtained at 274 11 kHz, provided that the values of µ0HCHyp have a lower error bar due to the higher range of magnetic field used. In Figure 5a, the coercive field µ0HCHyp deduced from the highest slope of the SAR(µ0Hmax) function is displayed for the four samples. The value for the sample 7 (d = 27.5 nm) is also shown for indication but is not taken into account in the quantitative analysis due to its multi-domain configuration. The size evolution of the coercive field is fitted using Equation ( 6) with τ0 and Keff as varying parameters and σS = 212 A m 2 kg -1 , T = 300 K and f = 274 kHz. The best fit of experimental data is obtained for Keff = 9.1×10 4 J m -3 and τ0 = 5×10 -11 s (see Figure 5a). These values have a correct order of magnitude, though τ0 is smaller than the values usually reported for individual MNPs (10 -9 -10 -10 s) and Keff slightly larger than the bulk value. These values were used for numerical simulations and the results compared with experimental data. In Figure 5b, the numerically calculated SAR obtained for the four samples using these parameters are shown. The coercive fields extracted from these curves are shown in Figure 5a

Conclusions

Extensive measurements have been performed on a model system consisting of iron MNPs in a size range going from 5.5 to 27.5 nm. The results have been analysed using state-of-the-art numerical simulations of hysteresis loops and analytical formula. Several features theoretically predicted have been observed for the first time in MH experiments: i) the increase of coercive field with the size of MNPs ii) the correlation between the amplitude of the coercive field and the amplitude of the SAR iii) the fact that the optimal size depends on the amplitude of the magnetic field. None of these features are predicted by the LRT, though this is the theory that has been the most used so far in MH. These features are however a natural consequence of theories deriving from the Stoner-Wohlfarth model and also appear clearly in numerical simulations. Furthermore the transition toward multi-domain MNPs is clearly detected experimentally by an abrupt drop of both coercive field and SAR, confirming that such MNPs are not the most efficient for MH, the optimized objects being single-domain MNPs.

The variation of coercive field with the diameter of the MNPs is quantitatively reproduced by a simple model of monodisperse non-interacting single-domain MNPs with size-independent τ0 and Keff values. However, such a model does not correctly reproduce the amplitude of the SAR and its magnetic field dependence. The quantitative analysis we have performed shows that this discrepancy is not due to the effect of size distribution but more likely to the effect of magnetic interactions.

Iron MNPs with optimized properties display hysteresis area values A = 11±1 mJ g -1 at µ0Hmax = 73 mT, which is the highest value of the literature so far. The study of their α parameter has shown that they reveal their potential only when large magnetic fields are applied. Keeping a µ0Hmaxf product similar to the one used during MH treatments, they 13 display a SAR of 290 W g -1 . We have seen that this value could be further enhanced since α and thus the SAR could be multiplied by a factor between 3 and 7. The next step of our study will be to reduce the importance of magnetic interactions in our system to check if MH experiments and numerical simulations would in quantitative agreement in system of noninteracting MNPs.

Once protected from oxidation, these high-magnetization nano-objects could be useful for several biomedical applications besides magnetic hyperthermia, such as magnetic labelling, drug targeting or contrast enhancement in MRI [30]. For this purpose, a silica layer can be used, as recently patented by our group [31]; this layer could also ensure an easy functionalization of the objects. However, the question of their eventual toxicity should unquestionably be addressed.

Experimental

Nanoparticles synthesis : All the preparations and syntheses are performed under either an inert atmosphere of argon or a reductive atmosphere of hydrogen using a glove-box and Fisher-Porter bottles. Mesitylene (99%, VWR Prolabo) was dried over alumina and degassed by three pump-freeze-thaw cycles. Hexadecylamine (HAD, 99%, Aldrich) was used without purification. Hexadecylammonium chloride was prepared according to published procedure. [32] The {Fe[N(SiMe3)2]2}2 dimer was brought to NanoMePS and used as-received. Ultrasmall iron(0) NPs were prepared according to published procedure [20].

Synthesis of samples 1 and 2:

A solution of palmitic acid (256 mg, 1 mmol for sample 1 and 358 mg, 1,4 mmols for sample 2; 384 mg, 1.5 mmols) in mesitylene (5 ml) was added to a green solution of {Fe[N(SiMe3)2]2}2 (376.5 mg, 0.5 mmol) in mesitylene (10 ml) in a Ficher Porter bottle. A solution of hexadecylamine (483 mg, 2 mmols) in mesitylene (5 ml) was then added to this mixture. Its color turned to black. The reacting medium was pressurized under 3 bars dihydrogen then heated at 150°C for 48 hours. Excess dihydrogen was evacuated and the solution was let to cool at room temperature. The magnetically collected black powder was washed three times with toluene (15 ml). It yielded a material that contained xx% and 24.1% of iron in samples 1 and 2 respectively. 14 Synthesis of samples 3 and 4: Hexadecylamonium chloride (277 mg, 1 mmol and 415,5 mg, 1,5 mmols) was added to a green solution of {Fe[N(SiMe3)2]2}2 (376,5 mg, 0.5 mmol) in mesitylene (15 ml) in a Ficher Porter bottle. The color changed immediately from green to yellow. A solution of hexadecylamine (483 mg, 2 mmols in both cases) in mesitylene (5 ml) was then added to this mixture. Its color turned to black. The reacting medium was heated at 90°C for 30 minutes, cooled down to room temperature, then pressurized under 3 bars dihydrogen and heated again at 150°C for 48 hours. Excess dihydrogen was evacuated and the solution was let to cool down at room temperature. The magnetically collected black powder was washed three times with toluene (15 ml). It yielded a material that contained 70.4% and 68.8% of iron for samples 3 and 4 respectively. Synthesis of samples 5, 6 and 7: Ultra-small iron NPs (1 mmol of Fe) were prepared in 10 mL of mesitylene. Hexadecylammonium chloride (304,5 mg, 1.1 mmols for 5; 277 mg, 1 mmol for 6; 554 mg, 2 mmols for 7) and hexadecylamine (483 mg, 2 mmols in the three cases) were successively added to the NPs solution. The reacting medium was heated at 90°C for 30 minutes, cooled to room temperature, pressurized under 3 bars dihydrogen and then heated again at 150°C for 48 hours. Excess dihydrogen was evacuated and the solution was let to cool down at room temperature. The magnetically collected black powder was washed three times with toluene (15 ml). It yielded a material that contained 81.1%, 75.6% and 91.2% of iron in samples 5, 6 and 7 respectively.

Magnetic measurements:

The magnetization and the coercive field of as-synthesized iron(0) NPs were measured on a Quantum Design Model MPMS 5.5 SQUID magnetometer. The absolute magnetization was deduced from the iron(0) total content determined by Inductively Coupled Plasma -Mass Spectrometry. These studies were carried out on powder samples that were prepared and sealed under an argon atmosphere to preserve the metallic character of iron. Extreme care was taken to avoid oxidation during the transfer to the apparatus.

Hyperthermia measurements: For hyperthermia measurements, a schlenk containing about 12 mg of powder of Fe NPs and 0.5 ml of mesitylene was filled under inert atmosphere to prevent any oxidation of the NPs. The mesitylene was used as solvent to ensure the solubility of the nanoparticles and because of his high boiling point. The schlenk is then placed in a calorimeter with 1.5 ml of deionised water, the temperature of which was measured. The measurement time was varied between 30 and 100 s, depending on the experimental 15 parameters, so that the temperature rise never exceeds 20°C. The temperature rise at the end of the magnetic field application was always measured after shaking the calorimeter to ensure the temperature homogeneity, which was checked by putting two probes at the top and the bottom of the calorimeter. The temperature rise is determined after this process from the mean slope of the ∆T/∆t function. Then the SAR values were calculated using the expression: for Fe NPs, Cp = 1750 J kg -1 K -1 for mesitylene, Cp = 4186 J kg -1 K -1 for water and Cp = 720 J kg -1 K -1 for glass), and mFe is the mass of the pure Fe MNPs. The mass of glass is put to 2930 mg. Dividing SAR values by the magnetic field frequency lead to the specific losses A. [a] α values for measurements at µ0Hmax = 60 mT and f = 54 kHz.

[b] α values for measurements at µ0Hmax = 73 mT and f = 274 kHz. Heating power of metallic iron magnetic nanoparticles whose diameter ranges from 5.5 nm to 28 nm is measured in an alternative applied magnetic field of frequency of 54 kHz and of amplitude ranging from 8 to 60 mT. The optimal size is visible as a maximum in the curve.

The table of contents entry

Optimized nanoparticles display losses of 8 mJ g -1 at this frequency.

  Figure 5c. For large MNPs, the size distribution has a negligible influence, while for small MNPs, the size distribution increases the SAR values. Thus, the size distribution effect can not be held responsible for the smaller experimental SAR values.

Where

  Cpi and mi are specific heat capacity and mass for each component (Cp = 449 J kg -1 K -1

Figure 1 .

 1 Figure 1. Transmission electron microscopy image of the different samples and corresponding size distribution. The histogram was fitted with a log-normal distribution to estimate the mean particle size d0 and the standard deviation σ.

Figure 2 .

 2 Figure 2. Hyperthermia measurement at frequency f = 54 kHz (left) and f = 274 kHz (right).

  (a) and (b) Evolution of the SAR as a function of magnetic field for different samples. (c) and (d) Normalized SAR value for NPs ranging between 8.9 and 27.5 nm. (e) and (f) Size dependence of SAR at different magnetic fields. (g) and (h) Data of (a) and (b) modified using Equation (1) to plot the magnetic field dependence of α.

Figure 3 .

 3 Figure 3. Numerical simulations of the hysteresis area using the parameters of bulk iron. Keff = 4.8×10 4 J m -3 , τ0 = 5×10 -11 s, MS = 1.7 kA m -1 , T = 300 K, f =100 kHz. (a) Evolution of the hysteresis loops as a function of the magnetic field for d = 16 nm. (b) Evolution of area as a function of magnetic field for different sizes. (c) Normalized hysteresis area for different size.

Figure 4 .

 4 Figure 4. (a) Evolution of the hysteresis loop as a function of size for µ0Hmax = 20 mT. (b) Evolution of hysteresis area as a function of size NPs at different magnetic field. (c)Comparison between the optimal diameters calculated analytically using Equation (2) (solid line), Equations (3), (4), (5) (dashed line) , and the ones extracted from numerical simulations (star points).

Figure 5 .

 5 Figure 5. Quantitative analysis of data at f = 274 kHz (a) Evolution of µ0HCHyp determined from experiments (round dots), from a fit using Equation (6) (solid line), and from the numerical simulations (triangles). Ms = 1.7×10 6 A m -1 , T = 300 K, f = 274 kHz, τ0 = 5×10 -11 s and Keff = 9.1×10 4 J m -3 . (b) SAR (µ0Hmax) calculated numerically for four diameters using these parameters. (c) Numerical calculations of the SAR at µ0Hmax = 73.5 mT as a function of the diameter when neglecting the effect of size distribution (dashed line), or using the size distribution deduced from TEM (round dots). Corresponding experimental values (square points). (d) Normalized SAR (µ0Hmax). Dots represents experimental values and lines are results from numerical simulations assuming no size distribution.

Table 1 .

 1 Summary of magnetic and hyperthermia properties of the different samples: mean diameter d0, standard deviation σ, saturation magnetization σS, remnant magnetization σR, ratio σR / σS, coercive field µ0HC, losses per cycle A at f = 54 kHz, 60 mT, and at f = 274 kHz, 73 mT, and corresponding α values calculated using Equation (1).

  should be fifty to sixty words long, written in the present tense, and refer to the chosen figure. Magnetic Nanoparticles, Biomedical Applications, Fe Nanoparticles By B. Mehdaoui, A. Meffre, J. Carrey*, S. Lachaize, L. M. Lacroix, M. Gougeon, B. Chaudret ,and M. Respaud Optimal size of nanoparticles for magnetic hyperthermia: A combined theoretical and experimental study ToC figure ((Please choose one size: 55 mm broad × 50 mm high or 110 mm broad × 20 mm high. Please do not use any other dimensions))

  is around one week, it is not reasonable to expect to vary both τ0 and Keff in numerical simulations until a good fit of experimental data is achieved. To make a quantitative analysis possible, a special method must be used and a few assumptions made. The method is the following: i) a "coercive field value deduced from MH experiments"

		µ 0	H	CHyp	is defined as the point of
	highest slope in SAR(µ0Hmax) functions. The interest of	µ 0	H	CHyp

Acknowledgements

This work was supported by the InNaBioSanté Foundation. Supporting Information is available online from Wiley InterScience or from the author.

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

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

Published online: ((will be filled in by the editorial staff))


L Lilensten

Y Danard

C Brozek

S Mantri

P Castany

T Gloriant

P Vermaut

F Sun

R Banerjee

F Prima

On the heterogeneous nature of deformation in a strain-transformable beta metastable Ti-V-Cr-Al alloy

Keywords: Titanium alloy, deformation twinning, martensitic phase transformation, strain hardening, deformation structures

Ti-10V-4Cr-1Al wt% (TVCA) is a new grade of titanium alloy, developed to combine twinning induced plasticity (TWIP) and transformation induced plasticity (TRIP) effects. The TVCA alloy exhibits a very high strain-hardening rate and an excellent balance between strength and ductility for great potential in aerospace applications.

Deformation mechanisms are investigated using in-situ techniques as synchrotron X-ray diffraction (SXRD) and in-situ electron backscatter diffraction (EBSD) analysis during

Introduction

The development of β-titanium alloys has attracted an increasing attention over the last few years for fulfilling the growing demand from the aerospace industry. Until recently, the interest brought to this family of titanium alloys was limited because of their rather weak work hardening and ductility, specific to bcc structures. However, proof was given that triggering twinning and/or martensitic phase transformation in the β-phase under

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

3 applied stress, even for dual-phase alloys, could greatly enhance both the work hardening capability as well as the elongation of the alloy, leading to materials with great potential for various applications [START_REF] Sun | A new titanium alloy with a combination of high strength, high strain hardening and improved ductility[END_REF][START_REF] Marteleur | On the design of new β-metastable titanium alloys with improved work hardening rate thanks to simultaneous TRIP and TWIP effects[END_REF][START_REF] Brozek | A β-titanium alloy with extra high strain-hardening rate: Design and mechanical properties[END_REF][START_REF] Kim | Martensitic transformation, shape memory effect and superelasticity of Ti-Nb binary alloys[END_REF][START_REF] Grosdidier | Deformation induced martensite and superelasticity in a β-metastable titanium alloy[END_REF][START_REF] Min | Enhancement of uniform elongation in high strength Ti-Mo based alloys by combination of deformation modes[END_REF][START_REF] Min | Mechanism of twinning-induced plasticity in β-type Ti-15Mo alloy[END_REF][START_REF] Talling | On the mechanism of superelasticity in Gum metal[END_REF][START_REF] Ahmed | The influence of β phase stability on deformation mode and compressive mechanical properties of Ti-10V-3Fe-3Al alloy[END_REF]. It is now well known that the transformation induced plasticity (TRIP) and the twinning induced plasticity (TWIP) effects are closely related to the chemical stability of the β-phase at room temperature [START_REF] Kolli | Phase Stability and Stress-Induced Transformations in Beta Titanium Alloys[END_REF]. The stability domains of the later can be predicted with good accuracy using the Bo-Md (where Bo stands for mean bond order and Md for mean d-orbital energy) design method [START_REF] Marteleur | On the design of new β-metastable titanium alloys with improved work hardening rate thanks to simultaneous TRIP and TWIP effects[END_REF][START_REF] Abdel-Hady | General approach to phase stability and elastic properties of β-type Ti-alloys using electronic parameters[END_REF], DFT [START_REF] Raabe | Theory-guided bottom-up design of β-titanium alloys as biomaterials based on first principles calculations: Theory and experiments[END_REF], but also Calphad approach. They allow to target the boundaries of the stable β domains, where TRIP and TWIP are more likely to occur. The understanding of the deformation mechanisms, taken independently, benefits from several studies. Concerning the TRIP effect, a good understanding of the formation of stress induced martensite was reached with the study of shape memory alloys, superelastic alloys and gum-metals, used for biomedical applications [START_REF] Kim | Martensitic transformation, shape memory effect and superelasticity of Ti-Nb binary alloys[END_REF][START_REF] Grosdidier | Deformation induced martensite and superelasticity in a β-metastable titanium alloy[END_REF][START_REF] Talling | On the mechanism of superelasticity in Gum metal[END_REF][START_REF] Besse | Mechanisms of deformation in gum metal TNTZ-O and TNTZ titanium alloys: A comparative study on the oxygen influence[END_REF][START_REF] Castany | In situ synchrotron X-ray diffraction study of the martensitic transformation in superelastic Ti-24Nb-0.5N and Ti-24Nb-0.5O alloys[END_REF]. As for the twinning effect, it has been investigated on various aspects, such as the influence of the grain orientation [START_REF] Bertrand | Twinning system selection in a metastable β-titanium alloy by Schmid factor analysis[END_REF][START_REF] Min | Deformation microstructural evolution and strain hardening of differently oriented grains in twinning-induced plasticity β titanium alloy[END_REF][START_REF] Zhou | Accommodative 332〈113〉 primary and secondary twinning in a slightly deformed β-type Ti-Mo titanium alloy[END_REF], deformation rate [START_REF] Ji | Strain-rate effect on work-hardening behavior in β-type Ti-10Mo-1Fe alloy with TWIP effect[END_REF] and chemical homogeneity [START_REF] Min | Heterogeneous twin formation and its effect on tensile properties in Ti-Mo based β titanium alloys[END_REF].

Twinning mechanisms for the bcc structure are still in hypotheses among various models, especially regarding the {332}<113> twinning, theoretically unfavorable yet widely observed in metastable β-Ti-alloys [START_REF] Crocker | Twinned martensite[END_REF][START_REF] Tobe | Origin of 332 twinning in metastable β-Ti alloys[END_REF][START_REF] Kawabata | Mechanical properties of TiNbTa single crystals at cryogenic temperatures[END_REF][START_REF] Lai | On the mechanism of 332 twinning in metastable β titanium alloys[END_REF]. Recent assumptions based on experimental observations suggest that the {332}<113> twinning process could be mediated by the formation of orthorhombic martensite, formed as a martensitic reversion product [START_REF] Lai | On the mechanism of 332 twinning in metastable β titanium alloys[END_REF][START_REF] Chen | Transitional structure of 332〈113〉β twin boundary in a deformed metastable β-type Ti-Nb-based alloy, revealed by atomic resolution electron microscopy[END_REF][START_REF] Castany | Reversion of a Parent {130}⟨310⟩α" Martensitic Twinning System at the Origin of {332}<113>β Twins Observed in Metastable β Titanium Alloys[END_REF]. This shows the complexity of the interactions between the

M A N U S C R I P T A C C E P T E D

ACCEPTED MANUSCRIPT 4 various phenomena, and how interdependent they could be. Thus, even more complex than single TWIP or single TRIP effects, the interplay of the two deformation modes in TWIP/TRIP alloys, although leading to highly desirable mechanical properties, remains far from being well understood; even if some first attempts to describe the deformation sequence and the interactions between the various features (including the role of the stress induced ω-phase, not observed in this study) have been reported [START_REF] Sun | A new titanium alloy with a combination of high strength, high strain hardening and improved ductility[END_REF][START_REF] Sun | Investigation of early stage deformation mechanisms in a metastable β titanium alloy showing combined twinning-induced plasticity and transformation-induced plasticity effects[END_REF][START_REF] Zhang | Complexion-mediated martensitic phase transformation in Titanium[END_REF][START_REF] Lai | Origin of shear induced β to ω transition in Ti-Nb-based alloys[END_REF]. Yet, no agreement on the deformation sequence was found, and microstructural mechanisms are derived mainly based on post-mortem imaging, limiting the understanding of the complex interactions which are very sensitive to stress-strain environment [START_REF] Sun | Investigation of early stage deformation mechanisms in a metastable β titanium alloy showing combined twinning-induced plasticity and transformation-induced plasticity effects[END_REF][START_REF] Zhang | Microstructural evolution of a ductile metastable β titanium alloy with combined TRIP/TWIP effects[END_REF][START_REF] Gao | Deformation mechanisms in a metastable beta titanium twinning induced plasticity alloy with high yield strength and high strain hardening rate[END_REF]. This study focuses on the in-situ investigation of the deformation mechanisms at the various length scales, and their interaction, using a Ti-V-Cr-Al alloy designed from d-electron approach [START_REF] Marteleur | On the design of new β-metastable titanium alloys with improved work hardening rate thanks to simultaneous TRIP and TWIP effects[END_REF][START_REF] Abdel-Hady | General approach to phase stability and elastic properties of β-type Ti-alloys using electronic parameters[END_REF]. The composition is designed for the formation of TRIP and TWIP upon deformation, using Bo (mean bond order) and Md (mean d-orbital energy) parameters of 2.785 and 2.353 respectively [START_REF] Marteleur | On the design of new β-metastable titanium alloys with improved work hardening rate thanks to simultaneous TRIP and TWIP effects[END_REF][START_REF] Abdel-Hady | General approach to phase stability and elastic properties of β-type Ti-alloys using electronic parameters[END_REF][START_REF] Brozek | Turbomachine part comprising an alloy based on titanium[END_REF]. In-situ analyses are of particular interest in this case, considering the partially reversible (when tensile stress is released) behavior of the martensite. A special focus is given to the understanding of the criteria ruling the selection of the deformation mechanisms, to try to add to the existing results on the "deformation sequence" in such alloys. The macro-scale data (in-situ synchrotron X ray diffraction SXRD) and the microstructural analyses (in-situ SEM and EBSD, with supplementary TEM) are used to investigate the evolution of the martensitic and the β-twins fractions.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

5

The triggering of one and/or the other mechanism -among TRIP and TWIP -for the various grains, using a Schmid factor (SF) analysis, is used to explain the exceptional tensile properties of the material.

Experimental

Ti-10V-4Cr-1Al alloy (TVCA) in wt% is processed by arc melting under a high Ar purity atmosphere. A 200g ingot is melted using industrial purity Ti sponges, Cr granules, a master alloy of AlV and Al shots. Plates of 3mm thickness are cut out, homogenized at 1173K under air for 900 s and then water quenched. Slice-cut samples are cold rolled down to 0.5mm sheets, corresponding to a reduction level of 80%. To restore a fully β state, cold rolled sheets are finally solution treated (ST) under high purity Ar atmosphere at 1033 K for 900 s, and then water quenched. Tensile tests are performed at room temperature using an INSTRON5966 machine with an extensometer, on specimen with dimensions of 50 mm x 5 mm x 0.5mm, at a strain rate of 10 -3 s -1 .

In-situ SEM tensile samples, with dimensions 35 mm x 2mm x 0.5mm, are prepared by electropolishing. Electron backscatter diffraction (EBSD) scans are collected using a field emission gun scanning electron microscope operating at 20kV, with step size ranging from 0.5 to 0.01µm. In-situ tensile tests are performed using a Proxima 100-Micromecha machine. Deformation is followed by secondary electron (SE) imaging during the loading state and EBSD scans are recorded at different strains.

Because of some surface relief related to the activated deformation mechanisms, preventing indexation, a few microns of the surface at the rupture stage were In situ SXRD is conducted at the high resolution beam line ID22 of the European Synchrotron Radiation Source (ESRF), Grenoble, France, from a tensile sample with gauge width of 4mm and 0.5mm in thickness. The incident X-ray beam had a diameter of 1 mm and a wavelength of 0.040002106 nm. Data collection of transmitted diffracted beams was performed over the 2θ angular range 6-24°, with a step size of 0.005°. Six scanning stages are carried out at room temperature, starting from the unloaded state, and then for increasingly loaded states, up to a strain of 0.20 with a strain rate of 10 -4 s -1 .

Analysis of the diffraction patterns consists in identifying the phases corresponding to the different peaks, and then estimating the lattice parameters owing to Pawley refinement using the TOPAS software. Estimation of the intensity of specific peaks is also carried out. Site-specific samples for transmission electron microscopy (TEM) were prepared using a FEI Nova NanoLab 200™ focused ion beam (FIB) system. TEM studies were carried out using a FEI Tecnai F20-FEG operating at 200 kV and FEI Tecnai G2 30 TWIN operating at 300 kV.

Results

Mechanical performances of Ti-V-Cr-Al alloy

The mechanical behavior is shown in form of a uniaxial tensile curve in figure 1a, in which a stable plastic flow starts at about 420MPa and is maintained until 1200MPa.

Initiation of localized necking is observed at about ε=0.35, yet the total engineering elongation exceeds 50%. Consistent with the Bo/Md design strategy, the TVCA alloy

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 7
displays a complex, non-monotonic, stress-strain profile previously observed in both Ti-12Mo and Ti-8.5Cr-1.5Sn alloys, especially due to the hump in the work-hardening rate curve (figure 1b). The work hardening rate of TVCA reaches a maximum of about 2.6 GPa at ε=0.12. The vertex divides the work hardening process into two stages: early stage -from yielding to ε=0.12 -and late stage from 0.12 to 0.35. The work hardening rate is characterized by an increase in the early stages before to slowly decrease during the late stage.

In situ SXRD investigation of deformation mechanisms

Before deformation (ε=0), the single β phase (body center cubic bcc) with cell parameter a=3.26 Å is confirmed by SXRD (Figure 2). SIM is, as expected, triggered at the end of the elastic regime. The orthorhombic α" martensitic phase is indeed detected by the first SXRD analysis in tension, for a deformation of ε=0.02 (Fig. 2a), with lattice parameters a=3.06 Å, b=4.88 Å and c=4.62 Å determined by Rierveld refinement. The evolution of the volume fraction of α" is estimated from the quantitative Rietveld refinement results for each diffraction diagram, assuming that the texture is random (Figure 2b). Although this last hypothesis is a bit strong, it allows us to access approximately the amount of α'', as well as to get the relative evolution of the α'' phase fraction with strain. The volume fraction increases rapidly to about 30 vol% at at ε≈0.08, and then reaches a plateau, during which the amount increases very little, to reach approximately 35 vol% at ε=0.2.

Local analysis of the deformation mechanisms: a heterogeneous behavior

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 8

At the macro-scale, an EBSD analysis of the microstructure (Figure 3) reveals a 100% β structure prior deformation, consistently with the SXRD information, and an average grain size of 46 microns. The color code refers to the calculated SF of the {332}<113> twinning (in red) and to the one of the SIM formation, in green, calculated using the habit plane of the martensite. One can see that the grains are grouped, defining some zones that are more likely, based on the SF, to deform either by twinning or by formation of martensite. It can also be noted that when the grains are well oriented for one mechanism, the other mechanism is rather unlikely to occur. Therefore, further analyses were performed to identify the respective sequence of deformation mechanisms in the two populations of grains. The two populations are referred to in the following as "type A" grains, displaying high SF for SIM (TRIP), and "type B" grains displaying high SF for mechanical twinning (TWIP effect),

Deformation induced martensitic phase transformation ("Type A" grain)

For further investigation, a model "type A" grain is selected, and called grain A in the following. Plate-like deformation bands (Fig. 4a) are observed by SEI (secondary electron imaging) after in-situ tensile loading to ε=0.03. EBSD mapping identifies the nature of the bands as α" martensite. It can be seen from Fig. 4a that one of the bands is composed of two orientations of α" that are in twinning relationship. The orientation relationships between the α" phase variants and β matrix are verified by pole figure projection (Fig. 4c and4d). The classical Burgers relationship (110) α" // [START_REF] Marteleur | On the design of new β-metastable titanium alloys with improved work hardening rate thanks to simultaneous TRIP and TWIP effects[END_REF][START_REF] Brozek | A β-titanium alloy with extra high strain-hardening rate: Design and mechanical properties[END_REF][START_REF] Kim | Martensitic transformation, shape memory effect and superelasticity of Ti-Nb binary alloys[END_REF][START_REF] Grosdidier | Deformation induced martensite and superelasticity in a β-metastable titanium alloy[END_REF][START_REF] Min | Enhancement of uniform elongation in high strength Ti-Mo based alloys by combination of deformation modes[END_REF][START_REF] Min | Mechanism of twinning-induced plasticity in β-type Ti-15Mo alloy[END_REF][START_REF] Talling | On the mechanism of superelasticity in Gum metal[END_REF][START_REF] Ahmed | The influence of β phase stability on deformation mode and compressive mechanical properties of Ti-10V-3Fe-3Al alloy[END_REF][START_REF] Kolli | Phase Stability and Stress-Induced Transformations in Beta Titanium Alloys[END_REF][START_REF] Abdel-Hady | General approach to phase stability and elastic properties of β-type Ti-alloys using electronic parameters[END_REF]

β is found M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 9
for the major SIM variant (green), but not found for purple one, confirming that the purple variant is obtained by twinning of the primary martensite (green variant). After releasing the tensile load, the martensite disappears (Fig. 4b), and a full β structure is restored, confirmed by EBSD mapping on the same zone. It shows that the martensitic phase transformation can be reversible in the early stage. In contrast, no mechanical twinning of beta phase was observed in the area for these deformation levels.

Stress induced twinning and secondary SIM ("Type B" grain)

In-situ SEM investigations are then performed to monitor the deformation of a grain, showing a high SF for mechanical twinning, denoted by Grain B hereafter.

SEI show that at yield (Fig. 5b), two straight parallel bands appear on the surface of Grain B, growing about 76° away from the tensile direction (vertical). At lower deformation (up to ε≈0.2), deformation seems to proceed by band multiplication, before thickening of the bands takes over.

Deformation strain is held at ε≈0.1 and ε≈0.2 for EBSD mapping, and a final EBSD is performed on the fractured sample (ε≈0.35). Phase maps are over imposed on the SE images (figure 5), and details of the orientation maps at the different stages can be found in the supplementary materials. Difficulties of indexation are noticed at the interfaces between bands and matrix, probably due to the shadow of the relief or to very fine and strained microstructures. The in-depth analysis of the nature of the deformation products, their orientation and surface fraction in the sample after rupture is then performed, and the results are given in Fig. 6a. For an easier identification, they are As can be seen in figures 5 and 6, the strain-induced deformation bands are mainly identified as β twins belonging to the {332}<113> twinning system (three twinning systems T1, T2 and T3, for a total surface fraction of 27.2%). A major twinning system T1 is observed, and one {112}<111> twin (T4) is formed at the late stage of the deformation (Figure 1 of supplementaty material), with low surface fraction of 0.7%. Some martensite is detected too, mostly on the right side of the grain. Three α" variants are found, with surface fractions of 2.3%, 1.3% and 6.7%. The a3 variant, it is the only one detected only after rupture (Figure 1 of supplementary materials), unlike the other variants which are identified for lower strains.

Considering the indexed parts, the {332} twins possess a higher volume fraction than the martensitic phase. The martensite surface area is also much lower than macroscopic SXRD estimation (Fig. 2b). Yet, this is a local result, and an analysis of this zone with lower magnification leads to 22.9% (surface fraction) of martensite at failure, closer to the SXRD results (Figure 2 of supplementary material). By plotting the stereographic projections of the different features (Figure 6c-e), classical crystallographic relationships are confirmed for the major T1 and the matrix, and for a1 variant and T1 but not for a1 and the matrix. This leads to the hypothesis that a1 is a secondary deformation product obtained from the deformation of T1 and not the matrix. Similar results are found for a2, that does not fulfill the Burgers relationship with the matrix but with T1 (not displayed here).

M A N U S C R I P T A C C E P T E D

Detailed microstructures of twins and martensite

An in-band microstructure is studied by TEM, to try to better understand the complexity of the secondary deformation mechanisms and the underlying phenomena. Grain B is selected, and the lift-out is performed after fracture, on a band with low confidence of index (black pixels in Fig. 7a). Figure 7 shows the cross section of the band, where a step is observed at the surface. It is evidenced from the SAED patterns (Fig. 7b) that the microstructure is composed of a {332} twin with internal secondary α" phase. However, the α" phase seems unstable, since SAED pattern on the left of the Figure 7b, taken on the same area after six months of "aging" at room temperature, shows that the α" phase is actually reversed into its parent phase, the β {332} twin.

Schmid factor analysis

Details of the SF, that led to the selection of grains A and B as model grains, are given in Figure 8. It is completed with analysis of the experimentally observed deformation features (highlighted) and details of T1 are also given, as this system was shown to display internal secondary deformation features. The SF for the three best-oriented variants (highest SF) of {332} twins, {112} twins and α'' martensite are given. In all the cases, there is a {112} twinning system which displays the highest SF of all systems.

However, such twins are not observed, except in grain B, for large strains (detected only in the last EBSD mapping). Apart from this twinning system, the mechanisms triggered in grains A and B fulfill the SF selection criterion: the formation of martensite is the favored mechanism in grain A and is indeed experimentally observed, similar to the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 12 {332}
twinning system for grain B, with the highest SF of 0.459 for the main variant T1 [131]. The fact that the SF of martensite is not the highest one for the secondary deformation mechanism in T1, and yet it is the mechanism observed, may indicate that the SF rule is not respected for secondary deformation mechanisms. This is in line with the idea that secondary deformation mechanisms accommodate the mechanical contrast inside the grains, without any direct link with the macroscopic strain.

Correlation between the mechanical data and the analysis of the deformation mechanisms

Finally, the comparison between strain hardening rate and strain transforming products at two different scales is shown in Figure 9. In Figure 9a, results from macroscopic measurements, strain hardening from tensile test and α'' fraction from in-situ SXRD, are plotted as a function of tensile strain. The plateau reached by the martensitic fraction happens to correspond to the hump in the work-hardening rate. Microscopic evolutions of the area fractions of β twins and α" martensite in grain B (TWIP/TRIP mode) are presented in Figure 9b as a function of tensile strain. The increase in the twinned fraction is monotonic up to about 25 area% at rupture. In contrary the α'' fraction is very low and increases slowly until 20% and then strongly at the fractured stage. The difference between macroscopic and microscopic observations may confirm the preferential selection of local deformation mechanism according to the grain orientation.

Discussion

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 13

The tensile behavior coupled with in-situ analysis of the deformation mechanisms of a β-metastable Ti-alloy are reported in this paper. As expected from the shape of the tensile curve, a signature of strain transformable alloys with deformation processes related to TRIP and TWIP, stress induced martensite and twins are formed upon deformation, as indicated by the SXRD (Figure 2), SEM (EBSD) (Figure 456) and TEM analyses (Figure 7). Yet, an EBSD analysis at the macro-scale prior deformation suggests that, depending on their orientation and so their Schmid factor, some grains may be more favorable for TRIP and some other for TWIP. This could indicate that, rather than the existence of a unique deformation sequence in the alloy, each type of grain undergoes a different pattern. Therefore, the following discussion will focus on the "type A" and the "type B" ones. The different roles of the martensite, in the earlier deformation stage and in the late one, will be also discussed. Finally, the impact of this contrast, between the two populations of grains, will be discussed in the light of the mechanical properties and possible hardening mechanisms.

Investigation of the different behaviors observed in type A and type B grains

Selection of the primary deformation mechanism

The textural influence on mechanical twinning has been highlighted in several papers, through Schmid factor (SF) analysis [START_REF] Bertrand | Twinning system selection in a metastable β-titanium alloy by Schmid factor analysis[END_REF][START_REF] Min | Deformation microstructural evolution and strain hardening of differently oriented grains in twinning-induced plasticity β titanium alloy[END_REF][START_REF] Zhou | Accommodative 332〈113〉 primary and secondary twinning in a slightly deformed β-type Ti-Mo titanium alloy[END_REF]. However, this particular parameter has not been evaluated in case of potentially combined TRIP and TWIP effects, with subsequent questions on the twins/martensite interactions. As can be seen in Figure 2, the β-grains at the initial state seem to be, based on the calculation of the SF, well

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

14 oriented for one or the other mechanism, either TRIP or TWIP. Experimentally, it is confirmed that SF predicts with good accuracy which grain will TRIP or TWIP (at least as primary deformation mechanism). "Type A" grains, as the one investigated in figure 4, shows that the first deformation mechanism observed is TRIP, in good agreement with the calculation of the SF (Figure 8). In this grain, up to ε=0.03, no twinning of the β-phase is found.

Similarly, in the grain investigated in figures 5-7 ("type B" grain), TWIP is the primary deformation mechanism, with the formation and multiplication of {332}<113> twins in the first stage. In grain B, the detected martensite is identified as a secondary deformation product (discussed later in section 4.1.2). At the end of the deformation process, up to 30% (surface fraction) of the grain is twinned, and less than 10% is indexed as martensite, well below the average value of 35% calculated on the macroscopic level by SXRD. This highlights the fact that different behaviors are found at the grain-size level, and the macroscopic data only cannot reflect on the variety and complexity of the local behaviors. Therefore, the history of the deformation, applying in every grain, cannot be treated as unique, but must be split into at least two stories, for the type A and type B grains in the alloy. The fact that both mechanisms can be triggered as primary deformation mechanisms also indicates that the critical resolved shear stresses (CRSS) required for TRIP and TWIP are similar, as is classically expected for TRIP/TWIP alloys. Indeed, if one of the mechanisms had a much lower CRSS, it would be triggered regardless of the grain orientation.

M A N U S C R I P T A C C E P T E D

Formation of secondary deformation products

When increasing the strain in grain B (figures 5-7), deformation products other than {332}<113> twinning are observed: {112}<111> twinning is found, as well as martensite variants. The fact that {112}<111> twinning is not observed earlier in the alloy (see supplementary material for details), in spite of its high SF (Figure 8) may be explained by a possible high CRSS, that prevented its earlier formation. Regarding the martensite, the indexed variants a1 and a2 do not fulfill the Burgers relationship with the matrix but with T1 (Figure 7d), therefore they can be identified as a secondary deformation mechanism.

Taking a closer look at the location where the martensite is formed, it seems to mostly be coupled with other features, such as a grain boundary or a {332} twin, with sometimes a very intricate structure for the latter case. This suggests that martensite is triggered at late deformation stages in the highly strained areas, to accommodate strong local strain fields.

Bands composed only of martensite (fully martensitic needles) are not found in this grain. Contrary to a1 and a2, detected at ε=0.1 and ε=0.2 in association with T1, a3 is observed in grain B only after rupture. This green variant, detected in several zones of grain B, has a twinning relationship with a1 and a2, and also a negative SF. This leads to the hypothesis that a3 is formed by twinning of a2 (previously indexed in the same zone in the grain), aiming at accommodating the backstress in the sample.

On the {332}<113> twinning

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 16

The formation of {332}<113> twins has been subjected to a lot of hypothesis, including shuffling, or martensitic mediation [START_REF] Lai | On the mechanism of 332 twinning in metastable β titanium alloys[END_REF][START_REF] Chen | Transitional structure of 332〈113〉β twin boundary in a deformed metastable β-type Ti-Nb-based alloy, revealed by atomic resolution electron microscopy[END_REF][START_REF] Castany | Reversion of a Parent {130}⟨310⟩α" Martensitic Twinning System at the Origin of {332}<113>β Twins Observed in Metastable β Titanium Alloys[END_REF]. Yet, the controversy remains, as no "universal mechanism" is reached, each alloy presenting its own characteristic features.

For this alloy, to get a better understanding of how the defects form and interact together, one of the band, that is not easily indexed, of the fractured sample, was looked at by TEM. One deformation step before, parts on the right side of the band was indexed as martensite (a1), and the left side of the band was indexed as T1.

Bright field images show that the band is actually composed of several components stacked on top of each other. The central and largest band was unambiguously indexed as orthorhombic martensite in (310) zone axis. It is surrounded by one band with {332}<113> twinning orientation relationship on each side, making the transition with the matrix. This result is the opposite of what was observed in several papers [START_REF] Lai | On the mechanism of 332 twinning in metastable β titanium alloys[END_REF][START_REF] Chen | Transitional structure of 332〈113〉β twin boundary in a deformed metastable β-type Ti-Nb-based alloy, revealed by atomic resolution electron microscopy[END_REF][START_REF] Castany | Reversion of a Parent {130}⟨310⟩α" Martensitic Twinning System at the Origin of {332}<113>β Twins Observed in Metastable β Titanium Alloys[END_REF], where the boundary between a {332}<113> twin and the matrix was made of orthorhombic martensite, more or less well defined. Results from the literature lead to the hypothesis that the formation of {332}<113> twins could be mediated by martensitic phase transformation. In the present case, several hints suggest that this scenario is not happening. First, the martensite is contained in the twin. Besides, the interface between the martensite and the twin is not parallel to the interface between the twin and the matrix. This indicates that the twin formed first, and that the martensite is a deformation sub-product, in good agreement with the EBSD data. Finally, the sample was analyzed

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

again after 6 months at room temperature ("aged" diffraction pattern in figure 7). The new TEM analysis of the "aged" sample shows that the martensite needle has transformed back into a bcc structure, the latter having the same orientation as the neighboring twin, and not the matrix. The martensite therefore reverted into its parent structure: the {332}<113> twin.

The deformation sequence observed here does not aim at being universal, since it is important to keep in mind that all the Ti β-metastable alloys with different compositions have a different sensitivity for the TRIP and the TWIP mechanism. But it shows that the mechanisms suggested, the martensitic mediation of twinning, although undeniable in the reported cases, also cannot be extended to the {332}<113> in general.

Understanding the deformation mechanisms

Correlation between the microstructural investigations and the tensile testing curve

Figure 1 shows that the tensile behavior can be seen as made of two different domains.

The first one goes up to about ε≈0.12, and the second one goes from that point on. In the first stage, the work-hardening rate increases a lot, and decreases in the second one.

The first deformation stage is also marked by a strong increase of the quantity of martensite, as indicated by the SXRD data (figure 2), showing a raise from 0% (initial state, undeformed sample) to about 30%. The shape of the tensile curve is characteristic of alloys forming SIM with the double yielding, observed between ε≈0.055 and ε≈0.12.

It is usually considered that the second yielding corresponds to the elastic deformation

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

18 of the martensite that was formed during the first stage (here up to ε≈0.055). This is in good agreement with the fact that according to SXRD data, the fraction of martensite in the material reaches a plateau at about 10% of deformation, which corresponds to the beginning of the hump in the work-hardening rate. Yet, as the present alloy is not purely superelastic but TRIP/TWIP, the mechanisms associated with this step combines elastic deformation of the martensite, but also twinning and dislocation glide.

The first deformation domain therefore seems to be controlled by two "independent" primary deformation mechanisms, twinning and SIM, each of them occurring in different grains (type A and type B).

During the second deformation stage, the amount of martensite increases very little.

Considering this slight increase and the discussion of the previous part, it can be assumed that the martensite formed during that stage is a secondary deformation product, and is formed to accommodate local stress fields, explaining the much smaller quantity. This also explains the rather heterogeneous repartition of the martensite in the grain, mostly localized in the twins or at the grain boundaries, which are locations expected to be under higher strain levels.

Towards a two-levels composite effect

The study of the grains A and B brought a lot of information on the various behaviors observed in the grains. The first major point is that there is a heterogeneity between the grains depending on their orientation. Some of them are more favorably oriented for TRIP and some others for TWIP, therefore grains deforming by TWIP are lying next to

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

19 grains deforming by TRIP. However, the consequence of one or the other mechanism is not the same from a mechanical point of view. In type-B grains, the refinement of the microstructure by twinning creates hardening via the dynamic Hall-Petch law [START_REF] Remy | Kinetics of f.c.c. deformation twinning and its relationship to stress-strain behaviour[END_REF].

Twin boundaries act as barriers to gliding dislocations, consequently reducing their mean free path, and these grains harden. In comparison, the role of martensite is still under discussion, as it is classically considered as a soft phase, with a modulus lower than the one of the bcc phase. Besides, as was shown in a recent study, martensite may not act as a strong barrier for gliding dislocations, in comparison with twin boundaries [START_REF] Lai | ω phase acts as a switch between dislocation channeling and joint twinning-and transformation-induced plasticity in a metastable β titanium alloy[END_REF]. It could however participate in the hardening by acting as dislocation source, and still contribute to hindering the dislocation glide. Overall, the hardening is expected to be different in type-A grains and type-B grains, leading to a first grain-level type of heterogeneities. The second point is that in type B grains (favorable for TWIP), martensite can form to accommodate the strong local strain fields, at twins or at grain boundaries, that may otherwise lead to earlier failure of the sample. Therefore, heterogeneities are observed at two different levels. The first one is the level of individual grains or groups of grains, and the second one is the intragranular scale. Thus, the distinct yet complementary roles of the two main deformation mechanisms TRIP and TWIP in type A and B grains may be used to explain the outstanding ductility and UTS in this alloy. It is well known that both mechanisms, taken individually or together, allow the beta structures to reach better mechanical properties, compared to simple dislocation glide [START_REF] Min | Mechanism of twinning-induced plasticity in β-type Ti-15Mo alloy[END_REF][START_REF] Castany | In situ synchrotron X-ray diffraction study of the martensitic transformation in superelastic Ti-24Nb-0.5N and Ti-24Nb-0.5O alloys[END_REF][START_REF] Min | Deformation microstructural evolution and strain hardening of differently oriented grains in twinning-induced plasticity β titanium alloy[END_REF][START_REF] Zhou | Accommodative 332〈113〉 primary and secondary twinning in a slightly deformed β-type Ti-Mo titanium alloy[END_REF][START_REF] Hanada | Effect of Plastic Deformation Modes on Tensile Properties of Beta Titanium Alloys[END_REF][START_REF] Yang | Characterization of the martensitic transformation in the superelastic Ti-24Nb-4Zr-8Sn alloy by in situ synchrotron X-ray diffraction and dynamic mechanical analysis[END_REF]. Two main effects are used to explain this

M A N U S C R I P T A C C E P T E D

ACCEPTED MANUSCRIPT 20 improvement: the dynamic Hall-Petch effect and the composite effect, which relates to the co-deformation of the various phases. This second effect seems of particular interest in this case, considering the two levels of heterogeneities discussed before. There is therefore a two-level composite effect observed in this alloy. The first one is an intergranular mechanical contrast between harder grains and softer grains. The second one is the intragranular mechanical contrast, seen between the matrix and the deformation products -harder for twins, softer for martensite. This architectural strain distribution may explain the excellent work-hardening capability and ductility of the alloy.

Conclusion

In-situ studies on the heterogeneous nature of deformation in a strain-transformable Ti-V-Cr-Al alloy are reported. The multi-scale analysis of the deformation mechanisms occurring in this alloy allowed to confirm the deformation mechanisms TRIP and TWIP, and the focus is put on the heterogeneity of the different zones of the sample, which depends highly on the grain orientation with respect to the tensile axis. This heterogeneous microstructure results in complementary co-deformation of the material at two different scales (intragranular and intergranular). Finally, with the help of in-situ techniques, the behavior of the martensite could be accurately investigated. An important result of this study is that, based on their orientation, grains seem to prefer one or the other deformation mechanism between TRIP and TWIP. Two different behaviors, twinning and martensitic transformation, are thus observed in the primary

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

21 stages of the deformation, each of them for a different group of grains. This coincides well with a fast increase of the martensite fraction, found in the form of needles in the grains favorably oriented, that reaches about 30% for ε≈0.12. In the latter stage of the deformation, the overall martensitic content stabilizes, even though it is still formed in small amounts in TWIP-favored grains as secondary mechanism, to accommodate the highly strained areas. These results suggest also that there is, for this alloy, no unique deformation path, but rather two parallel scenarios. It is suggested that the exceptional mechanical properties of this alloy, a yield strength of 420 MPa, a UTS of almost 1200 MPa and a uniform elongation of 35%, can be analyzed as the result of the multi-scale heterogeneities and co-deformation of the various deformation features of this material.

  the last EBSD scan.

  M A N U S C R I P T A C C E P T E DACCEPTED MANUSCRIPT 10 labelled T1 to T4 for the β-twins, and a1 to a3 for the martensitic variants.

1 Figure 1 : 2 Figure 2 :

 1122 Figure 1: Uniaxial mechanical properties of Ti-V-Cr-Al alloy in the ST state. (a) true stress-true strain curve. (b) Work hardening rate as a function of the true stress.

3 Figure 3 : 4 Figure 4 : 5 Figure 5 : 6 Figure 6 : 7 Figure 7 :

 3344556677 Figure 3: EBSD mapping of the in-situ traction specimen before deformation. (a) Inversed pole figure, IPF. (b) Image quality, IQ. (c) Schmid factor (0.4<SF max <0.5) distribution of SIM a" along tensile direction. (d) Schmid factor (0.4<SF max <0.5) distribution of {332}<113> twinning system along tensile direction. Examples of groups of Type A grains and Type B are marked by yellow dash lines in (c) and (d).

Acknowledgements:

This work was supported by the French Agence Nationale de la Recherche, on the ANR TiTWIP. The authors acknowledge the European Synchrotron Radiation Facility for provision of synchrotron radiation facilities and would like to thank Andy Fitch for assistance in using the beam line ID22.


M A Hadj-Youcef
email: amine.hadjyoucef@ias.u-psud.fr

F Orieux

A Fraysse

A Abergel

Spatio-Spectral Multichannel Reconstruction from few Low-Resolution Multispectral Data

Keywords: Inverse problems, Image reconstruction, Deconvolution, System modeling, Multispectral restoration

This paper deals with the reconstruction of a 3-D spatio-spectral object observed by a multispectral imaging system, where the original object is blurred with a spectralvariant PSF (Point Spread Function) and integrated over few broad spectral bands. In order to tackle this ill-posed problem, we propose a multispectral forward model that accounts for direct (or auto) channels and between (or cross) channels degradation, by modeling the imaging system response and the spectral distribution of the object with a piecewise linear function. Reconstruction based on regularization methods is proposed, by enforcing spatial and spectral smoothness of the object. We test our approach on simulated data of the Mid-InfraRed Instrument (MIRI) Imager of the James Webb Space Telescope (JWST). Results on simulated multispectral data show a significant improvement over the conventional multichannel method.

I. INTRODUCTION

Multispectral imaging systems are used in many fields, e.g. astrophysics [START_REF] Lionello | Multispectral emission of the sun during the first whole sun month: Magnetohydrodynamic simulations[END_REF], remote sensing [START_REF] David A Landgrebe | Signal theory methods in multispectral remote sensing[END_REF], medicine [START_REF] Troy O Mcbride | Multispectral near-infrared tomography: a case study in compensating for water and lipid content in hemoglobin imaging of the breast[END_REF] or microscopy [START_REF] Me Dickinson | Multispectral imaging and linear unmixing add a whole new dimension to laser scanning fluorescence microscopy[END_REF]. This paper deals with the inverse problem of joint restoration. Our goal is to reconstruct a discrete 3-D spatio-spectral object from a small number of 2-D Multi-Spectral (MS) observed data when this continuous 3-D object is degraded by the instrument that suffers the diffraction due to the limited size of the optical system. This physical degradation affects its spatial resolution (in the form of blur) accordingly to the wavelength. Moreover, before its spatial sampling, the blurred object is integrated by the detector over the different wide spectral bands, which results in low spectral resolution multispectral data. Therefore, the multispectral data are severely degraded and contain limited spectral information about the original object.

Multichannel restoration has been extensively studied in the literature. Multichannel forward models have been proposed in [START_REF] Hunt | Karhunen-loeve multispectral image restoration, part i: Theory[END_REF], [START_REF] Nikolas | Digital restoration of multichannel images[END_REF], where the system response is a block-diagonal matrix with circulant blocks. For instance, [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF], [START_REF] Bongard | 3d deconvolution of hyper-spectral astronomical data[END_REF], [START_REF] Song | Regularization parameter estimation for non-negative hyperspectral image deconvolution[END_REF] address multichannel 2-D deconvolution problem for hyperspectral image deconvolution. They take into account the within-channel degradation, but not the between channel (or cross-channel) degradation. Hence, this approach is not suitable for MS imaging, especially if spectral bands are broad and overlapping, which implies a strong correlation between channel. In [START_REF] Nikolas P Galatsanos | Least squares restoration of multichannel images[END_REF] a model is proposed that reduces these limitations since the system response is represented by a block matrix corresponding to within and between channel degradation. However, this model is mostly used when the number of channels and observations is the same, e.g. color image restoration [START_REF] Richard | Stochastic modeling and estimation of multispectral image data[END_REF], [START_REF] Wen | Efficient total variation minimization methods for color image restoration[END_REF], [START_REF] Yang | A fast algorithm for edge-preserving variational multichannel image restoration[END_REF].

In this paper we propose a multispectral forward model that accounts for within-and between-spectral channel degradations (or auto and cross-channel), where (1) the number of MS data is much lower than the number of spectral channels and (2) a set of low-resolution MS data are degraded by a spectral-variant PSF and integrated over broad spectral bands. Reconstruction of a spatio-spectral object is performed using regularization, by accounting for spatial and spectral quadratic regularization. Simulated results are provided with a comparison to multichannel 2-D deconvolution for an application to the MIRI Imager on board the JWST 1 .

This paper is organized as follows. In Sec. II we present the problem formulation. The imaging system response and the forward model are described in Sec. III. The reconstruction method is presented in Sec. IV. Simulation and results are presented in Sec. V including a brief description of the JWST/MIRI Imager. Conclusions and perspectives are in Sec. VI.

II. PROBLEM FORMULATION

The general form of the multispectral problem we are considering is the one proposed in [START_REF] Nikolas P Galatsanos | Least squares restoration of multichannel images[END_REF], [START_REF] Richard | Stochastic modeling and estimation of multispectral image data[END_REF]. It follows the following discrete linear forward model :

y = Hx + n, (1)
where x = x (1)T , x (2)T , . . . , x (M )T T is the stack of M spectral channels represented in a vector form, each channel contains N × N pixels, where x (m) i,j denotes the (i, j) th spatial position in the (m) th wavelength. The vector y = y (1)T , y (2)T , . . . , y (P )T T ∈ R P N 2 is the stack of multispectral observed data acquired via P broad spectral bands of the imaging system. n = n (1)T , n (2)T , . . . , n (P )T T ∈ R P N 2 represents an additive unknown noise. The full system response is a P N 2 by M N 2 block matrix,

H =    H 1,1 H 1,2 • • • H 1,M . . . . . . . . . . . . H P,1 H P,2 • • • H P,M    , P M, (2)
defined by a set of P × M Toeplitz sub-matrices H p,m ∈ R N 2 ×N 2 , which are approximated for computational ease by circulant blocks. All sub-matrices are ill-conditioned, meaning that H is also ill-conditioned, which leads to an ill-posed problem. Each multispectral observation y (p) depends on all spectral channels through the blocks of H, where H p,t , t = m represents the direct (or auto) observation, whereas the block H p,t , t = m accounts for between (or cross) degradation occurring between channels. We are interested in the particular case where we have few low-resolution multispectral observed data compared to spectral channels, i.e. P M , which means there is a lack of spectral information in the data. For instance, multispectral data observed by the JWST/MIRI imager, e.g. P = 9 and M = 1000, and the size of H is 9 × 256 2 by 1000 × 256 2 , for a 256 × 256 pixel detector.

III. MODELING OF WITHIN AND BETWEEN CHANNELS DEGRADATION

A. Object Model

We first define the 3-D spatio-spectral object of interest with φ(α, β, λ) : R 3 → R, having two spatial parameters (α, β) ∈ R 2 and one spectral parameter λ ∈ R + . In order to handle the lack of spectral information in the data, we propose to perform approximation of the M spectral channels by M channels, with M < M , by modeling the object spectral distribution with a piecewise linear function [START_REF] Ma Hadj-Youcef | Restoration from multispectral blurred data with non-stationary instrument response[END_REF]. Moreover, we are interested in reconstructing a discrete version of the object, hence, we define two basis functions, for spatial and spectral discretizations, b s and b λ , respectively. They are defined upon two grids,

G s = {α k , β l } N k,l=1 and G λ = λ (m) M m=1 , respectively. The object is modeled by φ(α, β, λ) = M m=1 N k,l=1 x (m) k,l b (k,l) s (α, β)b (m) λ (λ), (3) where b (m)
λ (λ) is a uniform piecewise linear function, for instance first-order B-spline function [START_REF] Thévenaz | Interpolation revisited [medical images application[END_REF]. The parameter M compromises between the sharpness of the spectral sampling of the modeled object and the unknown spectral channels to reconstruct x (m) .

B. Imaging System Response

In this section we provide the multispectral imaging system response by establishing an imaging system model that relates the input to the output. The instrument we are considering is composed of an optical system and a detector. Due to light diffraction of φ, the optical system response is modeled by a 2-D spatial convolution [START_REF] Joseph W Goodman | Introduction to Fourier optics[END_REF] with a spectral variant optical response, known as Point Spread Function (PSF) h(α, β, λ).

This blurs the object accordingly to the wavelength and limits its spatial resolution (as illustrated in Sec. V). The blurred object is integrated over broad spectral bands τ p (λ) and sampled pixel-by-pixel on the 2-D detector grid, G samp = {α i , β j } N i,j=1 : α i , β j being the 2-D angular positions of pixels (i, j) and N is the total number of pixel according dimensions α and β. We introduce a basis function b (i,j) samp (α, β) to carry out spatial sampling. It is defined on the pixel sensitive surface Ω pix . Moreover, a noise term n (p) i,j is added for each pixel (i, j) and band p, e.g. readout noise of the detector. Finally, the imaging system model is given by

y (p) i,j = R+ τ p (λ) Ωpix R 2 φ(α , β , λ) (4)
h(α -α , β -β , λ)dα dβ b (i,j) samp (α, β)dαdβ dλ + n (p) i,j
this model links the 3-D continuous input φ(α, β, λ) to the 2-D discrete output y (p) i,j through a complex instrument response, which includes spectral windowing and five sums, two for spatial 2-D convolutions, two for spatial sums and one for spectral integration. Note that the above model does not include any non-ideal characteristics of the detector, which are assumed to be corrected upstream.

C. Forward Model and definition of H p,m

The discrete forward model links the discrete spectral channels to the discrete multispectral data. It is obtained by substituting equation ( 3) in ( 4). This yields

y (p) i,j = M m=1 N k,l=1 H p,m i,j;k,l x (m) k,l + n (p) i,j , (5)
with

H p,m i,j;k,l = R+ τ p (λ)b (m) λ (λ) Ωpix R 2 b (k,l) s (α , β ) h(α -α , β -β , λ)dα dβ b (i,j) samp (α, β)dαdβ dλ. (6)
In addition, we consider for instance a rectangular impulse function [START_REF] Yaroslavsky | Digital holography and digital image processing: principles, methods, algorithms[END_REF] for the sampling function. i.e. b

(i,j) samp (α, β) = 1 α β Π α, β (α-α i , β-β j )
, with α, β are the sampling steps according to dimensions α and β, respectively. Thus, the system response becomes a convolution matrix H p,m i,j;k,l = H p,m i-k;j-l . Therefore, the vector-matrix representation of ( 5) is

y (p) = M m=1 H p,m x (m) + n (p) , (7)
where the p-th multispectral data y (p) is a sum of M discrete 2-D spatial convolutions between spectral channels and convolution matrices H p,m (blocks of the matrix H in (2)). Thus, it accounts for within and between channels degradation. The discrete multispectral forward model with the full imaging system H response takes the form in [START_REF] Lionello | Multispectral emission of the sun during the first whole sun month: Magnetohydrodynamic simulations[END_REF].

IV. RECONSTRUCTION

Reconstruction of the object of interest φ relies on the reconstruction of its spectral channels x (m) : m = 1, 2, . . . , M using a regularization method. The solution of reconstruction is obtained by minimizing an objective function J (x), x = argmin

x (1) ,...,x (M ) {J (x)} , (8)
this method is called regularized least squares [START_REF] Nikolaevitch Tikhonov | Numerical methods for the solution of ill-posed problems[END_REF], where we add regularization terms to J (x) in order to correct the illposedness of the problem. The objective function is defined by

J (x) = Q(x, y) + R s (x) + R λ (x). (9)
where the first term Q(x, y) = y -Hx

D λ ∈ R M N 2 ×M N 2
is a firstorder finite difference operator along the spectral direction. µ s and µ λ are regularization parameters, they are set to compromise between fidelity to the data and spatial smoothness, and spectral smoothness across channel, respectively.

By rewriting the objective function, (9) becomes

J (x) = y -Hx 2 2 + µ s D s x 2 2 + µ λ D λ x 2 2 .
(10) It is the sum of quadratic terms only, therefore, an explicit solution is available by canceling the gradient of J (x). This yields Qx = H T y, [START_REF] Richard | Stochastic modeling and estimation of multispectral image data[END_REF] with

Q = (H T H + D T s D s + D T λ D λ ).
It contains Toepltiz blocks Q (i,j) : i, j = 1, . . . , M of size N 2 × N 2 . However, Q (i,j) = Q (i+t,j+t) , hence Q is not a Toeplitz matrix. We propose to compute the solution without inverting Q, but by computing the solution iteratively using the following form:

xk+1 = xk -a Qx k -H T y , (12)
with xk=0 = 0 corresponds to the initialization and a is a convergence parameter of the algorithm. A conjugated gradient (CG) algorithm [START_REF] Richard | An introduction to the conjugate gradient method without the agonizing pain[END_REF] is implemented.

V. SIMULATION RESULTS

A. JWST/MIRI Imager

We apply the proposed reconstruction algorithm to multispectral data simulated using the model in (4) of the Mid-InfraRed Instrument (MIRI) Imager [START_REF] Bouchet | The Mid-Infrared Instrument for the James Webb Space Telescope, III: MIRIM, The MIRI Imager[END_REF] on-board the James Webb Space Telescope (JWST), the next flagship space telescope of NASA, ESA and the Canadian Space Agency (CSA) to be launched in 2020. The imager provides nine multispectral observations (P = 9) integrated over a broad range of spectral bands, from 5 µm to 28 µm [START_REF] Glasse | The Mid-Infrared Instrument for the James Webb Space Telescope, IX: Predicted Sensitivity[END_REF]. The nine bands are shown in Fig. 1. Note that overlapping of the spectral bands increase the between channels degradation. The MIRI Imager detector has a pixel pitch of 0.11 arcsecond, i.e. Ω pix = 0.11 × 0.11 arcsecond 2 . We use the official PSF simulator of the JWST mission, WebbPSF [START_REF] Marshall D Perrin | Simulating point spread functions for the James Webb Space Telescope with WebbPSF[END_REF], [START_REF] Marshall D Perrin | Updated point spread function simulations for jwst with webbpsf[END_REF], to simulate a realistic PSF image at different wavelengths, as shown in Fig. 2. The PSF is complex due to the segmented mirror of the JWST. We clearly observe an enlargement of the PSF according to the wavelength, i.e. the longer the wavelength the wider the PSF, as expected from diffraction theory [START_REF] Joseph W Goodman | Introduction to Fourier optics[END_REF]. Fig. 2. Monochromatic PSF of the JWST/MIRI imager simulated using WebbPSF [START_REF] Marshall D Perrin | Simulating point spread functions for the James Webb Space Telescope with WebbPSF[END_REF], [START_REF] Marshall D Perrin | Updated point spread function simulations for jwst with webbpsf[END_REF] and displayed in the same logarithmic scale. We clearly observe the dependency of the PSF to the wavelength as expected from diffraction theory.

B. Setup of the Experiment

The original 3D object φ is a simplified spatio-spectral model of the Horse Head nebula [START_REF] Abergel | Isocam and molecular observations of the edge of the horsehead nebula[END_REF]. A spatial region of 256 × 256 pixels (N = 256) is taken for the simulation with M = 1000 spectral samples uniformly distributed within 4-28 µm. Nine multispectral data are simulated using (4) with a zero-mean white Gaussian noise added in order to obtain a global Signal-to-Noise Ratio (SNR) of 30, 20, 10 dB. SN R (dB) = 10 log 10 1

P N 2 y 2 2 σ 2 n , (13)
where σ n is the standard deviation of the noise. P is the number of MS data and N 2 is the total number of pixels in the MS data.

Reconstruction results are summarized in Table I together with a comparison between the proposed method and the multichannel 2-D deconvolution method (MDec) [START_REF] Nikolas | Digital restoration of multichannel images[END_REF] (independent channel restoration using an averaged PSF per channel). The regularization parameters µ s and µ λ are adjusted by running the code for different values in a range [10 -4 , 10 -2 ] and keeping the pair that minimizes the objective function J (x( μs , μλ )). For a quantitative comparison between the original object f orig and the reconstructed f rec , we compute the relative reconstruction error as defined by Error(%) = 100 × f orig -f rec 2 / f orig 2 .

C. Discussion

Fig. 3 shows the spectral reconstruction result of one single pixel (100, 150), comparing the original spectrum φ orig , the reconstructed spectrum using our method φ rec , and the reconstructed spectrum using multichannel 2-D deconvolution φ M Dec . The original spectral distribution is complex with spectral features at short wavelengths (4 -13 µm) and continuum. Therefore, an accurate reconstruction using a few MS data is difficult, if not impossible, without using a strong prior knowledge of the spectrum of the object. The reconstructed spectrum computed with our method (using piecewise linear model) allows us to reconstruct an envelope-like spectral distribution which significantly increases the spectral resolution compared to multichannel 2-D deconvolution. Several values of M = {20, 40, 60} have been tested, and the reconstruction results for three wavelengths, 7.8, 16 and 21 µm, are reported in Table I. Increasing M improves the spectral resolution of the object model, but increases the between channels degradation and the number of unknowns. Moreover, we find that there is not much error improvement for M > 60. In any case, the proposed reconstruction shows smaller reconstruction errors compared to the multichannel 2-D deconvolution ; this is due to our model accounting for within and between channels degradation. Spatial reconstruction results at different wavelengths are illustrated in Fig. 4. As anticipated, a better reconstruction is obtained at λ = 16 µm and λ = 21 µm than at λ = 7.8 µm (see the fourth row of the figure) since within the integration windows at long wavelengths the spectrum of the object does not contain any feature.

VI. CONCLUSION

In this paper we address the reconstruction of 3-D spatiospectral object observed by a multispectral imaging system from a few low-resolution data. A discrete forward model is defined accounting for within and between channels degradation using a piecewise linear function to model the spectral distribution of the sought object. A quadratic reconstruction is proposed by considering spatial and spectral regularization terms. Results on simulated data applied to the JWST/MIRI Imager highlights the complexity of the problem. A clear increase of spatial and spectral distribution is achieved compared to multichannel 2-D deconvolution method.

2 2 2 2 2 .

 222 is the data fidelity. It enforces agreement of the solution with the data. The second term R s (x) = µ s D s x is a spatial regularization. It enforces spatial smoothness between neighboring pixels of x, with D s ∈ R M N 2 ×M N 2 is a second-order finite difference operator along the spatial dimension. The last term refers to spectral regularization R λ (x) = µ λ D λ x 2 It enforces the similarity between intensity values of corresponding pixels in neighboring channels, with

Fig. 1 .

 1 Fig. 1. The nine broad bands of the JWST/MIRI Imager [21] covering the mid-infrared wavelength range from 5 to 28 µm.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Comparison between one single pixel (100, 150) spectrum from the original object φ orig , the proposed reconstruction φrec (with M = 60) and the multichannel 2-D deconvolution φ M Dec . The nine MS data (P = 9) were corrupted with zero-mean Gaussian noise of 30 dB.

TABLE I PROPOSITION

 I AND COMPARISON OF HORSE HEAD NEBULA[START_REF] Abergel | Isocam and molecular observations of the edge of the horsehead nebula[END_REF] OF SIZE 1000 × 256 × 256

	SNR	λ		Error (%)
	(dB)	(µm)	Proposed Reconstruction	MDec
			M'=20	M'=40 M'=60
		7,8 49,44	42,37	41,42	52,85
	30	16,0	2,44	4,11	4,80	7,89
		21,0	1,87	3,82	4,26	11,92
		7,8 49,50	43,07	41,46	52,84
	20	16,0	7,41	7,66	8,98	8,02
		21,0	4,42	5,40	5,77	11,97
		7,8 50,71	43,71	42,38	52,84
	10	16,0 19,67	21,25	25,51	8,56
		21,0 10,85	11,31	13,38	12,13

https://jwst.nasa.gov/