Emmanuel Bertrand
email: emmanuel.bertrand@univ-nantes.fr

Philippe Castany

Thierry Gloriant

An alternative way to orient the parent phase in the cubic/orthorhombic martensitic transformation of titanium shape memory alloys

Keywords: EBSD, Phase reconstruction, Titanium, Shape memory alloys

Titanium-based shape memory alloys undergo a cubic to orthorhombic martensitic transformation (β to α'' martensite). A new easy way operating as a reconstruction of the parent microstructure is used from electron backscattered diffraction. It is shown that the electron backscattered patterns of the martensitic microstructure can be directly indexed as the parent phase in order to obtain the correct orientation of the high temperature microstructure. This method is however specific to the cubic/orthorhombic martensitic transformation in titanium shape memory alloys.

Electron BackScattered Diffraction (EBSD) measurements are commonly used in order to determine the Orientation Relationships (OR) between the parent phase and the product phase that exists in displacive transformations and enables to determine the crystallographic orientation of the high temperature phase before transformation. This technique which is named phase reconstruction was already used for the α/β transformation of titanium and zirconium alloys [1][2][3][4] or for the martensitic transformation in steels [5][6][7].

In titanium alloys, the OR between the stable low temperature α phase and high temperature β phase is known as Burgers relationship. However, the application of the OR on the α phase (hexagonal) is not sufficient to determine the orientation of the β (body-centered-cubic) parent phase because there are 6 equivalent ways to apply the OR on the α phase leading to 6 possible different β orientations [1,3]. It is thus not possible to determine exactly the orientation of the former β phase with this technique and the real orientation has to be determined at high temperature, on residual β-phase existing at room temperature or by other advanced techniques [4]. In Ti-based alloys, it also exists the metastable martensitic orthorhombic α'' phase, which can be obtained by quenching from the high temperature βphase domain. Indeed, in the case of metastable β titanium alloys, a sufficient amount of βstabilizing elements (such as Ta, Nb, Mo, V…) allows the reversible martensitic β/α'' transformation [8][9][10][11]. A major application of titanium alloys showing such a displacive transformation is the production of shape-memory or superelastic nickel-free biocompatible alloys for medical applications [10,[12][13][14] The martensitic α'' phase has a C-centered orthorhombic structure and the β/α'' OR is described in literature as a Burgers type except that the shear and atomic shuffle from the β to the hexagonal α phase do not go to completion [15]. The Burgers OR becomes thus for this β/α" transformation [9,16]: (-110)β // (001)α '' and

[111]β // [110]α '' . An illustration of the correspondence between structures of the β phase and the α'' martensite is given in supplementary material. By application of the OR, each β grain is decomposed during the quench into several possible α'' variants by application of 6 possible lattice correspondences, leading to a α'' selfaccommodating room temperature microstructure [17] and a shape memory effect [8,18].

The shape memory effect comes from the fact that there is one unique path to do the inverse transformation upon heating. From the crystallographic point of view, this means that each of the six α" variant can only transform into a same and unique parent crystal. Hence, the crystallographic reconstruction of the parent phase in an α'' microstructure is expected to be simpler than for other martensitic transformations, for which several paths to do the inverse transformations are possible [1,3]. A classical way to reconstruct the parent phase would be to apply the OR on each lath of α'' martensite. However, the microstructure in such shape memory alloys in usually very fine with nano-sized laths [15,17] leading to difficulties (or impossibility in some cases) to index individual lath by EBSD and thus reconstruct the orientation of the parent phase from martensite indexation. The objective of this paper is to show, in the specific case of the cubic to orthorhombic martensitic transformation, the possibility to directly determine the crystallographic orientation of the high temperature parent phase from the room temperature martensitic microstructure by indexing Electron BackScatter Patterns (EBSPs) without any further data processing and without indexation of the martensite.

The composition of the alloy for this study is Ti-25Ta-20Nb (mass%). It underwent an annealing at 850°C for 0.5 h, followed by water quenching in order to obtain the α''-phase microstructure. Preliminary experiments showed that this alloy exhibits a shape memory effect with austenite start and austenite finish temperatures of A s =205°C and A f = 232°C, respectively. The phases repartition is determined by X-ray diffraction (XRD) using a Bruker D8 Advance diffractometer with LynxEye detector and Cu Kα radiation. EBSD studies are carried out in a Merlin scanning electron microscope (SEM) from Carl Zeiss equipped with a Nordlys 2 camera. An accelerating voltage of 20 kV and a probe current of 3000 pA are used.

The samples were prepared by mechanical polishing to a "mirror finished" state. Chemical etching with a solution of 5% HF, 5% HNO 3 and 90% H 2 O (vol.%) was used to remove the residual deformation layer induced by polishing.

Fig. 1 shows EBSD maps of the water quenched Ti-25Ta-20Nb sample and the X-rays diffractogram. The diffractogram (Fig. 1d) shows that the sample is only composed of α'' martensitic phase which lattice parameters, known to be dependent on the chemical composition [9], are measured as: a=0.3174 nm, b=0,4826 nm and c=0,4638 nm. The Band Contrast (BC) map (Fig. 1a) is used to reveal the microstructure. Hence, grain boundaries and martensite/martensite interfaces have a lower crystallinity than the rest of the microstructure and appear darker. One can observe large domains delimited by dark lines corresponding to the lowest BC, which seem to correspond to grain boundaries of former β grains present at high temperature: these dark lines are thus a signature of the high temperature β microstructure. These ex-β grains will be called grains further in this paper. In these grains, a microstructure made of laths and needles is also visible. This thinner microstructure corresponds to a mix of the different possible variants of α". The BC map confirms the effect of the thermal treatment: the β microstructure at high temperature is decomposed into α'' martensite when the alloy is quenched. Fig. 1b shows the orientations represented by Euler angles (see color online) when the EBSPs are indexed using the α'' phase. The largest laths are unambiguously oriented while the orientation of some domains with a thinner microstructure is more imprecise. Indeed, the indexation rate is close to 60 % which is mainly due to a nanometric sized α'' microstructure in some domains [15,17]. However, five different α'' orientations are unambiguously detected in the largest laths in the grain at the center of the figure. As it will be shown further, all these five α" variants are formed from the same β grain at high temperature (delimited by the dark lines on Fig. 1a).

Usually, reconstruction of the parent phase orientation should be made by applying the Burgers OR on the orientation of each α" variant found previously, needing additional calculations. Another possible way to determine the orientation of the parent β phase would be to index the residual β phase by EBSD, but it is impossible in the present case since there is no β phase in the microstructure (see XRD result in Fig. 1d). In order to determine easily and directly the orientation of the β parent phase, the idea is now to force the indexation of the same EBSPs acquired previously with the β phase, even if, of course, there is no β phase in the microstructure. If the choice between the α" and β phases was let possible, the EBSD software would select automatically the α" phase to index the EBSPs due to a better matching than with the β phase. That is why only the β phase is proposed in the EBSD software as possible phase in order to force the indexation with the parent phase. However, all EBSPs cannot be indexed with the β phase, but relatively few can be indexed because of the similarity of both phases. As a consequence, the indexation rate drops to 25 %. Even if this value is dramatically low, there is no area wherein no pixel is indexed: indexed pixels are thus uniformly distributed in each ex-β grain (even when the α'' microstructure is very thin) making a noise reduction by neighbor extrapolation possible (Fig. 1c). Moreover, the asobtained orientations of the β phase are now delimiting grain boundaries that correspond to the expected ex-β grain highlighted in the BC map (Fig. 1a). Each ex-β grain is thus actually indexed as a unique β orientation from this method, i.e. each α" variant in a given ex-β grain is indexed as a same and unique β orientation, as theoretically expected. It will also be shown further that the obtained orientation of the β phase is accurate in spite of the very low indexation rate. In turn, although the sample is only composed of α'' phase (Fig. 1d), directly indexing the EBSPs with the β phase operates as a crystallographic reconstruction of the parent phase.

It is now necessary to verify if the orientation of the β phase found by this "direct reconstruction" corresponds to its real orientation. In order to verify this orientation, the OR with the α" martensite is first tested and, finally, a thermal treatment is performed to allow a direct indexation of the real β phase by EBSD.

The pole figures of the β orientation and the five α'' orientations obtained from EBSPs indexation are given in Fig. 2. It can be clearly seen that the (001)α '' pole of each α" variant corresponds to a {110}β pole of the β parent phase and each [110]α '' direction corresponds to a

[111]β direction, which is in accordance with the OR between these two phases. The maximum misfit between poles of the two phases is 2° but is mainly measured to be lower than 1° (Fig. 2). Moreover, each different α" variant gives a same and unique orientation of the β phase that is consistent with the expected microstructure of a shape memory alloy. The orientation of the β phase determined by direct reconstruction is thus in accordance with the OR applied with the five α" orientations that are really present in the microstructure: the direct indexation of the martensitic phase by the parent β phase respects the OR.

Once the crystallographic reconstruction is achieved, one technical difficulty is to evaluate the error made on the orientation of the reconstructed crystal. This is generally done by orienting the residual high temperature phase (β phase in titanium alloys or austenite in steels) [4]. With the alloy family that is studied in this paper, no residual β phase is present in the microstructure, but it is possible to obtain the high temperature parent phase at room temperature after a short annealing above A f temperature followed by air cooling. The same sample investigated previously (Fig. 1) is thus annealed at 300°C for 1 minute and then air cooled to obtain the real β microstructure. A short annealing duration is required to prevent significant formation of ω phase which is favored in this temperature domain. As previously, the sample is then prepared once again for EBSD characterization by mechanical polishing and chemical etching. Fig. 3 shows an EBSD mapping of the same area as in Fig. 1 and the associated X-ray diffractogram after annealing. The grain morphology is quite different to the Fig. 1 due to the thin metallic layer that is removed during the new polishing step. It is clear that the lath morphology inherited from the martensitic transformation is no more present, and the indexation rate is now excellent (over 99%) when EBSPs are indexed with the β phase.

The diffractogram also confirms that the sample is mostly composed of β phase (Fig 3b). A few amount of ω phase is also detected. But, as ω phase is known to be nanometric sized and fully coherent with the β phase [16,19], its presence has not effect on the EBSPs indexation.

The parent phase microstructure obtained at room temperature after annealing is thus suitable to be compared with the orientation that was obtained by direct crystallographic reconstruction from the martensitic microstructure. An accurate measurement of the misorientation between the ex-β grain in the center of Fig. 1 and the same β grain on Fig. 3 is 4.4°. This operation is done on each grain of the map and gives an average value of 5.2°.

These relatively high values of misorientation are mainly due to the experimental conditions (annealing, sample repreparation and sample repositioning in the SEM between the two mappings). In order to prove that, the classical phase reconstruction is used in order to reconstruct the ex-β grain at the center of the figure, i.e. the OR is applied on each variant of martensite which is detected in this grain. The difference between this classical method and the one which is proposed in this paper is of 1.1° in average, proving that the rest of the error is due to the experimental conditions. It confirms that the β orientation that is determined on Fig. 2c from the martensitic microstructure corresponds to the actual orientation of the parent phase. The "direct reconstruction" from the martensitic microstructure gives thus the real orientation of the former high temperature phase with a rather good precision in spite of a very low indexation rate.

Obtaining the orientation of the parent phase by indexing EBSPs of martensite as it was the parent phase can seem very peculiar. This orientation matching can be explained by the similarity of the EBSPs of the two phases. As an example, Fig. 4 shows a simulated EBSP based on the Euler angles obtained from the β indexation of Fig. 1c and a simulated EBSP obtained with one of the five α'' orientations determined on the same domain (Fig. 1b). One can observe that the main Kikuchi band intersections on Fig. 4a and Fig. 4b match very closely. Moreover, the (101)β pole matches exactly to the (001)α '' pole, confirming that the OR is directly applicable for the indexation of the EBSPs of α" phase by the β phase. This similarity of EBSPs is due to the similarity of both structures wherein positions of atoms are very near: a small displacement of some atoms of the β phase leads to the orthorhombic structure of the martensite.

As a summary, we demonstrate that the orientation of the high temperature parent β phase in titanium shape memory alloys can be easily obtained from EBSD experiments by directly indexing the α'' martensitic microstructure with the parent phase. The quality of indexation decreases but is sufficient to obtain the real orientation of the parent phase, even if this phase is not present in the microstructure. This new method does not need any further data processing after EBSD experiments and is easier than those commonly used. But this method works when martensite can only transform into a same and unique parent crystal with simple martensitic transformations such as body-centered cubic to orthorhombic. It would not work with more complex transformations such as bcc to hcp or fcc to bcc. As it does not need indexation of the martensite, it can also be applied when the martensitic microstructure is too thin to be indexed unambiguously by EBSD to allow classical reconstruction.

Fig. 1 .

 1 Fig. 1. EBSD maps of Ti-25Ta-20Nb alloy: band contrast (BC) only (a); BC superimposed to Euler angles obtained from indexation with the α'' phase (b); BC superimposed to Euler angles obtained from indexation with the β phase (c); X-rays diffractogram (d).

Fig. 2 .

 2 Fig. 2. Poles figure of the ex-β grain in the center of the Fig. 1. The {110}β and <111>β obtained from EBSD indexation are superimposed to the (001)α '' and [110]α '' of the five orientations of α" detected by EBSD.

Fig. 3 .

 3 Fig. 3. EBSD mapping of the same area as in Fig. 1 after annealing: band contrast and Euler angles from β indexation (a) and X-rays diffractogram (b).


Antoine Perasso
email: antoine.perasso@univ-fcomte.fr

Ulrich Razafison
email: ulrich.razafison@univ-fcomte.fr

Identifiability problem for recovering the mortality rate in an age-structured population dynamics model

Keywords: Mathematics Subject Classification: 35Q92, 35R30, 92D25, 93B30 Parameter identifiability, population dynamics, age-structured model, transport PDE, non-local boundary condition

published or not. The documents may come

Identifiability problem for recovering the mortality rate in an age-structured population dynamics model

Antoine Perasso, Ulrich Razafison

Introduction

When studying a mathematical problem that is modeling realistic and observed phenomena, an important issue is to check if the model satisfies a parameter identifiability property, i.e. whether the unknown parameters involving in the model can be uniquely determined from the observation and sometimes from the initial condition and an input of the model. This property, called parameter identifiability problem, is an inverse problem that consists to check if the parameter to output map of the model is into, initial condition and input being known. The main goal of this article is to study a parameter identifiability problem in the context of a classical age-structured population dynamics model, the McKendrick-Von Foerster model. This latter describes the evolution in time of a population structured according to the age of the individuals, leading to a problem formulated by a linear partial differential equation of transport type with a non local boundary condition. This model is classically used to describe population dynamics [START_REF] Murray | Mathematical Biology I, An introduction[END_REF] or also in epidemiology, where the age variable represents the age of infection [START_REF] Diekmann | Mathematical Epidemiology of Infectious Diseases[END_REF]. In this article, we wonder whether the observation of an age group of the population may recover the uniqueness of a fundamental population dynamics parameter, the mortality rate.

To our knowledge, the identifiability analysis of the mortality rate in such a model has never been considered before. The main difficulties here are inherent of the age-dependence of the mortality rate and of the boundary condition stated by a non local loopback on the state of the system.

There is a well-established theory for the parameter identifiability of controlled and uncontrolled dynamical systems described by ordinary differential equations [START_REF] Chappell | Global identifiability of the parameters of nonlinear-systems with specified inputs: a comparison of methods[END_REF][START_REF] Walter | Identifiability of State Space Models[END_REF][START_REF] Walter | of Parametric Models[END_REF]. Three main approaches have been used: (i) the state isomorphism method [START_REF] Denis-Vidal | Equivalence and identification analysis of uncontrolled nonlinear dynamical systems[END_REF][START_REF] Vajda | Similarity transformation approach to identifiability analysis of nonlinear compartmental-models[END_REF]; (ii) the Taylor series expansion method [START_REF] Pohjanpalo | System identifiability based on the power series expansion of the solution[END_REF]; (iii) the algebro-differential elimination method [START_REF] Ljung | On global identifiability for arbitrary model parametrizations[END_REF][START_REF] Margaria | Structural identifiability analysis of some highly structured families of statespace models using differential algebra[END_REF][START_REF] Saccomani | Parameter identifiability of nonlinear systems: the role of initial conditions[END_REF][START_REF] Saccomani | An effective automatic procedure for testing parameter identifiability of hiv/aids models[END_REF], aiming at obtaining and exploiting algebro-differential relations between the input and output of the system. In infinite dimension, identifiability results exist for fairly general classes of linear problems. Results concerning convolutive systems, which include the delay-differential equations, can be found in [START_REF] Belkoura | Identifiability of systems described by convolution equations[END_REF][START_REF] Orlov | On identifiability of linear time-delay systems[END_REF][START_REF] Zhang | Parameter identifiability of nonlinear systems with timedelay[END_REF]. Identifiability results derived from the use of spectral theory are given in [START_REF] Nakagiri | Review of Japanese work of the last ten years on identifiability in distributed parameter systems[END_REF] for the 1-D heat and wave equations with boundary observations as well as for abstract homogeneous evolution equations with whole state observation. Results on various classes of linear models with pointwise observation where obtained using Carleman estimates, for instance for a non-stationary particle transport equation (see [START_REF] Klibanov | Global uniqueness for a coefficient inverse problem for the non-stationary transport equation via carleman estimate[END_REF] and references therein). In the nonlinear case there are only few results dealing with transport equations using algebrodifferential approach in inifite dimension [START_REF] Perasso | Identifiability analysis of an epidemiological model in a population[END_REF] or parabolic equations using either Carleman estimates [START_REF] Egger | Global uniqueness and hölder stability for recovering a nonlinear source term in a parabolic equation[END_REF][START_REF] Kaltenbacher | An inverse problem for a nonlinear parabolic equation with applications in population dynamics and magnetics[END_REF][START_REF] Klibanov | Global uniqueness of a multidimensional inverse nonlinear parabolic equation by a carleman estimate[END_REF] or when the observation is the state of the system [START_REF] Handrock-Meyer | Identifiability of distributed parameters for a class of quasilinear differential equations[END_REF].

This article is structured as follows. Section 2 is dedicated to the statements of the parameter identifiability problem, where an abstract formulation is performed in the linear case. Then we illustrate throughout the Malthus model how the nature of the observation may affect the identifiability problem. Section 3 contains the mortality rate identifiability analysis of the Von Foerster-McKendrick model. We introduce in section 3.1 this age-structured population model and formulate the parameter identifiability problem in Section 3.2. The identifiability analysis is split into two cases, the homogeneous case, with a null boundary condition, and the inhomogeneous case, with a non local boundary condition. Finally, we conclude with perspectives in Section 4.

Statement of the identifiability problem

Theoretical framework

Throughout this paper, we use the following notations,

R + = [0, +∞); X , Y are Banach spaces; I d Z is the identity element of the set Z ;
L (X , Y ), resp. L (X ), the space of linear operators with domain X and range in Y , resp. X ; Consider the following first order linear dynamical system on X ,

     d u(t ) d t = Au(t ), u(0) = u 0 ∈ X , (1)
where A : D(A) ⊂ X → X is a differential operator. Given a set P , an unknown parameter p ∈ P involving in system (1) and an observation linked to the state u(t ) of the system, we consider the model system S p , reflecting the p-dependence in system (1), defined by

S p :            d u(t ; p) d t = A p u(t ; p), u(0) = u 0 ∈ X , y p (t ) = Bu(t ; p), (2)
where B ∈ L (X , Y ) is called the observation operator and function y p is the system output. In all that follows, we suppose the following assumptions on model system S p :

(i) The initial condition u 0 is known.

(ii) The operator B is known.

(iii) The parameter p is unknown.

(iv) The operator A p generates a strongly continuous semigroup of operators {T p (t )} t ≥0 of L (X ).

Remark 1.

From assumption (iv), it follows that system S p has a unique solution in X given by u(t ) = T p (t )u 0 , t ≥ 0, and consequently the observation in S p is given by y p (t ) = B T p (t )u 0 for every t ≥ 0. If the context is clear, the p-dependence of the solution u(t ; p) of system S p will be omitted.

From the given system (2), consider the system S p obtained by replacing p by p in S p with p satisfying assumption (iii) and A p satisfying assumption (iv). For J ⊂ R + , the difference

ǫ(t ; (p, p)) = y p (t ) -y p (t ), t ∈ J , (3)
is called the output error linked to S p and S p on J . We now define the parameter identifiability problem.

Definition 2.1. The system S p is said P -identifiable on J if for every p, p ∈ P the condition p = p follows from the zero output error,

ǫ(t ; (p, p)) = 0 in Y , t ∈ J . ( 4
)
When the equality (4) fails, we say that the parameter identifiability problem is ill-posed on (P, J ), or simply ill-posed on J if the set P is certain.

The parameter identifiability problem as stated above consists in studying if the parameter to output map of the problem is into. This study induces two equalities to work with, in one hand an equality between parameters and in other hand between the system outputs. These latter equalities require to determine the suitable functional spaces in which they are satisfied. In particular, one can check that equality (4) is satisfied in the Banach space Y which is linked to the space X through the relation y = Bu. As a consequence, performing a mathematical analysis of the well posedness of the abstract Cauchy problem is a prerequisite in the study of the parameter identifiability problem.

A short example: mortality rate identifiability in the Malthus model

In this section is studied an identifiability problem related to the mortality rate in the classical population dynamics model of Malthus. We recall that this model is given by the following ODE system,

   d u d t (t ) = βu(t ) -µu(t ), u(0) = u 0 ∈ R + , (5)
where the quantity u(t ) denotes the population at time t ≥ 0, µ ≥ 0 is the mortality rate and β ≥ 0 the birth rate. We refer to [START_REF] Murray | Mathematical Biology I, An introduction[END_REF] for the historic and a description of this model. The identifiability problem we are interested in consists in the determination of the uniqueness of parameters β and µ together from the observation either of the whole population or of the death flow.

To achieve that goal, one can check that the notations of Section 2 rewrite the Banach space X = R, the unknown parameter p = β µ ∈ P = (R * + ) 2 and the differential operator with domain R as A p = (βµ)I d R . Finally, the model system S p is linked to (5) with a given output y p and an observation operator B ∈ L (R) that are described below.

It is clear that the identifiability of the model system S p depends on the nature of the observation. A biological intuition should make us believe that the observation of the whole population u(t ) at any time gives a lot of information on the parametrization of the model. The following proposition shows the opposite. Theorem 2.2. Suppose that the operator B : R → R is given by 1. B = µ.I d R , then the model system S p is P -identifiable on every J ⊂ R + .

B = I d R , then the parameter identifiability problem is ill-posed on R + ;

Proof. To prove the point 2 of the theorem, we have to exhibit (at least) two different parameters p, p ∈ P such that ǫ(t ; (p, p)) is zero for every t ≥ 0. The integration of System [START_REF] Egger | Global uniqueness and hölder stability for recovering a nonlinear source term in a parabolic equation[END_REF] implies that the operator (A p , R) is generator of the semigroup {T p (t )} t ≥0 on R defined by

T p (t ) = e (β-µ)t , t ≥ 0. ( 6
)
Then for fixed (β, µ) ∈ P , the family of parameters {p λ } given by

p λ = λβ (λ -1)β + µ , satisfies for every λ, λ 2 > max 0, β-µ β , p λ ∈ P, ǫ(t ; (p λ , p λ 2 )) = 0, t ≥ 0.
As a consequence, the parameter identifiability problem is ill-posed on R + . We now prove the point 1 of Theorem 2.2. A direct consequence of the assumption (i) about the initial data combined with the zero output error valued at t = 0 implies that µ = μ. Then the semigroup expression given in ( 6) implies the identifiability on every observation period J ⊂ R + .

For the observation operator B = I d R , the output is the density u(t ) meaning that the observation consists of the total population whereas for B = µI d R , the observation is the death flow. It is clear that one could either observe the birth flow, given by B = βI d R , and the result of Theorem 2.2.1 would still hold. Theorem 2.2 states that the knowledge of the whole population at any time is not a sufficient contribution to recover the uniqueness of parameters µ and β. On the contrary, the observation of the death flow, respectively of the birth flow, is sufficient. To conclude, Theorem 2.2 shows that the observation of a part of the population can be more useful to determine some parameters than the observation of the entirety.

Remark 2.

One can check that the knowledge of u 0 , stated by assumption (i) in Section 2, is necessary to conclude the identifiability result. In other word, the model system is identifiable under the condition that the state of the system has to be known for at least one time, what will be supposed in all that follows.

Mortality rate identifiability in an age-structured population dynamics model

We start with several notations of functional spaces that will be used in all that follows.

L 1 (R + ) is the set of Lebesgue integrable functions on R + ; L 1 l oc (R + ) is the set of locally Lebesgue integrable functions on R + ; L ∞ (R + )
is the set of Lebesgue measurable functions on R + which are essentially bounded;

W 1,1 (R + ) = {ϕ ∈ L 1 (R + ), ϕ ′ ∈ L 1 (R + )};
and for any functional space Z of real valued functions, the subset Z + ⊂ Z will denote the cone of non negative functions of Z .

An age structured population dynamics model

We consider a population density u(t , a) at time t ≥ 0 and age a ∈ [0, a † ) where a † ∈ (0, +∞] is the upper bound of lifespan. Such a population is described by the following classical problem called Von Foerster-Mck Kendrick problem [START_REF] Murray | Mathematical Biology I, An introduction[END_REF][START_REF] Perthame | Transport equations in biology[END_REF],

               ∂u(t , a) ∂t + ∂u(t , a) ∂a = -µ(a)u(t , a), u(t , 0) = a † 0 b(a)u(t , a) d a, u(0) = u 0 , (7)
with the assumptions

• µ ∈ L 1 l oc,+ (0, a † ) and a † 0 µ(a)d a = +∞; • b ∈ L ∞ + (0, a † ); • u 0 ∈ W 1,1 + (0, a † ).
In this model, the total number of individuals at time t with age a between a 1 and a 2 is

a 2 a 1 u(t , a) d a.
Function µ represents an age-dependent mortality rate of the population. Assumptions made on µ mean that mortality may occur at any age and all the individuals leave the population before a finite age. As a consequence, the following quantity

π(a) = exp - a 0 µ(s)d s , (8)
which satisfies π(0) = 1 and π(a † ) = 0, represents the probability for an individual to survive at age a ≥ 0. Since µ ∈ L 1 l oc (0, a † ), the probablity π is clearly a continuous function on [0, a † ). The non local boundary condition in [START_REF] Iannelli | Mathematical Theory of Age-Structured Population Dynamics[END_REF], called renewal of the population, induces the birth rate b which gives the repartition of birth for individuals with age a ∈ (0, a † ). When b is zero in L ∞ (0, a † ), there is no renewal of the population. For more details concerning the biological significance of the assumptions of the model, we refer to [START_REF] Iannelli | Mathematical Theory of Age-Structured Population Dynamics[END_REF][START_REF] Webb | Dynamics of populations structured by internal variables[END_REF].

We now make clear the functional framework to study the abstract Cauchy problem linked to [START_REF] Iannelli | Mathematical Theory of Age-Structured Population Dynamics[END_REF]. To this goal, let us denote in all that follows

X = L 1 (0, a † ) and for b ∈ L ∞ + (0, a † ) let A b : D(A b ) ⊂ X → X be the differential operator defined by D(A b ) = ϕ ∈ W 1,1 (0, a † ), ϕ(0) = a † 0 b(a)ϕ(a)d a , A b : ϕ → -ϕ ′ -µϕ.
The abstract Cauchy Problem related to [START_REF] Iannelli | Mathematical Theory of Age-Structured Population Dynamics[END_REF] is then

     d u(t ) d t = A b u(t ), u(0) = u 0 ∈ W 1,1 (0, a † ). (9)
The following theorem is a classical result, see for instance [START_REF] Thieme | Mathematics in Population Biology[END_REF][START_REF] Webb | Dynamics of populations structured by internal variables[END_REF],

Theorem 3.1. (A b , D(A b )
) is an infinitesimal generator of a strongly continuous semigroup {T b (t )} t ≥0 on X . Theorem 3.1 has several consequences in the study of the Cauchy Problem [START_REF] Klibanov | Global uniqueness of a multidimensional inverse nonlinear parabolic equation by a carleman estimate[END_REF], such as existence and uniqueness of a non negative mild solution u ∈ C (R + , X ) of the problem, such that (t , a) → u(t , a) is continuous outside the maniflod a -t = 0. Moreover, this solution is defined in a strong sense in case where u 0 ∈ D(A b ).

The mortality rate identifiability problem

In this section is stated an identifiability problem related to system [START_REF] Klibanov | Global uniqueness of a multidimensional inverse nonlinear parabolic equation by a carleman estimate[END_REF]. For convenience, we suppose that a † = +∞. It is easy (but tedious) to prove that all the identifiability results of Subsection 3.3 and 3.4 can be extended for a finite value of a † . The parameter identifiability problem is the following: one wonders whether the observation of an age group of the population helps to recover the uniqueness of the mortality rate µ on the age domain R + , or at least on some subsets of the age domain. To study this parameter identifiability problem, we make explicit the notations of Section 2 in the context of the age-structured model [START_REF] Iannelli | Mathematical Theory of Age-Structured Population Dynamics[END_REF], then we perform an abstract formulation of the identifiability problem. In all that follows we denote, according to notations of Section 2, the set of unknown parameters

P = L 1 l oc,+ (R + ) the observation operator B ∈ L (X ) defined for 0 ≤ a ⋆ < a ⋆ ≤ +∞ by B : X → X , ϕ → ϕχ (a ⋆ ,a ⋆ )
and the observation period

J = [0, T ], T > 0.
As stated by the notations above, the aim of the identifiability problem is to recover the uniqueness of the mortality rate µ ∈ P from the observation of the population density of the age group (a ⋆ , a ⋆ ) ⊂ R + . Given µ, μ ∈ P the dependence on parameter µ of any variable Φ is denoted Φ µ , and the output error related to µ, μ is denoted ǫ(t ; (µ, μ)), accordingly to Section 2. Finally, S µ is the model system reflecting the dependence on the unknown parameter µ, as sated in equation (2).

Remark 3.

1. The observation operator is defined from X to X . This is a particular case where the functional space of the output is Y = X .

2. We can check that the divergence of the integral +∞ 0 µ = +∞ is not included in the definition of the set of parameters P . Indeed, this assumption has to be considered for its biological convenience more than as a technical need to the identifiability results.

Even if the identifiability problem is ill-posed on a time interval J , meaning the non uniqueness of parameter µ in X , it may be possible to recover the uniqueness of an unknown parameter on some subsets of the age domain R + . To take this case under consideration, we now extend the definition of the parameter identifiability. Let K ⊂ R + be a compact subset such that |K | > 0 where | • | denotes the Lebesgue measure. Then the application • L 1 (K ) defines a semi-norm on L 1 l oc (R + ). Consider now the quotient projection

θ K : L 1 l oc (R + ) → L 1 l oc (R + )/ ker • L 1 (K )
This projection has the following sense: for a function f ∈ L 1 l oc (R + ), the equivalence class is given by θ

K ( f ) = {g ∈ L 1 l oc (R + ), f (a) = g (a) f.a.e a ∈ K }.
From here, we have the following definition, Definition 3.2. Let K ⊂ R + a compact set. The model system S µ is P / K -identifiable on J if for every µ, μ ∈ P the equality θ K (µ) = θ K ( μ) follows from the zero output error

ǫ(t ; (µ, μ)) = 0, t ∈ J .
The latter definition, which allows to recover the uniqueness of the restriction of the unknown parameter µ on some compact subsets of R + , will be frequently used in the sequel.

Taking into account Problem [START_REF] Iannelli | Mathematical Theory of Age-Structured Population Dynamics[END_REF] we consider in all that follows the notations

A = {a ∈ R + , u 0 (a) > 0}, a = inf A , K ⋆ = [0, a ⋆ ] and K ⋆ = [0, a ⋆ ], A ⋆ = A ∩ K ⋆ and A ⋆ = A ∩ K ⋆ .
The study of the mortality rate identifiability problem is split into two cases, the homogeneous case, that describes a population without renewal, and the inhomogeneous case, that is a population subject to renewal.

Identifiability results in the homogeneous case

We suppose in this case that the population is not subject to renewal, meaning that the boundary condition in ( 7) is homogeneous, translated by b = 0 in L ∞ (R + ). In that case, a necessary and sufficient condition on the initial density is given to ensure the uniqueness of the mortality rate µ on K ⋆ , the largest reachable compact subset of the age domain R + .

Lemma 3.3. Let T > a ⋆ be arbitrary but fixed and K

= A ⋆ ∪ [a ⋆ , a ⋆ ]. If a = 0, then the model system S µ is P / K -identifiable on [0, T ].
Proof. Using the characteristics of the transport equation in [START_REF] Iannelli | Mathematical Theory of Age-Structured Population Dynamics[END_REF], we check that the expression of the semigroup {T 0 µ (t )} t ≥0 generated by (A 0 , D(A 0 )) is given for every ϕ ∈ L 1 (R + ) by

T 0 µ (t )ϕ(a) = π(a)(π(a -t )) -1 ϕ(a -t ) for t ≤ a, 0 for 0 ≤ a ≤ t . ( 10
)
A consequence is the following alternative expression for the output,

y µ (t ) = B T 0 µ (t )u 0
inducing that the output in model system S p satisfies for (t , a)

∈ [0, T ] × R + y µ (t , a + t ) = u 0 (a)π(a + t )(π(a)) -1 χ (a ⋆ ,a ⋆ ) (a + t ). ( 11
)
Consider now two parameters µ, μ ∈ P such that for t ∈ [0, T ] the following equality holds in

L 1 (R + ), ǫ(t ; (µ, μ)) = 0.
Using the expression [START_REF] Ljung | On global identifiability for arbitrary model parametrizations[END_REF] of the output we then have for every (t , a)

∈ [0, T ] × R + , u 0 (a) π(a + t )(π(a)) -1 -π(a + t )( π(a)) -1 χ (a ⋆ ,a ⋆ ) (a + t ) = 0, (12)
where π, resp. π, is a notation used for convenience instead of π µ , resp. π μ.

Performing the change of variables (t , a) → (a + t , a) = (ξ, a) in the latter equation yields on

{(ξ, a) ∈ (a ⋆ , a ⋆ ) × A , max(0, ξ -T ) ≤ a ≤ ξ}, π(a)( π(a)) -1 = π(ξ)( π(ξ)) -1 . ( 13
)
Since T > a ⋆ and π, π are continuous functions, equation [START_REF] Murray | Mathematical Biology I, An introduction[END_REF] implies that there exists a constant c(a) independent of a and ξ such that

π(a)( π(a)) -1 = c(a), ∀a ∈ A ⋆ , π(ξ)( π(ξ)) -1 = c(a), ∀ξ ∈ [a ⋆ , a ⋆ ].
If a = 0, then necessarily c(0) = 1 since π, π are probabilities, which implies π(a) = π(a), ∀a ∈ K , and consequently the model system S µ is P / K -identifiable on [0, T ].

Theorem 3.4. Let T > a ⋆ be arbitrary but fixed. Then the model system

S µ is P / K ⋆-identifiable on [0, T ] if and only if A ⋆ is a dense subset of K ⋆ . Proof. If A ⋆ is a dense subset of K ⋆ , then a = 0 and A ⋆ ∪ [a ⋆ , a ⋆ ] = K ⋆ ,
so the sufficient condition is a direct consequence of Proposition 3.3. We now aim at proving that the condition is necessary. To achieve that goal, we build two different parameters µ, μ ∈ P / K ⋆ such that the zero output error stands in Y . If A ⋆ is not dense in K ⋆ , there exists a non empty interval U = (a -, a + ) ⊂ K ⋆ \A ⋆ . Let us consider µ and μ ∈ P such that

(i) µ -μ L 1 (U ) > 0, (ii) π(a) = π(a), ∀a ∈ R + \U .
where π, resp. π, is the probability to survive defined in [START_REF] Kaltenbacher | An inverse problem for a nonlinear parabolic equation with applications in population dynamics and magnetics[END_REF] and linked to µ, resp. to μ. One can check that condition (i) clearly implies that µ and μ are two different parameters of P / K ⋆. In order to reconstruct the outputs y µ (t ) and y μ(t ) over time, the difference yielded by (i) is offset in (ii) by imposing that the probabilities to survive π and π are the same outside U . One can check that this also implies that µ(a) = μ(a) f.a.e. a ∈ R + \U .

For a fixed t ≥ 0, we now reconstruct the outputs y µ (t ) and y μ(t ) on some subsets of the age domain R + , in order to prove that the output error is 0 in X . Due to [START_REF] Klibanov | Global uniqueness for a coefficient inverse problem for the non-stationary transport equation via carleman estimate[END_REF], it is obvious to check that ǫ(t ; (µ, μ))χ (0,t ) = 0 in X .

We now prove that the zero value of the outpu is satisfied for a ∈ (t , +∞). Suppose that a ≥ t , and consider the three following cases.

Case a ∈ [t + a + , +∞) Then a -t ≥ a + and a ≥ a + so (ii) implies that π(a -t ) = π(a -t ) and π(a) = π(a). Equation [START_REF] Klibanov | Global uniqueness for a coefficient inverse problem for the non-stationary transport equation via carleman estimate[END_REF] then yields u(t , a; µ) = u(t , a; μ). ( 14)

Case a ∈ [t + a -, t + a + ] Then a -t ∈ U ⊂ K ⋆ \A ⋆ so u 0 (a -t ) = 0 and (14) still holds.
Case a ∈ [a + , t + a -] This latter case holds when t ≥ a + -a -. Then a -t ≤ a -and assumption (ii) implies that π(at ) = π(a -t ). Moreover, since a ≥ a + we also have π(a) = π(a). So equality ( 14) is satisfied.

To conclude the several cases, equality ( 14) is satisfied for a ∈ [t , +∞)∩[min(a + , t +a -), +∞). Finally, since U ⊂ K ⋆ , then either the observed age period satisfies (a 1. The result of Theorem 3.4 can be easily extended to the case where a ⋆ = +∞, meaning that every individual of the population with age larger than a ⋆ is observed. In this particular case of observation, the model system is P -identifiable on [0, T ] if and only if A is dense in R + .

⋆ , a ⋆ )∩[t , +∞) ⊂ [min(a + , t + a -), +∞) or (a ⋆ , a ⋆ ) ∩ [t , +∞) =
2. In order to get the mortality rate identifiability, one could believe that it is necessary to observe the population on an horizon time T larger than a ⋆ the maximum age of observation. Theorem 3.4 shows that T > a ⋆ is sufficient.

Identifiability results in the inhomogeneous case

We now suppose that the population is subject to reproduction, implying an inhomogeneous boundary condition in (7) since b is not zero in L ∞ (R + ). This section shows that the boundary condition yields a fundamental perturbation in the study of the mortality rate identifiability problem. Indeed, whereas the sufficient condition on the initial density u 0 stated in Theorem 3.4 still holds to prove the identifiability result, a necessary condition on the initial density and on the boundary condition via the birth function b is stated in Theorem 3.5. Moreover, to prove that this necessary condition is not a sufficient one, we exhibit a subclass of parameters of that satisfy an identifiability property whereas u 0 and b are such that the necessary condition does not hold. This is done in Proposition 3.6.

Extension of Theorem 3.4 in the inhomogeneous case Theorem 3.5.

1. Let T > a ⋆ be arbitrary but fixed. If A ⋆ is dense in K ⋆ then the model system S µ is P / K ⋆- identifiable on [0, T ].

If there exists an open set U ⊂ K

⋆ \A ⋆ such that b(a) = 0 f.a.e a ∈ U , then the parameter identifiability problem is ill-posed on [0, T ] for any T ≥ 0.

Proof. Contrary to homogeneous boundary condition case, where equation ( 10) gives an explicit expression of the semigroup {T 0 (t )} t ≥0 , it is not possible to make explicit such an expression for {T b (t )} t ≥0 when b is not zero in L ∞ (R + ), which is due to the non local loopback boundary condition in [START_REF] Iannelli | Mathematical Theory of Age-Structured Population Dynamics[END_REF]. However, a change of variables using the characteristics of the transport equation proves that the state of the system u(t , a) satisfies the following implicit formulation

u(t , a) = π(a)(π(a -t )) -1 u 0 (a -t ) for t ≤ a, π(a)u(t -a, 0) for 0 ≤ a ≤ t , (15)
that will be used in the following. From this expression, one can check that the proof of the point 1. of Theorem 3.5 holds by the same arguments that ones developed to prove the sufficient condition in Theorem 3.4.

So we now focus on the proof of the point 2 of the theorem. Without lack of generality, one can suppose that U = (a -, a + ) an open interval and consider two parameters µ, μ ∈ P that satisfy the assumptions (i)-(ii) as in the proof of Theorem 3.4. To prove that the output error is zero, the main difficulty is to get the zero value for a ∈ (0, t ). Indeed, [START_REF] Orlov | On identifiability of linear time-delay systems[END_REF] implies that it is necessary to reconstruct the boundary conditions u(t , 0; µ) and u(t , 0; μ) to reach that goal. Equation ( 15) obviously implies, as in the proof of Theorem 3.4, the following zero value of the output error,

ǫ(t ; (µ, μ))χ (t ,+∞) = 0 in X . (16)
Moreover, equation [START_REF] Orlov | On identifiability of linear time-delay systems[END_REF] implies that the difference between boundary conditions in [START_REF] Iannelli | Mathematical Theory of Age-Structured Population Dynamics[END_REF] rewrites for every t ≥ 0,

u(t , 0; µ) -u(t , 0, μ) = t 0 b(a) u(t -a, 0; µ)π(a) -u(t -a, 0; μ π(a)) d a + +∞ t b(a) u(t , a, µ) -u(t , a; μ) d a.
Let us prove that the latter integral in the equation above is zero. To achieve that goal, one can check that equality ( 14) can be extended on [t , +∞) ∩ [0, a -]. Indeed, in case where a ∈ [0, a -] one gets a -t ≤ a ≤ a -and assumption (ii) implies π(a -t ) = π(a -t ) and π(a) = π(a), so we finally get

u(t , a; µ) = u(t , a; μ), ∀a ∈ [t , +∞) ∩ E t , (17)
where

E t = [0, a -] ∪ [min(a + , t + a -), +∞). Then (17) yields +∞ t b(a) u(t , a, µ) -u(t , a; μ) χ E t (a) d a = 0, ∀t ≥ 0.
Moreover, from the definition of E t one deduces that t ≥0

(R + \E t ) ⊂ U , and since b(a) = 0 f.a.e. a ∈ U , one gets

+∞ t b(a) u(t , a, µ) -u(t , a; μ) χ R + \E t (a) d a = 0, ∀t ≥ 0.
Consequently, the difference between boundary conditions satisfies

u(t , 0; µ) -u(t , 0, μ) = t 0 b(a) u(t -a, 0; µ)π(a) -u(t -a, 0; μ π(a)) d a.
From here, since b(a) = 0 f.a.e. a ∈ U and π(a) = π(a) for every a ∈ R + \U consequently to assumption (ii), a change of variables implies that the difference f (t ) = u(t , 0; µ) -u(t , 0, μ) is a solution of the following Volterra integral equation

f (t ) = t 0 b(t -a)π(t -a) f (a) d a, ∀t ≥ 0. Since b ∈ L ∞ (R + ) and π is a probability, such a function f satisfies | f (t )| ≤ b L ∞ t 0 f (a) d a, ∀t ≥ 0,
From here, since U ⊂ K ⋆ the infimum a ⋆ of the age range of observation satisfies a ⋆ ≥ a + and one gets from assumption (ii) (proof of Th.3.4), π(a) = π(a) for a ∈ [a ⋆ , a ⋆ ], so [START_REF] Saccomani | An effective automatic procedure for testing parameter identifiability of hiv/aids models[END_REF] combined with [START_REF] Orlov | On identifiability of linear time-delay systems[END_REF] imply

ǫ(t ; (µ, μ))χ (0,t ) = 0 in X . ( 20
)
Equations ( 16) and [START_REF] Saccomani | Parameter identifiability of nonlinear systems: the role of initial conditions[END_REF] together end the proof.

Illustration of the ill-posed case A consequence of Theorem 3.5.2 is that for a population that is composed of juveniles and adults that is described by the age-structured model [START_REF] Iannelli | Mathematical Theory of Age-Structured Population Dynamics[END_REF], it is necessary to start the observation when birth occur to get the identifiability of mortality rate µ.

Let us illustrate that case by considering a population with the following characteristics: Hypotheses of Theorem 3.5.2 are satisfied with K := K ⋆ = K ⋆ = [0, a ⋆ ], A ⋆ = and U = (0, a ⋆ ) ⊂ K \A ⋆ . This latter theorem then states that there exists different mortality rates of P that lead to the same observations, or in other words that it is possible to find µ, μ ∈ P such that θ K (µ) = θ K ( μ) but y µ = y μ. We therefore aim at building such different mortality rates leading to the same observation. Consider the mortality rate µ ∈ P given by

µ(x) = c x a ⋆ χ [0,a ⋆ ] + 1 -a ⋆ 1 -x χ [a ⋆ ,1) , c > 0
To illustrate the ill-posed case, we build a parameter μ ∈ P by performing some perturbations of µ on K = [0, a ⋆ ] such that θ K (µ) = θ K ( μ) while y µ (t , a) = y μ(t , a) . In other word, we seek to maximize μ ∈ P → µμ L 1 (K ) under the constraint {π(a ⋆ ) = π(a ⋆ )}. To reach that goal, let us choose mortality rates as an affine perturbation on K of µ given for any α ∈ (0, a ⋆ 2 ] by

μα (x) = cx 2α χ [0,α] + c 2 χ [α,a ⋆ -α] + c x -a ⋆ 2α + 1 χ [a ⋆ -α,a ⋆ ] + µ(x)χ [a ⋆ ,1) . ( 21
)
One can easily check that θ K ( μα ) = θ K (µ) when α = a ⋆ 2 and θ K ( μα ) = θ K (µ) else. Figure 1 shows the mortality rate µ and some perturbations μα , where the numerical values are given below.

The simulations are made using an upwind scheme based on finite volume method with the following numerical values:

• T = 3, a ⋆ = 0.4, ǫ = 0.4, c = 4, α i = a ⋆ 3 i , i ∈ N 0 0.5 1
i = u(t ; µ) -u(t ; μα i ) L 1 (0,1) for i ∈ {1, 2, 3}. • u 0 (a) = (1 -χ [0,a ⋆ ] ) + 100(a -a ⋆ )χ [a ⋆ ,a ⋆ +0.1] + 10χ [a ⋆ +0.1,1] • b(a) = (1 -χ [0,a ⋆ ] ) + 200(-a 2 + (1 + a ⋆ )x -a ⋆ )
As we could expect, these simulations illustrate well the ill-posed case. Indeed, whereas Figure 2 shows that the L 1 -norms between u and u i are different, Figure 3 shows that, at final time t = T , the ouptuts are superposed and consequently nondifferentiable at any age that is observed, even if the densities u and u i are different before the age of observation a ⋆ = 0.4. Finally, Figure 4 shows that the superposition of the curves for ages upper a ⋆ holds for any time t ≤ T .

The sufficient condition in Theorem 3.5 is not necessary Contrary to Theorem 3.4 that states in the homogeneous case a necessary and sufficient condition dealing with the initial density to get identifiability, we can see in Theorem 3.5 that the loopback boundary condition yields a fundamental perturbation in the study of the parameter identifiability problem. Indeed, the following result shows that when dealing with a non homogeneous boundary condition, the density of A ⋆ in K ⋆ is not a necessary condition to get identifiability on P .

Let us suppose that there exists an open interval (a -, a + ) ⊂ K ⋆ such that Proof. Let T > a ⋆ + a -and consider µ, μ ∈ P such that ǫ(t ; (µ, μ)) = 0 in X for t ∈ [0, T ]. To prove the identifiablity result, we aim at proving that the probabilities to survive satisfy

π(a) = π(a), ∀a ∈ [0, a ⋆ ].
We recall that the population density u(t , a) satisfies the equality [START_REF] Orlov | On identifiability of linear time-delay systems[END_REF]. Then for a ∈ (a ⋆ , a ⋆ ) and t ≤ a the zero value of the output error implies

u 0 (t ) π(a)(π(t )) -1 -π(a)( π(t )) -1 = 0. From here, assumption (a 2 ) gives π(a)(π(t )) -1 = π(a)( π(t )) -1 on {(t , a) ∈ [0, T ]\(a -, a + )×(a ⋆ , a ⋆ ), a ≥ t }, which is non empty since T > a -. A direct
If a ⋆ ≤ a -then π = π on K ⋆ . Suppose now that a ⋆ > a -and let us prove that π = π on [a -, a + ]. prove this result, one can firstly check that the boundary conditions linked to µ and μ are identically on the horizon time [0, T -min(a ⋆ , a + )]. Indeed the zero value of the output error and the characteristic expression [START_REF] Orlov | On identifiability of linear time-delay systems[END_REF] imply that for a ∈ (a ⋆ , a ⋆ ) and t ≥ a, π(a)u(t -a, 0; µ) = π(a)u(t -a, 0; μ), and ( 22 From the definition of the boundary condition in Problem [START_REF] Iannelli | Mathematical Theory of Age-Structured Population Dynamics[END_REF] and assumption (a 2 ) on b, the latter equation yields (

) 23
Since T -min(a ⋆ , a + ) ≥ a -one can consider [START_REF] Walter | Identifiability of State Space Models[END_REF] for t ∈ [0, a -]. Then using the expression on characteristics [START_REF] Orlov | On identifiability of linear time-delay systems[END_REF], equation [START_REF] Vajda | Similarity transformation approach to identifiability analysis of nonlinear compartmental-models[END_REF] and assumption (a 1 ) on u 0 one gets the equality

Consider now the set E = {a ∈ [a -, a + ], π(x) = π(x) ∀x ≤ a}. The set E is clearly non empty since a -∈ E so one can consider a 0 = max E . Suppose by contradiction that a 0 < a + . Then there exists 0 < δ ≤ a + -a 0 such that π(a) -π(a) = 0 for every a ∈ (a 0 , a 0 + δ). Consider a fixed t ∈ (a 0 -a -, δ+a 0 -a -). Then we have t +a -≤ a + . Moreover, from assumption (a 3 ) we have δ + a 0 -a -≤ a -and so t < a -. Consequently equation ( 24) rewrites for t ∈ (a 0a -, δ + a 0 -a -) as The latter equality implies a contradiction since in one hand b(a)u 0 (a-t ) π(a-t ) > 0 f.a.e. a ∈ (a 0 , t + a -) and in other hand π(a) -π(a) = 0 for every a ∈ (a 0 , t + a -) ⊂ (a 0 , a 0 + δ). So a + = max E and π = π on [a -, a + ]. Finally, the equation ( 22) then implies that µ(a) = μ(a) f.a.e a ∈ K ⋆ which proves that the model system S µ is P / K ⋆-identifiable on [0, T ].

Perspectives

The present work, dealing with the mortality rate identifiability in an age-structured model, shows how the non local loopback boundary condition affects the study of the uniqueness of the parameter. The inverse problem is studied from the expression of the semigroup generated by the differential operator. One think that the present results could be adpated to the identifiability of functional parameters, such as growth function, in the case of size-structured models. Indeed, the latter models are described by transport equations with non constant velocities in the transport term, that may derive a semigroup formulation as in age-structured models.

  so one gets ǫ(t ; (µ, μ))χ (t ,+∞) = 0 in X , which ends the proof.

  and a standard Gronwall argument implies| f (t )| ≤ f (0) e t b L ∞ , ∀t ≥ 0. (18)Moreover, the assumption (i) stated in Section 2 on the initial condition yieldsf (0) = +∞ 0 b(a) u 0 (a; µ) -u 0 (a; μ) d a = 0,and (18) finally gives u(t , 0; µ) -u(t , 0, μ) = 0, ∀t ≥ 0.

Figure 1 :

 1 Figure 1: Mortality rate µ and perturbations μα i for i ∈ {1, 2, 3}.

(a 1 ) 3 Figure 3 :Proposition 3 . 6 .

 13336 Figure 3: Age-structured densities u = u(T ; µ) and u i = u(T ; μα i ), i ∈ {1, 2, 3}, at final time t = T .

  consequence of the continuity of the probability functions π, π is thenπ(a) = π(a), ∀a ∈ [0, a -] ∪ [min(a ⋆ , a + ), a ⋆ ].

Figure 4 :

 4 Figure 4: Difference between u(µ) and u( μα 3 ) over age and time.

  )(u(t , a; µ) -u(t , a; μ)) d a = 0, ∀t ∈ [0, T -min(a ⋆ , a + )].

  min(t +a -,a + ) a - b(a)u 0 (a -t ) π(a -t ) (π(a) -π(a)) d a = 0, ∀t ∈ [0, a -].

  )u 0 (a -t ) π(a -t ) (π(a) -π(a)) d a = 0.


Pablo Rauzy

Sylvain Guilley

Formal Analysis of CRT-RSA Vigilant's Countermeasure Against the BellCoRe Attack A Pledge for Formal Methods in the Field of Implementation Security

Keywords: RSA (Rivest, Shamir, Adleman), CRT (Chinese Remainder Theorem), fault injection, Bell-CoRe (Bell Communications Research) attack, formal proof, OCaml

In our paper at PROOFS 2013, we formally studied a few known countermeasures to protect CRT-RSA against the BellCoRe fault injection attack. However, we left Vigilant's countermeasure and its alleged repaired version by Coron et al. as future work, because the arithmetical framework of our tool was not sufficiently powerful. In this paper we bridge this gap and then use the same methodology to formally study both versions of the countermeasure. We obtain surprising results, which we believe demonstrate the importance of formal analysis in the field of implementation security. Indeed, the original version of Vigilant's countermeasure is actually broken, but not as much as Coron et al. thought it was. As a consequence, the repaired version they proposed can be simplified. It can actually be simplified even further as two of the nine modular verifications happen to be unnecessary. Fortunately, we could formally prove the simplified repaired version to be resistant to the BellCoRe attack, which was considered a "challenging issue" by the authors of the countermeasure themselves.

Introduction

Private information protection is a highly demanded feature, especially in the context of global defiance against most infrastructures, assumed to be controlled by governmental agencies. Properly used cryptography is known to be a key building block for secure information exchange. However, in addition to the threat of cyber-attacks, implementation-level hacks are also to be considered seriously. This article deals specifically with the protection of a decryption or signature crypto-system (called RSA [START_REF] Rivest | A Method for Obtaining Digital Signatures and Public-Key Cryptosystems[END_REF]) in the presence of hardware attacks (e.g., we assume the attacker can alter the RSA computation while it is being executed).

It is known since 1997 [START_REF] Boneh | On the Importance of Checking Cryptographic Protocols for Faults[END_REF] that injecting faults during the computation of CRT-RSA could yield to malformed signatures that expose the prime factors (p and q) of the public modulus (N = p • q). Notwithstanding, computing without the fourfold acceleration conveyed by the CRT is definitely not an option in practical applications. Therefore, many countermeasures have appeared that consist in step-wise internal checks during the CRT computation. Last year we formally studied some of them [START_REF] Rauzy | A Formal Proof of Countermeasures Against Fault Injection Attacks on CRT-RSA[END_REF]. We were able to formally prove that the unprotected implementation of CRT-RSA as well as the countermeasure proposed by Shamir [START_REF] Shamir | Method and apparatus for protecting public key schemes from timing and fault attacks[END_REF] are broken, and that Aumüller et al.'s countermeasure [ABF + 02] is resistant to the Bell-CoRe attack. However, we were not able to study Vigilant's countermeasure [START_REF] Vigilant | Countermeasure securing exponentiation based cryptography[END_REF][START_REF] Vigilant | RSA with CRT: A New Cost-Effective Solution to Thwart Fault Attacks[END_REF] or its repaired version by Coron et al. [CGM + 10], because the tool we developed lacked the ability to work with arithmetical properties necessary to handle these countermeasures. In particular, our framework was unaware of the easiness to compute the discrete logarithm in some composite degree residuosity classes. These difficulties were foreseen by Coron et al. themselves in the conclusion of their paper [CGM + 10]:

"Formal proof of the FA-resistance of Vigilant's scheme including our countermeasures is still an open (and challenging) issue." This is precisely the purpose of this paper. The state-of-the-art of formal proofs of Vigilant's countermeasure is the work of Christofi, Chetali, Goubin, and Vigilant [START_REF] Christofi | Formal verification of an implementation of CRT-RSA Vigilant's algorithm[END_REF]. However, for tractability reasons, the proof is conducted on reduced versions of the algorithms parameters. One fault model is chosen authoritatively (viz. the zeroization of a complete intermediate data), which is a strong assumption. In addition, the verification is conducted on a specific implementation in pseudocode, hence concerns about its portability after compilation into machine-level code.

Contributions. We improved finja (our tool based on symbolic computation in the framework of modular arithmetic [START_REF] Rauzy | A Formal Proof of Countermeasures Against Fault Injection Attacks on CRT-RSA[END_REF]) to enable the formal study of CRT-RSA Vigilant's countermeasure against the BellCoRe attack. The finja tool allows a full fault coverage of CRT-RSA algorithm, thereby keeping the proof valid even if the code is transformed (e.g., optimized, compiled, partitioned in software/hardware, or equipped with dedicated countermeasures). We show that the original countermeasure [START_REF] Vigilant | RSA with CRT: A New Cost-Effective Solution to Thwart Fault Attacks[END_REF] is indeed broken, but not as much as Coron et al. thought it was when they proposed a manually repaired version of it [CGM + 10]. We simplify the repaired version accordingly, and formally prove it resistant to the BellCoRe attack under a fault model that considers permanent faults (in memory) and transient faults (one-time faults, even on copies of the secret key parts, e.g., during their transit on buses or when they reside on register banks), with or without forcing at zero, and with possibly faults at various locations. Thanks to the formal analysis, we are able to simplify the countermeasure even further: two of the nine checks are unnecessary.

Organization of the paper. We recall CRT-RSA cryptosystem and the BellCoRe attack in Sec. 2. Vigilant's countermeasure and its variant by Coron et al. are exposed in Sec. 3. In Sec. 4 we detail the modular arithmetic framework used by our tool. The results of our analysis are presented in Sec. 5. Finally, conclusions and perspectives are given in Sec. 6. After that, Appx. A details the practical issues related to fault injection analysis, for single-and multiple-fault attacks.

CRT-RSA and the BellCoRe Attack

This section recaps known results about fault attacks on CRT-RSA (see also [START_REF] Koç | High-Speed RSA Implementation[END_REF], [START_REF]Introduction to Hardware Security and Trust[END_REF]Chap. 3] and [START_REF] Joye | Fault Analysis in Cryptography[END_REF]Chap. 7 & 8]). Its purpose is to settle the notions and the associated notations that will be used in the later sections (that contain novel contributions).

RSA

RSA is both an encryption and a signature scheme. It relies on the identity that for all message 0 ≤ M < N , (M d ) e ≡ M mod N , where d ≡ e -1 mod ϕ(N ), by Euler's theorem 1 . In this equation, ϕ is Euler's totient function, equal to ϕ(N ) = (p -1) • (q -1) when N = p • q is a composite number, product of two primes p and q. For example, if Alice generates the signature S = M d mod N , then Bob can verify it by computing S e mod N , which must be equal to M unless Alice is only pretending to know d. Therefore (N, d) is called the private key, and (N, e) the public key. In this paper, we are not concerned about the key generation step of RSA, and simply assume that d is an unknown number in 1, ϕ(N ) = (p -1) • (q -1) . Actually, d can also be chosen equal to the smallest value e -1 mod λ(N ), where λ(N ) = (p-1)•(q-1) gcd(p-1,q-1) is the Carmichael function (see PKCS #1 v2.1, §3.1).

CRT-RSA

The computation of M d mod N can be speeded-up by a factor four by using the Chinese Remainder Theorem (CRT). Indeed, figures modulo p and q are twice as short as those modulo N . For example, for 2, 048 bits RSA, p and q are 1, 024 bits long. The CRT-RSA consists in computing S p = M d mod p and S q = M d mod q, which can be recombined into S with a limited overhead. Due to the little Fermat theorem (special case of the Euler theorem when the modulus is a prime), S p = (M mod p) d mod (p-1) mod p. This means that in the computation of S p , the processed data have 1, 024 bits, and the exponent itself has 1, 024 bits (instead of 2, 048 bits). Thus the multiplication is four times faster and the exponentiation eight times faster. However, as there are two such exponentiations (modulo p and q), the overall CRT-RSA is roughly speaking four times faster than RSA computed modulo N .

This acceleration justifies that CRT-RSA is always used if the factorization of N as p • q is known. In CRT-RSA, the private key is a more rich structure than simply (N, d): it is actually comprised of the 5-tuple (p, q, d p , d q , i q ), where:

• d p . = d mod (p -1),
• d q .

= d mod (q -1), and

• i q . = q -1 mod p.
The CRT-RSA algorithm is presented in Alg. 1. It is straightforward to check that the signature computed at line 3 belongs to 0, p • q -1 . Consequently, no reduction modulo N is necessary before returning S.

Algorithm 1: Unprotected CRT-RSA Input : Message M , key (p, q, d p , d q , i q ) Output: Signature M d mod N

1 S p = M dp mod p // Signature modulo p 2 S q = M dq mod q // Signature modulo q 3 S = S q + q • (i q • (S p -S q ) mod p) // Recombination 4 return S

The BellCoRe Attack

In 1997, an dreadful remark has been made by Boneh, DeMillo and Lipton [START_REF] Boneh | On the Importance of Checking Cryptographic Protocols for Faults[END_REF], three staff of BellCoRe: Alg. 1 could reveal the secret primes p and q if the line 1 or 2 of the computation is faulted, even in a very random way. The attack can be expressed as the following proposition.

Proposition 1 (Original BellCoRe attack). If the intermediate variable S p (resp. S q ) is returned faulted as S p (resp. S q )2 , then the attacker gets an erroneous signature S, and is able to recover p (resp. q) as gcd(N, S -S).

Proof. For any integer x, gcd(N, x) can only take 4 values:

• 1, if N and x are coprime, • p, if x is a multiple of p, • q, if x is a multiple of q,
• N , if x is a multiple of both p and q, i.e., of N .

In Alg. 1, if S p is faulted (i.e., replaced by S p = S p ), then S -S = q • ((i q • (S p -S q ) mod p) -(i q • ( S p -S q ) mod p)), and thus gcd(N, S -S) = q. If S q is faulted (i.e., replaced by S q = S q ), then S -S ≡ (S q -S q ) -(q mod p) • i q • (S q -S q ) ≡ 0 mod p because (q mod p) • i q ≡ 1 mod p, and thus S -S is a multiple of p. Additionally, S -S is not a multiple of q. So, gcd(N, S -S) = p.

Vigilant Countermeasure Against the BellCoRe Attack

Fault attacks on RSA can be thwarted simply by refraining from implementing the CRT. If this is not affordable, then the signature can be verified before being outputted. Such protection is efficient in practice, but is criticized for two reasons. First of all, it requires an access to e; second, the performances are incurred by the extra exponentiation needed for the verification. This explains why several other countermeasures have been proposed.

Until recently, none of these countermeasures had been proved. In 2013, Christofi, Chetali, Goubin, and Vigilant [START_REF] Christofi | Formal verification of an implementation of CRT-RSA Vigilant's algorithm[END_REF] formally proved an implementation of Vigilant's countermeasure. However, for tractability purposes, the proof is conducted on reduced versions of the algorithms parameters. One fault model is chosen authoritatively (the zeroization of a complete intermediate data), which is a strong assumption. In addition, the verification is conducted on a specific implementation in pseudocode, hence concerns about its portability after compilation into machine-level code. The same year, we formally studied [START_REF] Rauzy | A Formal Proof of Countermeasures Against Fault Injection Attacks on CRT-RSA[END_REF] the ones of Shamir [START_REF] Shamir | Method and apparatus for protecting public key schemes from timing and fault attacks[END_REF] and Aumüller et al. [ABF + 02] using a tool based on the framework of modular arithmetic, thus offering a full fault coverage on CRT-RSA algorithm, thereby keeping the proof valid even if the code is transformed (e.g., optimized, compiled or partitioned in software/hardware). We proved the former to be broken and the latter to be resistant to the BellCoRe attack.

Vigilant's Original Countermeasure

In his paper at CHES 2008, Vigilant [START_REF] Vigilant | RSA with CRT: A New Cost-Effective Solution to Thwart Fault Attacks[END_REF] proposed a new method to protect CRT-RSA. Its principle is to compute M d mod N in Z N r 2 where r is a random integer that is coprime with N . Then M is transformed into M * such that

M * ≡ M mod N, 1 + r mod r 2 ,
which implies that

S * = M * d mod N r 2 ≡ M d mod N, 1 + dr mod r 2 .
The latter results are based on the binomial theorem, which states that:

(1 + r) d = d k=0 d k r k = 1 + dr + d 2 r 2 + higher powers of r, (1)
which simplifies to 1 + dr in the ring Z r 2 . This property can be seen as a special case of an easy computation of a discrete logarithm in a composite degree residuosity class (refer for instance to [START_REF] Paillier | Public-Key Cryptosystems Based on Composite Degree Residuosity Classes[END_REF]). Therefore the result S * can be checked for consistency modulo r 2 . If the verification S * ? = 1 + dr mod r 2 succeeds, then the final result S = S * mod N is returned. Concerning the random numbers used in the countermeasure, its author recommend using a 32 bits integer for r, and 64 bits for the R i s. The original algorithm for Vigilant's countermeasure is presented in Alg. 49.

Coron et al. Repaired Version

At FDTC 2010, Coron et al. [CGM + 10] proposed some corrections to Vigilant's original countermeasure. They claimed to have found three weaknesses in integrity verifications:

1. one for the computation modulo p -1;

2. one for the computation modulo q -1;

3. and one in the final check modulo N r 2 .

They propose a fix for each of these weaknesses. To correct weakness #1, lines 11 to 15 in Alg. 49 becomes:

Algorithm 3: Vigilant's CRT-RSA: Coron et al. fix #1

1 p minusone = p -1 2 d ′ p = d p + R 1 • p minusone 3 S pr = M ′ p d ′ p mod p ′ 4 p minusone = p -1 5 if d ′ p ≡ d p mod p minusone then 6 return error 7 end
Similarly, for weakness #2, lines 32 to 36 in Alg. 49:

Algorithm 4: Vigilant's CRT-RSA: Coron et al. fix #2 1 q minusone = q -1 2 d ′ q = d q + R 2 • q minusone 3 S qr = M ′ q d ′ q mod q ′ 4 q minusone = q -1 5 if d ′ q ≡ d q mod q minusone then 6 return error 7 end
These two fixes immediately looked suspicious to us. Indeed, Coron et al. suppose that the computations of p -1 and q -1 are factored as an optimization, while they are not in the original version of Vigilant's countermeasure. This kind of optimization, called "common subexpression elimination", is very classical but is, as for most speed-oriented optimization, not a good idea in security code. We will see in Sec. 5 that the original algorithm, which did not include this optimization, was not at fault here.

N = p • q 2 if p • q • (S -R 4 -q • iq • (R 3 -R 4 )) ≡ 0 mod N r 2 then 3 return error 4 end
What has changed is the recomputation of p • q instead of reusing N in the verification. The idea is that if p or q is faulted when N is computed then error will be returned. In the original countermeasure, if p or q were to be faulted, the faulted version of N would appear in the condition twice, simplifying itself out and thus hiding the fault.

Formal Methods

Our goal is to prove that the proposed countermeasures work, i.e., that they deliver a result that does not leak information about neither p nor q (if the implementation is subject to fault injection) exploitable in a BellCoRe attack. In addition, we wish to reach this goal with the two following assumptions:

• our proof applies to a very general attacker model, and

• our proof applies to any implementation that is a (strict) refinement of the abstract algorithm. Definition 1 (CRT-RSA). The CRT-RSA computation takes as input a message M , assumed known by the attacker but which we consider to be random, and a secret key (p, q, d p , d q , i q ). Then, the implementation is free to instantiate any variable, but must return a result equal to: S = S q + q • (i q • (S p -S q ) mod p), where:

CRT-RSA and Fault Injection

• S p = M dp mod p, and

• S q = M dq mod q.
Definition 2 (fault injection). An attacker is able to request RSA computations, as per Def. 1. During the computation, the attacker can modify any intermediate value by setting it to either a random value or zero. At the end of the computation the attacker can read the result.

Of course, the attacker cannot read the intermediate values used during the computation, since the secret key and potentially the modulus factors are used. Such "whitebox" attack would be too powerful; nonetheless, it is very hard in practice for an attacker to be able to access reliably to intermediate variables, due to protections and noise in the side-channel leakage (e.g., power consumption, electromagnetic emanation). Remark that our model only takes into account fault injection on data; the control flow is supposed not to be modifiable.

Remark 1. We notice that the fault injection model of Def. 2 corresponds to that of Vigilant [START_REF] Vigilant | RSA with CRT: A New Cost-Effective Solution to Thwart Fault Attacks[END_REF], with the exception that the conditional tests can also be faulted. To summarize, an attacker can modify a value in the global memory (permanent fault ), and modify a value in a local register or bus (transient fault ), but cannot inject a permanent fault in the input data (message and secret key), nor modify the control flow graph.

The independence of the proofs on the algorithm implementation demands that the algorithm is described at a high level. The two properties that characterize the relevant level are as follows:

1. The description should be low level enough for the attack to work if protections are not implemented.

Any additional intermediate variable that would appear during refinement could be the target of an

attack, but such a fault would propagate to an intermediate variable of the high level description, thereby having the same effect.

From those requirements, we deduce that:

1. The RSA description must exhibit the computation modulo p and q and the CRT recombination; typically, a completely blackbox description, where the computations would be realized in one go without intermediate variables, is not conceivable.

2. However, it can remain abstract, especially for the computational parts. For instance a fault in the implementation of the multiplication (or the exponentiation) is either inoffensive, and we do not need to care about it, or it affects the result of the multiplication (or the exponentiation), and our model takes it into account without going into the details of how the multiplication (or exponentiation) is computed.

In our approach, the protections must thus be considered as an augmentation of the unprotected code, i.e., a derived version of the code where additional variables are used. The possibility of an attack on the unprotected code attests that the algorithm is described at the adequate level, while the impossibility of an attack (to be proven) on the protected code shows that added protections are useful in terms of resistance to attacks.

Remark 2. The algorithm only exhibits evidence of safety. If after a fault injection, the algorithm does not simplify to an error detection, then it might only reveal that some simplification is missing. However, if it does not claim safety, it produces a simplified occurrence of a possible weakness to be investigated further.

How finja Works

Our tool3 , finja, works within the framework of modular arithmetic, which is the mathematical framework of CRT-RSA computations. The general idea is to represent the computation term as a tree which encodes the computation properties. Our tool then does symbolic computation to simplify the term. This is done by term rewriting it, using rules from arithmetic and the properties encoded in the tree. Fault injections in the computation term are simulated by changing the properties of a subterm, thus impacting the simplification process. An attack success condition is also given and used on the term resulting from the simplification to check whether the corresponding attack works on it. For each possible fault injection, if the attack success condition (which may reference the original term as well as the faulted one) can be simplified to true then the attack is said to work with this fault, otherwise the computation is protected against it. The outputs of our tool are in HTML form: easily readable reports are produced, which contains all the information about the possible fault injections and their outcome.

Computation Term

The computation is expressed in a convenient statement-based input language. This language's Backus Normal Form is given in Fig. 1. A computation term is defined by a list of statements finished by a return statement. Each statement can either:

• declare a variable with no properties (line 3);

• declare a variable which is a prime number (line 4);

• declare a variable by assigning it a value (line 5), in this case the properties of the variable are the properties of the assigned expression;

• perform a verification (line 6).

As can be seen in lines 9 to 15, an expression can be:

• zero, one, or an already declared variable;

• the sum (or difference) of two expressions;

• the product of two expressions; • the exponentiation of an expression by another;

• the modulus of an expression by another.

The condition in a verification can be (lines 18 to 23):

• the equality or inequality of two expressions;

• the equivalence or non-equivalence of two expressions modulo another (lines 20 and 21);

• the conjunction or disjunction of two conditions.

Optionally, variables (when declared using the prime or noprop keywords), expressions, and conditions can be protected (lines 3, 4, 7, and 16, mp stands for "maybe protected") from fault injection by surrounding them with curly braces. This is useful for instance when it is necessary to express the properties of a variable which cannot be faulted in the studied attack model. For example, in CRT-RSA, the definitions of variables d p , d q , and i q are protected because they are seen as input of the computation.

Finally, line 25 gives the regular expression that variable names must match (they start with a letter and then can contain letters, numbers, underscore, and simple quote).

After it is read by our tool, the computation expressed in this input language is transformed into a tree (just like the abstract syntax tree in a compiler). This tree encodes the arithmetical properties of each of the intermediate variable, and thus its dependencies on previous variables. Properties of intermediate variables can be everything that is expressible in the input language. For instance, being null or being the product of other terms (and thus, being a multiple of each of them), are possible properties.

An example of usage can be found in Appx. B, where the original version of Vigilant's countermeasure is written in this language. Note the evident similarity with the pseudocode of Alg. 49.

Fault Injection

A fault injection on an intermediate variable is represented by changing the properties of the subterm (a node and its whole subtree in the tree representing the computation term) that represents it. In the case of a fault which forces at zero, then the whole subterm is replaced by a term which only has the property of being null. In the case of a randomizing fault, by a term which have no properties.

Our tool simulates all the possible fault injections of the attack model it is launched with. The parameters allow to choose:

• how many faults have to be injected (however, the number of tests to be done is impacted by a factorial growth with this parameter, as is the time needed to finish the computation of the proof);

• the type of each fault (randomizing or zeroing);

• if transient faults are possible or if only permanent faults should be performed.

Attack Success Condition

The attack success condition is expressed using the same condition language as presented in Sec. 4.2.1. It can use any variable introduced in the computation term, plus two special variables and @ which are respectively bound to the expression returned by the computation term as given by the user and to the expression returned by the computation with the fault injections. This success condition is checked for each possible faulted computation term.

Simplification Process

The simplification is implemented as a recursive traversal of the term tree, based on pattern-matching. It works just like a naive interpreter would, except it does symbolic computation only, by rewriting the term using rules from arithmetic. Simplifications are carried out in the Z ring, and its Z N subrings. The tool knows how to deal with most of the Z ring axioms:

• the neutral elements (0 for sums, 1 for products);

• the absorbing element (0, for products);

• inverses and opposites;

• associativity and commutativity.

However, it does not implement distributivity as it is not confluent. Associativity is implemented by flattening as much as possible ("removing" all unnecessary parentheses), and commutativity is implemented by applying a stable sorting algorithm on the terms of products or sums. The tool also knows about most of the properties that are particular to Z N subrings and applies them when simplifying a term modulo N :

• identity:

-(a mod N ) mod N = a mod N , -N k mod N = 0; • inverse: -(a mod N ) × (a -1 mod N ) mod N = 1, -(a mod N ) + (-a mod N ) mod N = 0;
• associativity and commutativity:

-(b mod N ) + (a mod N ) mod N = a + b mod N , -(a mod N ) × (b mod N ) mod N = a × b mod N ; • subrings: (a mod N × m) mod N = a mod N .
In addition to those properties a few theorems are implemented to manage more complicated cases where the properties are not enough when conducting symbolic computations:

• Fermat's little theorem;

• its generalization, Euler's theorem;

• Chinese remainder theorem;

• a particular case of the binomial theorem (see Sec. 3.1).

The tool also implements a lemma that is used for modular (in)equality verifications (note that a or b

can be 0): a × x ? ≡ b × x mod N × x is simplified to the equivalent test a ? ≡ b mod N .

Analysis of Vigilant's Countermeasure

This section presents and discusses the results of our formal analysis of Vigilant's countermeasure.

Original and Repaired Version

The original version of Vigilant's countermeasure in our tool's input language can be found in Appx. B.

In a single fault attack model, we found that the countermeasure has a single point of failure: if transient fault are possible, then if p or q are randomized in the computation of N at line 52, then a BellCoRe attack works.

The repaired version by Coron et al. thus does unnecessary modifications, since it fixes three spotted weaknesses (as seen in Sec. 3.2) while there is only one to fix. Only implementing the modification of Alg. 5 is sufficient to be provably protected against the BellCoRe attack.

After that, we went on to see what would happen in a multiple fault model. What we found is the object of the next section.

Our Fixed and Simplified Version

Attacking with multiple faults when at least one of the fault is a zeroing fault, means that some verification can be skipped by having their condition zeroed. We found that two of the tests could be skipped without impacting the protection of the computation (the lines 38 and 57). Both tests have been removed from our fixed and simplified version, which can be found in Appx. C. We have proved that removing any other verification makes the computation vulnerable to single fault injection attacks, thereby proving the minimality of the obtained countermeasure.

Vigilant's original fault model did not include the possibility to inject fault in the conditions of the verifications. If we protect them from fault injection, we found that it is still possible to perform a BellCoRe attack exploiting two faults; Woudenberg et al. [START_REF] Jasper | Practical Optical Fault Injection on Secure Microcontrollers[END_REF] showed that two-fault attacks are realistic. Both faults must be zeroing faults. One on R 3 (resp. R 4 ), and one on S ′ p (resp. S ′ q ). It works because the presence of R 3 (resp. R 4 ) in the last verification is not compensated by its presence in S ′ p (resp. S ′ q ). This shows that even security-oriented optimizations cannot be applied mindlessly in security code, which once again proves the importance of formal methods.

Comparison with Aumüller et al.'s Countermeasure

Now that we have (in Appx. C) a version of Vigilant's countermeasure formally proven resistant to BellCoRe attack, it is interesting to compare it with the only other one in this case, namely Aumüller et al.'s countermeasure. Both countermeasures were developed without the help of formal methods, by trial-and-error engineering, accumulating layers of intermediate computations and verifications to patch weaknesses found at the previous iteration. This was shown, unsurprisingly, to be an ineffective development method. On one hand, there is no guarantee that all weaknesses are eliminated, and there was actually a point of failure left in the original version of Vigilant's countermeasure. On the other hand, there is no guarantee that the countermeasure is minimal. Vigilant's countermeasure could be simplified by formal analysis, we removed two out of its nine verifications. This proves that "visual" code verification as is general practice in many test/certification labs, even for such a concise algorithm is far from easy in practice. Moreover, the difficulty of this task is increased by the multiplicity of the fault model (permanent vs. transient, randomizing vs. zeroing, single or double).

Apart from these saddening similarities in the development process, there are notable distinctions between the two countermeasures. Vigilant's method exposes the secret values p and q at two different places, namely during recombination and reduction from modulo pqr 2 to pq. Generally speaking, it is thus less safe than Aumüller et al.'s method. In terms of computational complexity the two countermeasures are roughly equivalent. Aumüller et al.'s one does 4 exponentiations instead of 2 for Vigilant's, but the 2 additional ones are with 16 bits numbers, which is entirely negligible compared to the prime factors (p and q are 1, 024 bits numbers). Vigilant's countermeasure requires more random numbers, which may be costly but is certainly negligible too compared to the exponentiations.

Conclusions and Perspectives

We have formally proven the resistance of a fixed version of Vigilant's CRT-RSA countermeasure against the BellCoRe fault injection attack. Our research allowed us to remove two out of nine verifications that were useless, thereby simplifying the protected computation of CRT-RSA while keeping it formally proved. Doing so, we believe that we have shown the importance of formal methods in the field of implementation security. Not only for the development of trustable devices, but also as an optimization enabler, both for speed and security. Indeed, Coron et al.'s repaired version of Vigilant's countermeasure included fixes for weaknesses that have been introduced by a hasty speed-oriented optimization, (which is never a good idea in the security field), and two faults attacks are possible only because of random numbers that have been introduced to reinforce the countermeasure.

As a first perspective, we would like to further improve our tool. It is currently only able to inject faults in the data. It would be interesting to be able to take into account faults on instructions such as studied by Heydemann et al. [START_REF] Heydemann | Formal Verification of a Software Countermeasure Against Instruction Skip Attacks[END_REF]. It would be interesting as well to be able to automatically refine some properties of the variables, for instance to study attacks using a chosen message as input of the algorithm. Also, multiple fault analyses can take up to several dozens of minutes to compute. Since each attack is independent, our tool would greatly benefit (2 to 8 times speed-up on any modern multi-core computer) from a parallelization of the computations. Finally, it is worthwhile to note that our tool is ad hoc in the sense that we directly inject the necessary knowledge into its core. It would be interesting to see if general purpose tool such as EasyCrypt [BGZB09] could be a good fit for this kind of work. Injection). Some setups can feature the possibility of multiple fault injections. Simple, double and triple EM pulses are illustrated respectively in Fig. 2, 3 and4.

A.3 Practicality of the Attacks Identified by finja

The pulse amplitude (or "energy") shall be sufficiently high for the EMI to have an effect, and not too large to avoid chip damage. A typical method to set the adequate level of energy is to run the CRT-RSA computations iteratively, while increasing the energy from zero, with small increments (e.g., a few millivolts). After a given number of iterations, the CRT-RSA signature will be faulty or the circuit under attack will refuse to give an answer. From this value on, the EMI starts to have an effect on the chip. Perhaps this fault injection is probabilistic, meaning that it does not always have an effect. So, it is safe to increase the energy level of a few percents to raise the injection probability of success as close to one as possible.

The pulse duration is an important parameter, because it conditions the accuracy in time of the stress. As mentioned in Sec. A.1, there is a great variability in the duration of the operations. Therefore, some attacks will be more difficult to set up successfully. For exemple, the two corrections (Alg. 3 and 5) finally repair a flaw that is fairly difficult to fault, because the operations last:

• 1 clock cycle, i.e., 10 ns, and • 512 clock cycles, i.e., 5.12 µs respectively. Implementations of CRT-RSA are generally working in constant time, to prevent timing attacks [START_REF] Paul | Timing Attacks on Implementations of Diffie-Hellman, RSA, DSS, and Other Systems[END_REF]. Those attacks are indeed to be taken very seriously, as they can even be launched without a physical access to the device. Paradoxically, this helps carry out timely attacks on focused stages in the algorithm. So, with the setup described in Sec. A.2, such fault attacks remain possible. Also, multiple faults can be injected easily if the algorithm timing is deterministic, because the repetition rate of the setup is greater than once per clock cycle. Indeed, a new pulse can be generated every other period of 2 ns, i.e., at a maximal rate of 500 MHz. However, a bench of lower quality would deny some exploits, for instance those that require: (a) to fault an operation of small duration; (b) for high-order fault attacks, operations that are too close one from each other.

  The third fix proposed by Coron et al. consists in transforming lines 42 to 45 of Alg. 49 into: Algorithm 5: Vigilant's CRT-RSA: Coron et al. fix #3 1

First

  , we must define what computation is done, and what is our threat model.

Figure 1 :

 1 Figure 1: BNF of our tool's input language.

Figure 2 :

 2 Figure 2: First-order fault injection attack -proof of concept

Figure 3 :

 3 Figure 3: Second-order fault injection attack -proof of concept We detail below the technical characteristics of a fairly high-end bench 4 :

Figure 4 :

 4 Figure 4: Third-order fault injection attack -proof of concept

We use the usual convention in all mathematical equations, namely that the "mod" operator has the lowest binding precedence, i.e., a × b mod c × d represents the element a × b in Z c×d .

In other papers related to faults, the faulted variables (such as X) are noted either with a star (X * ) or a tilde ( X); in this paper, we use a hat, as it can stretch, hence cover the adequate portion of the variable. For instance, it allows to make an unambiguous difference between a faulted data raised at some power and a fault on a data raised at a given power (contrast X e with X e ).

http://pablo.rauzy.name/sensi/finja.html

We credit Laurent Sauvage for this information, that comes from the specification of one plug-in of the Smart-SIC Analyzer tool of Secure-IC S.A.S.

Acknowledgements

The authors wish to thank David Vigilant for his insightful answers to our questions. We also sincerely acknowledge the SIGPLAN PAC program of the ACM for a financial sponsor.

Algorithm 2: Vigilant's CRT-RSA

Input

: Message M , key (p, q, dp, dq, iq) Output : Signature M d mod N , or error 1 Choose random numbers r, R 1 , R 2 , R 3 , and R 4 .

A Practical Feasibility of the Identified Attacks

In this section, we discuss the key features that an attacker must control for his attacks to succeed. We first describe in Sec. A.1 the timing characteristics of a CRT-RSA computation implemented on a smartcard. Second, we analyze in Sec. A.2 the performance of a fault injection bench. Third, we analyse in Sec. A.3 which fault injection is realistic or, at the contrary, challenging practically speaking.

A.1 CRT-RSA Algorithm in a Smartcard

We consider a smartcard implementation of CRT-RSA. It consists in an arithmetic hardware accelerator, capable of computing multiplications and additions. It is thus termed a MAC (Multiplier-Accumulator Circuit). This MAC can manipulate data on 64 bits. In 65 nm and lower technologies, such a MAC can run at the frequency of 100 MHz (after logical synthesis in a low power standard cell library). Let us consider that N is a 2048 bits number, i.e., p and q hold on 1024 bits. The computations are carried out in the so-called Montgomery representation. For instance, when computing a modular multiplication, the multiplication and the reductions are interleaved. Thus, we derive some approximate timings for the unitary operations.

• One schoolbook addition on Z p (or equivalently on Z q ) of two large numbers takes 1024/64 = 16 clock cycles, i.e., 160 ns;

• One schoolbook squaring or multiplication on Z p of two large numbers takes 2 × (1024/64) 2 = 512 clock cycles, i.e., 5.12 µs; The factor two takes into the account the reduction steps.

• One schoolbook exponentiation on Z p of one large number by another large number takes 512 × (1024 × 3 /2) = 786432 clock cycles, i.e., about 7.8 ms; This estimation assumes the square-and-multiply exponentiation algorithm, in which there are as many squarings as bits in p, but only half (in average) multiplications.

Besides, we notice that some operations in Alg. 49 are hybrid: they involve one large number and one small number (e.g., 64 bits number, of the size of r, R 1 , R 2 , R 3 or R 4 ). Obviously, those operations unfold more quickly than operations with two large numbers. Some operations in Alg. 3 can even be implemented to be extremely fast. For instance, the subtraction on line 1 or 4 of Alg. 3 (resp. Alg. 4) can be imagined to be done in one clock cycle, i.e., 10 ns. Indeed, as p and q are odd (because they are large prime numbers), this "decrementation" consists simply in resetting the least significant bit of the operand p (resp. q). So, concluding, the operations in the CRT-RSA algorithm can be of diverse durations, ranging 6 orders of magnitude, from 10 ns for operations computable in one clock cycle, to about 8 ms for one modular exponentiation.

A.2 Fault Injection Setup

An efficient fault injection method consists in the production of a transient electromagnetic (EM) pulse of high energy in the vicinity of the circuit that runs the CRT-RSA. This kind of fault is termed EMI (EM d'q := dq + R2 * (q -1) ; 41 Sqr := M'q^d'q mod q' ; 42 43 if d'q !=[q -1] dq abort with error ; ---recombination ---51 S := S'q + q * (iq * (S'p -S'q) mod p') ; S := S'q + q * (iq * (S'p -S'q) mod p') ; _ != @ /\ ( _ =[p] @ \/ _ =[q] @ )