TOUFIC Mohamed Amine Cherif
email: mohamedamin.cherif@yahoo.fr

E L Arwadi
email: telarwadi@gmail.com

Jean-Marc Sac-Épée
email: ♯♯jean-marc.sac-epee@univ-lorraine.fr

DIRICHLET-TO-NEUMANN SEMIGROUP ACTS AS A MAGNIFYING GLASS

Keywords: Mathematics Subject Classification. 35J25, 47F05, 47D06. Key words and phrases. Dirichlet-to-Neumann operator, γ-harmonic lifting

The first aim of this paper is to illustrate numerically that the Dirichlet-to-Neumann semigroup represented by P. Lax acts as a magnifying glass. In this perspective, we used the finite element method for the discretization of the correspondent boundary dynamical system using the implicit and explicit Euler schemes. We prove by using the Chernoff's Theorem that the implicit and explicit Euler methods converge to the exact solution and we use the (P1)-finite elements to illustrate this convergence through a FreeFem++ implementation which provides a movie available online. In the Dirichlet-to-Neumann semigroup represented by P. Lax the conductivity γ is the identity matrix I n , but for an other conductivity γ, the authors of [3] supplied an estimation of the operator norm of the difference between the Dirichlet-to-Neumann operator Λ γ and Λ 1 , when γ = βI n and β = 1 near the boundary ∂Ω (see Lemma 2.1). We will use this result to estimate the accuracy between the correspondent Dirichlet-to-Neumann semigroup and the Lax semigroup, for f ∈ H 1/2 (∂Ω).

Introduction.

The semigroup character of a magnifying glass is intuitively quite obvious. From an optical point of view, the magnifying glass has to leave the viewing object unchanged if we look at it from a zero distance, while if we look at it from a distance t, and then from a distance s, it is as if we looked at it from a distance t + s. That is the semigroup character of a magnifying glass. Surprisingly, in his book [10, 36.2], P. Lax constructs an example of a (C 0 ) semigroup who acts in the same way. He takes the Banach space X = C(S n-1 ), the space of continuous functions f on the n-1 dimensional unit sphere or X the Hilbert space L 2 (S n-1 ). For each such function f there is a uniquely determined harmonic function u, called the harmonic lifting of f in the n-dimensional unit ball B which satisfies the following system: ∆u = 0, in B, u(ω) = f (ω), ω in S n-1 .

(1.1)

The Lax semigroup is defined by e -tΛ 1 f (ω) = u(e -t ω) f or ω ∈ S n-1 .

(1.2)

In [10, 36.2 Theorem 1] it is shown that e -tΛ 1 forms a compact (C 0 ) semigroup on X. Furthermore if u is the harmonic lifting of f , then lim t→∞ e -tΛ 1 f -u(0) X = 0.

(1.3)

In fact for any f ∈ X = C(S n-1 ), 1.3 is clear, since u is uniformly continuous. For f ∈ X = L 2 (S n-1 ), we can take for simplicity n = 2 and the Fourier series expansion of f . The unique harmonic lifting of f in S n-1 can be expressed as

u(r, θ) = a 0 2 + ∞ k=1
r k (a k cos kθ + b k sin kθ).

By replacing the Fourier coefficients by their integral expressions

a k = 1 π 2π 0 f (φ) cos kφdφ, b k = 1 π 2π 0 f (φ) sin kφdφ
and using the identity

1 + 2 ∞ k=1 r k cos k(θ -φ) = 1 -r 2 1 -2r cos(θ -φ) + r 2 we get the Poisson integral expression of u u(r, θ) = 1 2π 2π 0 1 -r 2 1 -2r cos(θ -φ) + r 2 f (φ)dφ. Since u(0) = 1 2π 2π 0 f (φ)dφ, we get 2π 0 |u(e -t , θ) -u(0)| 2 dθ = 2π 0 1 2π 2π 0 2e -t (cos(θ -φ) -e -t ) 1 -2e -t cos(θ -φ) + e -2t f (φ)dφ 2 dθ.
Now, for any ε > 0 and any (θ, φ) ∈ [0, 2π] 2 , we can take t large enough such that

2e -t (cos(θ-φ)-e -t )
1-2e -t cos(θ-φ)+e -2t < ε. Consequently, by using the Cauchy-Schwarz inequality, one has

e -tΛ 1 f -u(0) 2 X = 2π 0 |u(e -t , θ) -u(0)| 2 dθ ≤ ε 2π 2π 0 2π 0 |f (φ)|dφ 2 dθ ≤ 2πε 2π 0 |f (φ)| 2 dφ, which proves (1.3).
In each stage of the convergence (1.3), if we take f t (x) = u(e -t x) as a function on the boundary, then according to the semigroup law e -sΛ 1 f t (x) = u(e -(t+s) x), this is just like taking a magnifying glass and going to the center of the disc. For visualization of this phenomena compare the Figures 1, 2 and 3, which present the Lax semigroup at time t = 0, t = 0.5 and t = 2. At the end of the paper an internet address supplies the movie of this phenomena.

In [START_REF] Emamirad | An approximating family for the Dirichelet-to-Neumann semigroup[END_REF], the authors showed that the Lax semigroup is in fact the Dirichlet-to-Neumann (DtN) semigroup, when the domain Ω is the ball B and the conductivity matrix, which will be defined later, is the identity matrix. To see this let us define the DtN semigroup in its whole generality as in [START_REF] Escher | The Dirichlet-Neumann operator on continuous functions[END_REF]. We consider a bounded smooth domain Ω ⊂ R n , and γ an n × n symmetric matrix with smooth (for the simplicity of the exposition we take γ ∈ C ∞ (Ω)) real elements in Ω which has uniformly bounded positive eigenvalues, i.e., there exists 0 < c 1 < c 2 such that for every x ∈ Ω and ξ ∈ R n ,

c 1 ξ 2 ≤ ξ T γ(x)ξ ≤ c 2 ξ 2 .
This matrix is known as the electrical conductivity. Let X := L 2 (Ω) or C(Ω) and let the corresponding boundary space be ∂X := C(∂Ω) or L 2 (∂Ω). We solve the following Dirichlet problem

∇ • (γ∇u) = 0, in Ω, u = f, on ∂Ω. (1.4)
For any f ∈ ∂X, we write u = L γ f . Such a function is called the γ-harmonic lifting of f and L γ the γ-harmonic lifting operator. The function u represents the electrical potential where this potential on the surface of the substance is f and the substance is in electrical equilibrium. Now define the action of the generalized DtN operator on f as the normal outward derivative of u on the boundary, i.e.,

Λ γ f := (ν.γ∇u)| ∂Ω .
In other words, we define Λ γ := ν(y) • γ∇L γ where ν(y) is the unit outer normal vector at y ∈ ∂Ω, and Λ γ is called the DtN operator. In the simplest case where γ is the identity matrix, we denote these operators by L 1 , the harmonic lifting operator, and Λ 1 , the corresponding DtN operator.

We consider Λ γ as an unbounded operator on ∂X with the domain,

D(Λ γ ) = {f ∈ ∂X : Λ γ f ∈ ∂X}.
See [START_REF] Arendt | The Dirichlet-to-Neumann operator on rough domains[END_REF][START_REF] Engel | The Laplacian on C(Ω) with generalized Wenzell boudary conditions[END_REF][START_REF] Vrabie | C 0 -Semigroups and Applications[END_REF] for various properties of this domain. The operator -Λ γ generates an analytic compact semigroup of contractions on these spaces, C(∂Ω) and L 2 (∂Ω). This semigroup can be identified as the trace of the solution to the following problem with dynamical boundary conditions:

     ∇ • (γ∇u(t, •)) = 0, for every t ∈ R + , in Ω, ∂ t u + ν • γ∇u = 0, for every t ∈ R + , on ∂Ω, u(0, •) = f, on ∂Ω.
(1.5)

In fact, by taking the trace of this solution (see [START_REF] Arendt | The Dirichlet-to-Neumann operator on rough domains[END_REF] for uniqueness) and denoting it by u(t, x) ∂Ω , we can also define the DtN semigroup as, e -tΛγ f := u(t, x) ∂Ω , for every f ∈ ∂X.

It is easy to see that the L 2 (∂Ω) version of Λ γ is selfadjoint and nonpositive, and hence Λ γ generates also an analytic semigroup of maximal angle π/2. Further properties such as contraction, compactness, positivity, irreducibility and Markov character of e -tΛγ can be found in [START_REF] Vrabie | C 0 -Semigroups and Applications[END_REF] and [START_REF] Emamirad | An approximating family for the Dirichelet-to-Neumann semigroup[END_REF]. In [START_REF] Escher | The Dirichlet-Neumann operator on continuous functions[END_REF], J. Escher showed that the operator Λ γ generates also in C(∂Ω) an analytic semigroup of some positive angle θ.

In the sequel, we are particularly interested in the case where γ is the identity matrix I n , so u will satisfy the following system.

     ∆u(t, •) = 0, for every t ∈ R + , in Ω, ∂ t u + ν • ∇u = 0, for every t ∈ R + , on ∂Ω, u(0, •) = f, on ∂Ω. (1.6)
In this case K.-J. Engel [START_REF] Engel | The Laplacian on C(Ω) with generalized Wenzell boudary conditions[END_REF] showed that θ = π/2 even for the C(Ω) version of Λ γ . This proves that for Ω = B, we get (1.1) and in this case the DtN semigroup is the one given by (1.2). The main advantage of the Lax representation is that it is only necessary to solve the problem (1.1) in order to calculate the action of the semigroup in all times. In [5], it is shown that the Lax representation is optimal in the sense that for γ = I n , the domain Ω should be a ball in R n . The first aim of this paper is to illustrate numerically the action of the Lax semigroup and perhaps this permits to have an overview of the above open question. For achieving this, we used the finite element method for the discretization of (1.5). The program is performed by FreeFem++ language (see [START_REF] Hecht | A finite element software for PDE : FreeFem++[END_REF]). This discretization is accomplished by using the implicit Euler scheme for the general case; i.e. γ is not a constant multiple of the identity matrix but a function β times the identity matrix. Thus, we have also to resolve the problem of accuracy. Since for β = 1 in whole B we have an explicit representation of the semigroup, it is convenient to consider a net β α which converges to 1 as α → 0 and to check the numerical accuracy.

In [START_REF] Cornean | Towards a D-bar reconstruction method for threedimensional eit[END_REF], it is written that, although the following result may be considered wellknown, due to lack of a suitable reference we present a full proof for the reader's convenience. In this result the authors show that if the conductivity γ = βI n with β = 1 near the boundary then one can estimate Λ γ -Λ 1 L(H 1/2 (∂Ω),H s (∂Ω)) with respect to β -1 L ∞ (Ω) . In the next section we cite again this result and we use it in order to compare Λ 1 with the general Dirichlet-to-Neumann operator with β = 1 near the boundary; this estimation implies a strong estimation in ∂X norm for the difference between e -tΛγ and e -tΛ 1 which directly depends on β -1 L ∞ (Ω) . As a consequence of this result, theoretically, the closer β is to 1, the closer e -tΛγ is to e -tΛ 1 . This shows the geometrical effects of the conductivity on the structure of the Dirichlet-to-Neumann semigroup.

The third section is divided into three parts. The first part is devoted to the Euler implicit scheme for our problem. To prove the convergence we construct an approximating family {V (t)} t≥0 and we use the Chernoff's product formula in order to theoretically prove the convergence in L 2 (∂Ω). Contrary to the explicit formulation (see [5, section 4]), here we do not have an explicit representation of V (t). So we explicitly construct its left-inverse W (t) and we define V (t) as the right-inverse of W (t). After proving the convergence, we use the finite element method via its variational formulation in order to numerically solve this problem. The second subsection is devoted to the explicit Euler scheme and we show why this method cannot be resolved by its variational formulation. Finally, in section 3.3, we consider an initial-boundary condition example which has a known harmonic lifting and we construct via Lax representation the solution at time t = 0.5. Then we define a conductivity function γ α = β α I n such that β α → 1, when α → 0 and we numerically illustrate that in this case e -tΛγ → e -tΛ 1 .

Comparison of the semigroups.

In this section we show that if γ = βI n and β = 1 in a neighborhood of ∂Ω, then for 0 ≤ t ≤ T and any f ∈ H 1/2 (∂Ω) we have

e -tΛγ f -e -tΛ 1 f L 2 (∂Ω) ≤ C T β -1 L ∞ (Ω) ,
where the constant C T depends on T , f and the constant C in (2.1). The proof of this estimation is a straightforward consequence of the following lemma which has been proved in [START_REF] Cornean | Towards a D-bar reconstruction method for threedimensional eit[END_REF].

Lemma 2.1. Let β ∈ C ∞ (Ω) be a positive conductivity satisfying β = 1 near ∂Ω, then for all s ∈ R, the following estimate holds

Λ γ -Λ 1 L(H 1/2 (∂Ω),H s (∂Ω)) ≤ C β -1 L ∞ (Ω) , (2.1)
where the constant C depends on s, Ω, the ellipticity constant of γ and the support of β -1.

Now we are able to prove the main theorem of this section.

Theorem 2.1. For β = 1 near ∂Ω such that β ∈ C ∞ (Ω) and f ∈ H 1/2 (∂Ω), there exists a constant C T depending on the ellipticity constant of γ, f and T such that

e -tΛγ f -e -tΛ 1 f L 2 (∂Ω) ≤ C T β -1 L ∞ (Ω) , for all t ∈ [0, T ].
Proof. One can readily write

e -tΛγ f -e -tΛ 1 f = t 0 d ds (e (s-t)Λ 1 e -sΛγ f )ds = t 0 e (s-t)Λ 1 (Λ 1 -Λ γ )e -sΛγ f

ds

Since the semigroups e -tΛγ and e -tΛ 1 are the analytic and contraction semigroups, then for any f ∈ H 1/2 (∂Ω) we have e -tΛγ f ∈ H 1/2 (∂Ω) and we get :

e -tΛγ f -e -tΛ 1 f L 2 (∂Ω) ≤ t 0 (Λ 1 -Λ γ )e -sΛγ f L 2 (∂Ω) ds ≤ C T β -1 L ∞ (Ω)
where C T = CM T , with

M T = sup t∈[0,T ] t 0 e -sΛγ f H 1/2 (∂Ω) ds.

Numerical results

Here we are going to approximate the Dirichlet-to-Neumann semigroup by means of Chernoff's Theorem. Let us recall this Theorem which is proved in [START_REF] Chernoff | Note on product formulas for operator semigroups[END_REF]. Theorem 3.1 (Chernoff's product formula). Let X be a Banach space and {V (t)} t≥0 be a family of contractions on X with V (0) = I. Suppose that the derivative V ′ (0)f exists for all f in a set D and that the closure Λ of V ′ (0) D generates a (C 0 ) semigroup S(t) of contractions. Then, for each f ∈ X,

lim n→∞ V t n n f = S(t)f, uniformly for t in compact subsets of R + .
This procedure was achieved in [START_REF] Emamirad | An approximating family for the Dirichelet-to-Neumann semigroup[END_REF] by choosing the approximating family [V (t)f ](x) = u(e -tγ(x) x), where |x| = 1 and Ω = {x ∈ R n : |x| < 1}. In [5] the approximation family {V (t)} t≥0 is chosen for the Euler explicit scheme, with

V (t)f (x) = (1 -α)u(x) + αu(x -α -1 tγ(x)ν(x)), 0 ≤ t ≤ αT, V (αT )f (x), t > αT, (3.1)
where u = L γ f and x ∈ ∂Ω. By using Theorem 3.1, it is proved that the following Euler explicit scheme

(EES)      div(γ∇u m+1 ) = 0 in Ω, 1 ∆t (u m+1 -u m ) + ∂u m ∂νγ = 0 on ∂Ω, u 0 = f on ∂Ω.
converges to the solution of the problem (1.5). For the Euler implicit scheme one can provide a variational formulation programmable in finite element method of the problem.

For this purpose we replace the time derivative by

∂u ∂t (m∆t, x) ≃ 1 ∆t u m+1 (x) -u m (x) , for x ∈ ∂Ω, and m = 0, • • • , [T /∆t],
where u m (x) = u(m∆t, x). So, the Euler implicit scheme for the discrete boundary dynamical system can be written as

(EIS)      div(γ∇u m+1 ) = 0 in Ω, 1 ∆t (u m+1 -u m ) + ∂u m+1 ∂νγ = 0 on ∂Ω, u 0 = f on ∂Ω.
Since any x with |x| = 1 belongs to ∂Ω, we have

∂u m+1 ∂ν γ ≈ u m+1 (x) -u m+1 (x -∆xγ(x)x) ∆x . (3.2)
By replacing (3.2) in (EIS), we get

1 + ∆t ∆x u m+1 (x) - ∆t ∆x u m+1 (x -∆xγ(x)x) = u m (x). (3.3)
Now let 0 < α ≤ 1 be a parameter. Since β(x) ≥ C -1 everywhere x ∈ ∂Ω, there exists T > 0 such that for all t ≤ T , we have x -tγ(x)ν(x) ∈ Ω. We may construct a family of operators in L (L 2 (∂Ω)) as follows. Set

W (t)f (x) = (1 + α)u(x) -αu(x -α -1 tγ(x)ν(x)), 0 ≤ t ≤ αT, W (αT )f (x), t > αT, (3.4)
where u is the γ-harmonic lifting of f . By taking α = ∆t ∆x we can write (3.3) as W (∆t)u m+1 (x) = u m (x).

(3.5)

Now we define the family V (t) as the inverse of W (t) in the following sense

W (t)V (t)f (x) = f (x), (3.6)
and we show that V (t) satisfies all the assumptions of Chernoff's product formula in X = L 2 (∂Ω) with . = . L 2 (∂Ω) .

(1) V (0)f = f . Since for any x ∈ ∂Ω, we have

W (0)f (x) = u(x) = f (x), (3.6)
shows

f (x) = W (0)V (0)f (x) = V (0)f (x).
(2) V ′ (0)f = -Λ γ f . The derivative of W (t) with respect to t is

W ′ (t)f = ∇u(x -α -1 tγ(x)x)γ(x)x
and at point t = 0 we have W ′ (0)f = Λ γ f . Furthermore, the derivative of (3.6) gives W ′ (t)V (t) + W (t)V ′ (t) = 0 and its evaluation at t = 0, yields

Λ γ f + V ′ (0)f = 0. (3) V (t) is of contraction. Let us first retrieve a lower bound for W (t) in L 2 (∂Ω). In fact W (t)f 2 = ∂Ω u(x) + α u(x) -u(x -α -1 tγ(x)x) 2 dσ. Since (u(x) -u(x -α -1 tγ(x)x)) /α -1
is an approximation of the normal derivative ∂u/∂ν γ , so according to Hopf's Lemma (see [START_REF] Evans | Partial Differential Equations Amer[END_REF]), this term is positive and consequently

∂Ω u(x) + α u(x) -u(x -α -1 tγ(x)x) 2 dσ ≥ ∂Ω |u(x)| 2 dσ
which implies that W (t)f ≥ f . This implies that f = W (t)V (t)f ≥ V (t)f . So we have proved that V (t) satisfies all the assumptions of Chernoff's product formula and consequently the implicit Euler scheme converges to its exact solution.

The variational formulation of this problem can be obtained by multiplying both sides of the dynamic boundary condition by a test function v and by using the divergence theorem, we get

Ω γ∇u m+1 ∇vdx - ∂Ω γ ∂u m+1 ∂n vdσ = 0, that is Ω ∆tγ∇u m+1 ∇vdx + ∂Ω u m+1 v - ∂Ω u m vdσ = 0, (3.7)
which is of the form

a(u m+1 , v) = ℓ(v), where a(u m+1 , v) = Ω ∆tγ∇u m+1 ∇vdx + ∂Ω u m+1 vdσ
is the bilinear form with the unknown of the problem u m+1 and

ℓ(v) = ∂Ω u m vdσ
is a linear form. This variational formulation can be resolved by finite elements implemented in FreeFem++ and the result is illustrated in the following subsection.

3.1. Numerical illustration. We begin by taking the initial function h(x, y) = x 2 + y on the boundary. The Dirichlet problem

∆u = 0, in Ω u = h,
on ∂Ω admits a unique solution u(x, y) = h(x, y) -g(x, y), where g(x, y) = (x 2 + y 2 -1)/2. Since g = 0 on the boundary of Ω, ∆h = 2 and ∆g = 2. So the solution will be u(x, y) = (x 2 -y 2 )/2 + y + 1/2.

(3.8) The peculiarity of this function is that on the straight lines y = -x + 1 for (x, y) ∈ R + × R + and y = x + 1 for (x, y) ∈ R -× R + one has u = 1 (see Figure 1). Now we include the potential γ = 1. Define γ α as γ α (x, y) = 1 + α exp((-9 16 + x 2 + y 2 ) -1 ) if x 2 + y 2 < 9 16 , 1 if x 2 + y 2 ≥ 9 16 . We remark that 1 -γ α has its support in the disc D := {z ∈ C : |z| ≤ 3/4}, furthermore 1 -γ α L ∞ (Ω) = αe -16/9 , hence in order to justify our Theorem 2.1, we have to move α towards zero and observe the convergence of the solution to the exact solution given by (3.8). This is done in the subsequent Figures. Figure 4 to Figure 7 respectively illustrate the solution for α =300, 50, 10 and 0.1.

The calculations are performed via variational formulations described as above. We take 50 points on the circle ∂Ω and the domain inside is triangulated with Lagrange P 1 finite elements. The visualization is drawn indicating the isovalues of the finite element solution. We note the similarity of the exact solution (Figure 2) with the finite element solution for α = 0.1 (Figure 7). To illustrate this similarity we have computed at the end the difference between this solution and the exact solution in Figure 8. The exact solution at time t = 0.5 u(0.5, x, y) = ((exp(-0.5)x) 2 -((exp(-0.5)y) 2 )/2 + exp(-0.5)y + 1/2. The exact solution at time

IsoValue -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.
t = 2 u(2, x, y) = ((exp(-2)x) 2 -((exp(-2)y) 2 )/2 + exp(-2)y + 1/2.
Finally for illustrating the first part of this paper we have performed a movie with mp4 extension which is available online in http://www.math.univ-metz.fr/ jmse/Emamirad

Here we have taken the boundary function f (x, y) = x 4 + y 2 sin(2πy), and 30 seconds after starting the movie res.mp4, we observe that we make a close-up of the center and that we approach the values of u( 0

Figure 2 .

 2 Figure 2.The exact solution at time t = 0.5 u(0.5, x, y) = ((exp(-0.5)x) 2 -((exp(-0.5)y) 2 )/2 + exp(-0.5)y + 1/2.

Figure 3 .

 3 Figure 3.The exact solution at time t = 2 u(2, x, y) = ((exp(-2)x) 2 -((exp(-2)y) 2 )/2 + exp(-2)y + 1/2.

  cos4 θ + sin 2 θ sin(2π sin θ) dθ = 0.375 Figure 6. Solution for α = 10 Figure 7. Solution for α = 0.1 Error between u exa and u 0.1

	) =	1 2π	0	2π
	U_300(0.5,x,y)			IsoValue -0.0151183 0.0419633 0.0800177 0.118072 0.156126 0.194181 0.232235 0.270289 0.308344 0.346398 0.384453 0.422507 0.460561 0.498616 0.53667 0.574724 0.612779 0.650833 0.688887 0.784023	U_50(0.5,x,y)	0.868492 0.746339 0.697478 0.648617 0.599755 0.550894 0.502033 0.453172 0.404311 0.355449 0.306588 0.257727 0.208866 0.160004 0.111143 0.062282 0.0134208 -0.0354404 -0.0843016 -0.157593 IsoValue
	Figure 4. Solution for α = 300	Figure 5. Solution for α = 50
	U_10(0.5,x,y)			IsoValue -0.277674 -0.19028 -0.132017 -0.0737541 -0.0154912 0.0427718 0.101035 0.159298 0.21756 0.275823 0.334086 0.392349 0.450612 0.508875 0.567138 0.625401 0.683664 0.741927 0.80019 0.945847	U_0.1(0.5,x,y)	0.997316 0.836587 0.772295 0.708003 0.643711 0.57942 0.515128 0.450836 0.386544 0.322252 0.257961 0.193669 0.129377 0.0650853 0.000793497 -0.0634983 -0.12779 -0.192082 -0.256374 -0.352811 IsoValue
		erreur			0.00233312 0.0018511 0.00165829 0.00146549 0.00127268 0.00107987 0.000887062 0.000694254 0.000501446 0.000308637 0.000115829 -7.69788e-005 -0.000269787 -0.000462595 -0.000655403 -0.000848211 -0.00104102 -0.00123383 -0.00142664 -0.00171585 IsoValue
	Figure 8.


Salwa Mansour
email: salwa.mansour@inria.fr

Mohamad Muhieddine
email: mohamad.muhieddine@liu.edu.lb

Édouard Canot
email: edouard.canot@irisa.fr

Ramiro March
email: ramiro.march@univ-rennes1.fr

Identification of Thermophysical Properties of the Soil in 3D-axisymmetric Coordinate System Using Inverse Problem

Keywords: Inverse Problem, Saturated Porous Medium, Phase Change, Agricultural and Archaeological Soils, Gauss Newton Method

This paper is motivated by the studies of agricultural and archaeological soils. We introduce a numerical strategy in 3D axisymmetric coordinate system to estimate the thermophysical properties of a saturated porous medium (volumetric heat capacity (𝜌𝐶) 𝑠 thermal conductivity 𝜆 𝑠 and porosity 𝜙 ) where a phase change problem (liquid/vapor) appears due to strong heating. The estimation of these thermophysical properties is done by inverse problem knowing the heating curves at selected points of the porous medium. In order to solve the inverse problem, we use the least square criterion in which sensitivity coefficients appear and which leads to a system of ordinary differential equations (ODE). At the stage of numerical computations, we propose a global approach based on the method of lines, where time and space discretizations are considered separately. The space discretization is done using a vertex-centered finite volume method; the discretization in time is done via an ODE solver combined with a modified Newton's method to deal with the nonlinearity of the system of coupled equations.

Introduction

The work presented in this paper is motivated by the studies of agricultural and archaeological soils. A systematic application of numerical modeling in a particular field of agriculture and archaeology which is the study of seed germination and archaeological hearths is presented. The authors introduce a numerical strategy in 3D axisymmetric coordinate system to estimate the thermophysical properties of the soil (volumetric heat capacity(𝜌𝐶) 𝑠 , thermal conductivity 𝜆 𝑠 and porosity 𝜙 ) of a saturated porous medium where a phase change problem (liquid/vapor) appears due to intense heating from above. The advantage of our 3D-axisymmetry configuration is that it is well adapted to the heating of a semi-infinite medium by a circular heating plate, which is a very simple device. Usually ϕ is the porosity, however when the soil is not saturated (which should concern most cases), ϕ may be taken equal to the part of water in the pores. This is of course an approximation which is correct for the energy balance but which neglects the capillary forces and the migration flow of the liquid inside the porous media; a complete model of such an unsaturated model is out of the scope of this paper. *Salwa Mansour. Tel : +33 6 99 84 25 11 E-mail: salwa.mansour@inria.fr

The investigation of the thermal properties of the soil can have significant practical consequences such as evaluation of optimum conditions for plant growth and development and can be utilized for the control of thermal-moisture regime of soil in the field [START_REF] Oladunjoye | Variability of soil thermal properties of a seasonally cultivated agricultural teaching and research farm, university of Ibadan, south-western Nigeria[END_REF]. These properties influence how energy is partitioned in the soil profile so the ability to monitor them is a tool to manage the soil temperature regime that affects seed germination and growth. It can also provide information about the use of fire by ancient civilizations whether for cooking or heating ... [START_REF] Eckmeier | Detecting prehistoric fire-based farming using biogeochemical markers[END_REF].

The inverse problem, presented in this paper, consists of the estimation of thermophysical properties of the soil knowing the heating curves at selected points of the altered soil [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identification of the thermophysical properties of the soil by inverse problem[END_REF]. In general, the mathematical formulation of inverse problems leads to models that are typically ill-posed. According to Hadamard, a mathematical problem is ill-posed if one of the following properties is violated: o For all admissible data, a solution exists. o For all admissible data, the solution is unique and depends continuously on the data [START_REF] Engl | Nonlinear inverse problems: theoretical aspects and some industrial applicatIons[END_REF].

In such ill-posed inverse problems, we usually minimize a discrepancy between some experimental data and some model data [START_REF] Favennec | Lecture 7: Optimization methods for nonlinear estimation and function estimation[END_REF]. In our problem, we use the least square criterion in which the sensitivity coefficients appear and where we try to minimize the discrepancy function which is expressed as the norm of the difference between the experimental temperature and the synthetic ones [START_REF] Bjorck | Numerical methods for least squares problems[END_REF].

The system composed of the energy equation together with three boundary initial problems resulting from differentiating the basic energy equation with respect to the soil properties must be solved [START_REF] Majchrzak | Identification of substitute thermal capacity of solidifying alloy[END_REF].

At the stage of numerical computations, the Gauss Newton method is used to minimize the least square criterion; that requires the solution of a system of four highly nonlinear ordinary differential equations. We propose a global approach similar to that presented by [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identification of the thermophysical properties of the soil by inverse problem[END_REF]. Basically, our paper generalizes the 1D inverse problem presented by [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identification of the thermophysical properties of the soil by inverse problem[END_REF] to the 3D-axisymmetry configuration. It is important to note that in this new configuration, the solution was reached after taking into consideration the temperature history at different time steps which was not the case in [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identification of the thermophysical properties of the soil by inverse problem[END_REF] where the authors reached the solution by taking the temperature history at the final time only. This approach is based on the method of lines, where time and space discretizations are considered separately. The space discretization is done using a vertex-centered finite volume method; the discretization in time is done via an ODE solver that uses a BDF method and a modified Newton method to deal with the high nonlinearity. The code validation stage is based on the comparison between the numerical results and the synthetic data.

Mathematical formulation

Consider a porous medium in 3D axisymmetric coordinate system as shown in Fig. 1. 𝛺 is the bounded domain with boundary 𝑑𝛺 = 𝛤 𝐷 ∪ 𝛤 𝑁 . The medium is heated from above with 𝑇 𝑐 (temperature of the fire) greater than the 𝑇 𝑣 (evaporation temperature which is normally 100℃).

Resolution of heat diffusion problem by apparent heat capacity method.

In order to model the heat conduction transfer in the soil, we use the energy equation. We neglect the convection term of the heat equation so that the energy conservation equation for the unknown temperature 𝑇 is expressed as:

(
where 𝑇 represents the temperature, 𝑇 0 is the initial temperature at 𝑡 0 = 0 , 𝑇 𝐷 is 𝑇 𝑐 at the fire and 𝑇 0 elsewhere, 𝜌 is the density, 𝐶 is the specific heat capacity, 𝜆 is the conductivity, 𝜙 is the porosity, the subscripts 𝑒, 𝑓 and 𝑠 indicate the equivalent parameters of the medium, the properties of the fluid and the porous matrix properties respectively.

The effective volumetric heat capacity and the effective conductivity are defined by: (𝜌𝐶) 𝑒 = 𝜙(𝜌𝐶) 𝑓 + (1 -𝜙)(𝜌𝐶) 𝑠

(2) and

1 𝜆 𝑒 = 𝜙 𝜆 𝑓 + 1-𝜙 𝜆 𝑠 (3)
The effective conductivity in equation ( 3) is calculated using the harmonic mean to test the algorithm. In real situation, the harmonic mean should be replaced by some approximated model, such as this of Kunii & Smith [START_REF] Malherbe | Measurement of thermal conductivity of granular materials over a wide range of temperatures. Comparison with theoretical models[END_REF].

To avoid the tracking of the interface of the phase change problem (liquid/vapor) which appears when the water existing in the soil turns into gas, the apparent heat capacity (AHC) method is used because it allows a continuous treatment of a system involving phase transfer. In this method [START_REF] Bonacina | Numerical solution of phase change problems[END_REF][START_REF] Muhieddine | Various approaches for solving problems in heat conduction with phase change[END_REF], the latent heat is taken into account by integrating the heat capacity over the temperature, and the domain is considered to be treated as one region. The singularity in the AHC method which is presented in the formulation of thermo-physical properties defined by [START_REF] Bonacina | Numerical solution of phase change problems[END_REF] is treated in [START_REF] Muhieddine | Various approaches for solving problems in heat conduction with phase change[END_REF]. The Dirac delta function, representing the equivalent heat capacity, (see Fig. 2 continuous lines) can be approximated by the normal distribution:

𝑑𝜎 𝑑𝑇 = (𝜀𝜋 -1 2)exp [-𝜀 2 (𝑇 -𝑇 𝑣 ) 2 ] ( 4
)
where 𝜀 = √2 ∆𝑇

, ∆𝑇 is the phase change interval and 𝑇 𝑣 is the phase change temperature. The integration of equation ( 4) yields the error functions approximations for the initial phase fraction:

The smoothed coefficients (see Fig. 2 dashed lines) could be written as:

𝐶 𝑓 = 𝐶 𝑙 + (𝐶 𝑣 -𝐶 𝑙 )𝜎 + 𝐿 𝑑𝜎 𝑑𝑇 (6) 𝜆 𝑓 = 𝜆 𝑙 + (𝜆 𝑣 -𝜆 𝑙 )𝜎 (7) 𝜌 𝑓 = 𝜌 𝑙 + (𝜌 𝑣 -𝜌 𝑙 )𝜎 (8)
where L is the latent heat of phase change and the subscripts 𝑙 and 𝑣 indicate respectively the properties of the liquid phase (water) and that of the vapor. It is important to mention that the thermophysical properties of the fluid are temperature dependent and that is why the problem is highly nonlinear.

Inverse Problem

In order to solve the parametric inverse problem consisting of finding the volumetric heat capacity (𝜌𝐶) 𝑠 , the conductivity 𝜆 𝑠 and the porosity 𝜙 of the saturated soil, it is necessary to know the values of temperatures 𝑇 𝑔𝑖 𝑓 at selected points 𝑥 𝑖 of the porous medium domain for times 𝑡 𝑓 : 𝑇 𝑔𝑖 𝑓 = 𝑇 𝑔 (𝑥 𝑖 , 𝑡 𝑓 ) where 𝑖 = 1, 2, … , 𝑀 and 𝑓 = 1, 2, … , 𝐹.

We use the least squares criterion to solve this inverse problem so we try to find the soil parameters that minimize the error function 𝑆 which is defined by:

𝑆((𝜌𝐶) 𝑠 , 𝜙, 𝜆 𝑠 ) = 1 2 ∑ ∑ (𝑇 𝑖 𝑓 -𝑇 𝑔𝑖 𝑓 ) 2 𝐹 𝑓=1 𝑀 𝑖=1 (9)
where 𝑇 𝑖 𝑓 = 𝑇(𝑥 𝑖 , 𝑡 𝑓 ) are the temperatures being the solution of the direct problem for the assumed set of parameters at the points 𝑥 𝑖 , 𝑖 = 1, 2, … , 𝑀 for the time 𝑡 𝑓 , 𝑓 = 1, 2, … , 𝐹 and 𝑇 𝑔𝑖 𝑓 are the measured temperatures at the same points 𝑥 𝑖 for time 𝑡 𝑓 .

Method of resolution

To illustrate the method of resolution, we define the following vectors: where 𝑖 = 1, 2, … , 𝑁 and 𝑗 = 1, 2, 3. The error function 𝑆 defined by equation ( 9) can be written as

𝒑 (𝑘) = ( (
𝑆(𝑝 (𝑘) ) = 1 2
𝑟(𝑝 (𝑘) ) 𝑡 𝑟(𝑝 (𝑘) ) and its first derivative is given by: ∇𝑆(𝑝 (𝑘) ) = 𝐽(𝑝 (𝑘) ) 𝑡 𝑟(𝑝 (𝑘) ). A necessary condition for 𝑝 * to be a local minimum of 𝑆(𝑝 (𝑘) ) is that ∇𝑆(𝑝 * ) = 0 but it is not sufficient. It is required to show that the second derivative of 𝑆(𝑝 (𝑘) ) is positive.

For this reason, we choose to solve our non-linear least square problem by Gauss-Newton method because it doesn't require the computation of the second derivative and has advantage to solve equation ( 8) by a small number of iterations and to converge quickly toward the local minimum of the error function 𝑆. It proceeds by the iterations:

𝑝 (𝑘+1) = 𝑝 (𝑘) -[𝐽(𝑝 (𝑘) ) 𝑡 𝐽(𝑝 (𝑘) )] -1 𝐽(𝑝 (𝑘) ) 𝑡 𝑟(𝑝 (𝑘) ) (10)
𝐽(𝑝 (𝑘) ) is given by:

𝐽(𝑝 (𝑘) ) = ( 𝑊 1 1(𝑘) ⋮ 𝑊 1 𝐹(𝑘) 𝑊 2 1(𝑘) ⋮ 𝑊 2 𝐹(𝑘) ⋮ 𝑊 𝑀 1(𝑘) ⋮ 𝑊 𝑀 𝐹(𝑘) 𝑍 1 1(𝑘) ⋮ 𝑍 1 𝐹(𝑘) 𝑍 2 1(𝑘) ⋮ 𝑍 2 𝐹(𝑘) ⋮ 𝑍 𝑀 1(𝑘) ⋮ 𝑍 𝑀 𝐹(𝑘) 𝑅 1 1(𝑘) ⋮ 𝑅 1 𝐹(𝑘) 𝑅 2 1(𝑘) ⋮ 𝑅 2 𝐹(𝑘) ⋮ 𝑅 𝑀 1(𝑘) ⋮ 𝑅 𝑀 𝐹(𝑘) )

Governed Equations

In the following, we present the heat equation ( 1) in 3Daxisymmetric coordinate system (see figure 3) together with 3 sensitivity equations resulting from the differentiation of the heat equation ( 11) with respect to the soil parameters.

The heat diffusion equation in the 3Daxisymmetric coordinate system is given by: The three sensitivity equations are obtained by differentiating the heat diffusion equation with respect to the unknown parameters (𝜌𝐶) 𝑠 , 𝜙 and 𝜆 𝑠 respectively. In order to determine the sensitivity coefficients (𝑊, 𝑍 and 𝑅) appearing in the Jacobian, we should solve the three sensitivity equations (12), (15) and (18) :

3 Numerical Strategy

The obtained system of coupled equations (sensitivity equations + heat diffusion equation) is a nonlinear system of ordinary partial differential equations. To solve this system, we propose a global approach similar to that presented by [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identification of the thermophysical properties of the soil by inverse problem[END_REF]. This approach is based on the method of lines, where time and space discretizations are considered separately. The space discretization is done using a vertex-centered finite volume method; the discretization in time is done by an Euler implicit scheme. The system of coupled equations ( 11), ( 12), ( 15) and ( 18) is discretized in space using the vertex-centered finite volume method in 3D-axisymmetric coordinate system. The computational domain is divided into a grid or mesh of dimensions ∆𝑟 and ∆𝑧 as shown in Fig. 3. After spatial discretization, the system of coupled equations can be written as a linear system:

( 𝑀(𝑇) 𝐴 𝑊 (𝑇) 𝐴 𝑍 (𝑇) 𝐴 𝑅 (𝑇) 0 𝐼 0 0 0 0 𝐼 0 0 0 0 𝐼 ) ( 𝑑𝑇 𝑑𝑡 𝑑𝑊 𝑑𝑡 𝑑𝑍 𝑑𝑡 𝑑𝑅 𝑑𝑡 ) + ( 𝑁(𝑇) 0 0 0 0 𝐵 𝑊 (𝑇) 0 0 0 0 𝐵 𝑍 (𝑇) 0 0 0 0 𝐵 𝑅 (𝑇) ) ( 𝑇 𝑊 𝑍 𝑅 ) = 𝑏 (21)
which is an ordinary differential equation. By classical transformations and by letting 𝑌 = ( 𝑇 𝑊 𝑍 𝑅 ) , the system (21) can be written in the general form:

{ 𝑌 ′ = 𝐹(𝑇, 𝑌) 𝑌(𝑡 0 ) = 𝑌 0 (22)
We use an ODE solver to solve this system, this solver applies a time discretization scheme using Backward Differentiation Formula (BDF) because at each time step a system of nonlinear equations using a modified Newton method must be solved.

Algorithm

The aim of the inverse problem is the calculation of the parameters' vector 𝑝 which makes the values of the temperature 𝑇(𝑡), calculated by the forward problem, approach the experimental temperature 𝑇 𝑔 (𝑡). The algorithm is as follows:

1. Choose an initial value 𝑝 (0) of 𝑝 at iteration k=0.

2. Solve the system (21) formed of the coupled equations using 𝑝 (𝑘) to define the properties of the soil. 3. Deduce the values of 𝑇 𝑖 𝑓,(𝑘) ,𝑊 𝑖 𝑓,(𝑘) , 𝑍 𝑖 𝑓,(𝑘) and 𝑅 𝑖 𝑓,(𝑘) for 𝑖 = 1, 2, … , 𝑀 and 𝑓 = 1 ,2, … , 𝐹.

4. Calculate the value of 𝑟(𝑝 (𝑘) ) and that of the Jacobian 𝐽(𝑝 (𝑘) ). 5. Solve the linear system 𝑝 (𝑘+1) = 𝑝 (𝑘) -[𝐽(𝑝 (𝑘) )

𝑡 𝐽(𝑝 (𝑘) )] -1

𝐽(𝑝 (𝑘) ) 𝑡 𝑟(𝑝 (𝑘) ) for 𝑝 (𝑘+1) .

6. Calculate 𝑆(𝑝 (𝑘) ). If 𝑆(𝑝 (𝑘) ) < 𝜀, end (the criterion of convergence must be determined) If not, iterate: 𝑝 (𝑘) ← 𝑝 (𝑘+1) and go to 2.

Some validation examples

The code validation stage is based on the comparison between the numerical results and the synthetic data. First, we choose a plausible example where the soil properties (volumetric heat capacity, thermal conductivity and porosity) are given constant values. These values are used by the forward problem to calculate the temperatures at different points of the domain which in return plays the role of the experimental data in the inverse problem. Table 1 shows that the numerical results are in good agreement with the exact experimental ones. From table 2, we can notice that if the initial guess of (𝜌𝐶) 𝑠 is equal to its exact value then ϕ and 𝜆 𝑠 converge to their exact values. The final results of ϕ and 𝜆 𝑠 are shown in figures 4 and 5.

To study the influence of the initial guess of (𝜌𝐶) 𝑠 on the convergence, we consider table 3. The computations show that we have relatively a good convergence for ϕ and 𝜆 𝑠 . However, the convergence of (𝜌𝐶) 𝑠 is far from the exact solution. This is might be due to the fact that the heat equation is not sensitive to (𝜌𝐶) 𝑠 in comparison to the other parameters (ϕ and 𝜆 𝑠 ). In most of the cases, (𝜌𝐶) 𝑠 stays equal to its initial guess.

Conclusion

The idea of this paper is the inverse problem which consists of the estimation of thermophysical properties of the soil knowing the temperature at selected points of the domain. In order to solve this inverse problem, we used the least square criterion where we try to minimize the error function between real measures and synthetic ones. The coupled system composed of the energy equation together with the three sensitivity boundary initial problems resulting from differentiating the basic energy equation with respect to the soil properties must be solved. To overcome the high nonlinearity of the coupled system, we used Gauss Newton's method. According to Gauss Newton's method, the convergence rate should be linear [START_REF] Dennis | Numerical methods for unconstrained optimization and nonlinear equations[END_REF] but we obtained poor convergence ( see tables 1 and 3) except when we removed (𝜌𝐶) 𝑠 , (see table 2). The code validation stage is based on the comparison between the numerical results and the synthetic data which shows good agreement and convergence is obtained after about 10 iterations. The rate of convergence of our inverse problem is linear which agrees with the order of convergence of Gauss Newton method.

Our problem of interest is very difficult to be applied to real data because convergence cannot be obtained for relatively small ∆𝑇. In this case, the numerical solution will be far from the physical one. It appears that it is necessary to reduce both the temperature interval ∆𝑇 and the mesh size in order to get accurate results but this will lead to a high numerical cost.

Fig. 1 . 2 ( 1 +

 121 Fig.1. The computational domain in 3D axisymmetric coordinate system.

  and the 𝑁 × 1 residual vector is defined by: 𝒓(𝑝 (𝑘) ) = 𝒈(𝒑 (𝒌) ) -𝑻 𝑔 , where 𝑁 = 𝑀 × 𝐹 and 𝑘 is the iteration number. The Jacobian of the residual vector 𝒓(𝒑(𝑘) ) is defined by 𝐽(𝒑 (𝑘) ) 𝑖,𝑗 = 𝜕𝑟 𝑖 (𝒑 (𝑘) )𝜕𝑝 𝑗 (𝑘)

Fig. 2 .

 2 Fig.2. Equivalent thermo-physical properties in the AHC method.

  with the following initial and boundary conditions: 𝑇(𝑥, 0) = 𝑇 0 (𝑥) 𝑖𝑛 𝛺 𝑇(𝑥, 𝑡) = 𝑇 𝐷 (𝑥, 𝑡) 𝑜𝑛 𝛤 𝐷 × (0, 𝑡 𝑒𝑛𝑑 ] (Dirichlet) ∇𝑇(𝑥, 𝑡). 𝜈 = 0 𝑜𝑛 𝛤 𝑁 × (0, 𝑡 𝑒𝑛𝑑 ] (Neumann).



  Differentiating with respect to (𝜌𝐶) 𝑠 ,: (𝜆 𝑣 -𝜆 𝑙 ) 𝑑𝜎 𝑑𝑇 (𝜙𝜆 𝑠 +(1-𝜙)𝜆 𝑓 )-(1-𝜙)(𝜆 𝑣 -𝜆 𝑙 ) 𝑑𝜎 𝑑𝑇 𝜆 𝑓 𝜆 𝑠 ].𝑊 [𝜙𝜆 𝑠 +(1-𝜙)𝜆 𝑓 ] ) 𝑓 -(𝜌𝐶) 𝑠 + 𝜙 [𝜌 𝑓 [(𝑐 𝑣 -𝑐 𝑙 ) (𝜆 𝑣 -𝜆 𝑙 ) 𝑑𝜎 𝑑𝑇 𝑍(𝜙𝜆 𝑠 +(1-𝜙)𝜆 𝑓 )-(𝜆 𝑠 -𝜆 𝑓 +(1-𝜙)(𝜆 𝑣 -𝜆 𝑙 ) 𝑑𝜎 𝑑𝑇 𝑍)𝜆 𝑓 𝜆 𝑠 ] [𝜙𝜆 𝑠 +(1-𝜙)𝜆 𝑓 ] Differentiating with respect to 𝜆 𝑠 : 𝑠 (𝜆 𝑣 -𝜆 𝑙 ) 𝑑𝜎 𝑑𝑇 𝑅+𝜆 𝑓 )(𝜙𝜆 𝑠 +(1-𝜙)𝜆 𝑓 )-(𝜙+(1-𝜙)(𝜆 𝑣 -𝜆 𝑙 ) equations (12), (15) and (18) are completed with adequate initial and boundary conditions analogous to those in equation[START_REF] Dennis | Numerical methods for unconstrained optimization and nonlinear equations[END_REF] but in homogenous form. 𝑊, 𝑍 and 𝑅 are the unknowns of the

Fig. 3 .

 3 Fig.3. Discretization in 3D-axisymmetric coordinate system

Fig. 4 .

 4 Fig.4. Convergence of the porosity of the soil with respect to the number of iterations. The dashed red line represents the desired solution.

Fig. 5 .

 5 Fig.5. Convergence of the thermal conductivity of the soil with respect to the number of iterations. The dashed red line represents the desired solution.

Table 1 . Properties of the soil obtained by inverse problem

 1

		(𝝆𝑪) 𝒔	ϕ	𝝀 𝒔
	Exact	𝟏. 𝟗𝟓 × 𝟏𝟎 𝟔	0.20	0.756
	Initial guess	𝟐 × 𝟏𝟎 𝟔	0.12	0.4
	Calculated	𝟏. 𝟗𝟗 × 𝟏𝟎 𝟔	0.199	0.770

Table 2 . Properties of the soil obtained by inverse problem

 2

		(𝝆𝑪) 𝒔	ϕ	𝝀 𝒔
	Exact	𝟏. 𝟗𝟓 × 𝟏𝟎 𝟔	0.20	0.756
	Initial guess	𝟏. 𝟗𝟓 × 𝟏𝟎 𝟔	0.12	0.9
	Calculated	𝟏. 𝟗𝟓 × 𝟏𝟎 𝟔	0.20	0.756

Table 3 . Properties of the soil obtained by inverse problem

 3

		(𝝆𝑪) 𝒔	ϕ	𝝀 𝒔
	Exact	𝟏. 𝟗𝟓 × 𝟏𝟎 𝟔	0.20	0.756
	Initial guess	𝟏. 𝟕 × 𝟏𝟎 𝟔	0.12	0.9
	Calculated	𝟏. 𝟕 × 𝟏𝟎 𝟔	0.203	0.686


C Chandre
email: chandre@cpt.univ-mrs.fr

Incomplete Dirac reduction of constrained Hamiltonian systems

First-class constraints constitute a potential obstacle to the computation of a Poisson bracket in Dirac's theory of constrained Hamiltonian systems. Using the pseudoinverse instead of the inverse of the matrix defined by the Poisson brackets between the constraints, we show that a Dirac-Poisson bracket can be constructed, even if it corresponds to an incomplete reduction of the original Hamiltonian system. The uniqueness of Dirac brackets is discussed.

Introduction

We consider an N-dimensional Hamiltonian system with phase space variables z = (z 1 , z 2 , . . . , z N ) given by its Hamiltonian H(z) and its Poisson bracket {•, •}. We impose a set of K constraints Φ n (z) = 0 for n = 1, . . . , K on this dynamical system. We define the Dirac bracket as

{F, G} * = {F, G} -{F, Φ n }D nm {Φ m , G}, (1)
with implicit summation over repeated indices. Usually the matrix D whose elements are D nm is taken as the inverse of the matrix C whose elements are

C nm = {Φ n , Φ m },
if it is invertible [1,2,[START_REF] Hanson | Constrained Hamiltonian Systems[END_REF][START_REF] Sundermeyer | Constrained Dynamics[END_REF][START_REF] Henneaux | Quantization of Gauge Systems[END_REF][START_REF] Wipf | Hamilton's formalism for systems with constraints, in Canonical Gravity: From Classical to Quantum[END_REF]. In this case, it has been shown that the usual Dirac bracket (1) is a Poisson bracket [1], in particular, that it satisfies the Jacobi identity everywhere in phase space (and not just on the surface defined by the constraints). Dirac's theory of constrained Hamiltonian systems has been used in a wide variety of contexts [START_REF] Wipf | Hamilton's formalism for systems with constraints, in Canonical Gravity: From Classical to Quantum[END_REF]. Recently, it has been used to derived reduced models in fluid and plasma physics [7,8,9,10].

In the literature there are some tentative definitions of weak Dirac bracket [11,12,13,14] with the aim of defining a Dirac bracket even in situations where the matrix C is not invertible. The non-invertibility of C is linked to the existence of first class constraints, i.e., constraints that commute with all the other constraints (in the weak sense). The tentatives to define such a generalization of the Dirac bracket so far did not result in the definition of a well defined Poisson bracket, i.e., which satisfies the Jacobi identity everywhere in phase space. More specifically, in these tentatives, the Jacobi identity together with the commutation of the constraints with any function only happens on the surface defined by the constraints.

Here we generalize the Dirac bracket to cases where C is not invertible by taking D as the Moore-Penrose pseudoinverse of C. The conditions which D has to satisfy are

C = CDC, (2)
D = DCD, (3)
CD = DC. (4)
In finite dimensions, the pseudoinverse always exists, so a Dirac bracket of the form (1) can always be computed regardless of the constraints and the original Poisson bracket.

If the Poisson bracket is given by

{F, G} = ∂F ∂z • J(z) ∂G ∂z , ( 5
)
where J is the Poisson matrix, the Dirac bracket {•, •} * has the same expression as Eq. [START_REF] Henneaux | Quantization of Gauge Systems[END_REF] where J is replaced by

J * = J -J Q † D QJ, (6)
where Q has elements Qni = ∂Φ n /∂z i . The linear operator P * = 1 -Q † D QJ is a projector, called Dirac projector [10]. In order to prove this we need condition [START_REF] Hanson | Constrained Hamiltonian Systems[END_REF]. From this projector, we have the following identities :

J * = P † * JP * = JP * ,
which means that the Dirac bracket is the same as the original bracket with the exception that the derivatives ∂F/∂z in the Poisson bracket have to be replaced by the constrained derivatives defined by ∂F ∂z c = P * ∂F ∂z .

The constraints Φ n (z) are Casimir invariants of the Dirac bracket, i.e., {Φ n (z), G} * = 0 for all functionals G, if and only if J * Q † = 0. Therefore a complete reduction corresponds to the case where all constraints are Casimir invariants of the bracket (1). However even if the Dirac bracket can always be computed (at least for finite dimensional systems), J * Q = J Q(1 -DC) is non-zero in general. We will see below that the obstacle to the complete reduction corresponds to the primary constraints which are not Casimir invariants of the original bracket {•, •}.

It has been shown in Ref. [10] that if the reduction is complete, then the Dirac bracket ( 1) is a Poisson bracket, i.e., it satisfies the Jacobi identity everywhere in phase space. Here we prove that, even in the case of an incomplete reduction, the generalized Dirac bracket defined by Eq. ( 1) with the conditions (2)-( 4) is still a Poisson bracket, i.e., in Sec. 2 we prove that the bracket (1) satisfies the Jacobi identity everywhere in phase space. In Sec. 5, we apply the generalized Dirac bracket to several examples, finite and infinite dimensional ones.

Proof of the Jacobi identity

First, we recall that first-class constraints are constraints which Poisson-commute with all the other constraints for the original bracket {•, •}. In particular, constraints which are Casimir invariants of the original bracket are first-class constraints, but obviously first-class constraints are not restricted to Casimir invariants of the original Poisson bracket.

Using a local change of variables, we use the constraints as part of the variables, i.e., we assume Φ k (z) = z k for k = 1, . . . , K. We divide this set of constraints/variables into three groups:

• the first-class constraints which are not Casimir invariants (of the original bracket), for k = 1, . . . , k 1 ,

• the (first-class) constraints which are Casimir invariants, for

k = k 1 + 1, . . . , k 1 + k c , • the second-class constraints, for k = k 1 + k c + 1, . . . , k 1 + k c + k 2 ,
where

K = k 1 + k c + k 2 .
The arrangement of the constraints is done such that there is no linear combination of second-class constraints which is a first-class constraint. The rest of the variables are unchanged. Using this partition of the variables, we rewrite the Poisson matrix J as

J =      0 0 0 -c † 0 0 0 0 0 0 A -B † c 0 B J      . Since z k for k = 1, . . . ,
Q =    1 k 1 0 0 0 0 1 kc 0 0 0 0 1 k 2 0    ,
where 1 k is the k-dimensional identity matrix. The matrix C = QJ Q † is given by

C =    0 0 0 0 0 0 0 0 A    .
The pseudoinverse of C is given by

D =    0 0 0 0 0 0 0 0 A -1    .
The Poisson matrix associated with the Dirac bracket (1) is obtained from Eq. (6) as

J * =      0 0 0 -c † 0 0 0 0 0 0 0 0 c 0 0 J + BA -1 B †      . ( 7
)
From this expression, we readily see that the first-class constraints that are not Casimir invariants of the original bracket, are not Casimir invariants of the Dirac bracket. The only way to resolve these first-class constraints in a Hamiltonian way, i.e., to have these constraints as conserved quantities, is to modify the Hamiltonian (by removing the dependence of the Hamiltonian on the variables corresponding to first-class non-Casimir constraints).

In this section, the goal is to prove the Jacobi identity for J * given by Eq. ( 7) under the hypothesis that J satisfies the Jacobi identity, i.e.,

J il ∂ l J jk + (i,j,k) = 0, (8)
for all (i, j, k). In what follows we denote I 1 the set of indices k for which z k is a non-Casimir first-class constraints, I c the set of indices such that z k is a Casimir invariant of J, I 2 the set of indices such that z k is a second-class constraint, and I n all the other indices. From Eq. ( 8) there are seven non-trivial identities listed below:

• If (i, j) ∈ I 1 and k ∈ I n , c li ∂ l c kj -c lj ∂ l c ki = 0. (9)
• If i ∈ I 1 and j, k ∈ I 2 ,

c li ∂ l A jk = 0. ( 10
)
• If i ∈ I 1 , j ∈ I 2 and k ∈ I n c li ∂ l B kj + A jl ∂ l c ki -B lj ∂ l c ki = 0. ( 11
)
• If i ∈ I 1 and j, k ∈ I n -c li ∂ l Jjk -c kl ∂ l c ji + c jl ∂ l c ki -B kl ∂ l c ji + B jl ∂ l c ki -Jkl ∂ l c ji + Jjl ∂ l c ki = 0. ( 12
)
• If i, j, k ∈ I n c il ∂ l Jjk + B il ∂ l Jjk + Jil ∂ l Jjk + (i,j,k) = 0. ( 13
)
• If i, j ∈ I n and k ∈ I 2 c il ∂ l B jk -c jl ∂ l B ik + B il ∂ l B jk -B jl ∂ l B ik + Jil ∂ l B jk -Jjl ∂ l B ik -B lk ∂ l Jij + A kl ∂ l Jij = 0. ( 14
)
• If i ∈ I n and j, k ∈ I 2 c il ∂ l A jk + A kl ∂ l B ij -A jl ∂ l B ik + B il ∂ l A jk -B lk ∂ l B ij + B lj ∂ l B ik + Jil ∂ l A jk = 0. ( 15
)
• If i, j, k ∈ I 2 A il ∂ l A jk -B li ∂ l A jk + (i,j,k) = 0. ( 16
)
We denote

S ijk = J * il ∂ l J * jk + (i,j,k) . (17)
Using Eqs. ( 9)-( 16), we prove below that S ijk = 0 for all i, j, k. Several cases have to be envisaged, depending on which set i, j and k belong to. There are only two cases where the derivation of the Jacobi identity is non-trivial, one corresponding to i, j, k ∈ I n and the other one corresponding to i ∈ I 1 and j, k ∈ I n .

• If one of the indices i, j, k belong to I 2 or I c , then all the terms J * il ∂ l J * jk vanish individually and the Jacobi identity is trivially satisfied. Consequently i, j, k belong to either I 1 or I n .

• If i, j, k ∈ I 1 , then J * jk = 0 and hence the Jacobi identity is trivially satisfied.

• If i, j ∈ I 1 and k ∈ I n , S ijk = c li ∂ l c kj -c lj ∂ l c ki ,
which vanishes considering Eq. ( 9).

• If i ∈ I 1 and j, k ∈ I n , S ijk = -c li ∂ l ( J + BA -1 B † ) jk -c kl ∂ l c ji + c jl ∂ l c ki -( J + BA -1 B † ) kl ∂ l c ji + ( J + BA -1 B † ) jl ∂ l c ki . (18)
In order to show that S ijk = 0, we need Eqs. ( 10), ( 11) and ( 12). Using Eq. ( 12), S ijk is rewritten as

S ijk = -c li ∂ l (BA -1 B † ) jk -(BA -1 B † ) kl ∂ l c ji + (BA -1 B † ) jl ∂ l c ki + B kl ∂ l c ji -B jl ∂ l c ki .
From AA -1 = 1 k 2 , we have

A -1 βγ ∂ l A αβ + A αβ ∂ l A -1 βγ = 0.
Multiplying the previous identity by c li and using Eq. ( 10), we have

c li A αβ ∂ l A -1 βγ = 0.
Since A is invertible, this identity becomes equivalent to

c li ∂ l A -1 βγ = 0,
which is analogous to Eq. ( 10) for A -1 . In order to prove that S ijk = 0, we compute

c li ∂ l (BA -1 B † ) jk : c li ∂ l (BA -1 B † ) jk = c li ∂ l B jα A -1 αβ B kβ + c li ∂ l B kβ A -1 αβ B jα + c li B jα B kβ ∂ l A -1 αβ .
Using c li ∂ l A -1 αβ = 0 and Eq. ( 11) to rewrite c li ∂ l B jα and c li ∂ l B kβ , we obtain

c li ∂ l (BA -1 B † ) jk = (-A αl ∂ l c ji + B lα ∂ l c ji )A -1 αβ B kβ + (-A βl ∂ l c ki + B lβ ∂ l c ki )A -1 αβ B jα .
Since A is antisymmetric, we rewrite this expression as

c li ∂ l (BA -1 B † ) jk = B kl ∂ l c ji -B jl ∂ l c ki + (BA -1 B † ) lk ∂ l c ji -(BA -1 B † ) lj ∂ l c ki ,
from which we deduce that S ijk = 0.

• The case i, j, k ∈ I n is the most involved one. The expression of S ijk is given by

S ijk = c il ∂ l ( J + BA -1 B † ) jk + ( J + BA -1 B † ) il ∂ l ( J + BA -1 B † ) jk + (i,j,k) .
We start by inserting Eq. ( 13) into S ijk :

S ijk = c il ∂ l (BA -1 B † ) jk + Jil ∂ l (BA -1 B † ) jk +(BA -1 B † ) il ∂ l ( J+BA -1 B † ) jk -B il ∂ l Jjk + (i,j,k) .
Expanding the above expression leads to

S ijk = c il ∂ l B jα A -1 αβ B kβ + c il B jα A -1 αβ ∂ l B kβ + Jil ∂ l B jα A -1 αβ B kβ + Jil B jα A -1 αβ ∂ l B kβ + B iα A -1 αβ B lβ ∂ l B jγ A -1 γδ B kδ + B iα A -1 αβ B lβ B jγ A -1 γδ ∂ l B kδ + c il B jα B kβ ∂ l A -1 αβ + Jil B jα B kβ ∂ l A -1 αβ + B iα A -1 αβ B lβ B jγ B kδ ∂ l A -1 γδ + B iα A -1 αβ B lβ ∂ l Jjk -B il ∂ l Jjk + (i,j,k) .

Using a circular permutation and antisymmetry of

A -1 with which c il B jα A -1 αβ ∂ l B kβ is replaced by -c jl ∂ l B iα A -1 αβ B kβ , the expression for S ijk becomes S ijk = (c il ∂ l B jα -c jl ∂ l B iα + Jil ∂ l B jα -Jjl ∂ l B iα )A -1 αβ B kβ + B iα A -1 αβ B lβ ∂ l B jγ A -1 γδ B kδ + B iα A -1 αβ B lβ B jγ A -1 γδ ∂ l B kδ + c il B jα B kβ ∂ l A -1 αβ + Jil B jα B kβ ∂ l A -1 αβ + B iα A -1 αβ B lβ B jγ B kδ ∂ l A -1 γδ + B iα A -1 αβ B lβ ∂ l Jjk -B il ∂ l Jjk + (i,j,k)
, where we have also replaced B iα B jγ ∂ l B kδ by B kα B iγ ∂ l B jδ by a circular permutation of (i, j, k). Inserting Eq. ( 14) inside S ijk gives

S ijk = (-B il ∂ l B jα + B jl ∂ l B iα + B lα ∂ l Jij -A αl ∂ l Jij )A -1 αβ B kβ + c il B jα B kβ ∂ l A -1 αβ + Jil B jα B kβ ∂ l A -1 αβ + B iα A -1 αβ B lβ B jγ B kδ ∂ l A -1 γδ + B iα A -1 αβ B lβ ∂ l B jγ A -1 γδ B kδ + B iα A -1 αβ B lβ B jγ A -1 γδ ∂ l B kδ + B iα A -1 αβ B lβ ∂ l Jjk -B il ∂ l Jjk + (i,j,k) . (19)
Since

A αl ∂ l Jij A -1 αβ B kβ = -B kl ∂ l Jij , the term B kl ∂ l Jij cancels with -B il ∂ l
Jjk by using a permutation on the indices i, j, k. Similarly, the term

B kβ A -1 αβ B lα ∂ l Jij cancels with B iα A -1
αβ B lβ ∂ l Jjk with a permutation of the indices i, j, k and the antisymmetry of A -1 . Next, we rewrite the terms

U ijk ≡ B iα A -1 αβ B lβ ∂ l B jγ A -1 γδ B kδ + B iα A -1 αβ B lβ B jγ A -1 γδ ∂ l B kδ + (i,j,k) , into U ijk = B iα A -1 αβ A -1 γδ B kδ (B lβ ∂ l B jγ -B lγ ∂ l B jβ ) + (i,j,k)
, using a circular permutation of (i, j, k) and a relabeling of the silent indices (α, β, γ, δ) as well as the antisymmetry of A -1 . Inserting Eq. [START_REF] Marsden | Introduction to Mechanics and Symmetry[END_REF] gives

U ijk = B il A -1 γδ B kδ ∂ l B jγ + B iα A -1 αβ B kl ∂ l B jβ + B iα A -1 αβ A -1 γδ B kδ c jl ∂ l A γβ + B iα A -1 αβ A -1 γδ B kδ B jl ∂ l A γβ + B iα A -1 αβ A -1 γδ B kδ Jjl ∂ l A γβ + (i,j,k) . From the identity A -1 βγ ∂ l A αβ = -A αβ ∂ l A -1
βγ (which comes from differentiating AA -1 = 1 k 2 ), we notice the following cancellations using circular permutations of (i, j, k): The term

B iα A -1 αβ A -1 γδ B kδ c jl ∂ l A γβ cancels with c il B jα B kβ ∂ l A -1 αβ ,

and the same holds for

B iα A -1 αβ A -1 γδ B kδ Jjl ∂ l A γβ which cancels with Jil B jα B kβ ∂ l A -1 αβ .
The first two terms in U ijk cancel with the first two terms of S ijk in Eq. (19). It follows that 20) By differentiating AA -1 = 1 k 2 and using the antisymmetry of A, we have

S ijk = B iα B jγ B kδ B lβ A -1 αβ ∂ l A -1 γδ +B iα B kδ A -1 αβ A -1 γδ B jl ∂ l A γβ + (i,j,k) .(
∂ l A -1 γδ = A -1 γm A -1 δn ∂ l A mn .
Inserting this expression into the first term of S ijk in Eq. (20) gives

S ijk = B iα B jγ B kδ A -1 αβ A -1 δm A -1 nγ B lβ ∂ l A mn + (β,n,m) + B iα B kδ A -1 αβ A -1 γδ B jl ∂ l A γβ + (i,j,k) , since B iα B jγ B kδ B lβ A -1 αβ ∂ l A -1 γδ + (i,j,k) = B iα B jγ B kδ B lβ A -1 αβ ∂ l A -1 γδ + (α,γ,δ) .
Now we use Eq. ( 16) from which we obtain

S ijk = B iα B jγ B kδ A -1 γm A -1 δn ∂ α A mn + (α,γ,δ) + B iα B kδ A -1 αβ A -1 γδ B jl ∂ l A γβ + (i,j,k) , = B il B jγ B kδ A -1 γm A -1 δn ∂ l A mn + B iα B kδ A -1 αβ A -1 γδ B jl ∂ l A γβ + (i,j,k)
. By using a circular permutation of (i, j, k), a relabeling of the silent indices and the antisymmetry of A -1 , we show that the two terms in the last equation for S ijk cancel each other. Therefore S ijk = 0, and the Jacobi identity is satisfied.

Non-unicity of Dirac brackets

Conditions ( 2)-( 4) on the determination of the pseudoinverse D corresponds to a unique D for finite dimensional Hamiltonian systems. A natural question is whether or not one of these conditions can be relaxed with the requirement that the resulting matrix J * given by Eq. ( 6) still satisfies the Jacobi identity. In other words, what is the minimal set of equations which has to be satisfied by D such that the resulting bracket is a Poisson bracket and the second class constraints are Casimir invariants?

We consider a general (antisymmetric) matrix D written in the same coordinates as in the previous section :

D =    D (11) D (12) D (13) -D (12) † D (22) D (23) -D (13) † -D (23) † D (33)    .
From the computation of J * given by Eq. ( 7), we notice that the resulting Dirac bracket is independent of D (12) , D (22) and D (23) . The second line and second column correspond to the Casimir invariants of the original Poisson bracket. In what follows, we choose these matrices to be zero. In addition, a natural choice for D (33) is A -1 in order to resolve the second class constraints. As a result, the matrix J * is given by

J * =      0 0 0 -c † 0 0 0 0 0 0 0 -AD (13) † c † c 0 -cD (13) A J *     
, where J * = J+BA -1 B † +cD (11) c † +cD (13) B † -BD (13) † c † . This matrix J * does not satisfy the Jacobi identity in general. The goal is to find the conditions on D (11) and D (13) such that the resulting Dirac bracket satisfies the Jacobi identity regardless of the specific form of the original Poisson bracket (i.e., for all c, A, B and J satisfying Eqs. ( 9)-( 13)).

For this purpose we first consider S ijk given by Eq. ( 17) for i ∈ I 2 and j, k ∈ I n and look at the terms proportional to D

αβ . After some algebra, the condition S ijk = 0 reduces to D

(13) αβ c kα c jl -c jα c kl + c kα ( J + BA -1 B † ) jl -c jα ( J + BA -1 B † ) kl ∂ l A iβ = 0.
Since the coefficient of D

αβ is non-zero in general, it implies that D (13) = 0. The condition D (13) = 0 implies that the second class constraints are Casimir invariants of the Dirac bracket (1). In other terms it is necessary to have the second class constraints as Casimir invariants in order to satisfy the Jacobi identity for the bracket (1). This result echoes the one in Ref. [10] where it was proven that if all the constraints are Casimir invariants (which is the case for a complete reduction) then the Dirac bracket is a Poisson bracket, i.e., it satisfies the Jacobi identity everywhere in phase space.

Concerning D (11) , there are two conditions to be satisfied :

• The first condition is given by S ijk = 0 when i ∈ I 1 and j, k ∈ I n . This condition reduces to (11) c † ) jk -(cD (11) c † ) kl ∂ l c ji + (cD (11) c † ) jl ∂ l c ki = 0.

S ijk = -c li ∂ l (cD
The terms proportional to D (11) αβ vanish due to Eq. ( 9) and the antisymmetry of D (11) . What is left are terms proportional to the derivatives of D (11) :

S ijk = -c jα c k βc li ∂ l D (11) αβ .
Since S ijk should vanish regardless of c, it implies that D (11) does not depend on the variables z l for l ∈ I n (i.e., the variables which are neither constraints nor Casimir invariants).

• The second condition is obtained for i, j, k ∈ I n as S ijk = c il ∂ l (cD (11) c † ) jk + (cD (11) (11) c † ) jk + (cD (11) c † ) il ∂ l (cD (11) c † ) jk + (i,j,k) . This decomposes in three series of terms, the ones linearly proportional to the coefficients of D (11) , the ones proportional to the derivatives of these coefficients, and the ones which are quadratic in the coefficients of D (11) and their derivatives. The first series of terms is equal to

c † ) il ∂ l ( J + BA -1 B † ) jk + ( J + BA -1 B † ) il ∂ l (cD
S (1) ijk = D (11) αβ c il ∂ l (c jα c kβ ) + c iα c lβ ∂ l ( J + BA -1 B † ) jk + ( J + BA -1 B † ) il ∂ l (c jα c kβ ) + (i,j,k) .
Using Eq. ( 18) together with a circular permutation of the indices (i, j, k) and the antisymmetry of D (11) , all these terms vanish, i.e., S

(1) ijk = 0. For the second series of terms, we have

S (2) ijk = c jα c kβ c il + ( J + BA -1 B † ) il ∂ l D (11)
αβ + (i,j,k) . These terms are in general non-zero, except when D (33) does not depend on the variables z l for l ∈ I 1 and l ∈ I n . As for the third series of terms, it is written as

S (3) ijk = c iα c jγ c kδ c lβ D (11) αβ ∂ l D (11) γδ + (i,j,k) .

The conditions S

(3) ijk = 0 implies that D (11) does not depend on z l for l ∈ I n , which is again the condition found in the previous case i ∈ I 1 and j, k ∈ I n .

In summary, the Jacobi identity is satisfied when D (13) = 0 and when D (11) does not depend on the variables z l for l ∈ I 1 and l ∈ I n , i.e., D (11) could have a dependence on the Casimir invariants of the bracket {•, •} * (whether these invariants originate from the Casimir invariants of the original bracket {•, •} or from the reduction of the secondclass constraints). Of course, depending on the specific choice of Poisson matrix J, less restrictive cases can be considered for D (13) or D (11) . An example is provided in Sec. 5.

Here we provided the conditions on these two matrices such that the Jacobi identity is satisfied regardless of the specific form of the original Poisson matrix.

Consequently, given that there are other possible choices for D (11) , the Dirac bracket (1) as a Poisson bracket is not unique even if all the second class constraints are imposed as Casimir invariants. We notice that with a non-zero D (11) , the resulting matrix D satisfies the conditions ( 2) and ( 4), but the condition (3) is not satisfied. As a consequence, P * = 1 -Q † D QJ is not a projector (with the additional assumption that c † does not belong to the kernel of D (11) ).

Reduced dynamics

We consider the Dirac-Poisson bracket given by the matrix

J * =      0 0 0 -c † 0 0 0 0 0 0 0 0 c 0 0 J *     
, where J * = J + BA -1 B † + cD (11) c † . The equations of motion for z l for l ∈ I 1 and l ∈ I n are given by ż1 = -c † ∂H ∂z n ,

żn = J * ∂H ∂z n + c ∂H ∂z 1 ,
where H is the Hamiltonian which a priori depends on z 1 and z n (since z 2 and z c are Casimir invariants, they can be forgotten in the analysis). In order to have a complete reduction, the first-class constraints have to be conserved quantities (even if they are not Casimir invariants) which is obtained under the condition

c † ∂H ∂z n = 0. ( 21
)
Under this condition, we notice that the resulting system of dynamical equations is unique, i.e., it does not depend on D (11) since the additional term in the equations cD (11) c † ∂H/∂z n vanishes. The condition ( 21) is rewritten as {z 1 , H} * = 0 since {z 1 , z 2 } * = {z 1 , z c } * = 0 (where z 2 denotes z l with l ∈ I 2 and z c denotes z l with l ∈ I c ), which means that the Hamiltonian after the reduction has to commute with the first-class constraints. In other terms, the Hamiltonian has to be changed so that the first-class constraints are conserved quantities of the Dirac bracket.

Examples

There are two trivial examples and they correspond to the two limits of the Dirac reduction, one for which all the constraints are first-class and another for which all the constraints are second-class. For the case where all the constraints are first-class, the matrix C is zero and its pseudoinverse D is also zero. Hence the Dirac bracket is the same as the original bracket, which means that the reduction has failed. The second limit example corresponds to the case where the matrix C is invertible. The reduction is complete and the Dirac bracket is the usual one.

Between these two limit cases, there are examples where the reduction is incomplete. Several examples are given below.

Finite dimensional examples

If we consider only one constraint, it is necessary a first-class constraint since it Poisson commutes with itself. Therefore the Dirac bracket computed with the pseudoinverse is identical to the original Poisson bracket.

If we consider two constraints, there are two cases: one in which they do not commute and hence both constraints are second class, and another one where they commute and hence they are both first-class. In the first case, the Dirac bracket is identical to the usual one since the matrix C is invertible. In the second case, the Dirac bracket is identical to the original one since C = 0.

More interesting cases occur when there are three constraints, one of which is first-class. The simplest example is afforded by the Poisson bracket

J =      0 0 0 1 0 0 -1 1 0 1 0 -1 -1 -1 1 0     
, with the three constraints Φ 1 (z) = z 1 , Φ 2 (z) = z 2 and Φ 3 (z) = z 3 . We notice that Φ 1 is a first-class constraint whereas Φ 2 and Φ 3 are second-class. The matrix C is given by

C =    0 0 0 0 0 -1 0 1 0    ,
and its pseudoinverse by

D =    0 0 0 0 0 1 0 -1 0    .
The Dirac bracket computed from the pseudoinverse is given by

J * =      0 0 0 1 0 0 0 0 0 0 0 0 -1 0 0 0     
, where we notice that Φ 2 and Φ 3 are Casimir invariants of J * , but Φ 1 is not.

The second example in finite dimensions is given by

J =          0 -z 3 z 2 0 -w 3 w 2 z 3 0 -z 1 w 3 0 -w 1 -z 2 z 1 0 -w 2 w 1 0 0 -w 3 w 2 0 0 0 w 3 0 -w 1 0 0 0 -w 2 w 1 0 0 0 0         
, which corresponds to the Poisson matrix for the rigid body [START_REF] Marsden | Introduction to Mechanics and Symmetry[END_REF], and three constraints given by Φ 1 (z, w) = z 3 , Φ 2 (z, w) = w 2 and Φ 3 (z, w) = w 3 . We notice that Φ 3 is a first-class constraint. The matrix C is given by

C =    0 w 1 0 -w 1 0 0 0 0 0    ,
and its pseudoinverse by

D =    0 -1/w 1 0 1/w 1 0 0 0 0 0    .
The Dirac bracket computed from the pseudoinverse is given by

J * =          0 z 1 w 3 /w 1 -z 3 0 -w 2 w 3 /w 1 0 w 2 -z 1 w 3 /w 1 + z 3 0 0 w 3 0 -w 1 0 0 0 0 0 0 w 2 w 3 /w 1 -w 3 0 0 0 0 0 0 0 0 0 0 -w 2 w 1 0 0 0 0          , (22)
with the same remark as above, Φ 3 is not a Casimir invariant of the Dirac bracket.

In order to illustrate the remark above on the unicity of the Dirac bracket, we consider matrices D which only satisfy condition (2). For instance a possible matrix is

D =    0 -1/w 1 0 1/w 1 0 -w 2 0 w 2 0    ,
which satisfies C = CDC and D = DCD but not CD = DC. This choice leads to the Dirac bracket with the following Poisson matrix :

J * =          0 -w 1 w 2 w 3 + z 1 w 3 /w 1 -z 3 w 1 w 2 2 -w 2 w 3 /w 1 0 w 2 w 1 w 2 w 3 -z 1 w 3 /w 1 + z 3 0 -w 2 1 w 2 w 3 0 -w 1 -w 1 w 2 2 w 2 1 w 2 0 0 0 0 w 2 w 3 /w 1 -w 3 0 0 0 0 0 0 0 0 0 0 -w 2 w 1 0 0 0 0         
, which is different from J * given by Eq. (22).

Infinite dimensional example: two-dimensional Euler's equation

For an infinite-dimensional Hamiltonian system with dynamical field variables χ(x), the generalized Dirac bracket is written as

{F, G} * = {F, G} -d N x d N x ′ {F, Φ n (x)}D nm (x, x ′ ){Φ m (x ′ ), G},
where D satisfies

C nm (x, x ′ ) = d N y d N y ′ C nl (x, y)D lp (y, y ′ )C pm (y ′ , x ′ ), D nm (x, x ′ ) = d N y d N y ′ D nl (x, y)C lp (y, y ′ )D pm (y ′ , x ′ ), d N yC nl (x, y)D lm (y, x ′ ) = d N yD nl (x, y)C lm (y, x ′ ),
which generalize the conditions (2)-( 4). The existence of such a pseudoinverse D is more complicated than in the finite dimensional case, and is beyond the scope of the present work. On a practical basis it has to be assessed case by case. In the infinite dimensional case, it is very difficult to separate first class and second class constraints as we shall see it below using an example. However it should be noted that the principles of the definition of the Dirac bracket does not rely on the separation between first class and second class constraints.

Here we consider the two dimensional Euler equations for the density ρ(x), the fluid velocity v(x) and the entropy s(x) where x = (x, y) ∈ T 2 (the two-dimensional torus) given by ρ

= -∇ • (ρv), v = -v • ∇v -ρ -1 ∇P, ṡ = -v • ∇s.
The Hamiltonian of the system is a function of the dynamical field variables ρ, v and s, given by

H[ρ, v, s] = d 2 x ρ v 2 2 + ρU(ρ, s) , ( 23
)
where U is the internal energy such that the pressure is given by P = ρ 2 ∂U/∂ρ, and the Poisson bracket is

{F, G} = -d 2 x F ρ ∇ • G v + F v • ∇G ρ -ρ -1 (∇ × v) • F v × G v +ρ -1 ∇s • (F s G v -F v G s ) , (24)
where F ρ denotes the functional derivative of the observable F with respect to the field variable ρ (and equivalently for F v and F s ). The associated Poisson operator is

J =    0 -∇• 0 -∇ -ρ -1 (∇ × v)× ρ -1 ∇s 0 -ρ -1 ∇s 0    ,
from which the Poisson bracket is defined as

{F, G} = d 2 xF χ • JG χ ,
where

F χ = (F ρ , F v , F s ).
In order to find reduced systems defined from some constraints imposed on the above Hamiltonian system, a first obvious choice would be to impose ρ(x) = ρ 0 (x) with a prescribed density ρ 0 . Since {ρ(x), ρ(x ′ )} = 0, these constraints are all first class constraints and not of them are Casimir invariants. The matrix C is zero and hence its pseudo-inverse is also zero. The associated Dirac bracket is the same as the original bracket. The reduction has failed.

A second choice is to impose a local constraint as divergence-free velocity field which consists of an infinite number of constraints, Φ(x) = ∇ • v(x) = 0. The idea behind this calculation is to have a divergence free velocity field (incompressibility) but keeping the density as a dynamical field variable. It should be noted that none of these constraints are Casimir invariants of the bracket (24).

The operator C = QJ Q † where Q = (0, ∇•, 0), is given by

C = ∇ • (ρ -1 (∇ × v) × ∇).
We rewrite this operator using q(x, y) = ρ -1 (∇ × v) • ẑ :

C = -[q, •], where the bracket [•, •] is given by [f, g] = ∂f ∂x ∂g ∂y - ∂f ∂y ∂g ∂x .
The operator C is not invertible since any function of q belongs to its kernel. Therefore the constraints form a complicated admixture of first-class and second-class constraints.

In order to see whether or not C has a pseudoinverse, we consider the characteristics associated with the linear operator C which are fictitious trajectories with q(x, y) as Hamiltonian. In this fictitious dynamics, the variables are x and y. Locally, we change variables to action-angle variables (ϕ, A) so that q(x, y) = q 0 (A). Since this change of fictitious variables is canonical, the expression of the linear operator C becomes

C = q ′ 0 (A) ∂ ∂ϕ .
In Fourier series, it is straightforward to see that C has a pseudoinverse D, and the action of these two operators are given by Cf = q ′ 0 (A)

k ikf k e ikϕ , Df = 1 q ′ 0 (A) k =0 f k ik e ikϕ .
The Poisson operator J * associated with the Dirac bracket defined from the pseudoinverse D is given by

J * =    ∆D∆ -∇ • P -∆D∇ • (ρ -1 ∇s•) -P † ∇ -ρ -1 (∇ × v) × P P † (ρ -1 ∇s•) ρ -1 ∇s • ∇D∆ -ρ -1 ∇s • P -ρ -1 ∇s • ∇D∇ • (ρ -1 ∇s•)    ,
where • indicates where the operator acts, and P = 1 -∇D∇ • (ρ -1 (∇ × v) × •) is a projector with the following properties :

∇ × P = P, P∇ = ∇(1 -DC),
where 1 -DC is the projector onto the kernel of C. In fictitious action-angle variables, this projector is given by 1 -DC = dϕ.

From the Poisson matrix J * the Dirac bracket becomes :

{F, G} * = d 2 x (∇F ρ -ρ -1 ∇sF s ) • Ḡv -Fv • (∇G ρ -ρ -1 ∇sG s ) +ρ -1 (∇ × v) • Fv × Ḡv , (25)
where the constrained derivative Fv is defined from the Dirac projector P * = 1-Q † D QJ as Fv = PF v -∇D∇ • (∇F ρ -ρ -1 ∇sF s ).

For Hamiltonian (23), the constrained derivative Hv leads to the definition of a constrained velocity field v such that Hv = ρv with

v = v -ρ -1 ∇D∇ • [(v • ∇)v + ρ -1 ∇P ].
It should be noted that a priori ∇ • v is non-zero. Functionals C satisfying the following conditions :

∇C ρ -ρ -1 ∇sC s = 0, Not all the constraints have been fulfilled since

PC v =
∇ • v = -(1 -CD)∇ • (v • ∇)v + ρ -1 ∇P ,
which is in general non-zero since C is not invertible. In order to have a fully incompressible (in the sense of ∇ • v = 0) model, the pressure needs to be adjusted such that ∇ • v = 0 when ∇ • v = 0 as it is the case in the standard incompressible model. In summary, the Dirac procedure as outlined above provides a Dirac bracket which is a Poisson bracket, given by Eq. ( 25), but does not have all the constraints as Casimir invariants, contrary to the usual procedure where the matrix of the Poisson brackets between the various constraints is invertible. The origin of this incomplete reduction is the presence of first class constraints which are not Casimir invariants of the original bracket [21].

  0, are Casimir invariants of the bracket (25). The first condition leads to C = d 2 xρψ(s) which are also Casimir invariants of the original bracket (24). The reduction procedure has preserved the Casimir invariants as expected. The second equation states that the second-class constraints are Casimir invariants. However due to the intricacy of first and second class constraints, it is cumbersome to explicit those invariants. The resulting equations of motion are the following ones : ρ = -∇ • (ρv), v = -P † (v • ∇)v + 1 ρ ∇P , ṡ = -v • ∇s.

  k 1 are first-class constraints, the first column has three zeros, but the last element is non-zero (otherwise they would be Casimir invariants). The second column is zero since z k are Casimir invariants for k = k 1 + 1, . . . , k 1 + k c . The matrix A

	is invertible, since otherwise there would exist a linear combination of the second-class
	constraints which would be a first-class constraints, which is excluded by construction.
	The operator Q is given by

Acknowledgments

This work was supported by the Agence Nationale de la Recherche (ANR GYPSI) and by the European Community under the contract of Association between EURATOM, CEA, and the French Research Federation for fusion study. The views and opinions expressed herein do not necessarily reflect those of the European Commission. CC acknowledges fruitful discussions with P.J. Morrison and with the Équipe de Dynamique Nonlinéaire of the Centre de Physique Théorique of Marseille.


Édouard Canot
email: edouard.canot@irisa.fr

Renaud Delannay

Salwa Mansour

Mohamad Muhieddine

Ramiro March

Effective thermal conductivity of a wet porous mediumpresence of hysteresis when modeling the spatial water distribution for the pendular regime

Keywords: Wet porous media, Evaporation, Effective thermal conductivity, Hysteresis, Contact

This paper deals with the heat transfer between two spherical grains separated by a small gap; dry air is located outside the grains and a liquid water meniscus is supposed to be present between them. This problem can be seen as a micro-scale cell of an assembly of solid grains, for which we are looking for the effective thermal conductivity. For a fixed contact angle and according to the volume of the liquid meniscus, two different shapes are possible for the meniscus, giving a "contacting" state (when the liquid makes a true bridge between the two spheres) and a "non-contacting" one (when the liquid is split in two different drops, separated by a thin air layer); the transition between these two states occurs at different times when increasing or decreasing the liquid volume, thus leading to a hysteresis behavior when computing the thermal flux across the domain.

Introduction

The effective thermal conductivity is an important parameter for the models which describe the evaporation of water in porous media. Many publications deal with heat transfer through granular media, via two-component models or experiments ([1], [START_REF] Siu | Effective conductivity computation of a packed bed using constriction resistance and contact angle effects[END_REF], [8]). On the other hand, it is known that the presence of water affects a lot the thermal properties because the thermal conductivity of water is typically 25 times that of dry air. This has been mentioned by Chen [START_REF] Chen | Thermal conductivity of sands[END_REF] who explicitly emphasized the strong variation of the conductivity for very small quantity of water. Other papers ( [START_REF] Usowicz | Thermal conductivity modelling of terrestrial soil media -A comparative study[END_REF], [START_REF] Yun | Fundamental study of thermal conduction in dry soils[END_REF]) introduce the presence of liquid water as a component but without specifying its spatial localisation. The interesting paper of Mitarai & Nakanishi [START_REF] Mitarai | Simple model for wet granular materials with liquid clusters[END_REF] deals with the pendular regime (fig. 1) in which the liquid menisci act as a bridge between the grains but their study is restricted to the mechanical properties, not the thermal ones.

When a granular media is constituted by grains which are more or less spherical in shape, the liquid water is located around each contact zone between the grains due to the fact that liquid tends to minimize its surface energy. As stated above, we can expect that even a very small quantity of water increases the effective thermal conductivity by an important factor. No theoretical model for this three-component medium exists and it is difficult to be discovered by experimental results. However, the values of this effective property can be found by numerical simulations, using the Fourier phenomenological law.

In this paper, we focus on the heat transfer between two spherical solid grains of same radius with the presence of liquid water attached to the solid grains.

Surface tension is taken into account by computing the exact shape of the liquid meniscus (which must have a constant total curvature) and the role of the contact angle between the liquid-gas interface and the solid surface is also included in our model. Applications concern mainly the drying of materials and agriculture.

Micro-scale model

We want to solve the steady-state heat transfer in an elementary cell containing only two spherical solid grains, a few quantity of liquid water, and dry air. Actually, the liquid water should be in thermodynamic equilibrium with its vapor, so the surrounding gas should be a binary mixture of dry air and water vapor. In the present study, the water vapor is neglected and we assume that there is no adsorbed water on the grain surface.

Top and bottom sides are kept at constant temperature whereas the vertical ones are supposed to be isolated (this is related to the model used by [START_REF] Yun | Fundamental study of thermal conduction in dry soils[END_REF] in their study of thermal conduction in dry soils). The computational domain can be found using the symmetry planes of the problem knowing the main direction of the heat flux -for example, the bottom symmetry plane comes from the fact that the gravity is neglected. As stated in the introduction, the liquid meniscus is centered along the axis between the two spheres leading to an axisymmetric problem, so we can reduce the geometrical independent variables to the (r, z) cylindrical coordinates. Each component (solid, liquid water and dry air) has a constant thermal conductivity and the thermal contact between the components is supposed to be perfect.

The shape of the liquid meniscus must have both a constant total curvature (it is a minimum surface) and a prescribed contact angle with the solid. The percentage of liquid water, i.e. the ratio of its volume to that of the available space between the spheres is called humidity.

The input geometrical parameters are: R : the radius of the two solid spheres ε : the gap between the spheres V : the volume of the liquid water θ : the contact angle between the liquid interface and the solid. In this study, the radius R is equal to 1 (that is, all other geometrical lengths are expressed in terms of R), whereas ε and θ are kept constant (resp. 0.1 and 30 degrees), so the only input variable is V. The output variable of the problem is the (vertical) thermal flux Φ across the domain, which is related to the effective thermal conductivity of our composite medium. Note that in a complete model, the liquid volume V should be related to the gas pressure but this is out of the scope of this study.

Equilibrium shape of the liquid meniscus

The geometrical shape of the liquid meniscus is obtained by integrating a differential algebraic system of equations, stating that:

• the total curvature C is constant (but unknown); • the contact angle at the solid surface has the prescribed value of θ; • the volume of the liquid meniscus is V. Actually, the solution is obtained numerically by a shooting method over the two parameters C and the starting position of the point (r0 in the "contacting" state, or z0 in the "non-contacting" state); these two parameters are varied until the prescribed values (θ and V) are found. The family of the curves obtained, in both states, are drawn in the next figures.

Numerical computation of the heat flux

Solving

our problem numerically appears unavoidable: some authors (resp. [START_REF] Batchelor | Thermal or electrical conduction through a granular material[END_REF], [7]) have used analytical solutions but they are respectively restricted to an asymptotic behavior (so, adapted to a local geometrical zone) or to too crude approximations.

The steady-state heat transfer is solved using a Finite Volume scheme applied on a regular structured mesh of rectangular cells. The matrix of the linear system is stored in a sparse way and the UMFPack linear sparse solver (from SuiteSparse [ref]) is used. Some tries have shown that a 500 by 500 mesh is required to obtain a good accuracy (fig. 5). Each numerical computation takes about few seconds on a laptop (Intel Core i7 @ 2.7 GHz). During the evaporation/condensation of water in a wet granular medium, the liquid volume in each elementary cell changes with time. It is interesting to know the variation of the effective thermal conductivity with respect to the liquid volume. By using the Fourier law on our computational domain, the effective thermal conductivity is proportional to the heat flux. This heat flux depends of course on the water liquid volume because, as stated in the introduction, the thermal conductivity of the water is much more important than that of the dry air. Further, we expect a jump in the heat flux curve because there are two possible geometrical configurations for the liquid meniscus (the "contacting" state and the "noncontacting" one, as described in figures 4a and 4b).

Below, the heat flux is presented in figure 6-a, for humidity (or liquid volume) ranged from 0 to 20%. As expected, a hysteresis behavior is highlighted by our numerical computations and is more visible in figure 6-b which presents an enlarged view of the previous one.

Conclusion and perspectives.

A strong hysteresis behavior for the effective thermal conductivity has been shown when changing the humidity of a granular medium. It is due to a switch between two different geometrical configurations of the liquid meniscus attached to the two solid spheres. The presented results are preliminary and work is in progress to obtain more advanced ones especially when other parameters are varying (contact angle θ, gap between the grains ε).

Fig. 1 .

 1 Fig. 1. Real granular medium containing few quantity of water (pendular regime)..

Fig. 2 .

 2 Fig. 2. Simplification of the grains' assembly. The dashed rectangle is the computational domain used in this paper.

Fig. 3 .

 3 Fig. 3. Sketch of the cylindrical computational domain (taking into account all the symmetries). Boundary conditions are of homogeneous Dirichlet type on top and bottom sides, and of Neumann type on the external vertical side.

Fig. 4 -Fig. 4 -

 44 Fig. 4-a. This represents the family of the liquid meniscus when the volume varies, in the "contacting" state. Each curve is the boundary of the half of a bridge linking the two grains.

Fig. 5 .

 5 Fig. 5. The contour curves of the temperature obtained after a numerical computation on a 500x500 mesh.

Fig. 6 -Fig. 6 -

 66 Fig. 6-a. Thermal heat flux w.r.t. the liquid volume, for the "contacting" state (blue) and the "non-contacting" state (red). Note the great increase in the heat flux due to the presence of water: about a factor 5 when humidity is only 20%.