H T Diep

Quantum Theory of Helimagnetic Thin Films

Keywords: numbers: 75.25.-j, 75.30.Ds, 75.70.-i

We study properties of a helimagnetic thin film with quantum Heisenberg spin model by using the Green's function method. Surface spin configuration is calculated by minimizing the spin interaction energy. It is shown that the angles between spins near the surface are strongly modified with respect to the bulk configuration. Taking into account this surface spin reconstruction, we calculate self-consistently the spin-wave spectrum and the layer magnetizations as functions of temperature up to the disordered phase. The spin-wave spectrum shows the existence of a surface-localized branch which causes a low surface magnetization. We show that quantum fluctuations give rise to a crossover between the surface magnetization and interior-layer magnetizations at low temperatures. We calculate the transition temperature and show that it depends strongly on the helical angle. Results are in agreement with existing experimental observations on the stability of helical structure in thin films and on the insensitivity of the transition temperature with the film thickness. We also study effects of various parameters such as surface exchange and anisotropy interactions. Monte Carlo simulations for the classical spin model are also carried out for comparison with the quantum theoretical result.

I. INTRODUCTION

Helimagnets have been discovered a long time ago by Yoshimori [1] and Villain [2]. In the simplest model, the helimagnetic ordering is non collinear due to a competition between nearest-neighbor (NN) and next-nearestneighbor (NNN) interactions: for example, a spin in a chain turns an angle θ with respect to its previous neighbor. Low-temperature properties in helimagnets such as spin-waves [3][START_REF] Rastelli | Quantum fluctuations in helimagnets[END_REF][START_REF] Diep | Low-temperature properties of quantum Heisenberg helimagnets[END_REF][START_REF] Quartu | Phase diagram of bodycentered tetragonal Helimagnets[END_REF] and heat capacity [START_REF] Stishov | Magnetic phase transition in the itinerant helimagnet MnSi: Thermodynamic and transport properties[END_REF] have been extensively investigated. Helimagnets belong to a class of frustrated vector-spin systems. In spite of their long history, the nature of the phase transition in bulk helimagnets as well as in other non collinear magnets such as stacked triangular XY and Heisenberg antiferromagnets has been elucidated only recently [START_REF] Diep | Magnetic transitions in helimagnets[END_REF][START_REF] Ngo | Stacked triangular XY antiferromagnets: End of a controversial issue on the phase transition[END_REF][START_REF] Ngo | Phase transition in Heisenberg stacked triangular antiferromagnets: End of a controversy[END_REF]. For reviews on many aspects of frustrated spin systems, the reader is referred to Ref. [START_REF]Frustrated Spin Systems[END_REF].

In this paper, we study a helimagnetic thin film with the quantum Heisenberg spin model. Surface effects in thin films have been widely studied theoretically, experimentally and numerically, during the last three decades [START_REF]Ultrathin Magnetic Structures[END_REF][START_REF] Zangwill | Physics at Surfaces[END_REF]. Nevertheless, surface effects in helimagnets have only been recently studied: surface spin structures [START_REF] Mello | Magnetic surface phase of thin helimagnetic films[END_REF], Monte Carlo (MC) simulations [START_REF] Cinti | Exotic magnetic structures in ultrathin helimagnetic holmium films[END_REF] and a few experiments [START_REF] Karhu | Helical magnetic order in MnSi thin films[END_REF][START_REF] Karhu | Chiral modulation and reorientation effects in MnSi thin films[END_REF]. Helical magnets present potential applications in spintronics with predictions of spin-dependent electron transport in these magnetic materials [18][19][20]. This has motivated the present work. We shall use the Green's function (GF) method to study a quantum spin * Electronic address: diep@u-cergy.fr model on a helimagnetic thin film of body-centered cubic (BCC) lattice. The GF method has been initiated by Zubarev [21] for collinear bulk magnets (ferromagnets and antiferromagnets) and by Diep-The-Hung et al. for collinear surface spin configurations [START_REF] Diep-The-Hung | Effect of surface spin-waves and surface anisotropy in magnetic thin films at finite temperatures[END_REF]. For non collinear magnets, the GF method has also been developed for bulk helimagnets [START_REF] Quartu | Phase diagram of bodycentered tetragonal Helimagnets[END_REF] and for frustrated films [START_REF] Ngo | Effects of frustrated surface in Heisenberg thin films[END_REF][START_REF] Ngo | Frustration effects in antiferrormagnetic face-centered cubic Heisenberg films[END_REF]. In helimagnets, the presence of a surface modifies the competing forces acting on surface spins. As a consequence, as will be shown below, the angles between neighboring spins become non-uniform, making calculations harder. This explains why there is no microscopic calculation so far for helimagnetic films.

The paper is organized as follows. In section II, the model is presented and classical ground state (GS) of the helimagnetic film is determined. In section III, the general GF method for non-uniform spin configurations is shown in details. The GF results are shown in section IV where the spin-wave spectrum, the layer magnetizations and the transition temperature are shown. Effects of surface interaction parameters and the film thickness are discussed. Concluding remarks are given in section V.

II. MODEL AND CLASSICAL GROUND STATE

Let us recall that bulk helical structures are due to the competition of various kinds of interaction [1,2,[25][26][27]. We consider hereafter the simplest model for a film: the helical ordering is along one direction, namely the c-axis perpendicular to the film surface.

We consider a thin film of BCC lattice of N z layers, with two symmetrical surfaces perpendicular to the caxis, for simplicity. The exchange Hamiltonian is given by

H e = - i,j J i,j S i • S j (1)
where J i,j is the interaction between two quantum Heisenberg spins S i and S j occupying the lattice sites i and j.

To generate a bulk helimagnetic structure, the simplest way is to take a ferromagnetic interaction between NNs, say J 1 (> 0), and an antiferromagnetic interaction between NNNs, J 2 < 0. It is obvious that if |J 2 | is smaller than a critical value |J c 2 |, the classical GS spin configuration is ferromagnetic [3][START_REF] Rastelli | Quantum fluctuations in helimagnets[END_REF][START_REF] Diep | Low-temperature properties of quantum Heisenberg helimagnets[END_REF]. Since our purpose is to investigate the helimagnetic structure near the surface and surface effects, let us consider the case of a helimagnetic structure only in the c-direction perpendicular to the film surface. In such a case, we assume a non-zero J 2 only on the c-axis. This assumption simplifies formulas but does not change the physics of the problem since including the uniform helical angles in two other directions parallel to the surface will not introduce additional surface effects. Note that the bulk case of the above quantum spin model have been studied by the Green function method [START_REF] Quartu | Phase diagram of bodycentered tetragonal Helimagnets[END_REF].

Let us recall that the helical structure in the bulk is planar: spins lie in planes perpendicular to the c-axis: the angle between two NNs in the adjacent planes is a constant and is given by cos α = -J 1 /J 2 for a BCC lattice. The helical structure exists therefore if

|J 2 | ≥ J 1 , namely |J c 2 |(bulk)= J 1 [see Fig. 1 (top)].
To calculate the classical GS surface spin configuration, we write down the expression of the energy of spins along the c-axis, starting from the surface:

E = -Z 1 J 1 cos(θ 1 -θ 2 ) -Z 1 J 1 [cos(θ 2 -θ 1 ) + cos(θ 2 -θ 3 )] + ... -J 2 cos(θ 1 -θ 3 ) -J 2 cos(θ 2 -θ 4 ) -J 2 [cos(θ 3 -θ 1 ) + cos(θ 3 -θ 5 )] + ... (2)
where Z 1 = 4 is the number of NNs in a neighboring layer, θ i denotes the angle of a spin in the i-th layer made with the Cartesian x axis of the layer. The interaction energy between two NN spins in the two adjacent layers i and j depends only on the difference α i ≡ θ i -θ i+1 .

The GS configuration corresponds to the minimum of E.

We have to solve the set of equations:

∂E ∂α i = 0, for i = 1, N z -1 (3)
Explicitly, we have

∂E ∂α 1 = 8J 1 sin α 1 + 2J 2 sin(α 1 + α 2 ) = 0 (4) ∂E ∂α 2 = 8J 1 sin α 2 + 2J 2 sin(α 1 + α 2 ) +2J 2 sin(α 2 + α 3 ) = 0 (5) ∂E ∂α 3 = 8J 1 sin α 3 + 2J 2 sin(α 2 + α 3 ) +2J 2 sin(α 3 + α 4 ) = 0 (6) ∂E ∂α 4 = ...
where we have expressed the angle between two NNNs as follows:

θ 1 -θ 3 = θ 1 -θ 2 + θ 2 -θ 3 = α 1 + α 2 etc.
In the bulk case, putting all angles α i in Eq. 5 equal to α we get cos α = -J 1 /J 2 as expected. For the spin configuration near the surface, let us consider in the first step only three parameters α 1 (between the surface and the second layer), α 2 and α 3 . We take α n = α from n = 4 inward up to n = N z /2, the other half being symmetric. Solving the first two equations, we obtain

tan α 2 = - 2J 2 (sin α 3 + sin α 1 ) 8J 1 + 2J 2 (cos α 3 + cos α 1 ) (7)
The iterative numerical procedure is as follows: i) replacing α 3 by α = arccos(-J 1 /J 2 ) and solving (4) and ( 7) to obtain α 1 and α 2 , ii) replacing these values into (6) to calculate α 3 , iii) using this value of α 3 to solve again (4) and [START_REF] Stishov | Magnetic phase transition in the itinerant helimagnet MnSi: Thermodynamic and transport properties[END_REF] to obtain new values of α 1 and α 2 , iv) repeating step ii) and iii) until the convergence is reached within a desired precision. In the second step, we use α 1 , α 2 and α 3 to calculate by iteration α 4 , assuming a bulk value for α 5 . In the third step, we use α i (i = 1 -4) to calculate α 5 and so on. The results calculated for various J 2 /J 1 are shown in Fig. 1 (bottom) for a film of N z = 8 layers. The values obtained are shown in Table I. Results of N z = 16 will be shown later.

Some remarks are in order: i) result shown is obtained by iteration with errors less than 10 -4 degrees, ii) strong angle variations are observed near the surface with oscillation for strong J 2 , iii) the angles at the film center are close to the bulk value α (last column), meaning that the surface reconstruction affects just a few atomic layers, iv) the bulk helical order is stable just a few atomic layers away from the surface even for films thicker that N z = 8 (see below). This helical stability has been experimentally observed in holmium films [28].

Note that using the numerical steepest descent method described in Ref. [START_REF] Ngo | Effects of frustrated surface in Heisenberg thin films[END_REF] we find the same result.

In the following, using the spin configuration obtained at each J 2 /J 1 we calculate the spin-wave excitation and properties of the film such as the zero-point spin contraction, the layer magnetizations and the critical tempera-

III. GREEN'S FUNCTION METHOD

Let us define the local spin coordinates as follows: the quantization axis of spin S i is on its ζ i axis which lies in the plane, the η i axis of S i is along the c-axis, and the ξ i axis forms with η i and ζ i axes a direct trihedron. Since the spin configuration is planar, all spins have the same η axis. Furthermore, all spins in a given layer are parallel. Let ξi , ηi and ζi be the unit vectors on the local (ξ i , η i , ζ i ) axes. We write

S i = S x i ξi + S y i ηi + S z i ζi (8)
S j = S x j ξj + S y j ηj + S z j ζj (9)
We have (see Fig. 2)

ξj = cos θ ij ζi + sin θ ij ξi ζj = -sin θ ij ζi + cos θ ij ξi ηj = ηi
where cos θ ij = cos(θ i -θ j ) is the angle between two spins i and j. Replacing these into Eq. ( 9) to express S j in the ( ξi , ηi , ζi ) coordinates, then calculating S i • S j , we obtain the following exchange Hamiltonian from (1):

H e = - <i,j> J i,j 1 4 (cos θ ij -1) S + i S + j + S - i S - j + 1 4 (cos θ ij + 1) S + i S - j + S - i S + j + 1 2 sin θ ij S + i + S - i S z j - 1 2 sin θ ij S z i S + j + S - j + cos θ ij S z i S z j (10)
At this stage, let us mention that according to the theorem of Mermin and Wagner [29] continuous isotropic spin models such as XY and Heisenberg spins do not have long-range ordering at finite temperatures in two dimensions. Since we are dealing with the Heisenberg model in a thin film, it is useful to add an anisotropic interaction to create a long-range ordering and a phase transition at finite temperatures. We choose the following anisotropic interaction along the in-plane local spinquantization axes z of S i and S j :

H a = - <i,j> I i,j S z i S z j cos θ ij (11)
where I i,j (> 0) is supposed to be positive, small compared to J 1 , and limited to NN on the c-axis. The full Hamiltonian is thus

H = H e + H a . ξ ξ i j j ζ i ζ j S S i Q FIG. 2:
Local coordinates in a xy-plane perpendicular to the c-axis. Q denotes θjθi.

A. General formulation for non collinear magnets

We define the following two double-time Green's functions in the real space:

G i,j (t, t ′ ) = << S + i (t); S - j (t ′ ) >> = -iθ(t -t ′ ) < S + i (t), S - j (t ′ ) > (12) F i,j (t, t ′ ) = << S - i (t); S - j (t ′ ) >> = -iθ(t -t ′ ) < S - i (t), S - j (t ′ ) > (13)
We need these two functions because the equation of motion of the first function generates functions of the second type, and vice-versa. These equations of motion are

i d dt G i,j (t, t ′ ) = S + i (t) , S - j (t ′ ) δ (t -t ′ ) - H, S + i (t) ; S - j (t ′ ) , (14)
i d dt F i,j (t, t ′ ) = S - i (t) , S - j (t ′ ) δ (t -t ′ ) - H, S - i (t) ; S - j (t ′ ) , (15)
Expanding the commutators, and using the Tyablikov decoupling scheme for higher-order functions, for example

<< S z i ′ S + i (t); S - j (t ′ ) >>≃< S z i ′ ><< S + i (t); S - j (t ′ ) >> etc.
, we obtain the following general equations for non collinear magnets:

i dG i,j (t, t ′ ) dt = 2 < S z i > δ i,j δ(t -t ′ ) - i ′ J i,i ′ [< S z i > (cos θ i,i ′ -1) × × F i ′ ,j (t, t ′ ) + < S z i > (cos θ i,i ′ + 1)G i ′ ,j (t, t ′ ) -2 < S z i ′ > cos θ i,i ′ G i,j (t, t ′ )] + 2 i ′ I i,i ′ < S z i ′ > cos θ i,i ′ G i,j (t, t ′ ) ( 16
)
i dF i,j (t, t ′ ) dt = i ′ J i,i ′ [< S z i > (cos θ i,i ′ -1) × × G i ′ ,j (t, t ′ ) + < S z i > (cos θ i,i ′ + 1)F i ′ ,j (t, t ′ ) -2 < S z i ′ > cos θ i,i ′ F i,j (t, t ′ )] -2 i ′ I i,i ′ < S z i ′ > cos θ i,i ′ F i,j (t, t ′ ) (17)

B. BCC helimagnetic films

In the case of a BCC thin film with a (001) surface, the above equations yield a closed system of coupled equations within the Tyablikov decoupling scheme. For clarity, we separate the sums on NN interactions and NNN interactions as follows:

i dG i,j (t, t ′ ) dt = 2 < S z i > δ i,j δ(t -t ′ ) - k ′ ∈N N J i,k ′ [< S z i > (cos θ i,k ′ -1) × × F k ′ ,j (t, t ′ ) + < S z i > (cos θ i,k ′ + 1)G k ′ ,j (t, t ′ ) -2 < S z k ′ > cos θ i,k ′ G i,j (t, t ′ )] + 2 k ′ ∈N N I i,k ′ < S z k ′ > cos θ i,k ′ G i,j (t, t ′ ) - i ′ ∈N N N J i,i ′ [< S z i > (cos θ i,i ′ -1) × × F i ′ ,j (t, t ′ ) + < S z i > (cos θ i,i ′ + 1)G i ′ ,j (t, t ′ ) -2 < S z i ′ > cos θ i,i ′ G i,j (t, t ′ )] ( 18
) i dF k,j (t, t ′ ) dt = i ′ ∈N N J k,i ′ [< S z k > (cos θ k,i ′ -1) × × G i ′ ,j (t, t ′ ) + < S z k > (cos θ k,i ′ + 1)F i ′ ,j (t, t ′ ) -2 < S z i ′ > cos θ k,i ′ F k,j (t, t ′ )] -2 i ′ ∈N N I k,i ′ < S z i ′ > cos θ k,i ′ F k,j (t, t ′ ) + k ′ ∈N N N J k,k ′ [< S z k > (cos θ k,k ′ -1) × × G k ′ ,j (t, t ′ ) + < S z k > (cos θ k,k ′ + 1)F k ′ ,j (t, t ′ ) -2 < S z k ′ > cos θ k,k ′ F k,j (t, t ′ )] (19)
For simplicity, except otherwise stated, all NN interactions (J k,k ′ , I k,k ′ ) are taken equal to (J 1 , I 1 ) and all NNN interactions are taken equal to J 2 in the following. Furthermore, let us define the film coordinates which are used below: the c-axis is called z-axis, planes parallel to the film surface are called xy-planes and the Cartesian components of the spin position R i are denoted by (ℓ i , m i , n i ).

We now introduce the following in-plane Fourier transforms:

G i,j (t, t ′ ) = 1 ∆ BZ dk xy 1 2π +∞ -∞ dωe -iω(t-t ′ ) ×g ni,nj (ω, k xy ) e ikxy•(Ri-Rj) , (20)
F k,j (t, t ′ ) = 1 ∆ BZ dk xy 1 2π +∞ -∞ dωe -iω(t-t ′ ) ×f n k ,nj (ω, k xy ) e ikxy•(R k -Rj ) , ( 21
)
where ω is the spin-wave frequency, k xy denotes the wavevector parallel to xy planes and R i is the position of the spin at the site i. n i , n j and n k are respectively the zcomponent indices of the layers where the sites R i , R j and R k belong to. The integral over k xy is performed in the first Brillouin zone (BZ) whose surface is ∆ in the xy reciprocal plane. For convenience, we denote n i = 1 for all sites on the surface layer, n i = 2 for all sites of the second layer and so on. Note that for a three-dimensional case, making a 3D Fourier transformation of Eqs. ( 18)-( 19) we obtain the spin-wave dispersion relation in the absence of anisotropy:

ω = ± A 2 -B 2 (22)
where

A = J 1 S z [cos θ + 1]Zγ + 2ZJ 1 S z cos θ +J 2 S z [cos(2θ) + 1]Z c cos(k z a) +2Z c J 2 S z cos(2θ) B = J 1 S z (cos θ -1)Zγ +J 2 S z [cos(2θ) -1]Z c cos(k z a)
where Z = 8 (NN number), Z c = 2 (NNN number on the c-axis), γ = cos(k x a/2) cos(k y a/2) cos(k z a/2) (a: lattice constant). We see that ω is zero when A = ±B, namely at k x = k y = k z = 0 (γ = 1) and at k z = 2θ along the helical axis. The case of ferromagnets (antiferromagnets) with NN interaction only is recovered by putting cos θ = 1 (-1) [START_REF] Diep-The-Hung | Effect of surface spin-waves and surface anisotropy in magnetic thin films at finite temperatures[END_REF]. Let us return to the film case. We make the in-plane Fourier transformation Eqs. ( 20)- (21) for Eqs. ( 18)- (19). We obtain the following matrix equation

M (ω) h = u, (23)
where M (ω) is a square matrix of dimension (2N z × 2N z ), h and u are the column matrices which are defined as follows

h =              g 1,n ′ f 1,n ′ . . . g n,n ′ f n,n ′ . . . g Nz,n ′ f Nz,n ′              , u =        2 S z 1 δ 1,n ′ 0 . . . 2 S z Nz δ Nz,n ′ 0        , (24)
where, taking = 1 hereafter,

M (ω) =             ω + A 1 0 B + 1 C + 1 D + 1 E + 1 0 0 0 0 0 0 0 ω -A 1 -C + 1 -B + 1 -E + 1 -D + 1 0 0 0 0 0 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • D - n E - n B - n C - n ω + A n 0 B + n C + n D + n E + n • • • • • • -E - n -D - n -C - n -B - n 0 ω -A n -C + n -B + n -E + n -D + n • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0 0 0 0 0 0 D - Nz E - Nz B - Nz C - Nz ω + A Nz 0 0 0 0 0 0 0 -E - Nz -D - Nz -C - Nz -B - Nz 0 ω -A Nz             (25)
where

A n = -8J 1 (1 + d) S z n+1 cos θ n,n+1 + S z n-1 cos θ n,n-1 -2J 2 S z n+2 cos θ n,n+2 + S z n-2 cos θ n,n-2
where n = 1, 2, ..., N z , d = I 1 /J 1 , and

B ± n = 4J 1 S z n (cos θ n,n±1 + 1)γ C ± n = 4J 1 S z n (cos θ n,n±1 -1)γ E ± n = J 2 S z n (cos θ n,n±2 -1) D ± n = J 2 S z n (cos θ n,n±2 + 1
) In the above expressions, θ n,n±1 the angle between a spin in the layer n and its NN spins in layers n ± 1 etc. and γ = cos kxa 2 cos kya 2

.

Solving det|M| = 0, we obtain the spin-wave spectrum ω of the present system: for each value (k x , k y ), there are 2N z eigen-values of ω corresponding to two opposite spin precessions as in antiferromagnets (the dimension of det|M| is 2N z × 2N z ). Note that the above equation depends on the values of < S z n > (n = 1, ..., N z ). Even at temperature T = 0, these z-components are not equal to 1/2 because we are dealing with an antiferromagnetic system where fluctuations at T = 0 give rise to the so-called zero-point spin contraction [START_REF] Diep | Theory of Magnetism[END_REF]. Worse, in our system with the existence of the film surfaces, the spin contractions are not spatially uniform as will be seen below. So the solution of det|M| = 0 should be found by iteration. This will be explicitly shown hereafter.

The solution for g n,n is given by

g n,n (ω) = |M| 2n-1 |M| , (26)
where |M| 2n-1 is the determinant made by replacing the 2n -1-th column of |M| by u given by Eq. [START_REF] Ngo | Frustration effects in antiferrormagnetic face-centered cubic Heisenberg films[END_REF] [note that g n,n occupies the (2n -1)-th line of the matrix h].

Writing now

|M| = i [ω -ω i (k xy )] , (27)
we see that

ω i (k xy ) , i = 1, • • • , 2N z , are poles of g n,n
. ω i (k xy ) can be obtained by solving |M| = 0. In this case, g n,n can be expressed as

g n,n (ω) = i D 2n-1 (ω i (k xy )) [ω -ω i (k xy )] , (28)
where

D 2n-1 (ω i (k xy )) is D 2n-1 (ω i (k xy )) = |M| 2n-1 (ω i (k xy )) j =i [ω j (k xy ) -ω i (k xy )]
.

Next, using the spectral theorem which relates the correlation function S - i S + j to the Green's function [21], we have

S - i S + j = lim ε→0 1 ∆ dk xy +∞ -∞ i 2π g n,n ′ (ω + iε) -g n,n ′ (ω -iε) dω e βω -1 e ikxy•(Ri-Rj ) , ( 30
)
where ǫ is an infinitesimal positive constant and β = (k B T ) -1 , k B being the Boltzmann constant.

Using the Green's function presented above, we can calculate self-consistently various physical quantities as functions of temperature T . The magnetization S z n of the n-th layer is given by

S z n = 1 2 -S - n S + n = 1 2 -lim ǫ→0 1 ∆ dk xy +∞ -∞ i 2π [g n,n (ω + iǫ) -g n,n (ω -iǫ)] dω e βω -1 (31)
Replacing Eq. ( 28) in Eq. ( 31) and making use of the following identity

1 x -iη - 1 x + iη = 2πiδ(x) (32)
we obtain

S z n = 1 2 - 1 ∆ dk x dk y 2Nz i=1 D 2n-1 (ω i ) e βωi -1 (33)
where n = 1, ..., N z . As < S z n > depends on the magnetizations of the neighboring layers via ω i (i = 1, ..., 2N z ), we should solve by iteration the equations (33) written for all layers, namely for n = 1, ..., N z , to obtain the magnetizations of layers 1, 2, 3, ..., N z at a given temperature T . Note that by symmetry, < S z 1 >=< S z Nz >, < S z 2 >=< S z Nz-1 >, < S z 3 >=< S z Nz-2 >, and so on. Thus, only N z /2 self-consistent layer magnetizations are to be calculated.

The value of the spin in the layer n at T = 0 is calculated by

S z n (T = 0) = 1 2 + 1 ∆ dk x dk y Nz i=1 D 2n-1 (ω i ) (34)
where the sum is performed over N z negative values of ω i (for positive values the Bose-Einstein factor is equal to 0 at T = 0). The transition temperature T c can be calculated in a self-consistent manner by iteration, letting all < S z n > tend to zero, namely ω i → 0. Expanding e βωi -1 → β c ω i on the right-hand side of Eq. ( 33) where β c = (k B T c ) -1 , we have by putting S z n = 0 on the left-hand side,

β c = 2 1 ∆ dk x dk y 2Nz i=1 D 2n-1 (ω i ) ω i (35)
There are N z such equations using Eq. ( 33) with n = 1, ..., N z . Since the layer magnetizations tend to zero at the transition temperature from different values, it is obvious that we have to look for a convergence of the solutions of the equations Eq. ( 35) to a single value of T c . The method to do this will be shown below.

IV. RESULTS FROM THE GREEN'S FUNCTION METHOD

Let us take J 1 = 1, namely ferromagnetic interaction between NN. We consider the helimagnetic case where the NNN interaction J 2 is negative and |J 2 | > J 1 . The non uniform GS spin configuration across the film has been determined in section II for each value of p = J 2 /J 1 . Using the values of θ n,n±1 and θ n,n±2 to calculate the matrix elements of |M|, then solving det|M| = 0 we find the eigenvalues ω i (i = 1, ..., 2N z ) for each k xy with a input set of S z n (n = 1, ..., N z ) at a given T . Using Eq. (33) for n = 1, ..., N z we calculate the output S z n (n = 1, ..., N z ). Using this output set as input, we calculate again S z n (n = 1, ..., N z ) until the input and output are identical within a desired precision P . Numerically, we use a Brillouin zone of 100 2 wave-vector values, and we use the obtained values S z n at a given T as input for a neighboring T . At low T and up to ∼ 4 5 T c , only a few iterations suffice to get P ≤ 1%. Near T c , several dozens of iteration are needed to get convergence. We show below our results.

A. Spectrum

We calculated the spin-wave spectrum as described above for each a given J 2 /J 1 . The spin-wave spectrum depends on the temperature via the temperaturedependence of layer magnetizations. Let us show in Fig. 3 the spin-wave frequency ω versus k x = k y in the case on an 8-layer film where J 2 /J 1 = -1.4 at two temperatures T = 0.1 and T = 1.02 (in units of J 1 /k B = 1). There are 8 positive and 8 negative modes corresponding two opposite spin precessions. Note that there are two degenerate acoustic surface branches lying at low energy on each side. This degeneracy comes from the two symmetrical surfaces of the film. These surface modes propagate parallel to the film surface but are damped from the surface inward. As T increases, layer magnetizations decrease (see below), reducing therefore the spin-wave energy as seen in Fig. 3 (bottom). is known that in antiferromagnets, quantum fluctuations give rise to a contraction of the spin length at zero temperature [START_REF] Diep | Theory of Magnetism[END_REF]. We will see here that a spin under a stronger antiferromagnetic interaction has a stronger zero-point spin contraction. The spins near the surface serve for such a test. In the case of the film considered above, spins in the first and in the second layers have only one antiferromagnetic NNN while interior spins have two NNN, so the contraction at a given J 2 /J 1 is expected to be stronger for interior spins. This is verified with the results shown in Fig. 4. When |J 2 |/J 1 increases, namely the antiferromagnetic interaction becomes stronger, we observe stronger contractions. Note that the contraction tends to zero when the spin configuration becomes ferromagnetic, namely J 2 tends to -1.

C. Layer magnetizations

Let us show two examples of the magnetization, layer by layer, from the film surface in Figs. 5 and6, for the case where J 2 /J 1 = -1.4 and -2 in a N z = 8 film. Let us comment on the case J 2 /J 1 = -1.4:

(i) the shown result is obtained with a convergence of 1%. For temperatures closer to the transition temperature T c , we have to lower the precision to a few percents which reduces the clarity because of their close values (not shown).

(ii) the surface magnetization, which has a large value at T = 0 as seen in Fig. 4, crosses the interior layer magnetizations at T ≃ 0.42 to become much smaller than interior magnetizations at higher temperatures. This crossover phenomenon is due to the competition between quantum fluctuations, which dominate low-T behavior, and the low-lying surface spin-wave modes which strongly diminish the surface magnetization at higher T . Note that the second-layer magnetization makes also a crossover at T ≃ 1.3. Similar crossovers have been observed in quantum antiferromagnetic films [START_REF] Diep | Quantum effects in antiferromagnetic thin films[END_REF] and quantum superlattices [START_REF] Diep | Theory of antiferromagnetic superlattices at finite temperatures[END_REF].

Similar remarks can be also made for the case J 2 /J 1 = -2.

Note that though the layer magnetizations are different at low temperatures, they will tend to zero at a unique transition temperature as seen below. The reason is that as long as an interior layer magnetization is not zero, it will act on the surface spins as an external field, preventing them to become zero. The temperature where layer magnetizations tend to zero is calculated by Eq. ( 35). Since all layer magnetizations tend to zero from different values, we have to solve self-consistently N z equations (35) to obtain the transition temperature T c . One way to do it is to use the self-consistent layer magnetizations obtained as described above at a temperature as close as possible to T c as input for Eq. (35). As long as the T is far from T c the convergence is not reached: we have four 'pseudotransition temperatures' T cs as seen in Fig. 7, one for each layer. The convergence of these T cs can be obtained by a short extrapolation from temperatures when they are rather close to each other. T c is thus obtained with a very small extrapolation error as seen in Fig. 7 for p = J 2 /J 1 = -1.4: T c ≃ 2.313 ± 0.010. The results for several p = J 2 /J 1 are shown in Fig. 8.

D. Effect of anisotropy and surface parameters

The results shown above have been calculated with an in-plane anisotropy interaction d = 0.1. Let us show now the effect of d. Stronger d will enhance all the layer magnetizations and increase T c . Figure 9 shows the surface magnetization versus T for several values of d (other magnetizations are not shown to preserve the figure clarity). The transition temperatures are 2.091 ± 0.010, 2.313±0.010, 2.752±0.010, 3.195±0.010 and 3.630±0.010 for d = 0.05, 0.1, 0.2, 0.3 and 0.4, respectively. These values versus d lie on a remarkable straight line.

Let us examine the effects of the surface anisotropy and exchange parameters d s and J s 1 . As seen above, even in the case where the surface interaction parameters are the same as those in the bulk the surface spin-wave modes exist in the spectrum. These localized modes cause a low surface magnetization observed in Figs. 5 and6. Here, we show that with a weaker NN exchange interaction between surface spins and the second-layer ones, namely J s 1 < J 1 , the surface magnetization becomes even much smaller with respect to the magnetizations of interior layers. This is shown in Fig. 10 for several values of J s FIG. 7: (Color online) Top: Transition temperature is calculated at p = J2/J1 = -1.4 for d = 0.1, Nz = 8: at each temperature, using the self-consistent values of layer magnetizations at T < Tc, Eq.( 35) is solved to obtain Tcs for each layer. The convergence is reached when Tcs tend to a single value Tc. One has Tc ≃ 2.313 ± 0.010. Red circles, black void circles, blue squares and blue void squares are Tcs obtained from Eq. ( 35) for first, second, third and fourth layers, respectively, at different temperatures. Bottom: Extrapolation by lines to obtain Tc is shown for surface parameter J s 1 /J1 = 0.7. The precision for self-consistent convergence is 1% for layer magnetizations. See text for comments. at low T due to quantum fluctuations as discussed earlier.

The transition temperature strongly decreases with J s 1 : we have T c = 2.103 ± 0.010, 1.951 ± 0.010, 1.880 ± 0.010 and 1.841±0.010 for J s 1 = 1, 0.7, 0.5 and 0.3, respectively (N z = 16, J 2 /J 1 = -2, d = d s = 0.1). Note that the value J s 1 = 0.5 is a very particular value: the GS configuration is a uniform configuration with all angles equal to 60 • , namely there is no surface spin rearrangement. This can be explained if we look at the local field acting on the surface spins: the lack of neighbors is compensated by this weak positive value of J s 1 so that their local field is equal to that of a bulk spin. There is thus no surface reconstruction. Nevertheless, as T increases, thermal effects will strongly diminish the surface magnetization as seen in Fig. 10 (middle). As for the surface anisotropy parameter d s , it affects strongly the layer magnetizations and the transition temperature: we show in Fig. 11 the surface magnetizations and the transition temperature for several values of d s .

E. Effect of the film thickness

We have performed calculations for N z = 8, 12 and 16. The results show that the effect of the thickness at these values is not significant: the difference lies within convergence errors. Note that the classical ground state of the first four layers is almost the same: for example, here are the values of cosinus of the angles of the film first half for N z = 16 which are to be compared with the values for N z = 8 given in Table I At p = J 2 /J 1 = -2, the transition temperature is 2.090 ± 0.010 for N z = 8, 2.093 ± 0.010 for N z = 12 and 2.103 ± 0.010 for N z = 16. These are the same within errors. At smaller thicknesses, the difference can be seen. However, for helimagnets in the z direction, thicknesses smaller than 8 do not allow to see fully the surface helical reconstruction which covers the first four layers: to study surface helical effects in such a situation would not make sense.

At this stage, it is interesting to note that our result is in excellent agreement with experiments: it has been experimentally observed that the transition temperature does not vary significantly in MnSi films in the thickness range of 11-40 nm [START_REF] Karhu | Helical magnetic order in MnSi thin films[END_REF]. One possible explanation is that the helical structure is very stable as seen above: the surface perturbs the bulk helical configuration only at the first four layers, so the bulk 'rigidity' dominates the transition. This has been experimentally seen in holmium films [28].

F. Classical helimagnetic films: Monte Carlo simulation

To appreciate quantum effects causing crossovers of layer magnetizations presented above at low temperatures, we show here some results of the classical counterpart model: spins are classical XY spins of amplitude S = 1. We take the XY spins rather than the Heisenberg spins for comparison with the quantum case because in the latter case we have used an in-plane Ising-like anisotropy interaction d. Monte Carlo simulations have been carried out over film samples of 100 × 100 × 16. Periodic boundary conditions are applied in the xy plane. One million of MC steps are discarded to equilibrate the system and another million of MC steps are used for averaging. The layer magnetizations versus T are shown in Fig. 12 for the case where surface interaction J s 1 = 1 (top) and 0.3 (bottom) with J 2 /J 1 = -2 and N z = 16. One sees that i) by extrapolation there is no spin contraction at T = 0 and there is no crossover of layer magnetizations at low temperatures, ii) from the intermediate temperature region up to the transition the relative values of layer magnetizations are not always the same as in the quantum case: for example at T = 1.2, one has M 1 < M 3 < M 4 < M 2 in Fig. 12 M 1 < M 2 < M 4 < M 3 in Fig. 12 (bottom) which are not the same as in the quantum case shown in Fig. 6 (top) and Fig. 10 (top). Our conclusion is that even at temperatures close to the transition, helimagnets may have slightly different behaviors according to their quantum or classical nature. Extensive MC simulations with size effects and detection of the order of the phase transition is not the scope of this present paper.

V. CONCLUSION

We have studied in this paper surface effects in a helimagnet of body-centered cubic lattice with quantum Heisenberg spins. The classical bulk ground-state spin configuration is exactly calculated and is found to be strongly modified near the film surface. The surface spin rearrangement is however limited to the first four layers in our model, regardless of the bulk angle, namely the NNN interaction strength J 2 . The spin-wave excitation is calculated using a general Green's function technique for non collinear spin configurations. The layer magnetization as a function of temperature as well as the transition temperature are shown for various interaction parameters. Among the striking features found in the present paper, let us mention i) the cross-over of layer magnetizations at low temperatures due to the competition between quantum fluctuations and thermal effects, ii) the existence of low-lying surface spin-wave modes which cause a low surface magnetization, iii) a strong effect of the surface exchange interaction (J s 1 ) which drastically modifies the surface spin configuration and gives rise to a very low surface magnetization, iv) the transition temperature varies strongly with the helical angle but it is insensitive to the film thickness in agreement with experiments performed on MnSi films [START_REF] Karhu | Helical magnetic order in MnSi thin films[END_REF] and holmium [28], v) the classical spin model counterpart gives features slightly different from those of the quantum model, both at low and high temperatures.

To conclude, let us emphasize that the general theoretical method proposed here allows us to study at a microscopic level surface spin-waves and their physical consequences at finite temperatures in systems with non collinear spin configurations such as helimagnetic films. It can be used in more complicated situations such as helimagnets with Dzyaloshinskii-Moriya interactions [START_REF] Karhu | Chiral modulation and reorientation effects in MnSi thin films[END_REF].

FIG. 3 :

 3 FIG. 3: Spectrum E = ω versus k ≡ kx = ky for J2/J1 = -1.4 at T = 0.1 (top) and T = 1.02 (bottom) for Nz = 8 and d = 0.1. The surface branches are indicated by s.

FIG. 4 :

 4 FIG. 4: (Color online) Spin lengths of the first four layers at T = 0 for several values of p = J2/J1 with d = 0.1, Nz = 8. Black circles, void circles, black squares and void squares are for first, second, third and fourth layers, respectively. See text for comments.

FIG. 5 :

 5 FIG. 5: (Color online) Layer magnetizations as functions of T for J2/J1 = -1.4 with d = 0.1, Nz = 8 (top). Zoom of the region at low T to show crossover (bottom). Black circles, blue void squares, magenta squares and red void circles are for first, second, third and fourth layers, respectively. See text.

4 FIG. 6 :

 46 FIG. 6: (Color online) Layer magnetizations as functions of T for J2/J1 = -2 with d = 0.1, Nz = 8 (top). Zoom of the region at low T to show crossover (bottom). Black circles, blue void squares, magenta squares and red void circles are for first, second, third and fourth layers, respectively. See text.

FIG. 8 :FIG. 9 :

 89 FIG. 8: (Color online) Transition temperature versus p = J2/J1 for an 8-layer film with d = 0.1. See text for comments.

FIG. 10 :

 10 FIG.10:(Color online) Layer magnetizations as functions of T for the surface interaction J s 1 = 0.3 (top), 0.5 (middle) and 0.7 (bottom) with J2/J1 = -2, d = 0.1 and Nz = 16. Black circles, blue void squares, magenta squares and red void circles are for first, second, third and fourth layers, respectively.

FIG. 11 :

 11 FIG. 11: Top: Surface magnetization versus T for ds = 0.01 (circles), 0.1 (void circles), 0.2 (squares) and 0.3 (void squares), with J s 1 = 1, J2/J1 = -2 and Nz = 16. Bottom: Transition temperature versus ds for J s 1 =0.7, 0.5 and 0.3 (curves from up to down), with J2/J1 = -2, d = 0.1 and Nz = 16.

FIG. 12 :

 12 FIG.12:(Color online) Monte Carlo results: Layer magnetizations as functions of T for the surface interaction J s 1 = 1 (top) and 0.3 (bottom) with J2/J1 = -2 and Nz = 16. Black circles, blue void squares, cyan squares and red void circles are for first, second, third and fourth layers, respectively.

  , for p = J 2 /J 1 = -2

	(in parentheses are angles in degree):
	0.86737967 (29.844446),	0.41125694 (65.716179),
	0.52374715 (58.416061),	0.49363765 (60.420044),
	0.50170541 (59.887100),	0.49954081 (60.030373),
	0.50013113 (59.991325), 0.49993441 (60.004330).
	From the 4th layer, the angle is almost equal to the bulk
	value (60 • ).


Xixi Huang

Sandrine Bec

Jean Colombani

Influence of fine particles on the stability of a humid granular pile

Keywords: numbers: 45.70.Ht, 45.70.Mg, 62.20.Qp

come

Influence of fine particles on the stability of a humid granular pile

I. INTRODUCTION

The stability of humid granular media has been studied extensively during the last years and the influence of the liquid content [1][2][3], liquid surface tension [4], grain roughness [5,6], or pile age [7] have found explanations. In these studies, the working granular pile is a monodisperse set of spherical grains. Besides, the study of bidisperse granular media is devoted to the investigation of segregation and mixing [8]. But monodisperse granular media are hardly ever really monodisperse and contain dust due to the wear of the grains [9]. This is always true not only in sandcastle building on the beach, but also in industrial use of granular media (material grinding, cereal storage, pharmaceutical manufacturing, . . . ).

Adding a small controlled amount of micron-size particles to a granular pile, we show by rotating drum experiments that these fine grains have a dramatic effect on the stability of the pile. Furthermore, this influence is ambivalent: the fine particles fluidify the heap at low content, and solidify it at higher content. We propose a microscopic mechanism and a theoretical interpretation of this non-monotonic behavior, based on capillarity, jamming, lubrication, and friction considerations.

II. MATERIALS AND METHOD

The experimental setup consists of a stainless steel rotating drum of 10 cm radius and width, with glass windows on the sides. The windows have a hole in the center to allow air exchange with the outside atmosphere. The drum is partially filled with soda-lime glass beads of (250±100) µm diameter mixed with small beads of (5±4) * Electronic address: jean.colombani@univ-lyon1.fr µm diameter. The pile is illuminated from the front, observed with a camera, and its surface is detected by an automatic procedure. Its angle θ with the horizontal is measured with a 0.5 • accuracy. The drum is placed inside a closed vessel where a constant relative humidity H is maintained by the presence of a wide beaker of saturated aqueous inorganic salt solution. Each salt fixes the relative humidity in the atmosphere at a given value (6% with LiBr, 23% with KCH 3 CO 2 , 43% with K 2 CO 3 , 58% with NaBr). The values of H are measured with an accuracy of 1.5% by an hygrometer inside the drum. The temperature in the box is controlled by a heating resistor regulated by a temperature controller, and is always fixed at 29.0 • C. Homogeneity of temperature and hygrometry inside the vessel is provided by a fan. To enable the vapor to penetrate the pile, the drum is filled with the grains, placed in the closed vessel, and rotated at 1 rpm during 12 h. Just before the experiment, the drum is rotated at 5 rpm during a few seconds to provide an initial reference state.

As the waiting time and waiting angle before the pile destabilization have been found to be crucial parameters, due to capillary bridges nucleation, a fixed waiting time t w of 3 s, 1 • below the maximum stability angle, has been respected for each experiment [7,9]. To reduce the statistical uncertainty, every angle has been measured one hundred times for each of the 12 fine grains contents, 4 relative humidities and 23 rotation speeds. To check the absence of noticeable segregation, we have sampled the top and bottom of the heap at the end of some experiments, and verified by optical microscopy that they have similar content of fine particles. To test the presence of electrical charges in the heap, we have measured the voltage of the granular heap surface at a 20 mm distance with a fieldmeter before and just after a rotating drum experiment. We had not been able to measure an evolution, the voltage keeping a ∼50 V value, close to the sensitivity of the device, which can be considered as negligible, usual voltages in charged isolating materials being of the order of kV.

III. RESULTS

We have performed rotating drum experiments with a mixture of 250 and 5 µm glass beads. By varying the rotation rate ω, the percentage of fine particles φ (volume ratio of fine and coarse spheres) and the relative humidity H, we have built a stability diagram of the heap (see figure 1 and videos in Supplemental Material [25]).

At small concentrations, the rotating heap shows intermittent avalanches, i.e., the top surface remains flat until a maximum stability angle θ m is reached, where the grains flow to achieve a new stable flat surface state at the repose angle θ r (video #1). When the rotation speed is increased, the rotation time becomes too short compared to the avalanche duration, and a transition from the avalanche regime to a continuous flow regime is observed, where the grains at the top of the heap flow without rest, keeping a constant flat surface, characterized by the socalled dynamic repose angle θ d (video #2). An intermittent regime is observed during the avalanche/continuous flow transition (figure 2) [10]. When the fine particles fraction increases, the transition occurs for lower rotation rates (figure 1) and the avalanche angle decreases (figure 3), proof of a fluidization of the pile, or at least of its surface.

Up to that point, the granular assembly exhibits expected motion patterns [11]. But at larger fine particles concentration, the pile shows a stabilization of its surface. It remains now cohesive and, unlike standard rotating drum experiments, the motion of the drum is accommodated by a stick-slip behavior at the drum wall. When the drum rotates, the pile remains stuck to the wall until a maximum stability angle is attained (as for avalanching), where it detaches and slips downward to the repose angle, without any observable inner motion, except in a layer of grains in contact with the drum wall (video #3). When the rotation rate is enhanced, stickslip stops and a thin layer of grains at the drum wall rolls continuously, thereby leading to an overall sliding motion with a constant slope of the heap surface (video #4). In configurations where granular matter exhibits stick-slip at a wall contact, e.g. shearing in a Couette geometry, such a transition between stick-slip and continuous sliding is also observed and has been interpreted as a dilatancy transition [12].

IV. DISCUSSION

So, for low rotation rates, when the fine particles content increases (i) the heap flows at the surface during more and more fluid avalanches, until (ii) this flow is stopped, and the heap slips at the bottom behaving as a quasi-solid body. This non-monotonous nature of the evolution of the maximum stability angle with fine particles thus enables to bridge a gap between independent observations. Small particles are known to increase the flowability of powders [13] or the run-out length of landslides [14], in weakening the adhesive contact between the coarse ones. Besides, the introduction in a particulate assembly of particles of smaller diameter is known to promote jamming [15], although little is known quantitatively about the influence of polydispersity on jamming [16,17]. This contradictory nature of the influence of a small amount of fine grains on a granular assembly had not been reported yet. Indeed it had not been studied in one well-defined system, so a coherent view of the influence of dust was not possible. We mention that a change of behavior is expected when the fine particles begin to dominate the mixture, from φ ∼ 40% [18,19], but was unexpected for values of φ so low.

The locus of destabilization of a granular heap is known to be highly dependent on the cohesion of the particulate assembly, the more cohesive the assembly, the deeper the breakage [20]. When the fine particles volume fraction reaches the critical value φ c 0.2%, the locus of the breakage moves from the top to the bottom of the heap, proof of a change of cohesion. The value of φ c is therefore related to the ability of the fine particles to increase the cohesion of the grain assembly. No prediction of packing geometry in available for our experiments [15]. Indeed, for such a large bead radius ratio and low fine beads percentage, small particle usually percolate. So, in order to check the distribution of fine particles inside the network of large beads, we have carried out scanning electron microscopy observations of the mixture just after a rotating drum experiment. Figure 4 clearly indicates that, even if some fine grains are stuck at the surface of big beads, most of them populate the interstitial space between large particles. The most probable source of cohesion is therefore jamming, promoted by the adhesion of small beads on large ones and between small ones, which prevents them from percolating, and finally blocks the inner motions inside the heap The influence of this jamming on the mechanical behavior of the heap may be diverse. Whereas strong in the rotating drum, we have found it much less pronounced in a parallel plate rheometer (figure 5). In other words, stuck fine particles build kinds of solid bridges, that play the same role as liquid bridges, either fluidifying, or solidifying the heap, as observed by Chou & Hsiau [21].

For low fine grains content, the avalanche regime is characterized by three features: (i) the destabilization always occurs at the surface of the heap, (ii) the maximum stability angle is influenced by the fine particles content (so the cohesion) but (iii) not by the relative humidity (figure 3). This behavior is typical of the so-called granular regime [2]. In this regime, as the real contact surface between beads (so the adherence force between them) increases proportionally with the depth inside the heap, due to elastoplatic strains, the most unstable plane remains at the surface [20]. And for our relative humidity values, the equilibrium radius of curvature of the water menisci condensing between beads is smaller than their typical roughness (∼ 10 nm). So the capillary water is confined to the tip of asperities in contact between beads and the adhesion force between them is not influenced by humidity, as theoretically described [22] and experimentally observed [23] in monodisperse media. In this regime, the maximum stability angle obeys the simple law: tan(θ m ) = tan θ 0 (1 + C), with θ 0 the angle in absence of cohesion and C a function of the adhesion stress inside the heap [20]. The theoretical value of θ 0 , computed from pure geometrical considerations, is 23.4 • in a monodisperse system [24]. We see in figure 3 that the minimum experimental value of θ m , measured at φ c , corresponds to this value. We can deduce from this equality that, before jamming the heap at high content, the fine spheres begin by playing the role of ball bearings. By populating the intergrain space (figure 4), they contribute to separate the beads and destroy their cohesion, inducing granular lubrication, until recovering the behavior of a cohesion-free heap, when θ m reaches θ 0 . By fitting the experimental results with the above law, we can conjecture that the loss of cohesion, so the granular lubrication efficiency, due to the fine particles, follows a scaling law of the form C ∼ (φ-φ c ) n with n 2.3 (figure 3).

For high fine particles contents, the granular assembly now behaves as a solid body and its sliding at the drum wall is driven by the friction between the glass spheres and the steel cylinder. The maximum stability angle θ m can thus be deduced from a standard Coulomb analysis. At first order, we perform this analysis at the location of the drum wall with the tangent parallel with the pile surface. Here the humid granular heap loses its stability when the ratio of the shear stress P sin θ m (P the weight of the above grains) to the normal stress P cos θ m + c (c the pile-wall adhesion force per unit surface) exceeds the friction coefficient µ = tan θ s [5]. As P is ρgh, with ρ the density of the heap, g the acceleration of gravity and h the height of the heap above the considered point, the stability criterion can therefore be written: tan θ m = tan θ s (1 + c/(ρgh cos θ m )). The increase of adhesion with φ results from the capillary bridges nucleating between the fine grains and the wall, bonding them to the drum. The enhancement of the density of fine particles increases their number at the drum surface, and thereby increases the total number of capillary bridges and the pile-drum adhesion. If we consider that at the minimum value of the sliding angle, for φ m 0.3%, we are in the adhesionfree case, then at this point c = 0 and θ m = θ s . The adhesion stress c is linked to the number of fine particles per unit surface, so it can be written c = A(φ -φ m ) 2/3 , with A a factor linked to the beads size and roughness and to the relative humidity H. The dependence of A on H is A = B/ ln(1/H), B being now a purely geometric factor [7]. Therefore the stability of the pile writes:

tan θ m = tan θ s 1 + B ρgh cos θ s ln(1/H) (φ -φ m ) 2/3 (1)
Here, following Halsey and Levine, we have considered that cos θ m cos θ s to get a more tractable expression [5]. The fitting in figure 3 of all the experimental sliding angles with this law, using B as only fitting parameter, enables to catch satisfactorily the influence of humidity and the curved shape of the angle evolution. Thus considering that the enhancement of friction between the particulate assembly and the wall originates from the nucleation of capillary bridges gives the correct scaling for the influence of the fine particles content and humidity, thereby validating the assumption.

V. CONCLUSION

We have shown here that the presence of dust in a monodisperse granular assembly modifies its cohesion for a content as low as a fraction of percent. Furthermore this influence is ambivalent, the fine grains firstly lubricating the pile, and then jamming it. This modification of the inner dynamics of the pile has a huge impact on its mechanical stability. In our rotating drum configuration, it induces a change of destabilization mode, from avalanching to stick-slip. It may induce different changes in other configurations. This study helps to reconcile independent and contradictory knowledges about the fine particles influence. The implications of this decisive role of dust are direct. In industrial processes, the fine grains produced at the beginning of mineral grinding for instance appear here to be likely to lower the process efficiency. Beside, this system, despite simple, has shown to be an original way of investigating the nature of arrested states (what is the exact mechanism of jamming in a strongly bidisperse mixture?), the interaction inside polydisperse aggregates (what drives the granular lubrication efficiency?), or the tribology of humid granular matter (how does a bidisperse pile interact with the container walls?).

FIG. 1 :

 1 FIG. 1: Stability diagram of the granular heap. All regimes are shown as a function of the rotation velocity of the drum and of the percentage of fine particles for a relative humidity H = 23%. The well-defined regimes are separated by transition domains showing intermittent regimes. Each of the 276 dots corresponds to an average of 100 measurements. The increase of humidity induces exclusively a slight shift of the avalanche/continuous flow and stick-slip/continuous sliding transitions toward the high rotation rate.

Acknowledgments

We thank Elisabeth Charlaix, Lydéric Bocquet and Laurent Joly for enlightening discussions, Catherine Barentin and Loren Jorgensen for the rheometric measurements, François Gay and Matthieu Guibert for experimental help, and the Carnot Institute Ingénierie@Lyon for financial support.


J Colmars
email: julien.colmars@imelavi.fr

Frédéric Dubois
email: frederic.dubois@unilim.fr

Joseph Gril
email: joseph.gril@univ-montp2.fr

One-dimensional discrete formulation of Hygrolock model for wood hygromechanics

Keywords: hygrolock, wood, mechanosorption, moisture activated creep

a new 1D discrete formulation of hygrolock models is proposed for modeling wood time-dependent behavior. This discrete formulation is compared to an integral formulation presented in a previously published paper. Simulations with various humidity cycles are performed under constant stress (creep) or constant strain (relaxation) using first a hygrolock spring and then a generalized Kelvin-Voigt model with a time-spectral distribution of hygrolocks springs. The discrete formulation presented here, based on original idea of mixed series/parallel rheological model, is shown to be very practical for implementation in scientific numerical codes. In order to compensate a lack of complete data-set on wood material, a generic time and moisture-dependent material is proposed to compare the various models: this dataset could be re-used in other studies.

Finally the relevance of hygrolock models for wood time-dependent behavior is discussed: it appears that hygrolock models suits well to the upper part of wood hygroscopic domain, whereas further hypothesis should be tested for dry domain. A judicious choice in moisture activation of material parameters, combined with the blocking effect of hygrolock, agree with recent experimental results on simplifying the description of so-called ++ effect.

Introduction Hygromechanical couplings in wood

Wood is a hygroscopic material that exhibits significant mechanical couplings under combined load and humidity changes. Wood relative creep (ratio of timedependent to elastic strain) increases with moisture content, but during moisture changes under a given load within the elastic domain, wood creeps even more than at the highest moisture condition ever reached during loading. The terms "mechanosorptive effect" or "mechanosorption" introduced by [START_REF] Grossman | The use of Leicester's "rheological model for mechano-sorptive deflection of beams[END_REF] are classically used in literature to describe this additional creep (or relaxation) due to moisture variations under stress (resp. imposed strain). This was evidenced experimentally since the 40ies (see comments by [START_REF] Schniewind | Über den Einfluss von Feuchtigkeitsänderungen auf das Kriechen von Buchenholz quer zur Faser unter Berûcksichtigung von Temperatur und Temperaturänderungen[END_REF][START_REF] Gril | Une modélisation du comportement hygro-rhéologique du bois à partir de sa microstructure[END_REF]), and later by Armstrong and Kingston (1960). These couplings are of great importance for timber structure prediction, but researchers are still working on finding simple and accurate mathematical solutions to describe this complex behavior.

A large number of rheological models used for modeling mechano-sorptive effect have been developed in the past. The first models were based on a de-coupling between viscoelastic behavior and additional effects affected by the moisture content h and its variation rate t h d d . Moisture content of wood is given by the ratio of water mass within the sample (total mass minus mass of dried sample) to the mass of the dried sample. [START_REF] Leicester | A rheological model for mechano-sorptive deflections of beams[END_REF] proposed a model which takes a dashpot form whose viscosity  is a moisture content function as described by the following stress-strain expression (1):

      t h h d d ,  (1)
At constant stress, these models produce an increment of strain under both drying and moistening (because of term t h d d ); conversely, at constant strain the model produces relaxation under any moisture change. This model also used by [START_REF] Grossman | The use of Leicester's "rheological model for mechano-sorptive deflection of beams[END_REF] allowed describing longitudinal creep and relaxation tests on hoop pine timber. In the same way, [START_REF] Rybarczik | A Theoretical Description of the Swelling Pressure of Wood[END_REF] described swelling pressure of wood across the grain (i.e. swelling-induced stress under constant strain). Later [START_REF] Ranta-Maunus | The viscoelasticity of wood at varying moisture content[END_REF], using various dataset on softwoods, proposed a generic approach of mechanosorption in wood as:

h a d d    (2)
where a is a material parameters taking different values  a and  a according to the sign of h d . He also stated that specific value   a is needed to describe the socalled "++ effect" when moisture reach the highest level ever under load.

4 Similar approach was conducted on concrete [START_REF] Bazant | Microprestress-Solidification Theory for Concrete Creep. Part I : aging and drying effects[END_REF]) to describe additional creep caused by the decreasing of humidity in the concrete pores due to hydratation reaction, see Pickett (1941) , which results mathematically in infinite strain for continuous moisture variation between given upper and lower bounds. For very small but continuous variations of h , this mechano-sorptive behavior does not converge to creep at constant moisture content.

These models have been widely used in wood studies in addition to classical elastic and viscoelastic models, which mean that modeling mechanosorption usually requires more parameters than modeling wood mechanics at constant humidity. In contrast to this, one aim of the hygrolock models showed in this paper will be to rely only on constant moisture behavior to deduce complex behavior at variable humidity; similar approach has been attempted theoretically by [START_REF] Nakano | Analysis of creep of wood during water adsorption based on the excitation response theory[END_REF].

Leicester phenomenological model was used both in series for producing nonrecoverable mechano-sorptive strain, or associated in parallel to a linear spring to take into account a linear creep limit (see [START_REF] Hunt | Creep Trajectories for Beech during Moisture Changes under Load[END_REF], Hanhijärvi and Hunt (1998), [START_REF] Montero | Influence of hygromechanical history on the longitudinal mechanosorptive creep of wood[END_REF]). For a larger review on phenomenological models applied to wood material, we recommend following references: [START_REF] Chassagne | Three Dimensional Creep Model for Wood Under Variable Humidity-Numerical Analyses at Different Material Scales[END_REF], [START_REF] Randriambololona | Modélisation du comportement différé du bois en environnement variable[END_REF][START_REF] Bou Saïd | Contribution à la modélisation des effets différés du bois et du béton sous conditions climatiques variables[END_REF], Moutee (2006) and [START_REF] Montero | Influence of hygromechanical history on the longitudinal mechanosorptive creep of wood[END_REF] for recent advances in experimental characterization of mechanosorption.

Hygrolock models

Concept of hygrolock was first proposed by [START_REF] Gril | Une modélisation du comportement hygro-rhéologique du bois à partir de sa microstructure[END_REF] to model mechanosorption in wood. The idea was to create additional deformation by blocking strain (or stress) during drying by using various modeling levels corresponding to moisture content levels. [START_REF] Dubois | Creep in Wood Under Variable Climate Conditions: Numerical Modeling and Experimental Validation[END_REF] proposed a first numerical formulation of hygrolock spring by using two different behaviors during drying and wetting: during wetting (softening of wood) a simple Hooke law is used whereas during drying (hardening of wood) Bazant law (tangent Hooke law) is preferred; this produces blocking effect during drying as proposed by [START_REF] Gril | Une modélisation du comportement hygro-rhéologique du bois à partir de sa microstructure[END_REF], but [START_REF] Randriambololona | Modélisation du comportement différé du bois en environnement variable[END_REF] showed the difficulties involved when synchronizing these two constitutive laws for moisture content cycles.

Later [START_REF] Husson | Elastic response in wood under moisture content variations: analytic development[END_REF] developed a second version of hygrolock models involving the principle of internal (or fixed) stresses proposed by [START_REF] Liu | Thermomechanics of shape memory polymers: Uniaxial experiments and constitutive modeling[END_REF]. In this case total stress on the material is the sum of mechanical stress (that induces deformation) and fixed stresses (due to rigidity variation during drying) which are stored temporary within the material and show no effect on deformation. The calculations by [START_REF] Husson | Loi de comportement viscoélastique avec effet mémoire : application à la mécanosorption dans le bois[END_REF] for the hygrolock spring led to the following integral form:

         t e t t k t 0 min d . , 1     (3) in which   min , kt 
designates the rigidity minimum between times  and t .

Introducing this mechano-sorptive spring in a generalized Kelvin Voigt (GKV) representation, [START_REF] Husson | A finite element model for shape memory behavior[END_REF] defined a double integral of time and moisture which leads to a complete time-dependent formulation incorporating strain blockings during the drying period as well as memory effects during wetting phases after unloading. Each Kelvin-Voigt body admits the following behavior:

                      d . . d . exp . 1 0 d , min                  t t k t t (4)
By using incremental formulation deduced from equation (5), [START_REF] Dubois | Modeling of the viscoelastic mechano-sorptive behavior in wood[END_REF] gave the formulation of a discrete time-dependant behavior including hygrolock effect for implementation in finite element code.

In parallel to this work, we developed another model by using a different physical approach based on moisture levels initially proposed by [START_REF] Gril | Une modélisation du comportement hygro-rhéologique du bois à partir de sa microstructure[END_REF]: unless the basic assumptions are quite different, it will be shown that these different models behave qualitatively and quantitatively the same way under various load configurations.

With the specific aim of comparing various models on the same loading history, we also propose in this paper a method to calculate a complete set of parameters from two creep curves at low and high moisture content, in order to model creep of wood at varying humidity. Finally, we discuss the relevance of such models to describe moisture and time-dependent behavior of wood.

General principles of hygrolock spring

This paper focuses on a one-dimensional approach using a discrete scheme in time and moisture levels. Hygrolock basis is an elastic body (spring) that exhibits a blocking behavior during moisture changes under the highest moisture value ever under load. The schematic representation of the hygrolock unit proposed by [START_REF] Husson | Loi de comportement viscoélastique avec effet mémoire : application à la mécanosorption dans le bois[END_REF] is shown in Figure 1.

Figure 1: symbol for hygrolock spring as proposed by [START_REF] Husson | Loi de comportement viscoélastique avec effet mémoire : application à la mécanosorption dans le bois[END_REF] Hygrolock models have various objectives in wood material behavior, which will be detailed next. It should for example:

 simplify the ++ effect of mechanosorption, defined by Ranta-

Maunus (1975),

 avoid non-thermodynamic behavior of hygro-elastic element,

 be an alternative (or additional solution) to Leicester model,  generate complex behavior at varying moisture content using only parameters determined at constant humidity.

The following sections will give more details on each point.

Hygrolock for ++ effect

The concept of hygrolock proposed by [START_REF] Gril | Une modélisation du comportement hygro-rhéologique du bois à partir de sa microstructure[END_REF] reproduces the so-called "++ effect" of mechanosorption defined by [START_REF] Ranta-Maunus | The viscoelasticity of wood at varying moisture content[END_REF]. This effect is observable during a creep test under moisture content changes: it consists in a marked increase of the strain to moisture rate ( h d d

) when moisture content reaches the highest level ever reached under load. Expression ( 2) is used to describe creep of various softwoods by using 3 parameters  a ,  a and   a . The ++ effect is illustrated on Figure 2a, with the schematic representation proposed by [START_REF] Hunt | Creep Trajectories for Beech during Moisture Changes under Load[END_REF], where total strain (of creep test) minus free swelling strain (of matched samples) is plotted against moisture content.

The hygro-lock spring as proposed by [START_REF] Gril | Une modélisation du comportement hygro-rhéologique du bois à partir de sa microstructure[END_REF] blocks strain during drying (-) and moistening (+) phases in order to underline the (++) phases (Figure 2b). and moistening phases over the highest moisture level ever (++) under load.

+ - h ε-ε h ε viscous ε elastic couplings ε-ε h h + ++ ++ - + ++ ++ - + - (a) (b)

Hygrolock as a thermodynamic spring

Another objective of hygrolock spring has been already discussed in details by [START_REF] Husson | Loi de comportement viscoélastique avec effet mémoire : application à la mécanosorption dans le bois[END_REF] and Dubois (2004): it consists in proposing a hygro-elastic body that fulfills the second thermodynamic law during both drying and moistening.

Actually a problem is raised by the dependency of elastic rigidity of wood to moisture content. Because rigidity of wood decreases with increasing moisture content [START_REF] Guitard | Mécanique du matériau bois et composites[END_REF], Kollmann and Côté (1968)), the minimal hygro-elastic body should recover partly its deformation under drying and constant load: this consequence does not fulfill the second law of thermodynamics (see [START_REF] Chassagne | Three Dimensional Creep Model for Wood Under Variable Humidity-Numerical Analyses at Different Material Scales[END_REF], [START_REF] Randriambololona | Modélisation du comportement différé du bois en environnement variable[END_REF], [START_REF] Husson | Loi de comportement viscoélastique avec effet mémoire : application à la mécanosorption dans le bois[END_REF], Dubois (2004)) and is never observed experimentally.

The classical solution for modeling strain increase observed in such conditions consists in adding other models in series (Leicester model for example) to the basic linear spring. Another solution is to enhance the formulation of the hygroelastic elementary body. Hygrolock spring as shown next is thus an enhanced hygromechanical spring, i.e. a self-sufficient model for hygro-elastic behavior. In this paper we show that using the same Hygrolock spring into a Kelvin-Voigt body leads to a satisfactory time-dependent formulation.

Incremental formulation of hygrolock spring

In this paper we model a hygrolock spring as a combination of n Hooke springs. ) is thus divided into n boxes (see Figure 3) ordered by increasing h .

Box k is thus defined as:

    n k n h k n h k h h h n n k k , ... , 1 , , 1 , 1            (5)
Characteristic moisture content of each box is here chosen to be the mean value

  2 1 k k h h  
of its upper and lower bounds. Spring k is characterized by a constant compliance k J and exhibits classical linear behavior as defined by (eq.6). As assumed for all wood species and all orthotropic directions, compliance is an increasing function of h (

k k J J  1 ).
Each box means a range of moisture content, or a level of sorption energy. The repartition of the corresponding sorption sites in the molecular structure is complex an does not fit with either a parallel or a series disposition. Therefore, a mixture of both has been chosen. Total strain E of hygrolock spring is proposed to be the sum of all spring strains (eq.7). The original point is that total stress S will also be the sum of all stresses in each level (eq.8). This choice results in a mixing between series and parallel rheological assembly.

n k J k k k    1   (6)    n k k E 1  (7)    n k k S 1  (8)
As an alternative choice, we could have defined the stress as the mean box stress or the strain as the mean box strain. In any case, the concept of boxes does not require any physical mechanism to be applicable. Moisture content at current time will set the so-called "active box": determination and role of the active box are detailed next. This type of hygrolock model is thus based on the discretization of both time and hygroscopic domain; the latter will impose sometimes the choice of time steps depending on moisture history, which is explained in next section.

ε k = J k σ k J k+1 > J k E,S Active Box E,S

Moisture event in the hygro-lock

In the following, we use a discrete scheme in time, with time steps

  1 , ,    i i i t t t t
. All superscripts i refer to variables at time i t . Active box at i t is noted i k and is given by moisture content i h . At 1

 i t active box is 1  i k .
Thus during time step i t  , three different cases might occur: moistening ( The hygrolock effect presented in this paper is based on stress transfer to subsequent active box when moistening, and stress blocking in original box when drying. Strain variations in hygrolock are directly deduced from stress transfers using expression [START_REF] Leicester | A rheological model for mechano-sorptive deflections of beams[END_REF](8).

i i k k  1 ), drying ( i i k k  1 ) or stabilization of moisture content (statu quo, i i k k  1 ). 0. n a b c Time k 0.3 1 h(t) h k(i+1) h i+1 t i t i+1 h k h k-1 h i h box h k h k-1
Whatever the event, only two boxes are put under varying stresses: i k and

1  i k ,
whereas all other boxes remain in their previous stress-strain state We further assume that a moisture change does not skip any box, i.e. that 1

1    i i k k
(this is guaranteed in numerical code by initial interpolation of input moisture data, see case (c) in Figure 4). Table 1 sums up the three different cases: we define two binary variables,  and  , which will be used to write hygrolock spring equations whatever the event.

Table 1: description of active boxes for each moisture event during time step.

Case

  Active box   i t Active box at   1  i t Moistening 1 1 i k 1 1    i i k k Statu quo 0 0 i k i i k k  1 Drying 0 1 i k 1 1    i i k k
The total variation of strain and stress in the hygrolock spring during i t  is given for any case (moistening, drying or statu quo) by expression (9-10). Note that, in order to improve readability, i k will be written   i k when used as subscript.

    1       i k i k S   (9) ) 1 ( ) 1 ( ) ( ) (         i k i k i k i k J J E   (10)
The general behavior of hygrolock spring during i t  will be written using expression (11), where P is a hygro-strain and Q a compliance.

S Q P E     (11)
In next sections we will identify the expressions of P and Q for each moisture event. Finally, generic expressions of P and Q for any event will be proposed.

Moistening

All following expressions are written for

  1 ,   i i t t t
; i k will be simplified into k in case of no possible mix-up. During moistening phases, we assume that stress of box k (previous active box) is progressively (linearly) transferred to the new

active box 1  k :     1 1 1 ) ( 1 ) ( 1 ,                                     i i i i k i k k i i k k i j j t t t t t t S t t t t t k k j t        (12)
Here we supposed that all external stress S  occurring during time step is fully supported by box 1  k , in other words that 1  k is the active box during all time step. In case of moistening, this assumption has no consequence on elastic behavior: indeed even if k bears some part of ΔS, it will be transferred later

during t  to 1  k
with all other stresses. Strain transfer (eq.13) is directly deduced from linear elastic behavior of each box (eq.6).

            1 1 1 1 1 1 1 1 ,                                          i i i i k k i k k k k i i k k k k k i j j t t t t t t S J t J t t t t J t J t k k j t          (13)
Total strain in the hygrolock spring (eq.14) is obtained by adding strain changes in boxes i k and 1  i k , and all the possible internal strains stored in the other boxes.

            t t t S J J J E J t J t J t t t E i k i k k k i k k j i j j k k k k k k j i j k k                                     1 1 1 , 1 1 1 ; 1 ) (        (14)
By derivation of expression ( 14), we confirm that total strain rate is constant under the assumption of linear external stress evolution. Then we deduce the set of equations ( 17) giving E  during moistening:

                  1 1 k i k k k J Q J J P S Q P E  (15)
Expression ( 15) is rewritten in a more general form (eq.16), to be used whatever the moisture event.

                          1 1 i k i i k i k i k J Q J J P S Q P E  (16)

Drying

During drying hygrolock spring blocks the stress in the upper box. Then stresses are not transferred from box i k to box 1  i k . Thus the only variation in stress state is due to external increment S  (eq.17). Once again, we suppose that external stress is fully supported by 1

1    i i k k .     i i k k t t t S t        1 1   (17)
In expression (17), i k 1   is often equal to 0, because all stresses would have been transferred in i k  during previous moistening; exception occurs  

0 1   i k  if
simulation has started with an initial stress state.

Expression ( 17) is then injected into (6) to give the expression of strain in the new active box:

        i k i k i i k k k k k t t t S J t t t S J t J t                         1 1 1 1 1 1 1     (18)
Whereas former active box doesn't see any change, reproducing the blocking effect.

 

i k k k k k J t J t      ) ( (19)
The total strain is deduced from expressions (18-19):

      i k i k k j i j j k k k k t t t S J E J J J t E                     1 ; 1 1 1    (20)
Once again, total strain rate on the hygrolock spring is constant. As for moistening, we provide general behavior under the (eq.11) form.

               1 0 i k J Q P S Q P E (21)

Statu quo

Statu quo is the simplest case: the active box bears the total stress increment S  (eq.22) and no stress transfer occurs:

    i i k k t t t S t        (22)     i k i k k k t t t S J J t        (23)     i k i k j i j j k k t t t S J E J J t E                  (24)
which leads obviously to classical Hookean behavior:

                  1 0 i k i k J J Q P S Q P E (25)
Drying and statu quo are equivalent in terms of stress-strain relation; they are not from internal stress transfer point of view (see expression (27) below).

General behavior or hygrolock spring

To sum up we give hygrolock spring constitutive law whatever the event (eq.26), by using parameter  , defined in Table 1, which is equal to 1 in case of moistening, 0 otherwise.

                     1 1 i k i i k i k i k J Q J J P  (26)
We also give expression for the stress transfers (eq.27) using parameter П equal to 0 during statu quo, and 1 otherwise:

                               0 if 1 1 1 1 1 i i k i i k i i k i i k i i k S      (27)
Note that restriction of expression ( 27) in the case of statu quo needs to be correctly implemented in calculation program because

1   i i k k
.

Simulations of hygrolock spring

The previous equations have been implemented into the interpreted programming environment Scilab (www.scilab.org). A simulation test shown on Figure 5 is commented next. Stress and moisture content cycling are imposed, and strain is calculated. Significant events occur during the loading history, they are numbered from 1 to 9 on Figure 5:

1. HygroLock spring is wetted under load, compliance and strain increase, 2. Spring is dried, compliance decreases but strain is blocked, 3. Higher moistening leads to strain increase, 4. Spring is unloaded at low humidity, remaining strain has been blocked during the previous drying, 5. New moistening phases allow to recover strain completely, 6. Now spring is loaded wet, 7. Unloading at dry state ; again strain has been blocked in the upper boxes, 8. Small moistening allows recovering partly residual strain, 9. A new moistening to the condition when loaded allows complete recovery.

Incremental formulation of a hygrolock Kelvin-Voigt body

Here we explain the behavior of a hygrolock spring when assembled in parallel to a Newtonian dashpot, in order to form a hygrolock Kelvin-Voigt body (Figure 6).

The Kelvin-Voigt body is put under stress σ and strain ε ; E and S are still the global strain and stress of the hygro-lock spring. In this specific case

E  ε . σ σ S S ε = E
      t t S t ε σ     (28)
We assume that external stress σ varies linearly during time step

t  .     1 ,         i i i i t t t t t t t σ σ σ (29)
This allows solving easily the differential equation ( 28). The solution is detailed next for each moisture event.

Moistening

Stress transfers in the hygro-lock spring remain the same that in elastic case (see expression ( 27)). General equation ( 28) of the Kelvin-Voigt body is first multiplied by

1  k J to obtain:       t J t S J t k k σ ε 1 1       ( 30
)
where

  1   k J
is a characteristic time of the Kelvin-Voigt body.

Then expression of   t S is modified to make E  ε appear into the ordinary differential equation (ODE). According to characteristic equations of hygro-lock spring (eq.6-8), we can write:

                             1 , 1 1 1 1 1 k k j k k i j k n j j k k t t J t J t S J     (31) Whereas expression of   t E gives:               n j j k k j n j j j t J J J t J t E 1 1 1 1   (32)
Subtraction of expressions ( 31) and (32) gives:

        t J J J J J J J J t J J J t E t S J k k k k k k k j i j k k j k n j j k k j k    1 1 1 1 , 1 1 1 1 1 1 1 1                                      (33)
This last expression is used to modify the ODE (30):

              t t t J J S J t t t J J t t t J J J J t t t J J t t i i k k k i i k i k i k i i k k k n j i j j k i k i k                                     1 1 1 1 1 1 1 1 1 ε σ σ σ σ ε ε  (34)
which is rewritten into the general form:

                                 t J J t J D S J C t t D C t t i k k k k i i i k i   1 1 1 σ ε σ ε ε  (35)
Solution of equation ( 35) is given by:

         i i t t i i t t e t t D D C e t                   ε ε 1 (36) General expression of i i ε ε ε    1
is thus found to be:

                                t e t D C i t t i 1 1 ( 37
)
We finally get expression (38) to describe the general behavior of hygro-lock Kelvin-Voigt body in case of moistening.

                            1 1 i k i k i k i k i i J Q J J S P Q P     σ σ ε (38)

Drying (or statu-quo)

Again, expressions of stress transfers are taken from hygro-lock spring (eq.27).

General equation ( 28) of the Kelvin-Voigt body is now multiplied by

1  k J to obtain:      t J t S J t k k σ ε 1 1        (39) where   1    k J
is a new characteristic time. We assume next that     , i.e that viscosity varies together with compliance: further explanations are given in the following. Using the same calculation steps, the ODE is changed into:

                           t J D S J C t t D C t t k i i i k i σ ε σ ε ε 1 1   (40)
Solution of the ODE (eq.37) is used again to provide the drying behavior of hygrolock Kelvin-Voigt unit:

                  1 i k i i J Q S P Q P   σ σ ε (41)
The same expression is found for the statu quo case. However actualization of stresses within the boxes is slightly different (see expression ( 27)).

General behavior of hygrolock Kelvin-Voigt body

Expression (42) the general equation of the hygrolock Kelvin-Voigt unit to be implemented in numerical code.

                             1 1 i k i k i k i k i i J Q J J S P Q P     σ σ ε (42)
Note that internal stresses should be correctly updated during calculation, by using the same equations as for hygrolock spring (eq.27).

Theoretical material for simulations

Hygrolock simulations were performed with two different models: the formulation described in [START_REF] Husson | A finite element model for shape memory behavior[END_REF] (integral method) and the one presented in this paper (discrete method or box method). Wood modeling requires a large number of material parameters: variability of natural material, orthotropy, time-dependence and moisture related behavior make it very difficult to obtain consistent experimental parameters. In order to compare these simulations over a large range of humidity levels and times, a theoretical material is here proposed for simulation purpose; suitability of this data set is assessed by presenting a literature review on time-dependent effects in wood. Another method found in literature consists in using up-scaling homogenization methods to provide a generic orthotropic time-dependent behavior [START_REF] Eitelberger | Multiscale prediction of viscoelastic properties of softwood under constant climatic conditions[END_REF][START_REF] Eitelberger | Multiscale prediction of viscoelastic properties of softwood under constant climatic conditions[END_REF]).

Hygrolocks models need only to be fitted on viscous behavior at constant humidity (for at least two different moisture levels); complex behavior at varying humidity is directly deduced from the hygrolock formulation, without further parameters. Thus we propose a theoretical material defined by a series of creep curves given by power laws. Hygrolock model parameters is calculated by fitting the model on the creep curves over a time range of 1 to 8 of log(time in sec).

We give reference parameters for two anatomical directions, one is the longitudinal direction (L), i.e. fibers direction in wood, and the second one is a mean transverse (TR) behavior (for mean radial and tangential anatomical directions of wood). Elastic rigidity (one over compliance) at reference moisture content 12 . 0  h are chosen to be 10 GPa for L, and ten times lower in TR (1 GPa); these elastic parameters are moisture-dependent, according to expression (43) proposed by [START_REF] Guitard | Mécanique du matériau bois et composites[END_REF].

      12 . 0 1 1 1 % 12    h J h J X X X  (43)
where superscript X stands for the anatomical direction (L or TR), and We define two reference creep curves by using 2-parameters power laws:

  3 . 0 , . 0  h is
                      p t t J t J  1 0 0 ( 44
)
where  

0  t J
is the instantaneous compliance, p and  ~are 2 parameters defining the time-dependent behavior. Power laws were shown to be suitable for describing wood creep over a large range of directions, moisture content, and time (see [START_REF] Bardet | Modelling the transverse viscoelasticity of green wood using a combination of two parabolic elements[END_REF] 2 gives all parameters for the reference curves shown in Figure 8. As shown in the literature [START_REF] Gril | Using wood creep data to discuss the contribution of cellwall reinforced material[END_REF]), Dhoula et al. (2009)) value of parameters p and  ~decrease with moisture content.

 ~[sec] 1E+07 1E+05 1E+07 1E+05
The hygrolock model used next leans on a hygrolock Generalized Kelvin-Voigt The compliances of our GKV model are not independent variables: a good approximation of the power-law creep response is given by a GKV with a series of viscous compliances j J that fulfill expression (45).

         j j log . log J (45) with:   1 log 1   j j /τ  (46)
and where  and  are two linear functions of h, identified at a given humidity:

2 1      h (47) 2 1      h (48)
Here  and  were calculated at h = 0.05 and h = 0.25 by best fitting of GKV with reference power laws (by simulating simple creep test at constant humidity).

Then coefficients (time to double the elastic compliance) of equivalent power law decreases with increasing of moisture content (cf. Table 2 and Figure 8). Thus for a given Kelvin-Voigt body, rigidity and viscosity are equally modified by humidity. Alternative (equivalent) formulations with characteristic times modified by humidity would be possible, but this was not the choice made here.

To summarize: from two theoretical creep power-laws at low (0.05) and high Small corrections could be done on all parameters to fit better a given wood material, but we have here a consistent generic set of time-dependent parameters for L and TR directions to perform numerical simulations and compare different models. 2 andTable 3.

Simulations Stress-strain-moisture cycle

Various simulations have been performed using both the formulation presented in this paper and the integral method described in Dubois ( 2012), using same input data (Table 2-3) and constant time step of 300 seconds.

A theoretical creep test simulation at varying moisture content is shown on Mixing parallel and series configuration (eq.7,8) is here of great interest. The general equation ( 11) can be rewritten in a reciprocal form:

ε σ       1 1 Q P Q (49)
where P Q 1  is a hygroscopic pre-stress and 1  Q the rigidity, which allows to simulate also time-dependent relaxation test as shown in 11. Again, the two models give very similar results. This last result was not so obvious considering that discrete hygrolock spring has been built by reallocating stresses between the boxes. (2005)).

The discrete hygrolock model proposed here relies on the use of "boxes" corresponding to moisture levels, and on stress transfer between boxes according to external moisture events. The originality of considering both total strain and total stress in the hygrolock as the sum of strain (resp. stress) of each box led to a complete formulation for creep and relaxation. Although the fundamental hypotheses are quite different, the discrete model gives very similar results

to the previous formulation by [START_REF] Dubois | Modeling of the viscoelastic mechano-sorptive behavior in wood[END_REF].

In order to implement the present hygrolock model, one should focus on equations (26,27) for the hygro-elastic part, and (42,27) for the time-dependent GKV model.

The response of the discrete model slightly depends on the box number, but appears to be robust: unlike the Leicester model, hygrolock model tends to viscoelastic creep when the amplitude of moisture variations tends to zero.

The GKV model has also been used to simulate a time-dependent relaxation test with complex moisture variations. In the future this model could be used to simulate mixed conditions (alternate strain and stress control) likely to occur in wooden assemblies.

The hygrolock time-dependent model presented in this paper seems to be very interesting for wood at high humidity. Indeed, [START_REF] Montero | Influence of hygromechanical history on the longitudinal mechanosorptive creep of wood[END_REF] showed that samples which creep at high humidity level and samples subjected to previous moisture cycling converge to the same creep limit. In other words for high moisture levels, mechano-sorptive behavior doesn't give much additional deformation compared to stable viscoelastic response. However a "++ effect", i.e.

an increase of apparent strain rate during new moistening, can still be observed during experimentation on wood. We interpret this phenomenon as the combination of two effects: the first is due to blocking/releasing of hygrolock during new moistening, and the second is the dependency of delayed compliance to moisture content which has been shown to be non-linear with moisture content.

On the other hand it is quite evident that hygrolock model is not sufficient for low humidity domain where we still observe a significant amount of additional mechano-sorptive deformation compared to creep at constant humidity. This strong difference between the two domains (high and low moisture content) evidenced by [START_REF] Montero | Influence of hygromechanical history on the longitudinal mechanosorptive creep of wood[END_REF], is a real difficulty for models to fit both:

while Leicester model seems to reproduce the additional strain at low moisture content, hygrolock model seems to fit better the high moisture content domain.

Our interpretation is that wet state of wood is closer to the initial state within the tree, whereas in dry state wood exhibits specific behavior, maybe due to a specific initial state: the presence of high internal stresses initially blocked in dry wood could be a reason for significant additional creep to occur during the following moisture variations. Conversely, in the wet domain, wood is most closer to an initial free of stress state, for which our hygrolock model is qualitatively satisfactory.

In near future hygrolock formulation could be used to model the effect of high internal stresses as initial conditions for dry wood.

Figure 2 :

 2 Figure 2: Creep trajectories according to Hunt's representation, giving the total strain (  ) corrected for free swelling strain (

Figure 3 :

 3 Figure 3: Schematic representation of boxes within a hygrolock spring. The hygroscopic range of wood (from totally dry wood h 0 = 0. to fiber saturation point h n ≈ 0.3) is divided into n boxes. Current moisture level gives the "active box". As for all wood species and orthotropic directions, box compliance increases with moisture content.

Figure 4 :

 4 Figure 4: Discretisation of time and moisture along the moisture curve h(t). (a) Moistening with zoom; (b) drying; (c) jump of 2 boxes during dt. Thick black bars: active moisture content during Δt (which is moisture content at t i+1 ). Dotted black bars: characteristic moisture h k(i+1) of active box during Δt. Thick gray bars: new active box after interpolation of original data (to avoid jump (c)). Gray vertical bars: coloration from box 1 to active box.

Figure 5 :

 5 Figure5: simulation test of the hygrolock spring. Stress and moisture content cycling is imposed, strain is calculated. Various significant events are numbered from 1 to 9.

Figure 6 :

 6 Figure 6: Hygrolock Kelvin-Voigt body, ε and σ are external strain and stress on the body. E and S are the external strain and stress applied on the hygrolock spring.

  the moisture content of wood. The rigidity dependence to moisture parameters / % moisture content]), are also given in[START_REF] Guitard | Mécanique du matériau bois et composites[END_REF].

  , Dlouha et al.(2009)), with m values around 0.1-0.4 (Dhoula et al. (2008) (2009)). 20 We give two reference curves for each anatomical direction, one at low moisture content (0.05) and one at very high moisture content (0.25) to cover common values of moisture content. Table



  model (GKV), with parameters fitted on the previously given creep curves. The 1D GKV model used here consists of a series of: a hygrolock spring of equally spaced in log time scale from 10 sec to 3 years (1 to 8 in log scale) (see eq.46), associated with delayed compliances distribution GKV) is represented in Figure7.

Figure 7 :

 7 Figure 7: Hygrolock Generalized GKV: hygrolock spring of compliance distribution J 0 is associated in series with m hygrolock Kelvin-Voigt bodies of compliances distribution J j and characteristic time τ j .

( 0 .

 0 25) moisture content, a complete material behavior formulation based on the hygrolock GKV model is determined at any moisture content level. Creep response of such model at intermediate moisture levels is shown on Figure 8. In order to discuss the relevance of this theoretical material for creep modeling, we show on Figure 9 the relative delayed compliance     shown by experimental results of[START_REF] Schniewind | Über den Einfluss von Feuchtigkeitsänderungen auf das Kriechen von Buchenholz quer zur Faser unter Berûcksichtigung von Temperatur und Temperaturänderungen[END_REF],[START_REF] Nakano | Analysis of creep of wood during water adsorption based on the excitation response theory[END_REF] andMatar (2003), longitudinal relative compliance increases with humidity, and this effect is amplified by time. After 1 day, our material model exhibits relative creep of 1.4 at air-dry condition  

Figure 9 :

 9 Figure9: Relative creep compliance J(t)/J(t=0) versus moisture content (0.05-0.25) for longitudinal direction at different times; simulation results using parameters of Table2 and

Figure 10 .Figure 10 :

 1010 Figure10. The two models give very similar results, but small discrepancies might occur for hazardous choice of moisture boxes, as shown on strain curve from Figure10. For instance, if active box is not perfectly centered around the current value of moisture content at any time, then the solution given by discrete model slightly moves away from integral model. This is the consequence of regular division of hygroscopic range of wood into n boxes.

Figure 11 :

 11 Figure 11: relaxation tests at varying moisture content for a hygrolock GKV model; (plain curve)and integral method by Dubois (2012) (dashed line); maximum difference between the two curves less that 0.2%

Aknowledgments:

  The authors gratefully acknowledge the financial support of the Austrian Research Promotion Agency (FFG, project number 815234/12791) and the wood industry partnership "Building with Wood" within CEI-Bois for supporting this research work within project " MechWood ". COST (European Cooperation in Science and Technology) action IE0601 (Wood sciences for cultural heritage) is also greatly acknowledged for having supported cooperation among authors of this paper and others researchers. Special thanks to Cédric Montero for his fruitful comments.

  . In this case the viscosity depends

	on value h d d h d t . Modeling both drying and wetting in wood leads to use absolute t d

Table 2 :

 2 theoretical parameters for creep power law at two different moisture content level and anatomical directions.

		Longitudinal	TRansverse
	h	0.05	0.25	0.05	0.25
	J 1 [MPa] /	11050	8050	1210	610
	J [MPa -1 ] 9.0497E-05 1.2422E-04 8.2644E-04 1.6393E-04
	p	0.3	0.15	0.2	0.15

Table 3 :

 3 Parameters for moisture-time-dependent behavior of hygrolock GKV model at variable moisture content, obtained by fitting it on two reference creep curves (power law).

			Longitudinal TRansverse
		1	-0.9624	-0.3574

  The use of moisture content levels as proposed by[START_REF] Gril | Une modélisation du comportement hygro-rhéologique du bois à partir de sa microstructure[END_REF] allowed a formulation with internal variables: this is of particular interest for not having to store the entire moisture history for knowing the highest moisture content under load (see Dubois

	Conclusion
	For modeling mechanosorption in wood, models using varying compliance with
	d	h d	t	have been widely used: although they are quite convenient in terms of
	programming, they require additional parameters and they produce infinite strain
	during creep when  t h d d	tends to +∞.

The hygrolock model presented in this paper aims at giving a stable formulation for the mechanical behavior of wood at varying moisture content.