Thomas Rey
email: trey@cscamm.umd.edu

Giuseppe Toscani
email: giuseppe.toscani@unipv.it

LARGE-TIME BEHAVIOR OF THE SOLUTIONS TO ROSENAU TYPE APPROXIMATIONS TO THE HEAT EQUATION

Keywords: 2010 Mathematics Subject Classification. Primary: 35K05, 35Q20; Secondary: 35B40, 65M06 Large time behavior, Heat equation, Fourier metrics, Central differences scheme, Rosenau approximation, non-local model

, by showing that the solution to this approximation approaches the fundamental solution of the heat equation at a sub-optimal rate. The result is valid in particular for the central differences scheme approximation of the heat equation, a property which to our knowledge has never been observed before.

1. Introduction 1.1. Rosenau Regularization of the Chapman-Enskog Expansion. In his seminal work [START_REF] Rosenau | Tempered diffusion: A transport process with propagating fronts and inertial delay[END_REF], Rosenau has proposed a regularized version of the Chapman-Enskog expansion of hydrodynamics. This regularized expansion resembles the usual Navier-Stokes viscosity terms at low wave-numbers, but unlike the latter, it has the advantage of being a bounded macroscopic approximation to the linearized collision operator. The model is given by the scalar equation

(1.1) ∂f ∂t + ∂Ψ(f ) ∂v = -εξ 2 1 + ε 2 m 2 ξ 2 f (ξ) ∨ ,
where f (ξ) denotes the Fourier transform of f (v), while f (ξ) ∨ denotes the inverse Fourier transform.

The operator on the right hand side looks like the usual viscosity term εf vv at low wave-numbers ξ, while for higher wave numbers it is intended to model a bounded approximation of a linearized collision operator, thereby avoiding the artificial instabilities that occur when the Chapman-Enskog expansion for such an operator is truncated after a finite number of terms.

One of the advantages of this regularization is that it is equivalent in small frequencies to the heat equation (and then regularizes initial data), whereas high frequencies behave as an absorption term. In particular, information will travel at finite speed, and tools from hyperbolic equations can be used. Note that the right side of (1.1) can be written as

(1.2) -εξ 2 1 + ε 2 m 2 ξ 2 f (ξ) ∨ = ε (εm) 2 1 1 + ε 2 m 2 ξ 2 f (ξ) -f (ξ) ∨ = 1 mε [M ε * f -f ] ,
where ε = mε, * denotes convolution and

(1.3) M γ (v) = 1 2γ e -|v|/γ
is a non-negative function satisfying M γ L 1 = 1. In other words, the Rosenau approximation consists in substituting the linear diffusion equation 2 (v, t) with the linear kinetic equation

(1.4) ∂g ∂t (v, t) = σ 2 ∂ 2 g ∂v
(1.5) ∂g ∂t (v, t) = σ 2 ε 2 [M ε * g(v, t) -g(v, t)
] in which the "Maxwellian" M ε [START_REF] Cercignani | The Boltzmann Equation and its Applications[END_REF] is given by (1.3).

Equation (1.1) has then been studied by S. Schochet and E. Tadmor in [START_REF] Schochet | The Regularized Chapman-Enskog Expansion for Scalar Conservation Laws[END_REF] in the context of vanishing viscosity solutions of scalar conservation laws. H. Liu and E. Tadmor then studied in [START_REF] Liu | Critical Thresholds in a Convolution Model for Nonlinear Conservation Laws[END_REF] the Burgers equation with a Rosenau-like nonlocal viscosity term and proved the existence of a so-called critical threshold, namely the existence of a critical value of the total variation of the initial condition triggering or not a finite-time blow-up of the solutions. This nonlocal term was also used by C. Rohde to model capillarity effects close to fluid-vapour phase transitions in Navier-Stokes equations [START_REF] Rohde | On local and non-local Navier-Stokes-Korteweg systems for liquid-vapour phase transitions[END_REF]. Interestingly, this non-local approach seems more physically relevant than the original Korteweg's idea to use the heat operator for describing such phenomenons.

More recently, it has been replaced in a general context of system of hyperbolic balance laws with non-local source term by R. Colombo and G. Guerra in the two companion papers [START_REF] Colombo | Hyperbolic balance laws with a non local source[END_REF][START_REF] Colombo | Hyperbolic balance laws with a dissipative non local source[END_REF]. They showed in particular the well posedness in L 1 (globally in time) of the Cauchy problem for (1.1) together with uniform stability estimates, in the framework of viscosity solutions. Finally, C. Besse and T. Goudon obtained recently in [START_REF] Besse | Derivation of a Non-Local Model for Diffusion Asymptotics -Application to Radiative Transfer Problems[END_REF] a system of macroscopic equations with a nonlocal Rosenau term by computing the diffusion limit of a linear space inhomogeneous kinetic equation. A similar approach was also used by T. Goudon and M. Parisot in [START_REF] Goudon | On the Spitzer-Härm Regime and Nonlocal Approximations: Modeling, Analysis, and Numerical Simulations[END_REF] on a kinetic equation describing the behavior of a multi-species gas of charged particles to obtain a macroscopic model with a Rosenau-like nonlocal term.

Despite the previous studies, it is not completely clear and well established from a mathematical point of view if the correction proposed by Rosenau is a good approximation to the linear heat equation. In particular, it is not clear whether or not the large-time behavior of the solution to the Rosenau approximation agrees with the large-time behavior of the linear diffusion equation.

The aim of this article is to give an answer to the previous question, and to underline that the Rosenau approximation can be viewed as a particular case of a general approximation to the heat equation by means of a linear kinetic equation of type (1.5), provided the background density M ε is a probability density function of zero mean and variance proportional to ε 2 . In particular, it will be shown that, in a certain metric equivalent to the weak*-convergence of measures, the distance between the solution to the heat equation and the solution to the kinetic equation can be bounded uniformly in terms of ε and t, provided the background density has a sufficiently high number of moments (typically more than two).

The plan on the article is as follows. In Section 2, using tools of the kinetic theory of rarefied gases, we will introduce a possible derivation and the main features of the Rosenau approximation with a general kernel. This allows to describe well known models such as the central differences schemes for the heat equation. In Section 3, this kinetic formulation will be used to obtain explicit solutions to the Rosenau equation (1.5), using both Fourier transform and Wild sums. Last, in Section 4, we will investigate the large time behavior of the solutions to (1.5), showing that the convergence towards the fundamental solution occurs, in a suitable Fourier-based metric (see next subsection), at a suboptimal rate (compared to the heat equation). Finally, we will combine these results to show that strong convergence in L 1 towards the fundamental solution to the heat equation is obtained after a suitable regularization of the Rosenau equation, obtained by discarding its singular part.

1.2. Functional Framework. Before entering into the main topic of this paper, we list below the various functional spaces used in the following. For p ∈ [1, +∞) and q ∈ [1, +∞), we denote by L p q the weighted Lebesgue spaces

L p q := f : R → R measurable; f p L p q := R |f (v)| p (1 + v 2 ) q/2 dv < ∞ .
In particular, the usual Lebesgue spaces are given by

L p := L p 0 . Moreover, for f ∈ L 1
q , we can define for any k ≤ q the k th order moment of f as the quantity

m k (f ) := R f (v) |v| k dv < ∞.
For s ∈ N, we denote by W s,p the Sobolev spaces

W s,p :=    f ∈ L s ; f p W s,p := |k|≤s R f (k) (v) p dv < ∞    . If p = 2 we set H s := W s,2 .
Given a probability density f , we define its Fourier transform F v (f ) by

F v (f )(ξ) = f (ξ) := R e -i ξ v f (v) dv, ∀ξ ∈ R.
The Sobolev space H s can equivalently be defined for any s ≥ 0 by the norm

f H s := F v (f ) L 2 2s
.

The homogeneous Sobolev space Ḣs is then defined by the homogeneous norm

f 2 Ḣs := R |ξ| 2s f (ξ) 2 dξ.
Finally, we introduce a family of Fourier based metrics in the following way: given s > 0 and two probability distributions f 1 and f 2 , their Fourier based distance d s (f 1 , f 2 ) is the quantity

d s (f 1 , f 2 ) := sup ξ∈R\0 f 1 (ξ) -f 2 (ξ) |ξ| s .
This distance is finite, provided that f 1 and f 2 have the same moments up to order [s], where, if s / ∈ N, [s] denotes the entire part of s, or up to order s -1 if s ∈ N. Moreover d s is an ideal metric [START_REF] Carrillo | Contractive probability metrics and asymptotic behavior of dissipative kinetic equations[END_REF]. Its main properties are the following (1) For all probability distributions f 1 , f 2 , f 3 ,

d s (f 1 * f 3 , f 2 * f 3 ) ≤ d s (f 1 , f 2 );
(2) Define for a given nonnegative constant a the dilatation

f a (v) = 1 a f v a .
Then for all probability distributions f 1 , f 2 , and any nonnegative constant a

d s (f 1,a , f 2,a ) = a s d s (f 1 , f 2 ).
The d s -metric is related to other more known metrics of large use in probability theory [START_REF] Gabetta | Metrics for probability distributions and the trend to equilibrium for solutions of the Boltzmann equation[END_REF].

A Kinetic Description of the Rosenau Approximation

Let v ∈ R denote velocity, and let us assume to have at time τ > 0 a system of particles immersed in a background. Let us suppose that the number of particles in the system is sufficiently large to be studied by means of statistical mechanics, namely by giving the velocity distribution f (v, τ ) at time τ > 0. Moreover, let us assume that the main phenomenon which can modify particle's velocity is the interaction of particles with the background. Let M (w) denote the fixed in time (probability) distribution of the particles of the background, which we will assume of finite variance γ 2 . Assume that the collision process of a particle with velocity v with a background particle with a velocity w generates a post-collision velocity v * given by (2.1)

v * = v + w.
Then, in a suitable scaling [START_REF] Cercignani | The Boltzmann Equation and its Applications[END_REF], the effect of interactions (2.1) on the time-variation of the density f (v, τ ) can be quantitatively described by a linear Boltzmann-type equation, in which the variation of the density is due to a balance between a gain and loss terms. More precisely, for a given number v, we take into account all the interactions of type (2.1) which end up with the number v * (gain term) as well as all the interactions which, starting from the number v, lose this value after interaction (loss term). The balance equation for the density of particles can be fruitfully written in weak form. It corresponds to say that the aforementioned interaction process on particles modifies the solution f (v, τ ) according to

(2.2) d dτ R ϕ(v)f (v, τ ) dv = λ R 2 ϕ(v * ) -ϕ(v) f (v, τ )M (w) dv dw,
where the constant λ > 0 denotes the intensity of the variation process, and ϕ(v) is a smooth function. Note that choosing ϕ(v) = 1 shows that, independently of the background distribution, f (v, τ ) remains a probability density if it so initially

R f (v, τ ) dv = R f 0 (v) dv = 1.
This is in general the unique conservation law associated to equation (2.2). From now on, let us assume in addition that the probability distribution of the background is centered, and its variance depends on a small parameter ε > 0. To emphasize this dependence, we will denote this distribution by

M ε (w) = ε -1 M (ε -1 w). Then M ε (w) satisfies (2.3) R M ε (w) dw = 1, R w M ε (w) dw = 0, R w 2 M ε (w) dw = ε 2 γ 2 .
The weak formulation (2.2) yields immediately the time evolution of the moments of f . Taking

ϕ(v) = v one obtains d dτ R vf (v, τ ) dv = λ R w M ε (w) dw R f (v, τ ) dv = 0. Moreover, if ϕ(v) = v 2 (2.4) d dτ R v 2 f (v, τ ) dv = λ R 2 (v * ) 2 -v 2 f (v, τ ) M ε (w) dv dw = λ ε 2 γ 2 .
Thus, the second moment of f grows linearly with respect to time and depends on ε. One way to avoid this dependency is to scale the time properly. Setting t = ε 2 τ and introducing a new particles distribution function g ε such that g ε (v, t) = f (v, τ ) gives according to (2.4),

d dt R v 2 g ε (v, t) dv = 1 ε 2 d dτ R v 2 f (v, τ ) dv = λ γ 2 ,
and the second moment of g ε does not depend on ε.

The distribution g ε is a weak solution to

(2.5) d dt R g ε (v, t)ϕ(v) dv = λ ε 2 R 2 [ϕ(v + w) -ϕ(v)] g ε (v, t )M ε (w) dv dw, with g ε (v, t = 0) = g 0 (v) = f 0 (v). Note that, since R 2 ϕ(v + w)g ε (v, t) M ε (w) dv dw = R 2 ϕ(z) g ε (v, t) M ε (z -v) dv dz = R ϕ(z)M ε * g ε (z) dz, equation (2.5
) can be rewritten as

(2.6) d dt R g ε (v, t) ϕ(v) dv = λ ε 2 R ϕ(v) (M ε * g ε (v) -g ε (v)) dv,
which is the weak form of (1.2). Hence, if the background distribution is given by (1.3), equation (2.5) coincides with Rosenau's approximation (1.5), where σ 2 = λ. Of course, other choices of the background are possible, and, provided conditions (2.3) are satisfied, the evolution of the moments of the solution to (2.5) up to the second order do not depend on the background distribution.

It is then easy to see that, in the case in which M ε coincides with the "Maxwellian" (1.3), g ε is an approximate (at the leading order in ε) weak solution of the heat equation (1.4). Indeed, given that ε is small enough, one can Taylor expand ϕ(v * ), where v * = v + w, to obtain

ϕ(v * ) = ϕ(v) + w ϕ ′ (v) + w 2 2 ϕ ′′ (v) + 1 3! ϕ (3) (ṽ)w 3
for ṽ ∈ (v, v + w). Using this relation in (2.5), one obtains

(2.7) d dt R g ε (v, t)ϕ(v) dv = λγ 2 2 R g ε (v, t)ϕ ′′ (v) dv + R(ε)
where the remainder R(ε) satisfies

(2.8) |R(ε)| ≤ λ 3! ϕ (3) L ∞ 1 ε 2 R w 3 M ε (w) dw = ε λ 3! ϕ (3) L ∞ R w 3 M (w) dw.
Remark 1 . One can easily notice thanks to (2.7) that equation (2.5) is an approximation to the heat equation (1.4) (with diffusion coefficient σ 2 = λ γ 2 /2), provided that the remainder converges to zero as ε → 0. In order to have this convergence, it is enough that the distribution of the background is such that some moment of order greater than two remains bounded. We will use this in the following to obtain various approximations of the heat equation, just changing the distribution of the background.

Remark 2 . The Rosenau type kinetic equation (2.5) is such that mass, momentum and energy of its solution have the same evolution of the corresponding moments of the solution to the heat equation (1.4). A further interesting analogy with the heat equation is given by studying the evolution of convex functionals along the solution. Let Φ(r), r ≥ 0 be a (regular) convex function of r. Then, using equation (2.6) we obtain

d dt R Φ(g ε (v, t)) dv = R Φ ′ (g ε (v, t)) ∂g ε (v, t) ∂t dv = λ ε 2 R Φ ′ (g ε (v, t)) (M ε * g ε (v) -g ε (v)) dv.
Thanks to the convexity of Φ(•), for r, s ≥ 0

Φ ′ (s)(r -s) ≤ Φ(r) -Φ(s),
and one obtains

d dt R Φ(g ε (v, t)) dv ≤ λ ε 2 R (Φ(M ε * g ε (v)) -Φ(g ε (v))) dv.
Now, use the fact that M ε is a probability distribution, so that by Jensen's inequality

R Φ(M ε * g ε (v)) dv = R Φ R g ε (v -w))M ε (w) dw dv ≤ R 2 Φ(g ε (v -w))M ε (w) dv dw = R Φ(g ε (v)) dv,
and one concludes with

d dt R Φ(g ε (v, t)) dv ≤ 0.
Thus any convex functional is non-increasing along the solution to the Rosenau type kinetic equation (1.5), in agreement with the analogous property of the heat equation.

Remark 3 .

A leading example of Rosenau type approximation is obtained by assuming that the background distribution M ε is a balanced Bernoulli distribution

(2.9) M ε (v) = 1 2 [ δ 0 (v + εγ) + δ 0 (v -εγ)] ,
where as usual δ 0 denotes a Dirac mass concentrated at v = 0. Note that the "Maxwellian" (2.9) satisfies (2.3). Then, according to (2.6), the distribution g ε solves (2.10)

∂g ε ∂t (v, t) = λ 2ε 2 [ g ε (v -εγ, t) -2g ε (v, t) + g ε (v + εγ, t)]
. Now, let us fix λ = 2 and γ = σ. Then, given a parameter ∆v > 0 and a uniform grid

v i = i∆v, i ∈ Z on R, by setting g i (t) := g ε (v i , t) and ε = ∆v/σ, equation (2.10) reads g ′ i (t) = σ 2 g i+1 (t) -2g i (t) + g i-1 (t) ∆v 2 .
This is exactly the classical semi-explicit second order central differences scheme for the heat equation (1.4). It follows that this scheme furnishes an approximation to the heat equation which satisfies all the properties outlined in Remarks 1 and 2 , as well as the weak convergence towards the exact solution, which follows easily by applying the estimate (2.8).

Remark 4 . The previous example introduces into the matter a background with a non regular distribution (presence of point masses). This situation clearly differs from the standard Rosenau approximation, characterized by the regular density (1.3). Due to its wide applications in numerical simulations, the study of the properties of the Rosenau approximation given by this type of background turns out to be important. For this reason, in addition to the study of the classical Rosenau model, in what follows we will deal also with the Bernoulli type background, outlining when possible the principal differences.

Representations of the Solutions to the Rosenau Equation

The Rosenau approximation (1.5) of the heat equation is a linear kinetic equation of Boltzmann type. Therefore, we can resort both to linear kinetic theory and to the theory of linear diffusion equations to obtain explicit representations of the solution. We shall present here two equivalent ways to construct these solutions, each one giving a different insight on the behavior of g ε . In the rest of this Section, to avoid inessential difficulties, we will assume that the initial value g 0 (v) is a probability density with finite moments up to a given order (in general more than two).

Wild Sums.

The first method we introduce to obtain an explicit representation of the solution g ε is closely related to the so-called theory of Wild sum expansions of the Boltzmann equation. Wild indeed proved in [START_REF] Wild | On Boltzmann's equation in the kinetic theory of gases[END_REF] that one may represent solutions to the nonlinear Boltzmann equation for Maxwell molecules by using convergent power series. This method is particularly simple when one deals with linear convolution equations such as (1.5). To start with, let us set

h ε (v, t) := exp λ t ε 2 g ε (v, t).
Then the Cauchy problem for (1.5) can be rewritten as a fixed point problem (which is nothing but a Duhamel formula). Let h → Φ ε (h) define the map

(3.1) Φ ε (h) = g 0 + λ ε 2 t 0 M ε * h(s) ds. Then g ε (v, t) solves (1.5) exactly when Φ ε (h ε ) = h ε .
To find the fixed point, it is then sufficient to make a Picard iteration. Starting from h (0)

ε := g 0 one defines, for n ∈ N (3.2) h (n+1) ε = Φ ε h (n) ε .
By recurrence, for n ≥ 1 it yields

(3.3) h (n) ε (v, t) = h (n-1) ε (v, t) + λ t ε 2 n 1 n! M * n ε * g 0 (v),
where we use the shorthand

M * n := M * • • • * M (n times). Clearly, for n ≥ 0 h (n+1) ε (v, t) -h (n) ε (v, t) ≥ 0.
Hence the sequence h

(n) ε (v, t) n≥0 , being bounded in L 1 , converges towards h ε (v, t) ≥ 0 when n → ∞
, where according to (3.3)

h ε (v, t) = g 0 (v) + n≥1 λ t ε 2 n 1 n! M * n ε * g 0 (v).
By passing to the limit n → ∞ in (3.2), one obtains that h ε is a (nonnegative) fixed point for (3.1). This procedure allows to write the (nonnegative) solution g ε to the Cauchy problem for (1.5) as

(3.4) g ε (v, t) = e -λ t/ε 2 g 0 (v) + e -λ t/ε 2 n≥1 λ t ε 2 n 1 n! M * n ε * g 0 (v).
3.2. Fourier Transform. The standard approach to solve the heat equation on the whole space is to use Fourier transform. This is even more obvious for equation (1.5). Indeed, passing to Fourier variables in this equation shows that the Fourier transform of g ε is given by

(3.5) g ε (ξ, t) = g 0 (ξ) exp (-A ε (ξ) t) , ∀(ξ, t) ∈ R × R + ,
where g 0 is the Fourier transform of the initial density and

(3.6) A ε (ξ) := λ 1 -M ε (ξ) ε 2 . Since M ε is a probability density, | M ε (ξ)| ≤ 1 for all ξ ∈ R,
A ε (ξ) := λγ 2 2 ξ 2 + o(ξ 2 ),
where the o(ξ 2 ) term is considered for ξ → 0 (which will be the case until the end of the article). Finally, it is possible to invert the Fourier transform (3.5) and obtain in a distributional sense

(3.7) g ε (v, t) = g 0 * P ε (•, t)(v),
where we set Clearly, representations (3.4) and (3.7) have to coincide. Indeed, applying a Fourier transform to (3.4) yields

P ε (v, t) = exp {-A ε (ξ)t} ∨ (v).
g ε (ξ, t) = e -λ t/ε 2 g 0 (ξ) + e -λ t/ε 2 n≥1 λ t ε 2 n 1 n! M ε n (ξ) g 0 (ξ) = e -λ t/ε 2 e Mε(ξ)λ t/ε 2 -1 + 1 g 0 (ξ) (3.8) = exp (-A ε (ξ) t) g 0 (ξ),
where A ε is given by (3.6).

Remark 5 . According to relation (3.8), any Rosenau-type approximation to the heat equation is such that the smoothness of the solution is lost, independently of the regularity of the kernel M ε . Indeed, the fundamental solution to equation (1.5) (namely the solution to the Cauchy problem obtained from the initial value g 0 = δ 0 ) can be represented in Fourier variable by

G ε (ξ, t) = e -λ t/ε 2 e Mε(ξ)λ t/ε 2 -1 + e -λ t/ε 2 =: G 1 (ξ, t) + G 2 (t). (3.9) If M ε ∈ L 1 , then according to Riemann-Lebesgue Lemma, M ε (ξ) → 0 when |ξ| → ∞ and then |G 1 (ξ, t)| → 0.
The loss of smoothness comes from G 2 (t) which is the Fourier transform of a Dirac mass (of weight depending on time). According to (3.9), the singularity disappears exponentially when ε → 0 or time goes to infinity.

If one considers the

Rosenau kernel (1.3) with λ = σ 2 , M ε (ξ) = 1 1 + (εσ) 2 ξ 2 ,
and for all t, ε > 0,

G 1 (ξ, t) = e -σ 2 t/ε 2 exp σ 2 t ε 2 1 1 + (εσ) 2 ξ 2 -1 .
Thus G 1 (, ξ, t) decreases exponentially fast towards zero when |ξ| → ∞, and the loss of smoothness comes only from G 2 .

On the contrary, if one considers the central differences kernel (2.9) with λ = 2 and σ > 0, then M ε (ξ) = cos(εσξ) and for all t, ε > 0,

G 1 (ξ, t) = e -2t/ε 2 exp 2t ε 2 cos(εσξ) -1 ,
which does not converge at infinity. Then, if the kernel has some singularities, the solution to the Rosenau equation has at most the regularity of the initial condition. This enlightens a substantial difference of the Rosenau approximation with respect to the heat equation, in which the diffusion operator is such that the solution is instantaneously smoothed out. The representation (3.4) can be used to understand the behavior in time of the solution obtained from the central differences kernel. Indeed, each convolution by M ε consists in creating a weighted sum of each point (i.e. Dirac masses) of the grid (see Figure 3.1):

M ε * δ 0 (v) = 1 2 (δ 0 (v -εσ) + δ 0 (v + εσ)) , M * 2 ε * δ 0 (v) = 1 4 δ 0 (v -2 εσ) + 1 2 δ 0 (v) + 1 4 δ 0 (v -2 εσ), M * 3 ε * δ 0 (v) = 1 8 δ 0 (v -3 εσ) + 3 8 δ 0 (v -εσ) + 3 8 δ 0 (v + εσ) + 1 8 δ 0 (v + 3 εσ)
and so on. Thus, the fundamental solution of the central difference scheme is given by

g ε (v, t) = e -2t/ε 2 δ 0 (v) + e -2t/ε 2 n≥1 2t ε 2 n 1 n! g (n) ε (v),
where, for n ∈ N we defined

g (2n) ε (v) = 1 2 n n k=0 2n 2k [δ 0 (v + 2 kεσ) + δ 0 (v -2 kεσ)] , g (2n+1) ε (v) = 1 2 n+1 n k=0 2n + 1 2k + 1 [δ 0 (v + (2k + 1)εσ) + δ 0 (v -(2k + 1)εσ)] .

Asymptotic Behavior of Solutions to the Rosenau Equation

In Section 2 the similarities between the Rosenau-type kinetic equation and the linear diffusion equation have been enlightened. In particular, as well known from the classical literature on numerical approximation of partial differential equations, equation (1.5) provides a consistent approximation of the heat equation in a fixed time interval. It is not known, however, if equation (1.5) still realizes a good approximation for large times, and, in case of a positive answer, in which way the difference between the solutions of the Rosenau and diffusion equations can be estimated with respect both to time and the parameter ε.

In this Section, we will furnish a partial answer to this question. To clarify our intent, we will briefly resume various well-known facts about the large-time behavior of the solution to the heat equation posed in the whole space. Let g = g(v, t) be a solution to the heat equation (1.4) and define the heat kernel Ω σ by

(4.1) Ω σ (v, t) := 1 √ 4πσ 2 t exp - v 2 4σ 2 t , ∀(v, t) ∈ R × R + .
This fundamental solution of (1.4) represents an intermediate asymptotics of a large class of solutions to the heat equation. The recent review article [START_REF] Bartier | Improved intermediate asymptotics for the heat equation[END_REF] gives a precise state of the art on this topic. To make this concept more precise, we define the Boltzmann entropy of f as

H(f ) := R f (v) log f (v) dv.
Then it can be shown (see e.g. [START_REF] Toscani | Kinetic approach to the asymptotic behaviour of the solution to diffusion equations[END_REF]) that g(v, t) behaves as the heat kernel when t → ∞, provided that the initial condition g 0 is of finite kinetic energy and entropy:

g 0 ∈ L 1 2 and H(g 0 ) < ∞.
Moreover, the rate of convergence towards the fundamental solution can be computed in

L 1 norm (4.2) g(t) -Ω σ (t) L 1 ≤ C √ 1 + 2t ,
where C is an explicit constant. The bound (4.2) is sharp. A marked improvement of the constant in (4.2) has been recently obtained in [START_REF] Arnold | Improved entropy decay estimates for the heat equation[END_REF], by selecting well parametrized Gaussian functions, characterized either by mass centering or by fixing the second moments or the covariance matrix of the solution.

Most of these results follow by considering the equivalence between the heat equation and the linear Fokker-Planck equation

(4.3) ∂u ∂t = ∂ ∂v vu + σ 2 ∂u ∂v .
which has a stationary solution given by the Gaussian (4.4)

ω σ (v) := Ω σ (v, t = 1).
Setting T (t) := e 2t-1 /2, we can in fact consider the change of variables (4.5) u(v, t) := e t g e t v, T (t) .

Then u(v, t) is solution to the Fokker-Planck equation (4.3) as soon as g is a solution to the heat equation (1.4). The converse also holds true. It has been shown in [START_REF] Carrillo | Exponential convergence toward equilibrium for homogeneous Fokker-Plancktype equations[END_REF] that the solution to the linear Fokker-Planck equation converges for large times to the stationary solution ω σ and the following asymptotic behavior holds

u(t) -ω σ L 1 ≤ Ce -t , ∀t ≥ 0.
Using the self-similar scaling (4.5), it is clear that this result is equivalent to the rate (4.2), and once more, the bound is sharp. We shall use a similar rescaling in our study of the intermediate asymptotics of equation (1.5).

Our goal in the following is to prove that the same kind of result holds for the Rosenau-type approximation to the heat equation. Due to the generality of this approximation, which depends of the background distribution M , we will obtain a weaker convergence (with respect to the L 1 distance) to the fundamental solution, as discussed in Section 3, Remark 5 . The L 1 distance will be here substituted by the suitable weaker Fourier based metrics d s described in Section 1.2, which are particularly adapted to the convolution structure of the kinetic equation (2.5). This family of metrics has been introduced in the article [START_REF] Gabetta | Metrics for probability distributions and the trend to equilibrium for solutions of the Boltzmann equation[END_REF] to study the trend to equilibrium of solutions to the space homogeneous Boltzmann equation for Maxwell molecules, and subsequently applied to a variety of problems related to kinetic models of Maxwell type. For a more detailed description, we address the interested reader to the recent lecture notes [START_REF] Carrillo | Contractive probability metrics and asymptotic behavior of dissipative kinetic equations[END_REF].

These metrics can be easily and fruitfully applied to the study of the large-time behavior of the heat equation. Indeed, we have seen that in Fourier variables, the solution to the heat equation (1.4) with initial condition g 0 is g(ξ, t) = g 0 (ξ) exp -σ 2 ξ 2 t . Thus, if g 0 is of finite mass, the solution will converge pointwise towards 0 when t → ∞. A way to consider a nontrivial limit distribution is to make the change of variable ξ → V (t) ξ with V (t) := (1 + t) -1/2 . Then, the scaled distribution will converge towards the Gaussian distribution

g 0 (0) exp -σ 2 ξ 2 ,
and one can improve the rate (4.2) in Fourier distance, as soon as the initial datum possesses finite moments of certain order equal to those of the Gaussian [START_REF] Goudon | Fourier-based distances and Berry-Esseen like inequalities for smooth densities[END_REF]. From now on, let us set

(4.6) ω σ (ξ) = exp -σ 2 ξ 2 ,
the Fourier transform of the Gaussian distribution (4.4). We have the following result.

Proposition 4.1. Let V (t) := (1 + t) -1/2 . For a given s > 1, let 0 ≤ g 0 ∈ L 1 s a distribution of unit mass. Then, if g(v, t) is the unique solution to the Cauchy problem (1.4) with initial condition g 0 , the scaled distribution h(v, t) := V (t) -1 g V (t) -1 v, t is such that, for all t ≥ 0 (4.7)

d s (h(t), ω σ ) ≤ 1 (1 + t) s/2 d s (g 0 , ω σ ).
Proof. Thanks to the scaling properties of the Fourier transform one has

d s (h(t), ω σ ) = sup ξ∈R\0 | g(V (t)ξ, t) -ω σ (ξ)| |ξ| s = sup ξ∈R\0 1 |ξ| s g 0 (V (t)ξ) -exp -σ 2 ξ 2 V 2 (t) exp -σ 2 ξ 2 t 1 + t ≤ V (t) s sup ξ∈R\0 1 |ξ| s g 0 (ξ) -exp -σ 2 ξ 2 .
Remark 6 . Taking a Dirac mass as the initial condition in the Cauchy problem for equation 1.4 shows that these rates are optimal. 4.1. Approximate Solutions with Finite Energy. To recover the asymptotic behavior of the solution to the Rosenau-type approximation (2.6), we are going to apply a technique similar to that used in Proposition 4.1. Let us consider again the scaling ξ → V (t)ξ with V (t) := (1 + t) -1/2 , and define h(v, t)

:= V (t) -1 g V (t) -1 v, t and h ε (v, t) := V (t) -1 g ε V (t) -1 v, t .
We choose in this section an initial datum 0

≤ g 0 ∈ L 1 2 such that (4.8) R g 0 (v) dv = 1, R v g 0 (v) dv = 0, R v 2 g 0 (v) dv = E < ∞,
and λ γ 2 /2 = σ 2 . Then the quantity d 2 (g 0 , ω σ ) is bounded, and, according to (4.7) one obtains (4.9)

d 2 (h ε (t), ω σ ) ≤ 1 1 + t d 2 (g 0 , ω σ ) + d 2 (h ε (t), h(t)), ∀t ≥ 0.
Hence, we have to estimate the Fourier distance between the exact scaled solution h and its approximate counterpart h ε . We shall use for this the representation (3.5) of g ε in Fourier variables. Let us introduce a positive parameter R. Since the mass of the initial datum g 0 is equal to 1, one has

d 2 (h ε (t), h(t)) ≤ sup ξ∈R\0 1 ξ 2 exp -A ε ξ √ 1 + t t -exp -σ 2 ξ 2 t 1 + t ≤ 2 R 2 + sup |ξ|≤R 1 ξ 2 exp -A ε ξ √ 1 + t t -exp -σ 2 ξ 2 t 1 + t ≤ 2 R 2 + sup |ξ|≤R 1 ξ 2 exp -σ 2 ξ 2 t 1 + t 1 + t σ 2 ξ 2 A ε ξ √ 1 + t -1 -1 .
Using the elementary inequality |1 -e -x | ≤ |x|, one finally has (4.10)

d 2 (h ε (t), h(t)) ≤ 2 R 2 + sup |ξ|≤R D ε (ξ, t),
where

D ε (ξ, t) := t (εξ) 2 λ 1 -M ε ξ √ 1 + t - (εσξ) 2 1 + t .
In order to compute this last quantity and to obtain the rate of convergence, let us specify both the background distribution and the constants γ and λ.

Central Differences. As in Remark 3 , let us set λ = 2 and

M ε (v) = 1 2 [δ 0 (v + εσ) + δ 0 (v -εσ)].
Then M ε (ξ) = cos(εσξ), and thanks to Taylor's theorem, there exists θ ∈ (0, 1) such that

(4.11) D ε (ξ, t) = t 3 εσξ (1 + t) 3/2 sin θ εσξ √ 1 + t ≤ t (1 + t) 2
(εσξ) 2 3 .

Gathering inequalities (4.10) and (4.11) yields for R > 0 the inequality

d 2 (h ε (t), h(t)) ≤ 2 R 2 + t (1 + t) 2
(εσR) 2 3 .

It just remains to optimize on R to obtain (4.12)

d 2 (h ε (t), h(t)) ≤ 3σ 2 2 ε √ t 1 + t .
Inequality (4.12) shows that the second term in (4.9) has a slower rate of decay with respect to the first one. In contrast to the exact solution h, which is known to converge in d 2 towards the self-similar one at rate (1 + t) -1 [START_REF] Goudon | Fourier-based distances and Berry-Esseen like inequalities for smooth densities[END_REF], we found that the approximate solution converges towards the exact one at rate (1 + t) -1/2 . We just stated the following Proposition 4.2. Let 0 ≤ g 0 ∈ L 1 2 satisfy (4.8). If g ε is solution to the associated Cauchy problem for the Rosenau-type approximation (2.6) with the central difference kernel (2.9), then the scaled distribution

h ε (t, v) := V (t) -1 g ε V (t) -1 v, t for V (t) := (1 + t) -1/2 verifies for all t ≥ 0 d 2 (h ε (t), ω σ ) ≤ 1 1 + t d 2 (g 0 , ω σ ) + 3σ 2 2 ε √ t 1 + t .
Rosenau Regularization. In this case the Fourier transformed "Maxwellian" reads

M ε (ξ) = 1 1 + (εσ) 2 ξ 2 ,
and λ = σ 2 . We obtain

D ε (ξ, t) = (εσ) 2 t 1 + t ξ 2 1 + t + (εσ) 2 ξ 2 .
Inserting this expression in inequality (4.10) gives for R > 0

d 2 (h ε (t), h(t)) ≤ 2 R 2 + t (1 + t) 2
(εσR) 2 3 .

It just remains to optimize on R to obtain the rate of convergence in d 2 norm of h ε towards h:

(4.13) d 2 (h ε (t), h(t)) ≤ σ 2 2 ε √ t 1 + t .
Hence, this approximation has exactly the same order of convergence than the one given by the central differences kernel. We proved

Proposition 4.3. Let 0 ≤ g 0 ∈ L 1 2 satisfy (4.

8). If g ε is solution to the associated Cauchy problem for the Rosenau-type approximation (2.6) with the Rosenau kernel (1.3), then the scaled distribution

h ε (t, v) := V (t) -1 g ε V (t) -1 v, t for V (t) := (1 + t) -1/2 verifies, for all t ≥ 0 d 2 (h ε (t), ω σ ) ≤ 1 1 + t d 2 (g 0 , ω σ ) + σ 2 2 ε √ t 1 + t .
Remark 7 . Note that as a by-product of our analysis of the large-time behavior of the solution to the Rosenau-type approximation of the diffusion equation we proved that the d 2 -distance between the solution of the diffusion equation and its numerical approximation by the central difference scheme is uniformly bounded in time with respect to the parameter ε, as given by (4.12).

Approximate Solution with

Finite Fourth Order Moment. In Subsection 4.1, the convergence results required the boundedness of the d 2 -distance between the initial value of the Rosenau equation and the Gaussian density. This boundedness has been achieved by assuming that the initial datum is of finite energy. Indeed, moments play an important role with respect to convergence. We will improve here the results of Subsection 4.1, provided that the initial condition has finite moments up to order four. We will see that this allows to compute the rate of convergence (in time) towards self-similarity with very few assumptions on the kernel M ε . Without loss of generality (thanks to scaling and translational invariance), we can choose the initial distribution 0 ≤ g 0 ∈ L 1 4 to satisfy

(4.14) R g 0 (v) ϕ(v) dv = (1, 0, 1, 0, µ),
where

ϕ(v) = (1, v, v 2 /2, v 3 , v 4
), and µ is a positive constant. If g is solution to the Cauchy problem for equation (1.4) with initial condition g 0 , then on the one hand

(4.15) R g(v, t) ϕ(v) dv = (1, 0, 1 + tσ 2 , 0, µ + 24t + 12σ 2 t 2 ).
On the other hand, a solution g ε to equation (2.6) with same initial datum verifies

(4.16) R g ε (v, t) ϕ(v) dv = (1, 0, 1 + tσ 2 , 0, µ + (24 + B ε )t + 12σ 2 t 2 )),
provided that M ε follows assumption (2.3) and where

B ε := 2 ε 2 R M ε (v)v 4 dv.
We can see that the fourth order moment of the exact and approximate solutions differ at time t > 0 of a quantity equal to B ε t. Keeping the same notations of Subsection 4.1, we will use this fact to compute the d 3 distance between the scaled approximate solution h ε and the rescaled exact one h. Indeed, it is well know that if k ∈ N, the k th moment of a distribution correspond to the k th derivative of its Fourier transform. Then g(•, t) and g ε (•, t) are at least four times differentiable and according to Taylor's theorem (4.15)-(4.16) imply

g(ξ, t) = 1 + 1 + tσ 2 2 ξ 2 + 1 + 24t + 12σ 2 t 2 24 ξ 4 + O |ξ| 5 , g ε (ξ, t) = 1 + 1 + tσ 2 2 ξ 2 + 1 + (24 + B ε )t + 12σ 2 t 2 24 ξ 4 + O |ξ| 5 .
Thus, one has for R > 0

d 3 (h ε (t), h(t)) = sup ξ∈R\0 1 |ξ| 3 g ε ξ √ 1 + t , t -g ξ √ 1 + t , t ≤ 2 R 3 + B ε t 24(1 + t) 2 R. (4.17)
Optimizing (4.17) over R and gathering the result with inequality (4.7) (with s = 3) gives Theorem 4.1. Let 0 ≤ g 0 ∈ L 1 4 satisfy (4.14). If g ε is solution to the associated Cauchy problem for Rosenau approximation (2.6) with a kernel M ε verifying (2.3), then the scaled distribution

h ε (t, v) := V (t) -1 g ε V (t) -1 v, t for V (t) := (1 + t) -1/2 is such that, for all t ≥ 0 d 3 (h ε (t), ω σ ) ≤ 1 (1 + t) 3/2 d 3 (g 0 , ω σ ) + 13 √ 2 24 (B ε ) 3/4 √ t 1 + t 3/2
, where ω σ is the Gaussian distribution (4.6) and

B ε := 2 ε 2 R M ε (v)v 4 dv.
Remark 8 . The result of Theorem 4.1 indicates in a clear way that the rate of convergence towards the fundamental solution improves as soon as the initial value has a a higher number of moments equal to those of the Gaussian. In the present case, with four moments the difference in Fourier norm between solutions converges to zero at rate (1 + t) -3/4 . This suggests that the rate of convergence approaches the optimal rate when the number of finite moments of the initial value approaches infinity. This aspect of the problem is deeply connected to the analogous one studied in the central limit theorem, which is referred to as Berry-Esseen estimates. We remark that Berry-Esseen-like estimates for Fourier metric based metrics have been studied in [START_REF] Goudon | Fourier-based distances and Berry-Esseen like inequalities for smooth densities[END_REF].

Let us now specify the kernel M ε , to check if there are essential differences among the rate of convergence in ε.

The choice of the central differences kernel

M ε (v) = 1 2 [δ 0 (v + εσ) + δ 0 (v -εσ)]
gives B ε = 2ε 2 σ 4 . Then

d 3 (h ε (t), ω σ ) ≤ 1 (1 + t) 3/2 d 3 (g 0 , ω σ ) + C 1 ε 3/2 √ t 1 + t 3/2
, where C 1 is a positive constant depending on σ. Compared to (4.12), the rate in ε is improved by a power 1/2.

If we consider now the Rosenau kernel

M ε (ξ) = 1 1 + (εσ) 2 ξ 2 , B ε = Kε 3 for a nonnegative constant K. Then d 3 (h ε (t), ω σ ) ≤ 1 (1 + t) 3/2 d 3 (g 0 , ω σ ) + C 2 ε 9/4 √ t 1 + t 3/2
, where C 2 is a positive constant depending on σ. Remark 9 . We can see here that the Rosenau kernel gives a better convergence rate in ε towards the exact solution, compared to the central differences kernel. This improvement also depends of the number of finite moments. This shows that the number of finite moments plays an essential role also in connection with the order of convergence in ε. Indeed in formulas (4.12) and (4.13), which refer to the situation in which only moments up to order two have been considered, this difference between the approximations due to the two kernels was not evident. 4.3. Strong Convergence of a "Regularized" Approximate Solution. Let us assume that the kernel M ε belongs to L 1 (R). The analysis of Section 4 shows that, in a suitable scaling which allows to maintain the energy of the solution bounded, there is convergence in Fourier distance towards a Gaussian function. This result, however, cannot be directly used to conclude that the solution to the Rosenau kinetic equation converges towards the fundamental solution to the heat equation, as time goes to infinity. In fact, since the d s -metric is not scaling invariant, and the decay in time of the distance found in Theorem 4.1 is of order t -3/4 , by reverting to the original variables, the decay in time disappears. Nevertheless, these weak convergence estimates can be fruitfully employed to prove that the solution to Rosenau equation can be split into two parts, one of which is exponentially decaying to zero with respect to both time t and the parameter ε, while the other converges strongly in time towards the fundamental solution to the heat equation at a lower rate. This result clarifies the nature of the Rosenau-type approximation to the heat equation. This approximation produces a solution which essentially consists of two parts, which differs for their regularity. The part with lower regularity (essentially the first term in (3.4) with the same regularity of the initial datum) is rapidly decaying to zero. The regular part is shown to behave like the classical solution to the heat equation, and decays towards the self-similar Gaussian solution for large times. The result will be achieved by making a heavy use of the representations given in Section 3. In particular, we will write the fundamental solution (3.9) by splitting it in two parts

G ε (ξ, t) = e -λ t/ε 2 e Mε(ξ)λ t/ε 2 -1 -M ε (ξ) + 1 -M ε (ξ) e -λ t/ε 2 =: P ε,reg (ξ, t) + 1 -M ε (ξ) e -λ t/ε 2 . (4.18)
As we shall see, this splitting separates in a natural way the singular part of the kernel from its regularized part P ε,reg .

Let g (respectively g ε,reg ) a solution to the heat equation (1.4) (resp. a solution to the Rosenau equation (1.5) obtained by convolution with the regularized kernel). In other words

g(v, t) = g 0 (v) * Ω σ (•, t), g ε,reg (v, t) = g 0 (v) * P ε,reg (•, t). (4.19)
We remark that, resorting to formula (3.4), one recovers that g ε,reg (v, t) coincides with the Wild sum expansion in which the zero-order term is substituted by the regular term M ε * g 0 (v). Hence, g ε,reg (v, t) is nothing but a regularized version of the Wild sum expansion (3.4). Since the initial datum g 0 (v) is of unit mass,

g(t) -g ε,reg (t) L 1 = R R (Ω σ (v -w, t) -P ε,reg (v -w, t)) g 0 (w) dw dv ≤ Ω σ (t) -P ε,reg (t) L 1 . On the other hand, if λ = σ 2 P ε,reg (ξ, t) -Ω σ (ξ, t) = e -σ 2 t/ε 2 e Mε(ξ) σ 2 t/ε 2 -1 -M ε (ξ) -e -σ 2 |ξ| 2 t ≤ e -σ 2 t/ε 2 e Mε(ξ) σ 2 t/ε 2 -e -σ 2 |ξ| 2 t + 1 -M ε (ξ) e -σ 2 t/ε 2 ,
and this implies (4.20)

d 2 (P ε,reg (t), Ω σ (t)) ≤ d 2 (G ε (t), Ω σ (t)) + d 2 (M ε , δ 0 ) e -σ 2 t/ε 2 .
By means of the scaling property of the Fourier metric d s , formula (4.20) ensures that, after the scaling ξ → ξ/ √ 1 + t, the convergence results (in scaled variables) of Section 4 guarantee the convergence of the regularized kernel towards the Gaussian fundamental solution. Now, thanks to the regularity of P ε,reg (t), we can consider also convergence of the (regularized) approximate fundamental solution P ε,reg towards the heat kernel Ω σ in stronger norms, typically L 1 , which are invariant with respect to dilatation. Actually, it is enough to deal with L 2 norms. Indeed, it is well known [START_REF] Carlen | Propagation of Smoothness and the Rate of Exponential Convergence to Equilibrium for a Spatially Homogeneous Maxwellian Gas[END_REF][START_REF] Carrillo | Contractive probability metrics and asymptotic behavior of dissipative kinetic equations[END_REF] that if f ∈ L 2 (R) is of finite kinetic energy, there exists an explicit, nonnegative constant C 1 such that (4.21)

f L 1 ≤ C 1 f 4/5 L 2 R |v| 2 |f (v)| dv 1/5
.

Note that the second moment of the quantities involved is uniformly bounded only in the scaled variables used in Section 4. In original variables, in fact, the second moment of the solution, both of the heat equation and its Rosenau approximation, is increasing linearly with respect to time. Then, it is possible to interpolate the L 2 norm by some appropriate Fourier and homogeneous Sobolev norms. More precisely, it was shown in [START_REF] Carlen | Propagation of Smoothness and the Rate of Exponential Convergence to Equilibrium for a Spatially Homogeneous Maxwellian Gas[END_REF] that for all r ∈ (0, 1) there exists an explicit, nonnegative constant

C 2 = C 2 (r) such that (4.22) f -g L 2 ≤ C 2 d 2 (f, g) 2(1-r) f -g 2r ḢN ,
where N = 3 (1 -r)/r. Then, according to inequalities (4.21)-(4.22)-(4.20) and to the results of Subsection 4.2, we will obtain the strong convergence of P ε,reg (t, •) towards Ω σ (t, •) when t → ∞ if this function belongs to the homogeneous Sobolev space ḢN (R) for arbitrarily small N , and at the same time its Sobolev norm increases in time at a rate such that the right-hand side of inequality (4.22) decays to zero as time goes to infinity.

We will deal with the Rosenau kernel with λ

= σ 2 M ε (ξ) = 1 1 + (εσ) 2 ξ 2
which belongs to Ḣs (R) for all s < 1. Hence, in scaled variables, for all s < 1

e -σ 2 t/ε 2 M ε ξ √ 1 + t 2 Ḣs = e -σ 2 t/ε 2 √ 1 + t εσ 1+s M ε (ξ) 2 Ḣs ≤ C(ε) M ε (ξ) 2
Ḣs . According to (4.18), it remains to check the smoothness of G ε (ξ, t) := e Mε(ξ) σ 2 t/ε 2 -1, ∀t ≥ 0. Let us set s < 1. We have thanks to de l'Hospital rule

lim |ξ|→∞ ξ s+1 G ε (ξ, t) = lim |ξ|→∞ - σ 2 t (s + 1) ε 2 ξ s+2 M ε ′ (ξ) e Mε(ξ) σ 2 t/ε 2 = lim |ξ|→∞ 2 σ 4 t s + 1 ξ s+3 (1 + (εσ) 2 ξ 2 ) 2 exp σ 2 t ε 2 (1 + (εσ) 2 ξ 2 ) = 0.
Then we have proved that ξ s G ε (ξ, t) = o(ξ -1 ) which means according to Riemann criterion that G ε (•, t) ∈ Ḣs (R) for any t ≥ 0. Detailed computations, reported in Appendix, allow to conclude that, for s < 1,

(4.23) e -σ 2 t/ε 2 G ε ξ √ 1 + t Ḣs ≤ C(ε, s) (1 + t) δ
where the explicitly computable constant C(ε, s) depends on ε, s and M ε , and δ > 0 can be taken as small as we want. . On the other hand, the singular part of the fundamental solution G ε (ξ, t) converges to zero exponentially both in time and with respect to ε. Hence, after a transitory time, there is prevalence of the regular part of the fundamental solution, which decays to zero in L 1 at a suboptimal rate. Hence, the main difference between the large-time behavior of the heat equation and its Rosenau approximation consists in a slower rate of convergence of the latter towards the Gaussian fundamental solution. The method of proof also indicates that the rate of convergence in L 1 is strictly linked to the regularity of the distribution of the background. For this reason, the argument leading to Theorem 4.2 can not be applied to the central difference approximation. Maybe a different approach will help to clarify if, in some weaker norm, a result similar to that of Theorem 4.2 holds also for the central difference type approximation.

Conclusions

We studied in this article the validity of the approximation to the linear diffusion equation proposed by Rosenau [START_REF] Rosenau | Tempered diffusion: A transport process with propagating fronts and inertial delay[END_REF] as a regularized version of the Chapman-Enskog expansion of hydrodynamics. This approximation essentially is realized by substituting the heat equation with a linear kinetic equation of Boltzmann type, describing collisions of particles with a fixed background. This remark allows to consider the Rosenau approximation as a particular realization of a model Boltzmann equation, in which the background distribution is a general probability density with bounded variance. In addition to the Rosenau distribution, we considered in this article also a point masses background, which furnishes the central difference scheme to solve the linear diffusion equation. The main differences between the action of the two different backgrounds have been studied in some details. In particular, our analysis put into evidence that the approximation with a regular kernel is a good approximation, in that it possesses most of the typical properties of the heat equation, including the same large-time behavior, apart from a slower rate of decay towards the fundamental solution and a finite speed of propagation at high Fourier frequencies. We were not able to prove an analogous property for the non regular approximation, for which the rate of decay towards the equilibrium distribution has been only proven in scaled variables. The results of the present paper shall lead to new numerical approximations of the whole Rosenau approximation (1.1). For instance, Wild sums have been successfully used in connection with the approximation of the bilinear Boltzmann equation (see e.g. [START_REF] Pareschi | Time relaxed Monte Carlo methods for the Boltzmann equation[END_REF]) to develop a new family of Monte Carlo methods, which possess high order accuracy in time, and are asymptotic preserving (with respect to ε). Following the same strategy, our next goal is to develop an efficient Monte Carlo method to compute the Rosenau approximation using our analysis on the linear Wild sums. Since δ is arbitrarily small, inequality (4.23) follows.

Figure 3 . 1 .

 31 Figure 3.1. Schematic representation of the convolution steps of the central difference kernel. Each node is a Dirac mass, weighted by binomial coefficients.

2 dξ 1 / 2 .

 212 Appendix A. Growth in Time of the Sobolev Norm of the Regularized Kernel For any given s < 1, we investigate the growth in time of the quantityB s (t) = e -σ 2 t/ε 2 G ε ξ √ 1 + t Ḣs .To avoid inessential heavy notations, in what follows we fix ε = σ = 1. In this caseB s (t) = e -t R |ξ| 2s exp t 1 + (1 + t) -1 ξ 2 -1By a simple change of variable into the integral, we can rewrite B s (t) in the equivalent form (A.1) B s (t) = (1 + t) s+1/2 I s (t)

  Theorem 4.2 has various consequences. Coupling the decay of the d 2 metric obtained in (4.12) with inequalities (4.21) and (4.22), the regularity of P ε,reg implies a rate of decay of the L 1norm smaller or equal to (1+ t)1/15

	Thanks to inequalities (4.20)-(4.21)-(4.22)-(4.23) and Theorem 4.1 we obtain
	the following
	Theorem 4.2. Under the assumptions of Theorem 4.1, if g is a solution to the heat equation (1.4)

and g ε,reg is given by (4.19) for the Rosenau kernel (1.3), then one has lim t→∞ g(t) -g ε,reg (t) L 1 = 0. Remark 10 .

.

  1/2 = (1 + t) s+1/2The function I s (t) defined in (A.1) satisfies the differential equation -exp(-t) .If now |ξ| > (1 + t) -α we have A s (ξ, t) dξ ≤ t exp -t (1 + t) -2α /2 |ξ| 2s (1 + ξ 2 ) 2 dξ ≤ c s t exp -t (1 + t) -2α /2 , (A.4)where we denoted by c s the bounded constantc s = (1 + t) 2s+1 + c s t exp -t (1 + t) -2α /2 .Hence, considering that the second term in (A.5) decay faster at infinity than the first one, for any given constant a > 0 there exists a bounded constant C a such that if t > a(A.6) Ās (t) ≤ C a t -(2s+1-δ) . inequality (A.6) into (A.2), we obtain, for a given a > 0, and t > a I s (t) ≤ I s (a) e -2t + e -2t -(2s+1-δ) e 2τ dτ ≤ t -(2s+1-δ) , I s (t) ≤ I s (a)e -2t + C a t -(2s+1-δ) . B s (t) = (1 + t) s+1/2 I s (t) 1/2 ≤ I s (a)(1 + t) s+1/2 e -t + (1 + t) δ/2 C a .

	we conclude that, for t > a
	(A.2) where τ which implies d dt I s (t) = -2I s (t) + 2 A s (ξ, t) = |ξ| 2s 1 + ξ 2 exp -t ξ 2 1 + ξ 2 exp -t R A s (ξ, t) dξ, t e -2t a ξ 2 (A.7) 1 + ξ 2 For a given positive constant δ < 2p + 1, let us set Using (A.7) into (A.1) gives
	(A.8)							α = 1 -	2p + 1 δ
	Then we find
	(A.3)	0	(1+t) -α	A s (ξ, t) dξ ≤	0	(1+t) -α	|ξ| 2s dξ =	1 2s + 1	(1 + t) δ (1 + t) 2s+1 .
				ξ 2 1 + ξ 2 ≥	(1 + t) -2α 1 + (1 + t) -2α =	1 1 + (1 + t) 2α ≥	1 2	(1 + t) -2α ,
	and
		exp -t	ξ 2 1 + ξ 2 -exp(-t) =	exp -t ξ 2 1+ξ 2 -exp(-t) -t ξ 2 1+ξ 2 + t	t 1 + ξ 2 ≤	t 1 + ξ 2 .
	The previous inequalities imply
		+∞						+∞
		(1+t) -α						(1+t) -α
										R	|ξ| 2s (1 + ξ 2 ) 2 dξ.
	Putting together inequalities (A.3) and (A.4) we conclude with the upper bound
	(A.5) (1 + t) Substituting t Ās (t) := ξ∈R A s (ξ, t) dξ ≤ 1 2s + 1 Ās (τ ) e 2τ dτ
										a
					≤ I s (a) e -2t + C a e -2t	a t	τ -(2s+1-δ) e 2τ dτ.
	On the other hand, integration by parts shows that
	t a Choosing a in such a way that τ -(2s+1-δ) e 2τ dτ ≤	e 2t 2	   R t -(2s+1-δ) + |ξ| 2s exp -t 2s + 1 -δ 2a 1 + ξ 2 -exp(-t) t τ -(2s+1-δ) e 2τ dτ. 2 a ξ 2 2s + 1 -δ 2a ≤ 1 2 ,	dξ	  	1/2

δ

Acknowledgment

This work was initiated during a visit of Thomas Rey, which thanks the University of Pavia for its grateful hospitality. The research of TR was granted by the European Research Council ERC Starting Grant 2009, project 239983-NuSiKiMo. GT acknowledges support by MIUR project "Optimal mass transportation, geometrical and functional inequalities with applications".


Léa Fanchon

Florent Jacquemard

FORMAL TIMING ANALYSIS OF MIXED MUSIC SCORES

Interactive music systems coordinate in real-time an artificial perception of dynamics of human musicians with timely execution of machine reactions. As every human performance will differ from another, it is a challenging task to be able to predict the behavior or such systems in response to any possible performance, and prevent unwanted outcomes.

We present here the application of formal models and methods from the real-time systems verification literature to the static analysis of interactive music systems.We consider in particular the good parameters problem, which consists in synthesizing a set of timing parameter valuations (representing performances here) guarantying a good behavior of the system analyzed. The methods presented have been applied to the system Antescofo, and are general enough to apply to other interactive music systems.

INTRODUCTION

There are two dimensions in authoring time and interaction for interactive music systems, corresponding to two main processes [START_REF] Rowe | Interactive Music Systems: Machine Listening and Composing[END_REF]. The first process is the composition, or programming i.e. the specification of mixed scores containing both instrumental and electronic parts and instructions for their synchronization. The second is the performance: The real-time evaluation and synchronization of the parts. We consider here the case of the score follower Antescofo [START_REF] Cont | Antescofo: Anticipatory Synchronization and Control of Interactive Parameters in Computer Music[END_REF] which aims at addressing these two processes, enabling live interactive performances of written music between computers and human musicians. It relies on two subsystems: on the one hand, a listening machine [START_REF] Cont | A Coupled Duration-Focused Architecture for Realtime Music to Score Alignment[END_REF] tracks in real-time both the position of a musician in the score and his current tempo, and on the other hand a reactive engine, relying on a dedicated synchronous language, ensures the real-time coordination of the electronic part.

Such a system can be qualified as time critical. Like in other embedded systems e.g. for transportation or communication, time is not only a measure of efficiency but a real semantical issue: The objective is not to compute as fast as possible but to produce some actions at the right moment [START_REF] Lee | Computing needs time[END_REF]. Moreover, such a system is also reactive and open: during a performance, it has to react in realtime to triggering events (notes) coming from an unpredictable environment (the musicians). Indeed, whereas the duration of every note is specified precisely in a music score, it is well known that two human interpretations of the same piece can differ significantly, regarding timing values. Such a behavior is all fine artistically speaking and is considered as a part of the definition of musical performance [START_REF] Manoury | Considérations (toujours actuelles) sur l'état de la musique en temps réel. Etincelle[END_REF]. But wrt computer science unfortunately, it is called unpredictable and must be dealt with. For example, some slight changes in the delays of the instrumental events may change the curse of electronic actions, e.g. by swapping two actions whose ordering ought to be preserved (see Example 3 below). In this context, it is a particularly challenging task to be able to predict before a concert how an interactive system system such as Antescofo with react as expected to any human performance. Even intensive rehearsing cannot cover all cases and computer assistance is required for a complete timing analysis.

As part of the effort that has been made for the development of reliable embedded systems, several verification approaches have been developed in the past years for verifying statically that a system satisfies qualitative and quantitative temporal properties [START_REF] Baier | Principles of Model Checking[END_REF][START_REF] Berard | Systems and Software Verification: Model-Checking Techniques and Tools[END_REF]. Timed automata (TA) [START_REF] Alur | A theory of Timed Automata[END_REF] are a central abstract model of real-time systems in these approaches. They are finite state automata extended with clocks. These clocks are real-time variables increasing uniformly and which can be compared to fixed timing delays in guards of the automata's transitions. Moreover, several TA can be composed in a network (NTA), each component automaton modeling a building block of a complex system. Synchronization of automata is then performed by rendezvous on the symbols of transitions.

In this paper, we present the application of TA based formal methods to the timing analysis of mixed scores in a language close to the one used in Antescofo (presented in the next section). The similarity between augmented scores and NTA is clear: each part of a mixed score can be represented by a TA, and triggering events can be used for synchronization by rendezvous. A procedure of compilation of scores into NTA has been implemented and is presented in Section 3.

Uncertainty on the timing delays in a human performance can be modeled using parametric delays (i.e. unknown constants) in the transitions' guards of the TA modeling the instrumental score. This corresponds to the model of parametric timed automata (PTA) [START_REF] Alur | Parametric Real-Time Reasoning[END_REF]. In Section 4, we show how to apply to the analysis of mixed scores some techniques of parameter synthesis for PTA originally developed for the verification telecommunication protocols and asynchronous circuits, e.g. [START_REF] André | An Inverse Method for Parametric Timed Automata[END_REF]. The good parameters problem, consists, given a PTA A, in inferring a set of parameters values for which A behaves well [START_REF] Frehse | A Counterexample-Guided Approach to Parameter Synthesis for Linear Hybrid Automata[END_REF]. Such a set is presented by a linear constraint on parameters.

Finally, we propose in Section 5 another approach which can be used in a computer assisted composition process in other to find appropriate timing values for the delays associated to electronic actions. Assuming that the delays for instrumental events have been fixed, and that a constraint expresses the expected ordering of events and actions, a valuation of the actions' delays is chosen in a way that it maximizes a measure of the robustness of the mixed score to the timing variations in a musician's performance.

MIXED SCORES SPECIFICATION LANGUAGE

Antescofo uses for the specification of mixed scores a textual language which is similar to synchronous languages for programming real-time embedded systems such as Esterel [START_REF] Berry | The Esterel synchronous programming language and its mathematical semantics[END_REF]. Programs in this language describe both instrumental and electronic parts of the score and the instructions for their coordination during a performance. We present in the following grammar the syntax of a simplified fragment of this language, highlighting aspects relevant to the following automated verification techniques.

Mixed Score Syntax

A score is defined as an interleaving of instrumental events to be recognized by the listening machine and of the actions or group of actions triggered by the former. Antescofo can handle several kinds of events: notes, chords, trills, glissandi etc which are all abstracted here in symbols denoted e 0 , e 1 ,... in a an assumed predefined set E. Actions correspond to specific commands sent to the audio environment (a sound synthesis module or real-time audio processing system such as Max/MSP [START_REF] Puckette | Combining event and signal processing in the MAX graphical programming environment[END_REF] or Pure Data [START_REF] Puckette | Pure Data[END_REF]), and are represented as symbols a 0 ,... in a set A. We consider also a set L for group labels ( in the grammar). Every event (resp. action or group) is associated with a duration (resp. a preceding delay) d specifying a number of beats (relative time). This value can be converted into physical time (i.e. in seconds) using the tempo detected in real-time by the listening machine as explained in introduction. For now, let us assume rational valued delays d ∈ Q + . The extension to delay expressions with variables is discussed later.

Idealized Performance

Atomic actions can be sequentially composed into groups which can be compared to staves in traditional scores. Intuitively, a group g is triggered by the last event e occurring before g in the score S. If g occurs immediately after e in S, and is associated a delay d (denoted (d g) in the e 1 1.0 0.0 group g 1 (0.5 init) 0.25 group g 3 (0.5 msg) (0.5 off) e 2 1.0 0.0 group g 2 (0.5 on) e 2 1.0 score) then after the detection of e, the system will wait for d time units (relatively to the detected tempo) before launching the group g, which consist in playing the sequence of elements contained in the group in the specified order. If g is followed by a group g with delay d in the score, then after launching g, the system waits d time units and launch g concurrently. Groups can be nested, and the above interpretation of successive elements in the score is transposed recursively to nested sequences of elements in groups.

There are several ways to define the termination of loops in Antescofo language. Here, we assume for simplifying that a loop is associated to a terminating event or action κ. The loop's body will be iterated until it is killed, as soon as κ is received. The semantics of the execution of groups and loops will be made formal in Section 3.4. Let us illustrate this first in the following small example.

Example 1. We consider as a running example a trivial score containing 3 events e 1 , e 2 , e 3 and 3 groups g 1 , g 2 and g 3 (Figure 1). The groups g 1 and g 3 are triggered by the event e 1 and g 2 by e 2 . Like g 1 , g 2 is launched without delay, and contains one atomic action (an initialization message "init" for the former and a command for switching the lights on for the second). The group g 3 is launched concurrently to g 1 after a delay of 0.25 time unit, and contains two elements: a message "msg" and a command for switching lights off.

This behavior is illustrated in Figure 2, assuming the musician plays the events e 1 , e 2 and e 3 with delays conforming strictly to those specified in the score of Figure 1: δ 1 = δ 2 = 1.0 time units (we call this an idealized performance below). In this case, the lights are first switched off at time unit 1.25 (by group g 3 ), and switched back on 0.25 time units later (by g 2 ).

Group Attributes

In the language of Antescofo, the groups can be assigned attributes defining specific strategies of synchronization with events or for error handling, etc see [START_REF] Cont | Correct Automatic Accompaniment Despite Machine Listening or Human Errors in Antescofo[END_REF]. For the sake of conciseness, we shall shall consider here one option of synchronization strategy: if a group has attribute "tight", each of its elements is automatically realigned to the latest event scheduled before it (instead of the latest event occurring before it in the score, like for groups without this attribute).

Example 2. The score of Figure 4 is almost 1 equivalent semantically to the one of Example 1: the "tight" group g is played as it was split in the two groups g 1 and g 2 triggered respectively by e 1 and e 2 . Note that the delay of 0.5 associated to the action "on" in group g 2 corresponds to the remaining after the occurrence of e 2 of the original delay 1.0 associated to "on" in group g.

Using a "tight" attribute permits to express that a large group is implicitly split into smaller groups synchronized to the nearest detected events, which is useful e.g. in a context of automatic accompaniment, where the electronic (accompaniment) part has to follow closely the leading instrumental voice. This feature has been added to Antescofo in order to prevent the composer from the burden of having to manually segment the accompaniment part in small groups, keeping a high-level view of the grouping structure of the score [START_REF] Cont | Correct Automatic Accompaniment Despite Machine Listening or Human Errors in Antescofo[END_REF].

Musical Performances

By definition, in a musical performance (as opposed to the idealized performance presented in Example 1), the delays between events should always diverge from those specified in the score. This may lead to unexpected behaviors for the execution of actions, like e.g. unwanted inversions as illustrated in the next example.

Example 3. Let us come back to Example 1. If the musician plays e 2 earlier than expected, e.g. with δ 1 = 0.7 time units after e 1 in Figure 3, then the lights are first 1 The difference is explained in Example 6.

switched on (with no effect) and switched off 0.05 time units after, and the concerts finishes in the dark.

One could fix this by putting the actions "init" and "on" in a same group g triggered by e 1 , i.e. by the group of Example 2 without the attribute tight, or alternatively by giving an attribute tight to the group g 3 .

The purpose of the rest of the paper is to present formal methods for automatically analyze all possible behaviors of the system during musical performances. We use for this purpose the model of timed automata, presented in the next section, for a formal representation of mixed scores, and then apply methods based on this formalism.

TIMED AUTOMATA MODEL

Parametric Timed Automata

A timed automaton [START_REF] Alur | A theory of Timed Automata[END_REF] is a finite automaton extended with a finite set of real-valued variables called clocks. We assume given a set X of clocks denoted x, y,... Each clock can be independently reset to 0 on a transition and keeps track of the elapsed time since the last reset. Every transition of the automaton is guarded by a constraint on the clock values: the transition may be taken only if the current values of the clocks satisfy the associated constraint. Note that all the clocks evolve at the same rate. This corresponds exactly to the behavior of Antescofo, where all groups are synchronized on the tempo detected by the listening machine.

A timed automaton recognizes timed words: sequences of symbols associated with a timestamp (dates are expressed as real values). This formalism is well adapted to the expression of the timing constraints appearing in mixed scores. It is indeed a powerful model for describing both the logical ordering of the events and the timing between the events.

In order to be able to handle all possible musical performances, we are going to reason parametrically: the event's delays in a score are replaced by parameters δ 1 , δ 2 ,... from a given set ∆, as in Figures 2 and3. The instantiation of these parameters by real values defines then a performance. The appropriate formalism for modeling this are parametric timed automata (PTA) [START_REF] Alur | Parametric Real-Time Reasoning[END_REF].

Before stating the formal definition of PTA, let us describe the constraints that we shall use. A constraint over the clocks and the parameters is a conjunction of linear inequalities of the form t < t or t ≤ t , where t and t are linear terms over X and ∆:

i a i x i + i b i δ i + c with x i ∈ X, δ i ∈ ∆, and a i , b i , c ∈ Q + . The set of constraints over X and ∆ is denoted C(X, ∆). A PTA is a tuple A = L, 0 , F, Σ, inv , T
where L is a finite set of locations, 0 ∈ L is the initial location, F ⊆ L is the set of final (recognizing) locations, Σ is a finite alphabet of symbols, inv : L → C(X, ∆) assigns an invariant to each location (condition on clocks to be allowed to stay in this location), and T ⊆ L × C(X, ∆) × Σ × 2 X × L is a finite set of transitions. The transition , g, a, R, from the location to the location is guarded by the constraint g on clocks and parameters, and specifies an observed symbol a and the subset of clocks R ⊆ X to be reset when the transition is taken. Such a transition will also be denoted g,a,R ---→ . Examples of PTA can be found in Section 3.4. PTA without parameters correspond to standard timed automata (TA) [START_REF] Alur | A theory of Timed Automata[END_REF]. Given a PTA A and a valuation π : ∆ → Q + defined on the parameters of A, the TA A[π] is defined from A by replacing every parameter δ by its value π(δ).

Timed Traces and Languages

A run of a TA A = (L, 0 , F, Σ, inv , T ) is a finite sequence of the form

0 , v 0 d1,a1 ---→ 1 , v 1 d2,a2 ---→ . . . d k ,a k ---→ k , v k where for every 1 ≤ i ≤ k, i is a location in L, v i : X → Q + is a clock valuation, d i ∈ Q + is a delay and every a i ∈ Σ. Given d ∈ Q + , let us denote v i + d the clock valuation x → v i (x) + d.
The conditions on the above run are that v 0 satisfies inv ( 0 ), and for all

1 ≤ i ≤ k, • for all 0 ≤ d i ≤ d i , v i-1 + d i satisfies inv ( i-1 ),
• there exists a transition i-1 gi,ai,Ri -----→ i in T such that:

• v i-1 + d i satisfies g i , • for all x ∈ X, v i (x) = 0 if x ∈ R i and v i (x) = v i-1 (x) + d i otherwise, • v i satisfies inv ( i ).
The run is recognizing if k ∈ F , and v 0 (x) = 0 for all x ∈ X. In this case, we say that the timed word a 1 , t 1 , . . . , a k , t k ∈ (Σ × Q + ) * is recognized, where t i = i j=1 d i for all 1 ≤ i ≤ k. Its projection a 1 , . . . , a k over Σ * is called the associated untimed word. The (un)timed language of A is the set of recognized (un)timed words.

Automata Networks and Synchronization

It is convenient to use a parallel composition operator for PTA when modeling a system made of several subsystems for verification. This approach is particularly relevant for modeling polyphonic mixed scores, where voices (in our case the instrumental part and each of the electronic groups) evolve concurrently and each of them can be represented by a PTA.

A network of parametric timed automata (NPTA) is a set of PTA, denoted with the operator of parallel composition. The automata of the network are supposed to have disjoint sets of clocks and locations, but they may share alphabet's symbols which are used for synchronization. The synchronized product of a NPTA A 0

• • • A m , with A i = L i , 0,i , F i , Σ i , inv i , T i for each 0 ≤ i ≤ m, is the PTA m i=0 L i , 0,0 , . . . , 0,m , m i=0 F i , m i=0 Σ i , inv , T
where inv 0 , . . . , m = m i=0 inv i ( i ) and T is the set of all transitions 0 , . . . , m g,a,R ---→ 0 , . . . , m such that for all 0 ≤ i ≤ m, if a ∈ Σ i , then there exists g i and R i such that i gi,a,Ri ----→ i ∈ T i , and otherwise i = i , and g = a∈Σi

g i , R = a∈Σi R i .
The notions defined above for PTA (language, [π]...) are extended to NPTA using this synchronized product.

Compilation of Scores into Automata

The translation of the simplified mixed scores of Section 2 into NPTA is quite straightforward. Roughly, it returns one PTA representing the instrumental part of the score and one PTA for each group.

The events and actions become the symbols in transitions, and the guards are used to enforce the delays specified in the score. Instead of giving a complete definition of the translation process, we illustrate it on Example 4. A detailed translation procedure is described in [START_REF] Echeveste | Antescofo: A Domain Specific Language for Realtime Musician-Computer Interaction[END_REF] for a more complete fragment of the Antescofo language, in order to describe operational semantics for this language. Each of the three automata A 1 , A 2 , A 3 for electronic parts starts with a unguarded transition labeled by the event triggering the group. According to the semantics of the synchronized product of this NPTA presented in Section 3.3, this permits to synchronize the start of the groups with the transitions of the instrumental automaton. Let us consider for instance the PTA A 1 associated to the first group g 1 . It has three locations

g 0 1 , g 1 1 , g 2 1 .
Starting in the initial location g 0 1 , it waits for the detection of e 1 . Note that the guard of the transition between g 0 1 and g 1 1 is empty (i.e. equal to true): there are no condition on clock values to fire this transition. Once e 1 is detected, A 1 enters g 1 1 and waits for 0.5 time units (as expressed by the guard y 1 ≤ 0.5), sends the message "init", before entering the final state g 2 1 . The translation of a group g n nested in a group g follows the same principle, using the group label n of g n , instead of an event e i , as triggering symbol. In the outer group g, a transition labelled with n will trigger the nested group g n . The nested group g n is compiled in a separate PTA A n , which starts with a transition labelled with n and without guard, exactly like A 1 , A 2 , A 3 in Ex. [START_REF] André | Behavioral Cartography of Timed Automata[END_REF].

For a loop, we use just the same construction as for a group, except for an additional unguarded backward transition from the final location to the first one (actually the target location just after the triggering transition). Moreover, for stopping the loop when the terminating event or action κ is received, we add a transition labelled with κ from every non-final state to the final state.

Example 5. Figure 6 presents the PTA corresponding to the following loop, with triggering event e 1 and terminating event or action κ. The construction for groups with an attribute "tight" is illustrated in the next example. Example 6. The translation of the score of Example 2 (Figure 4) is presented in Figure 7. The group g with attribute "tight" is compiled into one PTA with 5 locations g 0 to g 4 . Note the transition labeled with e 2 between the transitions labeled with "init" and "on". The behavior of this PTA is be the following: after the detection of event e 1 , it enters location g 1 , waits for 0.5 time units, sends the message init, enters g 2 and waits for the detection of the event e 2 . After the detection of e 2 , it enters g 3 . However, if e 2 arrives before the expiration of the delay of 0.5, then init is skip and the PTA jumps to g 3 (transition with the guard y < 0.5). This is the difference between the "tight" group g of Figure 4 and g 1 , g 2 of Figure 1 (other options are possible in Antescofo [START_REF] Cont | Correct Automatic Accompaniment Despite Machine Listening or Human Errors in Antescofo[END_REF][START_REF] Echeveste | Antescofo: A Domain Specific Language for Realtime Musician-Computer Interaction[END_REF]). Once in g 3 , the PTA waits for 0.5 time units and send the message "on".

Example 7. The following is a run of the TA A 1 of Example 4 (Figure 5): Let A be the synchronized product of N and let π 0 = {δ 1 → 1.0, δ 2 → 1.0}. The following run corresponds to the timed trace of A[π 0 ] depicted in Figure 2 (with an additional delay of 0.1 before e 1 ). ----→ e 3 , g 2 1 , g 2 2 , g 3 3 , x=y2=1.0,y1=2.0,y3=1.25 . Let π = {δ 1 → 0.7, δ 2 → 0.9}. The following run corresponds to the timed trace of A[π] depicted in Figure 3. ----→ e 3 , g 2 1 , g 2 2 , g 3 3 , x=y2=0.9,y1=1.6,y3=0.85 . The NPTA obtained from a mixed score is the basis of verification procedures presented in the next section.

g 0 1 , y 1 = 0 0.1,e1 ----→ g 1 1 , y 1 = 0 0.5,init ----→ g 2 1 , y 1 = 0.5 g 0 g 1 y≤0.5
e 0 , g 0 1 , g 0 2 , g 0 3 , x=y1=y2=y3=0 0.1,e1 ----→ e 1 , g 1 1 , g 0 2 , g 1 3 , x=y1=y3=0,y2=0.1 0.5,init ----→ e 1 , g
e 0 , g 0 1 , g 0 2 , g 0 3 , x=y1=y2=y3=0 0.1,e1 ----→ e 1 , g 1 1 , g 0 2 , g 1 3 , x=y1=y3=0,y2=0.1 0.5,init ----→ e 1 , g

VERIFICATION OF MIXED SCORES

Given a NPTA model N of a mixed score, we will show how to infer a linear constraint K on the parameters of N guaranteeing an acceptable behavior of the automata network (referred as the good parameters problem in [START_REF] Frehse | A Counterexample-Guided Approach to Parameter Synthesis for Linear Hybrid Automata[END_REF]). In other words, K is such that if a valuation π : ∆ → Q + defined on the parameters of N satisfies K, then all the timed words recognized by A[π] belongs to an acceptable set (where A is the synchronized product of N ). The constraint K hence represents a range of human performances π for which the system is guaranteed to behave as expected, from a discrete point of view. We shall present below two notions of acceptable sets of behavior, and define first the crux tool for inferring the expected constraint K: Symbolic semantics of PTA.

Symbolic States and Traces

A constrained PTA is a pair, denoted as A(K), made of a PTA A and a constraint over parameters K ∈ C(∆).

The symbolic semantics of a constrained PTA A(K) with p clocks is a transition system over Σ denoted by A ς . It has a set of states of the form , C , with ∈ L and C ∈ C(X, ∆), initial state 0 , C 0 with

C 0 = K ∧ inv ( 0 ) ∧ p-1 i=1 x i+1 = x i ,
(x ) = ∃x, d C(x) ∧ g(x) ∧ xi∈R x i = 0 ∧ xi / ∈R x i = x i ∧ inv ( )(x ) ∧ inv ( )(x + d) where x = x 1 , . . . ,

Depth First Symbolic State Exploration

In this first approach, we assume given a finite state automaton B over Σ and say that a timed word w recognized by the NPTA N is acceptable if its projection on Σ * is recognized by B. Hence, B represents the acceptable discrete behavior of N , restricting e.g. the ordering on which some actions should be played. The automaton B can be seen as a special case of TA whose guards and invariants are all set to true. Therefore, it can be composed with N , to restrict (i.e. control) the language of the network.

The procedure for the inference of a constraint K as expected consists in a depth-first exploration of the symbolic states reachable by A ς , where A is the synchronized product of N B. The synchronized product is not computed a priori but on the fly, at each step of the exploration. Note that the tree T of partial symbolic traces of A ς is finite in the case of mixed scores without loops, or with loops terminated by an event or action occurring in the score. In this case, the depth of T is the length of the instrumental score. Its size can be exponential though. In order to avoid state space explosion, on can prune a search branch as soon as the constraint of the current symbolic state is unsatisfiable. It is correct since the sequence of symbolic constraints in a symbolic trace is monotonic. The expected constraint K is then defined as the disjunction of the constraints on the leaves of T whose state contains a final location. The untimed word e 1 , e 2 , init, msg, off leads to an unsatisfiable constraint (hence this branch is pruned). In-

c 0 c 1 c 2 Σ 0 off Σ 0 on Σ 0 Figure 8. Control automaton B (Example 8, § 4.2).
deed, the constraint reached after e 1 , e 2 contains δ 1 ≤ 0.5, because of the invariant y 1 ≤ 0.5 associated to the location g 1 1 in A 1 . Furthermore, it also contains 1.25 -δ 1 ≤ 0.5, because of the invariant y 2 ≤ 0.5 for location g 1 2 and the guards in the transition of A 3 with msg and off. Hence, we have a contradiction.

On the other hand, the branch corresponding to the idealized untimed word e 1 , init, msg, e 2 , off, on leads to the constraint 0.75 ≤ δ 1 ≤ 1.25 ∧ 0.5 ≤ δ 2 . This word can be obtained only for these values of parameters.

Despite the pruning, the state space remain large, and this naive approach (not implemented) should be tractable only for small excerpts of scores. It makes sense however for problems like the inversion of "off" and "on" in Ex. 1: The verification can be narrowed to a window delimited by the triggering events (e 1 and e 2 in the example) and the end of the groups (g 3 and g 2 ) containing these actions.

Inverse Method

A more efficient approach is to focus on the original score for defining the acceptable behavior. This corresponds to the purpose of the inverse method (IM), implemented in the system IMITATOR [START_REF] André | An Inverse Method for Parametric Timed Automata[END_REF]: Given with a NPTA N (with synchronized product A), it starts with a good parameter instantiation π 0 , and generalizes it. More precisely, it synthesizes a constraint K 0 ∈ C(∆) such that π 0 satisfies K 0 and for all parameter instantiation π satisfying K 0 , the untimed languages of A[π] and A[π 0 ] coincide. In other words, IM can be applied when the acceptable set contains the timed words whose projection are in the untimed language of A[π 0 ], i.e. the (unique) untimed word corresponding to an idealized performance in our case.

The algorithm IM can be described roughly as follows. Starting with K := true, it iteratively computes a set S of reachable states of A(K)

ς . When a π 0 -incompatible state , C is found in S (i.e. such that π 0 does not satisfy ∃x C(x)) then K is refined as follows: A π 0 -incompatible inequality on parameters J is selected in (the normal form of) ∃x C(x) and ¬J is added to K. This way, in the next iteration, A(K) ς will avoid the counterexample. The algorithm stops when the state set S reaches a fix point, and returns the intersection of the constraints of the states in the final set S. Note there are 2 non-deterministic choices in IM: the selection of the π 0 -incompatible state , C and of J, and in general IM is non-confluent: several application may lead to different results.

Example 9. When run on the NPTA of Example 4, Figure 5, IMITATOR reaches a fixpoint after 6 iterations and returns the constraint K = 0.75 < δ 1 ≤ 1.0 ∧ 1.25 < δ 1 + δ 2 , which defines a smaller domain than the naive approach. Moreover, the application of the algorithm is deterministic in this case.

Discussion

Reachability is undecidable for PTA hence in general, the algorithms presented in the two previous sections can diverge. The problem however becomes decidable on NPTA associated to scores without loops. IMITATOR has given good results for the verification of case studies including circuits and communication protocols [START_REF] André | Behavioral Cartography of Timed Automata[END_REF]. The associated models typically contained few (less than a dozen) parameters and locations. As we associate one parameter to every instrumental event in a mixed score, our models can be much larger, in term of number of parameters and locations. In counterpart, they have also a simpler structure. A traditional bottleneck for the verification of automata models equipped with clocks or registers is the occurrences of some control locations in more than one simple cycle [START_REF] Comon | Multiple Counters Automata, Safety Analysis and Presburger Arithmetic[END_REF]. This situation is generally avoided in mixed scores and the experiment conducted so far have given satisfying results. Another originality of our approach is that the PTA models are not constructed manually but obtained by compilation of mixed scores.

In the literature on the verification of embedded or cyber-physical systems, the uncontrollable environment is generally related to physical phenomena: stabilization time in circuits, temperature in a nuclear core, linear acceleration and angular velocity returned by sensors in a plane, GPS... The case of interactive music systems brings a human in the computation loop, possibly opening new application perspectives in this research field, see [START_REF] Cont | From authored to produced time in computer-musician interactions[END_REF].

COMPOSITION ASSISTANCE

We briefly describe in this last section another application of the parametric approach presented above. The objective is to provide assistance to a composer by proposing appropriate delays for some actions, according to a given ordering constraint on actions and events. More precisely, let us consider an unfinished mixed score S where all the delays of the instrument events e 1 , . . . , e n are set values d 1 , . . . , d n , but some actions' delays are still parameters α 1 , . . . , α m ∈ ∆ (we note α = {α 1 , . . . , α m }). Moreover, we assume that the expected behavior of S is specified by a linear constraint K ∈ C(α ∪ δ), where δ = {δ 1 , . . . , δ n } ⊂ ∆ is a set of parameters associated to the events' delays. We show below that we can propose a valuation π of α satisfying K[δ i := d i , 1 ≤ i ≤ n] and maximizing the robustness of S to the timing variations in the instrumental performance (K[π] denotes the replacement in K of parameters by their value defined by π).

Definition of Robustness. For 1 ≤ i ≤ n, let σ i be the valuation defined on δ\{δ i } by σ i (δ j ) = d j for all j = i, and for π defined on α, let inf [π] i and sup[π] i be the lower and upper bounds for δ i defined by K

[π ∪ σ i ]. We let rob S (K, π, e i ) = min(|d i -inf [π] i |, |d i -sup[π] i |)
and rob S (K, π) = min 1≤i≤n rob S (K, π, e i ) . Intuitively, the smaller this robustness value is, the closer to the delays d i the musician must play in order to satisfy K.

Total Ordering Constraint. When K is defined by a total ordering < on the events and actions of S (i.e. the expected untimed word), then a valuation π of α maximizing the above robustness can be defined, after some variable renaming, from the following instantiationof the delays between consecutive actions wrt < (see [START_REF] Fanchon | Temporal Analysis of Mixed Intrumental/Electronic Music Scores[END_REF]). Let 1 ≤ i ≤ n, and assume e i < a 1 i < . . . < a ni i < e i+1 , where there are no actions or events between successive elements in this sequence (wrt <).

For all 1 ≤ j < n i , if the actions a j i and a j+1 i are triggered by the same event, then the delay between them is set to the smallest delay supported by the system (let us call ε this value). Otherwise, the delay between a j i and a j+1 i is set to di+1-pi×ε ni+1-ci , where p i is the number of pairs of successive actions in a = {a 1 i < . . . < a ni i } triggered by the same event, and c i is the number of sets of successive actions in a triggered by the same event. Moreover, the delay between e i and a 1 i and between a ni i and e j+1 is also set to the same value. Intuitively, successive actions triggered by the same event must be as close as possible and other successive actions must be as far as possible.

Regular Constraint. The constraint K can also defined from a set of acceptable untimed words, e.g. presented by a FSA B as in Section 4.2. In this case, we can pick an untimed word w recognized by B and propose one optimal π for the total ordering defined by w, using the above technique. Finding one word w recognized by B such that the associated π maximizes robustness is an open question.

CONCLUSION

In general, the purpose of verification procedures such as model checking [START_REF] Baier | Principles of Model Checking[END_REF] is to ensure (e.g. with a formal proof destined to a certification authority), that a critical embedded system satisfies some correctness property in every possible execution. The approaches presented in Section 4 for the timing analysis of mixed scores aim rather at characterizing precisely the range of instrumental executions (performances) for which the reactive engine of Antescofo will behave as expected. This seems more appropriate in the context of interactive music systems. To our knowledge, it is the first time that such methods are applied to of interactive music systems. The outcome of the verification procedures can be exploited in several ways, for instance for warning performers about sensible parts of the score, or to assist the composer in adjusting the delays of electronic actions in these parts. The CAC approach proposed in Section 5 follows the same idea. A user friendly presentation of the outcome (currently linear constraints), directly in the score, is under study.

As possible extensions, one could use (linear) expressions over parameters for expressing delays in transitions' guards of PTA (instead of a single parameter), e.g. for specifying mixed scores with variables [START_REF] Vercoe | Real-time Csound: Software Synthesis with Sensing and Control[END_REF][START_REF] Mccartney | Supercollider: a new real time synthesis language[END_REF]. The generalization of IM to this case could be interesting. Another application target could be non linear scores (aka open scores), in which the instrumental part contains branches, the direction in these branches being chosen by the performer at real-time. Our framework could be extended to this case with a straightforward translation of the branches as non-deterministic transitions in the instrumental PTA (automaton A 0 in Example 4).

Other applications of the TA model of Section 3 could also be worth investigating. For instance reachability analysis or timed model checking [START_REF] Larsen | UPPAAL in a Nutshell[END_REF] over the unparametric TA model could be used for verifying temporal properties of the idealized performance. We are also planing to use approaches based on timed games theory [START_REF] Bouyer | Robust Reachability in Timed Automata: A Game-based Approach[END_REF] for giving a quantitative estimate of the robustness (upper bound on the time deviation guarantying expected behaviour) of a mixed score to timing variations in performances.

  score := (event d) | (event d) score | (d comp) score comp := action | group | loop group := group {list} loop := loop {list} until (event | action) list := comp | comp list

Figure 1 .

 1 Figure 1. A small mixed score.

Figure 2 .

 2 Figure 2. Ideal performance for the score of Fig. 1.

Figure 3 .

 3 Figure 3. Unexpected performance for the score of Fig. 1.

Figure 4 .

 4 Figure 4. Score equivalent to Fig. 1, with a tight group.

Example 4 .

 4 The NPTA N = A 0 A 1 A 2 A 3 associated to the mixed score of Example 1 is presented in Figure 5. It has an alphabet Σ = {e 1 , e 2 , e 3 , init, msg, on, off}. Note that the guards involve equalities constraints only. All the bounds of guards in the instrumental automaton A 0 are parameters, and they are real values for the other PTAs A 1 , A 2 , A 3 (which are hence all standard TA).

Figure 5 .Figure 6 .

 56 Figure 5. Automata network for Example 1, Fig. 1.

0. 0

 0 loop (0.1 a) (0.2 b) until κ Note the guards y 2 < 0.1 and y 2 < 0.2 in the transitions labeled with the terminating κ, aiming at the stop state.

Figure 7 .

 7 Figure 7. PTA for a tight group (Ex. 2, Fig.4).

  and transitions , C a -→ , C such that there exists a transition g,a,R ---→ of A and C

Example 8 .

 8 Let us come back to the NPTA N described in Example 4, and let Σ 0 = Σ \ {on, off}. The automaton B of Figure 8 recognizes the words in which both "off" and "on" occur exactly once, and in this order. The complete exploration of the symbolic traces of (N B) ς is presented in Appendix. Let us see two examples of traces.

δ 1 +δ 1 +e 1 init e 2 e 3 e 1 init e 2 msg e 3 δ 2 -( 6 )e 1 2 - 1 .25 + δ 1 0 δ 2 ≤ 0. 5 ( 7 )e 1 ( 8 )e 1 5 e 1 init msg e 2 e 3 δ 2 0 δ 2 5 e 1 init msg e 2 e 3 off δ 2 -δ 1 +e 1 init msg e 2 e 3 off on 0. 5 -2δ 2 +δ 1 -e 1 init msg e 2 off 1 .e 1 init msg e 2 off e 3 δ 2 - 1 .25 + δ 1 0 δ 2 ≤ 0. 5 ( 10 ) e 1 init msg e 2 off e 3 on 0. 5 -δ 2 0 e 1 init msg e 2 off on δ 1 -e 1 init msg e 2 off on e 3 δ 2 -Figure 9 .

 11326121257181525215112125105129 Figure 9. Symbolic traces (Example 8). Delay is from previous parent in the trace tree.

  x p , and x = x 1 , . . . , x p . The formula C (x ) can be transformed into a constraint of C(X, ∆) using for instance Fourier-Motztkin elimination of existentially quantified variables. It expresses the condition in the definition of runs in Section 3.2, applied to a TA A[π] such that the parameter valuation π satisfies K. A symbolic trace is a sequence of symbolic states corresponding to a run of A ς .

Acknowledgements.

The authors wish to thank Romain Soulat and Laurent Fribourg for valuable discussions and their precious assistance on using the tool IMITATOR.

Appendix

We present in Figure 9 a depth-first-search of the trace tree for the NPTA of Example 4, Figure 5 (score of Example 1, Fig. 1).

In this exploration, we give for each (partial) trace the constraint of the symbolic state reached, decomposed in four columns corresponding to the four clocks of the automaton: x, y 1 , y 2 and y 3 . Note that the constraints reached at the different prefix of a trace are accumulated (one should consider a conjunction of all the constraints reached at prefixes). The symbol indicates that the constraint reached is unsatisfiable and that a backtrack is necessary. The column "delay" presents the delay elapsed for performing the last step of the trace.

(1 -4) these partial traces lead to unsatisfiable constraints containing both δ 1 ≤ 0.5 and δ 1 ≥ 0.75.

(5 -6) these traces are possible only if δ 1 = 0.75, δ 2 = 0 (i.e. e 2 and e 3 are simultaneous) and off and on are simultaneous. Hence it is neither a realistic nor a desirable solution. (11) this corresponds to the ideal trace, defined by δ 1 = δ 2 = 1.0. This trace is possible under the constraints 0.75 ≤ δ 1 ≤ 1.25, 0.5 ≤ δ 2 and 1.25 ≤ δ 1 + δ 2 .