Valentin Jean
email: valentin.jean@univ-lorraine.fr

S Fumeron

K Termentzidis

X Zianni
email: xzianni@teihal.gr

David Lacroix
email: david.lacroix@univ-lorraine.fr

Monte Carlo simulations of phonon transport in Si nanowires with constrictions

published or not. The documents may come

Monte Carlo simulations of phonon transport in Si nanowires with constrictions

Introduction

During the past decade a broad attention has been paid to heat transport properties in nanostructured materials [START_REF] Cahill | Nanoscale thermal transport[END_REF][START_REF] Cahill | Nanoscale thermal transport. II. 2003-2012[END_REF][START_REF] Luo | Nanoscale heat transfer -from computation to experiment[END_REF][START_REF] Yang | Thermal transport in nanostructures[END_REF], and more specifically in semiconductors compounds that are used in several technologies.

Structures such as nanofilms, nanowires or nanoporous media have been extensively studied, experimentally and numerically. In order to understand the physical mechanisms that underlie heat transfer in these specific structures, one needs to study heat carriers properties such as phonon dispersion relations, the phonon lifetime and the density of states (DOS). Using these basic inputs, thermal properties such as the thermal conductivity (TC) or the inter-facial thermal resistance (Kapitza resistance) can be appraised. In many applications tailoring TC through a dedicated processing is considered as the most convenient way to elaborate materials with high, low and/or anisotropic TC. For example, in polymers carbon nanotube (CNT) inclusions can help increase the TC [START_REF] Han | Thermal conductivity of carbon nanotubes and their polymer nanocomposites: A review[END_REF]. Oppositely, in semiconductors, alloying, doping with substitution elements or nano-structuring are the best ways to reduce it [START_REF] Maldovan | Thermal energy transport model for macro-to-nanograin polycrystalline semiconductors[END_REF][START_REF] Wang | Thermal Conductivity of Nanocrystalline Silicon: Importance of Grain Size and Frequency-Dependent Mean Free Paths[END_REF]. In the latter case, several studies have pointed out that TC reduction can be addressed by increasing the phonon scattering. In this framework two regimes exist, the phonon wave confinement regime and the structure confinement. In the first case, the phonon wavelength is close to the characteristic length of the structure. For instance, in superlattices and phononic crystals, these wave effects lead to a modification of the dispersion properties (band folding) that reduces the material TC [START_REF] Ren | Thermal conductivity of superlattices[END_REF][START_REF] Yang | Partially coherent phonon heat conduction in superlattices[END_REF][START_REF] Zen | Engineering thermal conductance using a two-dimensional phononic crystal[END_REF]. In the second case, the structure characteristic size is shortened but remains large when compared to the phonon wavelength. In that case, phonon boundary scattering is an important issue that rules heat transport. In the well known case of Si nanowires, TC can be decreased by at least an order of magnitude [START_REF] Li | Thermal conductivity of individual silicon nanowires[END_REF][START_REF] Chen | Thermal Conductance of Thin Silicon Nanowires[END_REF].

Recently, nanowires were engineered in order to enhance phonon confinement through the modulation of their diameter [START_REF] Blanc | Phonon heat conduction in corrugated silicon nanowires below the Casimir limit[END_REF][START_REF] Maire | Reduced thermal conductivities of Si one-dimensional periodic structure and nanowire[END_REF]. Models and simulations [START_REF] Nika | Thermal conductivity inhibition in phonon engineered core-shell cross-section modulated Si/Ge nanowires[END_REF][START_REF] Termentzidis | Modulated SiC nanowires: Molecular dynamics study of their thermal properties[END_REF] were also developed to appraise the magnitude of the TC lowering assuming simple square cross-section modulated nanowires by means of kinetic theory or molecular dynamics. In these studies, cavities and constrictions are supposed to be the origin of the TC drastic decreases as the phonon group velocity is strongly lowered with the trapping of phonon modes. In this "abrupt" configuration (nanowire with rectangular modulation) heat flux, which is related to transported energy quanta along the nanowire axis, is strongly perturbed. Furthermore, heat transport in constrictions is also a major issue which has been considered to described the thermal resistance between a nano-sized object and a surface in the case of non perfect contact [START_REF] Shi | Thermal Transport Mechanisms at Nanoscale Point Contacts[END_REF][START_REF] Volz | Increase of thermal resistance between a nanostructure and a surface due to phonon multireflections[END_REF][START_REF] Singh | Phonon Transport Across Mesoscopic Constrictions[END_REF]. This is typically the case in scanning thermal microscopy (SThM) where the edge of the tip is often modeled as a constriction between two dissimilar materials. Constrictions have also been investigated theoretically to understand the diffusive/ballistic transport mechanism occurring through small apertures [START_REF] Prasher | Predicting the Thermal Resistance of Nanosized Constrictions[END_REF][START_REF] Prasher | Diffraction-limited phonon thermal conductance of nanoconstrictions[END_REF][START_REF] Sanjoy | Molecular dynamics simulation of thermal transport at a nanometer scale constriction in silicon[END_REF]. In the latter studies, phonon diffraction in the case of very small apertures was modeled. Constrictions are also present in nanojunctions and are often parts of micro electromechanical systems (MEMS). Recently, it has been experimentally shown that controlling their shape and size is important to adjust heat transport properties for several applications [START_REF] Jalabert | Ballistic Thermal Conductance of a Lab-in-a-TEM Made Si Nanojunction[END_REF].

Nanowires with realistic constrictions are difficult to simulate properly because of several numerical limitations. In molecular dynamics based studies, dealing with complex geometries and micrometer length scales cannot be simply done due to the currently available computational resources. On the other hand, modeling of phonon transport through the resolution of the Boltzmann transport equation (BTE) by finite volumes techniques is often limited to simple Cartesian or cylindrical geometries. In such an approach, the description of a realistic shape for a constriction might be challenging. Eventually, techniques based on the use of a transmission function to calculate heat flux through the constriction and the corresponding conductance, which are often faster, do not give information about temperature distribution in the nanostructures. Those issues or limitations are no longer existing when phonon motion and scattering is treated with the Monte Carlo (MC) simulation of the BTE. The present work is dedicated to the study of heat transfer through nanowires with different constriction shapes and lengths. Through MC numerical simulations, heat and temperature profiles in the nanostructures are recovered for different constriction profiles. Then, equivalent Fourier's conductance and TC are determined.

The paper is organized as follows. Following the introduction, theoretical models used to appraise the thermal conductance are recalled. Then, the constriction modeling is presented as well as the MC simulation tools used for numerical simulation of phonon transport. In section 3, conductance and TC for different constriction shapes are discussed and compared to existing theory. The particular case of a long constriction, similar to a nanowire between two heat reservoirs is addressed in subsection 3.4. Finally, conclusions and perspectives to this work are given.

Models and simulation tools

This section deals with models and simulations used to study constrictions in nanostructures. First, theoretical models giving expressions of the thermal conductance through apertures are recalled. Then, the model used in the framework of MC phonon transport is detailed.

Constriction models

Conductance or thermal resistance due to a "point contact" constriction can be theoretically assessed in two limiting cases, the diffusive limit and the ballistic limit. Those limits are characterized by the phonon scattering mean free path (mfp) Λ, compared to the characteristic size of the constriction. Usually for bulk materials, this mfp is appraised considering all the possible scattering mechanisms, namely the three phonon processes, Normal (N) and Umklapp (U), as well as the impurity or defect scattering. For bulk silicon, which is the material studied in the present work, at 300 K, Λ Si b is in the range of [170-250] nm according to the chosen phonon dispersion properties and the calculation method [START_REF] Chen | Thermal conductivity and ballistic-phonon transport in the cross-plane direction of superlattices[END_REF][START_REF] Minnich | Modified effective medium formulation for the thermal conductivity of nanocomposites[END_REF][START_REF] Jean | Monte Carlo simulations of phonon transport in nanoporous silicon and germanium[END_REF]]. Yet, for nanowires the phonon mfp is mainly ruled by phonon confinement and boundary scattering. Thus, it is commonly accepted that Λ nw = d, where d is the nanowire diameter.

Diffusive limit

In "the diffusive limit", the phonon mfp is smaller than the aperture of the constriction (minimum diameter) Λ < d min . In this case, the Maxwell's formula [START_REF] Sanjoy | Molecular dynamics simulation of thermal transport at a nanometer scale constriction in silicon[END_REF] holds and the thermal conductance G dif (or the thermal resistance R dif ) for a circular aperture reads

G dif = 1 R dif = kd min ( 1
)
where k is the thermal conductivity of the media on both sides of the constriction. Considering the Knudsen number Kn, defined as Kn = Λ/d min , the diffusive limit corresponds to small values of Kn. When the shape of the constriction changes, the diffuse conductance also changes [START_REF] Prasher | Predicting the Thermal Resistance of Nanosized Constrictions[END_REF].

Ballistic limit

In "the ballistic limit", the phonon mfp is much larger than the constriction aperture Λ > d min . In this case, the phonon transport across the constriction can be understood as the flow rate of gas molecules in the free molecular flow regime [START_REF] Prasher | Predicting the Thermal Resistance of Nanosized Constrictions[END_REF]. Two limit cases can be distinguished depending on the phonon dominant wavelength λ dom (T ) ≃ hv s /4.25k B T at a given temperature, where h is the Planck constant, v s refers to the speed of sound (here 6500 m/s for Si), and k B is the Boltzmann constant [START_REF] Blanc | Phonon heat conduction in corrugated silicon nanowires below the Casimir limit[END_REF]. First, if the size of the aperture d min is larger than λ dom the conductance can be derived from the net heat flux through the constriction.

Assuming constant group velocity and small temperature gradient on both sides of the constriction it reads

G bal = 1 R bal = 3kA 4Λ (2)
where A = πd 2 min /4 is the cross section of the aperture. G bal depends on the inverse of the phonon mfp. This expression was first derived from electron transport within constriction and named Sharvin conductance.

In the second case (which will be not considered here), d min is similar or smaller than λ dom and diffraction shall be accounted for. In the present study, simulations are mostly carried out at 300K, λ dom (T ) ≃ 0.25 nm and the smallest aperture diameter considered here is more than ten times larger: therefore, diffraction can be reasonably considered as negligible.

Monte Carlo simulation tool

Monte Carlo (MC) calculations of the phonon transport were performed using a dedicated home-built simulation tool. The basic specifications of the numerical method were extensively described in several papers [START_REF] Mazumder | Monte Carlo Study of Phonon Transport in Solid Thin Films Including Dispersion and Polarization[END_REF][START_REF] Péraud | Efficient simulation of multidimensional phonon transport using energybased variance-reduced Monte Carlo formulations[END_REF] for bulk material and successfully derived for nanowires and nanoporous thin films [START_REF] Lacroix | Monte Carlo simulation of phonon confinement in silicon nanostructures: Application to the determination of the thermal conductivity of silicon nanowires[END_REF][START_REF] Jean | Monte Carlo simulations of phonon transport in nanoporous silicon and germanium[END_REF][START_REF] Zianni | Scaling behavior of the thermal conductivity of width-modulated nanowires and nanofilms for heat transfer control at the nanoscale[END_REF]. The method lies on the resolution of the Boltzmann transport equation (BTE) in the frame of the relaxation time approximation (RTA). Bulk silicon dispersion properties have been used for the sampling of phonon frequency and group velocity, as well as for the calculation of the scattering relaxation times following the model proposed by Holland [START_REF] Holland | Analysis of lattice thermal conductivity[END_REF]. Relaxation time approximation is supposed to be valid as spatial discretization is chosen in order to ensure that sampled phonon do not cross several cells over a time step and thus have a non negligible probability to scatter. Besides, in the present work, most of the studied nanowire and constriction diameters are sufficiently large to ensure that bulk properties can be used for phonon transport. For the thinner nanowire (d=28 nm), the dispersion properties within the constriction are probably different. For instance, when looking at molecular dynamic studies dealing with this issue [START_REF] Huang | An investigation of the phonon properties of silicon nanowires[END_REF], for nanowires with diameters lower than 4 nm, supplementary modes at low frequencies with high wave-numbers appear. This means that phonon with lower energy and group velocity are dominant and thus heat transport should be lowered. Then in some peculiar cases of this study, the thermal conductivity and conductance might be a little bit overestimated.

The BTE for phonons is related to the variations of the distribution function f (T, r, p, K) which depends on time t, location r, polarization p and wave vector K. In the absence of external force, the BTE expression according to the RTA assumption is a classical transport equation in which time variation of f plus the advection term equal a collisional term [START_REF] Ashcroft | Solid state Physics[END_REF]

∂f ∂t + ∇ K ω • ∇ r f = f 0 -f τ (3)
where

f 0 (T, p, K) = 1/ [exp ( ω(p, K)/k B T ) -1]
is the Bose-Einstein equilibrium distribution with the reduced Planck's constant. τ (ω, p, T, L) is the total phonon scattering relaxation time, which depends on the phonon dispersion properties, on the temperature and possibly on a characteristic length L. In the framework of MC modeling, Eq. ( 3) is solved in two steps. First, the advective part of the BTE is addressed (left hand side term) assuming that there is no scattering. A statistically representative population of phonons is drawn (frequency, polarization and group velocity) from a cumulative distribution function based on the phonon energy at a given temperature [START_REF] Péraud | Efficient simulation of multidimensional phonon transport using energybased variance-reduced Monte Carlo formulations[END_REF][START_REF] Jean | Monte Carlo simulations of phonon transport in nanoporous silicon and germanium[END_REF]. Then each energy bundle is allowed to drift with respect to: i) its group velocity within the nanostructure, ii) the chosen time step ∆t, iii) the internal reflections at boundaries. Secondly, the phonon scattering step is achieved considering the right hand side term of Eq. ( 3). The characteristic relaxation time is given according to the "Matthiessen" summation rule. It takes into account: normal (N), umklapp (U) and impurity (I) scattering rates according to M.G. Holland [START_REF] Holland | Analysis of lattice thermal conductivity[END_REF][START_REF] Jean | Monte Carlo simulations of phonon transport in nanoporous silicon and germanium[END_REF] formalism. Knowing the energy distribution in the nanostructure one can compute in a control volume the energy and subsequently the temperature. Furthermore, phonon heat flux along the nanowire axis is given by

ϕ z = 1 V i=N b i=1 ω i × v gzi ( 4
)
where V is the control volume, N b the number of phonon energy bundles in V and v gzi the group velocity in the z direction (nanowire axis) of the phonon with frequency ω i . Then, considering the Fourier's formalism an equivalent thermal conductivity can be appraised k = -ϕ z ∆z ∆T . As all these points were already detailed in previous studies, technical issues concerning the sampling and the scattering procedure will not be recalled presently.

The major modification implemented in the MC simulation tool is about the handling of the geometry of constrictions. In the present work continuous function R(z) are used to set the radius of the constricted nanowire as a function of the z coordinate, in the Cartesian coordinates it satisfies

x 2 + y 2 -R 2 (z) = 0.
Using continuous functions makes easier the calculation of cross sections and control volumes, but also the boundary collision procedure. In the latter one, knowing the phonon location, its propagation direction and the shape of the domain one can easily extract from Newton-Raphson algorithm the intersection point where the phonon bundle crosses the external boundary of a constricted nanowire.

(x ini + v gx t I ) 2 + (y ini + v gy t I ) 2 -R 2 (z ini + v gz t I ) = 0 (5)
where t I is the phonon travel time between its initial position (x ini , y ini , z ini ) and the intersection point (x I , y I , z I ) with the external boundary. From this point, the normal direction n I at the surface is computed and the phonon is assumed to be diffusely reflected (sampling polar and azimuthal angles) in the hemisphere defined by n I .

n xI = -2x I n yI = -2y I (6)
n zI = 2R(z I ) dR(z) dz zI

Studied structures

Two geometries are presented in the paper, the "short" and the "long" constrictions, Figure [START_REF] Cahill | Nanoscale thermal transport[END_REF]. They are respectively defined by equations ( 7) and [START_REF] Ren | Thermal conductivity of superlattices[END_REF].

z 0 z 1 z 2 0 z z 0 d max d max d min d min L R(z) R(z)
R(z) = A 1 -B 1 - B 1 + C(z -z 0 ) 2 (7) R(z) = A + B {arctan[C ⋆ (z -z 2 )] -arctan[C ⋆ (z -z 1 )]} (8)
A and B are constants that can be calculated using geometric conditions;

R(0) = d max /2 and R(z 0 ) = d min /2
for Eq. ( 7), R(0) = d max /2 and R( z1+z2 2 ) = d min /2 for Eq. ( 8). L is the characteristic length of the "long constriction" (distance between z 1 and z 2 ). The parameters C and C ⋆ related to the "steepness" of the constriction, for a large value of C or C ⋆ the constriction is steep, while for a small value of C or C ⋆ it is smooth. In the case of a short constriction the radius as a function of z becomes

R(z) = d max 2 1 - 1 -R d 1 + C(z -z 0 ) 2 (9)
where R d = d min /d max stands for the diameter ratio.

Results and discussions

In this section phonon transport in short constrictions is addressed. Thermal conductivity of constricted nanowires and conductance of the constriction are studied as a function of: diameter ratio R d , phonon mfp Λ, constriction steepness C and temperature T . We divided our results in three subsections. along the nanowire length. First, without any constriction, the TC is really close to the experimental measurement [START_REF] Li | Thermal conductivity of individual silicon nanowires[END_REF], the heat flux is maximum and the temperature profile almost linear as the phonon transport in the structure is diffusive (except in first and last cells where ballistic exchanges with boundaries, at fixed temperatures, can occur). For the constricted nanowires, temperature profiles exhibit drops in the vicinity of the constriction. The latter one seems to be more continuous in the case of a smooth constriction while it is abrupt for the steep one. As a consequence, heat flux density increases at the constriction location.

The magnitude of the density heat flux peak shall be considered with caution as it depends on the spatial mesh choice. In the present calculations there are 40 control volumes with the same value of V to ensure reliable statistics in the MC procedure. It means that in the case of a steep constriction, no control volume is exactly inside the constriction, but rather on each side. On nanowires, it obviously decreases. This decrease is more important for the smooth one. This behavior can be explained from the phonon boundary scattering mfp calculation. When the constriction is smooth, phonon mfp for boundary scattering is getting smaller and thus the thermal conductivity. This point is addressed in the following.

Thermal conductivity

Steepness dependence

First, the dependence of the TC with respect to the phonon mfp in constricted nanowire is addressed through the steepness of the constriction. The latter parameter can be obtained numerically and analytically.

Using the MC simulation tool, a ray tracing (RT) experiment was carried out. For such a case, a very large amount of phonons is sampled in the whole structure and collisions with the nanowire surface envelope are counted during one small time step. The ratio of the number of collision N coll to the total number of drawn phonon N tot is equal to the collision probability. Thus, using the same technique as in Lacroix et al [START_REF] Lacroix | Monte Carlo simulation of phonon confinement in silicon nanostructures: Application to the determination of the thermal conductivity of silicon nanowires[END_REF], one can appraise Λ RT from numerical data as

Λ RT (ω, p, d, C) = v g ∆t ln Ntot(ω,p) Ntot(ω,p)-N coll (ω,p) (10)
where ω stands for the frequency and p the polarization. Simulation results show that there is no frequency or polarization dependence, but the mfp naturally varies according to the shape of the nanowire. Another way to express a particle mfp subject to diffuse transport was proposed by Blanco and Fournier [START_REF] Blanco | An invariance property of diffusive random walks[END_REF]. This analytical expression solely depends on the volume of the structure V tot and the surface S env of the contour envelope and it reads

Λ th (ω, p, d, C) = 4V tot S env (11)
Both expressions were used to calculate Λ for the previously studied nanowires using different values of the steepness parameter C. Results are reported in Table [START_REF] Cahill | Nanoscale thermal transport[END_REF]. The calculated mfp by both methods are in quite Table 1: Variation of the phonon mfp for boundary scattering as the function of the steepness parameter C using Eq. ( 10) and Eq. ( 11); d max =115 nm, R d = 0.5; Λ is given in nm.

C 10 13 10 14 5 10 good agreement. In the case of the analytic solution, for large values of C, the constriction is a very localized modification of the geometry (aperture of diameter d min in the cylinder of diameter d max ) and thus the mfp tends to d max . For smooth constrictions (i.e. small values of C), the mfp becomes smaller as the geometry looks like a "conical" constriction. On figure [START_REF] Han | Thermal conductivity of carbon nanotubes and their polymer nanocomposites: A review[END_REF], thermal conductivity results are reported as a function of Λ. Using this description obviously shows that TC is linearly varying with the mfp. Furthermore, for the largest value of C, where the constriction is similar to a single aperture in the nanowire, and for R d =0.5, the TC is reduced by 10% as compared to the cylindrical nanowire.

Diameter ratio dependence

On the basis of the previously adopted methodology, the TC is calculated for several nanowire diameters and constriction ratios R d . In Figure [START_REF] Maldovan | Thermal energy transport model for macro-to-nanograin polycrystalline semiconductors[END_REF] Moreover, it can be observed that the steepness of the constriction affects the calculated TC whatever the nanowire diameter is. This is due to the reduction of the phonon mfp for smooth constrictions discussed previously. The second point is that the constriction ratio R d lowering also induces the lowering of TC. This effect is more important for nanowires with large diameters for which phonon mfp is more affected by the constriction. Besides, it can be noted that there is a strong decrease of the TC when R d is small (here we have R d =0.1). In the latter case, the TC remains almost constant for all the studied diameters.

For such small values of R d , transport is ballistic through the aperture. Here, the phonons transmitted through the constriction are propagating along z axis. This point will be studied in more details in the next subsection dealing with the conductance. It shall be noted that simulations done for small d max and R d values, are in the limit of the BTE resolution with our technique since dispersion properties in the vicinity of the constriction are likely to be modified. Besides, several system lengths were simulated (here L z =2 µm) and no difference was found in the calculated TC. Here, the characteristic mfp is of the order of the nanowire diameter, that is much lower than its length, and thus thermal transport is diffusive.

Temperature dependence

Finally, the temperature dependence of TC in constricted nanowires is discussed. In Figure [START_REF] Wang | Thermal Conductivity of Nanocrystalline Silicon: Importance of Grain Size and Frequency-Dependent Mean Free Paths[END_REF] is plotted the TC versus temperature for a nanowire with d max =115 nm, C = 10 16 m -2 and two constriction ratios (0.5 and 0.3) over the temperature range 30 < T < 600 K. The temperature dependence for all the studied [START_REF] Li | Thermal conductivity of individual silicon nanowires[END_REF] geometries (with and without constriction) follows the same trends as in bulk Si, with a noticeable decrease of the TC due to boundary scattering. In the case of the nanowire with no constriction, simulation data are in quite good agreement with the one of Li et al [START_REF] Li | Thermal conductivity of individual silicon nanowires[END_REF]. Once again, TC of constricted nanowire is always lower than the one of regular nanowire. Nevertheless, the decrease is more important for temperatures between 50 K and 400 K. Above this temperature the phonon scattering due to phonon-phonon interactions is no more negligible as compared to the one of boundary scattering and a competition between both mechanisms occurs. Such dependence to the temperature paves the way to possible thermal rectification by adjusting the shape of the constriction and the temperature of thermal baths on both sides of the nanowire. For temperatures below 30 K, MC calculations are not carried out since phonon wavelength becomes larger and closer to the constriction aperture leading to potential diffraction effects.

Constriction conductance

In order to address more precisely the impact of the constriction itself on the phonon transport, the thermal conductance G = Φ ∆T is calculated using a simple model of thermal resistance in series. We suppose that the total thermal resistance of the constricted nanowire is the sum of two resistances R tot = R nw + R C , one for the regular nanowire plus one corresponding to the constriction (they are respectively the analogs of frictional and singular head losses in fluid mechanics). In terms of conductance it reads

G C = G tot × G nw G nw -G tot (12)
Using the Eq. ( 12) does not require to define the temperature gradient and the heat flux in the constrictions;

the exact location of the constriction boundaries being subject to different interpretation in the case of smooth and steep profiles, as it is shown in Figure [START_REF] Luo | Nanoscale heat transfer -from computation to experiment[END_REF].

In Figure [START_REF] Ren | Thermal conductivity of superlattices[END_REF] are reported the thermal conductance of the constriction for different values of the nanowire diameter and of the diameter ratio. Several points can be observed. First concerning the steepness of the constriction, it can be noted that in the case of the smooth one, there is no visible agreement with the theory for diffusive and ballistic limits. This is logical since Maxwell (diffusive) and Sharvin (ballistic) models are related to a point contact constriction (single aperture between two structures). This is not the case for the nanowire with C = 10 14 m -2 where the shape of constriction is hyperboloid like. In the case of a steep constriction, MC calculations of the conductance vary between both limits as the simulated geometry is similar to the one of the point contact model. Secondly, for a small diameter ratio (R d =0.1), the conductance trend is obviously close to the one in the ballistic limit. It is in agreement with remarks made in the previous section about the TC (Figure ( 6)). As the diameter of the nanowire gets larger, as well as the constriction ratio, the constriction conductance reaches the diffusive limit.

Long constrictions

We now discuss the case of a long constriction, defined as a nanowire with diameter m -2 ; comparisons to the ballistic limit Eq. (2) and diffusive the limit Eq. ( 1) sides of the constriction the length of d max nanowires is equal to 500 nm. The first observation concerns the length of the constriction. When L is larger than 500 nm, there is almost no change in the calculated mean TC of the nanostructure, whatever the diameter ratio R d is. The second important point is that TC of long constriction systems is always smaller than those corresponding to short constrictions with the same other geometrical characteristics (green diamonds in figure [START_REF] Yang | Partially coherent phonon heat conduction in superlattices[END_REF]). It can be recalled that in the present work, phonon transport is diffusive and consequently, difference in TC is due to the new internal thermal resistance. This is confirmed by the comparison of TC values of long constrictions with the pristine nanowire ones. For each [START_REF] Li | Thermal conductivity of individual silicon nanowires[END_REF] and attributed to nanowire roughness [START_REF] Chen | Thermal Conductance of Thin Silicon Nanowires[END_REF]. Therefore, the presence of heat reservoirs on both sides of the inner nanowire with the small diameter seems to affect the overall thermal conductivity and should be accounted for in thermal properties measurements at nanoscales. Such behaviour is also confirmed by calculations achieved on modulated nanowires [START_REF] Termentzidis | Modulated SiC nanowires: Molecular dynamics study of their thermal properties[END_REF][START_REF] Zianni | Scaling behavior of the thermal conductivity of width-modulated nanowires and nanofilms for heat transfer control at the nanoscale[END_REF], where TC of a structure that alternates large and small diameters is lower than the one of the corresponding small diameter.

Conclusions

Monte Carlo simulations of phonon transport in nanowires with a constriction were performed and discussed. Thermal conductivity as well as thermal conductance have been appraised and compared to available experimental data and literature models. A detailed study of the constriction geometry was carried out. Steepness of the constriction and the magnitude of diameter reduction were first addressed.

It has been shown that the TC behavior is ruled by the phonon mean free path related to boundary scattering. The latter one can be evaluated analytically. Concerning the magnitude of the constriction, two regimes exist. When the diameter of constriction is large, there is little resistance to phonon transport and the regime is "diffusive". Oppositely, when the diameter of constriction is small, thermal conductivity remains quite uniform and the transport is "ballistic". This behavior is also observed for the thermal conductance that varies between these two limit regimes when the constriction is steep. In the latter case, the constriction fit to the "point contact" model and the associated conductance can be simply appraised from the thermal conductivity of the nanowire. In such a case conductance is a linear function of the Knudsen number and of the thermal conductivity of the bounding material as suggested in former studies like [START_REF] Volz | Increase of thermal resistance between a nanostructure and a surface due to phonon multireflections[END_REF]. Otherwise, we show that in the case of smooth constriction, the conductance is no longer described by the Maxwell or Sharvin models.The temperature dependence was also investigated, as expected, the impact of the constriction on the TC is more or less pronounced as the the phonon mfp is dominated by boundary scattering or three phonons processes. TC being significantly reduced by constriction occurrence at low temperatures that at high. The last part of this work is about "long constrictions" that can be related to a nanowire between two heat reservoirs. For such structure, the TC is lower than in the case of a single nanowire. It puts forward the importance of thermal resistances induced by the constrictions. This behavior confirms previous studies on modulated nanowire and gives useful insights to understand thermal properties measurements experimentally achieved. This work will be continued considering modulations and asymmetric constrictions that could exhibit thermal rectification. Moreover, phonon transport within constrictions, and especially frequency and wavelength dependence will be deeper investigated.

Acknowledgements

People involved in this work acknowledge invited professorship of the University of Lorraine (LEMTA) for Pr X. Zianni and the project "SiMoNE", PEPS Mirabelle 2013 funded by CNRS and the University of Lorraine.

Figure 1 :

 1 Figure 1: "Short" (blue) and "long" (red) constrictions; description of characteristic lengths and coordinates

3. 1 .Figure 2 :

 12 Figure 2: Different shapes of nanowire constriction; L z =2 µm, d max =115 nm, d min =57.5 nm; top C = 10 13 m -2 , middle C = 10 14 m -2 and bottom C = 10 16 m -2

Figure ( 4 )Figure 3 :

 43 Figure 3: Heat flux density and temperature profiles, d max =115 nm, R d = 0.5; no constriction (black), smooth constriction C = 10 14 m -2 (red), steep constriction C = 10 16 m -2 (blue)

Figure 4 :

 4 Figure 4: Heat flux density as a function of the number of control volumes, d max =115 nm, R d = 0.5, C = 10 14 m -2 ; 40 cells (black), 180 cells (red)

Figure 5 :

 5 Figure 5: Thermal conductivity as a function of the mfp, d max =115 nm, R d = 0.5

  are reported the MC simulation results for six values of d max and three constriction ratio. Calculations are compared to the TC of a nanowire without constriction. First, the

Figure 6 :

 6 Figure 6: Thermal conductivity as a function of the nanowire diameter and the constriction ratio; d max =[37; 42; 56; 65; 86; 115] nm, R d =[0.1; 0.3; 0.5; 1.0]

Figure 7 :

 7 Figure 7: Thermal conductivity as a function of the temperature; d max =115 nm, C = 10 16 m -2 , R d =[0.3;0.5; 1.0]; yellow squares stand for experimental data taken from Li et al[START_REF] Li | Thermal conductivity of individual silicon nanowires[END_REF]

Figure 8 :

 8 Figure 8: Constriction thermal conductance as a function of the nanowire diameter and the constriction ratio; d max =[37; 42; 56; 65; 86; 115] nm, R d =[0.1; 0.3; 0.5]; black squares C = 10 16 m -2, open squares C = 10 14 m -2 ; comparisons to the ballistic limit Eq. (2) and diffusive the limit Eq.[START_REF] Cahill | Nanoscale thermal transport[END_REF]

RFigure 9 :

 9 Figure 9: Long constriction thermal conductivity as a function of the constriction ratio and the length of the constriction; d max =115 nm, R d =[0.25; 0.50; 0.75; 1.00], L=[500; 1000; 1500; 2000] nm; green diamonds steep constriction C ⋆ = 10 8 m -1 ; yellow hexagons pristine nanowires

  14 10 15 5 10 15 10 16 5 10 16 10 17

	Λ RT	96.6 106.4 108.3 110.1 110.6 110.9 110.9 111.0
	Λ th [34] 95.2 104.5 109.3 110.8 113.1 113.7 114.5 114.7

  14


Sébastien Roujol
email: sebastien.roujol@imf.u-bordeaux2.fr

Baudouin Denis De Senneville
email: baudouin@imf.u-bordeaux2.fr

Silke Hey
email: silke.hey@imf.u-bordeaux2.fr

Chrit Moonen
email: chrit.moonen@imf.u-bordeaux2.fr

Mario Ries
email: mario.ries@imf.u-bordeaux2.fr

M Roujol

Moonen Denis De Senneville

Robust Adaptive Extended Kalman Filtering for Real Time MR-Thermometry Guided HIFU Interventions

Keywords: Biomedical signal processing, Kalman filters, Real time systems, Magnetic resonance imaging. I

published or not. The documents may come   L'archive ouverte pluridisciplinaire

adapt the position of the focal point and the delivered power based on past measured temperature evolution. In addition, the evolution of the temperature (T (t), with t the time) can be used to obtain the thermal dose (T D) as follows:

T D = t 0 2 T (t)-43 , if T < 43 • C t 0 4 T (t)-43 , if T > 43 • C (1)
which provides an accurate and immediate prediction of tissue necrosis [START_REF] Sapareto | Thermal dose determination in cancer therapy[END_REF], which in turn allows to determine the therapy endpoint.

However, the precision of real time MR-thermometry is generally limited by both the available SNR [START_REF] Conturo | Signal to noise in phase angle reconstruction : Dynamic range extension using phase reference offsets[END_REF] and the influence of physiological motion. Consequently, algorithms used for retroactive control are liable to be biased or rendered unstable by the presence of noise on the temperature maps. Furthermore, the accuracy of the calculation of the thermal dose used for necrosis estimation deteriorates rapidly for low SNRs due to its exponential dependence on the temperature. To overcome these problems, temporal filtering has been proposed as a solution to improve the precision of the temperature maps using an infinite impulse response (IIR) filter [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF]. However, the use of this type of low-pass filters, or alternative designs such as finite impulse response (FIR) filters, introduces in general additional latency, leading to a reduced accuracy and limits the achievable temporal resolution of the observation process.

This can be alleviated by using more complex filter designs, which include physical knowledge of the observed system, such as Kalman filters [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]. Kalman filtering is based on the combination of both measured data and data derived from a forecast based on a physical model. Potocki & Tharp in [START_REF] Potocki | Concurrent hyperthermia estimation schemes based on extended Kalman filtering and reduced-order modelling[END_REF] proposed the bio heat transfer equation (BHTE) model [START_REF] Pennes | Analysis of tissue and arterial blood temperatures in the resting human forearm[END_REF] for this purpose. In their experiment, temperature information was only available at four different locations (obtained from optical fibers). The filter was essentially not designed for noise removal (since only four measurement points were available) but to estimate both blood perfusion and temperature information at unmeasured locations (spatial extrapolation). More recently, Ye et al [START_REF] Ye | Model-based ultrasound temperature visualization during and following HIFU exposure[END_REF] employed Kalman filtering for improved model-based ultrasound temperature visualization. The employed predictor model was based on the construction of isothermal ellipsoids around the heated area. Although the precision of the resulting temperature maps was clearly improved (due to a high confidence placed in the model), the accuracy of the method was not evaluated.

Here, we propose a novel spatio-temporal filter based on the Kalman filtering theory that aims to improve MR-thermometry precision while controlling its accuracy. For this, an extended Kalman filter is employed with the BHTE as a predictive model. However, the combination of the predicted data and measured data is not a simple problem. A too large influence of the model would lead to a high noise reduction, but may introduce severe bias on the output accuracy if the model is not properly configured. With increasing emphasis on measured data, the noise removal will not be efficient, however the filter will tolerate imprecise calibration without giving rise to systematic error. Therefore, the filter has to be tuned in function of the model accuracy which is a priori not known and may vary in time. Therefore a dynamic evaluation of the model accuracy was added to the filtering process in order to adjust in real time the confidence in the model and thus the balance between measured and predicted data. The proposed filter was evaluated on simulated datasets and its feasibility is demonstrated on MR-guided HIFU experiments on a agarose gel phantom and in-vivo on a porcine kidney. It showed significant improvement compared to non-adaptive temporal filter designs without giving rise to additional filtering latency.

II. MATERIAL AND METHODS

A. Temperature modeling using the Bio Heat Transfer Equation (BHTE) model

The BHTE model can be used to predict the temperature T from time t -1 to time t, based on the applied acoustic power P and a priori knowledge of the absorption rate (α), the heat diffusion coefficient (D) and the perfusion value (w) [START_REF] Pennes | Analysis of tissue and arterial blood temperatures in the resting human forearm[END_REF]. Note that this simplified model assumes these coefficients as spatialy and temporally invariant. The BHTE in the voxel of coordinates -→ r = (x, y, z) is defined as follows:

∂ ∂t T ( - → r ,t) = α.P ( - → r ,t) + D.∇T ( - → r ,t) -w.T ( - → r ,t) (2)
where ∇ denotes the lorentzian operator. In the absence of large vessels, perfusion effects are often neglected and thus w is set to 0. The BHTE is generally solved in the Fourier domain since the problem can be turned into a linear differential equation as follows:

∂ ∂t T ( - → k ,t) + (D.k 2 + w). T ( - → k ,t) = α. P ( - → k ,t) (3)
where T , P denotes the Fourier transform of T and P respectively and k denotes the frequency coordinates in the Fourier domain. The solution is computed based on the variation of constant as follows:

T ( - → k ,t) = T ( - → k ,t-1) .e -(D.k 2 +w)t + α. P ( - → k ,t)
1 -e -(D.k 2 +w)t D.k 2 + w (4) Since the BHTE models the physical processes of heat diffusion and absorption, it requires to be applied in 3D space in order to obtain unbiased results. Unfortunately, the available acquisition time in is typically too limited in interventional imaging to obtain full 3D temperature imaging, especially if physiological motion has to be resolved. As a consequence, the BHTE has in practice often to be applied to 2D or severely undersampled 3D datasets. This leads to a systematic underestimation of the heat evacuation and neglects potential heat inflow from adjacent slices, in particular for large ∆t .In order to evaluate the resulting bias on the filtered temperature data, the BHTE was implemented for both 2D and 3D temperature prediction and the results subsequently compared on simulated heating experiments.

Furthermore, although MR-thermometry provides coherent temperature information in areas of high signal (generally with a precision of few degrees), areas with very low signal level generally display randomized temperature values (in the range of several dozens of degrees of erroneous variations). Since the BHTE model prediction is based on preceding filtered temperature maps represented in Fourier space, such areas have to be excluded by manual ROI-based masking to prevent a bias due to undesired high frequency noise.

B. BHTE based extended Kalman filtering

Since the discrete solution of the BHTE represents a nonlinear model for data prediction, the original Kalman filtering theory, which requires a linear predictor, is not directly applicable. Non-linear predictors are addressed by the extendend Kalman filtering (EKF) formalism. An EKF can be seen as a two step process. In a first pass, the filter computes a data prediction (at time t) based on the last filtered data point (obtained at time t -1). In a second pass, the algorithm optimally combines both predicted data and measured data to obtain the final filtered data corresponding to time t. 1) First pass: time update equations: The first pass of the filter is often referred to as time update equations. The temperature prediction T - t at time t and the a priori estimate error covariance P - t are computed as follows:

T - t = f (T t-1 , u t-1 ) P - t = A t P t-1 A T t + Q (5)
where T - t denotes the model prediction at time t, T t-1 the filtered temperature at time t -1, f represents the BHTE model, u t-1 is the control input parameter which in our case corresponds to the HIFU delivered power at time t -1. A t is the Jacobian matrix of partial derivatives of f with respect to the temperature T at (T t-1 , u t-1 ), A T t denotes the transpose of A t , P t-1 is the a posteriori estimate error covariance at time t -1 and Q is the process noise covariance related to the model innacuracy.

2) Second pass: measurement update equations: In a second pass, the Kalman filter combines both predicted T - t and measured T m t data to obtain the final filtered data T t . The combination is weighted with a parameter K t , often referred to as the Kalman gain, while the difference T m t -HT - t is called the innovation S t . During the second pass, K t is first updated based on the new a priori estimate error covariance. Subsequently, the filtered temperature T t can be computed from the weighted innovation S t . The final step updates the a posteriori estimate error covariance P t . The measurement update equations are summarized as follows:

K t = P - t H T (HP - t H T + R) -1 T t = T - t + K t (T m t -HT - t ) P t = (I -K t H)P - t (6)
where H relates the model state to the measurement (here H is the identity since the measurement are directly obtained in state space) and R is the measurement noise covariance that relates to the MR-thermometry precision of the measure.

Note that R can be obtained prior to hyperthermia by evaluating the standard deviation of the measured temperature.

C. Autocalibrated extended Kalman filtering (AEKF)

The EKF combines measured and predicted data via the Kalman gain K t , that allows to optimally adjust the confidence between the employed model and the measured data. The Kalman gain is updated dynamically and predominantly influenced by two input parameters of the filter: The measurement noise covariance (R) and the process noise covariance (Q). While the measurement noise covariance can be determined from baseline data in absence of heating, the determination of Q is not straightforward. Q corresponds to the model accuracy and thus depends on the accurate knowledge of the physical parameters of the tissue (absorption coefficient and heat diffusion), which are a-priory not exactly known and thus in general only available as rough estimates.

In addition, heat absorption only occurs during the heating period when the HIFU system delivers acoustic energy, while heat diffusion is present during both heating and cooling. Therefore, the BHTE model may have a varying performance over time if not properly configured (for example using an incorrect absorption coefficient with the true diffusion coefficient). In such a case, Q would be better chosen as time variant.

This can be achieved using an adaptive EKF (AEKF), where Q is automatically adjusted over time based on a dynamic evaluation of the model accuracy for each new measurement. In the proposed implementation this is based on the assumption that temperature noise is white noise around the real value. The model is considered accurate at time t if ε t , the spatio-temporal sum of the difference between predicted and measured data, is below a predefined error threshold ε threshold , with

ε t = N i=t-N (x,y)∈ ( T - i (x, y) -T m i (x, y)). (7)
Here, the model is evaluated over the temporal window size N and over a voxel perimeter around the focal point, which is in the scope of this paper refered to as the spatial window size M . ε threshold is chosen as the maximum acceptable penality of the filtering process on the measurement accuracy.

D. Robust approach of AEKF

The computation of temperature maps requires an image processing pipeline, that can be rather complex especially in the case of MR-thermometry applied to mobile organs. Therefore, severe artifacts on temperature maps can be observed for a vaiety of reasons such as incorrect phase unwrapping or imperfect motion compensation. This type of error can have amplitudes much larger than the imprecision caused by low SNRs. Since the determination of the thermal dose corresponds to the integral over time of the temperature (see equation 1), temperature artifacts introduce a non reversible error bias on the thermal dose calculation. In addition, due to the exponential dependance on the temperature, even temporally sparse occurances of severe temperature artifacts lead frequently to an apparent thermal dose which is magnitudes off the true value. For HIFU interventions which use the thermal dose as a criteria for the therapy endpoint, this potentially directly affects the success of the intervention (overestimation) or patient safety (underestimation) and is thus highly undesireable. To detect such occurances, an outlier rejection based on the Chauvenet's criterion [START_REF] Chauvenet | A manual of spherical and practical astronomy[END_REF] was applied to the difference between measured and predicted temperature, i.e. the innovation S t . The outlier rejection considers an innovation S t for rejection if the probability to obtain its deviation from the mean St is less than 1/(2 • N S), N S being the number of measurement samples. Substituting St by ε t-1 as defined in equation 7, leads to the following formulation of Chauvenet's criterion:

T t = T - t , if |S t -ε t-1 | > σ(S t-1 ) * e max T - t + K t S t , if |S t -ε t-1 | ≤ σ(S t-1 ) * e max (8)
where e max is the ratio of maximum acceptable deviation to precision [START_REF] Chauvenet | A manual of spherical and practical astronomy[END_REF] and σ(S t-1 ) is the standard deviation of the innovation over the N × M sample points, excluding the last innovation S t . Note that a rejection of a measured temperature value T m t leads the filtered value T t to depend only on the prediction T - t , which is equivalent to temporarily setting R=∞ in equation 6.

E. Experimental set up

1) Simulations: A 3D reference temperature dataset was simulated using the BHTE model (absorption=0.02 KJ -1 , diffusion=0.1 mm 2 s -1 , delivered power=100 W (between 20 th and 70 th dynamics), field of view=32 × 32 × 16, voxel size=1 × 1 × 2 mm 3 , focal point size=1.23 × 1.23 × 7.88 mm 3 , dynamic scan time=1s). N (=100) datasets were derived with added Gaussian noise (σ=5 • C). The proposed AEKF was evaluated using a 3D implementation of the filter. Since in practical cases, MR-temperature information is often available in 2D, the impact of a 2D implementation of the filter was also analyzed.

A matched Kaiser Bessel FIR filter was used for comparison and configured to retain 90% of the integrated power spectrum of the original noise free temperature simulation (pass band cutoff frequency=12Hz, stop-band attenuation>21db, window size=15).

To quantify the performances of each filter, the mean square error (MSE) between filter output and reference data was used as quality criterion of the filter output accuracy.

2) MR-guided HIFU experiments: MR-imaging was performed on a Philips Achieva 1.5 T MRI (Philips, Best, the Netherlands) and HIFU heating was performed using the integrated Philips' Sonalleve MR-HIFU system (Philips, Helsinki, Finland). a) Heating experiment on an agarose gel phantom: Dynamic MR-temperature imaging of an agarose gel phantom was obtained with a gradient recalled single-shot echo-planar sequence (300 images, single slice, TR=106ms, TE=36ms, flip angle=40 • , voxel size=1.7×1.7×3mm 3 , matrix=64×64). HIFU heating protocol employed 50 W of electrical power to the HIFU transducer over a period of 40 s in a fix point P .

Physical filter parameters (absorption and diffusion) were estimated using [START_REF] Oliveira | Non-invasive determination of tissue thermal parameters from high intensity focused ultrasound treatment monitored by volumetric MRI thermometry[END_REF]. Since only 2D temperature maps were available, a 2D implementation of the EKF was employed. Finally, the EKF was compared with a matched Kaiser Bessel FIR.

A second HIFU heating was performed 20 mn later (once heat evacuation of the first experiment was completed) using a modified heating protocol. Here, the focal point position was electronically updated each TR/4 ms in four different locations (± 0.5 mm in each axis from the point P ). The AEKF was then applied, using the BHTE configuration obtained from the first fixed point experiment (which is thus not adapted to the present case), in order to observe its performance with a deliberately mis-configured BHTE model b) Heating experiment on an in-vivo porcine kidney: MRI guided HIFU heating was performed in vivo in the kidney of a pig under general anesthesia. Since the kidney was static during the experiment, no focal point position adjustment was required. Dynamic MR-temperature imaging was realized in the Philips Achieva 1.5 T MRI and employed sequence was configured as follows: 1500 dynamic sagittal images, one slice, TR=127 ms, TE=25 ms, flip angle=35 • , FOV=142.5×285 mm 2 , voxel size=3×3×6 mm 3 . Heating was performed using 250 Watts of electrical power during a period of 19 s. After completion of the experiments the animals were euthanized by intravenous injection of Pentobarbital. The animal experiments were conducted in agreement with the french law on animal experimentation and in compliance with institution's rules for animal care and use.

III. RESULTS

A. Simulations

Figure 1 shows the filter performances on a simulated dataset. The FIR filter (figure 1a) introduces a latency inherent to its design, while the EKF (figure 1b,c) avoids this effect. When configured with the exact physical parameters (figure 1b) the filter does not introduce any bias on accuracy and improves the output precision. An emphasis on the measured data (Q=10) leads to a little smoothing of the temperature. At the contrary, a high confidence in the model (Q=0.1) provides an efficient noise removal in this case. Figure 1c shows the influence of an imperfect configuration of the model (in this case with an 50% overestimated absorption coefficient). Here, a bias on accuracy is observed, especially for Q values giving a large influence to the model. On the other hand, using small values of Q reduces the penalty on the accuracy but limits the filter ability to reduce measurement noise.

The error bias introduced by an approximative configuration of the BHTE based model was then investigated using a 3D implementation of the EKF. With the employed absorption and diffusion coefficients and the chosen spatial resolution of the simulation, the choice of the absorption coefficient has a larger impact than the choice of the diffusion coefficient. Therefore, an illustration of the bias introduced by an approximative configurations of the BHTE based model, obtained with several incorrect aborption values, is reported in figure 2. For this, the filter was tested using different Q values over all N datasets with an absorption coefficient varying for -50% to +50% of its true value. Evaluation of the error bias introduced by a misconfiguration of the absorption coefficient are plotted in figure 2a,b. It can be observed that optimal Q values are different for the heating period (figure 2a,) and for the cooling period (figure 2b) that illustrates the time varying accuracy of the BHTE model. However, for all configurations, even in the case of a mis-configuration of the BHTE model with a bias of ± 50 % of the input absorption coefficient, a good choice of the Q value can guarantee a better filter performance compared to the FIR filter. Therefore, an adaptive tuning of Q during the intervention appears mandatory.

The performances of the presented adaptive EKF was investigated and the obtained filtered temperature with the time optimized Q values are plotted in figure 3. The AEKF used a 3D implementation of the model that was configured with an incorrect absorption coefficient of -50% of its true value and with the true diffusion coefficient. These configuration is thus similar to the EKF configuration employed in figure 1c and results are thus directly comparable. The initial mean square errors of 22.7 • C and 21.9 • C during the heat up and cool down periods were reduced to 8.1 • C and 0.5 • C using the mis configured AEKF. The same configuration of the EKF lead to a MSE of 26.9 • C and 1.2 • C with Q=0.1 and 8.4 • C and 8.5 • C with Q=10 (see figure 2c) for heating and cooling periods respectively. Therefore, AEKF performances are similar to EKF with Q = 10. during heating and similar to EKF with (a) (b) Fig. 3. Effects of an adaptive EKF using a 3D implementation of the BHTE model with a mis configured absorption coefficient (with α = 0.1 which represents -50% of its true value) and the true diffusion coefficient (similar configuration to the EKF in figure 1c). Filtered temperature curve obtained with the AEKF is shown in (a) and the time optimized Q values are reported in (b). Compared to the FIR ourput (figure 1a), the latency was substancially reduced. Then, before heating and during cooling, where the BHTE model was properly configured, the AEKF provided a better noise reduction compared to the FIR or the EKF with Q=10 (figure 1c) However, the AEKF depicted similar performances to the EKF with Q=0.01 (figure 1c, where an emphasis on the model is given). On the contrary, during the heating period, when the BHTE model was incorrectly configured, the mean sqaure error (MSE) obtained with the tested filters (FIR: 18.2 • C, EKF with Q=0.1: 26.9 • C) was reduced to 8 • C with the AEKF that provides comparable results to the EKF with Q=10 (MSE=8.4 • C, where an emphasis on measurement is given). Q = 0.1 during cooling confirming the benefit of an adaption of the Q value. Finally, the AEKF clearly outperforms the FIR that provides a MSE of 18.3 • C and 11.8 • C for each of the two periods (see figure 2a).

Since, the BHTE based model was only incorrect during the heat up period (since absorption is only relevant during power delivering). The time optimized Q values increased during this period in order to give more confidence in the measurements and converged into an interval of [log(0), log(2)] (figure 1b). This corresponds well to the optimum Q values (global minimum) obtained in figure 2a with the blue curve (-50%). Then, the model being true during the cooling down period, the optimized Q values decreased to improve the confidence in the model.

The potential of the proposed AEKF to automatically adapt and optimize the value of Q was tested over the same datasets and the results are plotted in figure 4. The influence of a mis configuration of both absorption (4a,b) and diffusion (4c,d) coefficients was evaluated during the heat-up (4a,c) and cooldown (4b,d) periods. The resulting accuracy obtained with the FIR filter was comparable to the accuracy without filtering, which is mainly due to the latency introduced by the filter and its inability to follow rapid temperature variation. However, during the cool down period, the FIR allows to improve the resulting accuracy by a factor of 3, due to the reduced temperature variation between two successive temperature points (compared to the heat up period). Using the AEKF, the resulting accuracy was reduced by a factor in the worst tested case by a factor 3 and 15 for the heat up and cool down period respectively. Therefore, the AEKF allows to outperform the FIR performances even in the case of a mis-configuration of the BHTE model with a bias of ± 50 % of the input parameters.

The influence of a 2D implementation of the filter was then evaluated over the same simulated datasets as for the 3D AEKF evaluation. Figure 5 shows the impact of a mis Here, the best performances of the AEKF are not achieved with the true configuration of the absorption and diffusion coefficients. Here, the focal point size in the third dimension is 7.88mm which is three time larger that the resolution of the simulation in the third dimension. In these conditions, all energy supposed to be delivered above and below the 2D single slice is not considered with this 2D implementation. Therefore, an AEKF with either an absorption coefficient over estimated of +25% of the true value or a diffusion coefficient under estimated of -25% of its true value allows to compensate for the lack of received energy and in this case outperform an AEKF with a true configuration. Similarly to the 3D AEKF evaluation, the resulting accuracy was reduced by a factor in the worst tested case by a factor 3 and 15 for the heat up and cool down period respectively. Therefore, the AEKF allows to outperform the FIR performances even The method was evaluated on an agarose gel phantom during a HIFU heating experiment (see figure 6). The AEKF provided a better reduction of the temperature noise and did not introduce large latency as obtained with the FIR (see figure 6a). A second heating experiment was then realized with a different heating protocol where the focal point position was electronically updated in four different locations. The AEKF was run using the same configuration as in the fixed point experiment, which is thus expected to be incorrect for this experiment. Results are reported in figure 6b and shows that the AEKF considers the model as incorrect during the heating period where the filter has a low confidence in the model. During the cooling period, the model is correct since only the heat diffusion is present (and is properly configured) and the filter stronly improves its confidence in the BHTE based model that result in a substancial reduction of the temperature noise. Again, the bias introduced by the latency of the FIR is also avoided with the AEKF. An in-vivo heating experiment was finally performed in a porcine kidney. A temperature elevation of 18 • C was reached in the focal point. A large temperature artifact (45 • C) was simulated at the 288 th dynamic in the focal point. This type of artifact induced a strong perturbartion of both the FIR filter (figure 7a) and the AEKF (figure 7b). The robust AEKF approach allowed to detect this artifacted datapoint and the resulting filtered temperature was not affected (figure 7c).

D. Real time benchmarking

The computation time of the EKF (using a dual processor, INTEL 3.1 GHz Penryn, two cores, 8 GB of RAM) was 0.71 ms for a 2D image of resolution 32×32 and 11.4 ms for a 3D image of resolution 32×32×16. In the AEKF, the additional time required for the autocalibration of the Q value is defined by the number of steps toward the optimum Q value. In our approach, an upper bound was imposed to the algorithm and was defined in order to ensure the real time constraint available in any case. For the case, of HIFU experiment where a TR of 100 ms was employed with a slide window of 10 images, the maximum number of iteration toward convergence was empirically fixed to 12 (since 12 × 10 × 0.71ms = 85.2ms < T R) to guarantee the real time condition.

IV. DISCUSSION

A. BHTE model based temperature prediction

The performance of the BHTE model relies primarily on the accurate representation of the physical heating process. For this, it requires a priori knowledge of the exact form of the focal point P ( - → r ,t) and the local absorption and diffusion coefficient α and D. In particular, the choice of the later two coefficients is in practice problematic, since it either requires empirical data or a suitable calibration experiment. This is complicated by the fact that the spatial resolution of typical MR-sequences for interventional imaging are of the same order than the observed diffusion phenomena and focal point sizes, which leads to partial volume effects [START_REF] Todd | The effects of spatial sampling choices on mr temperature measurements[END_REF]. Therefore, a calibration results only in apparent absorption and diffusion coefficients, which will vary with increasing voxel sizes. Furthermore, although both apparent coefficients can be expected to be sufficiently homogeneous within large organs such as the liver, this assumption could break down in smaller heterogeneous organs such as the kidney, or close to organ boundaries, where both become spatially variant. Another limitation of the predictor model arises from the negligence of the heat evacuation due to tissue perfusion. Although this is for the presented experiment on a gel of no consequence, future work on in-vivo tissue requires a careful assessment of the corresponding error bias and potentially to take the perfusion term in equation 4 into account.

In addition, since the achievable MR-acquisition time severely limits the available volume coverage, the BHTE has to be in practice applied to 2D or severely undersampled 3D datasets. As shown in figure 5 the influence of a 2D implementation of the model is of little consequence as long as the spatial resolution in slice direction remains much larger than the characteristic diffusion length at the given temporal resolution of the MR-sequence.

Finally, the temperature prediction model is also directly influenced by the choice of the size (in mm) of the prediction area. A large area could prevent the real time use of the method (especially for the autocalibration that requires several temperature prediction to adjust the Q parameter). On the other hand, a small area may be insufficient to encompass the whole heated area introducing undesirable high frequency in the Fourier domain that could be further reduced using apodization function. In the presented gel and in-vivo experiments under typical condition of clinical ablation, an area of 32×32 pixels with a pixel size in the range of a milimeter was found to be a good compromise.

B. Autocalibration performances of the EKF, Dynamic evaluation of the BHTE model accuracy

Although the BHTE model is an efficient way to predict the temperature, the uncertainty of the model is a priori not known and was shown to vary over time (especially between heating up and cooling down period). The autocalibration of the Q parameter allows to correctly handle the model accuracy in the proposed filter. Nevertheless, the autocalibration process is influenced by two parameters: the spatial window size (M ) and the temporal window size (N ) employed to evaluate the filter accuracy (see in equation 7).

The spatial window size has to cover voxels contained in the ablation area and thus subject to the influence of the heat absorption and the heat diffusion. In our experiment, 9 voxels around the focal point were selected for 2D experiments and 9 voxels of each of the two adjacent slices were included for 3D experiments (leading to a total of 27 voxels). Although the number of voxels prone to temperature increase can be limited in a fixed point heating experiment, this issue is expected to be clearly reduced for volumetric ablation (leading to a larger ablated area).

The temporal window size directly influences the frequency response of the proposed adaptation of Q. A large temporal window leads to stable convergence of the Q value, while introducing a latency on the Q adaptation. A small temporal window allows reactive adaptation of Q at the price of a reduced statistical sample contributing to equation 7. Consequently, this parameter has to be adjusted depending on the temporal dynamic of the temperature between two successive acquisitions. The auto-adaptation of Q allows to dynamically find the smallest Q value that maintains the introduced error bias of temperature below a predefined threshold defined by the ε parameter. This allows maximal measurement noise filtering under the condition of a guaranteed accuracy. For the case of a severely mis-configured filter, the auto-adaptation will gives a higher emphasis of the measured data (high value of Q) in order to limit the introduced accuracy bias. In practice this means the filter will be less efficient for noise removal, however, a severely mis-configured filter will not introduce systematic errors in the filtered data. ited to an upper bound defined by the ε parameter.

Contrary to autocalibration approaches directly based on the analysis of the innovation sequence such as in [START_REF] Mehra | On the identification of variances and adaptive kalman filtering[END_REF], the proposed adaptation allows a further reduction of the temperature noise at an accuracy cost, which is bounded by the predefined maximal accepted penatly of acuracy ε threshold .

C. Robust AEKF approach

In general, artifact detection is complicated by the requirement to find a robust criterion allowing to differentiate between a measurement artifact and changes due to the dynamic of the measured process. Since the proposed AEKF formalism employes the BHTE as a physical model for data prediction of the heating process, it is straightforward to seamlessly integrate a robust outlier detection without significant additional computational overhead. This is achieved by examining the innovation S t for unphysical changes and to replace unphysical measuements by the predictor value T - t . Therefore, contrary to temporal filtering, where the amplitude of the artifacted datapoint still influences the amplitude of the filtered data (figure 7 a,b), the proposed approach allows to replace the value entirely with a prediction based on the physical heating model (figure 7 c). Furthermore, the proposed AEKF filter also accounts for the variations of the average model error during the heating process and thus conviently allows to adjust the outlier rejection criterion by using the updated ε(t -1) and σ(S t-1 ) values.

D. Clinical aspects of temporal filtering of the temperature

The presented data has shown that due to the use of a physical model of the heating process, the AEKFs model based filtering performance surpasses the results of more simple filter designs. However, in order to evaluate the usefullness of the proposed AEKF for clinical applications, two other categories of requirements have to be considered: Patient safety aspects and real-world practicability. Since the proposed filter is intended for real-time MR-guidance of non-invasive HIFU ablations, additional latency and accuracy degradation can directly affect the success of the intervention and ultimately patient safety. Therefore, despite the fact that the BHTE model parameters are prior to the intervention only available as coarse estimates, the AEKF must be able to cope with a severe misconfiguration without impairing the accuracy of the filtered temperature values, as shown in figure 1 c for the EKF. The AEKF filtering process was designed so that even a worst case scenario merely results in an ineffective noise removal, without systematic errors introduced by the filter as shown in figure 3 a and 6 b. Furthermore, reliable thermal dose measurements are frequently hampered by the effects of occasional spontaneous motion events or instabilities/errors of the employed realtime data processing. The proposed integrated outlier detection reliably identifies and removes such temperature artifacts and thus allows to perfom sustained dosimetry over extended durations of several minutes. Since the artifacted datapoint is entirely replaced by a model based estimate, a systematic error bias of the final thermal dose value is prevented.

The proposed filter design is more complex than nonadaptive FIR or IIR filter designs, as a consequence the second important aspect is the practicality of it's use. The main disadvantage of the proposed BHTE based AEKF is that the filter requires an approximate estimate of the absorption coefficient α and the diffusion coefficient D prior to the intervention. Both can be provided on an individual basis based on a low-power test shot as shown in [START_REF] Oliveira | Non-invasive determination of tissue thermal parameters from high intensity focused ultrasound treatment monitored by volumetric MRI thermometry[END_REF]. Alternatively, since an adequate filter performance does not require the exact knowledge of either parameter, for standardised clinical scenarios such as uterine fibroid ablations, both values could be supplied based on averaged patient data. All other parameters of the filter are autocalibrated during the intervention, which reduces the risk of an arbitrary/empirical calibation by the user. Finally, although the AEKF implementation has a higher computational overhead than FIR or IIR designs, the benchmarking results show that the design is entirely compatible with the requirement of interventional guidance of high-framerate imaging associated with low processing latencies.

V. CONCLUSION

The proposed autocalibrated extended Kalman filter based on the bio-heat transfer equation was demonstrated to improve both precision and accuracy of MR-thermometry compared to simpler filter designs such as FIR-filters, without introducing undesired latency. However, the increased complexity of a model based filter design also increases the risk of systematic errors for the case of severe mis-configuration. For clinical applications in the field of therapy guidance this would represent an unacceptable risk, due to the fact that the filtered data is used for direct retroactive feedback control, and concequently errors can not be undone retrospectively. Since this directly influences patient safety, the proposed autocalibrated extended Kalman filtering formalism was used to control and limit the impact of the filter on latency and accuracy during the entire intervention. In conjunction to the additional outlier rejection algorithm, which replaces sparse temperature measurement artifacts with BHTE prediction based estimates, this represents a promising approach to improve the precision and the feasability of extended MR-thermometry and dosimetry for MR-guidance of HIFU ablations under clinical conditions.

Fig. 1 .

 1 Fig. 1. Examples of temporal filtering obtained in the focal point of the simulated dataset: results obtained using (a) a FIR filter (a) an EKF filter configured with the true absorption and diffusion values (α = 0.2 and D =1e-7), (c) an EKF filter badly configured (wrong absorption (α =0.1) with true diffusion D =1e-7).

Fig. 2 .

 2 Fig. 2. Mean squared error of the filtered signal obtained with a 3D implementation of the EKF filter in the focal point of the simulated dataset. Influence of absorption on the mean square error (MSE) of the filtered temperature curve during heat-up (a) and cool-down (b) periods.

Fig. 4 .

 4 Fig. 4. Temporal filtering performance obtain using several configurations of the model using a 3D implementation of the filter. Influence of absorption (a,b) and diffusion (c,d) on the mean square error (MSE) of the filtered temperature curve during heat-up (a,c) and cool-down (b,d). The AEKF achieves a better filter performance than the FIR even if the tissue absorption α or the thermal diffusion are deliberately miss-configured by ±50%. For each method, the MSE average and standard deviation over the N datasets are plotted in the empty box and the dashed box respectively.

Fig. 5 .

 5 Fig. 5. Influence of a 2D implementation of the filter: Similarly to figure 4 the impact of an incorrect absorption (a,b) and an incorrect diffusion (c,d) was evaluated during both heat up (a,c) and cool donw (b,d) periods.Again, the AEKF achieves a better performance than the FIR even with a 2D implementation of the filter and when the tissue absorption α or the thermal diffusion are deliberately miss-configured by ±50%. However, an over estimation of the absorption (by +25% of the true value) or an under estimation of the diffusion coefficient (by -25% of the true value) that help to compensate for the neglected inflow from adjacent slices lead to improved performances compared to a true configuration .

Fig. 6 .

 6 Fig. 6. Filtering performances of a FIR and the AEKF obtained on an agarose gel during heating experiments. (a) temperature evolution obtained during the first heating experiment in the focal point using a properly configured BHTE based model. (b) temperature evolution in the same point during a the second heating experiment where four different locations where iteratively heated. The AEKF follows measurements during heat up period since the model is incorrect but gains confidence during cool down period where the BHTE based model is correct.

CFig. 7 .

 7 Fig. 7. Temperature evolution of the HIFU experiment on an in-vivo porcine kidney using a (a) FIR filter, (b) the AEKF, (c) the robust AEKF. The robust AEKF allows to detect this strong temperature artifact as outlier and to replace it by the BHTE model prediction value.

This work was supported in part by Ligue Nationale Contre le Cancer, Conseil Régional d'Aquitaine, Diagnostic Molecular Imaging, Agence National de Recherche, Fondation InNaBioSanté and Philips Medical System. Asterisk indicates the corresponding author.