Meng Wu
email: meng.wu@hfut.edu.cn

Yicao Wang
email: yicwang@hhu.edu.cn

Bernard Mourrain
email: bernard.mourrain@inria.fr

Boniface Nkonga
email: boniface.nkonga@unice.fr

Changzheng

Changzheng Cheng
email: changzheng.cheng@hfut.edu.cn

Convergence rates with singular parameterizations for solving elliptic boundary value problems in isogeometric analysis

Keywords: singular parameterizations, convergence rates, elliptic boundary value problems, isogeometric analysis

come

Convergence rates with singular parameterizations for solving elliptic boundary value problems in isogeometric analysis

Introduction

Isogeometric analysis (IGA) [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] is a numerical method to solve geometric partial differential equations. The theory of approximation and error estimates has been developed in [START_REF] Bazilevs | Isogeometric analysis: approximation, stability and error estimates for h-refined meshes[END_REF][START_REF] Tagliabue | Isogeometric Analysis and error estimates for high order partial differential equations in fluid dynamics[END_REF], where the theory relies on the regularity of parameterizations, i.e., there are no singularities of the underlaying parameterizations.

Although singularities can be avoided by using some constraints [START_REF] Cohen | Analysis-aware modeling: Understanding quality considerations in modeling for isogeometric analysis[END_REF][START_REF] Xu | Constructing analysis-suitable parameterization of computational domain from CAD boundary by variational harmonic method[END_REF][START_REF] Xu | Parameterization of computational domain in isogeometric analysis: Methods and comparison[END_REF] in some cases, some singularities may stem essentially from the geometry or topology of the physical domain and cannot be avoided. For instance, in [START_REF] Lu | Circular element: Isogeometric elements of smooth boundary[END_REF], for a topologically circular domain, the parameterization is singular at 'a converging point'. In [START_REF] Wu | Hermite type spline spaces over rectangular meshes with arbitrary topology accepted by Communications in Computational Physics[END_REF], singularities appear at extraordinary vertices if the parameterization is asked to have C 1 continuity. In this paper, we present the convergence rates for solving elliptic boundary value problems with singular parameterizations.

Recently, there are works focused on singular parameterizations in the IGA. In [START_REF] Takacs | Existence of Stiffness Matrix Integrals for Singularly Parameterized Domains in Isogeometric Analysis[END_REF][START_REF] Takacs | H 2 regularity properties of singular parameterzations in isogeometric analysis[END_REF], the regularity of test functions on physical domains is discussed. Afterwards, the authors continue to work on this topic. One of the authors presents results on the approximation properties of isogeometric analysis function spaces in [START_REF] Takacs | Approximation properties of isogeometric function spaces on singularly parameterized domains preprint[END_REF] and construct smooth isogeometric function spaces [START_REF] Takacs | Construction of smooth isogeometric function spaces on singularly parameterized domains[END_REF] on a singularly parameterized domain. In addition, in enriched isogeometric analysis (EIGA) [START_REF] Jeong | Mapping techniques for isogeometric analysis of elliptic boundary value problems containing singularities[END_REF][START_REF] Oh | Enriched isogeometric analysis of elliptic boundary value problems in domains with cracks and/or corners[END_REF], the authors introduce a mapping technique which is an additional singular mapping patch used to generate singular functions in physical domains with cracks and/or corners. In [START_REF] Oh | Enriched isogeometric analysis of elliptic boundary value problems in domains with cracks and/or corners[END_REF], the error estimates of this mapping technique on an elliptic boundary value problem are presented. So far, singular parameterizations are discussed independently for solving elliptic boundary value problems with smooth weak solutions or non-smooth weak solutions. In the present work, the results are based on a weak solution space F(Ω) which contains all possible weak solutions of an elliptic boundary problem with smooth coefficients. The main contributions are listed in the following. 1. For functions in F(Ω), we present the approximation estimates with the L 2 (Ω)norm and the H 1 (Ω)-seminorm. Different from previous work, approximation errors are estimated locally, i.e., the impact of the singularities of the parameterization is considered. 2. Based on the approximation error estimates, we discuss the convergence rates of solving elliptic boundary problems with singular parameterizations. By introducing a controllable singular parameterization, optimal convergence rates are reached by setting singularities even if the optimal convergence rates cannot be reached by the finite element method (FEM). Moreover, convergence rates are discussed in the IGA with singular parameterizations.

The paper is structured as follows. Approximation error estimates are presented in Section 2. After that, in Section 3, we discuss convergence rates for solving elliptic boundary problems based on the error estimates of Section 2. Finally, in Section 4, we conclude this paper with a discussion of possibilities for future works.

Approximation Error Estimates of Singular parameterizations

In this section, as background, parameterizations and singularities of parameterizations are introduced in Section 2.1. Then, in Section 2.2, by Weyl's lemma, a set of functions F(Ω) called a weak solution space is introduced, where F(Ω) contains all the weak solutions of elliptic boundary value problems with smooth coefficients on a physical domain Ω. In addition, a Hermite projection Π Â H p h is introduced as a preparation for presenting the L 2 (Ω)-norm and H 1 (Ω)-seminorm approximation error estimates of singular parameterizations in Theorems 2.3 and 2.4 of Section 2.3.

Parameterizations and their singularities

Definition 2.1. A parameterization P(s, t) is a map from a parametric domain to a physical domain, i.e. P : Σ -→ Ω;

(s, t) -→ (x(s, t), y(s, t)).

For discussing the convergence rate for solving elliptic boundary value problems, the parameterization P(s, t) satisfies the following properties:

(1). P(s, t) is a one-to-one map.

(2). To solve the elliptic boundary problem, P(s, t) satisfies the H 1 integrability assumption [START_REF] Takacs | Existence of Stiffness Matrix Integrals for Singularly Parameterized Domains in Isogeometric Analysis[END_REF], i.e., for any f (s, t) in a given isogeometric function space, f If J(s 0 , t 0 ) = 0, then (s 0 , t 0 ) is called a singularity of P(s, t).

In the following, we suppose that J(s, t) ≥ 0 and the singularities are isolated. Furthermore, at a singularity (s 0 , t 0 ), P(s, t) satisfies ∂x ∂s (s 0 , t 0 ) = 0, ∂y ∂s (s 0 , t 0 ) = 0, ∂x ∂t (s 0 , t 0 ) = 0, ∂y ∂t (s 0 , t 0 ) = 0.

In Example 2.1, there is an H 1 -parameterization of a geometry with a smooth boundary.

Example 2.1 (A geometry with a smooth boundary). In this example, parameterizations of a geometry with a smooth boundary are presented. The parametric domain Σ is [0, 1] × [0, 1]. To fit the smooth boundary in Figure 1(b), P(s, t) is singular at the corner points P 1 (0, 0), P 2 (1, 0), P 3 (0, 1), P 4 (1, 1). Although the tangents of the s-curve and the t-curve are degenerate, the limits of the unit tangents of the s-curve and the t-curve at these points are collinear.

-∇ • (p(x, y)∇u) + q(x, y)u = f (x, y) in Ω, u = 0 on ∂Ω. (1)
Based on Weyl's lemma [START_REF] Strang | An analysis of the finite element method Prentice Hall[END_REF], if p(x, y), q(x, y), f (x, y) are analytic at the interior points of Ω, (i.e., on a neighborhood of an interior point P , their Taylor expansions converge to themselves), the weak solution u is analytic at these points. If ∂Ω is analytic, u is analytic on Ω. For Ω with non-smooth boundaries, non-analytic points are among these non-smooth boundary points. In the following, there is an example of a weak solution of (1).

Example 2.2. Let Ω be the physical domain shown in Figure 2. There is a corner P . The weak solution u of the elliptic boundary value problem (1) with p(x, y) = 1 and q(x, y) = 0 is non-analytic at P . If θ = απ, in [START_REF] Strang | An analysis of the finite element method Prentice Hall[END_REF], then the weak solution u around the corner P is

u(r, θ) = ∞ j=1 α j r v j φ j (θ) + ∞ j=1 ∞ =0 f j [( + 2) 2 -v 2 j ] -1 φ j (θ)r +2 , (2)
where

φ j (θ) = 2 απ sin v j θ, v j = j/α; f j (r) = απ 0 f φ j (θ)dθ, f j (r) = f j r , Ω 0 = {(r, θ)|0 < r < r 0 , 0 < θ < απ} ⊂ Ω, i.e.
, it has a singularity of type O(r 1/α ). Except for at the point P , the solution u is analytic. In this case, if 1/α / ∈ N, then u ∈ H σ (Ω) where 1 < σ < 2 and r = x 2 + y 2 . To describe all the weak solutions on Ω, we define a function space called a weak solution space,

F(Ω) = {f : Ω -→ R : f (x, y) is analytic on Ω o } ⊂ C ∞ (Ω o ).
The weak solutions on Ω of the elliptic boundary value problem are in F(Ω). Let P : Σ -→ Ω be a parameterization of Ω, where P(s, t) ∈ C (Σ). Thus,

F • P ⊂ C (Σ o ).
Let Q h be a regular mesh, i.e., it is a rectangular mesh and the valences of its interior vertices are 4, defined on a parametric domain Σ, where h is the maximum of the diameters of its cells. In Figure 3, there are three regular meshes. H p h (Q h ) is the spline space with bidegree (p, p) defined on Q h , and

H p h (Q h ) ⊂ C (Σ)
, where p = 2 + 1. In the following, we define an interpolation projection from

F(Ω) to H p h (Q h ). Associated to each vertex v k ∈ V(Q h ), we can construct ( + 1) 2 Hermite functions {H v k i,j (s, t) : i, j = 0, 1, • • • , }, where H v k i,j (s, t) satisfies ∂ r 1 +r 2 H v k i,j (v l ) ∂s r 1 ∂t r 2 = δ k,l δ r 1 ,i δ r 2 ,j , (3)
where

r 1 , r 2 = 0, 1, 2, • • • , , v k , v l ∈ V(Q h ), δ i,j = 0 if i = j, otherwise δ i,j = 1.
Following [START_REF] Deng | Dimensions of spline spaces over T-meshes[END_REF], it is easy to prove that

∪ v k ∈V(Q h ) {H v k i,j (s, t) : i, j = 0, 1, • • • , } is a set of basis functions of H p h (Q h ). Define a V(Q h )-based Hermite projection Π V(Q h ) H p h :C ( Σ) -→ H p h (Q h ) (4) such that Π V(Q h ) H p h (f (s, t)) = v∈V(Q h ) ∂ i+j f (v) ∂s i ∂t j H v i,j (s, t).
However, for the weak solution space F(Ω), F • P ⊂ C (Σ o ), in other words, the Hermite projection ( 4) is not well defined on the boundary vertices of Q h . To solve this problem, for each boundary vertex v b of Q h , there is a point q ∈ Σ o such that we can construct a set of linear independent splines {H q i,j (s, t) : i, j = 0, 1, • • • , } associated with q. In addition, {H q i,j (s, t)

: i, j = 0, 1, • • • , } is linearly equivalent to {H v b i,j (s, t) : i, j = 0, 1, • • • , }, i.e.,
for each H q i,j (s, t), it is a linear combination of {H v b i,j (s, t) : i, j = 0, 1, • • • , }; the reverse is also true. Here we present this theorem.

Theorem 2.1. Let Σ be a parametric domain. Q h is a regular mesh on Σ. For each boundary vertex

v b k ∈ V(Q h ), there is a point q k ∈ Σ o such that {H q k i,j (s, t) : i, j = 0, 1, • • • , } are linearly equivalent to {H v b k i,j (s, t) : i, j = 0, 1, • • • , }, and H q k i,j (s, t) satisfies ∀q l ∈ V(Q h ) ∪ {q k } -{v b k } ∂ r 1 +r 2 H q k i,j (q l ) ∂s r 1 ∂t r 2 = δ k,l δ ir 1 δ jr 2 .
Proof. The proof is divided into two parts. In the first part, we prove the existence of q k and construct H q k i,j (s, t). In the second part, prove the linear relationship between

{H q k i,j (s, t) : i, j = 0, 1, • • • , } and {H v b k i,j (s, t) : i, j = 0, 1, • • • , }. 1.
The existence of q k : Consider the matrix

H v b k (s, t) =   ∂ r 1 +r 2 H v b k i,j ∂s r 1 ∂t r 2 (s, t)   ( +1) 2 ×( +1) 2
.

Its rows and columns are arranged by the (i, j)-order and (r 1 , r 2 )-order, respectively, where the (i, j)-order and the (r 1 , r 2 )-order are (0, 0)

≺ (1, 0) ≺ • • • ≺ ( , 0) ≺ (0, 1) ≺ (1, 1) ≺ • • • , ≺ ( , 1), ≺ • • • ,≺ (0, ) ≺ (1, ) ≺ • • • ≺ ( , ). Then, on a cell C of Q h with v b k as one of its corner point, det H v b k (s, t)| C is a polynomial and det H v b k (v b k ) = 1 = 0 because H v b k (v b k ) is an identity matrix. Thus, there is a neighborhood of v b k , U (v b k ), such that for any q ∈ U (v b k ) ∩ C, det H v b k (v b k ) = 0. Take q k ∈ U (v b k ) ∩ C and construct the following functions in H p h (Q h ) H q k i,j (s, t) = k,l=0 C k,l H v b k k,l (s, t). (5)
and

C = (C 0,0 , C 1,0 , • • • , C ,0 , • • • , C 0, , • • • , C , ) T satisfy C = (H v b k (q k )) -1 e i,j , (6)
where

(H v b k (q k )) -1 exists because det H v b k (q k ) = 0; Additional, e i,j = (0, • • • , 0, 1, 0, • • • , 0)
T is a column vector, i.e., the (i, j)-element is 1 associated with the (i, j)-order. Other elements are zeros. Computing directly, we obtain

∂ r 1 +r 2 H q k i,j (q l ) ∂s r 1 ∂t r 2 = δ k,l δ i,r 1 δ j,r 2 , ∀q l ∈ V(Q h ) ∪ {q k } -{v b k }.
2. The linear relationship: Denote

H q k (s, t) = (H q k 0,0 (s, t), • • • , H q k i,j (s, t), • • • , H q k , (s, t)) H v b k (s, t) = (H v b k 0,0 (s, t), • • • , H v b k i,j (s, t), • • • , H v b k , (s, t))
By ( 6),

H q k (s, t) = H v b k (s, t)(H v b k (q k )) -1 , i.e., H q k (s, t)H v b k (q k ) = H v b k (s, t). Thus {H q k i,j (s, t) : i, j = 0, 1, • • • , } and {H v b k i,j (s, t) : i, j = 0, 1, • • • , } are linearly equivalent.
Based on this theorem, by replacing each boundary vertex

v b k ∈ V(Q h ) by q k ∈ Σ o , another set of interpolation points  of H p h (Q h ) is introduced and  ⊂ Σ o . Thus, an Â-based Hermite projection Π Â H p h can be defined such that Π Â H p h :F • P -→ H p h (Q h ),
where

Π Â H p h (f ) = q k ∈ Â ∂ i+j f (q k ) ∂s i ∂t j H q k i,j (s, t).
By the B-net method [START_REF] Deng | Dimensions of spline spaces over T-meshes[END_REF], the Bezier coordinates around q k are equivalent to the ones around v b k and they are independent of the Beźier coordinates around the other vertices of Â. Moreover, for each polynomial p(s, t) with bidegree (k, k)

Π V(Q h ) H p h (p(s, t)) = p(s, t),
where k ≤ p. Thus, we have Lemma 2.2.

Π Â H p h (p(s, t)) = p(s, t), p(s, t) ∈ P [k,k] .
where k ≤ p, P [k,k] is a polynomial space with bidegree (k, k).

Induced by H p h (Q h ) and P(s, t), one can define a projection

Π Â V h :F -→ V h u(x, y) → Π Â H p h (u • P(s, t)) • P -1 ,
where

V h = H p h (Q h ) • P -1 = span{H q i,j (s, t) : q ∈ Â, i, j = 0, 1, • • • , }
and {H q i,j (s, t) : q ∈ Â, i, j = 0, 1, • • • , } are the basis functions of H p h (Q h ). In the next section, we begin to measure the L 2 (Ω)-norm and the H 1 (Ω)-seminorm approximation error estimates of singular parameterizations.

2.3. L 2 (Ω)-norm and H 1 (Ω)-seminorm approximation error estimates of singular parameterizations To measure approximation error estimates of singular parameterizations, we introduce some symbols.

Classifying cells of Q

h : C s is a cell of Q h on Σ,
where α = (α 1 , α 2 ), ∂ α = ∂ α 1 1 ∂ α 2 2 , |α| = α 1 + α 2 . Especially, when p = 2, ||u|| L 2 (Ω) = ||u|| 0,2,Ω ; ||u|| H k (Ω) = ||u|| k,2,Ω (k ≥ 1); |u| k,2,Ω = |u| H k (Ω) .
In the following, the approximation error estimates with the L 2 (Ω)-norm and the H 1 (Ω)-seminorm are presented.

Theorem 2.3. Let Σ and Ω be a parametric domain and a physical domain.

P : Σ -→ Ω is a parameterization of Ω. u(x, y) ∈ F(Ω) ∩ L 2 (Ω) is a bounded function on Ω and u • P| Σr ∈ H k+1 (Σ r ), u • P| Σs ∈ H k +1 (Σ s ). Then, ||u -Π Â V h u|| L 2 (Ω) ≤ K 1 h k+1 |u • P| H k+1 (Σr) + K 2 h σ/2+k +1 |u • P| H k +1 (Σs) , ( 7
)
where h is the maximum diameter of the cells of the parametric mesh

Q h on Σ, V h = H p h • P -1 (p ≥ min{k, k }). J| Cs ∼ h σ Cs , σ = min{σ Cs }, K 1 = K 1 (P, k, Π Â H p h
) and

K 2 = K 2 (P, k , Π Â H p h
), i.e., K 1 and K 2 are independent of the mesh size h and u(x, y), but are dependent on the regularity of u • P, Π Â H p h and its parameterization P(s, t).

Proof.

||u -Π Â V h u|| 2 L 2 (Ω) = C C (u • P -Π Â H p h (u • P)) 2 |J|dsdt = Cr Cr (u • P -Π Â H p h (u • P)) 2 |J|dsdt + Cs Cs (u • P -Π Â H p h (u • P)) 2 |J|dsdt ≤ K 1 1 Cr Cr (u • P -Π Â H p h (u • P)) 2 dsdt + K 1 2 h σ Cs Cs (u • P -Π Â H p h (u • P)) 2 dsdt, = K 1 1 Cr ||u • P -Π Â H p h (u • P)|| 2 L 2 (Cr) + K 1 2 h σ Cs ||(u • P -Π Â H p h (u • P))|| 2 L 2 (Cs) ,
where K 1 1 and K 1 2 are only dependent on P. Based on Lemma 2.2,

Π Â H p h p(s, t) = p(s, t), ∀p ∈ P [k,k] and Π Â H p h p(s, t) = p(s, t), ∀p ∈ P [k ,k ]
because k, k ≤ p. Thus,based on Theorem 5 in [START_REF] Ciarlet | General lagrange and hermit interpolation in R n with applications to finite element methods[END_REF] on each cell,

||u • P -Π Â H p h || L 2 (C) ≤ K 2 1 |u • P| H k+1 (C) h k+1 , if C is a regular cell; ||u • P -Π Â H p h || L 2 (C) ≤ K 2 1 |u • P| H k +1 (C) h k +1 , if C is a irregular cell, where u • P| Cr ∈ H k+1 (C r ) and u • P| Cs ∈ H k +1 (C s ), K 2 1 , K 2 2 are only dependent on k, k , Π Â H p h
. Thus, we obtain the error bound [START_REF] Lu | Circular element: Isogeometric elements of smooth boundary[END_REF], where K 1 , K 2 are only dependent on

P, k, Π Â H p h and P, k , Π Â H p h , respectively. Theorem 2.4. Let u ∈ F(Ω) ∩ H 1 (Ω) and F • P ⊂ W k+1 ∞ (Σ) (for p ≥ k ≥ 0). Then |u -Π Â V h u| H 1 (Ω) ≤ Ch k |u • P| k+1,∞,Σ , ∀u ∈ F • P, where C = C(P, k, Π Â V h ) and V h = H p h • P -1 , (p ≥ 1).
Proof. J dsdt

|u-u h | 2 H 1 (Ω) = Ω [( ∂(u -Π Â V h u) ∂x ) 2 + ( ∂(u -Π Â V h u) ∂y ) 2 ]dxdy = Σ ( ∂E(
E(s, t) = u • P -Π Â H p h (u • P) ≤ Σ ( ∂E(s, t) ∂s ) 2 ( ∂y ∂t ) 2 + ( ∂E(s, t) ∂t ) 2 ( ∂y ∂s ) 2 + ( ∂x ∂s ) 2 ( ∂E(s, t) ∂t ) 2 + ( ∂x ∂t ) 2 ( ∂E(s, t) ∂s ) 2 J dsdt ≤ |E(s, t)| 2 1,∞,Σ Σ ( ∂x ∂s ) 2 + ( ∂y ∂s ) 2 + ( ∂x ∂t ) 2 + ( ∂y ∂t ) 2 J dsdt = |E(s, t)| 2 1,∞,Σ (|s • P -1 | H 1 (Ω) + |t • P -1 | H 1 (Ω) ) ≤ C 1 |u • P -Π Â H p h (u • P)| 2 1,∞,Σ ,
where

C 1 = C 1 (P) because |s • P -1 | H 1 (Ω) < ∞ and |t • P -1 | H 1 (Ω) < ∞ by the H 1 -parameterization P(s, t). Moreover, by k ≤ p and Lemma 2.2, Π Â H p h p(s, t) = p(s, t), ∀p ∈ P [k,k] .
Thus, by Theorem 5 in [START_REF] Ciarlet | General lagrange and hermit interpolation in R n with applications to finite element methods[END_REF],

|u • P -Π Â H p h (u • P)| 1,∞,Σ ≤ ||u • P -Π Â H p h (u • P)|| 1,∞,Σ ≤ C 2 |u • P| k+1,∞,Σ h k ,
where

C 2 = C 2 (k, Π Â H p h
). Thus, we obtain the result. From the results in Theorems 2.3 and 2.4, the regularity of u • P on a parametric domain determines the approximation order. By the multivariate Faà di bruno formula [START_REF] Constantine | A multivariate faa di bruno formula with applications[END_REF], which is a formula for derivatives of a composite function, we obtain the following theorem.

Theorem 2.5. If u ∈ F(Ω) ∩ H k (Ω) and on each cell of Q h P(s, t) ∈ C ∞ × C ∞ , then u • P ∈ H k (Σ). If u ∈ F(Ω) ∩ W k ∞ (Ω), then u • P ∈ W k ∞ (Σ), where k ≥ k.
For instance, in Section 3.2, u • P(s, t) becomes smooth even if u ∈ H σ (Ω) by controlling the singularities of P(s, t), where 1 < σ < 2.

Convergence Rates of Solving Elliptic Boundary Problems with Singular Parameterizations

In this section, we discuss convergence rates of solving elliptic boundary problems with singular parameterizations. An elliptic boundary problem is a model problem, and we introduce it as a background in Section 3.1. To discuss the impact of singularities, the convergence rates are presented by a controllable singular parameterization in Section 3.2. Then, in Section 3.3, we discuss the convergence rate with a singular parameterization within the IGA framework.

Model problem: an elliptic boundary problem

For elliptic boundary value problems,

-∆u = f in Ω, u = 0 on ∂Ω, (8)
where Ω is a physical domain and f is assumed to be analytic on Ω. Its weak form reads:

find u ∈ V : a(u, v) = l(v), ∀v ∈ V, (9)
where V = {ω ∈ H 1 (Ω) : ω| ∂Ω = 0} is a Hilbert space endowed with the norm of V. Moreover,

a(u, v) = Ω ∇v T ∇udΩ, l(v) = Ω f vdΩ
where a : V ×V -→ R is a continuous, (strongly) coercive bilinear form and l : V -→ R is a continuous linear functional. In the framework of the FEM, the approximate solution u h is obtained by solving the following problem:

find u h ∈ V h : a(u h , v h ) = l(v h ), ∀v h ∈ V h , (10)
where V h ⊂ V is a nontrivial and finite dimensional subspace.

Here, we collect the results of the a priori estimate from the FEM from [START_REF] Strang | An analysis of the finite element method Prentice Hall[END_REF]. Let Ω be a polygonal domain, f ∈ L 2 (Ω), and u h be the FEM solution.

If the exact weak solution

u ∈ H σ (Ω) (σ ≥ 2), then ||u -u h || L 2 (Ω) ≤ Ch l+1 ||u|| H l+1 (Ω) , ||u -u h || H 1 (Ω) ≤ Ch l ||u|| H l+1 (Ω) , (11)
where l = min(p, σ -1).

2. If u ∈ H σ (Ω) (1 < σ < 2), then the error bound is ||u -u h || L 2 (Ω) ≤ C(h l+σ-1 + h 2p )||u|| H l+1 (Ω) , l = min(p, σ -1), ( 12
)
where p is the degree of the piecewise polynomial basis functions.

In the IGA, Equation ( 10) is computed on the parametric domain Σ, i.e.,

a(u h , v h ) = Σ ∂(u h , y) ∂(s, t) ∂(v h , y) ∂(s, t) + ∂(x, u h ) ∂(s, t) ∂(x, v h ) ∂(s, t) J -1 dsdt, (13)
l(v h ) = Σ (f • P)(v h • P)Jdsdt, where ∂(f, g) ∂(s, t) = ∂f ∂s ∂g ∂t - ∂f ∂t ∂g ∂s .
The integral in ( 13) is a defective, if P(s, t) has singularities. However, it converges when P(s, t) is an H 1 -parameterization. Here, this integral is evaluated using the Gp × Gp points tensor-product Gaussian quadrature on each cell of Q h . To guarantee the accuracy of the numerical integral, we take as many Gaussian points as possible in the following experiments. Moreover, in these experiments, we take

V h = {H • P -1 : H(s, t) ∈ H 3 h (Q h ), H| ∂Σ = 0}
. By these a priori estimates [START_REF] Takacs | Approximation properties of isogeometric function spaces on singularly parameterized domains preprint[END_REF] and [START_REF] Takacs | Construction of smooth isogeometric function spaces on singularly parameterized domains[END_REF], in the FEM, there is positive correlation between the convergence rate of numerical solutions and the regularity of the weak solution. Especially, when the weak solution u ∈ H 4 (Ω), the optimal convergence rates with the L 2 (Ω)-norm and the H 1 (Ω)-seminorm are reached in the case of p = 3. In the following section, we find that the optimal convergence rates are also reached, even if u ∈ H σ (Ω), 1 < σ < 2, by choosing a proper singular parameterization.

Experiments with controllable singular parameterizations

In this section, we discuss the convergence rates with the L 2 (Ω)-norm and the H 1 (Ω)-seminorm by introducing a controllable singular parameterization at the corner points. The controllable singular parameterization P δ (s, t) is defined as:

P δ :Σ -→ Ω (s, t) → (s δ cos t, s δ sin t) (14)
where the Jacobian J = s 2δ-1 .

When δ = 1, P 1 (s, t) has been used to overcome the loss of accuracy encountered when applying the standard schemes to two-dimensional elliptic problems in the crack problems in [START_REF] Yosibash | Superelements for the finite element solution of twodimensional elliptic problems with boundary singularities Finite Elements in[END_REF]. For P δ (s, t), the singularity can be controlled by different values of δ. In addition, in Lemma 3.1, we prove that P δ is an H 1 -parameterization.

Lemma 3.1. If H(s, t) ∈ H 3 h (Σ) and H(s, t)| ∂Σ = 0, then

H • P -1 δ ∈ H 1 (Ω), ∀i = 1, 2, • • • , n.
Proof. Denote P δ (s, t) = (x(s, t), y(s, t)) = (s δ cos t, t δ sin t), where (s, t) ∈ [0, 1] × [0, 1]. Its Jacobian is J = s 2δ-1 . In the following, we will check that, for each H(s, t), the L 2 -norm and the H 1 semi-norm are bounded. On the one hand,

||H • P -1 δ || 2 L 2 (Ω) = Ω (H • P -1 δ ) 2 dxdy = Σ (H(s, t)) 2 s 2δ-1 dsdt.
Then, ||H • P -1 δ || L 2 (Ω) is bounded. On the other hand,

|H • P -1 δ | 2 H 1 (Ω) = Ω (( ∂H • P -1 δ ∂x ) 2 + ( ∂H • P -1 δ ∂y ) 2 )dxdy = [0,1]×[0,1] ( ∂H ∂s , ∂H ∂t ) Ĵ ĴT ( ∂H ∂s , ∂H ∂t ) T J dsdt = [0,1]×[0,1] s 2 ∂H ∂s + δ 2 ∂H ∂t s dsdt,
where Ĵ = s δ cos t s δ sin t -δs δ-1 sin t δs δ-1 cos t .

Moreover, s|H(s, t) because H(s, t)| ∂Σ = 0. So, s| ∂H ∂t . Thus,

|H • P -1 | H 1 (Ω) is bounded. In conclusion, H • P -1 δ ∈ H 1 (Ω)
. By this controllable singular H 1 -parameterization, we discuss the optimal convergence rates with L 2 (Ω)-norm and H 1 (Ω)-seminorm on a physical domain with corners such as the one shown in Figure 2. To discuss the regularity of weak solution cases differently, P δ is rewritten as

P δ,α :[0, 2] × [0, 2] -→ Ω α (s, t) → (( s 2 ) δ cos( απ 2 t), ( s
2 ) δ sin( απ 2 t)), (15)
where the included angle of Ω α is απ, Ω = P δ (Σ). In Examples 3.1 and 3.2, we discuss physical domains with α = 2 and α = 3/2, respectively.

Example 3.1. [α = 2] The elliptic boundary value problem: where the physical domain is shown in Figure 4(b) is

-∆u = 3 4 (x 2 + y 2 ) -5 4 (1 + 5 3 x 2 + y 2 )y in Ω 2 , u = 0 on ∂Ω 2 , (16)
Ω 2 = {(x, y) ∈ R 2 : x = r cos(θ), y = r sin(θ), 0 < r < 1, 0 < θ < 2π}.
The exact weak solution of this elliptic boundary value problem is u(x, y) = (x 2 + y 2 ) -1/4 (1 -x 2 + y 2 )y and

u • P δ,α (s, t) ( s 2 ) δ/2 [1 -( s 2 ) δ ] sin( απ 2 t).
Thus,

u • P δ,α ∈ H δ 2 +1-(Σ s ) (∀ > 0); u • P δ,α ∈ H 4 (Σ r ); u • P δ,α ∈ W δ 2 - ∞ (Σ s ).
By Theorems 2.3 and 2.4, we have

Corollary 3.2. The L 2 (Ω)-norm approximation order is not less than min{[ 3δ 2 + 1 2 - ], 4}; the H 1 (Ω)-seminorm approximation order is not less than min{[ δ 2 -1 -], 3}. Furthermore, when δ > 7 3
, the optimal L 2 (Ω)-norm approximation order is reached; when δ > 8, the optimal H 1 (Ω)-seminorm approximation order is reached.

In the following, the convergence rates between the numerical solution and the exact solution of (16) are presented. In Figure 5(a), we present the numerical result with the L 2 (Ω)-norm by taking δ = 2.4. The optimal L 2 (Ω) approximation order is reached by Corollary 3.2 because δ > 7/3. From the numerical result in Figure 5(a), the optimal convergence rate (i.e., 4) for solving the elliptic boundary value problem ( 16) is reached. Taking δ = 8.1, the optimal H 1 (Ω)-seminorm approximation order (i.e., 3) is obtained because δ = 8.1 > 8. The numerical results are presented in Figure 5(b). In Figure 5(b), the optimal convergence rate (i.e., 3) for solving the elliptic boundary value problem ( 16) is reached. -∆u = 2xy( 32 9 (x 2 + y 2 ) -5/3 -11 9 (x 2 + y 2 ) -7/6 ) in Ω,

u = 0 on ∂Ω, (17)
where the physical domain shown in Figure 6(b) is

Ω = {(x, y) ∈ R 2 : x = r cos θ, y = r sin θ, 0 < r < 1, 0 < θ < 3 2 π}.
The exact solution of this elliptic boundary value problem is

u(x, y) = 2xy(x 2 + y 2 ) -2/3 (1 -x 2 + y 2 ). u • P δ,α (s, t) = ( s 2 ) 2δ 3 (1 -( s 2 ) δ ) sin(απt) ∈ H σ (Ω) (1 < σ < 2)
Thus ∀ > 0,

u • P δ,α ∈ H 2δ 3 +1-(Σ s ); u • P δ,α ∈ H 4 (Σ r ); u • P δ,α ∈ W 2δ 3 - ∞ (Σ s ).
By Theorems 2.3 and 2.4, we obtain the following Corollary 3.3. The L 2 (Ω)-norm approximation order is not less than min{[ 5δ 3 +

1 2 ], 4}, and the H 1 (Ω)-norm approximation order is not less than min{[ 2δ 3 -1 -], 3}. When δ > 2.1, the optimal L 2 (Ω)-norm approximation order is reached. When δ > 6, the optimal H 1 (Ω)-seminorm approximation order is reached. From Corollary 3.3, when we take δ = 2.2, the L 2 (Ω)-norm approximation order is optimal (i.e., 4) because δ > 2.1. In Figure 7(a), we present the numerical results of errors with the L 2 (Ω)-norm by δ = 2.2. By this numerical result, the optimal L 2 (Ω)norm convergence rate is reached.

Taking δ = 6.1, the optimal H 1 (Ω)-seminorm approximation order (i.e. 3) is reached because δ > 6. In Figure 7(b), the numerical errors with the H 1 (Ω)-seminorm for solving the elliptic boundary value problem [START_REF] Constantine | A multivariate faa di bruno formula with applications[END_REF] are presented.

0 0.5 1 1.5 2 2.5 3 3.5 -18 -16 -14 -12 -10 -8 -6 -4
-2 log(1/h) log(||e|| L 2 (Ω) ) δ=2.2, ||e|| L 2 (Ω) ||e|| L 2 (Ω) , Gp=5 ||e|| L 2 (Ω) , Gp=8 ||e|| L 2 (Ω) , Gp=15 y=-4x-4 (a) δ = 2.2, L 2 (Ω)-norm 0 0.5 1 1.5 2 2.5 3 3.5 -12 -10 -8 -6 -4 -2 0 2 log(1/h) log(|e| H 1 (Ω) ) δ=6.1, |e| H 1 (Ω) |e| H 1 (Ω) , Gp=5 |e| H 1 (Ω) , Gp=8 |e| H 1 (Ω) , Gp=15 y=-3x+1 (b) δ = 6.1, H 1 (Ω)-seminorm

Experiments with singular parameterizations in isogeometric analysis

In this section, after introducing parameterizations in the IGA as background, we solve the elliptic boundary value problem [START_REF] Wu | Hermite type spline spaces over rectangular meshes with arbitrary topology accepted by Communications in Computational Physics[END_REF] and compare the results by the IGA with the ones by the FEM with H 3 h (Q h ) as test functions and trial functions. In addition, the FEM is a special case of the IGA method, where the parameterization is an identity map, which is a non-singular parameterization.

In the IGA, the geometric description and numerical analysis share the same basis functions. Thus, in this section, parameterizations are represented by splines in

H 3 h (Q h ), i.e., P(s, t) = (x(s, t), y(s, t)) T = n i=1 P i H i (s, t),
where

{H 1 (s, t), H 2 (s, t), • • • , H n (s, t)} is a set of basis functions of H 3 h (Q h ) and P i ∈ R 2 . Corollary 3.4. If u • P ∈ H 4 (Σ r ), u • P ∈ H k +1 (Σ s ), and J ∼ h σ on singular cells, then ||u -Π Â H 3 h u|| L 2 (Ω) ≤ C 1 h 4 |u • P| H 4 (Σr) + C 2 h k +1+σ/2 |u • P| H k +1 (Σs) ;
In the FEM, if u ∈ H 4 (Ω) the optimal convergence rate is reached. However, if u ∈ H 4-σ/2 (Ω s ) and u ∈ H 4 (Ω r ), then the approximation error order is optimal by using the singular parameterization P(s, t), where Ω s = P(Σ s ) and Ω t = P(Σ t ).

Numerical results

In the following, we present the numerical results in Examples 3.3 and 3.4 and analyze these results in Section 3.3.2. In these examples, the singular parameterization P(s, t) satisfies the integrability assumption 5.1 in [START_REF] Takacs | Existence of Stiffness Matrix Integrals for Singularly Parameterized Domains in Isogeometric Analysis[END_REF] and its construction method follows [START_REF] Wu | H1-parameterizations of planar physical domains with complex topology in isogeometric analysis[END_REF]. By its construction, its Jacobian on a singular cell is J ∼ h 2 , i.e. σ = 2 in Corollary 3.4.

u b (x, y) = y(4 -x 2 )(4 -y 2 ) (x 2 + y 2 ) 1/4 ,
then the numerical convergence rates within the IGA framework with the singular parameterization and FEM are shown in Figures 10(a

a (x, y) = (1 -x 2 )(1 -y 2 )x 3 y 2 ,
then the numerical convergence rates within the IGA framework with the singular parameterization and FEM are shown in Figures 13(a

Conclusions and the Future Work

In this paper, we discuss the convergence rates with singular parameterizations for solving elliptic boundary value problems in the IGA. To discuss all possible weak solutions of elliptic boundary value problems with smooth coefficients, this work is based on a weak solution space F(Ω). The approximation errors on F(Ω) with the L 2 (Ω)-norm and the H 1 (Ω)-seminorm are estimated locally, i.e., these errors reflect the impact of singularities of parameterizations. Based on this approximation error estimates, in Section 3, the convergence rates with the L 2 (Ω)-norm and the H 1 (Ω)-seminorm are discussed with a controllable singular parameterization and singular parameterizations in the IGA. By a controllable singular parameterization, although the optimal convergence rates can not be reached by the FEM, the optimal ones are reached by choosing proper singularities in the IGA framework. Later, we compare the convergence rates by singular H 1 -parameterizations in the IGA and the ones by the FEM. For smooth weak solutions by the FEM and the IGA with the singular parameterizations, both of them obtain the optimal convergence rates; for non-smooth weak solutions by these two methods, the results by the IGA with the singular parameterizations are better than the ones by the FEM. In the future, we will consider the behaviour of singular parameterizations for solving higher order partial differential equations.

Figure 1 :

 1 Figure 1: A parameterization of a geometry with a smooth boundary

Figure 2 :

 2 Figure 2: A physical domain Ω with a corner

Figure 3 :

 3 Figure 3: Regular Meshes Q h , where the bold edges are boundaries.

  where there is a singularity of P on it. Cells not of this type are called regular cells of Σ and are denoted by C r . Σ s = ∪C s , and Σ r = ∪C r . The definition of norms: ∀u ∈ F(Ω) ∩ W m,p (Ω), the norms and seminorms follow the classical definitions in W m,p (Ω), where W m,p (Ω) is the Sobolev space, i.e., for 1 ≤ p < ∞, ||u|| m,p,Ω = ( |α|≤m Ω |∂ α u| p dΩ) 1 p , |u| m,p,Ω = ( |α|=m Ω |∂ α u| p dΩ) 1 p , and for p = ∞, ||u|| m,∞,Ω = max |α|≤m {ess. sup|∂ α u(x, y)| : (x, y) ∈ Ω}; |u| m,∞,Ω = max |α|=m {ess. sup|∂ α u(x, y)| : (x, y) ∈ Ω},

2 Figure 4 :

 24 Figure 4: The parametric domain and physical domain of Example 3.1

Figure 5 :

 5 Figure5: Errors with the L 2 (Ω)-norm and the H 1 (Ω)-seminorm for solving the elliptic boundary value problem[START_REF] Deng | Dimensions of spline spaces over T-meshes[END_REF]

Example 3 . 2 .

 32 (α = 3/2) The elliptic boundary value problem:

2 Figure 6 :

 26 Figure 6: The parametric domain and physical domain of Example 3.2

Figure 7 :

 7 Figure 7: Errors with the L 2 (Ω)-norm and H 1 (Ω)-seminorm in Example 3.2

Example 3 .

 3 3 (α = 2). The parametric mesh is shown in Figure8(a). The parameterization is presented in Figure8(b). In the following, we present the numerical results. (a). If the exact solution is u a (x, y) = (4 -x 2 )(4 -y 2 ) sin y, then the numerical convergence rates within the IGA framework with the singular parameterization and FEM are shown in Figures9(a) and 9(b). (b). If the exact solution is

  ) and 10(b). (c). If the exact solution isu c (x, y) = y(x 2 + y 2 ) 1/4 (4 -x 2 )(4 -y 2 ),then the numerical convergence rates within the IGA framework with the singular parameterization and FEM are shown in Figures11(a) and 11(b).

Figure 8 : 3 0

 83 Figure 8: The parametric domain and physical domain of Example 3.3

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: The errors with the L 2 (Ω)-norm and the H 1 (Ω)-seminorm of Example 3.3 (a)

  u

  ) and 13(b).

Figure 12 :

 12 Figure 12: The parametric domain and physical domain of Example 3.4

Figure 13 :Figure 14 :

 1314 Figure 13: Errors with the L 2 (Ω)-norm and the H 1 (Ω)-seminorm in Example 3.4(a)

Figure 15 :

 15 Figure 15: Errors with the L 2 (Ω)-norm and the H 1 (Ω)-seminorm in Example 3.4(c)

3. 3 . 2 .

 32 Discussion of the numerical results in Examples 3.3 and 3.4: 1. For smooth weak solutions, the FEM and the IGA with singular parameterizations are reached optimal convergence rates. In Examples 3.3(a) and 3.4(a), u a (x, y) ∈ H 4 (Ω), 4 for the L 2 (Ω)-norm and 3 for the H 1 (Ω)-seminorm are reached in Figures 9(a),13(a) and 9(b), 13(b). 2. For non-smooth weak solutions, the convergence rates by IGA with singular parameterizations are better than the ones by the FEM. In Examples 3.3(b) and 3.3(c) (or Examples 3.4(b) and 3.4(c)), u b (x, y), u c (x, y) ∈ F(Ω) but u b (x, y), u c (x, y) are not in H 4 (Ω). Although u b (x, y) and u c (x, y) have different regularities, in Figures 10(a) and 11(a) (or Figures 14(a) and 15(a)), the convergence rates with the L 2 (Ω)-norm by IGA with singular parameterizations are better than the ones by the FEM. There is a similar result with the H 1 (Ω)-seminorm in Figures 10(b) and 11(b) (or Figures 14(b) and 15(b)).


V Dupuis

F Bert

J-P Bouchaud

J Hammann

F Ladieu

D Parker

E Vincent
email: evincent@drecam.saclay.cea.fr

Aging, rejuvenation and memory phenomena in spin glasses

Keywords: Spin glass, relaxation, aging, rejuvenation, memory PACS Nos 75.10.Nr, 75.40.Gb, 75.50.Lk

In this paper, we review several important features of the out-of-equilibrium dynamics of spin glasses. Starting with the simplest experiments, we discuss the scaling laws used to describe the isothermal aging observed in spin glasses after a quench down to the low-temperature phase. We report in particular new results on the sub-aging behaviour of spin glasses. We then discuss the rejuvenation and memory effects observed when a spin glass is submitted to temperature variations during aging, from the point of view of both energy landscape pictures and real-space pictures. We highlight the fact that both approaches point out the necessity of hierarchical processes involved in aging. Finally, we report an investigation of the effect of small temperature variations on aging in spin glass samples with various anisotropies which indicates that this hierarchy depends on the spin anisotropy.

Introduction

Spin glasses are magnetic systems in which the interactions are disordered and frustrated, due to a random dilution of magnetic ions. Below their glass temperature T g , their dynamic response to a small magnetic excitation is slow, and in addition it depends on the time spent below T g ('aging') [START_REF] Vincent | Complex behaviour of glassy systems[END_REF][START_REF] Bouchaud | Spin glasses and random fields[END_REF]. Contrary to structural and polymer glasses, aging in spin glasses is hardly influenced by the cooling rate. On the contrary, each cooling step leads to a restart of dissipation processes ('rejuvenation'), while the memory of previous agings achieved at different temperatures during cooling can be retrieved during re-heating [START_REF] Refregier | [END_REF]4]. It is possible to store independently the memory of several isothermal agings, which corresponds to imprinting the trace of various spin arrangements at well-separated length scales that are selected by temperature [5,6].

In §2 of this paper we first recall some basic features of isothermal aging in spin glasses. We discuss in particular the scaling laws used to describe the response of a spin glass to a field variation and report new results which raise the question of the link between the cooling time and the sub-aging behaviour observed in several spin glasses. In §3, we examine the effect of temperature variations on aging which reveals striking rejuvenation and memory effects and discuss these effects in terms of energy landscape/real-space pictures. Section 4 is devoted to a comparison between an Ising spin glass and several Heisenberg spin glasses with various degrees of random anisotropy, which shows that the aging properties are markedly different in these systems, the memory effect being more pronounced in the Heisenberg case [7,8]. These last results underline the role of anisotropy in the nature of the spinglass phase with Heisenberg spins and are discussed in the context of a cross-over expression for the relationship between time and length scales in an aging spin glass.

Isothermal aging in spin glasses

Basics features and scaling laws

A large variety of systems exhibit slow glassy dynamics and aging effects, i.e. history dependent properties. In all these systems, the response to an external excitation depends on the time t w (waiting time) that the system has spent in the low-temperature phase before the excitation is applied. In most cases, the response function can be typically expressed as the sum of a fast stationary part and a slow aging contribution which roughly scales as t/t w (figure 1a) [START_REF] Vincent | Complex behaviour of glassy systems[END_REF][START_REF] Bouchaud | Spin glasses and random fields[END_REF]:

R(t w + t, t w ) = R s (t) + R a (t/t w ), (1)
where t w is the time elapsed at the measurement temperature T m since the system has been cooled down below its glass temperature T g , and t is the time elapsed since the excitation has been applied. This scaling means that the effective relaxation time of the system below its freezing temperature is mainly set by its age t w + t. It is observed not only in the response to a field variation but also in the ac susceptibility (harmonic response) which obeys an equivalent ωt w scaling, where ω is the frequency of the ac field [START_REF] Vincent | Complex behaviour of glassy systems[END_REF]. The t/t w scaling is however only approximate. It appears that a more precise scaling can be obtained by assuming that the effective relaxation time goes as the full age (t + t w ) to a certain power µ with µ < 1 (figure 1b). This leads to an effective scaling variable ξ = λ/t µ w = 1/(1 -µ)[(t + t w ) 1-µ -t 1-µ w ] (see [START_REF] Vincent | Complex behaviour of glassy systems[END_REF][START_REF] Struik | Physical aging in amorphous polymers and other materials[END_REF]). µ < 1 means that the effective relaxation time grows more slowly than the age of the system, and this situation is thus called sub-aging. Many examples of outof-equilibrium systems such as spin glasses, glassy polymers, microgel pastes and other soft glassy materials exhibit this sub-aging behaviour. In all these cases, the exponent is smaller than 1, but remains close to 1 (of order 0.8-0.9) in a rather large domain of temperatures (0.4 < T /T g < 0.9). Such a behaviour is predicted in simple trap models [START_REF] B Rinn | [END_REF], in some mean-field models of spin glasses [START_REF] Bouchaud | Spin glasses and random fields[END_REF] and in a recent extended droplet model [11,12].

Intrinsic or extrinsic departure from full t/t w scaling

The important question of whether the experimental value of µ smaller than one may be due to experimental artefacts is however far from being settled and has given rise to much interest recently [13]. There are three main experimental limitations which may have a crucial effect on the results: the size limitations of the samples, the application of finite excitations, the long time needed to cool the sample down to the measurement temperature.

In a finite system composed of many subsystems (which could be grains in polycrystalline samples or powder), one can expect that the ergodic times needed to explore the whole configuration space of each subsystem remain also finite. In that case aging will be interrupted at some ergodic time and the exponent µ will decrease and eventually go to zero as the ergodic time is reached in an increasing number of subsystems. This is a natural mechanism for sub-aging which has been used in [14,15] to estimate, from the observed values of µ, the cut-off of the barrier height distribution or that of the size distribution of the considered samples. However, no experimental proof of the effect could be given yet.

The amplitude of the applied magnetic field also is obviously very important. Usually, only the minimum necessary excitations are used in an aging investigation and care is taken to ensure that the response remains in the linear regime. However, it is very difficult to understand in which way and how much the system is perturbed. A systematic study of the field effect was performed in spin glasses [16] and it is found that µ indeed decreases as the excitation field is increased. However, it seems to saturate at very low fields and stays at a value smaller than one as the field is reduced by large factors. In the case of the thiospinel compound, µ saturates around 0.85 below 10 G down to the smallest field probed (10 -3 G) [17].

Recently, some detailed investigations were devoted to the problem of the cooling time and the importance of the initial conditions in the measurement of the thermoremanent magnetization (TRM) in spin glasses [13,18,19]. Different cooling procedures were used and their effects on the TRM were considered. In the quoted investigations, three different cooling protocols were used: (A) direct cooling from the temperature above T g to the measuring temperature T m (usual procedure); (B) fast cooling from the temperature above T g down to a temperature smaller than T m and heating back up to T m ; (C) cooling down to a temperature larger than T m and waiting a certain time before the final cooling step down to T m . Some differences have been found using these various initial conditions. In [13] it was argued that a full aging (µ = 1) could almost be obtained in a Cu:Mn 6% spin glass sample with protocol (B). In [19], the difference between the two protocols (A) and (B) in a Ag:Mn 11% sample is only to reduce the initial age of the system for the shorter cooling time (protocol (B)).

Revisiting the link between cooling time and sub-aging

In order to clarify the situation, we have performed new experiments using protocols (A), (B) and (C) described above, to check the effects of different initial conditions on the scaling behaviour of the response function of three well-characterized spin glasses.

In order to precisely estimate the real cooling time associated with a given protocol, we recorded the field-cooled magnetization, which actually slightly depends on temperature and can thus be used as a thermometer giving the true temperature of the sample. In all cases, TRM measurements made just after the temperature stabilization at T m show a maximum of dM/d(log t) at a time of the order of the real cooling time. In other words, the cooling time in these experiments corresponds to some equivalent initial waiting time t ini at T m . In the scaling procedure, we have corrected t w by the addition of t ini which we use as a free parameter to account for extra aging during cooling.

Figure 2a shows a scaling of several TRM curves measured according to protocol (B) for various waiting times with the Heisenberg-like thiospinel CdCr 1.7 In 0.3 S 4 spin glass. The same sub-aging behaviour as in protocol (A) (figure 1b) can be observed with a slightly decreased value of the exponent µ (µ = 0.83 here, compared to µ = 0.87 in protocol (A)). However, an experiment using protocol (C) showed the same result (µ = 0.87) as with the usual cooling protocol (A), even though the sample was kept for 1 h at T /T g = 0.94 before being cooled down to the measuring temperature T m /T g (figure 2b). Other experiments on a Au:Fe 8% sample, another Heisenberg-like spin glass which is expected to behave similarly as the thiospinel [8], do not show any clear trend when slowing down the cooling procedure. The results are compatible with the thiospinel results, though less accurate due to a smaller signal, and point thus towards the same tiny cooling effects, if any.

We performed the same experiments on Fe 0.5 Mn 0.5 TiO 3 , an Ising spin glass in which stronger cooling effects can be expected [7]. The scaling of the TRM relaxation curves measured with fast (protocol (B)) and slow (protocol (C)) coolings is presented in figure 3. Fast cooling (figure 3a) yields µ = 0.84 and t ini = 75 s, with a good scaling. Slower cooling yields, as can be seen in figures 3b and 3c, (i) a downgrading of the scaling quality, as noted in [13] and (ii) a change in the scaling parameters which can be as well accounted for by an increase of t ini or a decrease of µ. From this analysis, we conclude that there is no evidence that µ would go to 1 for a very fast cooling. Our results are compatible with the idea of an 'intrinsic' nature of the sub-aging behaviour in spin glasses. In numerical simulations of (Ising) spin glasses, in which cooling can be instantaneous, it has been shown [20] that the effect of cooling rate on the value of µ very quickly saturates. µ was found to evolve from a value slightly above 1 for infinite cooling rate to a value less than 1 when the cooling rate is of the order of 10 -2 in microscopic units, and does not evolve significantly when the cooling is slower. Since a cooling rate of 10 -2 is still very much faster than any experimentally achievable cooling rate, this scenario observed in numerical simulations may be a plausible explanation for sub-aging in real materials.

Memory effects and hierarchy of states/length scales

In the previous section, we discussed what happens when one simply quenches a spin glass below its glass transition temperature T g and lets it evolve at constant temperature. Other protocols may be imagined in which one can probe the effects of temperature variations during aging and see, for example, how far heating an aging spin glass slightly can help it to reach its equilibrium [START_REF] Refregier | [END_REF]4].

Figure 4 shows the result of such a protocol on the imaginary part χ of the ac susceptibility measured at 0.1 Hz on the CdCr 1.7 In 0.3 S 4 spin glass. In this experiment, the sample was step-cooled through its glass transition temperature T g = 16.7 K down to 5 K (each plateau at constant temperature lasts approximatively 30 min) and reheated then continuously up to above T g . During each stop at constant temperature, aging is directly evidenced as a downward relaxation of χ . Looking carefully at the figure, we see that further cooling seems to rejuvenate the system as the susceptibility rises back: this is the rejuvenation effect. However, when heated back continuously (full symbols), χ shows marked dips at the plateau temperatures where the sample was allowed to age, revealing that the system has kept the memory of each stop that occurred during the cooling.

Energy landscape pictures: These striking rejuvenation and memory effects can be interpreted in terms of a hierarchical organization of the metastable states in which the system evolves during its search for equilibrium [START_REF] Vincent | Complex behaviour of glassy systems[END_REF][START_REF] Refregier | [END_REF]. This picture, originally inspired from the mean field solution of the SK spin glass model, was further developed in the past years in several toy models such as the 'traps on a tree' model [21] or the Sinai model (in which the hierarchy is put by hand) [22]. On the other hand, recent numerical simulations of the generalized random energy models [23] (roughly equivalent to the 'traps on a tree' model) proved unambiguously that the idea of a hierarchy of metastable states is able to capture the main features of experimentally observed aging.

Real-space pictures: While energy landscape pictures are illuminating, one cannot avoid trying to interpret the dynamics in an aging spin glass in real space, that is in terms of spin configurations. It is now widely accepted that aging in spin glasses corresponds to the growth of a certain coherence length (which would be the equivalent of a kind of out-of-equilibrium correlation length) [5,12]. In a ferromagnet, for example, this coherence length is simply the size of the ferromagnetic domains and the aging observed after a quench just reflects the growth of these ferromagnetic domains. In spin glasses however, and contrary to ferromagnets where the domain growth is directly related to cooling rate (the longer the time spent at higher temperature, the larger the domains), aging is hardly influenced by the cooling rate, as evidenced in figure 4 where each cooling step seems to entirely rejuvenate the system.

A simple real-space scenario which allows one to reconcile the idea of a growing coherence length with the observed rejuvenation and memory effects was discussed [5]. It was noticed that thermally activated dynamics over energy barriers B ∼ ∆(T ) ψ that vanish at T g as ∆(T ) ∼ (T -T g ) ψν (B is the barrier to overcome in order to flip a cluster of spins of typical size ) naturally lead to a hierarchy of the length scales probed in experiments. More precisely, and as described in [5], the temperature T at which aging is observed sets the length scale * (T ) at play in aging (often called active length scale). Larger length scales ( > * (T )) are frozen while smaller length scales ( < * (T )) are fully equilibrated. Decreasing the temperature from T to T -∆T reduces this typical length scale significantly as * (T -∆T ) < * (T ). This has two important consequences. Firstly the structure previously reached at length scale * (T ) freezes since the time needed to reconform details of the structure at this length scale has grown astronomically large. This will allow the frozen structure to be retrieved by heating back the system (memory effect). Secondly, the length scales such that * (T -∆T ) < < * (T ), which were equilibrated at T , are no longer equilibrated at T -∆T due to the change of the associated Boltzmann weights with temperature. They start to age giving rise to the rejuvenation effect. This scenario, which naturally leads to a hierarchy of embedded active length scales, is an appealing transcription of the hierarchical energy landscape pictures. An alternative scenario, based on 'temperature chaos', has been actively promoted in [24]. Although the mechanism leading to rejuvenation is different, this scenario also requires a rapid separation of time-scales as the temperature is reduced to preserve the memory.

Spin glass dynamics, anisotropy and growth of the coherence length

In order to go deeper into the understanding of the out-of-equilibrium dynamics of spin glasses and in particular the rejuvenation and memory effects discussed above, we recently performed a systematic study of the effect of temperature variations on the aging dynamics of various spin glasses with different degrees of anisotropy ranging from Ising to Heisenberg [7,8].

Temperature cycling experiments were performed in a monocrystalline Fe 0.5 Mn 0.5 TiO 3 Ising [25] sample #1 (T g = 20.7 K) and in Heisenberg SGs with decreasing random anisotropy [26]: an amorphous alloy (Fe 0.1 Ni 0.9 ) 75 P 16 B 6 Al 3 #2 (T g = 13.4 K), a diluted magnetic alloy Au:Fe 8% #3 (T g = 23.9 K), an insulating thiospinel CdCr 1.7 In 0.3 S 4 #4 (T g = 16.7 K) and the canonical Ag:Mn 2.7% SG #5 (T g = 10.4 K) (data from [27]).

In these experiments, an isothermal TRM procedure allows one to first define a set of reference curves that depend on the waiting time (such as the one in figure 1), and to which we will compare the relaxation curves obtained after more complicated histories. In the temperature cycle experiments, the sample is cooled to a temperature T < T g , kept at T for a time t 1 = 500 s, submitted to long wait at T -∆T during t 2 = 9000 s and finally re-heated to T for another short time t 3 = t 1 . The magnetic field is then switched off and the TRM is recorded. For each 'cycled' TRM curve, an effective waiting time t eff 2 is defined so as to superimpose the obtained TRM curve with that of the purely isothermal TRM recorded at temperature T , with waiting time t w = t 1 + t eff 2 (T, ∆T ) + t 3 . For ∆T = 0, the effective waiting time equals the actual waiting time, but becomes smaller for ∆T > 0: t eff 2 (T, ∆T ) < t 2 characterizes the effect at T of aging at T -∆T . The results for the different SGs and for two different temperatures, T /T g 0.72 and T /T g 0.85, are plotted in figure 5 (for sample #5, a similar procedure, with t 1 = t 3 = 0, was applied [27]).

The results of figure 5 which correspond to small ∆T allow us to probe the separation of length scale discussed in the previous section [5]. Indeed, when considering two temperatures T and T -δT such that * (T -δT ) ∼ * (T ), rejuvenation and memory effects are no longer observed and aging tends to accumulate from one temperature to another, which we described here in terms of an effective time t eff 2 . The stronger the decrease of this effective time with increasing δT , the sharper the separation of active length scales. Figure 5 thus shows that the separation of length scales becomes sharper and sharper when going from the Ising sample towards Heisenberg spin glasses with weaker and weaker anisotropy. Our results single out the role of the spin anisotropy in the nature of the Heisenberg spin glass phase.

A purely activated scenario over constant height barriers would correspond to a single straight line for all samples in figure 5 [8]. In order to account for the different observed behaviours, we introduced a cross-over expression for the relationship between time and length scales

t w = τ 0 z exp(∆(T ) ψ /T ). (2)
In this expression, which is simply an Arrhenius law with a temperaturedependent barrier, the usual attempt time τ 0 is replaced by a renormalized attempt time τ 0 z where z is the dynamic critical exponent. This modification has the advantage of giving the expected form for the critical dynamics close to T g , where due to the vanishing of the barrier height ∆(T ) = Υ 0 T g (1 -T /T g ) ψν the exponential term in eq. ( 2) tends to 1. In [8] the above expression was successfully used to explain consistently both the temperature cycling experiments and field variation experiments. This in turn gave, for each sample, estimates of the barrier exponents ψ and the energy scales Υ 0 of the barriers which are consistent with a sharper separation of the active length scales in Heisenberg-like spin glasses than in the Ising case (see [8] for details). The same trend has recently been observed in large-scale simulations of Ising and Heisenberg spin glasses [28].

Finally, let us mention an issue that concerns the value of the individual spin flip time τ 0 (arising, e.g., in eq. ( 2)), and that is usually assumed to be microscopic, of the order of 10 -12 s. This assumption is very important since it means, as mentioned above, that typical experimental cooling rates are extremely slow when compared with τ 0 . However, it might well be that the flip of a single spin is already slowed down by its environment (interactions, local anisotropy, etc.). In the Ising case, this slowing down might be somewhat large, and could significantly renormalize the value of τ 0 , bringing it much closer (especially at low temperatures) to the experimental time window. This could significantly affect the precise values of the parameters involved in eq. (2).

Conclusion

In this paper, we have briefly reviewed several important features of the slow out-ofequilibrium dynamics of spin glasses. We have first discussed the scaling laws that govern the isothermal aging of the response function of spin glasses, and showed experiments which suggest that the observed sub-aging behaviour may be intrinsic, although when compared to a similar effect observed in numerical simulations, laboratory cooling rates might still be too slow to conclude. We have then discussed the interpretation of the rejuvenation and memory effects observed in spin glasses submitted to temperature variations during their aging. In both energy landscape and real-space pictures, the need for a hierarchy of state/length scales seems unavoidable for the description of such effects. This strong separation of length scales probed at different temperatures is also needed in the 'temperature chaos' interpretation of these rejuvenation/memory effects (see [24] for a review and detailed arguments supporting this point of view, and [START_REF] Jimenez | Rejuvenation and memory in model spin glasses in 3 and 4 dimensions[END_REF] for very recent numerical simulations challenging this interpretation).

We have reported experiments probing the effect of very small temperature variations in spin glass samples with various degrees of anisotropy. The main results of this last investigation are: (a) a clear finger-print of a dynamical cross-over between a critical regime and an activated regime in spin glasses, with energy barriers vanishing at T c and (b) a trend to a sharper separation of active length scales with temperature in Heisenberg-like samples.

Figure 1 .

 1 Figure 1. (a) Thermoremanent magnetization relaxations at 0.7Tg on the thiospinel CdCr 1.7 In 0.3 S 4 spin glass (T g = 16.7 K); (b) scaling of the previous relaxation curves (see text for details).

Figure 2 .

 2 Figure 2. Scaling of thermoremanent magnetization relaxations measured on the CdCr 1.7 In 0.3 S 4 spin glass at 0.7T g . (a) Fast cooling (protocol (B)) (for this set of data, the exact vertical scale is unknown, but within the error bars A corresponds to A = 0.05); (b) slow cooling (protocol (C)) (see text for details).

Figure 3 .

 3 Figure 3. Scaling of thermoremanent magnetization relaxations measured on the Fe0.5Mn0.5TiO3 Ising spin glass at 0.7Tg: (a) fast cooling (protocol (B)); (b) slow cooling (protocol (C)), scaling with µ fixed to 0.84 and variable t ini ; (c) slow cooling (protocol (C), same data), but with t ini fixed to 70 s and variable µ.

Figure 4 .

 4 Figure 4. Series of dips imprinted on the ac susceptibility by successive stops at different temperatures while the system is cooled.

Figure 5 .

 5 Figure 5. Effective waiting times deduced from the temperature cycle experiments at (a) T 0.72T g and (b) T 0.85T g . The spin glasses are labelled from #1 to #5 by order of decreasing anisotropy.

Pramana -J. Phys., Vol. 64, No. 6, June 2005

Pramana -J. Phys., Vol. 64, No. 6, June 2005

Acknowledgements

This work was supported in part by the European Community's Human Potential Programme (contract HPRN-CT-2002-00307 referred to as DYGLAGEMEM).


Carlene Paula

Silva De Farias

Christopher Thomas

Catherine Pépin

Alvaro Ferraz

Claudine Lacroix

Sébastien Burdin

Carlene Farias

Spin liquid versus long range magnetic order in the frustrated body-centered tetragonal lattice

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

In most superconducting materials, magnetic ordering destroys superconductivity. However, magnetic correlations are essential for superconducting pairing in a large variety of unconventional superconductors. In the presence of frustration, classical magnetic order can become unstable against spin-liquid (SL) states, and superconductivity may emerge even with strong magnetic interactions. It is commonly believed that low dimensionality is essential for the formation of a SL, and as result this would also be a crucial ingredient for the realization of superconductors with higher critical temperature. Controverting this idea, we analyze the effect of frustration in the body-centered tetragonal (BCT) lattice, which is one of the 14 existing three-dimensional (3D) lattice structures [START_REF] Kittel | Introduction to Solid State Physics[END_REF]. This standard crystalline structure is realized in several strongly correlated electron materials with unusual magnetic and transport properties. Among the heavy fermion systems [START_REF] Stewart | [END_REF]3], different examples of materials with rare-earth atoms on a BCT lattice have been intensively studied for the last few decades: in URu 2 Si 2 , a still-mysterious hidden order (HO) phase appears below the critical temperature T H O ≈ 17 K close to a pressure-induced antiferromagnetic (AF) transition [4,5]; in YbRh 2 Si 2 and CeRu 2 Si 2 , non-Fermi-liquid properties are observed in the vicinity of the AF quantum phase transitions, which are still poorly understood [6][7][8][9]; CeCu 2 Si 2 was the first heavy-fermion material where unconventional superconductivity was discovered close to an AF transition [10]; CePd 2 Si 2 also exhibits unconventional superconductivity related to an AF transition [11,12]; multi-Q AF order has been observed in CeRh 2 Si 2 [13]. It is noticeable that the cuprate superconductors [14] can be included among the AF insulating parent compounds La 2 CuO 4 and Sr 2 CuO 2 Cl 2 in which the AF order originates from the Cu atoms that form a BCT crystal. However the relevant physics in their case is essentially two dimensional (2D) with the BCT structure being involved only in the formation of the square-lattice layers of Cu atoms that order antiferromagnetically. In this paper we suggest that the rich diversity of unusual physical properties observed in 3D materials with BCT structure is essentially associated with the underlying frustration. To do that we establish direct contact with important theoretical developments which were made in the past few years on the unconventional magnetic properties of the BCT lattice, following the pioneering study by Villain [15][16][17][18][19][20][21][22].

In this article, we analyze the ground states of an SU(n) generalization of the J 1 -J 2 -J 3 quantum spin S Heisenberg Hamiltonian H S n .

II. CLASSICAL SPIN AND QUANTUM SPIN-WAVE-FLUCTUATIONS APPROACH

The model Hamiltonian is introduced here, for simplicity, for n = 2:

H S n=2 = R,R J RR S R • S R , ( 1
)
where S R ≡ (S x R ,S y R ,S z R ) denotes the quantum spin S operators acting on site R of a BCT-lattice. The antiferromagnetic interaction J RR connects sites R and R , and can take three possible values J 1 ,J 2 ,J 3 > 0, as indicated in Fig. 1.

A. Classical ground state

Following the standard spin-wave (SW) approach, we start with the S = +∞ generalization of the Hamiltonian (1), which corresponds to the classical spin limit. For the sake of simplification, we consider only magnetic orderings that are characterized by a single wave vector Q = 2π (λ,μ,ν) identified as (λ,μ,ν) in reduced notation. Invoking the Fourier transform, this wave vector is used to minimize the classical dispersion,

J (q) ≡ 8J 1 γ q 1 + 2J 2 γ q 2 + 4J 3 γ q 3 , (2)
with

γ q 1 ≡ cos (q x /2) cos (q y /2) cos (q z /2), (3)
γ q 2 ≡ cos (q x ) + cos (q y ), (4)
γ q 3 ≡ cos (q x ) cos (q y ).

(

) 5
By tuning the dimensionless parameters p 2 ≡ J 2 /J 1 and p 3 ≡ J 3 /J 1 , we find that the ground state can be characterized by four kinds of possible wave vectors, as depicted in Fig. 2: Q I AF ≡ (1,1,1) and Q II AF ≡ (1/2,1/2,ν) corresponding to the regimes where the Weiss field is dominated by J 1 and J 2 , respectively. The ν degeneracy in the latter case indicates the underlying bidimensionality. The other possible ordering wave vectors are incommensurate and characterize two kinds of helical orders: Q III inc ≡ (0,ϒ 3 ,0) degenerate with (ϒ 3 ,0,0), (1,ϒ 3 ,1), and (ϒ 3 ,1,1) where

ϒ 3 = 1 π arccos -1 p 2 +2p 3 ; and Q II inc ≡ (ϒ 2 ,±ϒ 2 ,1) degenerate with (ϒ 2 ,1 ± ϒ 2 ,0), where ϒ 2 = 1
2π arccos 1-p 2 2p 3 . A different wave vector, Q III AF ≡ (0,1/2,ν) had been proposed [23] in a J 3 -dominated phase, which corresponds to the commensurate order characterizing a purely bidimensional square lattice. We find that J 1 = 0 corrections are relevant and Q III inc is energetically more stable than Q III AF . Not surprisingly these two vectors are asymptotically identical at large p 3 . Similarly,

Q II inc → Q II AF at large p 2 .
The transition lines separating these four classical phases are given by linear relations between p 2 and p 3 . The Q II inc -Q III inc transition is discontinuous, the other transitions are continuous.

B. Spin-wave fluctuations

Next, assuming a given classical ground-state wave vector Q, we investigate the large-S corrections. This expansion invokes a helical generalization of the Holstein-Primakov representation [16,19,24]: introducing boson annihilation (creation) a ( †) R operators, the spin operators are approximated as

S x R ≈ (S/2) 1 2 (a R + a † R ), (6)
and

S z R S y R ≈ cos θ R sin θ R -sin θ R cos θ R S -a † R a R -i(S/2) 1 2 (a R -a † R ) , (7) with θ R ≡ Q • (R -R 0 )
, where z is the easy axis characterizing a site R 0 chosen arbitrarily. The dispersions obtained for the Bogoliubov quasiparticles are

+ k ≡ J (k) -J (Q), (8)
- k ≡ J (k + Q) + J (k -Q) 2 -J (Q). ( 9
)
Order by quantum disorder. While the ground-state energy is proportional to S 2 J (Q) at the highest order, the first correction is proportional to

S BCT d 3 k √ + k - k
where the k integral runs over the first Brillouin zone of the BCT-lattice. Analyzing this correction for the classically degenerate vectors Q II AF , we find that the resulting order is stabilized by quantum disorder: the continuous degeneracy is lifted in favor of ν = 0, which is equivalent to ν = 1.

Fluctuation corrections to magnetization. We also studied the effects of fluctuations emerging from the linear SW corrections. Generalizing to the BCT structure the approach introduced in Ref. [25] for the square lattice model, the staggered magnetization is expanded around its classical value S z R 0 ≈ Sm(p 2 ,p 3 ). We find

m(p 2 ,p 3 ) = a † R 0 a R 0 = - 1 2 + BCT d 3 k 64π 3 + k + - k + k - k . ( 10
)
Unlike the 2D case [25], the fluctuation corrections do not diverge, which is not surprising for a 3D model. Frustration can relatively increase the critical value of S below which the linear spin-wave correction cancels the staggered magnetization, S c ≡ m(p 2 ,p 3 ). Indeed, for a fixed p 2 , we find that S c increases when p 3 approaches its critical value associated with the classical phase boundary. This maximal value is plotted in the inset of Fig. 2 as a function of p 2 for the continuous Q II AF /Q II inc transition and on each side of the discontinuous Q II inc /Q III inc transition. On each of these critical lines, we find S c ∼ √ p 2 at large p 2 . Furthermore, a logarithmic S c ∼ ln(p 2 /2p 3 ) and a power law S c ∼ 1/ √ p 3p 2 /2 are respectively obtained at large p 2 in the vicinity of the

Q II inc → Q III inc and Q III inc → Q II inc transitions.
This result is consistent with the square lattice spin wave analysis [25].

C. Discussion

The SW approach thus reveals some weaknesses of the classical magnetic orders, but the three-dimensionality protects these states against small fluctuations at the lowest, linear order. One may go further and study possible instabilities emerging from the next orders in the 1/S expansion, taking into account interactions between the spin-wave bosonic excitations. However, it has been shown for a two-dimensional lattice that the second-order correction can increase considerably the sublattice magnetization, although the first-order perturbation makes it decrease [26]. It is thus possible that a sign-oscillating series could also characterize the SW expansion from the orders characterizing our model. Therefore, stabilizing a spin-liquid state might require going beyond second order when invoking a perturbative approach from classical order. Hereafter, we bypass this open issue and we follow an alternative and simpler approach.

III. GENERALIZED SU(n) SYMMETRIC APPROACH

In this section, we analyze the possibility that the system forms a resonant valence bond state with fermionic excitations and SL correlations [27,28]. One of our physical motivations is driven by the physics of unconventional metallic systems with BCT structure: in several of these correlated systems, magnetic degrees of freedom seem to be "deconfined" into fermionic ones that may contribute to the formation of a Fermi surface, unlike weakly coupled bosons. Such a scenario, inspired by the physics of cuprate superconductors [29][30][31][32], could be easily strengthened by a coupling of the Heisenberg spins of the J 1 -J 2 -J 3 model to extra charge degrees of freedom. In the following we study this possible "fermionic deconfinement" of spin operators as an intrinsic property of the Heisenberg model. To achieve this goal, considering S = 1/2, the Hamiltonian (1) is generalized to its SU(n)-symmetric form:

H S=1/2 n = R,R J RR n σ σ χ † Rσ χ Rσ χ † R σ χ R σ , (11)
where

χ ( †)
Rσ (χ Rσ ) are annihilation (creation) fermionic operators with orbital degeneracy σ = 1, . . . ,n and satisfying the local constraints σ χ † Rσ χ Rσ = n/2. This is a standard SU(n) generalization [33] of the fermionic representation developed by Abrikosov for n = 2. The scaling factor 1/n ensures that the energy remains extensive, i.e., proportional to n, in the large-n limit.

A. Spin-liquid correlations

Using the Hubbard-Stratonovich decoupling as described in Refs. [34,35], the low-temperature phases of Hamiltonian (11) can be characterized by two kinds of order parameters: the local magnetization field

m σ R = χ † Rσ χ Rσ - 1 2 , ( 12
)
and the intersite spin-liquid fields

ϕ RR = - 1 n σ χ † Rσ χ R σ . ( 13
)
The purely magnetic classical mean-field theory characterized here by a staggered magnetization m σ R = ±S Q is equivalent to the one we analyzed for the S = +∞ limit. The corresponding ground-state phase diagram is thus given by Fig. 2. Hereafter, the stability of these classical magnetic orders is analyzed by testing various SL Ansätze as possible alternative ground states. Generalizing the modulated SL (MSL) order introduced elsewhere [35], we consider the nearest-neighbor intersite correlations:

ϕ 1 RR = 1 2 [ 1 + ie iQ I AF •( R+R 2 ) M ], (14)
ϕ 2 RR = 2 , ( 15
)
ϕ 3 RR = 3 , ( 16
)
with a bond index definition similar to that of Fig. 1. The free energy per spin component and per lattice site is expressed as

F = F 0 - k B T 32π 3 BCT d 3 k s=± ln (1 + e - E s k k B T ) - λ 0 2 , ( 17
)
where λ 0 denotes a Lagrange multiplier that minimizes F in order to satisfy the constraint for the fermionic occupation.

For the AF orderings,

F 0 = -2J (Q) n |S Q | 2 and the dispersion is E ± k = λ 0 ± J (Q) n S Q .
For the SL states we find

F 0 /J 1 = | 1 | 2 + | M | 2 + 2p 2 | 2 | 2 + 2p 3 | 3 | 2 and (E ± k -λ 0 )/J 1 = 2p 2 γ k 2 2 + 4p 3 γ k 3 3 ± 4 γ k 1 1 2 + γ k M M 2 .
Here, the non-BCT-periodic real term M with γ k M ≡ sin (q x /2) sin (q y /2) sin (q z /2) takes into account a possible spatial amplitude modulation of the SL field. Our numerical analysis is based on self-consistent calculations of the meanfield parameters for each AF or SL state. The most stable solution is then selected by free-energy minimization, invoking Eq. (17). Considering that the MSL is a BCT adaptation of the "kite" phase investigated in Ref. [33] for a square lattice, we also tested another nonhomogeneous "flux" phase SL characterizing a chiral state with complex

ϕ 1 RR = (ϕ 1 R R ) (a) (b) (c) p 3 0 0.5 1 1.5 0 0 .5 1 1 .5 2 Q I AF Q II AF Q III inc Φ 2 Q II inc p 3 0 0.5 1 1.5 0 0 .5 1 1 .5 2 Q I AF Q II AF Q III inc Φ 2 MSL Q II inc Q II inc p 3 p 2 0 0.5 1 1.5 0 0 .5 1 1 .5 2 Q I AF Q III inc Φ 2 MSL Φ 3 FIG. 3. Ground-state phase diagram of the J 1 -J 2 -J 3 model in coordinates (p 2 ,p 3 ) computed self-consistently for (a) n = 4, (b) n = 5, (c) n = 7.
that could be described within a very close formalism by simply replacing γ k M → sin (q x /2) sin (q y /2) cos (q z /2). This chiral SL was found to have a higher energy than the MSL.

B. Phase diagram

For n 3 we find purely AF ground states, and for n 10 the most stable states are SL characterized by finite values of either M , 2 , or 3 . At the mean-field level, these three SL parameters do not coexist. Furthermore, we remark that the transition between the MSL and the 2 -dominated SL phases is first order. Beyond mean field, we expect that only the SL critical temperature associated with a nonzero M still corresponds to a phase transition signaled by the translationsymmetry breaking. For 4 n 9 we find a rich phase diagram exhibiting AF and SL quantum phase transitions that are controlled by J 1 -J 2 -J 3 parameters, as illustrated in Fig. 3. Increasing n, the AF-SL instability shows up first within the Q II inc phase. Furthermore, comparing the large-S phase diagram (Fig. 2) with that obtained for n = 4 (Fig. 3), one observes that the S c = 1/2 spin-wave instability is located in the same region where the 2 -dominant state becomes stabilized. Also, beyond the specificity of the associated order parameters, Fig. 3 indicates that the SL instability "propagates" from large-p 2 ( 2 ) to smaller-p 2 (MSL) areas if we increase the value of n.

C. Modulated spin-liquid phase

An interesting feature also appears for the MSL solution: with a relatively high numerical accuracy the modulation field M is found to be always equal to the homogeneous field 1 . This leads to a very extreme situation for the interlayer field ϕ 1 RR = 1 2 [ 1 ± M ] which vanishes on half of the bonds while it keeps the finite value 1 = M on the other bonds. Introducing the probability p singlet RR that a given bond RR forms a singlet, the formation of the MSL state can be interpreted here as follows: First, for all the interlayer bonds such that

Q I AF • (R + R )/2 = π/2
, the interaction terms are effectively decoupled at mean-field level, leading to p singlet RR = 1/4 and S R • S R = 0. Then the SL with S R • S R = 0 is formed on the other interlayer bonds, with Q I AF • (R + R )/2 = -π/2, that remain effectively coupled. Using the numerical value 1 = M ≈ 0.45 computed at T = 0 in the MSL, we find that p singlet RR ≈ 0.60 on these effectively coupled bonds. This value has to be compared with the value ln(2) ≈ 0.69 that is predicted for a one-dimensional Heisenberg chain by using exact methods [36,37]. We may thus interpret the MSL as a crystal of interacting filaments formed by the connected effectively coupled bonds. In this picture, spin excitations are deconfined fermions moving along these filaments. This may generalize the usual concept of valence-bond crystal [START_REF]Introduction to Frustrated Magnetism[END_REF] where localized spin-1 excitations correspond to confined fermions.

IV. DISCUSSION

Here we considered a model with only localized spins. However, we know from previous works on heavy fermions and cuprates that charge fluctuations play a crucial role in destabilizing AF states. In the context of the cuprates, the AF phase of the insulating parent compounds corresponds to Q II AF . The SL phase introduced by Anderson et al. [28][29][30][31][32] corresponds to the homogeneous correlated state associated here with 2 = 0. One crucial specificity of this SL scenario for superconductivity in cuprates relies on the twodimensionality of the system. Stabilizing a SL state in 3D is commonly thought to be more tricky in view of the fact that the corresponding linear SW correction remains finite within a large-S approach. Nevertheless, we have shown how frustration in the BCT lattice can enhance the critical value S c that, in some sense, characterizes the weakness of AF order against SW fluctuations. In connection with this weakening, we identified various SL phases that can be stabilized when n is larger than a relatively small critical value. The possible emergence of stable SL-like phases in cubic lattices has also been proposed recently by invoking different approximations [START_REF] Laubach | [END_REF]40]. This opens new perspectives for the realization of unconventional electronic quantum orderings in 3D. In particular, our results suggest that the BCT-lattice structure can play a central role in crystalline materials such as the 122 and some cuprates in spite of numerous phases emerging in these systems. Indeed, the very rich phase diagram depicted in Figs. 2 and3 could provide a unifying framework for understanding and analyzing the intersite correlations in these compounds. The number n may be considered as an effective parameter related to the electronic orbital degeneracy, which could be phenomenologically increased or decreased by charge fluctuations or crystal-field effects. Considering a given compound, n might also be effectively decreased by applying an external magnetic field. Similarly, a tuning of the model parameters p 2 and p 3 may phenomenologically account for some effects of applied pressure [35]. For example, this scenario could explain two different AF instabilities of the HO phase that are observed experimentally in URu 2 Si 2 : assuming that the HO is a MSL order and applying a finite pressure this corresponds to increasing J 1 for a fixed n (see Fig. 3), leading to a commensurate AF instability characterized by Q I AF , as observed experimentally [4,5]. Alternatively, applying a magnetic field without pressure corresponds to lowering the effective value of n for fixed p 2 and p 3 : the MSL (HO) is destabilized to an incommensurate AF. Interestingly, by using different numerical values of J 1 , J 2 , and J 3 obtained from different fits of inelastic-neutron-scattering data, our scenario predicts an instability from MSL to Q III inc = (ϒ 3 ,0,0) with ϒ 3 ≈ 0.69 (from Ref. [41]), 0.66 (from Ref. [42]), 0.69 (from Ref. [23]), and 0.65 (from Ref. [43]). This scenario could be tested experimentally since it predicts that the AF order Q III inc could be continuously tuned to Q I AF by applying pressure on URu 2 Si 2 under a high magnetic field.

Invoking the SL instabilities described in this paper we may also generalize to 3D systems the spin-fluctuation pairing mechanism that was proposed in terms of gauge transformations in Refs. [44,[START_REF] Wen | Quantum Field Theory of Many-Body Systems: From the Origin of Sound to an Origin of Light and Electrons[END_REF]. Here, the link between the BCT lattice structure and the superconducting order parameter is natural. It can be tested experimentally since we predict that the symmetries of the resulting superconducting order parameters will result from the point-group symmetries of the SL, which may correspond to an in-plane pairing related to 2 or 3 or a fully 3D pairing associated with 1 . This SL mechanism driven by frustration on the BCT lattice may also be tested for the heavy-fermion superconductors CeRu 2 Si 2 and CePd 2 Si 2 , but in these systems valence-fluctuation effects need to be carefully included. The possible formation of a MSL could also give rise to a commensurately ordered pairing that would break the BCT symmetry down to simple tetragonal. Such a modulated pairing unconventional scenario could be tested with the superconducting instability observed in URu 2 Si 2 inside the HO phase. Alternatively, even if the chiral SL order was found here to be less stable than the MSL, an opposite result could occur by including charge fluctuations.

FIG. 1 .

 1 FIG. 1. Left: BCT lattice and the J 1 , J 2 , and J 3 interactions. In this work, the tetragonal lattice constants are set to a = b = c = 1. Right: Bold lines represent the three kinds of intersite SL correlations on the BCT structure.

2 FIG. 2 .

 22 FIG. 2. Classical ground-state phase diagram of the J 1 -J 2 -J 3 model in coordinates (p 2 ,p 3 ). The solid lines indicate transitions between the stable magnetic ordered states with modulating vectors Q I AF , Q II AF , Q II inc , and Q III inc . The SW instability region is computed for S = 1/2. Inset: Critical value S c as a function of p 2 when p 3 is bound to the critical lines. The three symbols indicate the corresponding location of critical value S c = 1/2 in each case. Numerical results obtained by using multidimensional Monte Carlo k integration over the BCT Brillouin zone.

ACKNOWLEDGMENTS

We acknowledge the financial support of Capes-Cofecub Ph 743-12. C.T. is bolsista Capes. This research was also supported in part by the Brazilian Ministry of Science, Technology and Innovation (MCTI) and the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). Research carried out with the aid of the Computer System of High Performance of the International Institute of Physics-UFRN, Natal, Brazil. The authors are grateful to Frédéric Bourdarot for useful discussions.