El-Houssaine Quenjel
email: el-houssaine.quenjel@centralesupelec.fr

Patrick Perré
email: patrick.perre@centralesupelec.fr.

Ian Turner
email: i.turner@qut.edu.au.

A 3D face interpolated discretisation method for simulating anisotropic diusive processes on meshes coming from wood morphology

Keywords: Diusive equations, 3D discrete gradient, anisotropy, second order accuracy, wood morphology, thermal conductivity

In this paper, a new 3D numerical discretisation method for solving the anisotropic, steady-state diusion problem is developed and analysed. The scheme is constructed on hexahedral meshes using the geometrical properties of the cells. Indeed, each cell provides a local basis formed using the centres of its lateral faces. The discrete cell gradient approximation is then obtained by using three discrete directional derivatives resolved in terms of this basis, and by invoking the consistent relationships between opposite faces of the same cell. The face degrees of freedom are interpolated to reduce the complexity of the numerical scheme, resulting in the main unknowns being entirely nodal based. The scheme is unconditionally coercive and admits a unique solution. Various numerical experiments are performed to highlight the accuracy and the robustness of the method with respect to the mesh and anisotropy. An important outcome is that second order convergence is observed for all problems considered, even for highly deformed meshes. After this validation process, the method is applied to the prediction of the eective thermal conductivity of wood from its real 3D morphology. The property 1 is estimated in the radial tangential and longitudinal directions. The solver is robust, ecient and yields to similar results compared to recent contributions.

Introduction

Partial dierential equations (PDEs) are present in an uncountable number of applications in science and engineering. As a consequence, there is an increasing demand for solving PDEs numerically. In particular, linear and nonlinear diusion terms are crucial components of many physical systems and are used to model, for example, complex ows in porous media, and heat and mass transfer in highly heterogeneous materials such as wood and plant stems [START_REF] Perré | A 3D version of Transpore: a comprehensive heat and mass transfer computational model for simulating the drying of porous media[END_REF]. The discretisation of various types of diusion problems (linear, nonlinear, transient or steady-state) on generic 3D meshes is a challenging topic because the reconstruction of the gradient from a nite set of points must be subject to consistency and stability rules. In this work we propose an original idea to devise a whole gradient across a hexahedral mesh that takes into account the shape and geometric properties of each cell. For this purpose, we will focus on a very standard linear, anisotropic Poisson equation arising in many applications.

Let Ω be an open bounded connected domain of R 3 . Let us consider the 3D elliptic model -div(Λ∇u) = f, in Ω, u = g, on ∂Ω,

where Λ is a symmetric positive-denite tensor dened from Ω to L ∞ (Ω) 3×3 that fulls the ellipticity condition

λ 0 ζ 2 ≤ Λζ • ζ ≤ λ 1 ζ 2 , ∀ζ ∈ R 3 ,
and λ 0 and λ 1 are some positive constants. The source term is a function of L 2 (Ω) and g ∈ H 1/2 (∂Ω). Under these assumptions, it is well-known that the model problem admits a unique weak solution u ∈ H 1 (Ω) whose trace is g and solves the variational formulation

Ω Λ∇u • ∇ϕ dx = Ω f ϕ dx, ∀ϕ ∈ H 1 0 (Ω). (2)
During the last two decades, it has been established that a good and proper approximation of the gradient can lead to an ecient and accurate resolution of the problem [START_REF] Aavatsmark | Discretization on unstructured grids for inhomogeneous, anisotropic media. Part II: Discussion and numerical results[END_REF]. Depending on the data, namely Λ, f , g, Ω and the mesh, many numerical methods were developed to conceive approximate directional derivatives as in the nite dierence approach. Their accuracy and robustness are observed through the FVCA benchmarks [START_REF] Eymard | 3D benchmark on discretization schemes for anisotropic diusion problems on general grids[END_REF][START_REF] Herbin | Benchmark on discretization schemes for anisotropic diusion problems on general grids[END_REF]. Many of them share several common points in terms of analysis, which renders them eligible to enter the theoretical framework called the gradient discretisation method [START_REF] Droniou | The gradient discretisation method[END_REF]. This paradigm includes shape function based classical methods such as the P 1 -nite elements [START_REF] Ciarlet | The nite element method for elliptic problems[END_REF]. It also contains modern methods and we mention among them the Hybrid Mimetic Mixed approaches [START_REF] Da Veiga | The mimetic nite dierence method for elliptic problems[END_REF][START_REF] Droniou | A unied approach to mimetic nite dierence, hybrid nite volume and mixed nite volume methods[END_REF][START_REF] Eymard | Discretization of heterogeneous and anisotropic diusion problems on general nonconforming meshes SUSHI: a scheme using stabilization and hybrid interfaces[END_REF], Multi-Point Flux Approximations [START_REF] Droniou | A unied approach to mimetic nite dierence, hybrid nite volume and mixed nite volume methods[END_REF], Diamond schemes [START_REF] Coudière | Convergence rate of a nite volume scheme for a two dimensional convection-diusion problem[END_REF][START_REF] Quenjel | Node-Diamond approximation of heterogeneous and anisotropic diusion systems on arbitrary two-dimensional grids[END_REF], and Discrete Duality Finite Volume methods [START_REF] Andreianov | On 3D DDFV discretization of gradient and divergence operators. I. Meshing, operators and discrete duality[END_REF][START_REF] Andreianov | Discrete duality nite volume schemes for LerayLionstype elliptic problems on general 2D meshes[END_REF][START_REF] Coudière | A 2D/3D discrete duality nite volume scheme. Application to ECG simulation[END_REF][START_REF] Domelevo | A nite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF][START_REF] Hermeline | A nite volume method for approximating 3D diusion operators on general meshes[END_REF]. In addition to the primary degrees of freedom (DOFs), many of the aforementioned methods require auxiliary unknowns in the reconstruction of the gradient as the mesh becomes distorted.

In this work, we rst propose and validate a new 3D numerical method for the resolution of [START_REF] Aavatsmark | Discretization on unstructured grids for inhomogeneous, anisotropic media. Part II: Discussion and numerical results[END_REF]. The motivation for this scheme results from the need to use distorted hexahedral meshes with larger sets of DOFs, without introducing auxiliary unknowns, as well as tackling the challenges of higher anisotropy ratios while still ensuring a straightforward implementation strategy. For instance, such constraints are naturally imposed during the prediction of macroscopic properties from real morphologies using computational homogenisation [START_REF] Carr | A dual-scale modeling approach for drying hygroscopic porous media[END_REF][START_REF] Louërat | Heat and moisture diusion in spruce and wood panels computed from 3-d morphologies using the lattice boltzmann method[END_REF][START_REF] Lux | Macroscopic properties of real brous materials: Volume averaging method and 3d image analysis[END_REF][START_REF] Perré | Determination of the material property variations across the growth ring of softwood for use in a heterogeneous drying model. part 2. use of homogenisation to predict bound liquid diusivity and thermal conductivity[END_REF][START_REF] Quenjel | Computation of the eective thermal conductivity from 3d real morphologies of wood[END_REF], or to track liquid imbibition in complex media such as natural wood [START_REF] Perré | A macroscopic washburn approach of liquid imbibition in wood derived from x-ray tomography observations[END_REF].

The contribution of the paper is twofold. The rst asset lies in the scheme construction. More importantly, the main idea of our work is to dene the numerical scheme using only the geometrical information of the mesh cells. Indeed, following the spirit of the diamond and the discrete duality schemes, one sets three discrete directional derivatives resolved in terms of a local basis, provided by the opposite faces for each hexahedron. This strategy allows the reconstruction of the gradient on the entire cell and therefore the mesh. Then, the numerical scheme stems from the discrete weak formulation. The proofs of the stability results, such as the discrete Poincaré's inequality or the derivation of the energy estimate, follows the guidelines of the literature. As a consequence, the stiness matrix is symmetric positive-denite and the scheme stencil is small and compact in the sense that any node is connected to at most 26 neighbouring nodes in the case of Cartesian meshes. The resulting linear scheme is solved using the conjugate gradient method. This method exhibits some ineciency when the mesh is severely distorted and the diusion tensor is highly anisotropic. We therefore implemented the SSOR (symmetric successive over-relaxation) preconditioner [START_REF] Saad | Iterative methods for sparse linear systems[END_REF] in order to improve the convergence behaviour of the linear solver.

The second contribution consists in the application of the 3D face interpolated scheme to the computation of the eective macroscopic thermal conductivity of wood. Thanks to image processing, the employed meshes are generated from the scanned real morphology of the sample. Then, the property is estimated following the orthotropic directions of wood. Compared to some recent works, the proposed numerical scheme is capable to recover similar ranges of the property in an excellent way and despite of the mesh deformation.

The rest of the paper is outlined as follows. In the next section, we survey the domain meshing and elaborate the approximation of the gradient operator. Then, we derive the 3D face interpolated discretisation method and study its stability together with solvability. In Section 3, we implement the numerical scheme and report the obtained results to conrm the eciency and the robustness of our methodology. In Section 4, we put into practice the 3D face interpolated scheme for the prediction of thermal conductivity of wood. The paper is then concluded in Section 5.

Numerical scheme

In this section we present the discretisation of the domain Ω. We introduce a consistent and accurate approximation of the gradient operator, and then derive the proposed 3D face interpolated scheme approximating the diusion model (2).

Mesh, notations and discrete gradient

The mesh M is given by a collection of conforming hexahedral cells {K} that cover the whole domain Ω. We assume that the volume faces are planar. In other words, the scheme does not handle generalised hexahedra with non-planar faces. The centre of the cell is determined by averaging the coordinates of its surrounding vertices. The volume of K is denoted by |K|. Each cell possess six faces that are expressed as F i , i ∈ {1, • • • , 6}. We denote by x Fi the barycentre of the face F i . Let T denote the set of vertices of the mesh, and let T D denote the set of boundary vertices. Let M s represent the cells sharing the same vertex s. For all s ∈ T , we assign a local sub-set ω K,s in the cell K, determined by the centre of the latter, the centres of the edges, and the faces sharing the same vertex, see Figure 1. The volume of ω K,s is computed using the identity

|ω K,s | = |Det -→ x s x e 1 K,s , -→ x s x e 2 K,s , -→ x s x e 3 K,s |,
where e 1 1,s , e 2 2,s , e 3 K,s are the edges of K whose extremity is s. Therefore, the total sub-volume associated with s and its measure are dened by

ω s = K∈Ms ω K,s , |ω s | = K∈Ms |ω K,s |.
The sum of these sub-volumes is equal to the volume of Ω. For i ∈ {1, 2, 3}, let us x the vector ------→ τ F2i-1F2i = -------→ x F2i-1 x F2i . The set T Fi contains the four vertices of the face F i .

Let h K denote the diameter of the cell K ∈ M. The mesh size is given by h T = max{h K , K ∈ M}. We denote by r K the radius of the largest sphere included in K. Consider

θ T = max K∈M h K r K , n T = max s∈T #M s , ξ T = max(θ T , n T ).
As T is rened, the mesh regularity assumption states that ξ T must be bounded with constant ξ 0 > 0 i.e., ξ T ≤ ξ 0 . In hexahedral meshes it is not evident how to construct a nodal (i.e., primary unknowns are nodal) accurate approximation of only one gradient without auxiliary unknowns, or without using the traditional shape functions of the well-known nite element framework.

The key point of our approach lies in the denition of a consistent whole discrete gradient per cell. The six faces of K yield three noncollinear directions forming a local basis in 3D, see Figure 2. First, in each cell K, we begin to set the approximation of the gradient to the formula

∇ K u T = 1 |K| u F2 -u F1 --→ N 12 + u F4 -u F3 --→ N 34 + u F6 -u F5 --→ N 56 , (3)
where An important observation is that the total number of faces of the mesh is larger than the number of vertices of the mesh. Consequently, the degrees of freedom are placed on the mesh vertices in order to reduce the complexity of the numerical scheme. For this purpose, because the faces are planar, the face unknowns can be eliminated by considering the second order interpolation

--→ N 12 = ---→ τ F3F4 ∧ ---→ τ F5F6 , --→ N 34 = ---→ τ F5F6 ∧ ---→ τ F1F2 , --→ N 56 = ---→ τ F1F2 ∧ ---→ τ F3F4 .
u Fi = 1 4 s∈T F i u s .
To simplify the exposition, we assume that the boundary datum g can be extended to a function that is still denoted g ∈ H 1 (Ω). This gives sense to the following integral approximation for the average of the function over the volume ω s

g s = 1 |ω s | ωs g(a) da.
Dene R #T to be the nite dimensional space containing vectors of the form u T = (u s ) t s∈T . We denote by R #T g (resp. R #T 0 ) the subspace for which u T is such that u s = g s (resp. u s = 0), for all s ∈ T D . For each u T ∈ R #T , we associate a unique function I h u T of L 2 (Ω) that is expressed as

I h u T = s∈T u s 1 ωs ,
where 1 ωs is the indicator function of the subset ω s , i.e. 1 ωs (x) = 1 for x ∈ ω s and is equal to 0 otherwise. Let ∇ h be the piecewise operator mapping 3 where

u T ∈ R #T into ∇ h u T ∈ L 2 (Ω)
∇ h u T = K∈M 1 K ∇ K u.
The tensor Λ is approximated by Λ h M where the latter is piecewise constant on the cells of the mesh and written

Λ h M = K∈M 1 K Λ K , Λ K = 1 |K| K Λ(x) dx.
Dene the norm • 0,2 and the semi-norm

• 1,2 on R #T such that v T 0,2 = s∈T |ω s |v 2 s 1/2 , v T 1,2 = K∈M |K||∇ K v T | 2 1/2
.

3D face interpolated scheme

Having the discrete gradient at hand, the 3D face interpolated scheme is naturally derived from the discrete variational formulation: i.e., nd

u T ∈ R #T g satisfying the relationship Ω Λ h M ∇ h u T • ∇ h ϕ T dx = Ω f I h ϕ T dx, ∀ϕ T ∈ R #T 0 . (4)
The pointwise version of the scheme can be equivalently given by

   K∈Ms |K|Λ K ∇ K u T • ∇ K e s = |ω s |f s , ∀s ∈ T \ T D , u s = g s ∀s ∈ T D , (5)
where {e s } s∈T denotes the canonical basis of R #T . For the sake of simplicity, assume that g = 0. The scheme (5) can then be recast in the matrix form

A T u T = F T , (6)
where

(F T ) s = |ω s |f s , if s ∈ T \ T D , 0 if s ∈ T D .
The matrix A T contains the stiness coecients. In view of the implementation perspective, the formulation ( 5) is not adequate. In practice, a possible way to ll A T is described in the sequel. For this purpose, let us consider the numbering strategy illustrated in Figure 3. Each vertex of the cell K belongs to three faces of the same cell. Using this orientation convention, the expression of the gradient already described in ( 3) is now rewritten as follows

∇ K u T = 8 j=1 u j a K j ,
where the vectors (a K j ) 1≤j≤8 are given by

a K 1 = 1 4|K| - --→ N 12 - --→ N 34 - --→ N 56 , a K 2 = 1 4|K| + --→ N 12 - --→ N 34 - --→ N 56 , a K 3 = 1 4|K| + --→ N 12 + --→ N 34 - --→ N 56 , a K 4 = 1 4|K| - --→ N 12 + --→ N 34 - --→ N 56 , a K 5 = 1 4|K| - --→ N 12 - --→ N 34 + --→ N 56 , a K 6 = 1 4|K| + --→ N 12 - --→ N 34 + --→ N 56 , a K 7 = 1 4|K| + --→ N 12 + --→ N 34 + --→ N 56 , a K 8 = 1 4|K| - --→ N 12 + --→ N 34 + --→ N 56 .
As a consequence, the stiness coecients of the matrix

A T read                          (A T ) sisi = K∈Ms i 8 j=1 |K|Λ K a K i • a K j , if s i ∈ T \ T D , (A T ) sisj = K∈Ms i ∩Ms j |K|Λ K a K i • a K j , if s i = s j ∈ T \ T D , (A T ) sisi = 1, if s i ∈ T D , (A T ) sisj = 0, if s i ∈ T \ T D , and s j ∈ T D , (A T ) sisj = 0, if s i ∈ T D , and s i = s j ,
where Λ K is approximated by its value at the centre of the cell for the numerical computations. The above strategy has been used in the implementation of the scheme in Section 3.

Remark 1. The stencil of the 3D face interpolated scheme is small and compact.

In other words, each row of the matrix consists of at least one non-zero element and at most 27 non-zero elements in the case of Cartesian meshes. In the latter case, the diagonal element corresponds to an internal vertex, which is connected to its 26 surrounding vertices.

As a result of the following discrete Poincaré's inequality, the semi-norm • 1,2 turns out to be a norm on the subspace R #T 0 .

Lemma 1 There exists a constant Cp depending only on the diameter of Ω and the mesh regularity such that

v T 0,2 ≤ Cp v T 1,2 , ∀v T ∈ R #T 0 .
Proof The guidelines of the proof are sketched below. Let s i ∈ K. Observe that

∇ K v T = 8 j=1,j =i a K j vs j -vs i .
The regularity of the mesh claim the existence of C 1 , C 2 depending only on θ T such that

C 0 h 2 K ≤ |a K i • a K j | ≤ C 1 h 2 K (7)
Then, it can be checked that

v T 2 1,2 ≥ C 3 si∈T K∈Ms i 8 j=1,j =i |K| h 2 K us i -us j 2 ,
holds true for some C 3 independent of the mesh. Adapting of the proofs established in [5, Lemma 3.3], [14, Lemma 3.1] (in the case of the vertex-centred TPFA scheme) to our setting, one nds

v T 2 0,2 ≤ C 4 si∈T K∈Ms i 8 j=1,j =i h 2 K |us i -us j |(|us i | + |us j |).
where C 4 depends only on θ T and the diameter of Ω. Following the arguments of [28, Lemma 2.1] and bearing in mind [START_REF] Ciarlet | The nite element method for elliptic problems[END_REF] we nally deduce that

v T 2 0,2 ≤ C 5 si∈T K∈Ms i 8 j=1,j =i |K| h 2 K us i -us j 2 ≤ C 2 p v T 2 1,2 ,
which concludes the proof.

In the rest of this section we prove that the numerical scheme is uniquely solvable. First, the following result is referred to as the coercivity property.

Proposition 2 Let u T be a solution to the scheme [START_REF] Andreianov | Discrete duality nite volume schemes for LerayLionstype elliptic problems on general 2D meshes[END_REF]. The norm of the discrete gradient of u T is bounded by a constant Ce depending only on, f , λ 0 , λ 1 , g and the mesh regularity i.e.,

u T 1,2 ≤ Ce.
Proof Let us consider g T ∈ R #T g such that (g T )s = gs for all s ∈ T . Let us also take [START_REF] Andreianov | On 3D DDFV discretization of gradient and divergence operators. I. Meshing, operators and discrete duality[END_REF]. Then, one writes

ϕ T = u T := u T -g T ∈ R #T 0 in
Z 1 = Z 2 + Z 3 ,
where

Z 1 = Ω Λ h M ∇ h u T • ∇ h u T dx, Z 2 = Ω Λ h M ∇ h g T • ∇ h u T dx, Z 3 = Ω f I h u T dx. Because Λ h M is coercive there holds Z 1 ≥ λ 0 u T 2 1,2 .
It is known that, there exists C 1 depending only on ξ 0 such that

g T 1,2 ≤ C 1 g H 1 (Ω) . (8)
Then, using (8) together with the Cauchy-Schwarz and Young's inequalities we nd

|Z 2 | ≤ λ 0 4 u T 2 1,2 + λ 1 (C 1 ) 2 λ 0 g 2 H 1 (Ω) .
Applying the Cauchy-Schwarz, Poincaré's and Young's inequalities allows us to obtain the estimation

|Z 3 | ≤ λ 0 4 u T 2 1,2 + C 2 p f 2 L 2 (Ω) λ 0 .
Consequently, one obtains u T 1,2 ≤ C 2 , for some constant C 2 depending only on the data and ξ 0 . We nally deduce

u T 1,2 ≤ g T 1,2 + u T 1,2 ≤ Ce.
The proof is then concluded.

Proposition 3 The numerical scheme [START_REF] Andreianov | Discrete duality nite volume schemes for LerayLionstype elliptic problems on general 2D meshes[END_REF] has a unique solution.

Proof It suces to prove that the kernel of

A T is trivial. Let X T ∈ R #T such that A T X T = 0.
Let us prove that X T = 0. By the denition of A T , one has (X T )s = 0 for all s ∈ T D . Thus, X T ∈ R #T 0 . Recall that the formulations ( 5)-( 6) are equivalent. It follows from the computations

0 = A T X T • X T = s∈T \T D K∈Ms |K|Λ K ∇ K u T • Xs∇ K es = Ω Λ h M ∇ h X T • ∇ h X T dx ≥ λ 0 X T 1,2 ,
where we have again used the fact that Λ h M is coercive. Hence, it must be that X T 1,2 = 0, implying X T = 0, because X T 1,2 is a norm on R #T 0 . This concludes the proof.

Mimicking the previous proof entails that the matrix A T is symmetric positivedenite.

Computational results

In this section, several test cases are simulated to assess the computational performance of the proposed numerical scheme. The rst objective is to test the accuracy of the developed 3D face interpolated scheme and investigate its robustness with respect to the anisotropy and the distortion of the mesh. The second aim is to illustrate the behaviour of the iterative solver in the presence, and the absence of the preconditioner. The computational domain Ω is the unit cube [0, 1] 3 . It is covered using four dierent partitions. They are referred to as the Cartesian, Kershaw, uctuated and sinusoidal meshes respectively, see Figure 4. The rst mesh is regular with cubic cells. The second mesh is taken from [START_REF] Eymard | 3D benchmark on discretization schemes for anisotropic diusion problems on general grids[END_REF][START_REF] Kershaw | Dierencing of the diusion equation in lagrangian hydrodynamic codes[END_REF]. The last two meshes are generated from the Cartesian mesh by transforming the old coordinates (x, y, z), that do not belong to the boundary edges, into the new positions (x , y , z ) dened as follows for the uctuated mesh. The numerical scheme is implemented in Fortran where the resulting linear system is solved using the conjugate gradient method. Its tolerance is xed to 10 -10 and the stopping criterion is taken on the norm of the successive iterates. We employ the SSOR preconditioner to improve the solver behaviour. Recall that the stiness matrix is denoted by A T . Let D be the diagonal component and E the strictly lower triangular part of the stiness matrix, so that A T = D + E + E T . We consider the SSOR preconditioner, with relaxation parameter set as 1, given as follows

P = (D + E)D -1 (D + E T ). (9)
We apply the preconditioner on the left of the stiness matrix. Numerically, this choice without parametrization was found to give good results after the extensive numerical investigation performed in this study. Some indicators are needed in order to quantify the performance of the solver. We denote by κ the condition number of A T , and #it the total number of iterations of the conjugate gradient method. The biggest eigenvalue is denoted by λ. Let CPUt denote the computation time measured in seconds, before preconditioning. The same quantities are denoted by κ , #it , CPUt , λ max after applying the preconditioner.

Five test cases are proposed below to test the performance of the numerical scheme in a progressive increase of the conguration severity: Case a : homogeneous and isotropic conguration, Case b : homogeneous conguration, but with an anisotropy ratio of 1000 along the z-axis, Case c : homogeneous and anisotropic conguration, Case d : anisotropic conguration, with material directions evolving in space, Case e : strongly discontinuous diusion matrix (piece-wise variations).

In all cases, an analytical solution is imposed and the function f is set for this analytical solution to be a solution of problem [START_REF] Aavatsmark | Discretization on unstructured grids for inhomogeneous, anisotropic media. Part II: Discussion and numerical results[END_REF]. The accuracy of the scheme is assessed using the corresponding manufactured analytical. For this purpose, we evaluate the discrete L 2 and H 1 relative errors:

u e,T -u T 0,2 u e,T 0,2 , u e,T -u T 2 0,2 + u e,T -u T 2 1,2 u e,T 2 0,2 + u e,T 2 1,2
,

where u e,T is the interpolate (value) of the exact solution on the vertices.

Test case a

In this rst test case involving a homogeneous medium, the diusion tensor is taken as the identity matrix. The Dirichlet boundary conditions are computed according to the analytical solution u e (x, y, z) = sin(2πx) sin(2πx) sin(2πz).

The right hand side f is f (x, y, z) = 12π 2 sin(2πx) sin(2πx) sin(2πz).

Figure 5 presents the log-log plots of the errors in terms of the mesh size. It is shown that the convergence rates of the method are second order for both norms.

Smaller convergence rates on the gradients are noticed on the rst elements of the Kershaw meshes. The solver behaviour is given in Table 2. The condition number as well as the total number of iterations are greatly reduced after the preconditioning. This translates into gains in the CPU time. Interestingly, the iterations increased after preconditioning for the Cartesian mesh. We now consider a global, highly anisotropic tensor in the z-direction

Λ =   1 0 0 0 1 0 0 0 1000   .
The exact solution is given by

u e (x, y, z) = 64x(1 -x)y(1 -y)z(1 -z).
The source function corresponding to this solution is given by

f (x, y, z) = 128 y(1 -y)z(1 -z) + x(1 -x)z(1 -z) + 1000x(1 -x)y(1 -y) .
A homogeneous Dirichlet boundary condition is prescribed. The obtained results are displayed in Figure 6. The order of convergence is quadratic regardless of the anisotropy. Table 3 summarises the solver behaviour in the presence of a highly anisotropic tensor. Signicant improvements in convergence rates of the iterative solver are observed after the application of the preconditioner, especially when the mesh is nonorthogonal and distorted.

Test case c

In this example, we look at the behaviour of the discretisation scheme in the case where diusion tensor is a fully heterogeneous and anisotropic matrix

Λ =   y 2 + z 2 + 1 -xy -xz -xy x 2 + z 2 + 1 -yz -xz -yz x 2 + y 2 + 1   .
This tensor is inspired from [START_REF] Lipnikov | The mimetic nite dierence discretization of diusion problem on unstructured polyhedral meshes[END_REF]. The analytical solution is set to u e (x, y, z) = x 2 yz 3 + x cos(πy) cos(πz),

where the right-hand side and nonhomogeneous Dirichlet boundary condition agree with this function. The numerical results are exhibited in Figure 7. The accuracy is again of second order for the solution and the gradient. This conrms the robustness of our scheme with respect to the mesh and anisotropy. Table 4 summaries the solver behaviour. Here again, the preconditioned system is better in terms of performance in the case of a heterogeneous tensor. In this experiment, the eigenvectors directions change over the domain. Then, the diusion tensor is decomposed as follows [START_REF] Guillaume | Numerical convergence of a parameterisation method for the solution of a highly anisotropic two-dimensional elliptic problem[END_REF][START_REF] Hermeline | A nite volume method for approximating 3d diusion operators on general meshes[END_REF]

Λ(x, y, z) = R θ (x) ×   1 0 0 0 ε 0 0 0 η(1 + x + y + z)   × R θ (x) t ,
where ε = 0.1, η = 10 and

R θ (x) =   cos(πx) -sin(πx) 0 sin(πx) cos(πx) 0 0 0 1   .
We impose the exact solution to be u e (x, y, z) = sin(πx) sin(πx) sin(πz).

The source term is routinely calculated from this solution. A homogeneous boundary condition is applied. The convergence results are plotted in Figure 8. Here we again observe the numerical optimal convergence rates in the L 2 -norm and H 1 -norm for the dierent meshes. The solver behaviour shown in Table 5 is consistent with the previous test case. We retain the impact of the preconditioner on reducing the CPU time, whereas the unpreconditioned solver requires considerably more iterations. This, for instance, is clearly apparent for the case of the Kershaw mesh.

Test case e

In this last example, we test our scheme in the case of a highly discontinuous diusion matrix. We divide the domain Ω into two sub-domains Ω 1 and Ω 2 where

Ω 1 = [0, 1] × [0, 0.5] × [0, 1], Ω 2 = [0, 1] × [0.5, 1] × [0, 1].
Following [START_REF] Aavatsmark | A new nite-volume approach to ecient discretization on challenging grids[END_REF], the anisotropy tensor is respectively dened on Ω 1 and Ω 2 by where ≥ 0. It can be checked that

Λ 1 =   + 1 + 1 + 1   , Λ 2 =   + 1 - - + 1 - - + 1   , Cartesian meshes h T κ κ λmax λ max #it #it CPU CPU 0.
u e (x, y, z) = 2y 3 -3y( + 1)(x -z) 2 ,
is the exact solution to the model with f = 0. The Dirichlet boundary condition is nothing more than the trace of u e on ∂Ω. The ratio of anisotropy is given by 3 + 1. In Figure 9, we exhibit the errors induced by the numerical scheme in the L 2 -norm for ∈ {33, 3333} leading to the ratio of anisotropy which is equal to 100 and 10000 respectively. One can observe that the convergence rate is of second order in both cases, independently of the last three meshes, and the scheme is exact on the Cartesian mesh.

In Table 6, a good performance of the new solver is observed. The reduction in the number of iterations as well as the CPU time is spectacular on the Kershaw mesh, which is a challenging mesh and accounts for a severe distortion. A similar behaviour is observed in the case where = 333.

4 Application to the prediction of wood thermal conductivity from its real 3D morphology

The objective here is to predict the eective thermal conductivity of wood using upscaling. One type of wood, namely poplar species is studied. For this purpose, we make use of high resolution scans provided by the Laboratory nano-tomograph and the proposed 3D face interpolated scheme to estimate the macroscopic property following the orthotropic directions of wood. For more details on the description of the sample preparation and the image processing of the scanned morphology, we refer to [START_REF] Louërat | Heat and moisture diusion in spruce and wood panels computed from 3-d morphologies using the lattice boltzmann method[END_REF][START_REF] Quenjel | Computation of the eective thermal conductivity from 3d real morphologies of wood[END_REF]. An illustration of the scanned sample is exhibited in Figure 10. It is referred to as Morpho . Let us stress that it is quite dicult to make predictions on this whole volume. A standard practical alternative is to extract sub-volumes from the original morphology. Such a subset is called the REV (Representative Elementary Volume). The representativeness of the volume with respect the thermal conductivity is underlined below according to the convergence test.

We denote by Ω = REV = (x 0 , x m )×(y 0 , y m )×(z 0 , z p ) ⊂ Morpho the extracted region of interest. Its centre and dimensions are specied by the user. Then, one performs the automatic segmentation so that one obtains the binary 3D image and generate the orthogonal Cartesian mesh made from cubes. The thresholded morphology contains only two phases consisting of the air and the solid matrix as the left side of Figure 11. Its right side shows the sinusoidal perturbation of the real morphology. This modication intends to test the eciency and the robustness of the 3D face interpolated scheme in predicting the thermal conductivity on a distorted sample. The next step is to discretise Fourier's law on the morphology mesh thanks to the proposed 3D face interpolated scheme. The model problem writes

-div α(x)∇T ) = 0, in Ω = (x 0 , x m ) × (y 0 , y n ) × (z 0 , z p ), T = T bd , on ∂Ω D , α(x)∇T • n = 0, on ∂Ω N ,
where ∂Ω = ∂Ω D ∪ ∂Ω N is the boundary and α = λ 0 /( c p ) is the local thermal coecient of each voxel (mesh vertex). It is composed of λ 0 the local thermal conductivity, the density and c p the heat capacity. For computation, we take c p = 1. The λ 0 is xed to 0.5 W.m -1 .K -1 for the solid phase and 0.023 W.m -1 .K -1 for the air phase.

As we are going to run several computations, we select an increasing sequence of REVs, namely REV i = (x i 0 , x i m ) × (y i 0 , y i m ) × (z i 0 , z i p ) where

x i 0 = y i 0 = z i 0 = 0, x i m = y i n = (32 × i) -1 × Resolution, z i p = ((64 × i) -1) -1 × Resolution, for i = 1, • • • , 5.
The resolution is given by 3.10 -6 m. For instance, the mesh of the rst REV is made from 32 × 32 × 63 vertices and so on. The resulting linear system is solved by the preconditioned conjugate gradient method. The behaviour of the temperature gradient is highlighted in Figure 12. On the both meshes, it is clearly seen that the heat follows naturally the cell walls. Accordingly, one can deduce the predicted thermal conductivity by assessing the ratio of the boundary ux Flux bd to the global imposed gradient δT /Long imposed by the taken boundary conditions. As a consequence, the radial, tangential and longitudinal directions conductivities are derived from

λ t R = Flux bd x × |x m -x 0 | δT , λ t τ = Flux bd y × |y n -y 0 | δT , λ t L = Flux bd z × |z p -x 0 | δT .
The results are provided in Table 7. It reports the predicted macroscopic property in terms of the REV size. This convergence test shows that the property oscillates on the rst volumes whereas it becomes to be stabilized as the REV size increases. Moreover, one can say that the REV5 is representative. On both meshes, the longitudinal value is bigger than the radial and tangential ones, which is in coherence with the wood structure. We also observe a slight reduction on the conductivity on the transformed morphology, which goes back to the impact of the waves induced by the sinusoidal perturbation. The results are similar to the ones obtained in [START_REF] Quenjel | Ecient prediction of the thermal conductivity of wood from its microscopic morphology[END_REF][START_REF] Quenjel | Computation of the eective thermal conductivity from 3d real morphologies of wood[END_REF]. However, in these references the method is the standard TPFA nite volume which works only on orthogonal meshes and does not support anisotropy. This ensures the robustness of the proposed 3D face interpolated method in capturing the ranges of the thermal conductivity even on a distorted mesh.

Mesh

Conclusion

In this paper, a new numerical method, referred to as the 3D face interpolated scheme, is proposed for the discretisation of highly anisotropic diusion equations on hexahedral meshes. The idea consists in dening an accurate and consistent approximation of the gradient using the opposite faces of the cells. The face unknowns are eliminated by employing a second order interpolation. This reduces greatly the complexity of the scheme by placing the unknowns at the mesh nodes. Because the scheme is derived from a discrete weak formulation, it is unconditionally coercive and therefore admits a unique solution. As a result, the stiness matrix is symmetric positive-denite. The numerical implementation makes use of the practical SSOR preconditioner in order to improve the behaviour of the iterative solver. The computational results show that the method is robust and attains second order accuracy independently of the chosen meshes and the anisotropy. After the validation step, the approach is applied for predicting the wood thermal conductivity in the species natural directions. Good results are obtained especially when the real morphology is perturbed using a sinusoidal transformation. In our future research we intend to extend this new methodology in a homogenisation framework for predicting the eective parameters of bres with more than two phases that arise in volume-averaged transport equations.

Figure 1 :

 1 Figure 1: Sub-set associated with the vertex s in the cell K.

Figure 2 :

 2 Figure 2: opposite faces in the cell K and their corresponding directions.

Figure 3 :

 3 Figure 3: Numbering of the vertices (left) and the faces (right) of the cell K.

x

  = x + 1 10 ϕ(x)ϕ(y), y = y + 1 10 ϕ(y)ϕ(z), z = z + 1 10 ϕ(x)ϕ(z), where ϕ(r) = sin(2πr) in the case of the sinusoidal mesh and ϕ(r) = 4r1 [0,0.5] + (2 -4r)1 [0.25,0.5] + (4r -2)1 [0.5,0.75] + (4 -4r)1 [0.75,1] ,

Figure 4 :

 4 Figure 4: Cartesian, Kershaw, uctuated and sinusoidal meshes.

Figure 5 :

 5 Figure 5: Test a : numerical relative errors in the L 2norm (left) and H 1 -norm (right).

Figure 6 :

 6 Figure 6: Test b : numerical relative errors in the L 2norm (left) and H 1 -norm (right).

Figure 7 :

 7 Figure 7: Test c : numerical relative errors in the L 2norm (left) and H 1 -norm (right).

Figure 8 :

 8 Figure 8: Test d : numerical relative errors in the L 2norm (left) and H 1 -norm (right).

Figure 9 :

 9 Figure 9: Test e: numerical relative errors in the L 2norm (left) and H 1 -norm (right).

Figure 10 :

 10 Figure 10: Scanned poplar sample using tomography with a resolution of 3µm.

Figure 11 :

 11 Figure 11: Original binary Cartesian morphology (left) and the corresponding sinusoidal perturbation (right)

Figure 12 :

 12 Figure 12: Gradient distribution of temperature on the Cartesian morphology (left) and on the sinusoidal one (right).

Table 1

 1 reports the number of degrees of freedom, dimension and sparsity of the stiness matrix for each mesh level.

	Mesh Nb Nb of DOFs	Nb NzCoefMat	Density
	1 2 3 4 5	729 4913 35937 274625 2097152	15625 117649 912673 7189057 57066625	2.94011 % 0.48741 % 0.07066 % 0.00953 % 0.00129 %

Table 1 :

 1 DOFs for each mesh and the corresponding number of non zero coecients of the stiness matrix.

Table 2 :

 2 Test a: solver statistics before and after the preconditioning.

	Cartesian meshes
	h T 0.125 8.30E+2 4.85E+1 1.00E+0 1.10E+0 18 κ κ λmax #it #it λ max 12 0.063 2.61E+4 7.58E+2 1.00E+0 1.12E+0 21 11 0.031 8.31E+5 1.19E+4 1.00E+0 1.13E+0 40 20	CPU CPU 0.0059 0.0044 0.0546 0.0417 0.4577 0.3482
	Kershaw meshes
	h T 0.178 2.91E+3 8.98E+1 5.43E+0 1.12E+0 98 κ κ λmax #it #it λ max 40 0.089 1.20E+5 5.13E+2 6.05E+0 1.06E+0 419 89 0.045 2.90E+6 1.39E+3 5.72E+0 1.15E+0 1158 134	CPU CPU 0.0129 0.0088 0.1864 0.0986 3.1635 0.9512
	Fluctuated meshes
	h T 0.163 8.76E+2 5.28E+1 1.00E+0 1.08E+0 49 κ κ λmax #it #it λ max 24 0.085 2.60E+4 8.33E+2 1.00E+0 1.10E+0 107 37 0.043 8.31E+5 6.76E+3 1.00E+0 1.12E+0 77 27	CPU CPU 0.0071 0.0100 0.0941 0.0611 0.5714 0.4546
	Sinusoidal meshes
	h T 0.156 7.75E+2 4.80E+1 1.00E+0 1.06E+0 50 κ κ λmax #it #it λ max 23 0.081 2.11E+4 7.93E+2 1.00E+0 1.12E+0 42 16 0.041 6.75E+5 7.86E+3 1.00E+0 1.13E+0 82 24	CPU CPU 0.0096 0.0072 0.0760 0.0574 0.5812 0.4406

3.2 Test case b

Table 3 :

 3 Test b: solver statistics before and after the preconditioning.

Table 4 :

 4 Test c: solver statistics before and after the preconditioning.

	Cartesian meshes
	h T 0.125 8.07E+2 4.57E+1 1.00E+0 1.11E+0 16 κ κ λmax #it #it λ max 9 0.063 1.62E+4 7.22E+2 1.00E+0 1.14E+0 23 8 0.031 5.15E+5 7.66E+3 1.00E+0 1.16E+0 29 12	CPU CPU 0.0274 0.0153 0.1933 0.1067 1.4337 0.9794
	Kershaw meshes
	h T 0.178 3.08E+3 8.41E+1 8.81E+0 1.07E+0 104 39 κ κ λmax #it #it λ max 0.089 1.22E+5 5.33E+2 1.00E+1 1.17E+0 441 81 0.045 3.71E+6 1.18E+3 1.03E+1 1.20E+0 352 47	CPU CPU 0.0399 0.0191 0.3149 0.1235 2.3797 1.0002
	Fluctuated meshes
	h T 0.163 1.17E+3 5.49E+1 1.26E+0 1.10E+0 48 κ κ λmax #it #it λ max 19 0.085 2.44E+4 8.41E+2 1.00E+0 1.10E+0 44 16 0.043 4.99E+5 4.42E+3 1.00E+0 1.14E+0 38 13	CPU CPU 0.0290 0.0111 0.2095 0.1008 1.6345 1.0540
	Sinusoidal meshes
	h T 0.156 1.10E+3 5.14E+1 1.37E+0 1.10E+0 37 κ κ λmax #it #it λ max 15 0.081 2.05E+4 8.28E+2 1.02E+0 1.14E+0 25 11 0.041 4.54E+5 6.22E+3 1.00E+0 1.17E+0 36 13	CPU CPU 0.0400 0.0260 0.1977 0.0981 1.5419 1.0680

3.4 Test case d

Table 5 :

 5 Test d: solver statistics before and after the preconditioning.

	125 2.47E+3 8.74E+1 1.22E+1 1.01E+0 66 0.063 2.69E+4 6.68E+2 7.43E+0 1.02E+0 46 0.031 4.53E+5 6.96E+3 4.17E+0 1.05E+0 74	27 18 14	0.0254 0.0153 0.1849 0.0905 1.4793 1.0899
	Kershaw meshes
	h T 0.178 3.35E+3 5.72E+1 2.42E+1 1.10E+0 145 κ κ λmax #it #it λ max 43 0.089 6.90E+4 2.22E+2 2.46E+1 1.05E+0 618 98 0.045 1.59E+6 1.08E+3 1.73E+1 1.12E+0 1017 111	CPU CPU 0.0321 0.0209 0.3600 0.1235 3.8786 1.1868
	Fluctuated meshes
	h T 0.163 3.49E+3 6.30E+1 1.66E+1 9.95E-1 κ κ λmax λ max 0.085 7.39E+4 6.56E+2 1.32E+1 1.11E+0 342 #it #it 131 37 64 0.043 1.41E+6 2.02E+3 8.68E+0 1.05E+0 284 44	CPU CPU 0.0377 0.0211 0.2631 0.1088 2.1575 1.0610
	Sinusoidal meshes
	h T 0.156 5.63E+3 6.63E+1 1.92E+1 1.08E+0 145 κ κ λmax #it #it λ max 32 0.081 9.43E+4 5.74E+2 1.53E+1 1.05E+0 152 31 0.041 1.69E+6 1.96E+3 1.01E+1 1.13E+0 162 26	CPU CPU 0.0294 0.0172 0.2138 0.1063 1.7939 0.9587

Table 6 :

 6 Test e: solver statistics before and after the preconditioning for = 33.

Table 7 :

 7 Predicted thermal conductivity of poplar with the real and the transformed morphology following the peculiar orthotropic directions of wood. Table8indicates the solver performance. The rst column lists the REV number. The second one informs the CPU time of the solver required to compute the three conductivities λ t R , λ t τ and λ t L . The last one provides the maximum of the errors committed in evaluating the latter values i.e.,Max residual = max Errorλ t R , Error λ t τ , Error λ tOne infers that the method is ecient and accurate on the Cartesian morphology, because the phases are quite well placed in parallel and ranged following the longitudinal direction compared to sinusoidal sample.

		λ t (W/(m.K))	REV1 REV2 REV3 REV4 REV5
	Cartesian	λ t R λ t τ	0.1906 0.1705 0.1430 0.1506 0.1522 0.1841 0.1466 0.1336 0.1254 0.1314
		λ t L	0.2507 0.2238 0.2009 0.2021 0.2046
	Sinusoidal	λ t R λ t τ	0.1814 0.1657 0.1317 0.1418 0.1441 0.1850 0.1511 0.1351 0.1266 0.1338
		λ t L	0.2386 0.2153 0.1864 0.1876 0.1920

L

.

Table 8 :

 8 Solver performance.

Acknowledgments: This study was carried out in the Centre Européen de Biotechnologie et de Bioéconomie (CEBB), supported by the Région Grand Est, Département de la Marne, Greater Reims and the European Union. In particular, the author would like to thank the Département de la Marne, Greater Reims, Région Grand Est and the European Union along with the European Regional Development Fund (ERDF Champagne Ardenne 2014-2020) for their nancial support of the Chair of Biotechnology of CentraleSupélec.


Jérôme Dubail

Jesper Lykke Jacobsen

Hubert Saleur

Conformal field theory at central charge c = 0: a measure of the indecomposability (b) parameters

A good understanding of conformal field theory (CFT) at c = 0 is vital to the physics of disordered systems, as well as geometrical problems such as polymers and percolation. Steady progress has shown that these CFTs should be logarithmic, with indecomposable operator product expansions, and indecomposable representations of the Virasoro algebra. In one of the earliest papers on the subject, V. Gurarie introduced a single parameter b to quantify this indecomposability in terms of the logarithmic partner t of the stress energy tensor T . He and A. Ludwig conjectured further that b = -5 8 for polymers and b = 5 6 for percolation. While a lot of physics may be hidden behind this parameter -which has also given rise to a lot of discussions -it had remained very elusive up to now, due to the lack of available methods to measure it experimentally or numerically, in contrast say with the central charge. We show in this paper how to overcome the many difficulties in trying to measure b. This requires control of a lattice scalar product, lattice Jordan cells, together with a precise construction of the state L-2|0 . The final result is that b = 5 6 for polymers. For percolation, we find that b = -5 8 within an XXZ or supersymmetric representation. In the geometrical representation, we do not find a Jordan cell for L0 at level two (finite-size Hamiltonian and transfer matrices are fully diagonalizable), so there is no b in this case.

Introduction

In the last twenty years or so since the seminal paper [1] of Belavin, Polyakov and Zamolodchikov, Conformal Field Theory (CFT) has proven amazingly successful. It is now an essential item in the toolbox of condensed matter and string theorists, and has had a profound impact on several sub-fields of modern mathematics.

Yet, despite this flurry of successes, some very fundamental questions have remained unanswered to this day. One of these questions concerns percolation, the very geometrical critical problem where CFT has obtained some of its most impressive results. In a nutshell, despite years of work and some progress (see below), we do not know any conformal field theory describing at least some of the geometrical observables (be it hulls, clusters or backbones) in a fully consistent way. As a result, many quantities involving bulk correlations functions-for instance the equivalent of the Binder cumulant [2] in this problem [3]-are, to this day, unknown analytically.

Another vexing question concerns the celebrated transition between plateaux in the integer quantum Hall effect. The evidence is strong that it corresponds to quantum critical points of the 2+1 dimensional electron gas. The physics of the transition is somewhat well understood (and involves an interplay between disorder, which tends to localize the electrons, and the kinetic energy, quenched by the strong magnetic field, which causes delocalization), and qualitatively well described by a 2D sigma model with topological term [4]. Nevertheless, a precise identification of the low energy effective field theory is still lacking to this day, despite the wealth of theoretical and numerical works on the topic, and the fact that this theory is expected to be conformally invariant (see [5] for an insightful review).

The origin of these difficulties is, ultimately, the fact that the CFTs describing these problems have to be non unitary, with vanishing central charge. Non-unitarity in a CFT can have rather mild consequences-like in the minimal theories such as the Yang Lee edge singularity. It took a while to realize that non-unitarity would most often (in particular, when c = 0) imply in fact indecomposability, leading to a very difficult kind of theory called a logarithmic CFT (LCFT).

LCFTs were probably first encountered in published form in a paper by Rozansky and Saleur [6]. In this paper, the authors studied a particular kind of c = 0 theory with U (1|1) supergroup symmetry, and stumbled upon four-point functions involving logarithmic dependence on the cross ratio. These authors correctly related this property to the indecomposability of the operator product expansions (OPE), and to the non-diagonalizability of the L 0 generator, inherited from the non semi-simplicity of the symmetry algebra. Shortly after, Gurarie [7] pointed out that these features were in fact necessary to have a consistent, non-trivial CFT at c = 0. Later, Gurarie [8] and Gurarie and Ludwig [9] built up a very attractive formalism within which CFTs at c = 0 must possess, in addition to their stress energy tensor T (z), an extra field whose holomorphic part, t(z), has conformal weight two. The singular part of the OPE between T (z) and t(z) is determined up to a new universal number, an "anomaly" usually denoted by b. This parameter is expected in [8,9] to play a very important physical role. It might obey a "c-theorem" [7] and thus indicate possible directions of RG flows within theories with vanishing central charge. Its value also dictates the existence of null vectors for conformal weights contained in the Kac table, and is thus profoundly related with the determination of four-point functions. Gurarie and Ludwig [9] suggested, based on knowledge of critical exponents from Coulomb gas calculations and some heuristic hypotheses, that b = - 5 8 for polymers and b = 5 6 for percolation.1 Despite the appeal of the ideas proposed in [8,9], only very little progress has happened to make them into a powerful theoretical tool. It is not clear for instance to what extent it is possible or helpful to study further the extension of the Virasoro algebra proposed in [9]. Considerable difficulties also arise when one tries to build non-chiral theories by combining the left and right moving sectors. Most examples worked out so far indeed concern only boundary LCFTs.

Meanwhile, the subject has matured with the understanding that, probably, rather than a single parameter b, LCFTs might be characterized by a complicated structure of indecomposable Virasoro modules. This structure involves maps between sub-modules, and several parameters to describe their precise action. This line of thought shows how close the problem is to the theory of non semi-simple Lie (super)algebras, whose representation theory is more often than not of the "wild" type. Progress has been steady in trying to understand indecomposable Virasoro modules [10,11,12,13,14] in general LCFTs. For c = 0 the b parameters appear indeed as some particular coefficient describing the embedding of a sub-module. Some simple algebraic arguments suggest then that b = - 5 8 for percolation and b = 5 6 for polymers, that is, the value of Gurarie and Ludwig [9] up to a switch.

In the last two years, more progress has occurred from the direction of lattice models. It had long been known that the representation theory of lattice algebras (mostly, the Temperley-Lieb algebra) bore some striking resemblance to that of the Virasoro algebra. In particular, it had been known that when q is a root of unity, the hamiltonian (the lattice discretization of L 0 , or L 0 + L0 depending on the boundary conditions) can sometimes be non-diagonalizable, and that Jordan cells appear, very much mimicking the ones expected in the continuum limit. This was explored in a much more systematic fashion in works by Pearce, Rasmussen and Zuber [15,16,17] and independently by Read and Saleur [18]. The outcome of these analyses was a rather coherent picture of the indecomposable Virasoro modules appearing in some families of boundary LCFTs, including percolation and polymers, together with some formal but potentially useful results about fusion. This picture was interpreted in interesting mathematical terms in [19].

Yet, large parts of these constructions are still speculative. For instance, the lattice analysis clearly shows which modules are mapped into each other under the action of the lattice algebras, but what this becomes precisely in the continuum limit-including which Virasoro generators and states are involved in these mappings, together with precise values of the corresponding coefficients-is largely conjectural. Moreover, it is clear that only the 'simplest indecomposable' modules have been encountered so far, and one is still far from being able to guess what would happen, say, for the quantum Hall plateau transition.

One of the crucial remaining obstacles is that the features that make LCFTs difficult are also very hard to observe or measure directly. It is quite striking for instance, that despite the 15 year old controversy around the value of the b parameter for percolation or polymers-problems which are usually rather easy to study numerically-no way to measure this parameter has been available, up to now. One of the only forays in this direction is the paper by Koo and Saleur [20], who attempted to define and study a straightforward regularization of the Virasoro algebra on the lattice. Their success was only partial, and hampered by the fact that the continuum limit of commutators is not the commutator of the continuum limits.

We shall report in this paper a method to measure the parameter b in lattice models. While it might not be general enough to work for all theories at c = 0, it certainly will allow us to assert what the values of this parameter for percolation and polymers actually are.

The paper is organized as follows. In the first section, we recall basic facts about the b number in LCFTs at c = 0, and define a general strategy to measure it. This strategy requires overcoming two difficulties. One is finding a scalar product on the lattice that goes to the Virasoro-Shapovalov form in the continuum limit. The other is finding a regularization of the L -2 operator-or more precisely, of its action on the vacuum. These difficulties are solved in the next two sections. In section 4, we finally gather all the pieces to study the Jordan cells for percolation and polymers, and measure b.

The b number

In a nutshell, the argument of Gurarie [7,8,9] went as follows. Conformal invariance fixes the coefficient of the stress tensor T in the OPE of an operator with itself (assuming there is a single field with h = 2, h = 0) to be of the form

φ(z)φ(0) ≈ a φ z -2h φ 1 + 2h φ c T (z) + . . . (1)
Here a φ is an amplitude, determining the normalization of the two-point function

φ(z, z)φ(0, 0) = |a φ | 2 1 z 2h φ z2 hφ (2)
If one want to keep a φ finite in a theory at c = 0 for an operator with h φ = 0, the 2h φ c factor poses problems as it diverges.

There are various ways to resolve the difficulty (this was also discussed in [21,[START_REF] Cardy | Logarithmic correlations in quenched random magnets and polymers[END_REF]. First, a divergence in the OPE coefficient (which would manifest itself in a physical quantity such as a four point function) might not be a problem after all, but have a physical meaning, related with some n → 0 limit. Second, it could be that demanding a φ finite is not physically meaningful, again in relation with an n → 0 limit. Recall that since we are dealing with non-unitary theories, having a vanishing two-point function does not mean that the corresponding operator is zero. The third way to resolve the difficulty is to admit that there are other operators with h = 2 appearing on the right-hand side of the OPE (1). These operators might be part of a supersymmetry multiplet, such as those that might occur in supergroup WZW models [START_REF] Schomerus | [END_REF] and other supersymmetric CFTs [24]. We note that in such cases, there is no need for a logarithmic dependence in the OPE at this order, although the theory will in general be logarithmic (see, e.g., [START_REF] Schomerus | [END_REF]. The other scenario proposed in [9] is of an n → 0 limit, with another operator whose dimension is generically different from 2 becoming degenerate with T . In this case-examples of which can be worked out in details for the models with central charge c = 1 -6

x(x+1) , in the limit x → 2-one expects the OPE of φ with itself to read instead

φ(z)φ(0) ≈ a φ z -2h φ 1 + 2h φ b z 2 [t(0) + log(z)T (0)] + . . . . (3)
The operator t is called 'logarithmic partner' of T , and forms with it a logarithmic pair, that is, a set of fields on which the L 0 operator is not diagonalizable but has the form of a Jordan cell

L 0 |T = 2|T , L 0 |t = 2|t + |T (4)
In terms of two-point functions one has

T (z)T (0) = 0 (5a) T (z)t(0) = b z 4 (5b) t(z)t(0) = -2b log z + a z 4 (5c)
where the last result comes from imposing global conformal invariance [9]. Note that the normalization of T is crucial in defining the numerical value of b, once a choice has been made that the L 0 Jordan cell has off-diagonal term unity (which translates into the relative normalization of the last two equations). There remains, on the other hand, a choice t → t + constant T , which does not affect the leading terms in the OPE, but only the constant a.

The universality of these OPEs, the value of the number b and its properties (e.g., under RG flows) have all been subject of intense debate. It is not our purpose to review this debate in detail. We only wish to recall that, when thinking of the limit of minimal models as c → 0, one is naturally led to two candidates which might become degenerate with T : the (3, 1) field and the (1, 5) field (in Kac table notations). The former choice leads to b = 5 6 , and the latter to b = - 5 8 . Using a variety of arguments, Gurarie and Ludwig [9] predicted that the former choice corresponds to percolation, and the latter to polymers, which are therefore profoundly different logarithmic CFTs.

While one can dwell at great length on the merits of the arguments in [9] we believe it is more constructive to ask whether one can, in fact, measure the parameter b. Some considerable difficulties arise when thinking of this possibility. One might, of course, try to identify t in a lattice model of percolation or polymers. After all, it is quite possible to get a handle on the meaning of T itself [25,[START_REF] Doyon | Conformal loop ensembles and the stress energy tensor[END_REF]. But measuring directly correlation functions such as those in (5) is a very daunting task, not to mention the difficulty of extracting properly a logarithmic dependence. We also note that the formalism in [9] is really a chiral one, where the job of glueing back together the left and right parts of the theory is left for further work, and might lead to considerable surprises. This chiral aspect naturally suggests turning to the physics of the boundary theory, which is in fact what is mostly considered in [9].

Short of measuring correlation functions on the lattice, the next best hope would be to access numerically b through a lattice analog of the Jordan cell structure in (4). Some progress in reproducing algebraic relations-e.g., the Virasoro commutation relations-in lattice models was reported in fact in [20], and it is a natural idea to try to push further the results obtained there. Introduce first the scalar product and the notion of adjoint from in and out states as usual [START_REF] Flohr | [END_REF], with L † n = L -n (of course now the form is not positive definite). Then the relations equivalent to (5) are

L 0 = 2 1 0 2 (6)
in the basis {T, t}, where

|T = L -2 |0 , T |T = 0 L 2 |t = b|0 , t|T = b . (7)
We see that a possible strategy to measure b would involve finding a Jordan cell for the hamiltonian (the lattice version of L 0 ) in finite size. One would then have to identify the states corresponding to T, t and normalize them properly. Finally one would have to calculate the scalar product t|T . We shall accomplish this program here, but first let us stress why it is difficult, and has not been done before.

The first obstacle is to find a Jordan cell in lattice models. This in fact is quite easy, as has been known for a while. The point is that the 'hamiltonians' in lattice models whose continuum limit is a CFT are often non-diagonalizable, a consequence of their deep relationship with lattice symmetries (quantum groups) or lattice algebras (Temperley-Lieb algebras) which have non semisimplicity properties fully emulating those of the continuum limit. Nevertheless, the question whether there exists a Jordan cell for L 0 in a given lattice model is not easily answered a prioriwe shall mention some surprises along these lines later.

Having obtained the Jordan cell, one needs to identify T . While it is natural to think that T is an eigenvector of the hamiltonian, a very important obstacle is its normalization. While crucial to determine the value of b, it cannot be done via a standard eigenvector normalization since the two-point function of T vanishes in the continuum, which manifests itself by the fact that |T -or rather its lattice equivalent-is a null state. To normalize T , we shall resort to a very indirect trick that we call the "trousers trick", which is really at the heart of our solution of the problem.

Finally, one needs to properly normalize t via the desired Jordan cell structure. This is quite involved with big transfer matrices, and requires some numerical tricks. We note in passing that, according to (5 c), |t turns out to have infinite norm square.

Lattice scalar products 2.1 Hamiltonians and Transfer Matrices

A 2D classical statistical model can be viewed as a 1+1D model evolving in imaginary time. At large length scales, the low-energy excitations of the latter are described by some quantum field theory. When the 2D model is at a critical point, these excitations are gapless, hence the field theory becomes conformal. A 1+1D system of width L is defined by its Hamiltonian H L , or by its transfer matrix T L when it comes from a 2D statistical model.

Let us start by considering the Ising model as an example. The critical Ising chain is defined by the Hamiltonian (say with periodic boundary conditions σ z L+1 = σ z 1 )

H L = - L i=1 σ z i σ z i+1 - L i=1 σ x i (8)
acting on the space of spin configurations |↑↓↓↑ . . . ↓ . H L is hermitian, therefore in the continuum limit is described by a unitary CFT. It is well-known that it corresponds to the minimal model M 4,3 with central charge c = 1/2. Now, our second example will turn out to be non-unitary. It is a dense loop model, given by the transfer matrix (in the picture L = 6)

T L = (9)
where each square is a sum over two terms + ≡ and one takes periodic boundary conditions (dotted lines are identified). Iterations of T L (from bottom to top) build the partition function of the dense loop model on a cylinder.

We have not specified yet on which configuration space T L is acting. Here the basis states will encode the planar pairings of L points, such as | or | for L = 6. A state in this model is a linear combination of such basis states. The action of T L on basis states is found by stacking the diagrams from bottom to top, e.g.:

| = = |
Each closed loop gets a Boltzmann weight n. For our formalism to be complete, we need to define a transposition operation (or equivalently a scalar product).

Transposition and scalar product

How can we define the forms .| ? In our geometrical formulation, the elements of the transfer matrix should act on such an object from top to bottom:

| = = n |
The loop scalar product is then defined for basis states by glueing the mirror image of the first state on top of the second one, giving a weight n per closed loop. For instance

| ≡ | = = n 2 (10)
This definition is extended by linearity to all the other states. Note that |l ≡ | -| has norm squared l|l = 2n2 (n -1), which is negative if n < 1. So for generic n the loop scalar product is not positive definite. There is no reason to expect otherwise, as the theories we are dealing with are not unitary. We note that this scalar product is well-known in the theory of the Temperley-Lieb algebra [28].

For the Ising chain (8), the scalar product is much simpler. It is the usual (Euclidian) one

L s 1 . . . s L |s ′ 1 . . . s ′ L L = δ s1,s ′ 1 . . . δ sL,s ′ L where s i , s ′ i ∈ {↑, ↓}.
In general, if there is a unitary representation of the model, one can find a basis in which the scalar product above is just the Euclidean one. For example, in minimal models the RSOS height representation provides an orthonormal basis |h L = |{h i } L , so writing the states as |a L = h a h |h L the scalar product is simply L a|b L = h a h b h . In non-unitary models, the simple loop model above shows that states and scalar product can be more subtle.

At a critical point the Hamiltonian can be related to the CFT translation operator along the cylinder L 0 + L0c/12 [29,30]

H L = E 0 L + 2πv F L L 0 + L0 - c 12 + O L -2 , ( 11
)
where v F denotes the Fermi velocity in the dispersion relation E ∼ v F |k|, and E 0 L is the (nonuniversal) extensive part of the ground state energy. A similar relation holds for log T L , since T L ≃ e -aHL , where a is the lattice spacing.

In finite size, when the Hamiltonian (or equivalently the transfer matrix) is diagonalizable, one can label the right eigenstates as

|0 L , |1 L , |2 L , . . . with E 0,L ≤ E 1,L ≤ E 2,L ≤ . . . .
The knowledge of these energies for different sizes L is a very useful way of extracting the conformal spectrum using (11) [29,30]. We claim in this paper that we can go further than the spectrum, and that scalar products can be measured as well. Let us introduce the left eigenstates of H L : L 0| , L 1| , L 2| , . . . . We expect these different states to correspond to the conformal states |0 , L -2 |0 , L-2 |0 , L -2 L-2 |0 , . . . (of course, descendants of other primary operators are also present in the spectrum, depending on the model we are looking at). We note that for different energies E i = E j (or even different momenta) one has L i|j L = 0, like in the continuum limit where scalar products between states at different levels are zero:

0|L k ′ 1 . . . L k ′ m L -kn . . . L -k1 |0 = 0 if m i=1 k ′ i = n i=1 k i .
In general, one can mimic on the lattice the construction of in and out states and scalar products of the continuum limit [START_REF] Ginsparg | Fields, Strings and Critical Phenomena[END_REF][START_REF] Flohr | [END_REF], and the scalar product we have defined appears as its only natural lattice regularization. We believe however that this might be too naive. In particular, an arbitrary state of the lattice model expands on an infinite sum involving scaling states (states associated with low energy eigenvalues) as well as non-scaling ones. Depending on the quantities one calculates, these non-scaling states might add up to non-vanishing contributions, even in the limit of a large system.

In the next section we give a strong numerical check that the lattice scalar products defined above for the Ising chain and for the dense loop model go over to the continuum limit ones as naively expected, for all quantities we shall be interested in.

Measure of Affleck-Ludwig boundary entropy 2

For a given CFT, different conformal boundary conditions exist [START_REF] Cardy | [END_REF]36]. When one perturbs a conformal boundary condition (CBC) by a relevant operator, it flows towards another CBC under the RG flow. These CBCs and their flows can be characterized by their boundary entropy [31]. A CBC can be encoded by a boundary state [36]. The scalar product of a boundary state |B with the ground state |0 of the conformal Hamiltonian on the cylinder L 0 + L0c/12 gives the boundary entropy s B =log B|0 . These numbers are universal and have been computed analytically for many CFTs and for many different CBCs. As a first application of the above considerations on finite-size scalar products, we are are going to recover them (at least numerically) from finite-size calculations in our two examples: the Ising chain and the dense loop model.

Affleck

f B s (L) = -log ( L B|0 L ).
When L → ∞ we expect this to scale as

f B s (L) ≃ f B s L + s B + O L -1 (12)
where f B s is the surface free energy per unit lenght, and s B is a universal term predicted by CFT. It is well-known [31,36] that s fixed =log √ 2/2 ≃ 0.34657 and s free = 0.

Numerically, we compute f B s (L) for L = 12, 14, 16, 18 spins. Then, writing

f B s (L) = 3 k=-1 c k L k ,
s B should be identified with c 0 according to (12). We find s fixed ≃ 0.3471 and s free ≃ 0.00008. This shows that the (universal part of the) finite-size scalar products L B|0 L can indeed be identified with the CFT result B|0 . This precise result for the Ising model is not new, and has deep links with calculations of Renyi entanglement entropies for critical wave functions [40].

Affleck-Ludwig term for loop models. To further check this method, we measure the boundary entropy of dense loops. When -2 < n ≤ 2 their continuum limit is a CFT with c = 1 -6(γ/π) 2 )/g, where g and γ are parameters related by g = 1γ/π, and n = 2 cos γ. We use the boundary condition shown in Fig. 1: any loop which touches the boundary (shaded) is given a weight n 1 instead of n. This is a CBC for any real n 1 > 0, with universal boundary entropy [37,38]

s n,n1 = -log   (2g) -1 4 sin rγ g sin rγ sin γ sin γ g 1 2   , (13)
where

n 1 = sin((r+1)γ) sin(rγ)
is parametrized by r [37,38]. Let |0 L be the ground state of T L , normalized such that L 0|0 L = 1 for the loop scalar product (10). The boundary free energy f

(n,n1) s (L) is obtained from |B = |
L , where L B|0 L is computed as in (10), but replacing n by n 1 because the loops touch the boundary (cf. Fig. 1).

The universal term is computed numerically from f

(n,n1) s (L) for L = 14, 16, 18, 20. The results (see Fig. 2) are in excellent agreement with (13).

The trousers trick

In the remainder of this paper we will always consider systems with boundaries, and boundary CFTs. As explained in section 1, we need to identify a lattice version of the conformal state L -2 |0 . To do this, we first note that L -2 generates a conformal transformation that maps the half-plane with an infinitesimal vertical slit onto the half-plane itself. This very basic idea is at the heart of the relation between the operator formalism of CFT and the celebrated Schramm-Loewner Evolution [41,42]. By turning this infinitesimal conformal transformation into a finite one, we will be able to build its lattice version.

We proceed as in [41,42]: consider the sequence of infinitesimal transformations z → z + dg(z) where dg(z) = dt/(2z). This is the infinitesimal version of the transformation

g t (z) = z 2 + t ( 14
)
in the sense that ∂ ∂t g t (z) = 1/(2g t (z)). The tranformation ( 14) maps conformally the half-plane minus a slit of height √ t onto the half-plane (figure 3). Taking z = g -1 t (w), where w is now a point in the upper half-plane, we have

∂g t ∂t • g -1 t (w) = 1/2 w (15)
so if F is an arbitrary function defined on the half-plane, it gets changed under the infinitesimal transformation dg by

dF = dt/2 w (∂ w F ) = - dt 2 l -2 F ( 16
)
where l -k = -w -k+1 ∂ w is a generator of the Witt algebra. If Gg : F → F • g, we see that equations ( 15) and ( 16) can be written G-

1 gt • ∂ Gg t ∂t = -(1/2)l -2
(the order comes from Gf•h = Gh • Gf ). We have thus obtained a differential equation which allows to compute Ggt . This operator maps a function F (defined on the upper half-plane) onto a function F • g t (defined in the upper half-plane minus a slit as in Fig. 3). In CFT, we want to act on the Verma modules rather than on simple functions. Thus we are looking for a linear operator G gt that maps a state |s of the theory in the upper half-plane geometry onto G gt |s , encoding the new geometry. Turning the Witt generators into the Virasoro ones, we get

G -1 gt . ∂G gt ∂t = - 1 2 L -2 (17)
and of course at t = 0 the operator is the identity G g0 = 1, so we end up with The fact that the action of L -2 encodes a small slit in the original geometry motivates the introduction of the "trousers" geometry (Fig. 3). Indeed, it is well-known that a CFT on the infinite strip of width L is equivalent to a CFT in the upper half-plane (the latter can be mapped on the former by the transformation z → w(z) = L π log z). Cutting a small slit in the half-plane is then equivalent to cutting a slit in the infinite strip, dividing it into two strips of width L/2 up to some height. This is the trousers geometry (Fig. 3). The ground state of the translation operator along the strip π L L 0 -c 24 is |0 . It is different from the ground state of the system with two legs of width L/2, which is a tensor product |Trousers ≡ |0

G gt = exp - t 2 L -2 . (18) |Trous
′ ⊗ |0 ′ if |0 ′ is the ground state of width L/2. Formula (18) shows that

|Trousers = e -1 2 L-2 |0 = |0 - 1 2 L -2 |0 + . . . (19)
We get a lattice version of the trousers geometry if we glue together two ground states of the transfer matrix T L/2 or Hamiltonian H L/2 , i.e., |Trousers L ≡ |0 L/2 ⊗ |0 L/2 (see Fig. 4 for example).

Measure of b

While we are mostly interested in polymers and percolation in this paper, there remains a certain variety of models one can consider. Indeed, one of the difficulties of the field is that geometrical problems are not defined in the usual terms of local degrees of freedom and hamiltonians, and the purported LCFT one is after might well depend on the kind of questions one wants to ask. For percolation for instance, one can decide to focus on the boundaries of clusters, or on the six-vertex model version, or on one of the supersymmetric versions [43,24,18]. We will discuss this point more in the conclusion. For now, we start with the six-vertex model (or equivalently the XXZ spin chain), largely for pedagogical reasons. We then move on to a purely geometrical set-up for the case of polymers.

4.1 Measure of b in the XXZ spin chain at q = e i π 3

We consider the U q (sl 2 ) symmetric XXZ spin chain with q = e i π 3 [39]. The Hamiltonian with open boundary conditions is

H L = L-1 i=1 σ x i σ x i+1 + σ y i σ y i+1 + q + q -1 2 σ z i σ z i+1 + q -q -1 2 (σ z 1 -σ z L ) . (20)
Although H L is not hermitian, its spectrum is real. When q = e i π 3 it has a Jordan cell structure [39] for L ≥ 4.

The low-energy spectrum of H L is described by a CFT with c = 0. H L and the Virasoro generator are related by

H L ≈ E 0 + πvF L L 0 + O 1 L 2 when L → ∞.
Here the Fermi velocity is v F = 3 √ 3 [39].

Scalar product for the XXZ chain. What is the right scalar product for the XXZ chain?

We claim that this is just the Euclidean scalar product in the spin basis, treating q as a formal parameter (i.e., without complex conjugation). For example, the q-singlet state |s 2 ≡ q -1/2 ↑↓ -q 1/2 ↓↑ has a norm squared 2 s|s 2 = q + q -1 (if q had been conjugated, one would have found

2 s|s 2 = |q| + |q| -1 = 2).
There are different ways of seeing that this is the right scalar product. One simple way is that the XXZ chain can be mapped on the dense loop model introduced in section 2, with a weight n = q + q -1 for each closed loop. This mapping consists in replacing each q-singlet by a half-loop . . . ⊗ q -1/2 ↑↓ -q 1/2 ↓↑ ⊗ . . . → . . . . . . . Thus the loop scalar product, as defined in section 2.2, is the same as the one we are considering now. In particular, this scalar product can be negative, and we have already seen that this is what makes it a good candidate for being a finite-size version of the Virasoro scalar product in the scaling limit. Another check that this scalar product is the right one is that H L becomes hermitian for this new scalar product. This is a property that should be expected, because H L is related to L 0 in the scaling limit, and L 0 is hermitian for the Virasoro scalar product.

A detailed example: L = 4. To explain our strategy, let us discuss in full detail the case of the lowest interesting even size. We consider the S z = 0 sector only (note that in general [H L , S z ] = 0). In the basis {↑↑↓↓, ↑↓↑↓, ↑↓↓↑, ↓↑↑↓, ↓↑↓↑, ↓↓↑↑}, the Hamiltonian is

H 4 =                   1 2 + i √ 3 2 0 2 -3 2 + i √ 3 2 2 0 2 -1 |3 4 → α
(H 4 -E 4,0 ) =   0 ∆ 4 1 ∆ 4  
where E 4,0 = -9 2 is the ground state energy, and

∆ 4 = 4 πvF 3 2 -E 4,0 = 4 √ 3 π . The operator L πvF (H L -E L,0
) should be viewed as a lattice version of L 0 , as follows from the boundary analog of (11).

We still have to build the state |Trousers 4 for L = 4. This is just a tensor product of two times the ground state of H 2 , which is a U q (sl 2 ) singlet: |s 2 ≡ q -1/2 ↑↓ -q 1/2 ↓↑.

|Trousers 4 ≡ |s 2 ⊗ |s 2 = 1 2 -i √ 3 2 ↑↓↑↓ -(↑↓↓↑ + ↓↑↑↓) + 1 2 + i √ 3 2 ↓↑↓↑ .
Note that it is normalized such that 4 Trousers|0 4 = 1, in order to match 0|0 = 1 when using (19).

General strategy.

We have to restrict to sizes L that are multiples of 4, because one needs to build the state |Trousers out of the ground state of H L/2 in the S z = 0 sector, so L/2 must be even.

If the spectrum of H L in the S z = 0 sector is E L,0 < E L,1 < E L,2 < . . . then E L,3 is always twice degenerate (for L ≥ 4). Let |3 L , 3
L be a basis of this equal-energy subspace, so that the lattice version of L 0 reads

L πv F (H L -E L,0 ) = ∆ L 1 0 ∆ L , (21)
where

∆ L ≡ L πvF (E L,3 -E L,0 ) → 2 as L → ∞.
Again, there is an undetermined overall normalization because L 3|3 L = 0 (with the same scalar product as above). The trousers trick is used again to define a normalization-independent lattice quantity b(L) as in the case L = 4. We normalize |0 L and |Trousers L such that

L Trousers|0 L = L 0|0 L = 1. Then define b(L) ≡ 4 L Trousers| 3 L 2 L 3| 3 L . ( 22
)
When L → ∞, we expect |3 L ∼ αL -2 |0 and 3

L ∼ αΦ 1,5 |0 . Using relation (19) we get b(L) -→ L→∞ 4 (α Trousers|Φ 1,5 |0 ) 2 α 2 0|L 2 Φ 1,5 |0 = 0|L 2 Φ 1,5 |0 = T |t = b. ( 23
)
Numerical results are given below (Tab. 1).

Measure of b in the dilute polymers model

We consider a dilute polymers model on the honeycomb lattice. It is defined by the tranfer matrix (in the picture L = 6)

T L =
where each losange is a sum over eight configurations with weight x per monomer:

1 x x x 2 x 2 x 2 x 2 x 2
The boundary triangles can appear with or without monomers, with weights x or 1. There are no closed loops in this model (i.e., the loop weight is n = 0). This model is the n → 0 limit of the O(n) model. Its critical point [35]

is x = (2 + √ 2) -1/2
, and it corresponds in the scaling limit to dilute polymers (or self-avoiding walks).

The transfer matrix T L acts on configuration states that contain half-loops and empty sites (marked with dots), such as or for L = 6. We also have to encode the fact that polymers are already in the system when one applies T L the first time, or in other words that they are "connected to the infinite past". We thus introduce strings (unpaired points interpreted as loop segments connecting to the infinite past) like for example or with two and three polymers coming from infinity respectively. The action of T L can lower the number of strings: it can turn out that two polymers coming from infinity were actually two pieces of the same polymer connected to infinity. But iterations of the transfer matrix cannot create new polymers connected to the past. This "irreversibility" is the origin of the Jordan cells of T L in this model. To illustrate this, consider the case of L = 2 sites. One has T 2 = . In the basis {| , | , | } it reads

T 2 =   1 0 x 2 x 4 x 4 0 0 0 x 4   . (24)
T 2 has a rank-two Jordan cell for the eigenvalue x 4 . The same structure (of course for different eigenvalues) is found in T L for higher L.

Scalar product. We proceed as in section 2.2 to establish which scalar product is relevant here. The scalar product has to be compatible with the action of the transfer matrix. For the basis states with only empty sites and half-loops, the scalar product is obvious. First, it has to be zero if the empty sites are not the same in the two states. Then one just uses the loop scalar product (section 2.2), which is actually zero as soon as there is a closed loop (recall that n = 0). For example | = 0 because the empty sites are not the same, and | = 0 because there is a closed loop. Actually, so far the only basis state that contributes to non-zero scalar products is | = 1. What about the strings connected to the infinite past? Since two strings can be connected by the transfer matrix, this should also be allowed by the scalar product. Thus we define the scalar product as 1 if one contracts a pair of strings against a half-loop, like

| = 1.
Detailed calculation of the lattice b(L) for L = 2. Like for the XXZ chain, we give here a concrete example of our strategy for the smallest interesting size. We start from the transfer matrix (24) which can be put in Jordan form in the basis

|0 2 , |1 2 , 1′ 2 T 2 =   1 x 4 1 x 4   ,
where

|0 2 = | + x 4 1-x 4 | |1 2 = -x 6 1-x 4 | 1′ 2 = | -x 2 1-x 4 | .
The state |0 2 is normalized so that |0 2 = 1. We also need the left eigenstates of T 2 (see section 2)

2 0| = | + x 2 1-x 4 | 2 1| = -x 6 1-x 4 | 2 1′ = | -x 4 1-x 4 | .
Note that this is not a simple transposition of the right Jordan basis, and that here the strings represent polymers connected to the infinite future (not past). There are two independent undetermined overall factors in the normalizations of the cells |1 2 , 1′ 2 and 2 1| , 2 1′ . They will be fixed later using the trousers trick.

Like in the XXZ chain, we introduce another normalization of the Jordan cell to relate our results to CFT. Let ∆ 2 ≡ -

√ 3 2 2 π log x 4 , and 1 2 ≡ -2 √ 3 π 2 x 4 1′ 2 . Then in the basis |1 2 , 1 2 T 2 = exp - 2 √ 3 π 2 ∆ 2 1 0 ∆ 2 .
The same normalization will be used for the left Jordan cell: This quantity is invariant under a rescaling of the right Jordan cell

2 1 ≡ -2
|1 2 → α |1 2 , 1 2 → α 1 2 . It is also invariant under a rescaling of the left one 2 1| → β 2 1|, 2 1 → β 2
1 , so it does not depend on the particular normalizations of the states we have chosen. With the above formulas, one can easily check that b(2) = 4π √ 3

x 4 1-x 4 = 0.680801 . . . (recall that x = (2 + √ 2) -1/2 ).
General strategy. Dilute polymers are described by a CFT with c = 0 [34]. T L is related to the (continuum) translation operator along the strip (π/L)L 0 when L → ∞ as follows: where λ L,0 > λ L,1 > λ L,2 > . . . are the eigenvalues of T L , and the factor √ 3/2 is due to the honeycomb lattice. One finds that λ L,1 is always degenerate twice. When L → ∞ the first three eigenstates of T L should be identified with the CFT states |0 , L -2 |0 and Φ 3,1 |0 of lowest conformal dimensions (0, 2 and 2). The operator Φ 3,1 inserts two strings at the infinite past [34]. One can check numerically that the lattice version of the conformal eigenvalue ∆ L ≡ -

T L ≃ λ L,0 e -2 √ 3 π L L0 + O 1 L 2 , (25)
T L T L/2 T L/2 1 L L/2 0| L/2 0| T L/2 T L T L/2 L 1 |0 L/2 |0 L/2
√ 3 2 L π log (λ L,1 /λ L,0 ) converges to 2 when L → ∞. In the two-dimensional subspace corresponding to λ L,1 , one can define |1 L , 1
L such that in this basis

T L = λ L,0 exp - 2 √ 3 π L ∆ L 1 0 ∆ L .
This defines a lattice version of the conformal states L -2 |0 ≃ α |1 L and Φ 3,1 |0 ≃ α |2 L up to some overall normalization constant α. We can do the same thing for the left action of T L , thus defining the left Jordan cell L 1| , L 1 , which should be viewed as 0| L 2 ≃ β L 1| and 0| Φ 3,1 ≃ β L 2|, where β is some unknown constant.

The trousers trick is used to get rid of α and β. The trousers states are again tensor products of two ground states of size L/2 (see Fig. 4). They are all normalized such that the component of the basis state with empty sites only (| or |) is 1. Gathering all the pieces of the puzzle, we can define

b(L) ≡ 4 L Trousers| 1 L L 1|Trousers L L 1| 1 L (26)
which is independant of all the different choices of normalization, etc. As in the XXZ chain case [START_REF] Schomerus | [END_REF]

Numerical results

We compute b(L) for dilute polymers and the XXZ spin chain at q = e iπ/3 (Tab. 1). The sizes we can access are relatively small, therefore our error bars are relatively large. However they are precise enough to show that b is indeed different for polymers and the XXZ chain, and compatible with the predictions of LCFT [9] Finally, we stress that fully comparable results are obtained for hamiltonians and for transfer matrices (for the XXZ case for instance, the alternative would involve the 6 vertex model transfer matrix).

There is no b in geometrical percolation

While hamiltonians and transfer matrices give the same results in the continuum limit, a profound difference can be observed if one switches representations. That is, instead of the XXZ or 6 vertex model, we could study directly the geometrical representation of the percolation problem, either in a loop or in a bond version. In this case we found that there is no Jordan cell at the level of conformal weight two. There are indeed two degenerate fields with energies (if one works with the geometrical hamiltonian) scaling to the conformal weight h = 2 (and they coincide with the energies of the XXZ hamiltonian), but the lattice hamiltonian (L 0 ) remains fully diagonalizable. This comes from the structure of the standard representations of the Temperley-Lieb algebra at n = 1 (see the discussion below). In the geometrical representation, both the Hamiltonian H L and the transfer matrix T L are built out of the geometrical generators of the Temperley-Lieb algebra e i 's. For size L = 4 these are

e 1 = e 2 = e 3 =
For generic n, they act on the 6-dimensional module (we use the same conventions as in previous sections)

{| , | , | , | , | , | } . (27)
In this basis they can be written

e 1 =        
n 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 n 1 0 0 0 0 0 0

        e 2 =        
0 0 0 0 0 0 1 n 0 0 0 0 0 0 0 0 0 0 0 0 1 n 1 1 0 0 0 0 0 0 0 0 0 0 0 0

        e 3 =        
n 1 0 0 1 0 0 0 0 0 0 0 0 0 n 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

       
For generic n = 0 and n = 1, the above representation (for L = 4) of the Temperley-Lieb algebra is reducible. Indeed, let us introduce

P =          1 0 -1 n 0 -1 n n (n+1)(n 2 -2) 0 1 0 -1 n 0 - 1 (n+1)(n 2 -2) 0 0 1 0 0 -n-1 n 2 -2 0 0 0 1 0 -n-2 n 2 -2 0 0 0 0 1 -n-1 n 2 -2 0 0 0 0 0 1         
which exists and is invertible as soon as n = 2 cos π k , with k an integer. In general (for larger L), one can construct P using Jones-Wenzl projectors, when they exist. Then

P -1 e 1 P =         n 1 0 0 0 0 0 0 0 0 0 1 n 0         P -1 e 2 P =         0 0 1 n 0 0 0 1 n 1 0 0 0 0         P -1 e 3 P =         n 1 0 0 n 1 0 0 0 0 0 0 0 0        
and one recognizes the three standard modules (of dimensions 2, 3 and 1) of the Temperley-Lieb algebra for L = 4. Of course, one expects that this generic construction fails whenever n is a Beraha number n = 2 cos π k with k an integer. The big surprise here is that the basis change encoded in P is not singular when n = 1. On the contrary, one gets two states that are annihilated by e 1 , e 2 and e 3 :

| -| and | -|
and this remains true for larger L. Then every operator built out of the geometrical e i 's (for example H L and T L ) will have the same structure in this representation: they have an eigenvalue which is twice degenerate, however it is still diagonalizable. This is the reason why there cannot be a Jordan cell at level two in geometrical percolation. Therefore we had to work in the spin 1/2 representation of the Temperley-Lieb algebra (namely the XXZ spin chain above), which in the case n = 1 (ie ∆ = 1/2) is not equivalent to the geometrical representation.

Deformed version of geometrical percolation: although the geometrical representation above does not give rise to Jordan cells, we can deform it slightly to get a geometrical version that is equivalent to the XXZ chain. This works as follows. For L = 4, the action of the TL generators (recall that n = 1) in the basis ( 27) is now

e 1 =        
1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0

        e 2 =        
0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 y 0 0 0 0 0 0 0 0 0 0 0 0

        e 3 =        
1 1 0 0 1 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

       
where y is now an arbitrary parameter. It is easy to check that this defines a representation of the TL algebra for any y, the usual geometrical representation corresponding to y = 1. y is a weight that is given when one contracts the second line (coming from the infinite past) with the third one. If the first one is contracted with the second one though (or the third with the fourth), the weight is 1. So y is somehow a way of keeping track of the parity of the lines one has contracted. We claim that for any y = 1 this is equivalent to the XXZ representation. One finds that, if

P y =         1 0 1 -1 0 -1 -1 1 0 0 -1 0 0 0 -1 1 0 0 0 0 y-2 y-1 0 1 0 0 0 -1 0 0 1 0 0 1 y-1 0 0 0         then P -1 y e 1 P y =         0 1 1 0 0 0 0 0 0 0 0 1 1         P -1 y e 2 P y =         0 1 1 0 0 0 0 1 1 1 0 0 0         P -1 y e 3 P y =         0 1 1 0 1 1 0 0 0 0 0 0 0         .
This clearly shows that for any y = 1 the deformed representation has the structure

(1) ⊕

where the integers are the dimensions of the irreducible representations, and the arrows indicate action of the algebra. This is exactly the structure expected for the XXZ chain on four sites at q = e iπ/3 (see the conclusion below). When y = 1, the arrow between the top (1) and the middle (1) is not there, so we have instead ( 1) ⊕

(1) ⊕

.

In this sense, the fact that geometrical percolation (as defined above, i.e., with y = 1) is actually diagonalizable is an accident.

Of course, one can use this deformed representation to compute a new parameter b as in section 4.1. The Hamiltonian of the XXZ spin chain (20) can be expressed in terms of the TL generators as

H L = L -1 2 -2 L-1 i=1
e i (see also [39]). Then, proceeding exactly as in section 4.1, one can define and compute a lattice b(L) (this involves choosing the scalar product that is compatible with the action of the algebra on the basis [START_REF] Flohr | [END_REF], so this scalar product involves the parameter y when one contracts the second and third lines against a half-loop). We find that this number b(L) does not depend on the choice of y (as soon as y = 1) and that it is exactly the same as for the XXZ chain.

Meanwhile one could ask what would happen if one studied polymers in a spin chain or vertex model representation. To answer this, it is time to move to a slightly more algebraic description.

Conclusion

Underlying both the XXZ spin chain, the 6 vertex model, and the geometrical percolation is the Temperley-Lieb algebra, generated by e i with i = 1, . . . , L -1, subject to the famous relations

with n = q + q -1 . The Jordan cell structure of the hamiltonian or the transfer matrix is determined by the particular representation of this algebra one is working with: the diagram (geometrical) representation and the vertex representation do not have to behave in the same way.

To state this in more details, we recall that generic irreducible (standard) modules of the TL algebra are labelled by a number j which is integer or half integer, and have dimensions

d j = L L/2 + j - L L/2 + j + 1 ( 29
)
with the restriction that L/2 + j is integer. In the diagram representation, these modules correspond to those where 2j strings (or "through lines") propagate. In the vertex representation, they correspond to a fixed value j of the U q (sl 2 ) spin.

For q a root of unity, and j such that [2j +1] q = 0 these standard modules remain irreducible. 3 For other values of j, these modules contain a proper submodule, but are indecomposable. Their structure is independent of j: the largest proper submodule is irreducible, and the quotient by this submodule is also. One can represent the structure of such a module by a diagram like

• ↓ • . ( 30
)
Here the top circle represents the states in the simple quotient module or "top" (or "head"), and the bottom circle the simple submodule or "foot". The arrow represents the action of the algebra; there is some element of the algebra that maps the top to the foot, but not vice versa, as well as elements that map the top into itself and the foot into itself. Equivalently, the diagram in (30) indicates that, in a basis ordered as (bottom,top), the Temperley Lieb generators take an upper triangular form.

Of course, since we are dealing with a non semi-simple situation, the particular indecomposables that appear in a physical realization-that is, a representation-can have a rather complicated structure. In the XXZ case, this structure involves further glueing of two standard modules to form a "diamond"

• ւ ց • • ց ւ • . ( 31
)
Here, each circle is a nonzero simple subquotient module, and arrows show the action of the algebra other than within the simple subquotients, with the convention that composites of arrows should also be understood as present implicitly.

To illustrate this, we borrow from [18] the structure of the XXZ spin chain for L = 6, as given in Fig. 5. The horizontal axis contains information about the U q (sl 2 ). The vertical axis encodes information about the TL algebra. Consider what happens above the point 0 of the horizontal axis. The module with top j = 2 and bottom j = 0, and the one with top j = 3 and bottom j = 2 are the standards, which have gotten glued by further action of the algebra: this is exactly the structure represented in figure 31 after a 90 o rotation. The two nodes at j = 2 will go over, in the continuum limit, to two irreducibles of the Virasoro algebra with character χ 1,5 (whose state with lowest energy corresponds to h = 2). The hamiltonian mixes them into rank-two Jordan cells.

The case L = 4 is similar, although the chain is too short to see the generic structure at j = 0 -one gets only three simple representations instead of the four ones forming a diamond, as illustrated at the end of the previous section.

In the continuum limit, it is expected that the diamond goes over a to diamond representation of the Virasoro algebra such as those in [12,14]

R 2 ւ ց R 0 R 3 ց ւ R 2 (32)
3 We define the q-analogs as [x]q = q x -q -x q-q -1 .

ÌÄ

U q (sl 2 ) The structure of the spin-1/2 chain for q = e iπ/3 and L = 6 sites, as a representation of U q (sl 2 )⊗TL L (q).

Here, R j corresponds to an irreducible of Virasoro with character χ 1,1+2j . In the U q (sl 2 ) spin chain (the situation would be similar in the SUSY case) there are more multiplicities associated with the extra symmetry. What this means is discussed in [18]. Now the point is that the diamond modules do not necessarily have to appear. In fact, in the diagram representation, they do not, at least for values up to L = 8 4 . In this case, there is no extra multiplicity due to the spin degrees of freedom, and what one gets is a collection of indecomposable standards with no action of Temperley-Lieb between them. This is represented schematically by the red dots and arrows in Fig. 5.

In the dilute polymer case, a similar discussion is possible, involving the dilute Temperley Lieb algebra and related geometrical representation, vertex model and spin chain. the decomposition of the Hilbert space for the latter is given in [18] and is qualitatively similar to Fig. 5. What happens now is that, in the geometrical transfer matrix, one already sees the diamonds, in contrast with the percolation case.

Our findings might have important consequences for the theory in [9], although one should be careful in jumping to conclusions. After all, b was introduced and its value conjectured within a discussion of the chiral sector of bulk properties, and this may not extend straightforwardly to the structure of Virasoro representations observed in the boundary case.

It is meanwhile a bit surprising that we do not found a Jordan cell for geometrical percolation. We are not sure how this affects the results in [15]. We however note that polymers and percolation are a bit different when considered from the n → 0 resp. Q → 1 limit point of view, as mentioned briefly in [START_REF] Cardy | Logarithmic correlations in quenched random magnets and polymers[END_REF]. We will discuss this in more details elsewhere.

  -Ludwig term for the critical Ising chain. The Hamiltonian H L defined in (8) acts on the 2 L -dimensional space of spin configurations |s L . The ground state |0 L = σ α s |s L is the eigenstate of H L with the largest eigenvalue. The α s are real and positive. Normalize |0 L such that L 0|0 L = s α 2 s = 1. We consider this system on a semi-infinite cylinder of circumference L. Two different CBC are obtained by taking the Ising spins on the boundary to be fixed, |B = | + • • • + , or free, |B = s |s L . The surface contribution to the free energy is

Figure 1 :Figure 2 :

 12 Figure1: Boundary conditions for dense loops: each loop touching the boundary (shaded) is given a weight n 1 . Bulk loops are given a weight n. The system is periodic in the horizontal direction (dashed lines are identified).

Figure 3 :

 3 Figure 3: The trousers geometry can be developed on the half-plane by the function g 1 (z) = √ z 2 + 1. This operation is encoded in the state |Trousers = e -1 2 L-2 |0 .

4 . 4 =

 44 Now we are ready to define a lattice quantity b(L = 4), that is invariant under a global rescaling of the Jordan cell |3 4 → α |3 4 and 3 4 → 3 4 and that will later correspond to b in the L → ∞ limit (see below): With the above expressions of the states one gets b(4) = -√ 3π -1.3603495 . . . .

√ 3 π 2 x 4 2 2 ) ≡ 4 2

 32422 1′ . Now we introduce the trousers state. It turns out that it is trivial for L = 2, because the ground state of the transfer matrix on L/2 sites, i.e. T 1 , is trivial. Then we have simply |Trousers 2 = | . It is normalized so that |Trousers 2 = 1. The same is true of course for the left trousers state 2 Trousers| = |. Now let us define b(Trousers|

Figure 4 :

 4 Figure 4: The different lattice geometries used in the computation of b(L) for polymers. The first one (left) is used represents the term L Trousers| 1L , the other one (right) the term L 1|Trousers L .

  , we expect that b(L) -→ L→∞ b. However b = 0|L 2 Φ 3,1 |0 now, and it should be different from the number b that one finds in the XXZ chain.

  the measures of b(L) and extrapolation for L → ∞. The results are compatible with b = 5/6 ≃ 0.833 for polymers and b = -5/8 = -0.625 for the XXZ chain at q = e iπ/3 .

e 2 i

 2 = n e i e i e i±1 e i = e i [e i , e j ] = 0 for |i -j| ≥ 2

Figure 5 :

 5 Figure5: The structure of the spin-1/2 chain for q = e iπ/3 and L = 6 sites, as a representation of U q (sl 2 )⊗TL L (q).

  |3 4 and 3′ 4 → α 3′4 , thus changing the scalar product 4 3| 3′ 4 . The whole point of the trousers trick will be to get rid of this undetermined factor α.To prepare the comparison with CFT, we introduce another normalization. Let 3 4 ≡

	πvF 4	3′	4 , then in the basis |0 4 , |3 4 , 3	4 we have
			4
			πv F

Table 1 :

 1 : 5/6 and -5/8. Numerical results for

	Polymers b(L)	XXZ b(L)
	Size L	b(L)	Size L	b(L)
	2	0.68080
	4	0.66431	4	-1.36035
	6	0.67032	8	-0.87027
	8	0.67893	12	-0.75399
	10	0.68753	16	-0.70564
	12	0.69551	20	-0.68012
	14	0.70273
	∞	0.79±0.08	∞	-0.61±0.02

We shall define b precisely below.

This section can be skipped on the first reading.

A general proof most likely exist; this will be addressed elsewhere.

Acknowledgments. We thank I. Affleck, A.M. Gainutdinov, V. Gurarie, A.W.W. Ludwig, V. Pasquier, J. Rasmussen and J.M. Stéphan for discussions. This work was supported by the Agence Nationale de la Recherche (grant ANR-06-BLAN-0124-03).

With the above scalar product (i.e., without complex conjugation), we see that

Note that 3′

4 is only defined up to some additional term 3′ 4 + λ |3 4 . The scalar product 4 3′ | 3′ 4 is not invariant under such a transformation, so this is certainly not a well-defined quantity. On the contrary, 4 3| 3′ 4 seems well-defined. However, since 4 3|3 4 = 0, we see that there is still some undetermination: the structure of the Jordan cell is the same if one rescales