André Galligo

Deformation of Roots of Polynomials via Fractional Derivatives

Introduction

Fractional calculus has attracted interest of engineers and researchers in several fundamental or applied fields, but it is often ignored by the researchers in Polynomial algebra and too seldom encountered in Computer algebra. In this paper we emphasize that Fractional calculus is based on simple rules which can be exploited in Computer algebra, then we build a novel geometric object associated to a real univariate polynomial. Fractional derivatives define a deformation of a univariate polynomial: the order of differentiation is the continuous parameter of deformation. This construction appeared naturally in our study of Budan tables (a 2 D diagram where the vertical axis represents the integer orders of differentiation) of polynomials with large degrees when we replace a discretized curve by a continuous one, see [START_REF] Bembe | A: Virtual roots of real polynomials and fractional derivatives[END_REF], ( 6), [START_REF] Galligo | Roots of the derivatives of some random polynomials[END_REF], and Section 4.

For a fixed polynomial f (x), we consider the order of differentiation q as a new indeterminate then factor the fractional derivative to get an associated bivariate polynomial that we denote by P f (x, q). For 0 ≤ q ≤ 1, P f (x, q) defines an homotopy between the polynomials f (x) and xf (x).

We associate to a degree n real polynomial f a plane spline curve, the stem of f . It is the union of n-1 real zero sets defined by P f (m) (x, q) = 0, 0 ≤ q ≤ 1, 0 ≤ m ≤ n -2 we called the layers of the stem. We classify the admissible configurations of layers into 4 types and proved that they are realizable by providing examples; they complement Rolle's theorem on the respective location of the real roots of a polynomial and of its derivative. Type 1 is the simplest and more common one while type 4 is more unusual. Stems of classical random polynomials exhibit regular patterns. Relying on many computational experiments we proposed conjectures to describe their geometry. In particular, we conjecture that in the layers of the considered random polynomials, almost surely types 3 and 4 never happen. The analysis of these patterns led us to extend the deformation process to the complex values of x. We observed experimentally that in the complex plane, P f (x, q) creates a correspondence between the complex roots and the critical points of our random polynomials f (x). It seems that this kind of pairing has not been reported before. In order to explain it, we consider a classical electrostatic interpretation: a critical point of f is viewed as an equilibrium of a particle submitted to a field of forces; and we suggest a mean field approach based on limit symmetries of the roots sets. However the analysis becomes complicated since there are two conflicting limit symmetries and we were not able to prove our claim, so we left it as a conjecture. In the conclusion we indicate several potential extensions of our constructions. We relied on the computer algebra system Maple to compute examples, and create pictures which illustrate our presentation.

The paper is organized as follows. Section 2 provides a brief introduction to Fractional calculus: history, definitions and properties. Sections 3 and 4 present our polynomial factor P f (x, q) of the q-th fractional derivative of a monic polynomial f , then describe two kinds of plane curves (FD-curve and stem) attached to f . Section 5 applies our constructions to some classical random polynomials, and set conjectures tested on many computational experiments. Section 6 concentrates on the homotopy between f (x) and xf (x) defined when 0 ≤ q ≤ 1; and on the correspondence between the complex roots of f (x) and xf (x). Then we sketch our mean field approach relying on limit symmetries of roots sets.

Notations:

We denote by R the field of real numbers, by f a monic polynomial of degree n, f := x n + i=n-1 i=0 a i x i ; and by f (m) the m-th derivative of f for an integer m, with f (0) := f .

Fractional derivatives

The attempt to introduce and compute with derivatives or anti derivatives of non integer orders goes back to the 17-th century. The traditional adjective "fractional", corresponding to the order of differentiation, is misleading since it need not be rational. An historical progression of the concept from 1695 to 1975, through a hundred citations of mathematicians, is provided in [START_REF] Oldham | The fractional calculus[END_REF]. More recently, fractional derivatives has attracted a renewed interest of engeniers and researchers in several fundamental or applied fields: applied mathematics, viscoelasticity, system theory, physics, chemistry, and even finance and social sciences. The paper [START_REF] Machado | Recent history of fractional calculus[END_REF] reports hundreds of specialized conferences, monographes, articles and softwares in the area of Fractional calculus during the last 40 years.

Let us emphasize that nowadays in Mathematics, fractional derivatives are mostly used in Functional analysis for the study of PDEs. A classical example is the formal factorization of the heat equation using the semi differentiation S such that SS = ∂ ∂t . Fractional derivatives are presented via Fourier or Laplace transforms, which are well adapted to the functions exponentials, sinus and their infinite sums. Moreover Fractional calculus can be considered as a laboratory for creating new representations of special functions and integral transforms. Unfortunately, Fractional derivatives are often ignored by the researchers in Polynomial algebra and too seldom encountered in Computer algebra.

An important property which remains true in Fractional calculus is that two fractional differentiations commute. However, for non-integer orders of differentiation, the fractional derivative at a point x of a function f does not only depend on the graph of f very near x; fractional differentiations do not commute with the translations on the variable x. Moreover, like for integrals a fractional derivative depends on the choice of an initial point.

For applications in engineering and science, the order of differentiation of fractional derivatives is generally either fixed or limited to few values; it will not be the case here since the order of differentiation will serve as a continuous parameter of deformation. This interpretation arrived naturally in our study of Budan tables of univariate polynomials with large degrees: The vertical axis represented the integer orders of differentiation, and we wanted to replace a discretized curve by a continuous one, see (1), ( 6), [START_REF] Galligo | Roots of the derivatives of some random polynomials[END_REF], and our Section 4. Notice that the idea of deforming a function via its fractional derivatives is not original since it already appeared in [START_REF] Larson | Continuous deformationof a polynomial into its derivatives[END_REF], and [START_REF] Miller | Insight into the Fractional Calculus via a Spreadsheet[END_REF]; the authors were motivated by mathematics education.

Several definitions

Fractional calculus can be introduced "pedagogically" by giving 3 easy formulae which apply respectively to the finite sums of monomials, exponentials and sinus. Each of them corresponds to the extension to the real orders of a usual formula of differentiation with integral orders.

Diff q denotes the fractional differentiation of order q ≥ 0. We have Diff

0 (f ) = f , Diff q (Diff r (f )) = Diff q+r (f ) and Diff 1 (f ) = f .
In the next sections, we will concentrate on the first one, attributed to Lacroix (1815) or to Peacock (1833), since our interest is in polynomial algebra.

(1) For monomials:

Diff q (x n ) := n! (n -q)!
x n-q ; f or q ≥ 0, n integer.

Where (n -q)! = Γ(n -q + 1); Γ is the gamma function, so:

Dif f 1/2 (x 2 -2x + 3) = ( 8 3 x 2 -4x + 3)x -1/2 1 √ π .
(2) For exponentials:

Diff q (e sx ) := s q e sx ; for q ≥ 0, s real positive.

(3) For sinus:

Diff q (sin(sx)) := s q sin(sx+q π 2 ); for q ≥ 0, s real positive.

An important task of Fractional calculus was finding a unique framework to interpret all special formulae or procedures. It was achieved by considering generalized functions (such as distributions) together with Laplace or Fourier transforms, and then rely on integrations (which correspond to fractional derivatives with negative orders) to extend the domain of validity of the proposed formulae. The theory now applies widely in Functional analysis and for numerical approximate computations. The general definition expresses clearly the fact that a fractional derivative is a non local concept. Hence Fractional calculus is necessarily more complicated than usual calculus.

The integral formulation shows that a good definition of fractional derivative of a function requires to remove one indetermination, e.g by choosing a lower limit (or origin) a ∈ R ; the corresponding notation is Diff q a . For instance in our first formula this origin is a = 0, in the second it is a = -∞, and in the third it has a more complicated expression since the function is periodic.

An unfortunate consequence of this generalization process is that fractional derivatives of polynomials which is elementary using the first above formula, seems very complicated through the computation of Fourier or other transforms using indefinite integrals.

In the sequel, we will only consider fractional derivatives with q ≥ 0 of polynomials and use the above formula.

Bivariate polynomial associated to f

A simple construction

We will consider a polynomial factor of the fractional derivatives of the polynomial f . We rely on Lacroix-Peacock rule for monomials, see above. Lemma 1. Let f (x) be a polynomial of degree n, then

P f (x, q) := 1 n! x q Γ(n -q)Diff q (f )
is a polynomial in x and q, of total degree n.

Notice that, for an initial value a ∈ R, the polynomial can be written in the basis formed by the powers of (x -a), then one considers the fractional derivative operator Diff q a and the corresponding polynomial P f,a (x, q). Here we set a = 0.

We decided that, since we focus on the polynomial roots of a fixed univariate polynomial f , it is convenient to choose for a the mean of the complex roots of f . Equivalently, thanks to the famous Tchirnaus transformation, after a real translation on the abscissa x we can assume that the coefficient of x n-1 vanishes, hence the mean of the roots of f is 0. This last property is also true for all derivatives with integer orders of f . To simplify the presentation of our study we will also assume that the other coefficients of f are not 0, a property generically satisfied. So we will consider only q ≥ 0 and concentrate on the following set. Definition 1. • We denote by E n , the subclass of monic polynomials f in R[x] of degree n such that:

f = x n + n-2 i=0 a i x i , a i ∈ R, a i = 0. • Let f ∈ E n .
We call bivariate polynomial associated to f , the polynomial P f (x, q) := x q (n-q)! n! Diff q (f, x)). It is a polynomial of total degree n in x and q, which may be written:

P f (x, q) := x n + n-2 i=0 ( n j=i+1 (1 - q j ))a i x i .
The polynomial P f (x, q) has the following property.

Lemma 2. Let f ∈ E n . Then, P f (0, q) = 0 if and only if q is an integer between 1 and n.

Hence it makes sense to consider the deformations of the roots of f (or of its derivatives) defined by P f (x, q) = 0 separately in the two half (x, q)-planes, x ≥ 0 and x ≤ 0.

For x > 0 and q > 0, this amounts to consider as in [START_REF] Miller | Insight into the Fractional Calculus via a Spreadsheet[END_REF] the deformation of (the graph of) f by fractional derivatives but skip fractional powers of x.

Remark 1. P f (x, q), considered as a polynomial in x, is an Hadamard product of two univariate polynomials, f (x) and a fixed polynomial S(q, x) defined by:

S(q, x) := x n + n-1 i=0 ( n j=i+1 1 - q j )   n i   x i .
Roots of Hadamard products have nice properties, see [START_REF] Rahman | Analytic theory of polynomials[END_REF] section 3.5 and chapter 5.

FD-curves

Definition 2. We call FD-curve of a polynomial f ∈ E n the real algebraic plane curve defined by the bivariate associated equation P f (x, q) = 0.

Figure 1 shows an example with f := (x -1)(x -2)(x -3)(x + 0.5)(x + 2.5)(x + 3), n = 6, an hyperbolic polynomial. The roots of f and its derivatives are represented by small green disks. Figure 2 shows another simple example with f := (x -2)(x -3)(x + 0.5)(x + 2.5)(x 2 + 1), n = 6.

FD-curves provide an interesting family of algebraic plane curves: each one is "naturally" attached to a univariate polynomial f and interpolates the real roots of the derivatives of f , hence gives geometric insights on the correlations between these roots. However we found that these correlations appear more neatly on another curve associated to f , the stem of f .

Budan tables and stems

Definitions

For a polynomial f (x) in E n , the real roots and the signs of all its derivatives with integer order are collected in a 2D diagram we called, in ( 1) and ( 6), a Budan table. The second coordinate indicates the degree of the derivative. Notice that if all roots are simple then in each row, the signs change at each root and the representation can be sketched to only indicate pairs formed by a root of a derivative and the degree of that derivative.

Example:

f := (x -5).(x 2 -x + 4) ; f = 3(x -1).(x -3) ; f = 6(x -2).
These polynomials have respectively 1, 2, 1 real roots, the Budan table BT can be represented by:

BT = {[5, 3], [1, 2], [3, 2], [2, 1]}.
Figures 3 and4 show Budan tables of two polynomials of degrees 64. Definition 3. Let f be a polynomial in E n . The stem of f is a C 0 spline formed by the union of n -1 portions of real curves, we called layers: For i from 0 to n -2 each layer is the zero set, of degree n -i, defined by P f (i) (x, q) = 0 and 0 ≤ q ≤ 1.

Stems were designed to study the roots of the derivatives of random polynomials of high degrees and exploit their symmetries.

To illustrate the differences between FD-curves and stems, Figures 5,6 show the stems corresponding to the two previous FD-curves: stems are much less curved and there is a C 0 continuity between the adjacent layers of a stem. However stems provide a good approximation of a Budan table when n is large and f is a random polynomial, compare Figures 7 and8 which correspond to polynomials of degree 128.

Configurations for

0 ≤ q ≤ 1
An important point in the study of stems is the analysis of the configurations of segments of curves (and ovals) appearing in one layer (i.e. between two integer values of the vertical coordinate); without loss of generality, in the deformations of the roots of f (x) into the roots of xf (x) when the parameter q in P f takes all the values between 0 and 1. To simplify the description of the configurations and since we will apply it to random polynomials, we will assume that f (x) and xf (x) have simple roots and that the curve C (the layer) defined by P f (x, q) = 0, 0 ≤ q ≤ 1 is smooth. This condition is generically true since it is implied by the non vanishing of the bivariate resultant Res (x,q) (P f , (P f ) x , (P f ) q ). Rolle's theorem implies that between two roots of f there are an odd number of roots of f . But we do not know a priori how they are connected through C. Proposition 1. With these assumptions, there are 4 types of configurations which happen in a layer C of the stem of f . They are enumerated here after.

(1) f has r real roots, X 1 < ... < X r and xf has also r real roots Y 1 < ... < Y r and each X i is connected, through C by a segment of curve, to Y i , for 1 ≤ i ≤ r.

(2) f has r real roots, X 1 < ... < X r and xf has r + 2u real roots Y 1 < ... < Y r+2u and each X i is connected, through C by a segment of curve, to some Y θ(i) , for 1 ≤ i ≤ r, with 1 ≤ θ(i) < θ(i + 1) ≤ r + u. The others Y j , not equal to any Y θ(i) , are connected pairwise through C by segments of curves which do not touch the line q = 0.

(3) f has r + 2v real roots, X 1 < ... < X r+2v and xf has r + 2u real roots Y 1 < ... < Y r+2u ; r roots of f , X j 1 < ... < X j r , are connected, through C by a segment of curve, to some Y θ(i) , for 1 ≤ i ≤ r, with 1 ≤ θ(i) < θ(i + 1) ≤ r + u. The other roots of f (x), not equal to any X j i , are connected pairwise through C by segments of curves which do not touch the line q = 1. And the others roots of xf (x), not equal to any Y θ(i) , are connected pairwise through C by segments of curves which do not touch the line q = 0.

(4) There is a loop in C which does not touch the line q = 0 nor the line q = 1.

It is enough to show that these configurations are effectively realized.

Example 2.

(1) All layers in Figures 5 and6 are of Type 1.

(2) Type 2 happens for all layers in Figures 8, where the curve admits an horizontal tangent. Another example is given in

Gaussian Random polynomials

Random polynomials is a classic and active subject in Mathematics and Statistics. It is at the core of extensive recent research and has also many applications in Physics and Economics; two books (2) and ( 5) are dedicated to it. Already in 1943, Mark Kac (9) gave an explicit formula for the expectation of the number of roots of a polynomial in a class that now bears his name (see below). The subject is naturally related to the study of eigenvalues of random matrices with its applications in Physics, see e.g. ( 4).

Here, for a fixed degree n, we consider a polynomial f in E n given in the monomial basis:

f := x n + n-2 i=0 a i x i
where the coefficients a i are instances of n -1 independent centered normal distributions N (0, σ). We are concerned with averaged asymptotic behaviors when n tends to infinity.

We consider 3 classes of random polynomials, specify their names and the corresponding variances σ 2 (i) of the coefficients

a i . • Kac polynomials: σ 2 (i) = 1. • SO(2)-polynomials: σ 2 (i) = n i . • Weyl-polynomials: σ 2 (i) = 1 i! .
The limit distributions when n tends to infinity, of the complex roots of these random polynomials is almost uniform in angles around the origin, this property is completed by an axial symmetry over the real axis due to complex conjugation. This observation can be quantified: (16) computed for Kac polynomials the density function h n (x, y) of the number of complex roots near a complex point x + iy, completing Kac's computation of the density function of the number of real roots near a real point x.

Our first experiments with Kac polynomials found that their Budan tables admit unexpected structured patterns, see Figure 3. The roots of the successive derivatives present almost dotted curves and alignments. To our best knowledge, this phenomenon has not been explored before.

We made more experiments with different instances of Kac polynomials and got very similar patterns, then we repeated the experiments with the different random polynomials defined above, see Figures 5 and6.

Based on our experiments, we propose the following conjectures on probabilistic behaviors when n tends to infinity.

Conjecture 1.

Let f be a Kac, SO [START_REF] Bharucha-Reid | Random Polynomials[END_REF] or Weyl random polynomial, as above. If the complex root with maximal module of f is real positive, then almost surely all the maximal real roots of f (m) for 0 ≤ m ≤ n -k with k << n tend to be aligned.

A similar statement for negative real roots is obtained replacing f (x) by (-1) n f (-x) .

Conjecture 2. Let f be a Kac, SO [START_REF] Bharucha-Reid | Random Polynomials[END_REF] or Weyl random polynomial, as above. In the stem of f almost surely all the layers are of types 1 or 2 (described in Lemma 4.1). The dotted line in Figure 11 suggests to examine the corresponding deformation of complex roots.

Complex critical points

There is an important bibliography on the location of the critical points of a polynomial (i.e. the roots of its derivative f ) with respect to the location of its roots, going back to Gauss with Gauss-Lucas theorem. Several recent works concentrate on the following conjecture of Sendov, which has been proved for small degrees and in several special cases. Their proofs rely either on the implicit function theorem or on extremal polynomials, or on refinements of Gauss-Lucas theorem. See the book (15), chapter 7.

Sendov Conjecture: Let f be a polynomial having all its roots in the disk D. If z is a root of f , then the disk z + D contains a root of f .

Observations

We performed experiments with the classes of random polynomials described in section 5, they exhibited interesting behaviors:

-for almost each root of f smaller disks z + D, with << 1, contain a critical point; in such a way that they suggest a correspondence between the roots of f and the roots of xf , -one can restrict these disks to small sectors, which indicates a direction towards the real axis or towards the origin.

Figures 13 to 14 illustrate the relation between roots and critical points for an SO(2) random polynomial of degree 32. Fractional derivatives are used to construct (as for real roots) a deformation between the zero sets of f and xf . The color chart is: the roots of f are blue, the roots of f are red and the roots of the fractional derivatives are green. Figure 15 shows the "top" part of a complex analog of the Budan table, notice the regular alignments towards the origin (it is slightly perturbed near the real axis).

Complex stem

By analogy with section 4.1 we attach to a polynomial f ∈ E n the C 0 spline formed by the n -1 layers γ i for 0 ≤ i ≤ n -2 in the (x, q)-space C × [0, 1]. Each layer γ i is a subset of a 3D algebraic curve of degree n -i, its "real" part C i is one of the layer of the (real) stem defined in section 4.1.

Since f has real coefficients, generically, the only double points of each layer γ i appear on C i , otherwise due to complex conjugation they would correspond to pairs of double points with the same value of q. These double points correspond to points of C i with an horizontal tangent. Hence, generically, the complex stem of f defines a continuation process which allows to distinguish the paths of each complex non real root of f until it arrives on the real axis.

Deformation process in the random setting

We view the non integer differentiation order q as a time parameter. For random polynomials all generically true conditions are almost surely satisfied.

Consider a random polynomial f as above and assume Conjecture 2 of section 5. Then when q varies, almost surely the deformation process in the complex plane is obtained when the roots of the derivatives of f move according to the following rules.

• The real roots of f (x) are deformed to real roots of xf (x) as described by a layer of type 1. • Some pairs of conjugated complex roots of f are deformed to real double roots of a fractional derivative and then to pairs of real roots of xf (x) as described by a layer of type 2. • The other complex non real roots of f (x) are deformed to complex non real roots of f (x). • The process is repeated inductively replacing f by each derivative f (m) for m = 0, ..., m -2.

Electrostatic attraction

The interpretation of the location of each critical point of f as an equilibrium of a logarithmic potential, where the roots of f are viewed as positively charged particles (or rods), goes back to F. Gauss. As reported in [START_REF] Rahman | Analytic theory of polynomials[END_REF], the following equality to zero provides a quick proof of Gauss-Lucas theorem.

Denote by x j the complex roots of the polynomial f assumed distinct from each other and distinct from z, another complex number. Then by logarithmic differentiation and conjugation we deduce:

f (z) = 0 ⇒ z -x j |z -x j | 2 = 0.
Thus a critical point z is the barycenter of the roots of f with positive weights 1 |z-x j | 2 , hence is in the convex hull of the set of roots.

We now sketch a mean field approach to explain the observed phenomena.

The vector in the complex plane -z+x j |z-x j | 2 is viewed as a force applied to z directed towards x j proportional to the reciprocal of the distance. Summing these forces, the point z (viewed as an electron) is attracted by the roots system of f (viewed as positively charged particles). When the limit distribution of the roots of f is uniform in angles, i.e. only depends on the radius, the resulting electrostatic force on a point z inherits a limit symmetry, and tends to be directed towards the origin.

Let us denote by L 1 the real line joining the origin to a root x k and by L 2 the real line orthogonal to L 1 through the origin. The number and distribution of roots below and above L 1 (respectively L 2 ) are asymptotically "almost" balanced. So, with a good probability, an equilibrium z k can be found "near" L 1 with the vector x k z k oriented towards the origin. One can expect as well that x k the farther is from the origin, the smaller the vector x k z k should be. This is what we observed in our experiments as illustrated with Figures 13.

The previous balanced count of forces is perturbed when we approach the real axis, because there is another axial symmetry due to complex conjugation, and a positive probability of real roots. This breaks the rotational symmetry, consequently the resulting electrostatic force is now also directed towards the real axis. Conjecture 3. For f a Kac, SO [START_REF] Bharucha-Reid | Random Polynomials[END_REF] or Weyl random polynomial, the presented continuation process realizes a bijection such that each non real critical point z k of f is attached to a root x k . Moreover in the limit distribution of (x k , z k ) when n tends to infinity, almost surely the vectors x k z k point towards the origin.

Remark 2. Conjecture 3implies that the zero set of f resembles the zero set of f but "shrunk". Hence the extremal module root of f will also be positive real. Then a simple computation will give the deplacement of the extremal root and inductively prove Conjecture 1.

Conclusion

In this article we considered deformations of the roots of a univariate polynomial f via a "natural" homotopy computed with a bivariate polynomial obtained as a factor of the fractional derivatives of f . We introduced geometric objects related to that construction, and initiated their study. We reported computational experiments performed with f in some classes of Gaussian random polynomials, then proposed 3 conjectures. Although we also suggested some insight and a mean field approach to better understand the underlying geometry, the proof of these conjectures will require advanced developments in Probability theory. We now list 3 other directions of research.

• Adapt our approach to fractional derivatives of trigonometric polynomials and study their roots. • Extend in 2 variables and experiment with deformations of algebraic plane curves.

• In this paper we put some limitations on the considered polynomial f , e.g. having distinct roots, study the more general case.

Fig. 1 .Fig. 2 .

 12 Fig. 1. A FD-curve

Fig. 3 .Fig. 4 .

 34 Fig. 3. Budan table of a Kac polynomial of degree 64

Fig. 5 . 1 Fig. 6 .

 516 Fig. 5. Stem of Figure 1

Fig. 7 .Fig. 8 .

 78 Fig. 7. Budan table of a SO(2) polynomial (detail)

Fig. 9 .Figure 9 ,

 99 Fig. 9. a type 2 layer

Figures 3 ,

 3 Figures 3, 5 and 6 illustrate these 2 conjectures.

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Roots of f and f

Preprint submitted to Elsevier

Acknowledgements

The author was partially supported by the European ITN Marie Curie network SAGA.


V Drean

G Perrot

G Couegnat

M Ricchiuto

G L Vignoles

Thème Num

Projet Bacchus

Image-based 2D numerical modeling of oxide formation in self-healing CMCs

Keywords: composite materials, self-healing materials, image based modelling, transversal crack model, liquid oxide formation, finite element discretization, oxidation mechanism

Excellent lifetimes make Self-healing Ceramic-Matrix Composites (CMCs) are promising candidates for jet engine hot parts. These composites are constituted of a 3D arrangement of SiC fiber tows infiltrated by a multilayer matrix. A pyrocarbon interphase acts as a crack deviator, SiC matrix layers bring stiffness, and boron-containing phases are able to produce above 450°C a liquid oxide preventing further oxidation by a diffusion barrier effect. This paper introduces an image-based numerical simulation of the self-healing mechanism under oxygen gas. Existing 0D or 1D models give the time evolution of the oxygen concentration in the weakest fiber and deduce from it a global lifetime through an oxygen-controlled slow crack growth rate law. We propose an approach in which the resolution domain is a 2D FE mesh directly obtained from transverse images of a tow containing the crack. Oxygen diffusion, carbon consumption around the fibers, and conversion of boron carbide into boron oxide are simulated. The model solves mass balances in terms of heights of oxygen (gaseous or dissolved), liquid oxide, pyrocarbon, and boron-containing phase. All the heights are considered perpendicular to the image plane (thin layer approximation), and represent the unit-volume (per square meter) occupied by each phase. Preliminary results on images containing several dozens of fibers and a multilayer matrix are discussed.

INRIA Research Report model

V. Drean3 , G. Perrot4 , G. Couegnat 2 , M. Ricchiuto 1 , G.L. Vignoles 2 Thème NUM Projet BACCHUS Rapport de recherche n° 7845 -Decembre 2011 -10 pages Résumé: Une importante durée de vie fait des composites à matrice céramique auto cicatrisante (CMCs) les candidats parfaits pour utilisation dans les composantes à haute température des moteur aéronautiques du future. Ce type de composites sont constitués par des rangement 3D de fibres de SiC infiltrées par une matrice multicouche. Une interphase de pyrocarbone protège les fibres en déviant la fissure, des couches de SiC dans la matrice la rendent plus rigide, et des phases qui contiennent du bore produisent à 450°C un oxyde liquide qui protège les fibres par un effet de barrière diffusive. Ce papier présente une modélisation numérique de ce phénomène de auto cicatrisation en présence d'oxygène. Les modèles 0D et 1D qui existent en littérature donnent une approximation de l'évolution de la concentration d'oxygène en correspondance de la fibre la plus faible et permettent d'en déduire une durée de vie globale via un mécanisme de croissance de la fissure contrôlée par l'oxygène. On propose une approche dans laquelle le domaine de calcul est un maillage éléments finis 2D, obtenu directement par des images de fils fissurés. La diffusion de l'oxygène, l'oxydation du pyrocarbone et la production d'oxyde liquide sont simulés. Le modèle résout des bilans de masse en termes de hauteur d'oxygène, d'oxyde de bore et de pyrocarbone. Ces hauteurs sont considérées comme étant perpendiculaires à la fissure et représentent le volume unitaire (par unité de surface) occupé par chaque phase. Résultats préliminaires sur des images qui contiennent des dizaines de fibres et une matrice multicouche sont présentés.

Mots clés: matériaux composites, matériaux auto cicatrisants, modélisation par images, modèles de fissures transverses, formation oxyde liquide, discrétisation éléments finis, mécanismes d'oxydation

Introduction

Self-healing Ceramic-Matrix Composites (SH-CMCs) are promising candidate materials for civil aircraft jet engine hot parts 1 . Tests in actual conditions of use have been performed, showing excellent lifetimes 2 . These materials are constituted of a woven/interlocked 3D arrangement of SiC fiber tows, infiltrated by a multilayer matrix. All the interphase and matrix layers have a specific role. First, there is a pyrocarbon interphase acting as a crack deviator and preventing premature fiber failure while the matrix undergoes multiple cracking, that is, progressive damage. Second, SiC matrix layers bring stiffness and are chemically inert at moderate temperatures. Finally, boron-containing phases produce above 450°C a liquid oxide which fills cracks and prevents further oxidation by a diffusion barrier effect 3 , thus increasing by large amounts the material lifetime under exposition to high-temperature oxidation 4 . The massive production of SH-CMC parts in engines can only be envisaged if there is a sufficient confidence in the material lifetime duration. The trouble is that the material can last as long as ten years; hence, a material development cycle based only on experimental characterization is not feasible. Modeling is mandatory in order to incorporate as much understanding of the material's physics as possible, and identify and explain the experiments that have been carried out in "accelerated degradation" conditions, or on material parts. The present work has been developed in this frame. The material lifetime results from a balance between degradation and self-healing 5 . In the SH-CMC, the network of matrix cracks brings the oxygen in contact with all the inner material constituents; they are transformed gradually into oxides, some of them being gaseous (e.g. carbon), liquid (e.g. boron) or solid (e.g. silica at moderate temperatures). The calculation of the chemical reaction rate everywhere in the material is a key point for lifetime prediction. Past reports have described models which give the time evolution of the oxygen concentration in, say, some key part of the composite (i.e. the weakest, most exposed fiber); they deduce from this quantity a global lifetime because the tensile strength of this fiber can be related to oxygen exposure through an oxygen-dependent slow crack growth rate law 6,7 . These models were developed in 0 8,9 or 1 [START_REF] Rebillat | Procs 5 th Intl. Conf on High Temperature Ceramic Matrix Composites (HT-CMC5)[END_REF] space dimensions. In order to gain further insight in the self-healing phenomenon at fiber and tow scale, we propose now an image-based approach, in which the resolution domain is a 2D FE mesh directly obtained from original mesh generation software. Though the whole history of the material includes liquid phase flow, we will deliberately discard this phenomenon here and focus on the first steps of the material evolution, that is, the formation of the liquid oxide and the plugging of the crack. The paper is organized as follows. First, we will recall the main physicochemical phenomena that are being accounted for in the model; then, model equations are built, and an analysis of their stiffness is presented. The next part describes the numerical implementation. Finally, some results are presented and discussed.

Material structure and physico-chemical processes

The SH-CMCs we model are depicted on fig. 1a. Several domains are easily recognized:

1. The SiC fibers, which are sensitive to oxidation through a subcritical crack growth behavior. Their diameter is in the range 8-15 µm.

2. The thin pyrocarbon (PyC) interphase coating the fibers. This layer acts a crack deviator, helping the development of multi-cracking, which gives the material its exceptional mechanical behavior [START_REF] Dréan | [END_REF]12 . In this work we will consider that cracks are perpendicular to the fiber and are limited by the fiber boundary. In particular, fig. 1b illustrates the fact that in our model the fibers are bridging the cracks.

3. The multi-layer matrix containing two constituents: silicon carbide, which is nearly inert at the temperatures typical of civil aircraft applications (it begins to oxidize appreciably above 1000°C 13,14 ), and boron carbide, yielding a liquid (glassy) phase upon oxidation 3,4 . INRIA 4. A number of cusp-shaped pores left by the production process. In particular, the PyC, SiC and boron carbide layers are deposited by Chemical Vapor Infiltration (CVI) 15,16 , a process that gives characteristic continuous layers growing in conformal shape with respect to the initial substrate (i.e. the fibers). The layer thicknesses are larger in the outer parts of the fiber bundles or tows, because of the competition between chemical deposition reactions and gas diffusion 17 . As a result, some pores remain between the fibers after infiltration. These imperfections may or may not be connected to the outer gas phase, depending on whether the crack network intercepts them.

Concerning the material evolution, the dominating phenomenon is oxidation. Accordingly, oxygen balance is the natural ingredient of the model, as well as balances for all oxidized materials, namely, pyrocarbon -undergoing active oxidation and giving gaseous CO 2 -and boron carbide, yielding a liquid B 2 O 3 oxide together with dissolved CO 2 . Silicon carbide may also produce a solid silica layer upon oxidation; however, the temperatures considered here (< 650°C) are not high enough to allow the development of a significant silica thickness. Moreover, here we only address the case in which the boron oxide liquid is viscous enough not to flow significantly away from the position where it has been created.

Ω f Ω i ∂Ω f Ω r Ω i ∂Ω e 2h p 2e f h l h m z r c) d) z r Fiber Inert matrix Inert matrix Glass- forming matrix Air Interphase Ω f Ω i ∂Ω f Ω r Ω i ∂Ω e 2h p 2e f h l h m z r Fiber Inert matrix Inert matrix Glass- forming matrix Air Interphase Ω f Ω i ∂Ω f Ω r Ω i ∂Ω e 2h p 2e f h l h m Fiber Inert matrix Inert matrix Glass- forming matrix Air Interphase Ω f Ω i ∂Ω f Ω r Ω i ∂Ω e 2h p 2e f h l h m z r c) d)
RR n° 7845

Model setup

In order to describe the formation of the protective layer, oxygen diffusion, carbon consumption around the fibers, and conversion of boron carbide into boron oxide are simulated. The underpinning model consists of mass balances of oxygen (gaseous or dissolved), boron carbide, boron oxide, and pyrocarbon. The liquid oxide obtained as a product is assumed viscous enough not to flow away at least in the initial phase of the healing process simulated here. As a consequence, its growth is vertical above the carbide/oxide interface (see fig. 1d). Using a thin layer approximation, the mass balances are written in terms of heights representing the unit-volume (per square meter) of certain face present locally. All the heights are considered perpendicular to the plane of the image, assumed to be the crack plane. In particular, the simulation domain consists of several subsets. As depicted by figs. 1c and 1d we can distinguish the following subdomains: -the inert matrix domain (silicon carbide) W i ; -the reactive matrix domain (boron carbide) W r ; -the fibers boundary ∂W f , modeled as an infinitely thin one-dimensional reactive layer of pyrocarbon ; -the pore boundaries ∂W e , subjected to exposure to external oxygen ; -the fibers W f , bridging the crack and, not part of the computational domain. The chemical balances of the oxidation processes are:

B 4 C 1.6 (s) + 4.6 O 2 (g) à 2 B 2 O 3 (l) + 1.6 CO 2 (g) (1)
C(s) + O 2 (g) à CO 2 (g) (2)
Note that the boron carbide in reaction (1) contains excess carbon with respect to the usual stoichiometry B 4 C, as found in CVD experiments 19,[START_REF] Rebillat | Procs 5 th Intl. Conf on High Temperature Ceramic Matrix Composites (HT-CMC5)[END_REF] -we will note it B x C in the following. The mass balance equations are considered as averaged over the direction perpendicular to the crack. This implies that several heights involved in the problem become variables defined over some part of the resolution domain. They are:

-the height h m of consumed boron carbide, defined only on W r ; -the height h l of liquid boron oxide, which is larger than h m because of the molar volume differences (i.e. B 2 O 3 has a larger molar volume than B x C), defined on W r ; -the height of consumed pyrocarbon h p , existing only on the fiber/matrix interface ∂W fi (considering that its thickness is so small that it does not need to be meshed). The balance equations for oxygen, boron carbide and boron oxide involve source terms describing diffusion-limited oxidation reactions. Pyrocarbon oxidation appears as a boundary condition on the fiber surface for the oxygen and pyrocarbon balances. The oxygen mass balance reads:

( ) ( ) 2 2 2 2 6 . 4 O l r O O g O g t C k C D h C h Φ - = ∇ - ⋅ ∇ + ∂ in W r and W i (3) 2 2 2 O p p fi O O g C e k n C D h - = ⋅ ∇ - on ∂W fi (4) ( ) ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = - - - = ⋅ ∇ - 0 2 2 2 2 2 2 2 ext O O ext O O g O ei O O C C or C C D n C D δ on ∂W ei (5)

INRIA

In eqs. (3-4) we have used the variable h g which is the available height for oxygen gas diffusion, and which is given by h g = e f + F r (h m -h l ), where F r indicates whether one lies or not in the reactive phase region. Eq. ( 5) displays two possible kinds of boundary conditions: either of the von Neumann type (similar to Fourier boundary condition in heat balances), or of the Dirichlet type (fixed assigned concentration). The time variations of the various heights involved in the system are given by the following ODEs:

( ) ( )

2 3 2 6 . 4 2 O l O B C B r g t C k h x υ υ - Φ = ∂ in W r (6) ( ) 2 3 2 6 . 4 2 O l O B r l t C k h υ Φ = ∂ in W r (7) ( ) 2 O p PyC p t C k h υ = ∂ on ∂W fi ( 8
)
The ! h term represents the possibility of a time variation in the crack opening, due for instance to variations in mechanical loading. The source terms for pyrocarbon and boron carbide oxidation depend on a reaction rate that follows a linear-parabolic law, as in the Deal-Grove model 20 : l /h l ), eq. ( 10) ensures a proper initialization of the numerical resolution. The initial conditions are:

1 2 1 - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = l O l BxC l D h k k (9)
( ) 1 ; 0 ; ; 0 2 2 2 = = = = = = g p l f O ext O O h h h e h C t C (11)
The remaining notation is explained in Table 1, where typical values for all coefficients are also reported.

Numerical resolution

The resolution domain is obtained from image processing routines. Two options are possible: (i) a routine allows extracting from a micrograph of a transverse slice of a given tow a finite element mesh, with a proper description of the domains as listed above; (ii) another routine allows extracting from a micrograph information on the fiber size and neighborhood distribution, then synthesizing a random artificial medium with the same fiber properties, and finally creating the multilayer matrix by mesh dilation operations. In this way, it is possible to work either on an actual material, or on a virtual material. The oxygen diffusion problem (3-5) has been discretized using a standard Galerkin approach on triangular P1 Lagrange finite elements. Time integration of the ODEs (6-8) is performed using an implicit backward differencing and solving the algebraic problem by means of standard Newton-Raphson iterations. Using dimensional analysis, an estimate of the relative importance of the different phenomena has been carried out, based on the data reported in Table 1. As a result of this analysis, one can the time scale associated to gaseous oxygen diffusion is so small that on the time scale of liquid oxide formation the oxygen diffusion equation (3-5) might be considered as a decoupled steady problem to be used as a starting point to integrate the height evolution equations (6-8). This RR n° 7845 decoupled approach allows to first solve for the oxygen concentration field with frozen the heights; then, larger time steps are taken for the computation of an increase of the heights, assuming a steady concentration field. After one time step, thenew values for the heights are used to obtain a solution of the oxygen diffusion equation. The procedure is carried out until h g reaches zero on some point of the mesh.

Results

Illustrative results are presented on an SH-CMC tow image containing approximately 50 fibers. The fibers were obtained by scanning and thresholding a transverse micrograph of a CMC; the matrix layers were generated by mesh dilations. The layer sequence is (from fiber to surface) PyC/SiC/B x C/SiC/B x C/SiC. Symmetry boundary conditions were applied everywhere the material was meeting the image border. All pores were considered as equally connected to the exterior oxygen source; a Neumann boundary condition was applied. Figure 2 gives shows the numerical results when plug formation is starting, at t = 112 s. Two-dimensional effects, such as competition in oxygen consumption between adjacent fibers, can be dramatic, as shown by the histograms of liquid height formation on Figure 3. As we can see, the spreading of the distribution is rather large at the beginning, then further spreads, before arriving at the maximal achievable height. When the surface barrier formation begins, i.e. when a first mesh point reaches the maximal achievable value for the liquid height (here, 1.35 µm), other surface points are still far from this state (minimum around 0.7 µm). Consequently, it takes a rather long time between the beginning and the end of barrier formation. The estimation in the present case is ~ 80 s, i.e. close to the beginning time itself (112 s); however it is certainly an overestimation, that can be corrected when taking into consideration the glass phase flow.

INRIA

Conclusion and outlook

This work focuses on the modeling of the first moments of the formation of a glassy boron oxide protective barrier in an SH-CMC tow. The originality is to carry out modeling in 2D in a crack plane supposed to be transverse with respect to fibers. The results show that it is possible to quantify the gradients of barrier formation throughout the section, and obtain an average and standard deviation of the self-healing behavior. Future work is aimed to account for the flow of the glass phase and of the diffusion of oxygen through it, in order to describe the "latency" period of the self-healing composite.

Figure 1 .

 1 Figure 1. a) Structure of a multilayer self-healing ceramic-matrix composite viewed in a transverse direction with respect to fibers b) Crack in the multilayer matrix with a bridging fiber 18 ; c) simulation domain viewed in elevation, d) detail of the simulation domain subsets and of several parameters.

  not very necessary on a physical point of view (because k BxC is very large in front of D O2

Figure 2 .

 2 Figure 2. Visualization of the numerical results in a base case (see Table 1), at the beginning of healing (t = 112 s). a) Computational mesh, b) Scaled oxygen partial pressure field, c) liquid oxide height, d) Consumed pyrocarbon height.

Figure 2b clearly evidences

  Figure2bclearly evidences the depletion of oxygen throughout the crack. Indeed, consumption by the glass-forming carbide and by carbon is not balanced by supply of gas by diffusion. So, only the boron carbide layers located close enough to the boundaries have reacted, as seen on fig.2c. In correlation with this, the consumed pyrocarbon heights are larger for the mesh points lying closer to the borders; a depletion effect is also seen for the points which lie close to a second fiber, as seen in fig.2d. In fig.3aare reported the consumed PyC heights as a function of the distance to the closest border, at t = 32s. Evidently, the mesh points lying closest to the surface are more oxidized; however one notes a large dispersion of the consumed height values everywhere in the material, as commented before on fig. 2d. To quantify this dispersion, the mesh points lying closest to the border (highlighted in fig.3a) are studied and their histogram is plotted at figure 3b. The standard deviation is 14% of the average value, a rather appreciable amount.

Figure 3 .

 3 Figure 3. Consumed PyC height results at t = 32 s. a) Values at mesh points vs. distance to the closest border. b) Values histogram for the mesh points lying closest to the border.

Figure 4 .

 4 Figure 4. Liquid height histograms at selected time values, for the mesh points lying closest to the border.

Table 1 . Meaning and units of the symbols used in the model.

 1

	Symbol	Meaning	Value	Unit
	C O2	Oxygen concentration		mol.m -3
	D O2	Oxygen diffusion coefficient		m 2 .s -1
	D O2	g	Oxygen diffusion coefficient	9.4•10 -5	m 2 .s -1
	D O2	l	Oxygen diffusion coefficient	5.4•10 -10	m 2 .s -1
	e f			Crack opening	10 -6	m
	e p			Pyrocarbon interphase thickness	10 -7	m
	h g			Available height for gaseous transfer		m
	h l			Liquid phase height		m
	h m			Consumption height of B x C		m
	h p			Consumption height of pyrocarbon		m
	k PyC	Rate constant for PyC consumption	1.7•10 -3	m.s -1
	k p			Effective rate constant for PyC consumption		m.s -1
	k BxC	Rate constant for B x C consumption	1•10 -3	m.s -1
	k l			Effective rate constant for B x C consumption		m.s -1
	u B2O3	Boron oxide molar volume	4.492 10 -5 m 3 .mol -1
	u BxC	Boron carbide molar volume	2.314 10 -5 m 3 .mol -1
	u PyC	Pyrocarbon molar volume	0.7059 10 -5 m 3 .mol -1
	d			Boundary layer thickness	2 10 -4	m
	F r			Reactive phase indicator (1 in W r , 0 elsewhere)		-

  1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 3.1 3.2

								a)				b)
											Number of mesh
											points (--)
												Consumed PyC height values (µm)
		700 700
		600 600
	Number of mesh points (-)	200 300 400 500 200 Number of mesh points (-) 300 400 500								t=016 t=036 t=056 t=096 t=016 t=036 t=056 t = 32 s! t = 72 s! t = 112 s! t=096 t = 192 s!
		100 100
		0	0
			0	0	0.2 0.2	0.4 0.4	0.6 0.6	0.8 0.8	1	1	1.2 1.2	1.4 1.4	1.6 1.6
							Liquid height (!m) Liquid height (!m)

Maximal achievable liquid height!

INRIA Bordeaux Sud-ouest, EPI BACCHUS -vdrean@inria.fr, mario.ricchiuto@inria.fr

LCTS, Lab. for ThermoStructural Composites, University Bordeaux Iperrot@lcts.u-bordeaux1.fr, couegnat@lcts.u-bordeaux1.fr, vinhola@lcts.u-bordeaux1.fr

INRIA Bordeaux Sud-ouest, EPI BACCHUS -vdrean@inria.fr, mario.ricchiuto@inria.fr

LCTS, Lab. for ThermoStructural Composites, University Bordeaux Iperrot@lcts.u-bordeaux1.fr, couegnat@lcts.u-bordeaux1.fr, vinhola@lcts.u-bordeaux1.fr

RR n° 7845

Acknowledgements

This work has been funded by GIS "Advanced Materials in Aquitaine" through a post-doc grant to V. D.. The authors also wish to thank Pr. F. Rebillat (LCTS, University Bordeaux) for fruitful discussion on the model.

INRIA


J.-C Walter

C Chatelain
email: chatelai@ijl.nancy-universite.fr

Relaxation at finite temperature in Fully-Frustrated Ising Models

Keywords: numbers: 05.10.Ln, 05.50.+q, 05.70.Jk, 05.70.Ln

We consider by means of Monte Carlo simulations the relaxation in the paramagnetic phase of the anti-ferromagnetic Ising model on a triangular lattice and of a fully-frustrated Ising model on a square lattice. In contradistinction to previous studies of the second model, we show that spin-spin correlation functions do not decay with a stretched-exponential law at low temperature but that both models display an exponential decay with logarithmic corrections that are interpreted as the signature of topological defects.

Introduction

The study of relaxation in frustrated systems is of particular interest because their dynamics may become anomalously slow at low temperature and eventually freeze below a certain temperature. The paradigmatic example of such a slow dynamics is given by spin glasses which combine both frustration and randomness. Frustration alone does not necessarily imply a glassy dynamics. Even though each one of its plaquettes is frustrated, the anti-ferromagnetic Ising model on a triangular lattice [START_REF] Wannier | Antiferromagnetism. the triangular ising net[END_REF] (AFIM) was shown to display the same dynamics as unfrustrated systems [START_REF] Zeng | Zero-temperature phase transitions of an antiferromagnetic ising model of general spin on a triangular lattice[END_REF]. Interestingly, if elastic deformations of the lattice are allowed and coupled to the spin degrees of freedom, the relaxation is well described by a (stretched-exponential) Kohlrausch-Williams-Watts law like in glasses [START_REF] Chen | Elastic antiferromagnets on a triangular lattice[END_REF][START_REF] Gu | Monte carlo study of a compressible ising antiferromagnet on a triangular lattice[END_REF][START_REF] Yin | Slow dynamics and aging in a nonrandomly frustrated spin system[END_REF]. This behavior is observed experimentally for example with closely-packed colloidal spheres [START_REF] Han | Geometric frustration in buckled colloidal monolayers[END_REF].

In this paper, we consider both the AFIM and the Fully-Frustrated Ising Model (FFIM) defined on a square lattice by the Hamiltonian -βH =

x,y [σ x,y σ x+1,y + (-1) f (x,y) σ x,y σ x,y+1 ], σ x,y = ±1

where f (x, y) = x + y in the so-called zigzag model and f (x, y) = x in pileup-domino configuration. In both cases, each plaquette of the square lattice contains an odd number of anti-ferromagnetic bonds and, as a consequence, is frustrated. The pileup-domino configuration allows for exact diagonalization of the transfer matrix by free fermion techniques [START_REF] Villain | Spin glass with non-random interactions[END_REF]. Both AFIM and FFIM are known to belong to the same universality class. At the critical temperature T c = 0, spin-spin autocorrelation functions decay algebraically

C(t, s) = σ(t)σ(s) ∼ (t -s) -η/z (2)
where η = 1/2 and z = 2 is the dynamical exponent. In the case of the AFIM, evidences have been given of the existence of topological defects interacting through a logarithmic Coulombian potential in the paramagnetic phase [START_REF] Moore | Vortex dynamics and entropic forces in antiferromagnets and antiferromagnetic potts models[END_REF][START_REF] Yin | Effective field theory of the zero-temperature triangularlattice antiferromagnet: A monte carlo study[END_REF]. While no signature of these defects was observed in earlier simulations of the AFIM [START_REF] Kim | Nonequilibrium critical dynamics of the triangular antiferromagnetic ising model[END_REF], we have shown that such defects manifest themselves in the FFIM [START_REF] Walter | Logarithmic corrections in the ageing of the fully frustrated ising model[END_REF]. When the system is initially prepared in the paramagnetic phase and then quenched at T c = 0, they pin the domain walls and slow down their motion so that aging takes place in the same way as in homogeneous systems but with logarithmic corrections. The correlation length grows as ξ(t) ∼ (t/ ln t) 1/z and spin-spin autocorrelation functions behave as

C(t, s) ∼ s -η/z t ln s s ln t -λ/z (3)
where both z and λ are compatible with 2.

This result suggests that at finite temperature, the relaxation should be exponential, like in unfrustrated systems, but with logarithmic corrections due to the presence of topological defects. However, on the basis of Monte Carlo simulations, a stretchedexponential relaxation ‡

C(t) ∼ e -(t/τ ) κ (4)
with a temperature-dependent exponent κ has been reported in the FFIM below the temperature T p ≃ 1.701 at which Kasteleyn-Fortuin clusters start to percolate [START_REF] Fierro | Percolation transition and the onset of nonexponential relaxation in fully frustrated models[END_REF][START_REF] Franzese | Precursor phenomena in frustrated systems[END_REF]. In this work, we reconsider the AFIM and the FFIM and show that Monte Carlo data are in better agreement with an exponential decay when taking into account logarithmic corrections than with a stretched exponential. In the first section, a dynamical scaling hypothesis is set up to predict the behavior of equilibrium spin-spin autocorrelation functions in fully-frustrated Ising models. The expression is then compared with Monte Carlo data for the AFIM in the second section and the FFIM in the third one.

Dynamical scaling of spin-spin autocorrelation functions

We analyze the two-time spin-spin correlation functions in the framework of dynamical scaling [START_REF] Hohenberg | Theory of dynamic critical phenomena[END_REF]. Upon a dilatation with a scale factor b, the equilibrium correlation C( r, t, T ) = σ 0 (0)σ r (t) at temperature T is assumed to satisfy the homogeneity relation

C( r, 1/t, |T -T c |) = b -2xσ C(r/b, b z /t, |T -T c |b 1/ν ) (5)
where x σ is the scaling dimension of magnetization density with 2x σ = η for twodimensional systems and z is the dynamical exponent. The motivation for the last two arguments of the scaling function in equation ( 5) comes for the behavior of the correlation length either with time, ξ ∼ t 1/z , or with temperature, ξ ∼ |T -T c | -1/ν . Letting b = t 1/z in equation ( 5), we obtain

C( r, t) = t -η/z C(r/t 1/z , |T -T c |t 1/νz ) (6)
The algebraic prefactor corresponds to the critical behavior while the scaling function includes all corrections to it. The characteristic time

τ ∼ ξ z ∼ |T -T c | -νz (7)
appears as the relaxation time of the system. In the following, we are interested only in autocorrelation functions, i.e. r = 0. Moreover, we expect an exponential decay of the scaling function C(t/τ ) in the paramagnetic phase. Therefore, the autocorrelation function generally reads at equilibrium

C(t, T ) ∼ e -t/τ t η/z (8)
Two modifications need to be made to apply this hypothesis to the AFIM and the FFIM. First, as already mentioned, the existence of topological defects slows down the ‡ The exponent in the stretched-exponential law is usually denoted β. We will use the notation κ instead in order to avoid confusion with the inverse temperature β = 1/k B T . motion of domain walls and thus the growth of the correlation length, i.e. ξ ∼ (t/ ln t) 1/z . A logarithmic correction has to be included in the scaling hypothesis ( 6)

C(0, t) = t -η/z C(0, |T -T c |(t/ ln t) 1/νz ) ∼ e -t/τ ln t t η/z (9)
The algebraic decay in front of the scaling function is not affected by logarithmic corrections because it describes the critical behavior for which topological defects are paired. In contradistinction, the scaling function C corresponds to the deviation to this behavior caused by non-vanishing scaling fields. The latter brings the system out-of-criticality, i.e. into the paramagnetic phase where free topological defects are encountered. C should thus involve logarithmic corrections. Note that the same kind of behavior was assumed during aging (3). The second modification concerns the relaxation time. Since the correlation length is known [START_REF] Forgacs | Ground-state correlations and universality in two-dimensional fully frustrated systems[END_REF][START_REF] Lukic | Finite-size scaling in villain's fully frustrated model and singular effects of plaquette disorder[END_REF] to diverge exponentially with temperature, i.e. ξ ∼ e 2/T and not algebraically, equation ( 7) has to be replaced by

τ ∼ ξ z ∼ e 2z/T . ( 10
)

Relaxation of the AFIM

We have studied the AFIM and FFIM by means of large-scale Monte Carlo simulations for a two-dimensional lattice with 192 × 192 sites. The dynamics is the heat-bath local Markovian process introduced by Glauber [START_REF] Glauber | Time-Dependent statistics of the ising model[END_REF]. We shall consider first the AFIM. Inverse temperatures β = 1/k B T in the range [0.75; 3.00] have been considered. Data have been averaged over 30, 000 independent histories. Error bars on the average correlation C(t, s) have been estimated at each times t and s as the standard deviation of the data produced during the different independent histories.

To check that the system had thermalized when the measurements were started, we monitored the two-time spin-sin correlation functions C(t, s) which are expected to depend only on t -s in the stationary state. On figure 1, spin-spin autocorrelation functions of the AFIM are plotted at different temperatures with respect to t -s. For sufficiently large values of s, one can observe the collapse of the curves, indicating that the system has reached equilibrium. In the example of the data presented in figure 1, one is led to the conclusion that below the inverse temperature β = 1/k B T = 1.75, equilibrium has already been reached at time s = 1000. At β = 1.75, equilibrium is reached only at time s = 2000 while for larger values of β, no collapse is observed indicating that equilibrium is not reached yet at time s = 4500. For this reason, only inverse temperatures β ≤ 1.75 with s ≥ 2000 will be considered in the following.

We have tested three possible scenarii: exponential decay of the spin-spin autocorrelation functions [START_REF] Hohenberg | Theory of dynamic critical phenomena[END_REF], exponential decay with logarithmic corrections [START_REF] Kim | Nonequilibrium critical dynamics of the triangular antiferromagnetic ising model[END_REF], and stretched exponential (4). On figure 2, the scaling function C(t, s)(t -s) versus t -s with a waiting time s = 2000 (ensuring equilibration as discussed above).

In the scenario [START_REF] Hohenberg | Theory of dynamic critical phenomena[END_REF], this function is expected to decay exponentially. This behavior is indeed observed over a large range of times t -s for all inverse temperatures β ≤ 1.75. We estimated the relaxation time by interpolation of the scaling function as e -t/τ with a sliding interpolation window. As can be seen in the inset of figure 2, the dependence of τ on β is well described by an exponential, as expected for the AFIM (see equation 10). A fit gives τ ∼ e 4.31(2)β , a behavior which is close to the expected one [START_REF] Lukic | Finite-size scaling in villain's fully frustrated model and singular effects of plaquette disorder[END_REF], thought significantly outside error bars.

On figure 3, the second scenario ( 9) is tested. The scaling function C(t, s)(t -s) 1/4 is plotted versus (t-s)/ ln(t-s). The interpolation with an exponential gives the values represented in the inset. As before, the relaxation time behaves exponentially with the inverse temperature. A fit gives the law τ ∼ e 3.92(2)β , a behavior which is much closer to [START_REF] Lukic | Finite-size scaling in villain's fully frustrated model and singular effects of plaquette disorder[END_REF] than without logarithmic corrections, even though still outside error bars.

Finally, the stretched exponential scenario is tested on figure 4. The expected power-law behavior in the long-time regime (according to equation 4) is observed only for a narrow interval of times, much narrower than in the two previous scenarii. We have nevertheless interpolated the data with equation (4) to extract the relaxation time. The result is plotted in the inset of figure 4. An exponential growth is observed but with a quite different factor: τ ∼ e 2.79 (2) is plotted versus the inverse temperature β. The straight line is the interpolated behavior τ ∼ e 2.79(2)β .

Relaxation of the FFIM

The FFIM has been studied in the zig-zag bond configuration for a 192 × 192 square lattice. We restricted ourselves to inverse temperatures smaller or equal to β = 1. As a consequence, smaller times and fewer histories were necessary. Monte Carlo data have been averaged over 1, 000 independent histories. On figure 5, spin-spin autocorrelation functions C(t, s) are presented for different temperatures. The collapse of the different curves indicates that thermalization is achieved very rapidly. Like for the AFIM, we have compared the three scenarii: exponential decay of the spin-spin autocorrelation functions (8), exponential decay with logarithmic corrections (9), and stretched exponential (4).

First, we present the test of the stretched-exponential scenario for the FFIM (figure 6). The logarithm -ln C(t, s) displays a behavior similar to the AFIM. At first sight, the power-law behavior that is expected according to equation (4) seems to be observed over a larger range of times t-s than for the AFIM but one should keep in mind that the inverse temperatures are different on figures 4 and 6. We have estimated the exponent κ of the stretched exponential. Our data confirm the general trend observed in references [START_REF] Fierro | Percolation transition and the onset of nonexponential relaxation in fully frustrated models[END_REF][START_REF] Franzese | Precursor phenomena in frustrated systems[END_REF]. The exponent κ indeed strongly depends on the temperature. However, we do not recover a purely exponential decay above the percolation temperature of the Fortuin-Kasteleyn clusters. We do not observe two regimes (κ = 1 for T > T p and κ decreasing below T p ) as previously found but instead, an exponent slowly approaching a value κ = 1 (see figure 7). The discrepancy may be explained by the larger lattice size used in this study (L = 192 instead of L = 64) and by the difficulty to identify a sufficiently large power-law regime, especially at low temperature. As shown in the inset of figure 6, the relaxation time τ does grow exponentially over the whole range of inverse temperatures β considered, in contradistinction to the theoretical prediction [START_REF] Lukic | Finite-size scaling in villain's fully frustrated model and singular effects of plaquette disorder[END_REF]. However, our data are compatible with two distinct regimes of exponential growth with a prefactor estimated to be 3.60 [START_REF] Glauber | Time-Dependent statistics of the ising model[END_REF] for β ≥ 0.6 and 2.50(2) for β < 0.6 (note that β t ≃ 0.59).

On figure 8, the second scenario is tested. As expected, the scaling function C(t, s)(t -s) 1/4 displays an exponential decay with t -s over a large range of times t -s for different temperatures. However, the relaxation time does not grow exponentially over the whole range of temperatures considered (inset of figure 8). Again, two regimes can be distinguished: the prefactor in the exponential ( 10) is estimated to be 4.78(2) for β ≥ 0.6 and 3.72(5) for β < 0.6. These values should be compared with the theoretical prediction 4.

On figure 9, the third scenario involving logarithmic corrections is tested. The scaling function C(t, s)(t -s) 1/4 decay exponentially with (t -s)/ ln(t -s) over a large range of times t -s for different temperatures. The relaxation time τ now displays an exponential growth over the whole range of temperatures considered and our estimate 4.097(26) of the prefactor in the argument of the exponential is much closer to the theoretical prediction 4 according to equation [START_REF] Lukic | Finite-size scaling in villain's fully frustrated model and singular effects of plaquette disorder[END_REF]

Conclusions

We have analyzed the decay of the equilibrium two-time spin-spin correlation functions in the paramagnetic phase against three scenarii: exponential decay of the spin-spin autocorrelation functions (8), exponential decay with logarithmic corrections (9), and stretched exponential (4). To distinguish between them, we tested the thermal behavior [START_REF] Lukic | Finite-size scaling in villain's fully frustrated model and singular effects of plaquette disorder[END_REF] of the relaxation time in the different scenarii. Our results are summarized in table 1. The best agreement with the theoretical prediction 4 is obtained with an exponential decay with logarithmic corrections (9) for both the AFIM and the FFIM. Error bars are unfortunately smaller than the deviation from the theoretical prediction. A great care has been taken in the computation of error bars (from the data production to the analysis) so we can only invoke finite-size effects or other systematic deviations. The deviation for the two other scenarii being much larger, we do not expect them to get closer to 4 than the exponential decay with logarithmic corrections. These results bring further evidences of the existence of topological defects in the paramagnetic phase of the AFIM and the FFIM.

Figure 5 .

 5 Figure 5. Relaxation of the spin-spin correlation function C(t, s) of the FFIM versus t -s at six different inverse temperatures β = 0.50, 0.60, 0.70, 0.80, 0.90 and 1.00. The different colors correspond to different times s as indicated in the legend.

Figure 6 .

 6 Figure 6. Logarithm -ln C(t, s) of the spin-spin correlation of the FFIM at different inverse temperatures β = 0.50, 0.60, 0.70, 0.80, 0.90, and 1.00 (from top to bottom). Only the time s = 160 is presented. In the inset, the relaxation time τ obtained by interpolation of the correlation function as e -t-s τ

Figure 7 .

 7 Figure 7. Exponent κ of the stretched exponential versus the inverse temperature β for both the AFIM (black) and the FFIM (red).

Figure 9 .

 9 Figure 9. Scaling function C(t, s)(t-s) 1/4 at time s = 1600 versus (t-s)/ ln(t-s) for the AFIM at different inverse temperatures β = 0.50, 0.60, 0.70, 0.80, 0.90, and 1.00 (from bottom to top). In the inset, the relaxation time τ obtained by interpolation of the scaling function as e -(t-s)/τ ln(t-s) is plotted versus the inverse temperature β. The straight line is the interpolated behavior τ ∼ e 4.097(26)β .

  1/4 is plotted

		1 1 1 1 1 1
								s=1000
								s=1500
								s=2000
								s=2500
		0.1 0.1 0.1 0.1 0.1 0.1						s=3000 s=3500
								s=4000
								s=4500
	C(t,s) C(t,s) C(t,s) C(t,s) C(t,s) C(t,s)	0.01 0.01 0.01 0.01 0.01 0.01						β=2.00
								β=1.75
		0.001 0.001 0.001 0.001 0.001 0.001		β=1.25	β=1.50		β=1.65
		0 0 0 0 0 0	β=0.75	2000 2000 2000 2000 2000 2000	4000 4000 4000 4000 4000 4000	t-s t-s t-s t-s t-s t-s	6000 6000 6000 6000 6000 6000	8000 8000 8000 8000 8000 8000

Figure 1. Relaxation of the spin-spin correlation function C(t, s) of the AFIM versus t -s at six different inverse temperatures β = 0.75, 1.25, 1.50, 1.65, 1.75, and 2.00. The different curves correspond to different times s. Error bars are represented but they are hardly visible on the figure because they are of order O(10 -5 ), i.e. much smaller than the symbols.

  β .

	9 9 9 9 9
				1000
			1 1 1 1 1
				100
	C(t,s)(t-s) 1/4 C(t,s)(t-s) 1/4 C(t,s)(t-s) 1/4 C(t,s)(t-s) 1/4 C(t,s)(t-s) 1/4 3 -ln C(t,s) 3 -ln C(t,s) 3 -ln C(t,s) 3 -ln C(t,s) 3 -ln C(t,s)	0,1 0,1 0,1 0,1 0,1	0,6 10	0,8 0,6 100	0,8	1	1	1,2 β	1,2 β	1,4 1,4	1,6 1,6	1,8 1,8
	1 1 1 1 1			10
		0,01 0,01 0,01 0,01 0,01
	C(t,s) (t-s) 1/4 C(t,s) (t-s) 1/4 C(t,s) (t-s) 1/4 C(t,s) (t-s) 1/4 C(t,s) (t-s) 1/4	0 0 0 0 0 0 0,01 0,1 1 0 0,01 0,1 1 0 0,01 0,1 1 0 0,01 0,1 1 0 0,01 0,1 1 Figure 3. Scaling function C(t, s)(t -s) 10 2000 4000 t-s 2000 4000 t-s 2000 4000 t-s 2000 4000 t-s 2000 4000 t-s 200 400 (t-s)/ln(t-s) 200 400 (t-s)/ln(t-s) 200 400 (t-s)/ln(t-s) 200 400 (t-s)/ln(t-s) 200 400 (t-s)/ln(t-s) 0,6 0,8 1 10 100 1000 τ 100 1000 600 600 600 600 600 1 t-s 10 100 1000 t-s 10 100 1000 t-s 10 100 1000 t-s 10 100 1000 t-s	6000 6000 6000 6000 6000 1,2 β		8000 8000 8000 8000 8000 1,6 1,8 10000 800 800 800 800 800 1,4 10000 10000 10000 10000

τ

Figure 2. Scaling function C(t, s)(t -s) 1/4 at time s = 2000 versus t -s for the AFIM at different inverse temperatures β = 0.75, 1.25, 1.50, 1.65, and 1.75 (from left to right). In the inset, the relaxation time τ obtained by interpolation of the scaling function as e -(t-s)/τ is plotted versus the inverse temperature β. The straight line is the interpolated behavior τ ∼ e 4.31(2)β . 1/4 at time s = 2000 versus (t -s)/ ln(t -s) for the AFIM at different inverse temperatures β = 0.75, 1.25, 1.50, 1.65, and 1.75 (from from left to right). In the inset, the relaxation time τ obtained by interpolation of the scaling function as e -(t-s)/τ ln(t-s) is plotted versus the inverse temperature β. The straight line is the interpolated behavior τ ∼ e 3.92(2)β . τ Figure 4. Scaling function -ln C(t, s) at time s = 2000 versus t -s for the AFIM at different inverse temperatures β = 0.75, 1.25, 1.50, 1.65 and 1.75 (from top to bottom). In the inset, the relaxation time τ obtained by interpolation of the correlation function as e -t-s τ κ

Table 1 .

 1 than without logarithmic corrections. Prefactor in the argument of the exponential growth[START_REF] Lukic | Finite-size scaling in villain's fully frustrated model and singular effects of plaquette disorder[END_REF] of the relaxation time for the two models AFIM and FFIM in the different studied scenarii. The theoretical prediction is 4.

				100
				10
		0,1
	1/4 C(t,s) (t-s)	0,01		1	0,2	0,4	0,6 β	0,8	1
		0,001
		0	100	200 t-s			300		400
	Model Stretched exponential Exponential		Log. corrections
	AFIM 2.79(2)		4.31(2)		3.92(2)
	FFIM 3.60(5) (β ≥ 0.6)	4.78(2) (β ≥ 0.6) 4.097(26)

τ Figure

8

. Scaling function C(t, s)(t -s) 1/4 at time s = 160 versus t -s for the FFIM at different inverse temperatures β = 0.50, 0.60, 0.70, 0.80, 0.90, and 1.00 (from left to right). In the inset, the relaxation time τ obtained by interpolation of the scaling function as e -(t-s)/τ is plotted versus the inverse temperature β. The two straight lines are the interpolated behaviors τ ∼ e 4.78(2)β for β ≥ 0.6 and τ ∼ e 3.72(5)β for β < 0.6.