Walid Belkhir
email: walid.belkhir@inria.fr

Alain Giorgetti
email: alain.giorgetti@femto-st.fr

Michel Lenczner
email: michel.lenczner@utbm.fr

A Symbolic Transformation Language and its Application to a Multiscale Method

Keywords: Symbolic transformation, term rewriting, strategies, multiscale modeling

The context of this work is the design of a software, called MEMSALab, dedicated to the automatic derivation of multiscale models of arrays of micro-and nanosystems. In this domain a model is a partial differential equation. Multiscale methods approximate it by another partial differential equation which can be numerically simulated in a reasonable time. The challenge consists in taking into account a wide range of geometries combining thin and periodic structures with the possibility of multiple nested scales.

In this paper we present a transformation language that will make the development of MEM-SALab more feasible. It is proposed as a Maple TM package for rule-based programming, rewriting strategies and their combination with standard Maple TM code. We illustrate the practical interest of this language by using it to encode two examples of multiscale derivations, namely the two-scale limit of the derivative operator and the two-scale model of the stationary heat equation.

Introduction

The context of this work is the design of microsystem array architectures, including microcantilevers, micromirrors, droplet ejectors, micromembranes, microresistors, etc., to cite only a few. A model for such arrays is a Partial Differential Equation (PDE). The numerical simulation of whole arrays based on classical methods like the Finite Element Method (FEM) is prohibitive for today's computers (at least in a time compatible with the time scale of a designer). The calculation of a reasonably complex cell of a threedimensional microsystem requires at least 1000 degrees of freedom which lead to at least 10 000 000 degrees of freedom for a 100 × 100 array. Fortunately there is a solution consisting in approximating the model by a multiscale method. The resulting approximated model is again a PDE. It can be rigorously derived from the exact one through a sequence of mathematical transformations, but these transformations differ for each case.

We are currently developing a software, called MEMSALab, for "MEMS Arrays Laboratory", dedicated to multiscale and multiphysics modeling of arrays of micro-and nanosystems. Unlike traditional software that is based on models built once and for all, MEMSALab is a software that constructs models. The challenge consists in taking into account a wide range of geometries combining thin and periodic structures with the possibility of multiple nested scales. One should also consider PDEs representing multiphysics systems with high contrast in equation coefficients.

Simulation software available in the market offers specialized tools for large arrays of micro-and nanosystems, but the construction of new models raises many problems. Firstly the time required for a new design varies from some weeks for a specialist to several months for a beginner. Secondly the mathematical machinery is too sophisticated to be manually applied to complex systems. Finally the resulting models require specific numerical simulation methods, that have to be implemented case by case.

The software MEMSALab we design aims at addressing these problems. It is based on multiscale models, especially on those derived by asymptotic methods. Such asymptotic models are derived from a system of PDEs when taking into account that at least one parameter is very small, such as thickness for a thin structure or the small ratio of a cell size to the global size for a periodic structure. The resulting models are other systems of PDEs, obtained by taking the mathematical limits of the nominal models, in a well-suited sense, when the small parameters tend toward zero. This approach provides a reasonably good approximation. It also offers the advantages and factors of reliability to be rigorous and systematic. The resulting PDEs can be implemented in a simulation software such as the finite element based simulator COMSOL (Multiphysics Finite Element Analysis Software, official site http://www.comsol.com), and simulations turn to be fast as needed.

The literature in this field is vast and a large number of techniques have been developed for a large variety of geometric features and physical phenomena. However, none of them have been implemented in a systematical approach to render it available to engineers as a design tool. In fact, each published paper focus on a special case regarding geometry or physics, and very few works are considering a general picture. By contrast our software will treat the problem of systematic implementation of asymptotic methods by implementing the construction of models rather than the models themselves. This approach will cover many situations from a small number of bricks. It combines mathematical and computer science tools. The mathematical tool is the two-scale transform originally introduced in [START_REF] Lenczner | Homogénéisation d'un circuit électrique[END_REF][START_REF] Casado-Díaz | Two-scale convergence for nonlinear Dirichlet problems in perforated domains[END_REF][START_REF] Cioranescu | Periodic unfolding and homogenization[END_REF] to model periodic and thin structures, and also referred as the unfolding method. We have extended its domain of application to cover in the same time homogenization of periodic media, see for instance [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF], and methods of asymptotic analysis for thin domains, see Ciarlet (1988). The computer science tools include term rewriting, λ-calculus, and type systems (Cirstea and Kirchner, 2001;[START_REF] Marin | Rule-based programming with mathematica[END_REF]Cirstea et al., 2001;[START_REF] Geuvers | Introduction to type theory[END_REF]. The software is written in the symbolic computation language Maple TM .

Compared to other techniques, our multiscale method requires more modular calculations, avoids any non-constructive proof and intensively relies on equational reasoning. The classical way to automate equational reasoning is to consider mathematical equalities as rewrite rules. The rewrite rule t → u orients the equality t = u from left to right and states that every occurrence of an instance of t can be replaced with the corresponding instance of u. Consequently symbolic computation with equalities is reduced to a series of term rewritings. Algebraic computation and term rewriting are two research domains with strong similarities. Both are separately well-studied but there are only few works about the combination of algebraic computation and term rewriting [START_REF] Fèvre | Combining algebraic computing and term-rewriting for geometry theorem proving[END_REF][START_REF] Bündgen | Combining computer algebra and rule based reasoning[END_REF]. Term rewriting provides a theoretical and computational framework which is very useful to express, study and analyze a wide range of complex systems. It is characterized by a repeated transformation of data objects such as words, terms or graphs. Transformations are described by a combination of rules which specify how to transform an object into another one in the presence of a specific pattern. Rules can have further conditions and can be combined by specifying strategies. The latter control the order and the way the rules are applied. Term rewriting is used in semantics in order to describe the meaning of programming languages as well as in program transformations. It is used to perform symbolic computations like in Mathematica, and also to perform automated reasoning. It is central in systems where the notion of rule is explicit such as expert systems, algebraic specifications, etc.

The computer algebra system Maple TM is widely used in the symbolic computation community. It is also used by members of our project for a prototypal implementation of MEMSALab. Maple TM is a suitable language for combining function-based and rulebased symbolic transformations. Unfortunately, it is only equipped with a limited rewrite kernel, namely the function applyrule(rule,expr). The main drawback of applyrule is that it iterates the application of the rule everywhere in the given expression until no matching sub-expression remains. Therefore there is a lack of flexibility: the user cannot express how and where rules must be applied. The other problem is that the Maple TM matching function patmatch(expr,pattern,sub) has two main drawbacks. Firstly, when the user defines new operators, the matching must be done modulo the properties of these operators, e.g. associativity and commutativity. Actually, patmatch does not provide this feature. Secondly the matching function computes just one solution (i.e. a substitution) of the matching problem of pattern with expr, whereas both supporting associative-commutative operators and implementing conditional rewrite rules require a matching that returns all the possible solutions.

Contributions

In this work we present a transformation language named symbtrans (for "symbolic transformation"). It extends Maple TM with conditional rewriting, strategies, and patternmatching modulo associativity and commutativity. All the ingredients of term rewriting with strategies are made explicit. In particular terms, patterns, rules, strategies, pattern matching and application of rules and strategies to terms are represented by Maple TM expressions. Rewrite rules and rewriting strategies are deterministic functions that raise an exception when they are not applicable. Such functions can be combined with Maple TM functions. We illustrate the practical interest of this transformation language by using it to write two formal derivations for the MEMSALab software. The first one states the weak two-scale limit of the derivative operator. The second one computes the two-scale model of the stationary heat equation.

Related work

Our transformation language is an adaptation for Maple TM of popular strategy languages such as ρ-log [START_REF] Marin | Rule-based programming with mathematica[END_REF] or Tom [START_REF] Balland | Tom: Piggybacking rewriting on Java[END_REF]). A conceptual difference with Tom is that the latter extends a host language with an additive syntax, whereas our transformation language smoothly integrates with standard Maple TM functions. The first work that have considered term rewriting from a functional point of view is Elan, see [START_REF] Borovansky | Rewriting with strategies in ELAN: a functional semantics[END_REF], within a non-deterministic framework. Our transformation language is comparable with the deterministic fragment of the rewriting calculus (Cirstea and Kirchner, 2001).

Paper outline

Section 2 introduces the two-scale transform and necessary term rewriting concepts and notations. Section 3 defines the transformation language symbtrans in a detailed way and illustrates transformations by several examples. Section 4 presents more advanced features of symbtrans: the combination of term rewriting with procedural programming, a delayed procedure evaluation mechanism, pattern-matching modulo associativity and commutativity, and conditional rewriting. Section 5 shows the transformation language at work on two realistic examples of multiscale derivations. Section 6 presents the theoretical bases of symbtrans and compares the present work with related ones. Section 7 draws conclusions.

A mathematical proof of the weak two-scale limit property of the derivative operator is reproduced in Appendix A. Its formal proof with symbtrans is reproduced in Appendix C. A mathematical derivation of the two-scale model of the stationary heat equation is reproduced in Appendix B. Its formal counterpart with symbtrans is given in Appendix D. Rules and transformations corresponding to mathematical properties for these formal proofs are reproduced in Appendix E.

Preliminaries and notations

Term, substitution, matching and rewriting

Term rewriting systems (in the classical sense) are defined by specifying a set of terms and a set of rewrite rules. Rewrite rules are applied to reduce terms. In general the process of reduction continues until no more rules can be applied, or forever in the case of non-terminating systems. A term which cannot be reduced to another term is called a normal form. If there is always a unique normal form then the system is said to be confluent.

Let F be a countable set of function symbols, each symbol having a fixed arity. Let X be a countable set of variables. The set of terms, denoted by T (F, X ), is inductively defined as the smallest set containing the elements of X and f (t 1 , . . . , t n ), for any symbol f with arity n in F and any terms t 1 , . . . , t n in T (F, X ). When n = 0 the symbol f is called a constant and the corresponding term is denoted f instead of f ().

The set of variables occuring in a term t is denoted by Var(t). If Var(t) = ∅ then t is said to be ground. A substitution is a function σ from X to T (F, X ) such that σ(x) = x only for finitely many variables x in X . If x 1 , . . . , x n are these variables, then σ is denoted by {x 1 → t 1 , . . . , x n → t n }, where t i = σ(x i ) (1 ≤ i ≤ n). If t is a term then σ(t) is the term that results from the application of σ to t.

A rewrite rule is a pair (l, r) of terms l and r in T (F, X ) s.t. Var(r) ⊆ Var(l). This pair is usually written l → r. Definition 1. For two terms t and t ′ in T (F, X ) the problem of finding substitutions σ such that σ(t ′ ) = t is called a matching problem and is denoted t ′ ≪ t. A substitution σ such that σ(t ′ ) = t is called a solution of this matching problem. The application of a rewrite rule l → r to a term t, denoted by [l → r](t), is defined by σ(r) where σ is any solution of the matching problem l ≪ t. It is undefined when this matching problem has no solution. It is worth mentioning that there is a difference between the notions of variables and constants in the mathematical sense and in the sense of term rewriting when mathematical expressions are viewed as terms. Let us give a concrete example. Consider the mathematical expression Ω f (x)dx and its representation by the term Integral(Omega, f(x), x). In this term, Omega and the mathematical variable x are viewed as function symbols of arity 0 (i.e. constants), f is viewed as a function symbol of arity 1, and Integral as a function symbol of arity 3. Further clarifications on this topic can be found in Section 3.

Running example

We illustrate the transformation language with the homogenization problem of the stationary heat equation. We consider a stationary distribution of temperature in a region Ω ⊂ R n (where n ∈ {1, 2, 3}) with an internal heat source and with imposed vanishing temperature along the boundary. The diffusion coefficient a ε : R n → R is assumed to be periodic on Ω in the n directions, with a small period ε. In other words, there is a function a : R n → R which is (0, 1) n -periodic and such that a ε (x) = a(x/ε) for x ∈ Ω. With a view to derive the so-called homogenized model, the parameter ε is considered as small, and we are interested in finding an approximation of the stationary heat equation when ε decreases to zero. In this mathematical asymptotic process, the distributed internal heat source f ε : R n → R can be considered as depending on ε, and the temperature distribution u ε : R n → R, vanishing on the boundary ∂Ω of Ω, is the unique solution to the stationary heat equation

n i=1 Ω a ε ∂ xi u ε ∂ xi v dx = Ω f ε v dx (1)
written in its variational form, as explained in [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF]. Here, v : Ω → R is any test function, i.e. a sufficiently regular function vanishing on ∂Ω.

To conduct the asymptotic process ε → 0 while keeping as much information as possible in the solution u ε at the small scale, the region Ω is unfolded into the Cartesian product Ω × Y with Y = (0, 1) n through the so-called two-scale transform or unfolding, where Ω ⊂ R n is called the macroscopic domain associated to Ω, and Y is called the microscopic domain associated to Ω. This sort of change of variables is applied to the sequence u ε : Ω → R which yields another sequence of functions T u ε : Ω × Y → R. We can show that the latter converges to a limit u 0 (x, y), so we say that the sequence u ε is two-scale convergent to u 0 . Similarly, we show that a ε , f ε and ∂ xi u ε are two-scale convergent towards some limits a 0 , f 0 and ∂ xi u 0 + ∂ yi u1 , and in the same time that u 0 is independent of the so-called microscopic variable y, it is vanishing on ∂ Ω and u 1 is Y -periodic, a concept that is explained later. Then, applying the two-scale transform in the variational formulation (1) of the heat equation, passing to the limit ε → 0, we arrive to the two-scale model satisfied by the pair (u 0 , u 1 ),

n i=1 Ω×Y a 0 (∂ xi u 0 + ∂ yi u 1 ) (∂ xi v 0 + ∂ yi v 1 ) dxdy = Ω×Y f 0 v 0 dxdy, (2)
where the pair (v 0 (x), v 1 (x, y)) is sufficiently regular, v 0 is vanishing on ∂ Ω and v 1 is Y -periodic. The last step consists in showing that u 1 is a function of u 0 which allows its elimination and yields the homogenized model

1 n j,k=1 Ω a H jk ∂ xj u 0 ∂ x k v 0 dx = Ω f H v 0 dx (3)
for all test functions v 0 . Here a H is the n× n matrix of effective diffusion coefficients, and f H is the effective heat source. We observe that in this example the macroscopic domain Ω is identical to Ω. However, we prefer to distinguish them for sake of generalization, e.g. to thin domains.

The statement of the approximate model (3) derived with the two-scale transform was announced in [START_REF] Lenczner | Homogénéisation d'un circuit électrique[END_REF]. Then several proofs have been published, in [START_REF] Lenczner | Homogenization of electrical networks including voltage-to-voltage amplifiers[END_REF], [START_REF] Casado-Díaz | Two-scale convergence for nonlinear Dirichlet problems in perforated domains[END_REF], [START_REF] Cioranescu | Periodic unfolding and homogenization[END_REF], [START_REF] Lenczner | A two-scale model for an array of AFM's cantilever in the static case[END_REF], and [START_REF] Cioranescu | The periodic unfolding method in homogenization[END_REF]. Here we follow the proof in [START_REF] Lenczner | A two-scale model for an array of AFM's cantilever in the static case[END_REF] where an effort has been made to formulate proofs in an algebraic way and to avoid any abstract reasoning in the sequences of formal transformations. The steps in this formal method are rigorously specified at a high level of generality, that make them independent of the domain geometry and applicable to other equations.

Before entering in more details, we introduce some notations and definitions. For any region Θ ⊂ R n , L 2 (Θ) denotes the set of square integrable functions on Θ, that is the set of functions

v : Θ → R with a bounded L 2 (Θ)-norm ||v|| L 2 (Θ) = ( Θ v 2 (x) dx) 1/2 . Then a sequence u ε ∈ L 2 (Θ) is said to be convergent (or strongly convergent) in L 2 (Θ) towards a limit u 0 if ||u ε -u 0 || L 2 (Θ) = O(ε)
where O(ε) is the Landau notation representing any sequence of ε tending to zero when ε → 0. Another concept of convergence used for asymptotic models is the concept of weak convergence. A sequence

u ε ∈ L 2 (Θ) is said to be weakly convergent in L 2 (Θ) towards a limit u 0 if Θ (u ε -u 0 )v dx = O(ε) for any v ∈ L 2 (Θ).
A by product of this definition is that a sequence u ε ∈ L 2 (Ω) is said to be two-scale weakly convergent towards a limit

u 0 ∈ L 2 ( Ω × Y ) if T u ε is weakly convergent towards u 0 in L 2 ( Ω × Y ). Finally, a function v : Y → R is said to be Y -periodic if v(y + ) = v(y -)
for any pair of opposite points y + and y -of the boundary of Y .

Now we specify the assumptions needed to build the model (3) as the limit of ( 1). The heat source f ε is assumed to be uniformly bounded in the L 2 (Ω)-norm, that is, there exists a constant C independent of ε such that

||f ε || L 2 (Ω) ≤ C. (4)
Then, the proof is divided into five parts. (i) We establish that the fields of temperature distribution u ε and its derivatives are uniformly bounded in the L 2 (Ω)-norm,

||u ε || L 2 (Ω) ≤ C and ||∂ xi u ε || L 2 (Ω) ≤ C for i = 1, ..., n. (5)
Given these results, it is assumed that u ε has an asymptotic expansion on the form

T u ε = u 0 + ε u 1 + εO(ε), (6)
where here O(ε) denotes a function that tends to zero weakly in L 2 ( Ω × Y ). We observe that this assumption is not necessary to get the desired result, however, we find that on the one hand it is not very strong once the a priori estimates are known and on the other hand it allows the proof to be entirely computational, i.e. without steps of abstract reasoning. (ii) The next step consists in deducing of ( 5) that the two-scale weak limit u 0 of T u ε is independent of y ∂ yi u 0 = 0 for all i.

(7) (iii) It comes to show that there exists u 1 (x, y) such that

∂ xi u ε is weakly two-scale convergent towards ∂ xi u 0 + ∂ yi u 1 in L 2 ( Ω × Y ), (8)
together with the relation between u 1 and (u 0 , u 1 ),

u 1 = u 1 + n j=1 y j ∂ xj u 0 (9)
as well as the fact that u 1 is Y -periodic. (iv) Once the two-scale limit of ∂ xi u ε has been computed, it is used in the variational formulation (1) in a manner that yields the two-scale model ( 2). (v) Finally, the homogenized model ( 3) is deduced by expressing each microscopic derivatives ∂ yi u 1 as a linear function of the macroscopic derivatives (∂ xj u 0 ) j=1,...,n , the linear operator between them being solution of a partial differential equation in the cell Y . It would take too long to present in detail all of the above proof. We chose to illustrate the transformation rules only on the third and fourth steps because they require the implementation of most of the useful concepts for the complete proof. In addition, we do not show the periodicity of u 1 nor the relation (9) between u 1 and (u 0 , u 1 ). Thus, the two examples discussed are to prove the following propositions. The mathematical proofs are reported in Appendices A and B respectively when their counterpart formalized by rewriting rules and strategies are in Appendices C and D.

Proposition 3. For a sequence of functions

u ε : Ω → R such that u ε and ∂ xi u ε are bounded in the L 2 (Ω)-norm, if T u ε has a formal expansion of the form T u ε = u 0 + εu 1 + ε n j=1 y j ∂ xj u 0 + εO(ε), ( 10
)
where u 0 is independent of y, the partial function y → u 1 (., y) is Y -periodic and O(ε) tends to zero in L 2 ( Ω × Y ) weak, then (8) holds.

Proposition 4. Assuming that the data of the stationary heat equation ( 1) satisfy

T a ε = a 0 and T f ε = f 0 + O(ε), (11)
where O(ε) tends to zero in L 2 ( Ω × Y ) and the solution u ε satisfies the assumptions and the conclusion of Proposition 3, then the pair (u 0 , u 1 ) is solution to the two-scale model [START_REF]LHSterm := Outermost(ApproximationT2Context)(LHSterm); LHSterm := Outermost[END_REF].

A number of mathematical tools are required to carry out the proofs. Some of them are refered in the body of the paper so they are described in the end of this section while the others are used only in the detailed proofs and thus are presented in the beginning of Appendix A.

In order to formalize the convergence concepts of sequences of numbers and of functions, and to handle them within a computational framework, we recall the usual set of computation rules of the Landau notation O(ε) and we also add redundant rules for sums and integrals. As usual the computation rules on the Landau notation apply from left to right only.

O(ε) + O(ε) = O(ε), -O(ε) = O(ε), O(ε) = O(ε), O(ε) dx = O(ε), O(ε) * O(ε) = O(ε), α * O(ε) = O(ε) if α is independent of ε, and ε = O(ε). (12)
This small system of "axioms" defines what we call the convergence calculus for the rewriting rules. We observe that it does not include a mean for distinguishing the various types of convergences, so for now that distinction is left to the user. We shall also repeatedly use the additional property of O(ε) whereby

Ω g O 1 (ε) dx = O(ε) (13)
as long as g is a function uniformly bounded in L2 (Ω) and lim

ε→0 ||O 1 (ε)|| L 2 (Ω) = 0. It results from the Cauchy-Scwharz inequality Ω g O 1 (ε) dx ≤ ||g|| L 2 (Ω) ||O 1 (ε)|| L 2 (Ω) ≤ C O(ε) = O(ε).
The usual operations on variational formulation require to use the extension of the rule of integration by parts to multidimensional domains. The so-called Green formula holds for sufficiently regular functions u and v defined in a domain Θ

⊂ R n , Θ u ∂ xi v dx = ∂Θ u v (n x ) i ds(x) - Θ v∂ xi u dx (14)
for any i = 1, ..., n where n x represents the outward pointing unit normal of the hypersurface volume element ds(x) on the boundary ∂Θ.

Finally, functions such as u and v are usually considered as elements of a vector space such as L 2 (Ω). As such, it is possible to define linear operators that apply to them. We recall that an operator L defined on a vector space is said to be linear if for any vectors v, w and any scalar α,

L(α v) = α L(v) (15) and L(v + w) = L(v) + L(w).
(16)

Transformation language

This section defines a transformation language based on the three notions of rule, strategy and transformation. By a rule we mean a classical rewrite rule. A strategy is a way to control how rules are applied. Strategies can be combined to define strategies with a finer control or a more powerful effect. We propose easy-to-remember names for the most popular -and indeed most useful -strategy constructors and combinators. The user can also extend the language with other strategies.

It is not obvious that any natural transformation of mathematical expressions and models can be concisely expressed as a rewriting-based strategy, in a natural way. Moreover what is exactly a strategy is not completely clear from the literature, and the name of "strategy" for a formal transformation may lead the user of our language to confusion. Therefore we consider an a priori independent notion of transformation. Transformations have the following three features: (i) a transformation is reproducible, (ii) a transformation may not progress, and (iii) a transformation may not terminate. The reproducibility feature (i) means that each application of a given transformation to a given expression produces the same effect: either it does not terminate each time, or it produces each time the same expression. In particular, this property excludes non-determinism from the notion of transformation. By contrast a strategy may be non-deterministic (see, e.g. One[s] in [START_REF] Balland | Tom: Piggybacking rewriting on Java[END_REF]). It would be enriching to develop a complete theory of transformations but it exceeds the scope of the present paper on rule-based transformations.

We propose an implementation of this transformation language as a new package for the Maple TM computer algebra system. The package is named symbtrans, for "symbolic transformations". 2 We strongly rely on the functional features of Maple TM by providing rules, strategies and transformations as Maple functions, possibly through higher-order functions constructing them from another representation. Functions faithfully provide the expected feature (i) of transformation reproducibility. Whether a transformation output differs or not from its input (feature (ii)) is controlled by the Maple TM mechanism of exceptions. Feature (iii) is left under the responsability of users that can interrupt execution with the Maple TM function timelimit.

Top rewriting

The Maple TM statement ruleName := [l,r] declares the rewrite rule l → r as a pair and assigns it the name ruleName. The function Transform associates to any such pair the function applying the corresponding rewrite rule at the top of any term: Given a term t, the function application Transform(ruleName)(t) applies the rewrite rule l → r to t, as defined in Definition 1. If the rule cannot be applied i.e. if t does not match its left side l, then the exception "Fail" is raised. This is the standard rewriting at the top or top rewriting strategy.

Example 5. Consider the property v + w dx = v dx + w dx of linearity of the integral. The rewrite rule corresponding to its application from left to right can be defined with symbtrans as the pair IntegralLinearity := [ Integral(A_ + B_, C_), Integral(A, C) + Integral(B, C)]; A convention in the package is that variable names end with " " in order to distinguish them from constants. In order to apply the IntegralLinearity rule at the top to the term t := Integral(v(x)+w(x),x); we write Transform(IntegralLinearity)(t). The resulting term is Integral(v(x),x) + Integral(w(x),x)

Elementary transformations

The two elementary transformations Identity and Fail are defined as follows.

Identity(t) = def t Fail(t) = def error "Fail"
The first one has no effect, since it transforms t into itself. The second one always fails and raises the exception "Fail". This exception is raised each time a transformation fails transforming a term. What is a failure for a transformation has to be defined for each transformation, as previously done for top rewriting.

It is sometimes useful to consider the non-progress of a transformation as a failure with the aim to handle this exception and enable other transformations. This feature is realized by the transformation combinator IdentityAsFail defined by IdentityAsFail(s)(t) = def if s(t) = t then Fail(t); else s(t); for any transformation s and any term t.

Conversely, it is also convenient to hide at some higher level the exception "Fail" raised by a transformation s. This is the purpose of the FailAsIdentity combinator: FailAsIdentity(s)(t) = def try s(t); catch "Fail" : t;

Transformation combinators

In this section we define three transformation combinators. They take transformations as parameters and control their order of application. Thus they help defining complex transformations by combination. They are generic in the sense that they do not depend on the nature and structure of the terms they are applied on. They are defined as follows.

LeftChoice([s 1 , . . ., s n ])(t) = def try s 1 (t); catch "Fail":LeftChoice([s 2 , . . ., s n ])(t) LeftChoice([ ])(t) = def Fail(t); Comp([s 1 , . . ., s n ])(t) = def Comp([s 2 , . . ., s n ])(s 1 (t)); Comp([ ])(t) = def t; STNormalizer(s)(t) = def if s(t) = t then t; else STNormalizer(s)(s(t));
The transformations LeftChoice([s]) and Comp([s]) are defined by induction on n ≥ 0 for any sequence of transformations s = (s i ) i=1,...,n . If n = 0, both do nothing. Otherwise, the application of LeftChoice([s]) to the term t returns the result s i (t) of the application of the first transformation in the sequence s which succeeds on the term t. It reports a failure if no one succeeds. When n ≥ 1 the transformation Comp([s 1 , . . ., s n ]) applies s 1 , s 2 , etc in sequence, until application of s n or a previous failure. The transformation STNormalizer(s) iterates the application of the transformation s until the latter fails or a fixed point is reached. The transformation STNormalizer(s) fails if and only if the transformation s fails during these iterated applications. To avoid the failure of the transformation s when computing the normal form of a term t with respect to s, one should write STNormalizer(FailAsIdentity(s))(t). The transformations involving STNormalizer presents a risk of non termination and should therefore be carefully employed. Notice that the above definition of STNormalizer is just a specification. For more efficiency, the implementation computes s(t) only once.

Traversal transformations

This section introduces transformations called traversal or term transformations because they explore the structure of the term they are applied on. We provide three traversal transformation constructors: Some, Outermost and Innermost. By extension of a classical terminology in rewriting theory, a redex of a term t for a transformation s is a subterm of t that can be transformed by s, i.e. on which the transformation s does not fail.

The transformation Some(s) tries to apply the transformation s to all the immediate subterms of a term t. It fails if all these applications fail, or if t is a constant or a variable. It is defined by Transform(IntegralLinearity)(t+2) raises the exception "Fail" because the rewrite rule cannot be applied at the top of the expression t+2. A solution is to replace it with the statement Some(Transform(IntegralLinearity))(t+2) which produces v(x) dx + w(x) dx + 2 since the expression t+2 is viewed as the term +(t, 2). The transformation Some is not very useful in practice. Its main purpose is to shorten the definition of the other traversal transformations, namely Outermost and Innermost.

Some(s)(t) = def if t = f(t 1 ,. . .,t n ) then if ∀i ∈ {1, . . . ,
The transformation Outermost(s) is very common in symbolic computation. It applies the transformation s once to all the redexes of t for s that are the closest ones to the root of t, i.e. to the largest subterms of t on which s succeeds. In other words the transformation Outermost traverses the term t down from its root and tries to apply s to each traversed subterm. If the transformation s succeeds on some subterm t ′ of t, then it is not applied to the proper subterms of t ′ . In particular, Outermost(s) fails if and only if s fails on all the subterms of t. It can be formally defined by

Outermost(s) = def LeftChoice([s,Some(Outermost(s))]) Example 7. Let t = 2 + v(x) + 3( w(x) + g(x) dx r1
) dx r0 be a term with two redexes r 0 and r 1 for the rule of integral linearity. Then Outermost(IntegralLinearity)(t)

gives the expression 2 + v(x) dx + 3( w(x) + g(x) dx) dx, since the rule IntegralLinearity is only applied to the outermost redex r 0 of t.

The strategy Innermost(s) works similarly, but in the opposite direction, i.e. it traverses a term t up from its smallest subterms and tries to apply the strategy s once to the smallest redexes of t for s. It is formally defined by Innermost(s) = def LeftChoice([Some(Innermost(s)),s]) Example 8. For the term t of Example 7, the expression Innermost(IntegralLinearity)(t) applies the rule IntegralLinearity only to the innermost redex r 1 of t and gives the expression 2 + v(x) + 3( w(x) dx + g(x) dx) dx.

We provide additional traversal transformations, namely All, TopDown and BottomUp. They are less useful in practice, but we include them in symbtrans because they exist in other strategy languages, e.g. the Tom strategy language [START_REF] Balland | Tom: Piggybacking rewriting on Java[END_REF].

The transformation All(s) applies the transformation s to all the immediate subterms of any term t, and it fails if and only if one of the applications to the immediate subterms fails. It is defined by

All(s)(t) = def if t = f(t 1 ,. . .,t n ) then f(s(t 1 ),. . .,s(t n ));
else // t is a constant or a variable t; fi;

The main purpose of the transformation All is to simplify the definition of the transformations TopDown and BottomUp.

The transformation TopDown(s) tries to apply the transformation s to all the subterms of any term t, at any depth, by starting with the root of t. It fails when there is a subterm of t where s fails. It is defined by

TopDown(s) = def Comp([s,All(TopDown(s))])
The transformation BottomUp(s) behaves similarly, but works in the opposite direction, i.e. it starts from the leaves (the smallest subterms) and goes up.

BottomUp(s) = def Comp([All(BottomUp(s)),s])

Advanced features

In this section we introduce some technical features of symbtrans, namely, the ability to combine rewriting and procedural programming, the matching modulo associativity and commutativity and the conditional rewriting.

Procedural programming and strict evaluation

It is possible to combine rewrite rules and procedural programming in symbtrans. That is, the right-hand side of the rewrite rules may contain calls of Maple TM predefined functions or of functions defined by the user. This feature is not available in pure rewriting languages, e.g. in Maude [START_REF] Clavel | All About Maude -A High-Performance Logical Framework, How to Specify, Program and Verify Systems in Rewriting Logic[END_REF], but it is available in rewriting languages built upon a host language, e.g. ρ-log [START_REF] Marin | Rule-based programming with mathematica[END_REF] which is built on Mathematica and Tom [START_REF] Balland | Tom: Piggybacking rewriting on Java[END_REF] which is built on Java.

When we declare a rewrite rule whose right-hand side contains some function calls a problem occurs: since Maple TM is a strict evaluation language, it completely evaluates all the sub-expressions of an expression before evaluating the expression itself. Therefore it evaluates the function calls present in the right-hand side of the rewrite rules at the declaration of the rewrite rule, whereas the expected behaviour is most often evaluation of function calls at the application of the rewrite rule.

Example 9. Let us consider the following rule declaration:

e := [L_,nops(L)]; where nops is the Maple TM function that computes the number of arguments of L. When L is a list nops(L) computes the length of L. The problem is that Maple TM immediately evaluates nops(L) and the rule e becomes [L_,1]. The expected behavior is that the function nops is applied after the rule application, i.e. after instantiation of the variable L by the substitution that arises from the matching of L with a given term.

The solution we propose is to write DelayEval -> r instead of r for the right-hand side of a rewrite rule that contains function calls. In the example, the declaration of e would be e := [L_, DelayEval -> nops(L)];

The Maple TM expression x -> r denotes the λ-term3 λx.r, and this solution works for two reasons. First, the evaluation of r in λx.r is delayed until an argument is provided to this λ-term. Second, Maple TM accepts the application of substitutions to λ-terms. With these two ingredients we can correctly implement the delayed evaluation in the rewrite rules as follows. The implementation of the application of the rule [l, DelayEval -> r] to the term t is done via the following steps:

(1) Compute a solution σ to the matching problem l ≪ t, (2) Compute the application σ(r) of σ to r, (3) Return the application (DelayEval -> σ(r))(DelayEval) of the λ-term to the protected variable DelayEval.

Since DelayEval is a protected variable in symbtrans and thus never appears in σ(r), the expression (DelayEval -> σ(r))(DelayEval) is β-reduced4 by Maple TM to σ(r).

It is worth mentioning that enclosing with unevaluation quotes '...' each function present in the right-hand side r of a rule, or using the parameter modifier ::uneval, does not give the required solution. If an unevaluated function enclosed by the quotes appears as an argument of a term, then Maple TM eliminates the quotes and evaluates this function. As a consequence, applying a substitution to the term r provokes an early evaluation of the functions in r enclosed by quotes.

Extension to associative and commutative function symbols

When applying a rewrite rule [l → r] to a term t, and terms l and/or t contain associative and commutative function symbols, such as + or * , then the matching problem l ≪ t and the rule application have to be done modulo associativity and commutativity.

The following definition generalizes Definition 1. It defines pattern-matching and rule application modulo a theory as in (Cirstea and Kirchner, 2001).

Definition 10. Let T be an equational theory. For two terms t and t ′ in T (F, X ) the problem of finding substitutions σ such that the equality σ(t ′ ) = t is a logical consequence of the axioms of T is called a matching problem modulo T and is denoted t ′ ≪ T t. Such a substitution σ is called a solution of the matching problem t ′ ≪ T t. The application of the rewrite rule l → r to a term t, denoted by [l → r](t), is the (possibly empty) set {σ 1 (r), σ 2 (r), . . .} where each σ i is a solution of l ≪ T t.

If the theory T axiomatizes the associativity and commutativity of some function symbols, then the set of solutions of the matching problem is finite. The associative and commutative symbols of symbtrans presently are * , +, ∪ and ∩. The corresponding theory is denoted by AC.

To deal with associative and commutative symbols the package symbtrans is equipped with two matching functions. The function Matching(t')(t) returns only one solution of the matching, always the same one, because transformations are expected to be reproducible. The function MatchingAll(t')(t) returns all the solutions. The first function is sufficient in many practical cases, but the second one is sometimes useful. In particular it is necessary for the application of conditional rules (see next section).

The present implementation of these two functions in symbtrans does not include optimizations suggested in the literature. In [START_REF] Belkhir | Lazy AC-pattern matching for rewriting[END_REF], the authors suggested a more flexible approach. They provided a lazy matching algorithm modulo associativity and commutativity. Lazy means that the algorithm produces only a first solution and a way to get the other ones. We plan to integrate this lazy algorithm in a future version of symbtrans.

Conditional rewriting

In symbtrans rewrite rules can be conditional. A conditional rule is a rewrite rule l → r and a condition c on the variables in l. In symbtrans a conditional rule is declared as a list [l,r,c] of three elements, where the third element is the condition. Notice that the DelayEval mechanism is usually required in the condition when it contains function calls.

Example 11. The linearity property (15) of an operator B can be encoded by the following conditional rewrite rule replacing B(αx) by αB(x) if α is a scalar:

[B(Alpha_*X_), Alpha*B(X), DelayEval -> IsScalar(Alpha)].

The application Transform([l,r,c])(t) of the conditional rule [l,r,c] to a term t, where c is a Boolean condition, gives a term σ(r) where σ is a solution of the matching problem l ≪ AC t such that σ(c) holds, or raises the exception "Fail" if there is no such solution.

Formal proof examples

As typical examples, we consider the mathematical proofs reproduced in Appendices A and B. They are respectively formalized in Appendices C and D as sequences of transformations with the symbtrans language. Appendix E shows the collection of rules and transformations used in the proofs. The present section explains this formalization.

Mathematical operators

The algebraic properties of the integration, derivation, and general summation operators are encoded with symbtrans in Appendix E, Section E.1. Two versions of the Green rule are given in Appendix E, Section E.2: the usual rule GreenRule, and a parametrized conditional one CondGreenRule. The pattern patt is used to ensure that the Green rule is applied in the right way.

About linearity

We recall that an operator L is said to be linear if ( 15) and ( 16) hold. The linear operators used in the proofs are T , T * , B, Ω , ∂ x and i . At the present level of formalization, scalars are not distinguishable from other symbolic expressions. As a consequence Eq. ( 15) cannot be turned into a general rewrite rule: otherwise this rule could also produce the unexpected term v L(α) from L(α v). We presently address this issue by writing one rewrite rule for each term α of interest in the proofs under consideration. It only concerns Steps 2 and 5 in Appendix B.

On the contrary Eq. ( 16) can be safely expressed by a rewrite rule. However the twoscale transform manipulates many linear operators and it is tedious to define a rewrite rule that expresses the linearity property (16) for each operator. Therefore we provide the generic constructor Linearity(n,fun,t), where fun is a function and t is a term with n (underscored) variables or more. This constructor generates a rewrite rule that states that the operator t is linear with respect to its n th variable, in the sense of fun. Notice that fun is usually the + function provided by Maple TM . However it is possible to have different + functions for different vector spaces. For instance, the Maple TM expressions Linearity( 2 For a given operator, the two rewrite rules that correspond to ( 15) and ( 16) can be merged together using transformation combinators. The collection of linear operators as well as their linearity properties is stored in an array named LinearityOf and defined in the file integral.mpl reproduced in Appendix E, Section E.1.

Convergence calculus

The notion of convergence has been introduced in Section 2. The theory of convergence ( 12) and ( 13) corresponds to the following rewrite rules:

O(ε) + O(ε) → O(ε), (17)
n i=1 O(ε) → O(ε), (18)
Ω O(ε) dx → O(ε), and (19)
z * O(ε) → O(ε) (20)
for a term z bounded with respect to ε. These rules are combined in a strategy defined by ConvergenceStrategy := STNormalizer( LeftChoice ([ Outermost(OEpsilonSum), Outermost(OEpsilonSUM), Outermost(OEpsilonIntegral), Outermost(OEpsilonConst) ]) ); where OEpsilonSum, OEpsilonSUM, OEpsilonIntegral and OEpsilonConst respectively encode the rewriting rules ( 17), ( 18), ( 19) and (20). They are defined in the file convergence.mpl reproduced in Appendix E, Section E.3. The result ConvergenceStrategy is a powerful strategy that reduces O(ε) terms as much as possible. In the present case it can be shown that it always terminates and thus can be systematically applied after each tranformation step.

We now present and address a problem arising when embedding this strategy within a computer algebra system with strict evaluation. If the notion of a function that tends to 0 when ε tends to 0 is represented by the Maple TM expression O(ε), then Maple TM simplifies any expression O(ε) -O(ε) to zero, whereas O(ε) -O(ε) should be simplified to O(ε). The solution we suggest consists in considering the term O(i, ε) instead of O(ε), where i is a fresh index. The term "fresh" means that the same index has never been produced before to construct such a term. This solution is natural because it basically relies on the mathematical semantics of the term O(ε). That is, two occurrences of O(ε) are two different functions, and it is natural to distinguish them using two different indexes. Technically, we provide a function FreshIndex() that returns a new index at each call. Moreover each occurrence of O(ε) in the right-hand side of a rewrite rule is replaced by O(FreshIndex(), ε).

Two-scale calculus

The algebraic properties of the two-scale operators have been stated and mathematically proved in Appendix A of [START_REF] Lenczner | A two-scale model for an array of AFM's cantilever in the static case[END_REF]. They are formulated as rewriting rules in Appendix E, Section E.4.

The two-scale limit of the gradient operator corresponds to Eq. ( 74) of [START_REF] Lenczner | A two-scale model for an array of AFM's cantilever in the static case[END_REF]. This equality is proved with symbtrans in Appendix C. The mechanism of the formal proof is to reduce the left-and the right-hand sides of the equality A.7, and then to show that their reduced forms are equal up to O(ε).

The formal derivation of the two-scale model of the stationary heat equation with symbtrans is reproduced in Appendix D. There, the transformation steps are applied to the term LHST -RHST, where LHST (resp. RHST) is the left-(resp. right-) hand side term of the heat equation ( 1). The two-scale model derivation uses the weak convergence property of the derivative operator (Step 6 in Appendix B) as a lemma. In the formal proof this lemma is written "by hand" as a rewriting rule after being formally proved.

Final remarks

Some proof steps require to apply some equation A = B from left to right and from right to left. However, including the two rewriting rules A → B and B → A in the same strategy may induce non-termination. This is typically the case for the linearity properties. In this case, non-termination is avoided by the introduction of a more specialized version of the rules corresponding to the second orientation, such as: AdhocSimplify := [Integral(Omega_,C_*SUM(F_,J_,D_),X_), SUM(Integral(Omega,C*F,X),J,D)]; Notice that the proofs use two Maple TM functions: has(expr1,expr2) that returns true if expr2 is a subexpression of expr1, and Evala(Expand(expr)) that expands products and powers of rational functions with algebraic coefficients.

We have also developed a formal proof of the property (7) that the two-scale weak limit u 0 of T u ε is independent of y. This proof is not reproduced here. Instead, ( 7) is encoded as a lemma (See Appendix E.5) used in the main proofs. The hypotheses of these proofs are gathered in the file hypothesis.mpl reproduced in Appendix E, Section E.6.

Theoretical basis and related work

The theoretical basis for this work is the deterministic fragment of the ρ AC -calculus (Cirstea and Kirchner, 2001), where AC is the theory axiomatizing the associativity and commutativity of the symbols +, * , ∪ and ∩. The ρ T -calculus, where T is an equational theory, is an extension of the λ-calculus; one abstracts on a pattern rather than on a single variable. The abstraction mechanism is based on the rewrite rule l → r, also viewed as a ρ-term. Notice that when l is just a variable x, this ρ-term corresponds to the λ-term λx. r. Moreover, the ρ T -calculus considers higher-order terms, i.e. terms that may contain abstractions and rule applications.

When an abstraction l → r is applied to a ρ-term t, which is denoted by [l → r]t, the matching mechanism is based on the binding of the free variables of l to the appropriate subterms of t. This matching is done modulo the theory T. The latter is often expressed by algebraic axioms such as associativity and/or commutativity. In symbtrans both the left-hand side term of a rule and the term under rule application corresponds to firstorder ρ-terms, i.e. they contain neither abstractions nor rule applications. However the right-hand side of a rule may be a higher-order ρ-term, since it may contain function calls. Those functions are nothing but λ-terms. Handling the priority between rule application and β-reduction is explained in the steps ( 1), ( 2) and (3) in Section 4.1. Despite the fact that strategies can be encoded with the ρ T -calculus [START_REF] Cirstea | Rewrite strategies in the rewriting calculus[END_REF] by means of some constructors, we preferred encoding the strategies in symbtrans by means of Maple functions for the sake of efficiency. Finally, we notice that in the symbtrans language, if a rule cannot be applied to a term then the exception Fail is raised. This makes a subtle difference with the semantics of the ρ T -calculus that consists in returning an empty set in this case. The problem of the ρ T -calculus approach is that we can not distinguish between an empty set which is a mathematical term that could arise from the symbolic transformations, and the empty set which denotes the failure of the application of a rule.

The proposed transformation language does not claim for originality. It is deliberately an adaptation for Maple of popular strategy languages such as ρ-log [START_REF] Marin | Rule-based programming with mathematica[END_REF] or Tom [START_REF] Balland | Tom: Piggybacking rewriting on Java[END_REF]. But, departing from Tom which extends an host language with an additive syntax, our transformation language smoothly integrates with standard Maple functions. Consequently, the Maple programmer learns it quickly, and is free to mix function-and rule-based programming styles. Moreover all the features of her development environment (such as refactoring, code completion, dependency analyses, etc) are preserved for free.

The closest implementation is ρ-log, a package developed upon the advanced rewriting kernel of Mathematica. It supports non-deterministic and conditional rewriting. The main drawback of ρ-log is that it considers the non-applicability of a rule as the identity. Technically speaking, the strategy FailAsIdentity is implicitly applied to all the transformations. However, when a transformation returns the same term given as an input, we do not know if this transformation fails or it performs some modifications and then returns the same term. Moreover in ρ-log it is not possible, at least in a straightforward way, to do higher-order rewriting, since the rewriting rules are not directly accessible to the user: They are declared by means of the constructor DeclareRule.

Conclusion

Our main motivation for the development of a transformation language in Maple TM was to facilitate the design of the MEMSALab software dedicated to the automatic derivation of multi-scale models. However symbtrans is a general tool that can be used by Maple TM programmers and mathematicians in the formalization of equational reasoning. It makes it possible to express rule-based symbolic computations in a concise and natural way, thus providing a good guarantee of the correctness of the formal proofs with respect to their hand-written counterparts. Since the symbtrans package is written in Maple TM , it obviously does not extend the expressivity of the Maple TM language, but it clearly increases readability and conciseness. Although this paper presents an implementation in Maple TMTM , the transformations presented here could easily be developed in a similar way in any functional language.

The transformation language symbtrans allows the derivation of the weak two-scale limit of the derivative operator and the two-scale model of the stationary heat equation at the same "level" as the mathematical derivations. The word "level" covers three aspects: the formal and the hand-written proofs have almost the same size, they follow the same steps, and the strategy term written at each step of the formal proofs is a natural formalization of its mathematical counterpart. The symbtrans package is used in (Yang et al., 2011a,b) to formally derive the two-scale model of the stationary heat equation in a region composed of a thin part and a part with periodically distributed holes.

For a more scalable treatment of linearity we plan in a near future to detect the scalar nature of terms by assigning a type to each expression. More generally a type system for mathematical expressions is under way. We plan to transform each proof into a module whose execution produces a parametrized rewrite rule. The latter can be instantiated and applied in other proofs.

where n y stands for the unit outward normal vector to the boundary ∂Y of Y . Remarking that ∂ yi u 0 vanishes5 and that v vanishes on ∂Y , the second subterm vanishes and we obtain Ψ = - By linearity of ∂ and B, and since v 0 does not depend on y, we get after application of (13) and simplifications,

Example 2 .

 2 Let f , g , + and * (resp. a and b) be function symbols of arity 2 (resp. 0) in F. Let x and y be variables in X . • The application [x → b](f (a, b)) of the rule x → b to the term f (a, b) is the term b. The substitution σ = {x → f (a, b)} is an obvious solution of the matching problem x ≪ f (a, b). • The application [f (x, a) → g(a, x)](f (b, a)) of the rule f (x, a) → g(a, x) to the term f (b, a) is the term g(a, b), with the substitution σ = {x → b}. • The application [a → b](a) of the rule a → b to the term a is the term b, with the empty substitution σ = {}.

  ,x->y->x+y,Integral(Omega_,_,Z_)); and Linearity(1,x->y->x+y,T(_)); respectively express the linearity property (16) of the Integral and T operators with respect to +. Their evaluation respectively produces the rule [Integral(Omega_,X_+Y_,Z_), Integral(Omega,X,Z)+Integral(Omega,Y,Z)] and [T(X_+Y_), T(X)+T(Y)]. a

∂∂∂∂

  xj u 0 ∂ yi v i dxdy + xj ∂ yi v i dxdy + O(ε).• Step 9. Similarly, repeating the Green formula application until no derivative is left on the test function v, i (nx ) i ds(x)dy -Ω×Y ∂ xi u 0 v i dxdy + Ω×∂Y u 1 v i n yi dxds(y) -Ω×Y ∂ yi u 1 v i dxdy + n j=1 Ω×∂Y y j ∂ xj u 0 v i n yi dxds(y) -n j=1 Ω×Y ∂ yi (y j ∂ xj u 0 ) v i dxdy + n j=1 ∂ Ω×∂Y u 0 v i n yi y j n xj ds(x)ds(y) yi (y j u 0 ) n xj v i ds(x)dy n j=1 Ω×∂Y y j ∂ xj u 0 v i n yi dxds(y) yi (y j ∂ xj u 0 ) v i dxdy + O(ε).Since v vanishes on all boundaries, xi u 0 v i dxdy+ Ω×Y ∂ yi u 1 v i dxdy + n j=1 Ω×Y ∂ yi (y j ∂ xj u 0 ) v i dxdy n j=1 Ω×Y ∂ yi (y j ∂ xj u 0 ) v i dxdy + O(ε), xi u 0 v i dxdy + Ω×Y ∂ yi u 1 v i dxdy + O(ε).Finally, thanks to linearity of integral and by factoring v i , u 0 + ∂ yi u 1 )v i dxdy + O(ε).B. A mathematical two-scale transformation of the heat equationWe detail the proof of Proposition 4. We start with test functionsv 0 ∈ C ∞ 0 ( Ω) and v 1 ∈ C ∞ ( Ω × Y ) that is Y -periodic. • Step 1. We choose v = B(v 0 + εv 1 )as a test function in the weak formulation (1) of the model, n i=1 Ω a(∂ xi u) ∂ xi B(v 0 + εv 1 ) dx = Ω f B(v 0 + ǫv 1 ) dx. • Step 2. Applying the rule (A.3) of partial derivatives of B(v) yields n i=1 Ω a(∂ xi u) (B(∂ xi (v 0 + εv 1 )) + 1 ε B(∂ yi (v 0 + εv 1 )) dx = Ω f B(v 0 + ǫv 1 ) dx.

  xi u) (B(∂ xi v 0 ) + B(∂ yi v 1 )) dx = Ω f B(v 0 ) dx + O(ε).

From# 2 .

 2 xi u) B(∂ xi v 0 + ∂ yi v 1 ) dx = Ω f B(v 0 ) dx + O(ε).• Step 3. The zero-order approximation (A.5) of B byT * implies n i=1 Ω a(∂ xi u) T * (∂ xi v 0 + ∂ yi v 1 ) dx = Ω f T * (v 0 ) dx + O(ε).• Step 4. Now, we apply the definition of the adjointT * of T , n i=1 Ω×Y T (a∂ xi u) (∂ xi v 0 + ∂ yi v 1 ) dxdy = Ω×Y T (f ) v 0 dxdy + O(ε).• Step 5. From the identity (A.2) and the assumptions (11) we getn i=1 Ω×Y a 0 T (∂ xi u) (∂ xi v 0 + ∂ yi v 1 ) dxdy = Ω×Y f 0 v 0 dxdy + O(ε).• Step 6. From the approximation (A.7) of the derivative operator applied to the test function a 0 (∂ xi v 0 + ∂ yi v 1 ) we get the wanted two-scale modeln i=1 Ω×Y a 0 (∂ xi u 0 + ∂ yi u 1 ) (∂ xi v 0 + ∂ yi v 1 ) dxdy = Ω×Y f 0 v 0 dxdy + O(ε). C. Aformal proof of the derivative weak convergence property # File gradient.mpl # Contributors: Walid Belkhir, Alain Giorgetti and Michel Lenczner # 1. Library loading new_lib_dir := "../lib": libname := new_lib_dir, libname: with(stmodule); Rewrite systems loading read 'integral.mpl': read 'Green.mpl': read 'twoscale.mpl': read 'convergence.mpl': read 'lemmas.mpl': read 'hypothesis.mpl':

  n} s(t i ) fails, then Consider again the rewrite rule IntegralLinearity of Example 5, encoding the linearity of the integral, and the term t encoding v(x) + w(x) dx. Notice that the statement

	Fail(t);
	else
	f(FailAsIdentity(s)(t 1 ),. . .,FailAsIdentity(s)(t n ));
	fi;
	else // t is a constant or a variable
	Fail(t);
	fi;
	Example 6.

We mention that we use the minimal scope convention of the derivative operator, i.e. ∂x i f g means (∂x i f )g.

The package is available as an archive file upon request to the authors. It can be executed on any machine running Maple TM .

We recall that a λ-term is either a variable x, an abstraction λt.u, or an application (u v), where u and v are two λ-terms.

The β-reduction is the reduction rule (λt.u)v ❀ β u[v/t], where u[v/t] is the replacement of t with v in u.

In the formal proof this simplification is done after the successive applications of Green rule and the elimination of the boundary terms, i.e. at the end of step 9.

Acknowledgements

The authors would like to thank the anonymous reviewers of a previous version of this paper for their helpful comments, and R. N. Dhara and B. Yang for the feedback provided by their development activity using symbtrans.

⋆ This work has been partially funded by the Labex ACTION (ANR-11-LABX-01-01), the INTERREG IV project OSCAR and the PPF MIDi.

A. A mathematical proof of the derivative weak convergence property

This section is devoted to the proof of Proposition 3. We start with some reminders of mathematics that complement those in Section 2.2.

Here the two-scale transform T can be viewed as a linear continuous operator from L 2 (Ω) into L 2 ( Ω×Y ), as such its adjoint T * is a linear continuous operator from L 2 ( Ω×Y ) into L 2 (Ω) defined by (A.1)

The two-scale transform can also be defined on integrable functions and it satisfies the property T (u v) = T (u) T (v) for any u, v ∈ L 2 (Ω). (A.2)

Then, we define the so-called regularized inverse two-scale transform B :

x/ε). It can be easily checked that B is a linear operator. The partial derivatives of B(v) for any sufficiently regular function v can be derived by applying the chain rule

It is also useful to know how the null condition of a function v(x, y) on the boundary ∂( Ω × Y ) is transfered to its range by B:

In the following lemma (admitted), O(ε) denotes any function that vanishes in the L 2 (Ω)norm when ε tends to zero.

Lemma 12. The operator B is a zero-order approximation of the adjoint operator T * in the sense that

) for any sufficiently regular and Y -periodic function v. Moreover, B(v) can be approximated at the first-order by

In the following, for simplicity we write u instead of u ε . We shall prove Proposition 3 or equivalently, by the density of the set C ∞ 0 ( Ω×Y ) n of infinitely continuously differentiable functions will all derivatives vanishing on ∂(

•

Step 1. Applying the definition of T * to the left-hand side Ψ of (A.7) yields

• Step 2. From the approximation (A.5) of T * (v i ) by B(v i ), the linearity of integral, the boundedness of ||∂ xi u|| L 2 (Ω) and the property (13) of O(ε) we get

• Step 3. Then, we apply the Green formula ( 14) and get

and the terms on the boundary are removed thanks to Property (A.4),

• Step 4. From the expression (A.3) applied to the partial derivatives of B(v i ) and by linearity of integral,

and

• Step 5. We apply (A.5) and (A.6) to approximate B(∂ xi v i ) at the zero-order and B(∂ yi v i ) at the first-order together with the rule (13) and thus get

and

Thanks to the linearity of T * ,

Grouping Ψ 1 and Ψ 2 ,

• Step 6. From the definition of the dual operator

After expanding and applying linearity of integral,

• Step 7. We use the zero-order approximation

and the first-order approximation, both deduced from (10), respectively in the first and third integrals and in the second integral to get

After simplification

• Step 8. We apply the following instance of the Green formula to the second subterm,


Alois Würger

Thermodiffusion in binary liquids: The role of irreversibility

Keywords: numbers 66.10.C-, 05.70.Ln, 82.70.-y numbers:

We study thermal diffusion in binary mixtures in the framework of non-equilibrium thermodynamics. Our formal result displays the role of partial enthalpies hi and Onsager's generalized mobilities Ai. The ratio A1/A2 provides a measure for the irreversible character of thermal diffusion. Comparison with experimental data on benzene, cyclohexane, toluene, and n-alkanes shows that irreversibility is essential for thermal diffusion, and in particular for the isotope effect.

Introduction. Thermal diffusion, or the Soret effect, describes the mass flow induced by a temperature gradient in complex fluids [1][2][3]. Together with the Dufour effect that acccounts for heat flow in a concentration gradient, it constitutes a classical example of Onsager's reciprocal relations for non-equilibrium systems [4]. Thermal diffusion is one of the mechanisms governing the compositional grading in the Earth's petroleum reservoirs [START_REF] Galliéro | Society of Petroleum Engineers EUROPEC/EAGE Annual Conference[END_REF] and the isotope fractionation in silicate melts [START_REF] Dominguez | [END_REF]. In recent years, thermal diffusion has proven a versatile means for confining colloidal suspensions [7,8], manipulating DNA through nanopores [9], or self-propelling Janus particles [10,11].

The Soret effect of macromolecules and colloidal particles is to a large extent determined by viscous effects and thus can be treated in the framework of macroscopic hydrodynamics [3]. Exploiting the reciprocal laws for heat and mass flows, Derjaguin related thermally driven transport to the solute-solvent interaction enthalpy [START_REF] Derjaguin | Surface Forces[END_REF]. Much recently much progress has been made for charged solute, confirming Derjaguin's picture of enthalpy flow as driving mechanism and pointing the out the role of the electrolyte Seebeck effect [START_REF] Putnam | [END_REF][14][15].

Much less is known, however, on thermal diffusion in molecular mixtures. In spite of the many experimental [16][17][18][19][20][21][22][23] and molecular dynamics data [24][25][26] on the composition, temperature, and mass dependencies, available models diverge on the drivinng mechanism [27][28][29][30][31][32][33][34]. Although the observed isotope effect [18] indicates irreversible behavior, there is at present no general agreement whether the Soret effect can be described by equilibrium thermodynamics, or on the contrary reflects the irreversible character of the underlying diffusion process.

The present Letter discusses these questions on the basis of formally exact expressions derived by de Groot & Mazur half a century ago, yet which so far have been given no attention. As a main purpose, we separate equilibrium and irreversible factors of the Soret coefficient. Comparison with recent experimental findings for binary mixtures, provides a simple physical picture for the driving mechanism and a clear signature of irreversibility.

Consider a two-component molecular liquid with volume fractions c 1 and c 2 = 1 -c 1 . Non-uniform composi-tion and temperature induce a current

J 1 = -D∇c 1 -c 1 c 2 D T ∇T. (1)
The first term or gradient diffusion tends to smoothen composition inhomogeneities, whereas the second one, or thermal diffusion, pushes one component to the cold and the other to the hot, and thus favors separation. The steady state J 1 = 0 is characterized by a spatial composition gradient ∇c 1 = -c 1 c 2 S T ∇T , where S T = D T /D is the Soret coefficient. Even for simple systems such as benzene-cyclohexane, both its sign and its dependence on composition and temperature lack a rationale so far.

Non-equlibrium thermodynamics. The coefficients D and D T are related to the thermodynamic forces introduced by Onsager and worked out in detail by de Groot and Mazur. As a generalized state function, the Planck potential µ i /T of component i is given by the ratio of its chemical potential µ i (or partial Gibbs energy) and the absolute temperature T . Just like the gradient of gravitational or electric potentials give mechanical forces, that of the Planck potential defines thermodynamic forces acting on each molecular species. For an otherwise homogeneous system one has

∇ μi T = - ĥi T 2 ∇T + μc ii T ∇c i , (2)
where the first term arises from the Gibbs-Helmholtz relation d(μ i /T )/dT = -ĥi /T 2 and the second one from applying the Gibbs-Duhem relation to μc ij = dμ i /dc j . For incompressible liquids, it is convenient to use volume specific quantities, e.g., μi = µ i /v i the volume per molecule v i ; thus ĥi is the corresponding enthalpy density.

The thermodynamic force on species j induces a current a ij ∇(μ j /T ) of species i, where the mobility matrix a ij is symmetric and positive definite. Eliminating the related heat flow and accounting for J 2 = -J 1 results in [4]

J 1 = -(c 2 a 11 -c 1 a 12 ) ∇ μ1 T + (c 1 a 22 -c 2 a 21 ) ∇ μ2 T . (3)
Inserting (2) and separating thermal and concentration gradients, one readily identifies D and D T in (1),

D = c 1 (c 1 a 22 -c 2 a 12 ) + c 2 (c 2 a 11 -c 1 a 12 ) c 1 c 2 T c 1 μc 11 , (4)
D T = (c 1 a 22 -c 2 a 12 ) ĥ2 -(c 2 a 11 -c 1 a 12 ) ĥ1 c 1 c 2 T 2 . ( 5
)
These relations for the diffusion and thermodiffusion coefficients have been obtained by de Groot and Mazur, albeit for mass instead of volume fractions [35]. They reveal several remarkable features. First, the coefficient D T is proportional to the partial enthalpies of the components; this confirms the central role of enthalpy for thermo-osmosis pointed out by Derjaguin [START_REF] Derjaguin | Surface Forces[END_REF].

Second, the factors a ij accounting for irreversibility, appear in both D and D T in the form

B 2 = c 1 a 22 - c 2 a 12 and B 1 = c 2 a 11 -c 1 a 12 .
In view of a two-paramter model discussed below, we absorb the chemical potential derivative in these coefficients and separate the dynamic viscosity η. Thus we put

A i = B i (η/c 1 c 2 T )c 1 μ11 and have for the diffusion coefficient D = c 1 A 2 + c 2 A 1 η . (6)
Note thato in general both a ij and A i are functions of composition and temperature. Inserting the A i in the thermal diffusion coefficient we have

D T = 1 c 1 μc 11 ηT A 2 h 2 v 2 -A 1 h 1 v 1 . (7)
The derivative of the chemical potential can be split in two factors,

c 1 μc 11 = k B T c 1 v 2 + c 2 v 1 Γ, (8)
where the first one arises from the mixing entropy; the remaining factor Γ is related to composition dependence of the "activity coefficient" γ and thus accounts for nonideal behavior [4]. The remainder of this paper deals with the Soret coefficient S T = D T /D. For its derivation it was advantageous to consider volume fractions c i = n i v i , where n i is the concentration. Since Soret data and thermodynamic excess quantities of binary mixtures are usually given in terms of mole fractions x i = n i /(n 1 + n 2 ), we change variables and obtain

S T = 1 Γk B T 2 h 2 A 2 /v 2 -h 1 A 1 /v 1 x 1 A 2 /v 2 + x 2 A 1 /v 1 . (9)
This formally exact expression reveals the intricate equilibrium and irreversibiltiy aspects of S T : The state function enthalpy is an equilibrium property, whereas the mobilities A i account for the irreversible nature of thermal diffusion. The data for benzene-cyclohexane are from Debuschewitz and Köhler [18], those for toluene-hexane from Zhang et al. [17].

The solid lines are calculated from Eq. ( 12) with parameters from the tables. For protenated benzene we have used A1/A2 = 0.765, and for deuterated benzene A1/A2 = 0.785.

It turns out instructive to consider the case where the mobility coefficients are constant and identical to each other. Then the Soret coefficient takes the form

S T = 1 Γk B T 2 h 2 /v 2 -h 1 /v 1 x 1 /v 2 + x 2 /v 1 (A 1 = A 2 ), (10)
which is proportional to the difference in partial enthalpy per volume. In other words, the molecules with the more negative enthalpy density migrates to the cold, similar to gravity-driven sedimentation where the denser component accumulate at the bottom. Remarkably, the mobilities have dropped out, such that S T depends on equilibrium properties only. This expression was first obtained by Haase for gas mixtures, with mole fractions weighted by molecular mass instead of volume [27]; later on it has been refined by several authors [28,29]. Yet the equilibrium hypothesis is not generally valid, and in particular fails for molecular liquids. With the enthalpy and volume parameters given in the Tables, Haase's expression (10) results in a strong positive Soret coefficient for benzene in cyclohexane, whereas the data [18] plotted in Fig. 1 rather show the opposite behavior. Even more strikingly, it misses the isotope effect: Since protonated and deuterated benzene, C 6 H 6 and C 6 D 6 , hardly differ in their enthalpies and molecular volumes, (10) gives the same Soret behavior, whereas the measured values are higher for the heavier isotope.

In the remainder of this paper, we discuss the parameters h i and A i in detail, and then compare Eq. (9) to measured Soret data. The partial enthalpies h i in general depend on composition. In the follow- ing we consider "regular mixtures" where the mean molar value is a parabolic function of the mole fractions, h = x 1 h 0 1 +x 2 h 0 2 +x 1 x 2 h E , with the pure-body enthalpies of molecular interactions h 0 i and the change upon mixing h E . Then the partial enthalpies take the form

h 1 = h 0 1 + x 2 h E , h 2 = h 0 2 + x 1 h E . (11)
Most molecular liquids have positive excess enthalpy of the order of a few kJ/mole, with the noticeable exception of alcohols in water [38]. Equilibrium thermodynamics is mainly concerned with enthalpy differences with respect to some reference state. Eq. ( 9) on the contrary requires absolute values, which are neither easily measured nor easily calculated. In the following we identify h 0 i with the vaporization enthalpies given in Table I. Indeed, the enthalpy change upon removing one molecule from the liquid to the gas state is essentially determined by molecular interactions; typical values correspond to ten times the thermal energy and by far exceed conformational and volume expansion contributions.

The partial molar volumes v i change upon mixing and have been discussed as one source for the composition dependence of S T [34]. For the systems considered here, the excess volume v E is rather small and thus has been neglected from the beginning. For example, for benzenecyclohexane mixtures one has v E = 0.65 cm 3 /mol, which corresponds to less than one percent of the molar volumes [38], whereas the excess enthalpy h E amounts to about 10 % of h 0 i . The thermodynamic factor Γ defined in ( 8) describes the non-ideal behavior of the chemical potential. In analogy to the enthalpy, the excess chemical potential µ E is exptected to be much smaller than the ideal values µ 0 i . This implies that the activity coefficient γ and the quantity Γ vary only weakly with composition. Since the latter appears as an overall factor in S T , its modulation is of little relevance and we may take Γ as a constant. On the contrary, the ideal enthalpies in the numerator of ( 9) largely cancel each other, thus enhancing the weight of the excess h E .

The mobilities A i account for the irreversible character of thermal diffusion. There is no thermodynamic theory for calculating these quantities; they have to be extracted from molecular dynamics simulations or experiments, or to rely on model assumptions. Here we adopt the simplest two-parameter model that consists in taking A 1 and A 2 constants. The experimental values of Table II are determined from data for the mutual diffusion coefficients at low dilution, e.g.,

A 1 = (η 2 /T )D 1 for x 1 → 0.
Comparison with experiment. Inserting the above partial enthalpies and defining ξ = A 1 /A 2 , we obtain

S T = h 0 2 /v 2 -ξh 0 1 /v 1 + h E (x 1 /v 2 -ξx 2 /v 1 ) Γk B T 2 (x 1 /v 2 + ξx 2 /v 1 ) . (12)
In Fig. 1 we compare this expression with data for toluene-hexane [17] and benzene-cyclohexane [18], which show an almost linear variation with the mole fraction of the first component. The theoretical curves are calculated with the parameters given in the Tables; the fit values for the irreversibility parameter A 1 /A 2 and the excess enthalpy h E agree rather well with diffusion and thermometry data. Both the numerator and the denominator of ( 12) contribute equally to the composition dependence.

Haase's expression (10) would result in positive and much too large value of S T . On the other hand, when including the irreversibiltiy parameter Eq. ( 12) provides a good description for the data. The experimental and fit values of A 1 /A 2 agree rather well, yet one should be aware of the uncertainties of the diffusion coefficients D i . As an important result, these fits show that the Onsager mobilities take a ratio which significantly differs from unity. This leads us to the conclusion that irreversibility is crucial for thermal diffusion.

This statement is confirmed by the isotope effect. Since the enthalpy density is insensitive to the molecular mass, the equilibrium expression (10) does not differentiate between protonated and deuterated benzene. On the other hand, the kinetics and thus the mobilities A i depend on the molecular mass; possible mechanisms are the collision rate as discussed in a hard-sphere model [32] and the jump rate of activated diffusion [START_REF] Dominguez | [END_REF]. From thermal diffusion in gases, it is well known that the isotope effect is not an equilibrium property. Finally we emphasize that the slight change from A 1 /A 2 = 0.765 for protonated benzene to 0.785 for the deuterated species results in a composition independent offset of S T , in good agreement with the data. Fig. 2 shows the temperature dependence of Soret data reported by Polyakov et al. for benzene-heptane [20]. An increase of 20 K reduces the Soret coefficient by one quarter, which corresponds to a relative change of about 1.2 % per K. Comparison with the thermal expansion coeffi-cient d ln v/dT ∼ 10 -3 K -1 , suggests that the change of the molecular volume is of little relevance. The temperature derivative of the molar enthalpy h 0 is determined by the heat capacities C P ; with the numbers of Table I one finds that a relative change of about 0.005 K -1 , which accounts for almost half of the variation of the Soret coefficient. The fit curves in Fig. 2 are obtained by varying in addition both A 1 /A 2 and h E , suggesting that the temperature dependence of S T has not a single cause.

In the present work we have discarded viscous effects. This is appropriate for molecular components of similar size, where the two terms of the thermodynamic force (2) carry the same mobilities. On the contrary, the motion of large particles or polymers is accompanied by viscous flow; the Einstein coefficient D carries a size-dependent friction factor, whereas that the thermal diffusion coefficient D T is constant [33]. Mixtures of normal alkanes would provide a model system for the emergence of viscous effects, because of their rather simple mixing properties and the available Soret data [23].

In summary, our comparison with experiments provides strong evidence for the irreversible character of thermal diffusion, even in systems where viscous effects are absent. The simple two-parameter model used here could be improved by considering three mobilities A ij ∼ a ij /c i c j , and adjusting the diffusion and thermal diffusion coefficients independently. We conclude with a remark on the irregular alcohol-water mixtures. In accordance with [START_REF] Derjaguin | Surface Forces[END_REF], the measured Soret coefficient [19] roughly follows the excess enthalpy as a function of composition [38], and in particular reflects the cusp at x eth = 0.15; because of the strongly irregular behavior, a more rigorous comparison would require to go beyond the linear law (11) for the partial enthalpies.

12 FIG. 1 :

 121 FIG. 1: Composition dependence of the Soret coefficient.The data for benzene-cyclohexane are from Debuschewitz and Köhler[18], those for toluene-hexane from Zhang et al.[17]. The solid lines are calculated from Eq. (12) with parameters from the tables. For protenated benzene we have used A1/A2 = 0.765, and for deuterated benzene A1/A2 = 0.785.

FIG. 2 :

 2 FIG.2: Soret coefficient of benzene-heptane mixtures at three temperatures. Data are from Polyakov et al.[20]. The solid lines are calculated from[START_REF] Derjaguin | Surface Forces[END_REF] with parameters given in the Tables.

TABLE I :

 I Pure-component parameters at 25 • C: vaporization enthalpy hvap ; molar volume v; viscosity η; specific heat at normal pressure[38, 39]. Throughout this paper we identify hvap with the interaction enthalpy h 0 .

		hvap ≡ h 0	v	η	CP
		(kJ/mol)	(cm 3 /mol)	(mPa.s)	(J/molK)
	benzene	-33.9	89	0.61	135
	cyclohexane	-33.3	108	0.88	156
	toluene	-38	106	0.59	156
	hexane	-30.7	132	0.31	198
	heptane	-35.4	147	0.39	225

TABLE II :

 II Parameters of protonated benzene (B), cyclohexane (CH), toluene (T), hexane (C6), and heptane (C7). The coefficient Di describes tracer diffusion of component i[START_REF] Funazukuri | [END_REF]. The parameters Ai for tracer diffusison are calculated with the pure solvent viscosity, e.g. A1 = η 2 D1/T . Since ST depends on the ratio A1/A2 only, we don't give absolute values. The viscosity of heptane at 40 • C is 0.33 cP. Data at 25 and 40 • C [37]; d [38]; e Data for B-C 6 at 25 and 40 • C [38]; f [34].

			D1	D2 A1/A2 A1/A2 hE (kJ/mol) Γ
			(10 -9 m 2 /s) (exp) (fit) (exp) (fit)
	B-CH 25 • 1.90 a 2.09 a 0.76 0.765 3.2 d	2.8 1.4
	T-C6 25 • 4.0 b 2.4 b	0.91	0.70		1.0 0.7 f
	B-C7	20 • 30 • 40 •	3.92 c 4.74 c	1.79 c 2.28 c	0.71 c 0.73 c	0.71 0.695 0.67	3.2 e 2.5 e	2.6 2.3 2.0	1 1 1
	a [36]; b [16]; c