Magalie Thomassin
email: magalie.thomassin@ims-bordeaux.fr

Rachid Malti
email: rachid.malti@ims-bordeaux.fr

Subspace method for continuous-time fractional system identification

Keywords: Fractional state-space representation, continuous-time identification, subspace method, multivariable model

The aim of this paper is to develop a subspace method for state-space identification of continuous-time systems using fractional commensurate models. As compared to the classical state-space representation, the commensurate differentiation order must be estimated besides the state-space matrices. The latter are estimated with conventional subspace-based techniques using QR and singular value decompositions, whereas the commensurate order is estimated using nonlinear programming. This is the first method developed for multi-input multi-output system identification of fractional models. The performances are demonstrated by simulations at various signal-to-noise ratios assuming a known then an unknown commensurate order.

INTRODUCTION

Fractional models have witnessed a growing interest during the last years. Many diffusive phenomena can be modeled by fractional transfer functions. In electrochemistry for instance, diffusion of charges in acid batteries is governed by Randles models [START_REF] Sabatier | Fractional system identification for lead acid battery sate charge estimation[END_REF]] that involve Warburg impedance with an integrator of order 0.5. Electrochemical diffusion showed to have a tight relation with derivatives of order 0.5 [START_REF] Oldham | Diffusive transport to planar, cylindrical and spherical electrodes[END_REF]]. In thermal diffusion of a semi-infinite homogeneous medium, [START_REF] Battaglia | Solving an inverse heat conduction problem using a noninteger identified model[END_REF] have shown that the exact solution of the heat equation links thermal flux to a half order derivative of the surface temperature on which the flux is applied. Time-domain system identification using fractional models was initiated in the late nineties. [START_REF] Oustaloup | Identification of non integer order system in the time-domain[END_REF] developed a method based on the discretization of the fractional differential equation using Grünwald definition and on the estimation of its coefficients using least squares. [START_REF] Trigeassou | Modeling and identification of a non integer order system[END_REF] based their identification method on the approximation of a fractional integrator by a rational model. Then, they deduced the fractional model after estimating its rational approximation. [START_REF] Cois | Fractional state variabe filter for system identification by fractional model[END_REF] proposed several extensions of equation error methods, such as the state variable filters and the instrumental variable (IV), to fractional system identification. [START_REF] Aoun | Synthesis of fractional Laguerre basis for system approximation[END_REF] synthesized fractional orthogonal bases generalizing various bases (Laguerre, Kautz,...) to fractional differentiation orders for identification issues. Recently, [START_REF] Malti | An optimal instrumental variable method for continuoustime fractional model identification[END_REF] have extended the concept of optimal IV methods to fractional systems. For an overview of these identification methods refer to Malti et al. (2008a).

In this paper, we consider the problem of identifying a continuous-time (CT) fractional system in its state-space form. Only few papers deal with system identification us-ing fractional state-space representation [START_REF] Cois | Fractional state variabe filter for system identification by fractional model[END_REF]; [START_REF] Poinot | Identification of fractional systems using an output-error technique[END_REF]]. They are based on the minimization of an output error criterion by nonlinear programming. These methods are well suited for single-input single-output (SISO) systems, and are generally difficult to apply in the multi-input multi-output (MIMO) case because the number of parameters to estimate becomes large. Here, a subspace method to estimate the matrices of the CT fractional state-space representation is proposed. It is an extension of the methods presented in the literature for rational systems [START_REF] Haverkamp | Identification of continuous-time MIMO state space models from sampled data, in the presence of process and measurement noise[END_REF]; [START_REF] Johansson | Stochastic theory of continuous-time state-space identification[END_REF][START_REF] Johansson | Stochastic theory of continuous-time state-space identification[END_REF]] to the fractional ones. Other subspace techniques for identifying CT systems with rational models can be found in [START_REF] Bastogne | A PMF-based subspace method for continuous-time model identification. application to a multivariable winding process[END_REF]; [START_REF] Ohsumi | Subspace identification for continuous-time stochastic systems via distribution-based approach[END_REF]; [START_REF] Mercère | Subspace based methods for continuous-time model identification of MIMO systems from filtered sampled data[END_REF]. So, the proposed method inherits the advantages of subspace methods which stem from the reliability of the numerical algorithms using the QR and the singular value decompositions [START_REF] Katayama | Subspace methods for system identification[END_REF]]. Thus, it does not involve nonlinear optimization to obtain state-space matrices. In addition, no canonical form (such as modal or companion) of the state-space representation is required. Finally, the proposed subspace algorithms can be applied to the identification of both SISO and MIMO fractional systems. As will be seen later, the state-space representation of a fractional commensurate system involves an additional parameter which is the commensurate order. This parameter is the only one computed by minimizing an output error criterion with a nonlinear optimization technique.

In section 2, some recalls about fractional systems are given. Section 3 presents the method proposed to estimate the matrices of the CT fractional state-space representation, followed by simulation examples in Section 4. Finally, Section 5 is devoted to the estimation of the fractional commensurate order. Monte Carlo simulations are made to check the statistical properties of the estimator.

FRACTIONAL SYSTEMS

The concept of differentiation to an arbitrary non-integer order α, with α ∈ R + * (set of strictly positive real numbers), D α ∆ = d dt α , was defined in the 19 th century by Riemann and Liouville. The α-order fractional derivative of f (t) is defined as an integer derivative of order ⌊α⌋ + 1 (⌊.⌋ stands for the floor operator) of a non-integer integral of order α -⌊α⌋ [START_REF] Samko | Fractional Integrals and Derivatives : Theory and Applications[END_REF]]:

D α f (t) = d ⌊α⌋+1 dt ⌊α⌋+1 1 Γ(⌊α⌋ + 1 -α) t 0 f (τ ) (t -τ ) α-⌊α⌋ dτ
where the Euler gamma function Γ(β) is defined for every β ∈ R + * by: Γ

(β) = ∞ 0 z β-1 e -z dz.
The Laplace transform, denoted by L , is a more concise algebraic tool which allows to write, in case of zero initial conditions, [START_REF] Oldham | The Fractional Calculus: Integrations and Differentiations of Arbitrary Order[END_REF]]:

L D α f (t) = s α F (s), with F (s) L f (t)
where s is the Laplace variable.

A SISO fractional system is governed by a fractional differential equation:

y (t) + a 1 D α1 y (t) + • • • + a mA D αm A y (t) = b 0 D β0 u (t) + b 1 D β1 u (t) + • • • + b mB D βm B u (t)
where (a j , b i ) ∈ R 2 , and the differentiation orders α 1 < α 2 < . . . < α mA , β 0 < β 1 < . . . < β mB are allowed to be non-integer positive numbers. State space representation was extended by [START_REF] Matignon | Some results on controllability and observability of finite-dimensional fractional differential systems[END_REF] to commensurate fractional systems, where all the differentiation orders are multiple integers of α. The extension was done by allowing the differentiation order of the state vector to be any commensurate order α ∈ R + * . The fractional state space representation is presented in a MIMO case as:

D α x(t) = Ax(t) + Bu(t), (1) y
(t) = Cx(t) + Du(t)
(2) where x ∈ R n is the state vector, u ∈ R m the input vector, y ∈ R p the output vector, A ∈ R n×n , B ∈ R n×m , C ∈ R p×n , D ∈ R p×m are constant matrices. Zero initial conditions are considered: x(t) = 0 for t ≤ 0. Matignon (1998) proved that the fractional system (1)-( 2) is stable if and only if:

0 < α < 2 and | arg(λ k )| > α π 2 ∀k = 1, . . . , n where λ k is the k th -eigenvalue of A and | arg(λ k )| ≤ π.
The conversion of ( 1)-( 2) to the MIMO transfer function form is obtained as for the rational systems by:

G(s) = C(s α I -A) -1 B + D.
In the following, assume that (A, B) is reachable and (C, A) is observable. The controllability and the observability conditions of a state space representation of a commensurate fractional system are the same as for rational systems [START_REF] Matignon | Some results on controllability and observability of finite-dimensional fractional differential systems[END_REF]].

One of the main difficulties with fractional models is the time-domain simulation. This problem has been extensively studied and an overview of the principal methods can be found in [START_REF] Aoun | Numerical simulations of fractional systems: an overview of existing methods and improvements[END_REF]. The most commonly used approximation of fractional operators is the recursive distribution of zeros and poles, proposed by [START_REF] Oustaloup | La dérivation non entière. Théorie, synthèse et applications[END_REF], which approximates the frequency behavior of s α in the frequency range [ω A , ω B ]. Nevertheless, this approximation has null asymptotic behaviors at low and high frequencies, which can introduce a static error between the fractional model and its approximation. To avoid this drawback, [START_REF] Trigeassou | Modeling and identification of a non integer order system[END_REF] suggested to use the conventional integrator outside the frequency range [ω A , ω B ]:

s -α [ωA,ωB ] = G α s Nc k=1 1 + s/ω ′ k 1 + s/ω k (3)
where:

• N c is the number of cells (directly related to the quality of the approximation), • G α is fixed so that s -α has the same gain as s -α

[ωA,ωB] in the middle of the interval [ω A , ω B ],
• ω ′ k and ω k are respectively zeros and poles recursively distributed in the frequency range [ω b , ω h ] = [σ -1 ω A , σω B ] where σ is generally set to 10 to minimize border effects. They are defined by the following relations:

ω ′ k = γω k , ω k+1 = ηω ′ k , α = 1 - log γ log γη .
This approximation is used to simulate the fractional systems presented in this paper with the parameters: N c = 20, ω A = 10 -5 and ω B = 10 5 .

SUBSPACE ALGORITHMS FOR FRACTIONAL STATE-SPACE IDENTIFICATION

Consider the linear CT fractional state-space representation (1)-( 2). The problem in this section is to estimate the system matrices A, B, C, D from sampled input-output data 1 {u k } N -1 k=0 and {y k } N -1 k=0 . The commensurate order α is assumed to be known. The case where α is unknown will be discussed in section 5.

The Laplace transform of (1)-(2) gives:

s α X(s) = AX(s) + BU (s) (4) Y (s) = CX(s) + DU (s).
(5) The problem in these relations is the use of the s α -operator which amplifies noise in the time domain, especially at high frequencies. To avoid this problem, the low-pass filter:

λ = 1 1 + s ω f α = 1 1 + τ s α with τ = (1/ω f ) α (6) is introduced. Then, (4) is transformed into: X(s) = A λ [λX(s)] + B λ [λU (s)] (7) with A λ = I +τ A and B λ = τ B.
Application of the inverse Laplace transform leads to the following system of linear equations:

x(t) = A λ [λx(t)] + B λ [λu(t)] (8) y(t) = Cx(t) + Du(t)
(9) where λx(t) and λu(t) correspond to the states and the inputs prefiltered by λ in (6). Then, from (9), it is found by recursion that:

y(t) = Cx(t) + Du(t) = CA λ [λx(t)] + CB λ [λu(t)] + Du(t) . . . = CA k λ [λ k x(t)] + k j=1 CA k-j λ B λ [λ k-j+1 u(t)] + Du(t)
for k ∈ N * , where λ k x(t) denotes the signals obtained from x(t) by filtering through a series of k low-pass filters λ. In the same way, it is found for l ∈ N * :

λ l y(t) = Cλ l x(t) + Dλ l u(t) = CA λ [λ l+1 x(t)] + CB λ [λ l+1 u(t)] + Du(t) . . . = CA k-l λ [λ k x(t)] + k-l j=1 CA k-j-l λ B λ [λ k-j+1 u(t)]
+ Du(t) with k ≥ l. As a consequence, the input-output data can be formulated as the following extended linear model:

Y(t) = O i X (t) + Ψ i U(t) (10
) with input-output and state variables:

Y(t) =       λ i-1 y(t) λ i-2 y(t)
. . .

λ 1 y(t) y(t)       , U(t) =       λ i-1 u(t) λ i-2 u(t) . . . λ 1 u(t) u(t)       , X (t) = λ i-1 x(t)
and:

O i =     C CA λ . . . CA i-1 λ     , Ψ i =       D 0 • • • 0 CB λ D . . . . . . . . . . . . . . . 0 CA i-2 λ B λ • • • CB λ D      
where O i ∈ R ip×n is the extended observability matrix and Ψ i ∈ R ip×im is a block Toeplitz matrix. Now, from N input-ouput available samples observed at discrete times t k = kT s for k = 0, . . . , N -1, the extended linear model (10) can be rewritten as:

Y N = O i X N + Ψ i U N (11)
where

U N =       [λ i-1 u] 0 [λ i-1 u] 1 • • • [λ i-1 u] N -1 [λ i-2 u] 0 [λ i-2 u] 1 • • • [λ i-2 u] N -1 . . . . . . . . . [λu] 0 [λu] 1 • • • [λu] N -1 u 0 u 1 • • • u N -1       ∈ R mi×N
by using the following notation for sampled filtered data:

[λ j u] k = λ j u(t k ). ( 12
)
The matrices Y N ∈ R pi×N and X N ∈ R n×N are constructed in a similar way. Equation ( 11) enables to use subspace identification algorithms as in their original nonfactional discrete-time version. We have chosen to use the ordinary MOESP (MIMO Output-Error State Space) algorithm [START_REF] Katayama | Subspace methods for system identification[END_REF]]. The principle of this algorithm is as follows:

(1) Compute an LQ decomposition of the data matrix:

U N Y N = L 11 0 L 21 L 22 Q T 1 Q T 2 ( 13
)
where

L 11 ∈ R im×im , L 21 ∈ R ip×im , L 22 ∈ R ip×ip
with L 11 , L 22 lower triangular, and

Q 1 ∈ R N ×im , Q 2 ∈ R N ×ip are orthogonal.
(2) Compute a singular value decomposition (SVD) of the L 22 matrix approximating the column space of O i :

L 22 = [U 1 U 2 ] Σ 1 0 0 0 V T 1 V T 2 (14)
where U 1 ∈ R ip×n and U 2 ∈ R ip× (ip-n) . The state order n can be estimated from the SVD since n = dim Σ 1 in the noiseless case. (3) Estimate the extended observability matrix:

Ôi = U 1 Σ 1/2 1 . (15)
(4) Estimate the C matrix: Ĉ = Ôi (1 : p, 1 : n).

(5) Estimate A λ by solving the linear equation: Ôi 1 :

p(i -1), 1 : n A λ = Ôi p + 1 : ip, 1 : n . ( 16
)
(6) Estimate the B λ and D matrices. For that purpose, it can be shown that:

U T 2 Ψ i = U T 2 L 21 L -1 11
(17) which is a linear equation with respect to B λ and D. Define:

U T 2 [L 1 L 2 • • • L i ] (18) U T 2 L 21 L -1 11 [M 1 M 2 • • • M i ] ( 19
)
with L k ∈ R (ip-n)×p and M k ∈ R (ip-n)×m for k = 1, . . . , i. Thus, from (17):

L 1 D + L 2 ĈD + • • • + L i Ĉ Âi-2 λ B λ = M 1 L 2 D + L 3 ĈD + • • • + L i Ĉ Âi-3 λ B λ = M 2 . . . L i-1 D + L i ĈD = M i-1 L i D = M i . Define Lk = [L k . . . L i ] ∈ R (ip-n)×(i+1-k)p , k = 2, .
. . , i, and get the following overdetermined system of linear equations:

       L 1 L2 Ôi-1 L 2 L3 Ôi-2 . . . . . . L i-1 Li Ô1 L i 0        D B λ =       M 1 M 2 . . . M i-1 M i       (20)
where the block coefficient matrix in the left-hand side is of dimensions i(ipn) × (p + n). Estimates of B λ and D are estimated by finding the least-squares solution of (20).

The matrices of the fractional CT state-space representation (4)-( 5) are then deduced as follows: Â = 1 τ ( Âλ -I), B = 1 τ Bλ . The Ĉ and D matrices do not change.

EXAMPLE

The algorithm is applied to input-output data of length N = 1023 generated by simulating the linear system (1)-( 2) with α = 0.9,

A = 0 -0.1 1 -0.2 , B = 1 0 , C = 0 0.1 0.5 -0.1 , D = 0 0 ,
and zero initial conditions. The input signal is a pseudorandom binary sequence (PRBS) with maximum length. The sampling period is T s = 0.05s. The outputs (Fig. 1) are corrupted by white noise with a signal-to-noise ratio (SNR) of 20 dB. The algorithm parameters are set to i = 8 and ω f = 6 and the commensurate order α is set to its true value. In these conditions, a second-order state-space representation is estimated with a good accuracy with the following companion form: Â = 0 -0.0984 1 -0.1853 , B = 1 0 , Ĉ = 0.0002 0.1015 0.5114 -0.0962 , D = 0.0004 0.0024 .

The estimated outputs ŷ are plotted in Fig. 1. The normalized prediction error norm ( ŷy 2 / y 2 ) equals 0.1109 (-19.1 dB) for the first output and 0.1082 (-19.3 dB) for the second one. The singular values of L 22 are plotted in Fig. 2. The state-space order was easily deduced since we observe an abrupt change at n = 2. Different order estimation criteria may also be used [START_REF] Bauer | Order estimation for subspace methods[END_REF]].

Fig. 3 shows the influence of the method parameters (the number of block rows i and the filter frequency ω f ) on the normalized error norm ŷy 2 / y 2 for the first output (similar results are obtained with the second output).

The estimation is considered incorrect if the normalized error norm is greater than or equal to one (i.e. if its logarithm is greater than or equal to zero). So, in Fig. 3, all normalized errors greater than one are set to one. These figures also show the sensitivity to stochastic disturbance by considering various SNRs: 20, 15 and 10 dB. The contour levels indicate that ω f can be chosen in a suitable range over one decade, between 1 and 10 for a SNR greater than 15 dB. In the case of a higher noise level, this range is slightly reduced and the choice of i acts more on the results: fine results are obtained from i = 6 with w f ∈ [3, 10] for SNR=10 dB.

To analyze the estimator statistical properties, 500 data sets, each with different realization of the noise, are generated for three SNRs (20, 15 and 10 dB). The means of the normalized error norms, obtained with i = 8 and ω f = 6, are given in table 1 and Fig. 4 shows the estimated poles.

The normalized mean squared error (MSE) of the poles are indicated in table 2. It can be seen that the estimator is biased. Indeed, the persistent excitation of the input, obtained with a PRBS for example, may not be enough to guarantee the consistency of the estimates for MOESP [START_REF] Bauer | Analysis of the asymptotic properties of the MOESP type of subspace algorithms[END_REF]].

Table 1. Means of the normalized error norms over 500 runs (α known). Table 2. Normalized MSE of poles (α known).

SNR=20 dB -0.0676.10 -4 ± i 0.3257.10 -4 SNR=15 dB -0.0684.10 -3 ± i 0.2475.10 -3 SNR=10 dB -0.0834.10 -2 ± i 0.2665.10 -2

FRACTIONAL ORDER ESTIMATION

In this section, we assume that the fractional derivative order α ∈]0, 2[ is unknown and has to be estimated by minimizing a quadratic criterion:

α = arg min α∈]0,2[ 1 2 ŷc (α) -y c 2 2 , ( 21
)
where y c is the vector (of length pN ) resulting from the concatenation of the p system outputs and ŷc (α) is the vector resulting from the concatenation of the p outputs estimated with the MOESP method presented in section 3 for a given α. The value of α in ( 21) is obtained through a nonlinear optimization technique. The iterative algorithm used is based on a subspace trust-region method and on the interior-reflective Newton method. As a consequence, the proposed MOESP algorithm is executed at each iteration of the optimization technique. The number of state variables n is assumed to be known. The normalized prediction error norm is equal to 0.1050 (-19.57 dB) for the first output and 0.1008 (-19.93 dB) for the second one. As expected, the normalized errors obtained here are better than those obtained in the previous section for α fixed to its true value because α is optimized so as to minimize the prediction error. In Fig. 5, the criterion to be minimized to estimate the commensurate order α is plotted. We observe that the minimum can be easily reached by the optimization algorithm. Fig. 5. Criterion to minimize with respect to α for a data set with SNR=20 dB (true commensurate order: 0.9) Now, in the context of commensurate order optimization, we present results obtained using a Monte Carlo simulation with 500 runs for three values of the SNR (20, 15 and 10 dB). The histograms of the estimated order α are given in Fig. 6. We observe a bias which increases with the noise level. This is obviously due to the bias introduced by the MOESP algorithm. Table 3 gives the means of the normalized error norms. Again, they are slightly better than the ones obtained in table 1 when α is fixed to its true value. Nevertheless, concerning the estimates of the poles, the previous results are better than those obtained here (see Fig. 7). Indeed, the normalized MSEs of the poles presented in table 4 are larger than the ones in table 2. Note that similar results are obtained by considering systems with an other commensurate order.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Output data y and output of the estimated model ŷ for a PRBS as input (SNR=20 dB).

Fig. 4 .

 4 Fig. 4. Estimated poles over 500 runs for various SNRs,where + denotes the true poles (α known)

Fig. 3 .

 3 Fig. 3. Normalized error norm ŷy 2 / y 2 in dB versus the number of block rows i and the filter frequency ω f for various SNRs. This method is applied on the same data set used in the first example, with SNR=20 dB. The algorithm parameters are also the same. The estimation results are: Â = 0 -0.1008 1 -0.1757 , B = 1 0 , α = 0.8894, Ĉ = 0.0012 0.1012 0.4997 -0.0923 , D = 1.3782 10 -3 -0.6763 10 -3 .

Fig. 7 .

 7 Fig. 7. Estimated poles over 500 runs for various SNRs, where + denotes the true poles (α estimated)

Table 3 .

 3 Means of the normalized error norms over 500 runs (α estimated).

		output SNR=20 dB SNR=15 dB SNR=10 dB
	mean of the normalized error norm	1 2	0.1026 (-19.78 dB) (-14.73 dB) (-9.56 dB) 0.1834 0.3327 0.1005 0.1797 0.3282 (-19.95 dB) (-14.91 dB) (-9.68 dB)

Table 4 .

 4 Normalized MSE of poles. This paper focuses on the identification of CT systems using fractional state-space models. Thanks to an adapted data filtering, we have shown that subspace-based algorithms can be used to estimate the state-space matrices. Simulation results have shown that the estimation results are not very sensitive to small variations of the tuning parameters. The commensurate differentiation order is estimated by using nonlinear programming. Simulation examples have shown that the estimators are biased. So, a Fig.6. Histograms of the estimated commensurate order α for various SNRs future work consists in developing stochastic methods to eliminate this bias.

	SNR=20 dB	-0.4952.10 -4 ± i 0.8875.10 -4
	SNR=15 dB	-0.2530.10 -3 ± i 0.6315.10 -3
	SNR=10 dB	-0.1208.10 -2 ± i 0.5740.10 -2
	6. CONCLUSION

The discrete-time variables are denoted by x k and correspond to the time sampling with a constant sampling period Ts of the CT variable x(t): x k = x(kTs).


Abdelhamed Mohamed

Stefan Nemanja

Member, IEEE Marco Di Perović

Fellow, IEEE Renzo

Intelligent Omni-Surfaces (IOSs) for the MIMO Broadcast Channel

In this paper, we consider intelligent omni-surfaces (IOSs), which are capable of simultaneously reflecting and refracting electromagnetic waves. We focus our attention on the multiple-input multiple-output (MIMO) broadcast channel, and we introduce an algorithm for jointly optimizing the covariance matrix at the base station, the matrix of reflection and transmission coefficients at the IOS, and the amount of power that is reflected and refracted from the IOS. The distinguishable feature of this work lies in taking into account that the reflection and transmission coefficients of an IOS are tightly coupled. Simulation results are illustrated to show the convergence of the proposed algorithm and the benefits of using surfaces with simultaneous reflection and refraction capabilities.

I. INTRODUCTION

Reconfigurable intelligent surfaces (RISs) have emerged as a promising approach to improve the wireless communication channel quality and to extend the network coverage [START_REF] Renzo | Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and road ahead[END_REF], [START_REF]Smart radio environments empowered by reconfigurable AI metasurfaces: An idea whose time has come[END_REF]. However, the vast majority of works consider surfaces that can only reflect the incident signals, which limits the coverage capabilities offered by RISs [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF], [START_REF] Huang | Reconfigurable intelligent surfaces for energy efficiency in wireless communication[END_REF]. To tackle this problem, the recently proposed concept of intelligent omnisurface (IOS) provides 360 • coverage thanks to surfaces that can simultaneously reflect and refract the incident signals [START_REF] Zhang | Intelligent omni-surfaces: Ubiquitous wireless transmission by reflective-refractive metasurfaces[END_REF]. The reflection and refraction capabilities of the incident electromagnetic waves are controlled through the optimization of two interlinked sets of reflection and transmission coefficients. In general, in other words, it is not possible to control the reflection and transmission coefficients independently.

Recently, a few research works have analyzed the performance of IOSs. In [START_REF] Perera | Sum rate maximization in STAR-RIS assisted fullduplex communication systems[END_REF], the authors investigate the weighted sum-rate maximization under quality of service requirements and unit modulus constraints for the IOS elements, by utilizing the successive convex approximation method. In [START_REF] Yang | Energy-efficient wireless communications with distributed reconfigurable intelligent surfaces[END_REF], the energy efficiency maximization problem is studied, and an optimization algorithm is proposed for jointly optimizing the transmit power and the passive beamforming at the IOS. In [START_REF] Niu | Simultaneous transmission and reflection reconfigurable intelligent surface assisted MIMO systems[END_REF], the weighted sum-rate of an IOS-aided multiple-input multiple-output (MIMO) system is maximized by using the alternating optimization method. The precoding matrices are obtained by the Lagrange dual method, while the reflection and A. Mohamed, N. S. Perović, M. Di Renzo are with Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des Signaux et Systèmes, 3 Rue Joliot-Curie, 91192 Gif-sur-Yvette, France. Email: marco.di-renzo@universite-parissaclay.fr. This work was supported in part by the European Commission through the H2020 PAINLESS project under grant agreement number 812991, the H2020 SURFER project under grant agreement number 101030536, the H2020 ARIADNE project under grant agreement number 871464, and the H2020 RISE-6G project under grant agreement number 101017011. transmission coefficients are obtained by the penalty concaveconvex method. Interested readers can consult [START_REF] Xu | STAR-RISs: Simultaneous transmitting and reflecting reconfigurable intelligent surfaces[END_REF], [START_REF] Liu | STAR: Simultaneous transmission and reflection for 360°c overage by intelligent surfaces[END_REF], [START_REF] Mu | Simultaneously transmitting and reflecting (STAR) RIS aided wireless communications[END_REF], [START_REF] Xu | Simultaneously transmitting and reflecting (STAR) intelligent omni-surfaces, their modeling and implementation[END_REF] for further information on IOSs.

In contrast to the available works, we focus our attention on optimizing the IOS in the MIMO broadcast channel. Also, for the first time, we explicitly take into account the dependence between the reflection and transmission coefficients, based on a recently implemented IOS prototype [START_REF] Zhang | Intelligent omni-surfaces for full-dimensional wireless communications: Principles, technology, and implementation[END_REF]. Specifically, the main contributions of this paper are as follows:

• We utilize the duality between the broadcast channel and the multiple access channel to maximize the achievable sum-rate. We formulate a joint optimization problem for the users' covariance matrices, the reflection and transmission coefficients, and the power ratio between the reflected and transmitted power. We analyze the case studies with continuous-valued and discrete-valued phase shifts for the reflection and transmission coefficients, and we assume that they are not independent of each other. • Due to the non-convexity of the formulated optimization problem, and the coupling between the optimization variables in the objective function, we propose an alternating optimization algorithm to solve the aforementioned problem. The optimal users' covariance matrices are obtained by applying the dual decomposition and the block coordinate maximization (BCM) method, while the optimal phase shifts of the reflection and transmission coefficients of the IOS elements are formulated in a simple expression. In addition, the power ratio between the reflected and refracted power is computed iteratively by utilizing the gradient ascent method. • Simulation results show that the proposed algorithm converges relatively fast (i.e., within a few iterations) to a local optima. Moreover, we quantify the impact of discretizing the reflection and transmission coefficients for a two-state IOS testbed platform.

The paper is organized as follows. In Section II, we present the system model for the IOS-aided MIMO broadcast channel. In Section III, we formulate the optimization problem to maximize the achievable sum-rate. In Section IV, we describe the optimization algorithm to solve the optimization problem. In Section V, we provide simulation results that illustrate the achievable sum-rate. Conclusions are drawn in Section VI.

Notation: Bold upper and lower case letters denote matrices and vectors, respectively. C m×n denotes the space of m × n complex matrices. H T , H H , |H| and Tr(H) represent the transpose, Hermitian transpose, the determinant and the trace

II. SYSTEM MODEL

We consider an IOS-assisted multi-user MIMO broadcast channel, where a multi-antenna base station (BS) with N t antennas communicate with K users, and each user is equipped with N r antennas. The total number of users in the reflection and transmission sides of the IOS are denoted by KR and KT , respectively, with K = KR + KT . We assume that the BS and the users' antennas form a uniform linear array (ULA), and the inter-antenna separations are s t and s r , respectively. In addition, we assume that the IOS comprises N elements that can simultaneously reflect and refract the incident signals.

Since the users can be located either in the reflection or the transmission sides of the IOS, a single user receives either the signal reflected or the signal refracted from the IOS.

To serve the users located in the reflection and transmission sides simultaneously, each element of the IOS applies a complex-valued reflection and transmission coefficient. The matrix of N reflection or transmission coefficients is denoted by

F(θ) = diag(θ) ∈ C N ×N , where θ = {θ r , θ t }, and θ r = √ ρ βr 1 , βr 2 , . . . , βr N T ∈ C N ×1 and θ t = √ 1 -ρ βt 1 , βt 2 , . . . , βt N T ∈ C N ×1
. Specifically, the parameter ρ ∈ [0, 1] is an optimization variable that accounts for the amount of power that the IOS directs towards the users in the reflection side (1 -ρ is the amount of power directed towards the transmission side). Also, βr n and βt n are the reflection and transmission coefficients of the IOS. These reflection and transmission coefficients are interlinked and are optimized in pairs for each IOS element, i.e., βr n , βt n ∈ ( βr s1 , βt s1 ), ( βr s2 , βt s2 ), . . . , ( βr sQ , βt sQ ) = Ψ, where Q is the number of (reflection, transmission) coefficient pairs.

Due to the presence of the IOS, the end-to-end channel for the k r ∈ [1, . . . KR]-user located in the reflection side is:

H kr = D kr + G kr F(θ r )U = D kr + G kr √ ρ × diag βr 1 , βr 2 , ..., βr N U = D kr + √ ρ N n=1 βr n g n,kr u n (1)
and the end-to-end channel for the k t ∈ [1, . . . KT ]-user located in the transmission side is:

H kt = D kt + G kt F(θ t )U = D kt + G kt 1 -ρ × diag βt 1 , βt 2 , ..., βt N U = D kt + 1 -ρ N n=1 βt n g n,kt u n
(2) In a compact form, we have:

H k (θ) = D k + G k F(θ)U (3)
where D k ∈ C Nr×Nt denotes the channel matrix between the BS and the k-th user, U ∈ C N ×Nt denotes the channel matrix between the BS and the IOS, and G k ∈ C Nr×N denotes the channel matrix between the IOS and the k-th user.

To simplify the notation, we write H k instead of H k (θ), where the dependence on θ is implicit. Thus, the received signal at the k-th user is written as:

y k = H k x k + K j=1,j =k H k x j + n k (4)
where H k ∈ C Nr×Nt is the channel matrix for the k-th user, x k ∈ C Nt×1 is the transmitted signal intended for the k-th user, and x j ∈ C Nt×1 for j = k are the transmitted signals intended for other users, which act as interference for the detection of x k . The noise vector n k ∈ C Nr×1 consists of independent and identically distributed (i.i.d.) elements whose distribution is CN (0, N 0 ), where N 0 is the noise variance.

III. PROBLEM FORMULATION

We are interested in maximizing the achievable sum-rate of the considered IOS-aided wireless communication system. To accomplish this, we exploit the fact that the achievable rate region of a Gaussian MIMO broadcast channel can be achieved by dirty paper coding (DPC) [START_REF] Weingarten | The capacity region of the Gaussian multiple-input multiple-output broadcast channel[END_REF]. DPC enables us to reduce the interference in a communication system, i.e., to perfectly eliminate the interference term j<k H k x j for the k-th user. In this regard, the ordering of the users clearly matters. Let π be an ordering of users, i.e., a permutation of the set {1, 2, . . . , K}. Then, for this ordering, the achievable rate for the k-th user can be computed as [15, Eq. ( 3)]:

R π(k) = log 2 I + H π(k) j≥k S π(j) H H π(k) I + H π(k) j>k S π(j) H H π(k) (5)
where

S k = E x k x H k
0 is the input covariance matrix of user k. We assume a sum-power constraint at the BS, i.e.:

K k=1 Tr S k ≤ P ( 6
)
where P is the maximum total power at the BS. Therefore, the achievable sum-rate optimization problem for the RIS-assisted MIMO broadcast channel can be expressed as:

max S, βl ,ρ K k=1 log 2 I + H π(k) j≥k S π(j) H H π(k) I + H π(k) j>k S π(j) H H π(k) (7a) s.t. K k=1
Tr S k ≤ P ; S k 0, ∀k, (7b)

βr n , βt n ∈ Ψ; ∀n ∈ [1, 2, . . . , N ], (7c)
0 ≤ ρ ≤ 1 (7d)
where βl = { βr , βt } and l ∈ {r, t},. To solve this problem, we exploit the duality between the MIMO broadcast channel and the multiple access channel, as recently done in [START_REF] Perović | On the achievable sum-rate of the RIS-aided MIMO broadcast channel[END_REF]. Accordingly, we reformulate the optimization problem in [START_REF] Yang | Energy-efficient wireless communications with distributed reconfigurable intelligent surfaces[END_REF] as follows:

max S, βl ,ρ f ( S, β, ρ) = log 2 I + K k=1 H H k SH k s.t. K k=1
Tr Sk ≤ P ; Sk 0, ∀k,

βr n , βt n ∈ Ψ; ∀n ∈ [1, 2, . . . , N ], 0 ≤ ρ ≤ 1, (8)
where H H k represents the dual multiple access channel corresponding to H k and Sk is the dual multiple access input covariance matrix of the k-th user. Once the input covariance matrices ( Sk ) K k=1 in the dual multiple access channel are found, the corresponding covariance matrices (S k ) K k=1 in the broadcast channel are obtained from [15, Eq. ( 11)], as:

S k = B -1/2 k F k G H k A 1/2 k Sk A 1/2 k G k F H k B -1/2 k ( 9
)
where

A k = I + H k ( k-1 i=1 S i )H H k , B k = I + K i=k+1 H H i Si H i , and F k Λ k G H k is the singular value de- composition of B -1/2 k H H k A -1/2 k .

IV. PROPOSED OPTIMIZATION METHOD

To solve the formulated problem, we propose an alternating optimization method. The users' covariance matrices are optimized by exploiting the dual decomposition method in [START_REF] Perović | On the achievable sum-rate of the RIS-aided MIMO broadcast channel[END_REF]. The reflection and transmission coefficients of the IOS are obtained by generalizing the method in [START_REF] Zhang | Capacity characterization for intelligent reflecting surface aided MIMO communication[END_REF], which is applicable to reflecting surfaces, under the assumption of continuous-valued coefficients. The corresponding discretevalued reflection and transmission coefficients are obtained by projecting the obtained solutions onto the feasible set of possible discrete values. Moreover, we present a gradientbased method for optimizing the power ratio.

A. Covariance Matrix Optimization

For given values of the power ratio and the reflection and transmission coefficients, the achievable sum-rate optimization problem in ( 8) is simplified as follows:

max S log 2 I + K k=1 H H k Sk H k (10a) s.t. S ∈ S (10b)
where S ( Sk ) K k=1 for S = { S | K k=1 Tr Sk ≤ P ; Sk 0 ∀k}. As described in [START_REF] Perović | On the achievable sum-rate of the RIS-aided MIMO broadcast channel[END_REF], [START_REF] Perović | On the maximum achievable sum-rate of the RISaided MIMO broadcast channel[END_REF], the optimization problem in ( 10) is solved by using the dual decomposition Lagrangian and the accelerated block coordinate maximization methods.

B. IOS Optimization

For given Sk K k=1 and ρ, the optimization problem in [START_REF] Niu | Simultaneous transmission and reflection reconfigurable intelligent surface assisted MIMO systems[END_REF] with respect to βl n , l ∈ {r, t}, can be written as follows:

max βl n f ( βl n ) = log 2 I + K k=1 H H k SH k s.t. βr n , βt n ∈ Ψ; ∀n ∈ [1, 2, . . . , N ] (11)
To tackle this problem, we first optimize both βr n and βt by assuming that they are continuous-valued and independent coefficients, and we then project the obtained solutions onto the set of feasible discrete phase shifts, by taking into account that the reflection and transmission coefficients are interlinked.

For ease of writing, we define the variables

G k = [g 1 , . . . , g n ], U = [u T 1 , ..., u T n ] T , and
Dk = D r i ∀i ∈ {1, ..., KR} D t j ∀j ∈ {1, ..., KT } (12)
ḡn,k = g r n,i ∀i ∈ {1, ..., KR} g t n,j ∀j ∈ {1, ..., KT }

Under the assumption that βr n and βt n can be optimized independently (this is removed next by projecting on the feasible set), the objective function can be rewritten as follows:

f ( βl n ) = log 2 I +     KR k=1 Dk + N n=1 √ ρ βr n ḡn,k un H ×S k Dk + N n=1 √ ρ βr n ḡn,k un     +     KT k=1 Dk + N n=1 √ 1 -ρ βt n ḡn,k un H ×S k Dk + N n=1 √ 1 -ρ βt n ḡn,k un     (14)
Thus, for optimizing βr n , we have:

f ( βr n ) = log 2 A r n + βr n B r n + βr * n B rH n ( 15
)
and, for optimizing βt n , we have:

f ( βt n ) = log 2 A t n + βt n B t n + βt * n B tH n ( 16
)
where:

B r n = √ ρ KR k=1 ( DH k + N l=1,n =l √ ρ βr * l u H l ḡH l,k )S k ḡn,k u n (17) B t n = KT k=1 ( 1 -ρ DH k + N n =l (1 -ρ) βt * l u H l ḡH l,k )S k ḡn,k u n ( 18
)
Algorithm 1: AO algorithm to solve [START_REF] Niu | Simultaneous transmission and reflection reconfigurable intelligent surface assisted MIMO systems[END_REF] Input: S(0) , βr(0) , βt(0) and ρ (0) with feasible values Apply the phase shift projection procedure, if required, based on (22)

1 Set i ← 0 2 repeat 3 Compute ( S (i+1)
7 end 8 ρ (i+1) = PD ρ (i) + υ (i) ∇ρf (ρ)| ρ=ρ (i) 9 i ← i + 1; 10 until a stopping criterion is met Output: S • = S (i) , ρ • = ρ (i) and βr•(i) , βt•(i) A r n = I + KT k=1 ( Dk + N n=1 √ 1 -ρ βt n ḡn,k un) H ×S k ( Dk + N n=1 √ 1 -ρ βt n ḡn,k un) + KR k=1 ρu H n ḡH n,k S k ḡn,k un + KR k=1 ( DH k + N l=1,n =l √ ρ βr * l u H l ḡH l,k ) ×S k ( Dk + N l=1,n =l √ ρ βr l ḡl,k u l ) (19)
A t n = I + KR k=1 ( Dk + N n=1 √ ρ βr n ḡn,k un) H ×S k ( Dk + N n=1 √ ρ βr n ḡn,k un) + KT k=1 (1 -ρ)u H n ḡH n,k S k ḡn,k un + KT k=1 ( DH k + N l=1,n =l √ 1 -ρ βt * l u H l ḡH l,k ) ×S k ( Dk + N l=1,n =l √ 1 -ρ βt l ḡl,k u l )
(20) The optimal solution of the optimization problem in ( 14) is then given by [START_REF] Perović | On the achievable sum-rate of the RIS-aided MIMO broadcast channel[END_REF]:

β * r n = exp(-j arg(σ r n )) (21a) β * t n = exp(-j arg(σ t n )) (21b)
where σ r n and σ t n are the only non-zero eigenvalues of (A r n ) -1 B r n and (A t n ) -1 B t n , respectively. The obtained continuous-valued solution of the reflection and transmission coefficients is denoted by β * r n , β * t n = SC. The corresponding discrete-valued solution in the considered feasible set is obtained by projecting SC on the discrete set Ψ. As a case study, we consider the IOS prototype in [START_REF] Zhang | Intelligent omni-surfaces for full-dimensional wireless communications: Principles, technology, and implementation[END_REF], where each IOS element can be configured in two states, where each state is identified by a pair of reflection and transmission coefficients, i.e., ( βr S1 , βt S1 ), ( βr S2 , βt S2 ) = {S1, S2} = Ψ. In this case, the projection can be formulated as follows:

βr• , βt• = (β r S1 , β t S1 ) , SC -S1 2 < SC -S2 2 (β r S2 , β t S2 ) , otherwise (22)
where βr• , βt• denotes the reflection and transmission coefficients that minimize the distance between the solution of the continuous-valued case and the feasible discrete set.

C. Power Ratio Optimization

For given values of S * k K k=1 and βr• , βt• , the optimization problem in (8) can be explicitly rewritten as:

max ρ f (ρ) = log 2 I + K k=1 H H k Sk H k (23a) s.t. 0 ≤ ρ ≤ 1 (23b)
The objective function f (ρ) can be rewritten as:

f (ρ) = log 2 I + K k=1 H H k S k H k (24) = log 2 I + 2 √ ρ KR k=1 N n=1 βr * n u H n ḡH n,k S k Dk +2 √ 1 -ρ KT k=1 N n=1 βt * n u H n ḡH n,k S k Dk + KR k=1   N n=1 N n=1 ρ βr * n βr n u H n ḡH n,k S k ḡn,k un + DH k S k Dk   + KT k=1   N n=1 N n=1 (1 -ρ) βt * n βt n u H n ḡH n,k S k ḡn,k un + DH k S k Dk  
(25) In order to find the optimal value of the power ratio ρ, we solve the equation f (ρ) = 0, where f (ρ) is the first-order derivative of f (ρ). For ease of notation, let us define:

X 1 = KR k=1 N n=1 βr * n u H n ḡH n,k S k Dk (26)
The first-order and second-order derivatives, f (ρ) and f (ρ), respectively, can be formulated as:

f (ρ) = ∇ ρ f (ρ) = Tr X-1 X1 √ ρ -X2 √ 1-ρ + (Y 1 -Y 2 ) (31) f (ρ) = -1 ln 2 IN V ( X) 2 × 0.5 -3/2 √ ρX 1 + -3/2 √ 1 -ρX 2 X + X1 √ ρ -X2 √ 1-ρ + (Y 1 -Y 2 ) 2 (32)
where ∇ ρ f (ρ) is the gradient of f (ρ) with respect to ρ. The objective function in (23a) is a concave function with respect to ρ, since the second derivative in (32) is no greater than zero. Therefore, the optimization problem in (23a) has a single optimum value. However, it is hard to find a closed-form solution by solving f (ρ) = 0. Thus, we utilize the projected gradient (PG) method. Specifically, ρ in the (i+1)-th iteration is updated as follows: where υ i in (33) is the step size that is updated by using the backtracking line search method in [START_REF] Perović | On the achievable sum-rate of the RIS-aided MIMO broadcast channel[END_REF], and the projection operator P D is defined as:

ρ i+1 = P D ρ i + υ i ∇ ρ f (ρ)| ρ=ρ i (
ρ = P D (ρ) = ρ min ρ • < ρ min ρ • ρ min ρ • ρ max ρ max ρ > ρ max (34)
where ρ • and ρ are the solutions obtained before and after implementing the projection, and ρ min and ρ max are the minimum and maximum values for ρ, respectively.

V. SIMULATION RESULTS

In this section, we evaluate the achievable sum-rate of the proposed algorithms with the aid of Monte Carlo simulations. Specifically, we compare the sum-rates under the assumption that the reflection and transmission coefficients are continuousvalued and independent values, and under the assumption that they are interlinked and belong to a discrete set.

The simulation setup is the following: the carrier frequency is f = 2 GHz (the wavelength is λ = 15 cm), s t = s r = s ris = λ/2 = 7.5 cm, the network topology is given in Fig. 1, N t = 8, P = 1 W, and N 0 = -110 dB. The IOS consists of N ris = 15 × 15 = 225 elements placed in a 15 × 15 square formation. The users are equipped with N r = 2 antennas and are randomly distributed within the disks shown in Fig. 1. The results are averaged over 100 independent channel realizations.

The achievable sum-rate is reported in Fig. 2 and Fig. 3. The proposed algorithm converges in a few iterations. By closing analyzing the sum-rate of the users located in the reflection and transmission sides of the IOS, we observe that the sum-rate loss in Fig. 3 as compared to Fig. 2 is due to the non-unitary value of the reflection and transmission coefficients for the considered feasible set in [13, Table 1].

VI. CONCLUSION

In this paper, we have proposed optimization algorithms for maximizing the sum-rate in IOS-aided MIMO broadcast channels. We have optimized the covariance matrices at the transmitter, the reflection and transmission coefficients of the IOS, and the ratio between the reflected and refracted power. The simulations results demonstrate that the presented algorithm provides a rapid convergence rate in a few iterations.

Fig. 1 :

 1 Fig. 1: Aerial view of the considered communication system. of H, respectively. The rank of H is denoted by rank(H) and λ max (H) is the largest singular value of H, while (•) * denotes the complex conjugate. The identity matrix is denoted by I and the matrix inversion operation is denoted by IN V (X).

  -j arg(σn)) using (21) 6

  S k ḡn,k u n (28) S k ḡn,k u n (29)

Fig. 3 :

 3 Fig. 3: Discrete-valued coefficients.


Marc Chaperon

François Laudenbach

Bernard Teissier

Publication of the mathematical works of René Thom in the collection Documents mathématiques of the French Mathematical Society *

Marc Chaperon and François Laudenbach

Volume II has just appeared. We present the whole of the project, initiated by André Haefliger.

René Thom in brief

Thom was one of the first members of the mathematical school created around Henri Cartan at the end of the second world war, and one of the most glorious, but from the start singular: provincial, born outside university circles, he had acquired very young an intimate knowledge of the differential calculus as its founders conceived it; unusual for his time, he did not distrust geometry, where he had developed his intuition to the point of "seeing" in four dimensions; finally, having followed Henri Cartan to Strasbourg as a young researcher in the CNRS, 1 he stayed there after the departure of his master, benefiting in particular from the influence of Charles Ehresmann. Thus blessed with a vision of the world complementary to that of the "Parisians," Thom resolved fundamental questions which others would doubtless not even have thought of. Some landmarks:

-From 1949 to 1956, Thom worked in algebraic topology and elaborated "a completely new way of studying differentiable manifolds" (Milnor), about which he obtained thus definitive results, giving birth to the theory of cobordism, for which he was awarded the Fields medal in 1958.

-From 1956 on he concentrated on the singularities of differentiable maps, which appear naturally in his vision of differential topology. Following on from Whitney, he then studied stratifications and introduced the "natural stratification of mapping spaces" made precise later by John Mather. 2 -From the middle of the 1960's, aware that singularities of maps and transversality aid in understanding all sorts of natural phenomena, he developed a catastrophe theory; this met with an extraordinary mediatic period after the appearance of his book Stabilité structurelle et morphogénèse in 1972.

-At the end of the 1970's, the sometimes delirious infatuation for catastrophes met with a sudden halt, without doubt hardly more justified than the ups and downs of fashions. Thom moved away then from mathematics in favour of philosophy and a fertile return to Aristotle. These changes of orientation must not mask a profound unity of thought: already on his entry to the École Normale Supérieure, Thom was tempted by the philosophy of science, from which the director of the École had to dissuade him. His great mathematical results have a strong philosophical connotation, far from the "modern" pretension to separate mathematics from the question of meaning.

Birth of the project

The admiration of André Haefliger for Thom goes back to their common Strasbourg years (1954-1958). 3 Their relations only ceased after the disappearance of Thom in October 2002.

Already possessing documents from a near half-century of exchanges, Haefliger consecrated much time from the end of 2010 to the study of the Thom archives, sorted, inventoried and classified from mid-April 2011 in the basements of IHES by its librarian Aurélie Brest, with, at the start, help from Herminia and André Haefliger. He discovered "veritable treasures," for example two quite amazing unpublished mathematical manuscripts, that we will speak of again. He associated Marc Chaperon to this discovery in September 2011, with the idea of putting in place with his help and that of Bernard Teissier a "classical" paper publication, annotated and commented on, of the mathematical works of Thom.

Teissier and Étienne Ghys would have liked this edition to be accompanied by the placing online of the complete works (mathematical or not) of Thom, published in the form of a CD-ROM at the beginning of 2003 by IHES. 4 Alas! This quite natural idea was halted by problems of copyright. 5 As to the CD-ROM, perfectly usable despite the obsolescence of its search engine, it is no longer for sale.

A first editorial committee, constituted in October 2011, met at the beginning of the following month; as well as Haefliger, Teissier and Chaperon, it included Alain Chenciner. The project was submitted at the end of November by Teissier to Pierre Colmez and immediately accepted in the collection Documents mathématiques newly created by the French Mathematical Society (SMF). François Laudenbach, Jean Petitot, David Trotman and, for volume I, Jean Lannes and Pierre Vogel quickly joined the editorial committee.

Overview

This publication does not pretend to be a substitute for the CD-ROM but to complete it:

-It is concentrated on the mathematical articles (or at least classified as such in Mathematical Reviews), which here are provided with mathematical or historical commentaries, justified by later developments and the continued relevance of this often visionary work. Certain commentaries are to be found following the article which they relate to, others serve as introduction to several articles. After each article one finds shorter notes concerning precise details.

-We have chosen to reproduce the originals instead of transcribing into modern mathematical typography, thus avoiding introducing errors and changing the page numbering.6 An exception comprises unpublished texts (absent from the CD-ROM) whose typed version was not legible enough.

As well as these unpublished texts, we publish a certain number of documents for their historical interest, for example large extracts from the correspondance with Cartan which led to Thom's thesis7 and fragments of letters written by him to his wife Suzanne during his stay in Princeton in 1951, they contain much information on the genesis of the work which followed.

-The bibliography, like the bibliographical notice, completes, corrects and enriches that of the CD-ROM, using earlier versions of Michèle Porte, Jean Petitot and Aurélie Brest, it covers all Thom's writings, mathematical or not.

Volume I

This volume, of 573 pages, appeared at last in April 2017 as n • 15 of the collection Documents mathématiques. It covers the publications of Thom from 1949 to 1959 except one, delayed to the beginning of the following volume. Given the dates, it is not surprising that it contains the major elements of Thom's work. We will evoke them a little more precisely.

The first article of Thom is a note in the Comptes Rendus de l'Académie des Sciences de Paris, 2 pages, entitled Sur une partition en cellules associée à une fonction sur une variété. It concerns a Morse function, that is to say a function at least of class C 2 whose critical points have Hessians of maximal rank. This note has had a very rich descendance that we have made explicit by a commentary in volume I. Thom is here curiously mute about what will become, starting in 1954, one of his paradigms: the notion of genericity. There lacks indeed the supplementary hypothesis of genericity ensuring that the said partition into cells be what is called a cellular decomposition, with its specific properties concerning the attachment of the cells. It is Stephen Smale, ten years later, who will make explicit this hypothesis on the gradient of the function in question, called today the Morse-Smale condition.

Next comes Thom's work in algebraic topology that we have preceded by a note of presentation written by J. Lannes and P. Vogel. This covers essentially the work of Thom for his thesis (1951) and his founding article on the theory of cobordism (1954). The thesis was published in the Annales de l'École Normale Supérieure (1952), 73 pages.

The article of 1954, 69 pages, has a title turned to the future, Quelques propriétés globales des variétés différentiables. It is here that appears for the first time the notion of transversality, a property which is generically satisfied, "most of the time," by a sufficiently differentiable map. This article is important in two aspects. On the one hand it opens the way to a new branch of differential topology which will become homotopic topology: this starts from the following problem that Thom solves, under what condition on a homology class of a manifold is it realisable by a submanifold. On the other hand, merely by the apparition of transversality-a few pages of this article-it opens the way to the study of singularities. With some exceptions, Thom is going to plunge into this second way. Let us see what happens.

Thom will realise that the way in which he proved his transversality theorem for a submanifold applies at once, without change, to certain situations that one now calls transversality under constraint. The first of these situations is the transversality to a submanifold in a jet space.

Jet bundles were discovered by Ehresmann. For a function f , let us say real for simplicity, k times differentiable and defined on a manifold M , its jet of order k at a point a of M is, in coordinates, the Taylor development of f at a. Of course, this polynomial depends on the coordinates but that two functions have the same Taylor development to order k does not depend on coordinates. One can then collect all the k-jets of functions at all points of M and we have thus created a new manifold J k (M, R) which fibres over M by indicating the point a where the k-jet is taken. A section of the bundle is rather rarely the collection of Taylor polynomials of a single function at all the points of M : there is said to be an integrability condition.

Despite this constraint, the theorem of transversality to a submanifold S of J k (M, R) holds, not in the space of all sections but among those which are integrable-and of course, one can replace R by any other manifold. However, the choice of S is decisive for the study of singularities of real functions. This is what Thom does in his article, Les singularités des applications différentiables, which appeared in the Annales de l 'Institut Fourier (1956). Before this paper, a "cap" written by Haefliger and based on his own archives, presents the work of Thom on singularities in the period 1956-1957. Moreover, with the agreement of the American Mathematical Society we have reproduced the long report of W. S. Massey in Math. Reviews; this provides a useful light on Thom's article.

Because we are in 1957, let us say a word about two unpublished papers. The first is entitled: Une démonstration d'un théorème de Lefschetz, and the second L'homologie des variétés de Stein. Thom writes at the very beginning of the second: "As a result we have the proof of a conjecture of J.-P. Serre: the homology groups H i (V, Z) of a Stein manifold of complex dimension n are zero for i > n." Thom presented these in a seminar in Chicago in February 1957. Why did they remain unpublished? We will never have the answer. In any case A. Andreotti and T. Frankel published an article, The Lefschetz theorem on hyperplane sections, Annals of Math. (1959), 14 pages. In their introduction one may read:

Recently Thom has given a proof (unpublished) which, as far as we know, is the first to use Morse's theory of critical points. We present in §3, in a slightly more general setting, an alternate proof inspired by Thom's discovery.

Let us finish this look at volume I with a last historical evocation. Thom gave a talk in a CNRS meeting in Lille in 1959. This talk was published in the Bulletin SMF with the title Remarques sur les problèmes comportant des inéquations différentielles globales. One cannot avoid relating this title with that of the book by M. Gromov, Partial Differential Relations, Springer-Verlag, 1986; moreover, Thom's article is cited there. Thom gives in this article, next to general considerations, a precise statement on the homology of an open set in a jet space. Here we are ten years before the thesis of Gromov (1969) who proves an analogous statement but, with an important difference, Gromov speaks of homotopy. As later, the statement of Gromov will be called an h-principle after the initial letter of homotopy, one could say that the statement of Thom is an h-principle (only homological) before its time.

The paradox of this affair is that Thom, at the time that S. Smale wrote his (brief) report for Math. Reviews, only believed in his theorem for the jet spaces of order one as shown by Smale's commentary. However, the commentary of D. Spring in our volume I, p. 562, makes precise that he himself knew how to give a proof using the holonomic approximation theorem of Eliashberg and Mishachev (L'Enseignement Math., 2001). Finally, there is a strong chance that the idea of Thom to introduce small sawtooths in his simplexes can be carried out to the end.

Volume II

As in volume I, the texts of Thom are preceded here by a complete bibliography of his works. 8 It begins with the course on singularities given at Bonn in 1959, as recorded by Harold Levine. Annotated and commented by Haefliger, this lecture course is followed by a translation of the preface and the table of contents of the Russian edition by V. I. Arnold, then by a letter to Haefliger of February 1959; despite a casual error, this shows that Thom had quickly noticed the modules (moduli) which complicate the theory.

The volume gathers then the articles of Thom published between 1962 and 1971,9 in chronological order of their appearance. The editorial committee was solicited more here than in volume I.

-Many of the long articles contained in the latter, written under the gaze of Henri Cartan and his school, required in fact no other commentary than the mention of their whys and wherefores. After his Fields medal, Thom is more alone but does not lose his tranquil audacity, which results in often prophetic (and, very rarely, badly written) mathematical papers requiring completion10 and, in this volume, commented on in more detail.

-This audacity led him to "leave the frame" with his catastrophe theory, of which the foundational articles figure in this volume II. Containing sometimes very new mathematics, they are part of another story, notably in biology, of which it was important to give an idea: Sara Franceschelli and Jean Petitot took this on with much talent.

Singularities.

A large part of this volume is formed of essential papers on singularities of mapsnotably on their topological stability and on stratifications-annotated and commented upon by Teissier and Trotman, experts in the field. This includes La stabilité topologique des applications polynomiales, which sketches already the celebrated isotopy theorems of Thom-Mather, then Propriétés différentiables locales des ensembles analytiques (d'après H. Whitney), excellent exposition of fundamental results of Whitney on the stratification of analytic sets and a draft of the theory of stratified sets; this is further developed in Local topological properties of differentiable mappings, astonishingly clairvoyant, and followed by On some ideals of differentiable functions, focusing on the difference between the differentiable and the analytic.

There next comes a "big part," Ensembles et morphismes stratifiés which, after inspiring Mather and a pleiad of other experts, continues after fifty years to be the source of active research. The bifurcation subset of a space of maps introduces the absolutely fascinating idea of a natural stratification of function spaces, concerning which much remains to be done.

The volume concludes, in a way closer to the two unpublished papers of volume I, by Un résultat sur la monodromie, commented on by Norbert A'Campo with Teissier, and by an unpublished manuscript on the monodromy11 which, despite a fatal error, contains important and beautiful ideas. For the anecdote, Thom's article had been accepted in a prestigious journal, which had previously refused at the same time the paper where the young A'Campo gave a counterexample! Happily, Thom withdrew his text and the work of A'Campo, which surprised the most distinguished, established his reputation.

Catastrophes. Preceded by a knowledgeable introduction of Petitot, this part contains five articles.

A dynamic theory for morphogenesis, "lost" but recovered by Tadashi Tokieda, presents the theory for the first time, insisting on its mathematical aspects-this is in particular the first appearance of universal unfoldings; the commentary led us to interrogate Mather, whose response, of great interest, is partially reproduced.

Thom presents his ideas in more detail in Une théorie dynamique de la morphogénèse, a fundamental article followed by correspondence with the great biologist C. H. Waddington; the whole is commented on by Sara Franceschelli and Petitot, who begin with the relation with The chemical basis of morphogenesis of Alan Turing (1952).

Then comes Topological models in biology, then A mathematical approach to morphogenesis: archetypal morphologies and Topologie et linguistique, which illustrates the scale of the project.

Varia. The volume contains other "concrete" applications of the theory of singularities: Sur la théorie des enveloppes, written in 1960 to clean up a rather "dirty" domain, is not very readable-we have tried to remedy this in part without departing from Thom's ideas. Sur les variétés d'ordre fini sketches notably the proof that, save for rare exceptions, a compact submanifold M of dimension n in R n+µ cuts every affine µ-plane P in a finite number of points, bounded as P varies.

More ambitious, the article Les symétries brisées en physique macroscopique et la mécanique quantique is confronted by Valentin Poénaru with the later progress in physics. Like the talk Travaux de Moser sur la stabilité des mouvements périodiques on "KAM theory," it bares witness to the variety of interests of Thom and to his exceptional insight.

The article Sur l'homologie des variétés algébriques réelles, commented on by Ilia Itenberg, is devoted to what is often called since the Smith-Thom inequality, very important in real algebraic geometry. This work dedicated to Marston Morse depends on Morse theory, of which a "foliated" version is proposed in Généralisation de la théorie de Morse aux variétés feuilletées.

Finally Jets de Liapunov, commented on by Krzysztof Kurdyka, examines the implications of the existence of a Liapunov function for a vector field in the neighbourhood of a point.

Volume III

Respecting the name of the collection, this last volume12 essentially contains mathematical writings, a small minority of Thom's publications after 1971, but often remarkable. Let us point out some of them. The first text, the Fermi lectures Modèles mathématiques de la morphogénèse13 is a clear exposition of the stakes of catastrophe theory and a lot of new Mathematics, which have had an important progeny. A critical rereading was therefore in order.

Then comes Sur le cut-locus d'une variété plongée, a paper which is very rich mathematically (it takes the functional viewpoint of The bifurcation subset of a space of maps) and even beyond Mathematics-one finds there a quite up-to-date model of visual perception.

Phase transitions as catastrophes which is mathematically substantial14 proposes in particular a proof of the Gibbs phase rule based on the generic structure of Maxwell sets. 15 This notion was at the heart of the cut locus paper.

Sur les équations différentielles multiformes et leurs intégrales singulières provides among other things a marvellous introduction to singularities, contact geometry and Pfaffian systems.

Symmetries gained and lost analyses symmetry breakings, in connexion with the work of the physicist Louis Michel and as a continuation of the paper of Volume II.

Introduction à la dynamique qualitative, a mathematico-historico-philosophical text à la Thom, evidences his long standing interest for this domain which was for a time not so popular in France, a situation which his seminar around 1970 had remedied, contributing to the orientation of future "leaders" in the subject, such as Michel Herman.

Gradients of analytic functions, absent from the CD-ROM and unearthed by our Iranian colleague Massoud Amini, notably states the "gradient conjecture" on the limits at a singular point of the tangents to integral curves of an analytic gradient. 16 Enthusiastic commentary by Kurdyka.

Tectonique des plaques et théorie des catastrophes is the trace of a daring incursion into Claude Allègre's territory.

The Les Houches lectures Mathematical concepts in the theory of ordered media is a follow-up to the work of Kléman and Toulouse on dislocations in crystalline media.

The note with Yannick Kergosien Sur les points paraboliques des surfaces, completed and corrected with Thomas Banchoff, deals with the apparent contours of a "generic" surface in euclidean space when the direction of projection varies.17

The two notes with Peixoto on Le point de vue énumératif dans les problèmes aux limites pour les équations différentielles ordinaires, illustrate in particular the very enlightening vision that Thom had of "well posed" problems: they are those in which the d-dimensional manifold of solutions of the differential system is transversal in the function space to the d-codimensional variety defined by the limit conditions of the problem.

Quid des stratifications canoniques? evidences again their extreme importance in the eyes of Thom.

Finally, the obituaries of Morse and Whitney written for the Académie des sciences are obviously of great interest, as well as La théorie des jets et ses développements ultérieurs published with Ehresmann's complete works.

We also write a little about "the catastrophe of catastrophe theory", which took place in 1978 and has without any doubt contributed to Thom distancing himself from mathematics, to its detriment.

Catastrophes, around 1975 c Archives de l'IHÉS -Droits réservés Appendix: a modest overview of transversality according to Thom

Transversality plays an essential role in Thom's work, both in mathematics and in catastrophe theory where a key idea is that one can only observe the phenomena that are stable under perturbation. The transversality lemma in jet spaces18 amazed its first "guinea pig" Whitney who at first could not believe this statement which put under the same umbrella almost all general position arguments. That is why we have thought it would be useful to present an overview of it.19

First elementary version.

The following facts seems to correspond to our intuition. If, in R 3 , a point p lies on a surface S a jiggling20 will push p outside of S. If a (compact) curve C is tangent to S at one or several of its points, a jiggling will result in that either C and S are disjoint or, at the possible common points, the tangent to C and the tangent plane to S will intersect transversally. Finally, if S is a (compact) surface which at some points shares a tangent plane with S, a jiggling of S will put it in a transversal position with respect to S; at the possible common points the tangent planes will be secant; in fact-and this is less intuitive-if this last condition on the tangent planes is satisfied at each point of the intersection S ∩ S , this intersection will be a smooth curve (this last fact is guaranteed by the implicit function theorem).

Thom's first statement on transversality is a generalization to all dimensions and all manifolds of those three examples. But one has still to specify what a "jiggling" is. It is a diffeomorphism g of the ambient manifold M , arbitrarily close to the identity (in the C ∞ topology) which one applies to a (compact) submanifold N in order to make it transversal to another submanifold N of M : at any point of g(N ) ∩ N the two tangent spaces T x g(N ) and T x N generate T x M . In what follows all manifolds will be assumed to be C ∞ (we shall say "smooth").

Elementary transversality theorem (Thom).

The property that g(N ) is transversal to N is generically satisfied for g in the group Diff

k (M ) of C k diffeomorphisms of M if k > max(0, dim N - codim N ).
Here, a property is said to be generic if it is satisfied at all points of a countable intersection of dense open subsets of Diff k (M ).

Remarks. 1) If N is compact, transversality is an open property in the C 1 topology and therefore in all finer topologies such as the C ∞ on the group Diff ∞ (M ). Moreover, under this assumption, genericity is also true in the subgroup Diff ∞ c (M ) of diffeomeorphisms with compact support. 2) The same statement is still valid if N is singular (presence of double points, lack of tangent space) or even the image in M of a manifold N by a C k map.

Proof of Thom's theorem based on theorems of A. Morse and A. Sard.

Theorem (Morse-Sard

). Let X be a manifold. If f : X → R q is C k , k > max(0, dim X -q),
then, almost every y in R q , in the sense of the Lebesgue measure, is a regular value of f , which means that for every x ∈ f -1 (y) the differential of f at x is of rank q.

For instance, if dim X < q, almost no y is in the image of f . By the implicit function theorem, the inverse image of a regular value is a C k submanifold of codimension q.

The inequality on k in the statement above is necessary. For example, Whitney has built a function f : R 2 → R of class C 1 and an arc A in the plane where f is not constant but along which the differential of f is identically zero (the arc A is continuous but of course not C 1 ); the interval f (A) contains no regular value of f although it is of positive length.

Here is how Thom uses Sard's theorem to prove his first transversality theorem. With its notations, let us consider the case k = ∞ and N compact and concentrate on the density part. Since Diff ∞ c (M ) is a group one easily shows that it suffices to prove that id M can be approximated in the C ∞ topology by g ∈ Diff ∞ c (M ) such that g(N ) is transversal to N . Thom's idea is to build, for a certain integer r, a map

G : R r × M → M , (v, x) → G(v, x) of class C ∞ , such that:
(1) G(0, x) = x for all x ∈ M ;

(2) For all small enough fixed v, g v := G(v, -) is a diffeomorphism of M with compact support;

(3) ∂ v G(0, x) is of rank n = dim M for all x ∈ N .

We shall think of G as a family {g v } v∈R r depending on r real parameters. Let us assume we have such a family and prove Thom's theorem. Since the rank condition is open, condition (3) implies that ∂ v G(v, x) is of rank n for all (v, x) in a neighborhood U of {0} × N . Let us still denote by G its restriction to U . With the rank condition, the implicit function theorem tells us that W :

= G -1 (N ) is a submanifold of U of class C ∞ .
Let Π : W → R r be the projection onto the parameter space. One now applies Sard's theorem to Π. Almost every (r-dimensional) parameter v is a regular value of Π. An easy lemma of linear algebra shows that for such a v, the diffeomorphism g v pushes N to a position transversal to N .

There remains only to build the family G. Let us begin with the case M = R n . We take as parameter space the vector space R n of translations of R n . By a well chosen partition of unity, one extends the translation of the compact N by the vector v into a diffeomorphism g v with compact support which is id M if v = 0. Conditions (1) -(3) are satisfied.

In the case where M is a manifold, using the compactness of N we cover it by a finite collection {B i } s i=1 of balls, each of which is contained in a coordinate chart O i of M . One chooses the translations v i in the coordinates of O i small enough for B i + v i to remain in a fixed compact subset of this chart. One easily extends this translation into a diffeomorphism g i,v i of M with compact support which is id

M if v i = 0. Let v := (v 1 , ..., v s ). Set g v = g 1,v 1 • • • • • g s,vs
. One easily verifies that conditions (1) -( 3) are satisfied, which ends the proof.

Transversality and jet spaces.

It is obviously in the jet spaces of order r ≥ 2 that Thom will give the deepest applications of the transversality theorem to singularities of differentiable mappings. However, one can already realize the innovative aspect of this theorem with jets of order 1 (or 1-jets) of real functions of one real variable.

In order to start with a compact manifold, let us look at the 1-jets of functions on the interval [0, 1], which is not exactly a manifold, but a manifold with boundary. The boundary creates no difficulty for understanding the space of real valued C ∞ functions.

If

f : [0, 1] → R is C k , k ≥ 1, we define its jet of order one j 1 f (a) = (a, f (a), f (a)) at any point a ∈ [0, 1]. 21 The map j 1 f is a section (that is, a right inverse) of the projection (x, y, z) → x of J 1 ([0, 1]) to [0, 1]; thus, it identifies with its image which is a C k-1 submanifold of J 1 ([0, 1]). Given a ∈ [0, 1], any P ∈ R[T ] of degree one is equal to j 1 a f , where f (x) = P (x -a). However a global section s(x) = (x, y(x), z(x)) is of the form j 1 f (x) if and only if z(x) = dy dx (x) (integrability condition).
This condition is quite rarely satisfied. Thus, if Σ is a submanifold of J 1 ([0, 1]), the proof given for the transversality theorem (elementary version) will produce from j 1 f at best a section transversal to Σ, with practically no chance to produce an integrable section.

We remember that everything lies in the choice of the family G of maximal rank with respect to the parameters. It is Thom's clairvoyance which made him see that the same proof worked in the subspace of integrable sections if we take as parameter space the (finite dimensional) space of degree one real polynomials, here in one variable. At the same time we shall replace the diffeomorphism group of J 1 ([0, 1]) by its subgroup, which one might call the gauge group consisting of those diffeomorphisms which preserve not only each fiber J 1

x ([0, 1]) but also the distribution of contact planes, i.e., the kernels of the differential form dy -zdx. 22 If P is a polynomial function of degree one, the translation (in the fibers) (x, y, z) → (x, y + P (x), z + P (x)) in other words j 1 f (x) → j 1 (f + P )(x) is a gauge transformation, sending each integrable section j 1 f to the integrable section j 1 (f + P ). This family of translations verifies conditions (1) -(3), once truncated to obtain (2). The same approach works in any dimension of base manifold and for any jet order. 23 Thom thus obtains the following theorem:

Transversality theorem in a jet space (Thom).

Let N be an n-dimensional manifold, let J r (N ) be its space of r-jets of real functions and let Σ be a submanifold of J r (N ) of codimension q and of class C , > max{0, n -q}. For any f ∈ C k (N, R), k -r > max{0, n -q}, generically for g in the C k-r gauge group, the submanifold g(j r f (N )) is transversal to Σ.

Here is a first application to Morse functions. At the beginning of the 1930's, Marston Morse showed the importance of the functions which now bear his name for the understanding of manifolds. Recall that a Morse function f : N → R is a function, of class at least C 2 , whose critical points (points where df = 0) are all non degenerate (in coordinates, the matrix of second partial derivatives is of maximum rank).

Corollary. Morse functions of class C k on N are dense in C k (N, R) for k ≥ 2. Indeed, a simple computation shows that f is a Morse function if and only if j 1 f is transversal to the zero section Σ := {j 1 ϕ(x) : d x ϕ = 0} of the projection j 1 ϕ(x) → j 0 ϕ(x) := (x, ϕ(x)). According to the theorem (here with q = n and r = 1), generically j 1 f is transversal to Σ, hence the density.

Numerous discussions between them about foliated manifolds had then led to the "concrete" part (analytic foliations) of the still famous thesis on which Haefliger worked under the direction of Charles Ehresmann, rarely present in Strasbourg during this period.

Thanks to the efforts of its director, Jean Pierre Bourguignon, and to the unfailing enthusiasm of Michèle Porte, director of this project begun in 1996 and to which Thom had actively collaborated.

Those for the books reproduced were ceded by their publishers only for a limited number of copies.

This choice, even if it has contributed to much delay in the publication of volume I, was imposed all the more since Thom often published in journals or collections of articles which have become almost impossible to find.

We have unfortunately not been able to take into account in volume I an inestimable document which is part of the archives given since then by the Thom family to the library of the École Normale Supérieure: chapter 6 on cobordism, written by Thom and excluded by Cartan, for whom it seemed to not yet make sense.

Mildly corrected with respect to volume I, which sometimes changes the numbering.

Except for the Fermi Lectures Modèles mathématiques de la morphogénèse given at the Scuola Normale Superiore of Pisa in April 1971, kept for volume III as the lectures in Bonn were kept for volume II.

Mathematicians of the calibre of John Mather and Vladimir Arnold worked on this, but many exciting problems remain open.

In a modern typography of Duco van Straten.

In preparation, we would like to see it published at the end of 2020.

It provides the first three chapters of the book with the same title published in 1974 in the pocket collection 10-18.

and bold, sometimes a little too much so!

The set of values of parameters for which a function depending on those parameters attains its minimum in several points (if one counts them with multiplicities).

A version of it has been proved by Kurdyka, Mostowski and Parusiński in a paper published in 2000 in the Annals of Math.

The case of a single projection was very well treated in the paper Sur les équations différentielles multiformes et leurs intégrales singulières.

See Thom's paper Un lemme sur les applications différentiables (1956), reproduced and commented in volume I.

Marc Chaperon had treated it very differently in issue 65 of the Gazette de la SMF, with many other applications.

A generic small movement; we borrow this term from W. Thurston.

In other wordsj 1 f (a) = (a, j 1 a f ), where j 1 a f ∈ R[T ] is the Taylor polynomial f (a) + f (a)T .

A "contactologist" tradition exchanges y and z, giving dz -ydx.


Sarah Visage

Pauline Souloumiac

Bertrand Maillot

Matthias Rosenau

Ehsan Kosari

Data Description

This dataset provides friction data from ring-shear tests (RST) on twice broken rice used in the GEC Laboratory in CY Cergy Paris University in stick-slip experiments. They were obtained by Sarah Visage as part of her doctoral training (funded by the ANR DISRUPT programme) during an invitation at the Helmholtz Laboratory for Tectonic Modelling (HelTec) at the GFZ German Research Centre for Geosciences in Potsdam. Like any granular material, the twice broken rice is characterized by several internal friction coefficients μ and cohesions C, classicaly qualified as dynamic, static, and reactivation coefficients. In adition, since the rice exhibits a stick slip behaviour, the various shear -velocity or shear-displacement curves exhibit high frequency oscillations and we therefore define maximum, minimum, and mean values corresponding respectively to the curve peaks, curve troughs and smoothed curve.

Material tested

The tested material consists of twice broken rice. The rice originates from Cambodia and is sold by the brand Perle d'Asie (Fig. 1).

Measurement procedure

The data presented here are the outcomes of ring shear tests (RST) using a SCHULZE RST-01.pc [START_REF] Schulze | Time-and Velocity-Dependent Properties of Powders Effecting Slip-Stick Oscillations[END_REF][START_REF] Schulze | Powders and bulk solids[END_REF] at the Helmholtz Laboratory for Tectonic Modelling (HelTec) of the GFZ German Research Centre for Geosciences in Potsdam. The RST-01 is designed to measure friction coefficients μ and cohesions C in loose granular materials at the low confining pressures and shear velocities found in sandbox experiments. In this apparatus, a granular layer is sheared internally at constant normal stress σN and shear velocity v, while the shear force and the lid displacement normal to the shear plane are measured continuously. For more details see [START_REF] Klinkmüller | Properties of granular analogue model materials: A community wide survey[END_REF] and [START_REF] Ritter | Scaling the sandbox-Mechanical (dis) similarities of granular materials and brittle rock[END_REF].

Sample preparation and test conditions

Each sample is prepared by the same person and measured following the same protocol, rice has been poured into the shear cell and not densified thereafter. The resulting bulk density has been 911 ± 11 kg/m3. The measurements presented here correspond to the internal friction, i.e. shearing inside the material. Normal force, shear force, shear velocity and lid displacement are measured at 100 Hz. Air conditionning ensures constant temperature and humidity during all the measurements (temperature: 23•C, humidity: 45%). The rice samples are labbeled GFZ-ID 498-01 to 498-18 for the Ring-Shear Tests (RST), 499-01 for the Velocity Stepping Tests (VST), and 500-01 for the Slide-Hold-Slide tests (SHS).

RST (Ring-shear test) procedure

During the RST a shear velocity of v = 30 mm.min-1 and a normal load as defined below are imposed while the shear force and the lid displacement are monitored. 18 individual tests are done at normal stresses of σN = 500, 1000, 2000, 4000, 8000, and 16000 Pa (with 3 repetitions per stress level). At the beginning of the test, the material is sheared during 5 minutes. During this period the shear stress τ reaches a peak (= peak friction) and then drops to a plateau indicating the shear has localized into a permanent shear zone (= dynamic friction). The sample is then unloaded by shortly reversing the rotation and immediately re-sheared for 5 minutes during which the shear stress τ is expected to reach a second peak (= reactivation friction) simulating reactivation of an existing shear zone before returning to the plateau (Fig. 2).

VST (Velocity stepping test) procedure

To determine the dependence of dynamic friction on the shear velocity v, a velocity stepping test (VST) is performed. Shear velocities ranging from 0.1 to 30 mm.min -1 and a normal stress of σN = 2000 Pa are imposed. The velocity is systematically increased in logarithmic steps. Durations at each values of the velocity are chosen to reach a displacement of 30 mm in each step (Table 1). This displacement is sufficient to ensure that the plateau of the dynamic friction is reached in each step. In order to initiate the shear zone, a first step at 30mm for 3min is performed (Fig. 4).

Shear velocity [mm.min -1 ]

Duration [hh:mm:ss] 30 00:03:00 0.1 05:00:00 0.3 01:40:00 1 00:30:00 3 00:10:00 10 00:03:00 30 00:01:00

SHS (Slide-Hold-Slide) procedure

To determine the dependence of reactivation friction on the holding time, a slide-hold-slide test (SHS) is performed. The shearing phases, or slide steps, last for three minutes, at the constant shear velocity of 30 mm.min -1 . During the holding phases, the shear velocity is zero. Durations of holding start at 5 hours and are reduced by a factor three down to ten seconds (

Analysis method

RST analysis: Friction coefficients and cohesion

After converting forces to stresses, by dividing the force in Newton by the area of the shear plane (226,19 cm 2 ) including correcting shear stress for the lever geometry and time to displacement (with the velocity of 30 mm.min -1 ), three characteristic values have been picked manually from the resulting shear stress curves: peak, dynamic, and reactivation values:

1) The shear strength τ at peak friction corresponding to the first peak in the shear curve reflecting hardening-weakening during strain localization 2) the shear strength τ at dynamic friction corresponding to the plateau after localization and representing friction during sliding 3) the shear strength τ at reactivation friction corresponding to the second peak and representing static friction during reactivation of the shear zone. However, in the twice broken rice we do not observe a reactivation peak.

Since the curve exhibits high frequency oscillations -because of the stick-slip behaviour of the rice -, we distinguish maximum, minimum and average values, as shown in Figure 2. We performed a linear regression of these friction data through all data pairs of shear strength τ and normal stress σN. The slope of the linear regression corresponds to the friction coefficient μ and the y-axis intercept, to the cohesion C (Fig. 3). This method assumes that the material behaves strictly as a Mohr-Coulomb material, i.e. has an affine failure envelope.

VST analysis: Rate-dependencies of dynamic friction

To calculate "a-b", we first pick the peaks and troughs of the apparent coefficient of friction (i.e., neglecting cohesion) μ=τ/σN calculated from the velocity stepping test (VST) (with the matlab code "velocity-weakening Fitting" by Ehsan Kosari ( 2021)). These values represent the maximum and minimum dynamic frictions. The slope of a linear regression through the data approximates the dependency of the dynamic friction on the shear velocity v (i.e., its rate-dependency, parameter (a-b) [START_REF] Scholz | Earthquakes and friction laws[END_REF]). We have plotted peaks, troughs, and as well as their average values as a function of log(v) and fitted three natural logarithm lines to the data. The slop of the fitted line to the average values is considered as "a-b". To derive a-b the slope derived in log(v) space has to be divided by ln(10) (in the table 7 the values are divided by ln( 10)).

SHS analysis: Frictional healing

From the SHS time series data (Fig. 6), the coefficient of apparent reactivation friction μ should be plotted as a function of log(THold), then the slope of a linear regression through the data gives an approximation of the dependence of static friction on holding time or frictional healing (i.e. the parameter b). However in the rice we don't observe a reactivation peak, so we decided to use the dynamic friction in the slip phase which shows a similar dependence to the previous hold phase which we plotted as a function of log(THold) (Fig. 7, with an adapted version of the matlab code "velocityweakening Fitting" by Ehsan Kosari ( 2021)). In the same way as for the RST and VST tests, we take the maximum, minimum and average values of the friction during the sliding phases. To derive b, the slope derived in log(v) space must be divided by ln(10) (in the table 7 the values are divided by ln( 10)).

File description

The following files exist in the zip folder "2022-001_Visage-et-al_Data":

(i) RST raw data ("ID [f][date time].tdms") (ii)
RST shear curve data ("ID ts.txt"; example Table 3)

(iii)
RST smoothed shear curve plot ("ID ts.pdf")

RST friction data ("ID peak.txt", "ID dynamic.txt", "ID reactivation.txt"; example Table 4)

(v)
RST friction plot and linear least-squares method data ("ID linregr.pdf", "ID fricstd.txt"; example Figure 3)

(vi)
VST raw data ("ID VST.tdms")

VST plot ("ID vst.pdf"; example Figure 4)

(viii)
SHS raw data ("ID SHS.tdms")

(ix)
SHS plot ("ID SHS.pdf"; example Figure 6)

(x)
MATLAB code "velocity-weakening Fitting"

RST shear curve data

RST shear curve data are derived from RST raw data (i) and given as (ii) time series (ts) data in .txtformat ("ID ts.txt", table 3) and visualized as (iii) shear stress τ versus shear displacement d plots ("ID ts.pdf", figure 2).

RST friction data and analysis

RST friction data are given as (iv) data pairs (normal stress σN and shear strength τ; Table 4) for peak, dynamic and reactivation friction in .txt-format ("ID peak.txt", "ID dynamic.txt", "ID reactivation.txt").

They are visualized by (v) plotting into Mohr Space (normal stress σN vs. shear stress τ) including a linear regression ("ID linregr.pdf"; figure 3).

Normal stress [

VST data and analysis

VST raw data (vi) are visualized (vii) by plotting shear velocity v and dynamic friction (simplified as shear stress τ divided by normal stress σN) against the shear displacement d ("ID vst.pdf"; 4). For analysis, dynamic friction is plotted against the shear velocity v including a logarithmic curve fit, the slope divide by ln(10) of which reflects the rate dependency of dynamic friction (the parameter a-b). An example of the structure of VST data is given in Table 5.

SHS data and analysis

SHS raw data (viii) are visualized (ix) by plotting time against static friction ("ID shs.pdf"; Fig. 6). For analysis, static friction is plotted against the holding time including a logarithmic curve fit, the slope divide by ln(10) of which reflects the paramater b of static friction. An example of the structure of SHS data is given in

Results

Parameter

  Material properties of analogue velocity-weakening material used for seismotectonic analogue modeling of strike-slip seismic cycle: the case of twice-broken rice Sarah Visage, Pauline Souloumiac, Bertrand Maillot, Matthias Rosenau, Ehsan Kosari To cite this version: Sarah Visage, Pauline Souloumiac, Bertrand Maillot, Matthias Rosenau, Ehsan Kosari. Material properties of analogue velocity-weakening material used for seismotectonic analogue modeling of strikeslip seismic cycle: the case of twice-broken rice. Geosciences data publisher, 2022, 10.5880/fidgeo.2022.001. hal-03838426 Table of content 1. Licence ............................................................................................................................................. 1 2. Citation ............................................................................................................................................ 1 3. Data Description .............................................................................................................................. 2 3.1. Material tested ........................................................................................................................ 3 3.2. Measurement procedure ........................................................................................................ 3 3.2.1. Sample preparation and test conditions ....................................................................

Figure 1 :

 1 Figure 1: Photo of twice broken rice grains with binoculars.

Figure 2 :

 2 Figure 2: Shear curve plot (498-04) at a normal load of 8000Pa and at a shear velocity of 30 mm min-1 with the location of each values analysed.

Figure 3 :

 3 Figure 3: Friction plot (498-01 to 18). Plot of all data pairs in the Mohr space (shear stress τ vs. normal stress σN ) including curves of the corresponding linear regression (shear velocity = 30 mm min-1). The results are shown in table 1.

Figure 4 :

 4 Figure 4: The VST data (499-01). The shear velocity v (red curve) and the coefficient of friction (shear stress τ/normal stress σN) (blue curve) at normal load of σN =2000Pa.

Figure 5 :

 5 Figure 5: Dynamic friction (peaks = blue, mean = yellow and troughs = red) as a function of the logarithm of the shear velocity v, with the linear regressions. The results are shown in table 1.

Figure 6 :

 6 Figure 6: The SHS data (499-01). The static friction (shear stress τ /normal stress σN ) at a normal load of 2000 Pa and a velocity of 30mm.min-1 during the sliding.

Figure 7 :

 7 Figure 7: Static friction (peaks = blue, mean = yellow and troughs = red) as a function of the holding time on a logarithmic scale with the linear regression. The results are shown in table 1.

  Ring-shear test) procedure ...................................................................................... 4 3.2.3. VST (Velocity stepping test) procedure ........................................................................... 4 3.2.4. SHS (Slide-Hold-Slide) procedure .................................................................................... 4 3.3. Analysis method ...................................................................................................................... 5

		..... 3
	3.2.2.	RST (7
	4.4. SHS data and analysis .............................................................................................................. 8
	5. Results ............................................................................................................................................. 8
	6. References ..................................................................................................................................... 12

3.3.1. RST analysis: Friction coefficients and cohesion ............................................................. 5 3.3.2. VST analysis: Rate-dependencies of dynamic friction ..................................................... 6 3.3.3. SHS analysis: Frictional healing........................................................................................ 6 4. File description ................................................................................................................................ 6 4.1. RST shear curve data ............................................................................................................... 7 4.2. RST friction data and analysis .................................................................................................. 7 4.3. VST data and analysis ..............................................................................................................

Table 1 :

 1 Logarithmic steps of the shear velocity v in a VST including the duration of each step.

Table 2 )

 2 (Fig.6). A normal stress of σN = 2000 Pa is imposed throughout the test.

		Duration [hh:mm:ss]
	Shear	00:03:00
	Hold time	05:00:00
	Shear	00:03:00
	Hold time	01:40:00
	Shear	00:03:00
	Hold time	00:30:00
	Shear	00:03:00
	Hold time	00:10:00
	Shear	00:03:00
	Hold time	00:03:00
	Shear	00:03:00
	Hold time	00:01:00
	Shear	00:03:00
	Hold time	00:00:30
	Shear	00:03:00
	Hold time	00:00:10
	Shear	00:03:00

Table 2 :

 2 Steps of the SHS test including the duration of each step.

Table 3 :

 3 Example of shear curve time series data (498). First line is header. First column is time (in s). Columns 2-19 are measured shear forces (in N) for corresponding normal stresses as specified in the header of the respective columns (6 stress values from 500 to 16 000 Pa, with three repetitions at each stress value).

	Time [s]	Normal stress [Pa]:	1000	2000	…
		500
	0.1	46.98	81.60	88.11	…
	0.2	…	…	…	…
	…	…	…	…	…

Table 4 :

 4 Example of friction data (498, peak). First column is normal stress σN (in Pa). Column 2-4 is shear strength τ (in Pa) for the high, mean and low values.

Table 5

 5

	: Example of VST data (499). First column is time (in s). Second column is shear velocity (in mm s-1). Third and fourth
	columns contain normal and shear force (in N)

Table 6 :

 6 Table 5 and the holding time is given in Table 2. Example of VST data (500). First column is time (in s). Second and third columns contain normal and shear force (in N)

	Time [s]	Normal force [N]	Shear force [N]
	0.1	1909.4	2.6074
	0.2	…	…
	…	…	…

Table 7 :

 7 Summary of RST, VST and SHS data.