J G Mattei

F Pelletier

Diana Ciuculescu-Pradines

Pierre Lecante

Jean-Charles C Dupin

N Yaacoub

Joachim Allouche

Jean-Marc M Greneche

Danielle Gonbeau

Catherine Amiens

D Ciuculescu

¶ N Yaacoub

¶ C Amiens

M J Casanove
email: casanove@cemes.fr

† Cemes

Formation of bimetallic FeBi nanostructured particles: Investigation of a complex growth mechanism

Keywords: Bimetallic nanoparticles, bismuth, iron, growth mechanisms, structure, soft chemistry

come   L'archive ouverte pluridisciplinaire

Introduction

Nanosciences and related technologies have recently focused an important research effort on the fabrication and study of multicomponent and hybrid particles. 1 The presence of discrete domains of different materials in a same particle is indeed highly desirable for many potential applications: on the one hand, it offers the remarkable advantage of combining different properties in a single object, on the other hand it opens a route to new properties and novel functions. 2-4 Among multicomponent particles, the ones with core-shell geometry are highly valuable because the core material is effectively protected from the environment by the shell, while the shell material can provide enhanced properties, either optical 5 or magnetic, 6 to the core-shell system. A great attention has been paid to the combination of iron with a non ferromagnetic metal in a single particle.

In these materials, the use of a metal presenting a large miscibility gap with iron is particularly appealing to create core-shell chemical distribution. Several experimental or theoretical studies have thus been dedicated to the Fe/Au system in which Fe is expected to remain at the core of the particle. [7][8][9][10] In the present paper, we focus on the iron-bismuth system in which both metals display very different properties and present a large miscibility gap. While α-Fe is ferromagnetic, bismuth behaves as a semi-metal with unique electronic properties. In particular, remarkable transport and thermoelectric properties have been reported in bismuth nanoscaled materials or nanocomposites. [11][12][13][14] Very few studies have been dedicated to the immiscible Fe-Bi system and all of them concerned thin films, either Fe/Bi multilayered films 15 or Fe-Bi bimetallic films grown by cosputtering of elemental targets or dual electron-beam evaporation or even through ion mixing in previously grown Fe/Bi multilayered films. [16][17][18] In a preliminary study, we developped a successful synthesis of FeBi bimetallic nanocomposites, although the assembly of particles was not fully homogeneous. 19 In the present paper, we developped further the synthesis in the presence of hexamethyldisilazane (HMDS), which allowed the formation of bimetallic nanocomposites, so as to control the formation of an assembly of identical FeBi particles, required to perform a thorough structural and chemical analysis of the synthesized particles. Typically, we used a Fe/Bi atomic ratio as small as 7/1 instead of 1/1 as in, 19 with the aim to avoid the nucleation and growth of pure bismuth nanocrystals. Excess HMDS (namely 10 equivalent versus metal content) was used to compensate for the high lability of this amine and ensure stabilization of the colloidal solution. As our samples grow in mild conditions of temperature and pressure, and are surrounded by ligands, their structure and morphology can be very different from those expected from thermodynamics. On the contrary, kinetics can play a rather important part, as we observed in the Fe-Rh system. 20 Thorough examination of the structural and chemical features, from the core to the surface of the particles, is thus mandatory if one aims at explaining the properties of such complex compounds. For these reasons, we combined in this study wide-angle x-ray scattering (WAXS), transmission electron microscopy (TEM) with related analytical techniques such as energy dispersive spectroscopy (EDS), extended x-ray absorption fine structure (EXAFS), x-ray photoelectron spectroscopy (XPS) and transmission 57 Fe M össbauer spectrometry. The magnetic properties of the Fe-Bi nanospheres were also measured using a superconducting quantum interference device (SQUID).

Experimental section

Synthesis

General methods

All operations (material preparation, sampling and packaging) were carried out using standard Fischer-Porter bottle techniques and glove-box under argon. Fe(N(SiMe 3 ) 2 ) 2 (>99%) was purchased from Nanomeps and anhydrous anisole from Sigma-Aldrich. Anisole was purified by distillation under argon (over sodium) and further degassed by the freeze-pump-thaw technique and finally dried over molecular sieve (0.4 nm) before use. Bi(N(SiMe 3 ) 2 ) 3 21 and iPr 2 NH -BH 3 22

were prepared and purified according to published procedures. Dihydrogen was purchased from Air Liquide and contained less than 3 ppm of H 2 O and 2 ppm of O 2 . To determine the Fe and Bi contents in the materials, elemental analyses were performed at Antellis 23 by ICP-MS (Inductively Coupled Plasma-Mass Spectrometer), after digestion of the samples in concentrated HNO 3 .

Sample preparation

Fe(N(SiMe 3 ) 2 ) 2 (80 mg, 0.21 mmol), Bi(N(SiMe 3 ) 2 ) 3 (20mg, 0.03 mmol), HMDS (343 mg, 2.1 mmol) were dissolved in anisole (18 mL, H 2 O < 3ppm) in a Fischer-Porter bottle in a glove box.

The homogeneous yellow-green solution was freezed under liquid nitrogen and 1.1 equivalent of diisopropylamine-borane with respect to iron (27 mg, 0.24 mmol) in anisole (4 mL) were canula-transferred in the Fischer-Porter bottle. Slow warm-up of the reaction medium up to room temperature afforded a dark brown solution. This homogeneous solution was then pressurized under 3 bar H 2 and heated 4h at 110

WAXS and EXAFS experiments

The samples for WAXS experiments consisted in small amounts of the pure powders sealed in 1 mm Lindemann glass capillaries. Measurements were performed using a dedicated two-axis diffractometer equipped with a graphite flat monochromator in order to select the molybdenum K α wavelength (0.07107 nm

TEM and EDS experiments

Samples for TEM experiments were dispersed in toluene or anisole previous to deposition on ultrathin carbon films supported by a lacey carbon film on a 400 mesh copper grid. They were exposed to air just before their introduction in the microscope. The samples were analyzed using a Cs-corrected Tecnai F20 and a CM20 FEI microscopes. Both microscopes are fitted with a field emission gun. Elemental mapping of isolated particles was performed by energy dispersive xray spectroscopy (EDS) to address compositional variations from the core to the surface of the particles. The EDS spectra were recorded in STEM mode on the CM20, which is fitted with a Bruker QUANTAX EDS system.

XPS experiments

Samples for XPS experiments were prepared in an ultra-pure atmosphere (within the argon glovebox directly connected to the spectrometer) and mounted on a stainless sample holder using conductive copper tape. Surface elemental and chemical state analyses were performed on a Kα

x-ray photoelectron spectrometer (Thermo Scientific). This instrument is equipped with a microfocusing monochromator (Al-Kα x-ray source at 1486.6 eV), which was operated at a spot size of 400 µm diameter, and a pass energy of 20 eV. The spectrometer was calibrated using the photoemission lines of Au (Au4 f 7/2 = 83.9 eV, with reference to the Fermi level) and Cu (Cu2p 3/2 = 932.5 eV); for the Au4 f 7/2 line, the full width at half maximum (FWHM) was 0.86 eV in the recording conditions. Charge effects were compensated by the use of a dual-beam charge neutralisation system (low energy electrons [typically 1.85 eV] and Ar + ions) which had the unique ability to provide consistent charge compensation. All the neutraliser parameters remained constant during analysis. For the study, charge referencing was done against adventitious carbon (C1s binding energy = 284.6 eV). Spectra were analyzed and mathematically described using CasaXPS software. The XPS signals were analysed by using a least squares algorithm and a non-linear baseline. The fitting peaks of the experimental curves were defined by a combination of gaussian (70 %) and lorentzian (30 %) distributions. [START_REF]Different standard compounds were analyzed as references in order to precise some attributions: B 2 O 3 material from Aldrich (CAS number 1303-86-2)[END_REF]

M össbauer experiments

The samples for M össbauer analyses consist of a thin layer containing about 5 mg of Fe per cm 2 of FeBi powders dispersed in PPO that was sealed in glove box with inert atmosphere to avoid oxidation. 57 Fe M össbauer spectra were recorded using a 57 Co/Rh γ-ray source mounted on a conventional electromagnetic drive using a triangular velocity form. They were obtained at 77K, 300 K without external magnetic field and at 10 K in an 8 T field oriented parallel to the γ-beam.

The M össbauer spectra were analysed using the 'MOSFIT'. [START_REF] Teillet | Mosfit software[END_REF] The hyperfine structure was modelled by a least-square fitting procedure involving Zeeman sextets composed of either lorentzian or gaussian lines. The isomer shift values (IS) were referred to that of the α-Fe at 300 K.

SQUID experiments

Hysteresis loops were recorded at 2 K, between -5 T and +5 T using a MPMS Quantum Design SQUID. The as-prepared powder was sealed in a gelatin capsule to prevent any oxidation during the transfer of the sample into the magnetometer. Elemental analysis was carried out directly on the capsule to enable normalisation of the magnetization with respect to the Fe content in the sample.

Results and discussion

A core-shell chemical distribution

Structure and morphology

TEM observations give evidence for the presence of only one type of particles: nanospheres with a mean diameter 150 ± 30nm, as determined from the analysis of about one hundred particles. A typical TEM image is displayed in Figure 1. The enlarged view of a nanosphere shows that it is composed of numerous small grains about 2 to 3 nm large, at least at the surface region thin enough to present electron transparency. The core-shell contrast observed in such large particles cannot be readily attributed to a difference in chemical composition as it is also the expected contrast for a spherical object, the thicker region of the core being much more absorbent for the incident electron beam than the surface region.

Figure 1: TEM image of the Fe-Bi nanoparticles. An enlarged view of one particle, revealing the presence of numerous nanograins, is displayed in inset.

In order to check the compositional variation across the nanospheres and in particular from core to surface, several nanospheres were examined using spatially resolved STEM-EDS with an electron probe size of 2 nm. Spectrum images were recorded using a drift correction and further processed for quantitative analysis of the different elements. Figure 2 presents the elemental mapping data for an ensemble of three particles, which provides clear evidence for a compositional variation from shell to core: the shell seems to be mainly composed by iron while most of the bismuth content seems concentrated in the particle core. The quantitative analysis of the displayed group of particles gives a Fe/Bi ratio of Fe : Bi = 0.875/0.125 = 7, remarkably close to the nominal ratio introduced at the beginning of the synthesis. Analyses performed on spectrum images of different particles gave very close results, the ratio remaining in the range 7 ± 0.25.

The compositional gradient from surface to core of the nanospheres is better evidenced when the atomic composition is plotted along a line profile taken across a nanosphere diameter. Such a plot is displayed as supporting information (figure SI1). Clearly, the ratio of the different elements strongly differ in the surface and in the core of the nanosphere. Of course, while the displayed surface region is effectively representative of the shell, the spectra recorded in the central part of the nanospheres include both the contribution from the particles core and the contribution from top and bottom regions of the shell. The shell composition, analyzed from small regions taken in the surface region, is thus estimated to Fe/Bi = 93.6/6.4 = 14.625. The composition of the core was then deduced from the two known compositions, i.e. mean composition and shell composition, using a simple model of a spherical particle fully occupied by Fe and Bi atoms with different atomic volumes. Such estimation gives a Fe/Bi ratio slightly less than 1. One can conclude from these analyses that the observed nanospheres are fully representative of the whole sample as they seem to gather all the metallic atoms introduced in the synthesis. The core of the nanospheres, with a volume about 5 times smaller than the shell, has a composition close the Fe/Bi = 1/1 (i.e; around 50% of Bi) while the shell includes only 6.4% of Bi. However, nothing can be said at this step about the distribution of the different metals in these regions, and in particuler about the possible formation of metastable Fe-Bi alloys at this nanometer scale. The observed core-shell distribution is not the one expected from both thermodynamic and kinetic considerations. Bismuth has a much lower surface energy than iron (0.537J.m -2 for Bi compared to 2.222J.m -2 for Fe) 27 as well as a much larger atomic size. These two parameters should promote a segregation of bismuth at the surface of the nanospheres owing to energy minimization considerations. Moreover, the synthesis procedure itself was designed to ensure the sequential decomposition of the two different precursors. Indeed, control experiments (detailed in SI) revealed that the bismuth precursor is very stable and requires a temperature as high as 150 • C to be decomposed under dihydrogen, while the reduction of the iron precursor in the presence of amine-borane is quasi-instantaneous at the very beginning of the process carried out close to -40 • C. In order to shed some light on this unexpected chemical distribution, the close atomic environment of the two metals was investigated by different complementary techniques.

A chemical overview of the Fe-Bi compound was first given through XPS analyses. In addition to iron and bismuth, which were both found in a metallic state, the XPS analysis highlights the presence of some traces of light elements coming from metal-organic precursors/HMDS decomposition as well as boron coming from the amine-borane used as reducing agent and easily dehydrogenated. The presence of B-Fe bonds was clearly revealed at the B-1s core peak (presented in supporting information, figure SI2). Indeed this peak displayed a small component at 188.2 eV associated with a low atomic percentage of boron (1.2%). This component corresponds to a specific fingerprint of Fe-B bonds recorded from a reference FeB alloy in agreement with a previous work. 28 In light of these results, one can conclude that boron, coming from the amineborane reducing agent, is present in the iron close environment. [START_REF] Glavee | Boron inclusion in Fe nanoparticles is well attested for when borohydride reducing agents are used (see for example[END_REF] X-ray absorption experiments were also carried out at the Fe K (7112eV) and Bi L 3 (13419eV) edges. Figure 3 presents the uncorrected Fourier transforms (FT) of the EXAFS functions at both edges together with the references from pure Fe and Bi foils. The FT at Fe edge is consistent with metallic iron, however it significantly differs from the bcc pattern observed for the reference foil:

the peak including bonding distances is broad and shifted to shorter distances. Its amplitude is reduced tenfold compared to bulk iron. This points towards an environment for Fe atoms characterized by a strong dispersion, both in bond lengths and chemical species (i.e. including Fe atoms but also lighter elements like C, N, O and B), however already observed for small iron particles also synthesized by a soft chemical route. 30 Such a disorder is not observed at the Bi edge, but the pattern clearly evidences a mixed environment for Bi atoms. Even if FTs of heavy elements like Bi must be considered carefully (e.g., the broad asymetrical peak for the Bi foil is actually related to a single well defined bonding distance), there is no real ambiguity : the longer distance (0.3 nm, uncorrected) is related to a metallic Bi environment, but the shorter one (0. For easier comparison, the reference foils were adjusted using a scaling factor (0.10 and 0.80 for Fe and Bi, respectively).

The local environment of Fe was further analyzed through M össbauer spectroscopy experiments. The zero field 57 Fe M össbauer spectrum collected at 77K displays a pure magnetic sextet composed of asymmetrical broadened lines (figure SI3 in Supporting Information). The broadening of the lines probably results from the chemical inhomogeneity. Different fitting models were attempted to describe the hyperfine structure, assuming either one or two magnetic components with distribution of hyperfine fields and lorentzian or gaussian lines. The best agreement, in par-ticular for the fit of the inner lines, was obtained when considering 2 magnetic components with gaussian lines and independent linewidths. However the solution involving two discrete distributions of hyperfine fields was preferred as the gaussian profile is not physically appropriate.The mean hyperfine field and isomer shift values are estimated at about 26 and 33T, 0.23 and 0.18 mm/s, respectively, while the quadrupolar shift values remain close to zero (see table 1). The results agree with the occurrence of two different environments for iron (consistently with XPS).

The low and broadened hyperfine field component can be attributed to Fe atoms surrounded by B atoms: indeed, investigation of crystalline and amorphous Fe borides reveals that the hyperfine field at 57 Fe sites highly depends on the presence of boron atoms as first neighbours, with values ranging from 20 to 25T at 77K. 32 In addition, the increase of the isomer shift up to about 0.25-0.30 mm/s at 77K is also nicely consistent with the presence of B species in the close vicinity of Fe. 33 On the other hand, the high hyperfine field component is most probably due to rather pure Fe disordered phase taking account of the observed line broadening (let us recall that the hyperfine field value in the well-ordered α-Fe phase is 33.7 T at 77K). The present situation appears to be rather consistent with that observed in grain boundaries of nanocrystalline Fe materials resulting from in situ consolidated Fe clusters generated by inert gas condensation. 34 Consequently, the hyperfine structure of present Fe-Bi nanospheres should result from combined effects including the presence of boron, the disordered topological structure and the surface or interface effects which contribute also to disturb the local environment of 57 Fe probes. Such a complexity prevents from an estimation of the boron content diffused into the iron clusters. The lack of reference on the Fe-Bi immiscible system prevents further analysis on the presence of Bi in the Fe environment.

Properties

Despite the presence of iron in the shell, the nanospheres display good magnetic properties and resistance to oxidation during air exposure. Figure 4 shows the hysteresis loops recorded by SQUID experiments on both fresh (directly after synthesis) and aged (after 1 day of air exposure) Fe-Bi samples, the magnetization being normalized with respect to the iron content. The fresh sample displays an enhanced magnetization compared to bulk iron (234 A.m 2 /kg to be compared with 218

A.m 2 /kg). Such an enhancement was already observed in 2nm-wide Fe nanoparticles stabilized in a polymer displaying a poorly ordered polytetrahedral atomic arrangement, and was attributed to strong surface effects in so small particles. 35 After air exposure, the magnetization remains high, typically at 92.2% of its value in bulk iron, which represents a decrease of 14% compared to the fresh specimen. The comparison of M össbauer spectra obtained at 77 K on the as-prepared powder and after about 4 months ageing however revealed significant changes, which can be attributed to an evolution of the Fe-Bi nanospheres. The difference essentially comes from the narrowing of lines which can be assigned a priori to some local atomic reordering. Further M össbauer spectra were performed at 300K (without external field) and at 10K (under 8T external field) on aged samples to better understand the atomic scale Fe surroundings and the magnetic orientation of magnetic Fe moments (figure 5). The in-field M össbauer spectrum consists of only 4 broadened and asymmetrical lines: the disappearance of the intermediate lines indicates unambiguously that the Fe magnetic moments are aligned parallel to the γ-beam, making thus easier the hyperfine structure.

As the applied field is parallel to the γ-beam, in the cryomagnetic device, one should conclude to a ferromagnetic structure; in addition, the reduction of the magnetic splitting compared to that of the zero-field spectrum suggests that the hyperfine field is opposite to the magnetic moment. The same fitting model involving distribution of magnetic components was applied at both 300K and 77K, and at 10K under external magnetic field, in order to ensure a suitable determination of the mean value of the canting angle. The effective field (equivalent to the vectorial sum of both applied field and hyperfine field) distribution, established from the in-field spectrum, allows us to estimate the hyperfine field distribution taking into account both the canting angle and the applied field values.

The hyperfine field distributions are compared in figure 5 (right) : one observes clearly a refinement of both components, particularly the low field one, with decreasing temperature. A summary of the values of hyperfine parameters is given in table 2. The refined values of isomer shift are consistent with the presence of metallic Fe species. In addition, the in-field M össbauer spectrum fairly agrees with the presence of a ferromagnetic structure with reference to both the early observations by Hanna et al. 36 and the low canting angle. All these results suggest an evolution of the sample with time. As we cannot exclude some air exposure during the samples transfer and residual oxygen exposure during the long stay in the glove box, the observed evolution (atomic reordering, ferromagnetic metallic iron) could be attributed to the migration of some interstitial (or substitutional) boron atoms from the volume of the iron nanograins to their surface so as to form a protective boron oxide. This scenario is also supported by the observation of boron oxide in XPS analyses of aged specimens. The role of bismuth in the observed properties remains unclear.

However, enhanced resistance to oxidation and good magnetic properties were reported in Fe-Bi alloys examined by Forester and co-workers 16 and our EXAFS experiments strongly suggest the formation of an Fe-Bi amorphous alloy, at least in the core of the nanospheres. Further analyses will be necessary to gain full understanding of the different parts played by these elements but the good resistance to oxidation of these bi-metallic core@shell compounds together with the low biological toxicity of bismuth made these compounds particularly valuable for further applications in biology. 37 Figure 5: 57 Fe M össbauer spectra obtained at 300 K, 77 K (without external magnetic field) and 10 K in 8 T magnetic field applied parallel to the γ-beam on powder after 4 months (left side). Spectrum and sub-spectra result from distribution of magnetic sextets based Lorentzian lines which are represented on right side.

Origin of the formation of Bi@Fe nanospheres

In order to understand the origin of the unexpected core@shell distribution, and thus to provide a better control of the synthesized samples, we analyzed the different stages of the formation of the Fe-Bi nanospheres by collecting samples at the main steps of the synthesis sequence, as sketched in Figure 6. The different samples are numbered from S1 to S5, the previously observed Fe-Bi nanospheres corresponding to sample S4 in the sequence.

The XPS analysis of the two precursors, Fe(N(SiMe 3 ) 2 ) 2 and Bi(N(SiMe 3 ) 2 ) 3 was undertaken in order to precise the binding energy (BE) of Fe-2p and Bi-4 f core peaks in these compounds. For bismuth, besides a Bi-4 f doublet characteristic of an oxygenated environment of bismuth (BE Bi -4 f 7/2 159.0eV ) a small Bi -4 f 7/2 component appears at 158.0 eV. This BE is in line with the value reported in a recent paper for Bi-N bonds in Bi-Fe(N) layer. 38 This result shows that the signature of Bi in the precursor is clearly different from Bi metal (BE Bi -4 f 7/2 at 159.0 eV). For

Fe(N(SiMe 3 ) 2 ) 2 it was not possible to precise the signature of Fe according to the well known reactivity of this precursor with traces of water favoring the occurrence of FeO nanoparticles. 39 Figure 7 presents high resolution spectra of Bi-4 f and Fe-2p ionization peaks at different steps of the synthesis (samples S1, S2 and S4). The Bi-4 f peak displays, on all the spectra, a well-defined 4 f 7/2-5/2 doublet at 157.0eV-162.2eV, significant of the decomposition of the Bi(N(SiMe 3 ) 2 ) 3

precursor. The decomposition of the bismuth precursor at the very beginning of the synthesis (sample S1) while Bi was expected to be reduced by dihydrogen at 150 • C only, is certainly an important information. In parallel, the iron ionization peak shows a slow evolution from S1 to S4.

Iron is essentially in an oxidized state at room temperature while a strong metallic iron component arises only during step 2, i.e. after 4 hours under dihydrogen at 110 • C. Sample S3 (not shown here) still contains an important proportion of iron in an oxidized state. Finally, the Fe-2p peak displays a 2p 3/2-1/2 doublet at 707.1eV-720.3eV, clearly assigned to a non-oxidized environment, after 19h at high temperature in the S4 sample. The oxidized state of iron observed in the previous steps of the synthesis is probably related the adventitious hydrolysis of the Fe amido-complex, still present in S1 to S3 samples during sample preparation for XPS. Finally, XPS analyses of the S5 specimen (not shown) are similar to the S4 ones.

These results allow to conclude that the simultaneous presence of the two precursors in the solution definitely favors the bismuth precursor decomposition. In other words, the zero-valent Fe atoms resulting from the fast reduction of the iron precursor in the presence of amine-borane, induce the subsequent reduction of the bismuth precursor according to a redox process which can be summarized as follows, 3Fe 0 + 2Bi(N(SiMe 3 ) 2 ) 3 -→ 3Fe(N(SiMe 3 ) 2 ) 2 + 2Bi 0 clearly regenerating the Fe precursor in the reacting medium as long as there is still some bismuth precursor to reduce. In this first step of the reaction, the Fe amido complex acts as a catalyst for the reduction of the bismuth amido one. The bismuth precursor is soon completely decomposed (the Fe/Bi ratio is equal to 7.6) and the iron molecular species released in the solution become stable and can grow to form nanoparticles. The same sequence of samples was analyzed using WAXS and EXAFS. WAXS could be applied to all steps of the synthesis and provides insight to the average structure. On the left of EXAFS was also performed on all samples excepted S5, as it provides the unique ability to investigate short range order for both metals. Figure 9 presents the evolution of the uncorrected Fourier Transforms (FT) of the EXAFS functions at both edges. Study at Fe K edge provides limited information: at all steps, only one broad peak can be observed. One major change can however be observed: for S1, the broad peak points to still shorter distances, which indicates that it includes a higher amount of bonding distances between Fe and light elements, in agreement with the fact that iron could be mostly in a molecular state, as expected if a transmetalation process occurs in the first stages of the reaction. Study at the Bi edge clearly evidences the process: the peak at 0.22 nm (uncorrected), already noticed in the previous section, can be observed from S1

and already points to Fe in Bi environment. At this step, there is no evidence of pure Bi. The minor peak at 0.17 nm (uncorrected) points to a distance between Bi and a light element, which could correspond to some remaining bismuth precursor although this is not supported by the XPS analyses. At S2, we observe the growth of the Bi alloy and a strong decrease of the peak at 0.17 nm. At S3, this Fe-Bi material is the only phase including bismuth in the sample. A least-square refinement was performed on this spectrum and very good agreement was obtained with a single Bi-Fe distance at 0.274 nm (figure SI4 in supporting information). At S4, we can see that the Bi-Fe peak decreases and the Bi-Bi one finally appears, in agreement with WAXS investigations.

This evolution clearly indicates that the mixed environment observed for Bi atoms in the initial sample S4 is actually related to a mixture of pure metallic Bi with an amorphous Fe-Bi alloy. The precise composition of this Fe-Bi alloy cannot be estimated from our analyses, despite the large range of techniques implemented in this study. However, in their study of Fe-Bi films grown at room temperature by electron-beam deposition, Forester et al. 16 showed that the most stable alloys formed for a Bi concentration around 20%. Besides, the magnetic moment per iron atom in this alloy was found to reach at least 95 % of its value in bcc iron. Solid solutions with higher or lower Bi content were seen to undergo phase separation with the formation of pure iron or bismuth nanocrystals. It seems therefore reasonable to assume a 20% Bi content in the observed Fe-Bi alloy. According to the STEM-EDS analyses, the nanospheres core contains 50% bismuth, which strongly suggests that all the pure Bi nanograins are concentrated in this region together with the iron rich amorphous alloy. The weak Bi content (around 6%) evidenced in the shell region could remain in between Fe particles which were found to include some boron interstitial or substitutional atoms, at least in fresh specimens.

Conclusion

We have developped the synthesis of bimetallic Fe-Bi particles in the non-miscible Fe-Bi system.

A clear description of these particles and their properties was obtained through the use of complementary techniques as STEM-EDS, WAXS, EXAFS, XPS and M össbauer spectrometry. The particles, 150 nm of mean diameter, present a core@shell chemical order with a bismuth rich core and a shell formed by an assembly of 3 nm large iron nanoparticles. Such a chemical order, which was fully unexpected from both thermodynamic and kinetic considerations, has been explained through the analyses of the compound by XPS, WAXS and EXAFS experiments at different stages of the synthesis. Evidence has thus been given for the occurrence of a transmetalation process in these compounds, which occurs at the very first stages of the synthesis and favours the formation of an Fe-Bi amorphous alloy. When the Fe-Bi particles are formed, this alloy remains concentrated at their core together with pure bismuth grains. To our knowledge this is the first evidence for the formation of a Fe-Bi alloy in nanomaterials and using a soft chemistry route. Finally, we demonstrated the enhanced resistance to oxidation in these particles, despite their shell mainly formed by Fe nanoparticles. Good magnetic properties and resistance to air exposure as well as a relative innocuity of bismuth in biology give to these particles a very good potential for multifunctional technological applications .

Figure 2 :

 2 Figure 2: Local elemental mapping of Fe-Bi nanospheres giving evidence for an iron-rich shell and a concentration of Bi atoms in the core region (intensity corresponds to counts).

  22 nm) points to a different backscattering atom. Considering both atom size and phase shift, this indicates that in the whole sample, a large part of Bi atoms have Fe neighbours.[START_REF]) displayed the same kind of environment for Bi while the activation of Si[END_REF] Here also, very light elements like C, N, O and B cannot be excluded but should only marginally contribute to the left side of the main peak. Metallic bonding between Bi and Fe atoms should also be observed at Fe edge, and actually the right side of the main peak may include such a contribution, however quite low for two reasons: the low percentage of Bi in the material, but also the high static disorder in Fe neighbouring already reported. Since the relatively short Fe-Fe bond is already largely suppressed, an even more dramatic effect is expected in these conditions for the longer Fe-Bi one. The EXAFS results thus agree with a mixed Bi and Fe environment for Bi atoms. As most of the iron atoms are located in the bismuth poor shell region of the nanospheres (the shell volume is 5 times larger than the core volume and contains nearly 94% of iron), the observed Fe-Bi bonds at the Bi edge most probably originate from the core region whose composition was previously estimated close to 1/1.The fact that both Fe-Bi and Bi-Bi bonds exist gives a strong indication for the presence of pure Bi nanograins in the core region together with Fe-Bi alloy and possibly also pure Fe grains.

Figure 3 :

 3 Figure 3: Magnitudes of uncorrected Fourier transforms at Fe K edge (left) and Bi L 3 edge (right).For easier comparison, the reference foils were adjusted using a scaling factor (0.10 and 0.80 for Fe and Bi, respectively).

Figure 4 :

 4 Figure 4: Hysteresis loops of the Fe-Bi nanospheres recorded just after the synthesis (black solid line) and after 24 hours of air exposure (red dotted line).

Figure 6 :

 6 Figure 6: Sketch of the synthesis sequence and samples collected at the different steps, numbered from S1 to S5.

Figure 7 :

 7 Figure 7: Evolution of the XPS spectra in the Bi-4 f and Fe-2p regions through the different steps of the synthesis. Top: sample S1 in the presence of amine-borane, middle: sample S2 after 4h at 110 • C under dihydrogen, bottom: sample S4 after 19h at 150 • C under dihydrogen-underlined symbols were chosen for Bi • and Fe • to include Bi-Fe, respectively Fe-Bi, bonds in the pure metals peaks.

Figure 8 (

 8 Figure 8 (reciprocal space), the prominent feature is a broad amorphous pattern from S1 to S4. S5 is dominated by the rhombohedral Bi pattern, which gives evidence of the formation of extended Bi crystals (several nanometers) in this last step. The Radial Distribution Function (RDF) on the right of Figure 8 provides more information: from S1 to S3, well defined distances point to a phase quite different from both iron and bismuth. Since amplitudes in the RDF are strongly related to the atomic numbers of atoms engaged in the pair, we can expect distances related to Bi atoms to

Figure 8 :

 8 Figure 8: Evolution from bottom to top S1 to S5; left: in reciprocal space -bottom axis, bcc Fe pattern, reverse top axis Bi pattern; right: in real space -on top, simulation from bulk Bi.

Figure 9 :

 9 Figure 9: Evolution of the magnitudes of uncorrected Fourier transforms from S1 to S4 at Fe K edge (left) and Bi L3 edge (right) -for each edge, all curves are plotted at the same scale.

  • C and 19h at 150 • C without any stirring. A black magnetic solid was recovered by filtration, washed with anisole (5 mL) and dried under vacuum (recovered mass: 16mg; Fe: 11.18% ; Bi: 5.52 % i.e. atomic composition Fe 7.6 Bi 1 ). The same procedure was repeated and the reaction stopped at different reaction times to investigate the evolution of the material: before pressurization under H 2 , after 4h at 110 • C, after 4 and 72h hours at 150 • C.

Table 1 : -Comparison between refined values at 77K of hyperfine parameters (mean iso- mer shift

 1

< δ > and mean hyperfine field < B hyp >) obtained

by means of different fitting procedures involving gaussian and lorentzian profiles, G and L respectively.

	Sub-spectrum < δ > (mm.s -1 ) < 2ε > (mm.s -1 ) < B hyp > (T ) Area (1)/(2)
		±0.01	±0.01		±0.01
	(1) G / L	0.19 / 0.18	-0.02 / -0.00	32.2 / 33.0	67 / 57
	(2) G / L	0.32 / 0.23	-0.01 / 0.1	25.5 / 26.0	33 / 43

Table 2 : -Summary of refined values of hyperfine parameters obtained at 300 K, 77 K and 10 K (mean isomer shift

 2 < δ >, mean quadrupolar shift < 2ε >, mean canting angle < θ > and mean hyperfine field < B hyp >). δ > (mm.s -1 ) < 2ε > (mm.s -1 ) < θ > < B hyp > (T )

	Area (1)/(2)

<

Acknowledgement

This work has been supported by the ANR, the french National Research Agency (Mag@M project ANR-09-BLAN-0002-01). The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement nř 312284. The EXAFS experiments were carried out at the light source DORIS III at DESY, a member of the Helmholtz Association (HGF). We would like to thank Roman Chernikov for assistance in using beamline C and Edmund Welter in using beamline X1, as well as A. Mari and L. Rechignat for their help in SQUID measurements.

Supporting information available

The four figures referred in the text as SI1, SI2, SI3 and SI4 are available as supporting information, together with details on the control experiments with the iron and bismuth precursors. This information is available free of charge via the Internet at http://pubs.acs.org.

Notes and References


domain: math spi

text:
Olivier Pironneau

Pierre-Henri Tournier

Reflective Conditions for Radiative Transfer in Integral Form with Compressed H-Matrices

Keywords: MSC classification 85A25, 37N30, 31A10, 35Q30, 68P30, 74S05, Radiative Transfer, Reflective boundaries, Integral equation, H-Matrix, Finite Element Methods

In a recent article the authors showed that the radiative Transfer equations with multiple frequencies and scattering can be formulated as a nonlinear integral system. In this article, the formulation is extended to handle reflective boundary conditions. The fixed point method to solve the system is shown to be a monotone. The discretization is done with a P 1 Finite Element Method. The convolution integrals are precomputed at every vertex of the mesh and stored in 2K compressed hierarchical matrices, using Partially Pivoted Adaptive Cross-Approximation, where K is the number of different values of absorption. Then the fixed point iterations involve only matrix vector product. The method is O(N 3 √ N ln N ), with respect to the number of vertices, when everything is smooth. A numerical implementation is proposed and tested on a simple example. As there are some analogies with ray tracing the programming is complex.

Introduction

The Radiative Transport Equations (RTE) describe the behavior of electromagnetic radiation in a domain Ω as it interacts with matter [START_REF] Pomraning | The equations of Radiation Hydrodynamics[END_REF]. It is used to model a wide range of physical phenomena, including the propagation of light through plasma, tomography [START_REF] Tarvainen | Hybrid radiative-transfer-diffusion model for optical tomography[END_REF], atmospheric media [START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF], etc. The RTE is derived from the basic principles of quantum and statistical mechanics; it is a partial differential equation (PDE) that describes the distribution of radiation intensity in space, time and frequencies, coupled with a budget balance equation (BBE) for the electronic temperature. The PDE takes into account both absorption and scattering of radiation by matter, as well as emission of radiation by sources, which, in the present case, will be restricted to the boundaries of the emitting material. In [START_REF] Golse | Radiative transfer in a fluid[END_REF], [START_REF] Golse | Stratified radiative transfer in a fluid and numerical applications to earth science[END_REF] the authors have shown that the PDE can be converted into an integral equation for the total radiation at each point in the domain and that the coupling with the BBE can be handled by fixed point iterations. The method leads also to a general proof of existence, uniqueness and regularity of the solution. The difference with earlier studies such as [START_REF] Fan | Fast algorithms for integral formulations of steady-state radiative transfer equation[END_REF] is in the coupling with the equation for the temperature, BBE or even the full temperature PDE. In [START_REF] Golse | Radiative transfer for variable 3d atmospheres[END_REF] the authors have presented an implementation of the method using H-Matrix compression, a crucial ingredient which makes the evaluation of the integrals O(N 3 √ N ln N ) with respect to the number of vertices N in the 3D mesh which discretizes the domain Ω; N ln N is the complexity of the H-Matrix part but each element of the matrix requires an integral along a line in the domain. Compared with a brute force solution of the equations as in [START_REF] Jolivet | Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning[END_REF], the integral method has still a usable computing time for problems with frequency dependent physical parameters. However it did not handle reflective boundary conditions [START_REF] Siewert | On radiative-transfer problems with reflective boundary conditions and internal emission[END_REF]. H-Matrix compression [START_REF] Hackbusch | A sparse matrix arithmetic based on h-matrices. part i: Introduction to h-matrices[END_REF], [START_REF] Bebendorf | Approximation of boundary element matrices[END_REF], [START_REF] Boerm | Hybrid cross approximation of integral operators[END_REF], is a mathematical technique used to efficiently represent and manipulate large matrices that arise in a variety of applications. The technique uses a hierarchical structured representation of the matrices which allows a fast and accurate numerical computations when the integrals have a convolution type integrand which decays with the distance. H-Matrix compression works by partitioning a large matrix into smaller submatrices, and then approximate these submatrices by hierarchical low-rank approximations. The resulting representation allows for efficient matrix-vector multiplication, among other operations of linear algebra. The technique is particularly important and popular for computational electromagnetics in integral form such as boundary element methods. With the Partially Pivoted Adaptive Cross-Approximation (ACA) [START_REF] Boerm | Hybrid cross approximation of integral operators[END_REF] only the needed coefficients of the matrices are computed. However the theory requires geometrical smoothness [START_REF] Bebendorf | Hierarchical Matrices[END_REF].

We have extended the implementation done in [START_REF] Golse | Radiative transfer for variable 3d atmospheres[END_REF] using FreeFEM [START_REF] Hecht | New developments in freefem++[END_REF] and htool1 ; htools is a C++ toolbox for boundary element methods with electromagnetism in sight. FreeFEM is a popular open-source software package used for solving PDE systems with the finite element method (FEM). FreeFEM provides a wide range of pre-built FEM, as well as tools for mesh generation. It has a dedicated high level programming language that allows users to meet their specific needs. FreeFEM also supports parallel computing with mpi. One of the main advantages of FreeFEM for the present study is its ability to handle complex geometries and boundary condition, especially because of its powerful interpolations from volume to surface meshes.

Adding reflective conditions (RC) to the FreeFEM code presented in [START_REF] Golse | Radiative transfer for variable 3d atmospheres[END_REF] turned out to require solving the following difficulties:

• Integrate the RC into the integral formulation of the problem • Show that the fixed point iterations are still monotone.

• Find a formulation compatible with the use of H-Matrices

• Implement the method in the FreeFEM code. This paper presents the solutions found to resolve these 4 difficulties. It ends with a numerical test proposed in [START_REF] Kobayashi | 3d radiation transport benchmark problems and results for simple geometries with void region[END_REF].

The Radiative Transfer Equations

The problem is formulated in a domain Ω ⊂ R 3 with boundary Γ. The unit sphere in R 3 is called S 2 . One must find the radiation (called light from now on) intensity I ν (x, ω) at all points x ∈ Ω, for all directions all ω ∈ S 2 and all frequencies ν ∈ R + , satisfying:

ω ⋅ ∇I ν + κ ν I ν = κ ν (1 -a ν )B ν (T ) + κ ν a ν J ν , J ν ∶= 1 4π S 2 I ν dω , (1)
∞ 0 κ ν (1 -a ν )(J ν -B ν (T ))dν = 0 , (2)
I ν (x, ω) = R ν (x, ω)I ν (x, ω -2(n ⋅ ω)n) + Q ν (x, ω), on Σ ∶= {(x, ω) ∈ Γ × S 2 ∶ ω ⋅ n(x) < 0 }, (3)
where

B ν (T ) = ν 3 e ν T -1
is the (rescaled) Planck function. In the RC (3), R ν is the portion of light which is reflected and Q ν is the light source; n(x) is the outer normal of Γ at x. κ ν and a ν are the absorption and scattering coefficients; in general they depend on ν and x.

Example 1 If an object O inside a box B radiates because it is at temperature

T 0 , then, Ω = O B, Q ν = Q 0 [ω ⋅ n] -B ν (T 0 )
on O and zero elsewhere and Σ ⊂ ∂B × S 2 .

An Integral Formulation

For clarity we drop the subscript ν on κ, a and I. Assume that Ω is bounded and convex (see remark 1). Let

S ν (x) = κ(1 -a)B ν (T ) + κaJ ν , (4)
For a given x and ω, let τ x,ω be such that (x Σ (x, ω) ∶= xτ x,ω ω, ω) ∈ Σ; the method of characteristics tells us that

I(x, ω) = I(x Σ (x, ω), ω)e -∫ τx,ω 0 κ(x-ωs)ds + τx,ω 0 e -∫ s 0 κ(x-ωs ′ )ds ′ S ν (x -ωs)ds. (5)
Notice that τ x,ω = x Σ -x (see Figure 1), therefore, let

J ν (x) ∶= 1 4π S 2 I(x, ω)dω = S E ν (x) + J [S ν ](x) with S E ν (x) ∶= 1 4π S 2 I(x Σ (x, ω), ω)e -∫
τx,ω 0 κ(x-ωs)ds dω,

J [S](x) ∶= 1 4π S 2 τx,ω 0 e -∫ s 0 κ(x-ωs ′ )ds ′ S(x -ωs)dsdω = 1 4π Ω S(y) e -∫[x,y] κ y -x 2 dy. (6)
where, ω ′ (ω) ∶= ω -2(n ⋅ ω)n and

[x,y] f ∶= y -x 1 0 f (sy + (1 -s)x)ds.
To justify the last formula we refer to the following lemma with Ψ(x, y) = S(y)e -∫[x,y] κ . Again, for clarity, we drop the first argument x.

Lemma 1 Let Ω be a convex bounded open set of R 3 ; let Γ be its boundary. Let

Ψ ∶ Ω ↦ R be continuous. Let τ x,ω ≥ 0 be such that x -τ x,ω ω ∈ Γ, x ∈ Ω. Then S 2 τx,ω 0 Ψ(x -ωs)dsdω = Ω Ψ(y) y -x 2 dy. Proof : Denote Ψ the extension of Ψ by zero outside Ω. Let ω = (cos θ sin ϕ, sin θ sin ϕ, cos ϕ) T , θ ∈ (0, 2π), ϕ ∈ (-π 2 , π 2
). Consider a partition of the semi infinite line starting at x in direction -ω into segments of size δs and denote

x n = x -nδsω. Then S 2 τx,ω 0 Ψ(x -ωs)dsdω = lim δs→0 n>0 δs 2π 0 π 2 -π 2 Ψ(x n ) cos ϕdϕdθ = lim δs→0 n>0 2π 0 π 2 -π 2 Ψ(x n ) x -x n 2 x -x n 2 x n+1 -x n cos ϕdϕdθ (7)
We note that xx n 2 x n+1 -x n cos ϕdθdϕ is the elementary volume in the sector dθdϕ between the spheres centered at x and of radii x-x n and x-x n+1 . Therefore the right handside is an integral in

x ′ ∈ R 3 of Ψ(y) x-y 2 x -x n 2 . ◻ Remark 1
When Ω is not convex, on may apply the lemma to its convex closure Ω with κ extended to +∞ in Ω Ω.

Remark 2 When R ν ≡ 0, S E is known. As (2)defines a map T ∶ J ↦ T , T (x) = T [J ν ](x), ∀x ∈ Ω,
then, (4), ( 6) is a nonlinear integral formulation for J:

J ν (x) = S E ν (x) + J [κ(1 -a)B ν (T [J ν ]) + κaJ ν ](x), ∀x ∈ Ω. ( 8
)
The following fixed point method was shown in [START_REF] Golse | Radiative transfer in a fluid[END_REF] to be monotone and convergent:

J k+1 ν (x) = S E ν (x) + J [κ(1 -a)B ν (T [J k ν ](x)) + κaJ k ν ](x), k = 0, 1, . . . (9)
Let us extend these properties to the RTE with RCs.

For clarity let x Σ be short for x Σ (x, ω) and let ω ′ be

ω ′ (ω) = ω -2ω ⋅ n(x ′ ) n(x ′ ), with ω = x -x ′ x -x ′ .
Let us insert ( 5) and ( 3) in [START_REF] Golse | Radiative transfer in a fluid[END_REF]. Then,

S E ν (x) = S E ν,1 + S E ν,2 + S E ν,3 with S E ν,1 (x) ∶= 1 4π S 2 Q ν (x Σ , ω)e -∫
τx,ω 0 κ(x-ωs)ds dω,

S E ν,2 (x) ∶= 1 4π S 2 R ν (x Σ , ω)Q ν (x Σ (x Σ , ω ′ ), ω ′ ) e -∫ τ x Σ ,ω ′ 0 κ(x Σ -ω ′ s)ds e -∫ τx,ω 0 κ(x-ωs)ds dω, S E ν,3 (x) ∶= 1 4π S 2 R ν (x Σ , ω)e -∫ τx,ω 0 κ(x-ωs ′ )ds ′ τ x Σ ,ω ′ 0 e -∫ s 0 κ(x Σ -ω ′ s ′ )ds ′ S ν (x Σ -ω ′ s)ds dω.
Hypothesis 1 Let us rule out multiple reflections and focal points

1. If R ν (x Σ (x, ω), ω) > 0, then R ν (x Σ (x Σ (x, ω), ω ′ ), ω) = 0.
2. Given x and y, there is only a finite number

N of x ′ n ∈ Γ such that [x ′ n , y] is the reflected ray of [x, x ′ n ].
Note that x ′ n depends on x and y.

Proposition 1 Under Hypothesis 1

S E ν,3 (x) ∶= M n=1 1 4π Ω R ν (x ′ n , x-x ′ n x-x ′ n ) e -∫[x,x′ n ]∪[x ′ n ,y] κ ( x -x ′ n + x ′ n -y ) 2 S(y)dy Proof Let x(s) ∶= x Σ -ω ′ s.
By Lemma 1, (7) can be written as

S 2 τ x Σ ,ω ′ 0 R ν (x Σ , ω)e -∫[x,x Σ ]∪[x Σ ,x(s)] κ S(x(s))ds dω = Ω R ν (x Σ , ω)S(y) e -∫[x,x Σ ]∪[x Σ ,y] κ ( x -x Σ + x Σ -y ) 2 dy, provided that [x Σ , y] is reflected from [x, x Σ ]
. Now, by hypothesis, if x and y are given in Ω there are only a finite number of x Σ ∈ Γ for which [x Σ , y] is reflected from [x, x Σ ], (see Figure 1). ◻ Proposition 2 Let Hypothesis 1 hold. Then the source terms from the boundaries are

S E ν,1 (x) = 1 4π Γ Q ν (y, y-x y-x ) [(y -x) ⋅ n(y)] - y -x 3 e -∫[x,y] κ dΓ(y), (10)
S E ν,2 (x) = M n=1 1 4π Γ R ν (x ′ n , x-x ′ x-x ′ )Q ν (y, x ′ n -y x ′ n -y ) [(x ′ n -y) ⋅ n(y)] -e -∫[x,x′ n ]∪[x ′ n ,y] κ x ′ n -y ( x -x ′ n + x ′ n -y ) 2 dΓ(y). ( 11
)
Recall that x ′ n depends y.

Proof : Recall that a solid angle integral at x of a surface Σ is

S 2 f (x, x ′ )dω(x ′ ) = Σ f (x, x ′ ) [(x -x ′ ) ⋅ n(x ′ )] - x -x ′ dΣ(x ′ ) x -x ′ 2 ,
Hence from the definition of S E ν,2 above we see that (10) holds. To prove [START_REF] Kobayashi | 3d radiation transport benchmark problems and results for simple geometries with void region[END_REF] we start from the definition of S E ν,2 above. For clarity let us assume that Q ν and R ν do not depend on ω. Observe that if a ray from x in the direction -ω does not hit, after reflection at x ′ on some Γ R , a boundary Γ Q at y where Q ν (y) is non zero then that ω does not contribute to S E ν,2 . So we can use the solid angle of Γ Q . However the solid angle is not seen from x but from x, the symmetric of x with respect to the tangent plane of Γ R at x ′ . As the distance from x to y is also xx ′ + x ′ -y , we obtain [START_REF] Kobayashi | 3d radiation transport benchmark problems and results for simple geometries with void region[END_REF].

◻

Corollary 1 J ν (x) = SE ν (x) + J [S ν ](x), (12)
with SE ν (x) ∶= S E ν,1 (x) + S E ν,2 (x) given by Proposition 2 and

J [S](x) = 1 4π Ω ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ e -∫[x,y] κ y -x 2 + M n=1 e -∫[x,x′ n ]∪[x ′ n ,y] κ ( x -x ′ n + x ′ n -y ) 2 R ν (x ′ n , x-x ′ n x-x ′ n ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ S(y)dy (13)

Example

Assume that Γ = Γ Q ∪ Γ R and Q ν (x, ω) = [ω ⋅ n(x)] -Q 0 with Q 0 > 0 on Γ Q and 0 on Γ R . Assume R ν (x, ω) = R 0 with R 0 > 0 on Γ R and 0 on Γ Q . Assume that
there is never more than one reflection point on Γ R , i.e. N = 1. Then

SE ν (x) = Q 0 4π Γ Q ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ [(y -x) ⋅ n(y)] - y -x 2 2 e -∫[x,y] κ + R 0 ([(x ′ 1 -y) ⋅ n(y)] -) 2 e -∫[x,x′ 1 ]∪[x ′ 1 ,y] κ x ′ 1 -y 2 ( x -x ′ 1 + x ′ 1 -y ) 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ dΓ(y), J [S](x) = 1 4π Ω ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ e -∫[x,y] κ y -x 2 + R 0 e -∫[x,x′ 1 ]∪[x ′ 1 ,y] κ ( x -x ′ 1 + x ′ 1 -y ) 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ S(y)dy. x x ω x Σ ∶= x Σ (x, ω) ω ′ = ω -2(n(x Σ ) ⋅ ω)n(x Σ ) n(x Σ ) x ΣΣ ∶= x Σ (x Σ , ω ′ ) n(x ΣΣ )
Figure 1: In this configuration the source Γ Q is the upper square. An RC is imposed on the lower plane Γ R . S E ν has an integral of the solid angle of the upper square seen from x plus an integral of the solid angle of the upper square seen from x, the symmetric of x with respect to Γ R .

Fixed Point Iterations

Consider the fixed point iterations initialized with T 0 and J 0 = 0.

Algorithm For k = 0, 1, . . . : Set S k ν (x) = κ(1 -a)B ν (T k ) + κaJ k ν , Set J k+1 ν (x) = SE ν (x) + J [S k ν ](x). Compute T k+1 by solving (using Newton algorithm) for each x ∈ Ω ∞ 0 κ ν (1 -a ν )(J k+1 ν -B ν (T k+1
))dν = 0.

(

) 14
Then

S ν,j ∶= aJ ν,j + (1 -a)B ν (T j ), J ν,i ∶= SE ν,i + j G ij κ S ν,j
where

G ij κ = 1 4π Ω ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ κ e -∫[xi-y] κ x i -y 2 dy + M n=1 R 0 (x ′ n ) e -∫[xi,x′ n ]∪[x ′ n ,y] κ ( x i -x ′ n + x ′ n -y ) 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ŵj (y)dy and where SE ν,i = 1 4π Γ Q 0 (y) ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ [(x i -y) ⋅ n(y)] - x i -y 2 2 e -∫[xi,y] κ + M n=1 R 0 (x ′ n ) ([(x ′ n -y) ⋅ n(y)] -) 2 e -∫[xi,x′ n ]∪[x ′ n ,y] κ x ′ n -y 2 ( x i -x ′ n + x ′ n -y ) 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ dΓ(y)
The integrals are approximated with quadrature at points {x j q } Mq 1 . The points are inside the elements; consequently x ix j q is never zero. A formula of degree 5, with M q = 14, is used when x iy is small and of degree 2, with M q = 4, otherwise; the results do not change when higher degrees are used. Luckily when x i is closed to Γ an analytical formula can be used [START_REF] Golse | Stratified radiative transfer in a fluid and numerical applications to earth science[END_REF].

To compute x ′ n such that [y, x ′ n ] is the reflected ray of [x ′ n ,
x i ] a loop on all the elements of the reflecting boundaries is necessary. This can be expensive, but in the case of planar reflective boundaries the symmetric point xi is easy to compute and so is the intersection of [x i , y] with the reflective boundary. Finally, to the vector { SE ν,i } N i=1 we associate a matrix { SE i,l } N,L i,l=1 by replacing Q 0 (y) above by wl (y). Then:

Q 0 (y) = L 1 Q 0 l wl (y) ⇒ SE ν,i = L 1 SE i,l Q 0 l .

Compression

So, for each ν we have two very large matrices, { Ḡi,j } N,N i,j=1 and { SE i,l } N,L i,l=1 . Remark 4 Note that for each value of ν two matrices are needed. However on close inspection it is really two matrices for each value of κ ν . Very often, less than ten values are sufficient to represent a general κ ν .

Both matrices can be compressed as H-matrix [START_REF] Bebendorf | Hierarchical Matrices[END_REF], [START_REF] Rjasanow | The Fast Solution of Boundary Integral Equations[END_REF], [START_REF] Sauter | Boundary Element Methods[END_REF] (and the references therein) so that the matrix-vector product has complexity O(N ln N ). The method works best when the kernel in the integrals, decays with the distance between x i and y. In both matrices the kernel decays with the square of the distance. The H-matrix approximation views G as a hierarchical tree of square blocks. The blocks correspond to interactions between clusters of points near x j and near x i . A far-field interaction block can be approximated by a low rank matrix because its singular value decomposition (SVD) has a fast decaying singular value. We use the Partially Pivoted Adaptive Cross-Approximation (ACA) [START_REF] Boerm | Hybrid cross approximation of integral operators[END_REF], to approximate the first terms of the SVD of the blocks, because only r-rows times r-columns are needed instead of the whole block, where r is the rank of the approximation. The rank is a function of a user defined parameter connected to the relative Frobenius norm error. Another criteria must be met: if R 1 (resp. R 2 ) is the radius of a cluster of points centered at x 1 (resp. x 2 ), then one goes down the hierarchical tree until the corresponding block satisfies max(R 1 , R 2 ) < η x 1 -x 2 where η is a user defined parameter. If the end of the tree is reached, the block is not compacted. The precision is not guaranteed on account of [(x -y) ⋅ n(y)] -if the jump from one triangular face to another is large. A similar singularity caused by normals is analyzed for a double layer potential formulation in [START_REF] Bebendorf | Hierarchical Matrices[END_REF] (Example 3.38, p148) and a remedy is proposed. To check whether this remedy is needed here we ran two cases, one without compression and one with 97% compression. No difference was observed. Convergence versus iteration number n. When the scaled temperature is initialized to T 0 = 0.001 at n = 0 the convergence is monotonously increasing. When T 0 = 0.44 the convergence is monotonously decreasing.

A Test Case

In [START_REF] Kobayashi | 3d radiation transport benchmark problems and results for simple geometries with void region[END_REF] a semi-analytic solution of the RTE is given for a geometry shown on Figure 3. In this test a = 0 and κ is not a function of ν. Hence the grey formulation can be used. Define Ī = ∫ ∞ 0 I ν dν. By averaging [START_REF] Rjasanow | The Fast Solution of Boundary Integral Equations[END_REF] in ν and due to the Stefan-Boltzmann relation, the following holds (for clarity a ν = 0):

∞ 0 B ν (T )dν = σT 4 with σ = π 4 15 ⇒ Jk+1 (x) = SE (x) + J [κσ(T k ) 4 ](x), κσ(T k+1 ) 4 = κ Jk+1 . (15)

The Geometry

The outer container is D = (0, 60) × (0, 100) × (0, 60), in cm. A cube C = [0, 10] 3 radiates with intensity Q 0 = 0.1. A rectangular cylinder prolonging the radiating cube (0, 10) × (10, 100) × (0, 10) has a low absorption κ = 10 -4 while the rest has κ = 0.1. In Kobayaski's test case 1A there is no scattering, and the three planes containing the origin reflects the radiations: (O, x, z), (O, x, y), (O, y, z).

Unfortunately the present method cannot handle volumic radiating region. Consequently we have kept the geometry but only the 3 faces of C inside Ω radiate in all directions ω with intensity Q 0 [ω ⋅ n] -, where n is the normal to the cube's face pointing inside the cube. The domain is Ω = D C. To assert the precision of the method we consider first only one reflective plane, Γ R = (0, y, z). So, R ν = 1 on that plane and zero elsewhere and Q 0 is non zero only on the faces of the C. We also take κ = 0.1 everywhere. First we verify that the convergence is monotone increasing if T 0 is small and monotone decreasing if T 0 is large (Figure 2). Note that the monotone increasing sequence converges faster.

Results

Next, we compare the results with a computation on a domain D = (-60, 60) × (0, 100) × (0, 60) which is D plus the symmetrized of D with respect to the plane (0, y, z). This is because reflection on a plane is equivalent to extending the domain by symmetry with respect to that plane. Figure 4 shows level surfaces of J computed on the symmetrized domain (but restricted to the original domain) and compared with the same level surfaces but computed with the RC. Surfaces with similar colors should be near each other. In fact the difference is not viible except near z = 0 . Figure 5 shows level surfaces of J computed with the RC, κ = 10 -4 or 0.1, and the same level surfaces but computed without any RC on the (O, y, z) plane. It is seen that surfaces with similar color are far from each others. By comparing Figure 4 with Figure 5 we see that the RC does almost the same as symmetry and that no condition at all is a non viable approximation for this problem. Similarly, figure 6 shows x ↦ J(x, 15, 15), computed on the symmetrized domain, or with the radiative condition or without it.

Finally, Figure 7 shows x ↦ J(x, 15, 15), computed with the RC on 3 meshes, coarse, medium and fine. The same 3 meshes are used in Table 1 where the theoretical complexity N

Kobayashi's Test 1A

Test 1A of [START_REF] Kobayashi | 3d radiation transport benchmark problems and results for simple geometries with void region[END_REF] has been computed, i.e. non constant κ and 3 reflective planes. Surface levels of J are shown on Figure 10. The comparison with the data in [START_REF] Kobayashi | 3d radiation transport benchmark problems and results for simple geometries with void region[END_REF] on the line (5, y, 5) is shown on Figure 8. Finally the L 2 error is computed using values on the line (x, 15, 15). The reference solution to compute the error is the solution on a mesh with N=26189 .

The results are displayed on Figure 9. It shows an L 2 -error on the line O(h 3 ). [START_REF] Kobayashi | 3d radiation transport benchmark problems and results for simple geometries with void region[END_REF]. A scaling is applied to each curve so that they coincide at y = 15 (because [START_REF] Kobayashi | 3d radiation transport benchmark problems and results for simple geometries with void region[END_REF] is given for volumic source and the present method handles only surface sources). ). The reference solution is computed on a mesh with 26189 vertices. The plotted points are computed on meshes with N as in Table 1.

Conclusion

Compressed H-matrices is an ideal tool for RTE in integral form because the complexity of the method is O(N 3 √ N ln N ) where N is the number of vertices in the 3D mesh and because it can handle frequency dependent absorption and scattering coefficients at the expense of a finite number of compressed matrices and a finite number of matrix-vector products. The integral nonlinear formulation of RTE proposed in [START_REF] Golse | Radiative transfer for variable 3d atmospheres[END_REF] has been extended to cases with reflective boundary conditions. The monotonicity property of the iterative solver is kept. The discretization with the finite element method is the same. However it is hard to write a general computer code because of the complexity of potential multiple reflections, as in Ray Tracing. Hence in this article the numerical validation has been done only for a finite number of plane reflective boundaries and for an academic geometry.

Figure 2 :

 2 Figure 2: Values of J along the y axis at x = z = 15 computed with a RC.Convergence versus iteration number n. When the scaled temperature is initialized to T 0 = 0.001 at n = 0 the convergence is monotonously increasing. When T 0 = 0.44 the convergence is monotonously decreasing.

Figure 3 :

 3 Figure 3: A small cube (colored blue and red on the figure) radiates normally to its faces in a medium which has a very small absorption coefficient κ = 10 -4 in the cylinder prolonging the cube and κ = 0.1 elsewhere.

Figure 4 :

 4 Figure 4: Level surfaces of J using a log scale computed with κ = 0.1 and only one reflective plane, (0, y, z) facing us, slightly to the left. Comparison between a computation done with the RC and a computation done on a symmetrized domain, double in size.Surfaces of equal colors are so near each other that it is hard to distinguish them.

Figure 5 :

 5 Figure 5: Same as in Figure 4 but the RC is not used in one computation. Surfaces of equal colors are far from those using the RC, indicating the absolute necessity of a RC.Here κ is not constant.

3 √Figure 6 :Figure 7 :

 367 Figure 6: Values of J along the x axis at y = z = 15 computed by different methods on the coarse mesh.

Figure 8 :

 8 Figure 8: Values of J versus y ≥ 15 at x = z = 5 and comparison with the values given in[START_REF] Kobayashi | 3d radiation transport benchmark problems and results for simple geometries with void region[END_REF]. A scaling is applied to each curve so that they coincide at y = 15 (because[START_REF] Kobayashi | 3d radiation transport benchmark problems and results for simple geometries with void region[END_REF] is given for volumic source and the present method handles only surface sources).

5 Figure 9 :

 59 Figure 9: Log-log plot of L 2 error versus N . The linelog N +1.5 indicates an error O( 1 N). The reference solution is computed on a mesh with 26189 vertices. The plotted points are computed on meshes with N as in Table1.

Figure 10 :

 10 Figure 10: Kobayashi's test: Level surfaces of J using a log scale. The reflective planes are the (O, x, z), (O, y, z), (O, x, y). The origin O is the lower right corner.

Table 1 :

 1 CPU time, complexity and compression for 3 meshes with one reflective plane and non constant κ. But since the radiative sources are different (volumic in Kobayashi's and surfacic in our case) we have scaled the result with Kobayashi's value at x = 5, y = 15, z = 5 by their values at y = 15.

	N vertices Surf. Space Saving Vol. Space Saving. CPU	10 5 CP U N √ 3 N log N
	2758	0.41	0.60	5.2"	1.6971
	8003	0.59	0.77	13.9"	0.97
	26189	0.67	0.89	60.7"	0.77
	84042	0.95	0.74	312"	0.75

https://github.com/htool-ddm/htool

Acknowledgement

We would like to thank warmly Frédéric Hecht, for his constant good will to adapt FreeFEM to our needs.

Proposition 3 Let {J * ν , T * } be the solution. If T 0 (x) > T * (x), ∀x ∈ Ω then the iterations are monotone decreasing: T k (x) > T k+1 > T * (x), ∀x ∈ Ω. Conversely if T 0 (x) < T * (x), ∀x ∈ Ω then the iterations are monotone increasing: T k (x) < T k+1 < T * (x), ∀x ∈ Ω.

Proof : Le us prove it for the monote increasing sequence. By subtracting the definition J k ν from that of J k+1 ν and using the linearity of J , we obtain

for all x then J k+1 ν (x) > J k ν (x). The equation for T k+1 is also monotone in the sense that

One sure way to impose it is to choose T 0 = 0 and J 0 = 0.

To prove that T k < T * we observe that

Henceforth, convergence and uniqueness can probably be proved as in [START_REF] Golse | Stratified radiative transfer in a fluid and numerical applications to earth science[END_REF], but there are technical difficulties of functional analysis which may not be appropriately discussed here.

FEM discretization and Compressed H-Matrices

For clarity consider exemple 1.2. As the values of Q 0 and R 0 depend on the boundary name, we write Q 0 (x) and R 0 (x).

The domain Ω is discretized by a tetraedral mesh; the boundary Γ is discretized by triangular mesh, not necessarily conforming with the volume mesh. Let {x j } N 1 be the vertices of the tetraedra of Ω and {x l } L 1 the vertices of the triangles of Γ . A continuous P 1 interpolation of J on the tetraedral mesh is:

where ŵj is the P 1 -Finite Element hat function of vertex x j .

Declarations

• As to the specific input of each author, P.-H. Tournier provided the Hmatrices part and the interface between htool and FreeFEM; the rest of the program and the theory were done together.

• The authors have no relevant financial or non-financial interests to disclose.

• The authors have no conflicts of interest to declare that are relevant to the content of this article.

• All authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript.

• The authors have no financial or proprietary interests in any material discussed in this article.