Valérie Labas

Ana-Paula Teixeira-Gomes

Frédéric Andersson

Sébastien Ménigot

Martine Batailler

Hans Adriaensen

Martine Migaud

Elodie Chaillou

Representations Combining in vivo 3D MRI and ex vivo 2D MSI Approaches to Study the Lipid Distribution in the Whole Sheep Brain

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


Alexandre Ern
email: alexandre.ern@enpc.fr

Iain Smears
email: iain.smears@inria.fr

Martin Vohralík
email: martin.vohralik@inria.fr

Discrete p-robust H(div)-liftings and a posteriori estimates for elliptic problems with H -1 source terms *

We establish the existence of liftings into discrete subspaces of H(div) of piecewise polynomial data on locally refined simplicial partitions of polygonal/polyhedral domains. Our liftings are robust with respect to the polynomial degree. This result has important applications in the a posteriori error analysis of parabolic problems, where it permits the removal of so-called transition conditions that link two consecutive meshes. It can also be used in a the posteriori error analysis of elliptic problems, where it allows the treatment of meshes with arbitrary numbers of hanging nodes between elements. We present a constructive proof based on the a posteriori error analysis of an auxiliary elliptic problem with H -1 source terms, thereby yielding results of independent interest. In particular, for such problems, we obtain guaranteed upper bounds on the error along with polynomial-degree robust local efficiency of the estimators.

Introduction and main results

We study in this paper two different but connected problems that we introduce separately.

Discrete p-robust H(div)-liftings

First, we are interested in the problem of finding liftings of piecewise polynomial data into piecewise polynomial subspaces of H(div), that are robust with respect to the polynomial degree (p-robust). More precisely, let Ω ⊂ R d , d ∈ {2, 3}, be a bounded Lipschitz polygonal/polyhedral connected open set. We consider L 2 (Ω) := L 2 (Ω; R d ), and H(div, Ω) := {v ∈ L 2 (Ω), ∇•v ∈ L 2 (Ω)}. Consider a partition of the boundary Γ of Ω into two connected components Γ D and Γ N . Let T be a given conforming, simplicial, possibly locally refined mesh of Ω. We assume that T matches Γ D and Γ N in the sense that every boundary face of the mesh T is fully contained either in Γ D or in Γ N . Let f be a given scalar function and let ξ be a given vector field, such that f and ξ are piecewise-polynomials with respect to T . We consider the question of finding a piecewise polynomial vector field σ h in the Raviart-Thomas-Nédélec subspace of H(div, Ω) over the mesh T , such that ∇•σ h = f in Ω, σ h •n = 0 on Γ N , and such that σ h + ξ is quasi-minimal 1 in the sense of satisfying a bound of the form

σ h + ξ min v∈H(div,Ω) ∇•v=f in Ω v•n=0 on ΓN v + ξ . (1.1)
The notation a b means that a ≤ Cb, with a constant C that can only depend on the shaperegularity of T and on the space dimension d, but is otherwise independent of the domain Ω, of the size of the mesh elements in T , and crucially of the polynomial degree. Note that the converse bound in (1.1) holds trivially with constant 1, since σ h is a member of the minimization set considered in the right-hand side. Problem (1.1) is also known as the problem of finding a stable right-inverse of the divergence operator, and it plays an important role in numerical analysis. We remark that many classical approaches based on projection operators that possess commuting diagram properties are not sufficient for p-robustness. Recently, building on [START_REF] Costabel | On Bogovskiȋ and regularized Poincaré integral operators for de Rham complexes on Lipschitz domains[END_REF][START_REF] Demkowicz | Polynomial extension operators. Part III[END_REF], Braess, Pillwein, and Schöberl [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF] showed the existence of discrete p-robust H(div)-liftings of piecewise-polynomial data on a patch of triangular elements sharing a common vertex, see [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF]Thm. 7]. The extension to three space dimensions is given in [START_REF]Stable broken H 1 and H(div) polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions[END_REF]Thm. 2.2]. These results imply that equilibrated flux a posteriori error estimators for elliptic problems are p-robust. In the present context, the results of [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF][START_REF]Stable broken H 1 and H(div) polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions[END_REF] establish the existence of a discrete p-robust H(div)-lifting satisfying (1.1) when T is a set of elements sharing a common vertex, with Ω being the patch composed of these elements. The existence of H(div)-liftings on locally refined meshes (and not just on element patches around mesh vertices) has been recently studied in [START_REF] Ainsworth | Discrete extension operators for mixed finite element spaces on locally refined meshes[END_REF], where the authors considered vanishing interior source terms but nonzero boundary data and studied liftings for a fixed polynomial degree, with constants typically depending on it.

We now present the first main contribution of this work on problem (1.1). Let

H 1 * (Ω) := H 1 ΓD (Ω) be the subspace of functions in H 1 (Ω) with vanishing trace on Γ D if Γ D is nontrivial, and otherwise (that is, if Γ N = ∂Ω), let H 1 * (Ω) := H 1 (Ω)
/R be the space of functions in H 1 (Ω) with mean-value zero. For an integer p ≥ 0, let RT N p (T ) ⊂ L 2 (Ω) denote the space of piecewise Raviart-Thomas-Nédélec vector fields of order p with respect to T ; note that in the present notation, we do not impose H(div, Ω)-conformity on RT N p (T ). For further details on the notation, see Section 2 below. Our first main result on discrete p-robust H(div)-liftings is the following theorem.

Theorem 1.1. Let p ≥ 1. For any f ∈ P p-1 (T ) and ξ ∈ RT N p-1 (T ), satisfying (f, 1) = 0 if Γ N = ∂Ω, we have min v h ∈H(div,Ω)∩RT Np(T ) ∇•v h =f in Ω v h •n=0 on ΓN v h + ξ min v∈H(div,Ω) ∇•v=f in Ω v•n=0 on ΓN v + ξ = max v∈H 1 * (Ω)\{0} (f, v) -(ξ, ∇v) ∇v . (1.2)
We emphasize that the constant in the stability bound (1.2) does not depend on the polynomial degree p. The last equality in (1.2) follows from classical equivalence results between primal and dual mixed formulations of elliptic problems. This identity is actually important in the proof of Theorem 1.1. The proof, as detailed below, is constructive and consists of the following two steps: first, we pose a primal problem using the data f and ξ and approximate it using H 1 -conforming finite elements of degree p = 1. Then we use this approximate solution to build equilibrated fluxes around each vertex of T by posing local minimization problems using Raviart-Thomas-Nédélec spaces of order p, and we use the discrete p-robust H(div)-liftings on patches that were described above. Theorem 1.1 requires that the discrete H(div)-conforming lifting be of one polynomial degree higher than the piecewise-polynomial data f and ξ. At present, we do not know if equal-order H(div)-liftings retain p-robustness for general data. Nevertheless, the next theorem shows that equal-order p-robust H(div)-liftings are possible when the data take a certain specialised form, and when the norm on the right-hand side of the stability bound is slightly strengthened. Let h Ω stand for the diameter of Ω and let • ∞ denote the L ∞ -norm.

Theorem 1.2. Let p ≥ 1. Let ψ † ∈ H 1 (Ω)∩P 1 (T ) be a continuous piecewise affine function with respect to T . Let F † denote the set of boundary faces F of the mesh T such that ψ † | F = 0, and let Γ N, † = ∪ F ∈F † F . Then, for any f ∈ P p-1 (T ) and ξ

∈ RT N p-1 (T ), with (f, ψ † ) = (ξ, ∇ψ † ) if Γ N, † = ∂Ω, we have min v h ∈H(div,Ω)∩RT Np(T ) ∇•v h =ψ † f -∇ψ † •ξ in Ω v h •n=0 on Γ N, † v h + ψ † ξ C(Ω, ψ † ) min v∈H(div,Ω) ∇•v=f in Ω v + ξ (1.3a) = C(Ω, ψ † ) max v∈H 1 0 (Ω)\{0} (f, v) -(ξ, ∇v) ∇v , (1.3b)
where

C(Ω, ψ † ) = ψ † ∞ + C P, † h Ω ∇ψ † ∞ and C P, † is the Poincaré constant of the space H 1 † (Ω) := H 1 (Ω)/R if Γ N, † = ∂Ω, and H 1 † (Ω) := H 1 ∂Ω\Γ N, † (Ω) otherwise, i.e. v ≤ C P, † h Ω ∇v for all v ∈ H 1 † (Ω).
Note that Theorem 1.2 is indeed optimal with respect to the polynomial degrees of the data, as ψ † f ∈ P p (T ) and ψ † ξ ∈ RT N p (T ). Note also that the multiplicative factor ψ † has been somehow factored out in the right-hand side of (1.3a) which, in particular, entails that the infinite-dimensional minimization set does not enforce any prescription on the normal component of the flux at the boundary. The proof of Theorem 1.2 is again constructive and uses as a key idea the link to a primal formulation as indicated by the right-hand side of (1.3b).

Theorem 1.2 has two important immediate applications. The first arises in the context of parabolic problems with mesh adaptation between time-steps. Previous a posteriori error analyses (see, e.g. [START_REF] Verfürth | A posteriori error estimates for finite element discretizations of the heat equation[END_REF]) required the so-called transition condition, which restricts the extent of mesh-adaptation between time-steps, since the constants in the efficiency of the estimators typically depended on the ratio of the sizes of the elements between the different meshes. However, in [START_REF] Ern | Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems[END_REF], we were able to remove this restriction for the first time, thanks to Theorem 1.2 which enables the construction of equilibrated flux a posteriori error estimates with efficiency bounds that do not depend on the mesh adaptation between time-steps. We refer the reader to [START_REF] Ern | Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems[END_REF] for further details. A second application of Theorem 1.2 arises in the context of nonmatching meshes: following [START_REF] Dolejší | hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF], Theorem 1.2 enables the construction of equilibrated flux a posteriori error estimates without any restriction on the number of levels of hanging nodes. We give the details in Appendix A.

A posteriori error analysis of problems with H -1 source terms

The study of discrete H(div)-liftings is connected to the a posteriori error analysis of elliptic problems with source terms in H -1 (Ω), as we now explain. We start by noting that the righthand side of (1.2) above corresponds to the H 1 -seminorm of the solution of the model problem

-∆u = f + ∇•ξ in Ω, u = 0 on Γ D , ∇u•n = -ξ•n on Γ N .
(1.4)

We recall that for a general ξ ∈ L 2 (Ω), the source term f + ∇•ξ and the Neumann boundary condition in (1.4) must be interpreted in the sense of distributions. The weak formulation of (1.4) then looks for u ∈ H 1 * (Ω) such that

(∇u, ∇v) = (f, v) -(ξ, ∇v) ∀ v ∈ H 1 * (Ω). (1.5)
In this work, we show that the a posteriori error analysis of problem (1.4) leads to an inherently constructive, practical, and efficient way of computing discrete p-robust H(div)-liftings, thereby justifying Theorems 1.1 and 1.2. Furthermore, we study the a posteriori error analysis of (1.4) independently, in the most general case f ∈ L 2 (Ω) and ξ ∈ L 2 (Ω), because elliptic problems with distributional source terms arise in many other important applications, such as the computation of scalar potentials in the Helmholtz decomposition of vector fields.

In comparison to the literature on a posteriori error estimates for problems with source terms in L 2 (Ω), there are comparatively few works treating the case of distributional data. Cohen, DeVore, and Nochetto [START_REF] Cohen | Convergence rates of AFEM with H -1 data[END_REF] propose a posteriori error estimates involving the sum of localized negative norms of the source term over the patches of the mesh and weighted norms of the jumps in the finite element solution over the faces of the mesh. However, it is known from examples [7, p. 704] that this estimator can significantly overestimate the error in some cases; this is due to the splitting of the residual into the source term and the jumps of the numerical solution. Recently, Kreuzer and Veeser [START_REF] Kreuzer | Oscillation in a posteriori error estimation[END_REF] derived a posteriori error estimates based on low-pass filters that are both reliable and efficient. Furthermore, the a posteriori error analysis of problems with Dirac delta source terms, which do not lie in [START_REF] Agnelli | A posteriori error estimates for elliptic problems with Dirac measure terms in weighted spaces[END_REF][START_REF] Araya | A posteriori error estimates for elliptic problems with Dirac delta source terms[END_REF][START_REF] Gaspoz | A posteriori error estimates with point sources in fractional Sobolev spaces[END_REF].

H -1 if d ≥ 2, is treated, for instance, in
The second main contribution of this work is to extend the results of [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF][START_REF]Stable broken H 1 and H(div) polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions[END_REF], where locally efficient and p-robust equilibrated flux error estimators for elliptic problems with sources in L 2 (Ω) are derived, to problems with source terms in H -1 (Ω) such as (1.4). In particular, we prove the following result (detailed notation is given in Section 2):

Theorem 1.3. Let f ∈ L 2 (Ω) and ξ ∈ L 2 (Ω), with (f, 1) = 0 if Γ N = ∂Ω, and let u ∈ H 1
* (Ω) be the weak solution of problem (1.4) defined by (1.5). Let p and p be positive integers with 1 ≤ p ≤ p, and let u h ∈ V h := H 1 * (Ω) ∩ P p (T ) be the finite element approximation of u such that

(∇u h , ∇v h ) = (f, v h ) -(ξ, ∇v h ) ∀ v h ∈ V h . (1.6)
Let the equilibrated flux reconstruction σ h ∈ H(div, Ω) ∩ RT N p (T ) be defined by (3.1) and (3.4) below. Then, we have the guaranteed upper bound on the error

∇(u -u h ) 2 ≤ K∈T σ h + ξ + ∇u h K + h K π f -Π hp f K 2 . (1.7)
Furthermore, for each K ∈ T , we have the local efficiency bound

σ h + ξ + ∇u h K a∈V K ∇(u -u h ) ωa + η a osc , (1.8)
where the local data oscillation η a osc defined by (4.2) below. Finally, the global efficiency can be summarized as

σ h + ξ + ∇u h ∇(u -u h ) + a∈V [η a osc ] 2 1 2
.

(1.9)

We consider in Theorem 1.3 the polynomial degrees 1 ≤ p ≤ p; in the a posteriori analysis of the model problem (1.5), one is typically interested in the situation where p = p. However, in the context of discrete p-robust H(div)-liftings, the approximation u h only serves as a tool in the analysis, and the choice p = 1 turns out to be sufficient for our purposes.

The rest of this paper is organized as follows. We detail the notation in Section 2. Section 3 then presents the equilibrated flux reconstruction in the setting of H -1 source terms. Sections 4, 5, and 6 respectively prove Theorems 1.3, 1.1, and 1.2, and Appendix A illustrates an application of our results to a posteriori error estimation on meshes with an arbitrary number of levels of hanging nodes.

Setting

We summarize here briefly the notation used in this paper. For an arbitrary open subset ω ⊂ Ω, we use (•, •) ω to denote the L 2 -inner product for scalar-or vector-valued functions on ω, with associated norm • ω . In the special case where ω = Ω, we drop the subscript notation, i.e.

• := • Ω . For each mesh element K ∈ T and for a fixed integer p ≥ 1, let P p (K) denote the space of polynomials of total degree at most p on K. Let P p (T ) ⊂ L 2 (Ω) denote the space of scalar piecewise-polynomials of degree at most p over T and let RT N p (T ) ⊂ L 2 (Ω) denote the piecewise Raviart-Thomas-Nédélec space, defined by RT N p (T ) :

= {v h ∈ L 2 (Ω), v h | K ∈ RT N p (K)}, where RT N p (K) := P p (K; R d ) + P p (K)x. Let Π RT N
hp denote the vector L 2 (Ω)orthogonal projection operator from L 2 (Ω) onto RT N p (T ). Let Π hp : L 2 (Ω) → P p (T ) denote the scalar L 2 -orthogonal projection operator from L 2 (Ω) onto P p (T ). Finally, let Π hp denote the vector L 2 (Ω)-orthogonal projection operator from L 2 (Ω) onto P p (T ; R d ); note that Π hp is simply obtained by application of Π hp component-wise. We also emphasize that all these projections are elementwise, in particular since H(div, Ω)-conformity is not imposed on the space RT N p (T ).

Let F denote the set of faces of the mesh, with F ext denoting the set of all boundary faces of the mesh. For each element K ∈ T , h K stands for the diameter of K. Let V denote the set of vertices of the mesh T . For each a ∈ V, the function ψ a is the hat function associated with a, and the set ω a is the interior of the support of ψ a , with associated diameter h ωa . Furthermore, let T a denote the restriction of the mesh T to ω a . In the case where Γ D is nontrivial, a vertex a ∈ V is said to belong to V int , the set of interior and Neumann boundary vertices, if a ∈ Ω ∪ (∂Ω \ Γ D ). Otherwise a belongs to V ext , the set of Dirichlet boundary vertices. In the case where Γ D = ∅ and Γ N = ∂Ω, all vertices are considered to be interior vertices and V int := V. Finally, for each element K ∈ T , we collect in V K the set of vertices of V belonging in K.

Flux equilibration for elliptic problems with H -1 source terms

The construction of the flux equilibration is based on independent local mixed finite element approximations of residual problems over the patches of elements around mesh vertices, in generalization of [START_REF] Ainsworth | A framework for obtaining guaranteed error bounds for finite element approximations[END_REF][START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF][START_REF] Braess | Equilibrated residual error estimator for edge elements[END_REF][START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF][START_REF]Stable broken H 1 and H(div) polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions[END_REF]. For each a ∈ V, let P p (T a ), respectively RT N p (T a ), be the restriction of P p (T ), respectively RT N p (T ), to the patch T a around the vertex a ∈ V. The local spatial mixed finite element spaces V a h and Q a h are defined by

V a h := {v h ∈ H(div, ω a ) ∩ RT N p (T a ), v h •n = 0 on ∂ω a } if a ∈ V int , {v h ∈ H(div, ω a ) ∩ RT N p (T a ), v h •n = 0 on ∂ω a \ Γ D } if a ∈ V ext , Q a h := {q h ∈ P p (T a ), (q h , 1) ωa = 0} if a ∈ V int , P p (T a ) if a ∈ V ext .
For each a ∈ V, let σ a h ∈ V a h be defined by

σ a h := arg min v h ∈V a h ∇•v h =g a h v h + ψ a (ξ + ∇u h ) ωa , (3.1)
where

g a h := Π hp ψ a f -∇ψ a •(ξ + ∇u h ) | ωa = (Π hp (ψ a f ) -∇ψ a • (Π hp ξ + ∇u h )) | ωa , (3.2)
where the last equality follows from Π hp (∇ψ a •ξ) = ∇ψ a •Π hp ξ, since ∇ψ a is piecewise constant, and since u h is a piecewise polynomial of degree at most p ≤ p. It is important to note that g a h satisfies the Neumann compatibility condition (g a h , 1) ωa = 0 for all a ∈ V int , i.e. g a h ∈ Q a h , thereby guaranteeing that σ a h from (3.1) is well-defined. Indeed, this is found by choosing the test function

v h = ψ a in (1.6) when Γ D = ∅, as here ψ a ∈ V h ⊂ H 1 * (Ω). When Γ D = ∅, ψ a / ∈ V h ⊂ H 1 * (
Ω) due to the mean-value zero condition on H 1 * (Ω), but the compatibility condition (f, 1) = 0 implies that (1.6) also holds for ψ a .

It is well-known that the Euler-Lagrange conditions for (3.1) are: find

σ a h ∈ V a h and r a h ∈ Q a h
(the Lagrange multiplier of the divergence constraint) such that

(σ a h , v h ) ωa -(∇•v h , r a h ) ωa = -(ψ a (ξ + ∇u h ), v h ) ωa ∀ v h ∈ V a h , (3.3a)
(∇•σ a h , q h ) ωa = (ψ a f -∇ψ a •(ξ + ∇u h ), q h ) ωa ∀ q h ∈ Q a h . (3.3b)
After extending each σ a h by zero in Ω \ ω a , we define the equilibrated flux reconstruction σ h ∈ RT N p (T ) by

σ h := a∈V σ a h . (3.4)
Lemma 3.1. Let the equilibrated flux reconstruction σ h ∈ RT N p (T ) be defined by (3.4). Then, σ h belongs to H(div, Ω) and satisfies

∇•σ h = Π hp f in Ω, (3.5a) σ h •n = 0 on Γ N . (3.5b)
Proof. The proof follows closely the arguments in [START_REF] Braess | Equilibrated residual error estimator for edge elements[END_REF][START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF][START_REF] Dolejší | hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]; we sketch it here for the sake of completeness. First, for any a ∈ V, the zero extension of σ a h belongs to H(div, Ω) as a result of the boundary conditions in V a h . Thus σ h ∈ H(div, Ω). Consider now any element K ∈ T having a face F contained in Γ N . Then, for each vertex a ∈ V K , the definition of V a h requires that σ a h •n = 0 on F , thus implying that σ h •n = 0 on F . Since Γ N = ∪ F ⊂ΓN F by hypothesis, we deduce (3.5b). Finally, to show (3.5a), we employ (3.1) and (3.2): thus, on each

K ∈ T a , ∇•σ h | K = a∈V K ∇•σ a h | K = a∈V K Π hp ψ a f -∇ψ a •(ξ + ∇u h ) | K = Π hp f | K , since the hat functions {ψ a } a∈V K form a partition of unity over K.
For a given vertex a ∈ V, let the space H 1 * (ω a ) be defined by

H 1 * (ω a ) := {v ∈ H 1 (ω a ), (v, 1) ωa = 0} if a ∈ V int , {v ∈ H 1 (ω a ), v| ∂ωa∩ΓD = 0} if a ∈ V ext .
We also set τ a h := Π RT N hp (ψ a ξ) + ψ a ∇u h . Then we have the following crucial stability result.

Lemma 3.2. For each a ∈ V, let σ a h be defined by (3.1). Then

σ a h + τ a h ωa min σ∈V a ∇•σ=g a h σ + τ a h = max v∈H 1 * (ωa)\{0} (g a h , v) ωa -(τ a h , ∇v) ωa ∇v ωa , (3.6)
where V a denotes the set of all vector fields σ ∈ H(div, ω a ) such that σ•n = 0 either on ∂ω a if a ∈ V int , or on ∂ω a \ Γ D if a ∈ V ext .

Proof. It follows from the definitions of the projectors

Π RT N hp that (τ a h , v h ) ωa = (ψ a (ξ + ∇u h ), v h ) ωa for all v h ∈ V a h . Therefore, (3.3a) implies that (σ a h , v h ) ωa -(∇•v h , r a h ) ωa = -(τ a h , v h ) ωa for all v h ∈ V a h . We deduce that (3.1) is equivalent to σ a h = arg min v h ∈V a h ∇•v h =g a h v h + τ a h ωa .
Since g a h ∈ P p (T a ) and since τ a h ∈ RT N p (T a ), we obtain (3.6) from [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF]Thm. 7] in the case of two space dimensions (up to straightforward adaptations for boundary vertices), and [13, Thm. 2.2] in the case of three space dimensions.

Proof of Theorem 1.3

The proof follows essentially the arguments in [START_REF] Braess | Equilibrated residual error estimator for edge elements[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF].

Proof of the guaranteed upper bound (1.7)

It is straightforward to see from the fact that u h ∈ H 1 * (Ω) and from (1.5) that

∇(u -u h ) = max v∈H 1 * (Ω)\{0} (f, v) -(ξ, ∇v) -(∇u h , ∇v) ∇v .
Consider v ∈ H 1 * (Ω) such that ∇v = 1. Then, by addition and subtraction and the facts that ∇•σ h = Π hp f and that σ h •n = 0 on Γ N by Lemma 3.1, we have

(f, v) -(ξ, ∇v) -(∇u h , ∇v) = (f -Π hp f, v) -(σ h + ξ + ∇u h , ∇v). (4.1)
Next, we note that since f -Π hp f has mean-value zero over each K ∈ T , we obtain the bound

|(f -Π hp f, v) K | ≤ h K π f -Π hp f K ∇v K for each K ∈ T

Proof of the local efficiency (1.8)

Let us define the data oscillation as

[η a osc ] 2 := K∈T a h 2 K p 2 ψ a f -Π hp (ψ a f ) 2 K + ξ -Π hp ξ 2 K + ψ a ξ -Π RT N hp (ψ a ξ) 2 K . (4.2)
Recall the definition τ a h := Π RT N hp (ψ a ξ) + ψ a ∇u h and note that, for each K ∈ T ,

σ h + ξ + ∇u h K ≤ a∈V K σ a h + ψ a (ξ + ∇u h ) K ≤ a∈V K [ σ a h + τ a h K + ψ a ξ -Π RT N hp (ψ a ξ) K ] ≤ a∈V K [ σ a h + τ a h ωa + η a osc ].
Employing Lemma 3.2, it is enough to bound the right-hand side (3.6). Fix a ∈ V K and consider v ∈ H 1 * (ω a ) such that ∇v ωa = 1. We then write (g a h , v) ωa -(τ a h , ∇v) ωa = 4 i=1 E i , where

E 1 := (f, ψ a v) ωa -(ξ + ∇u h , ∇(ψ a v)) ωa , E 2 := (Π hp (ψ a f ) -ψ a f, v) ωa , E 3 := ([ξ -Π hp ξ]•∇ψ a , v) ωa , E 4 := (ψ a ξ -Π RT N hp (ψ a ξ), ∇v) ωa .
Extend ψ a v by zero outside of ω a ; then two situations may arise. Either Γ N = ∂Ω, where we use the fact that the weak formulation (1.5) holds for all test functions in H 1 (Ω) by the compatibility (f, 1) ωa = 0. Or Γ D = ∅, where

ψ a v ∈ H 1 * (Ω) = H 1 ΓD (Ω) for any v ∈ H 1 * (ω a ).
In both cases, we conclude from (1.5) that (∇u, ∇(ψ a v)) ωa = (f, ψ a v) ωa -(ξ, ∇(ψ a v)) ωa . Therefore, we see that

E 1 = (∇(u -u h ), ∇(ψ a v)) ωa and thus |E 1 | ≤ ∇(u -u h ) ωa ∇(ψ a v) ωa .
Next, we recall that there is a constant depending only on the mesh shape-regularity such that ∇(ψ a v) ωa ∇v ωa for all v ∈ H 1 * (ω a ) for every a ∈ V, see [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]. Hence, using the hypothesis that ∇v ωa = 1, we find that

|E 1 | ∇(u-u h ) ωa . Next, we have E 2 = (Π hp (ψ a f )-ψ a f, v-Π hp v) by orthogonality
of the L 2 -projection, and thus we find that

|E 2 | 2 K∈T a h 2 K p 2 ψ a f -Π hp (ψ a f ) 2 K by using the approximation bound v -Π hp v K h K p ∇v K for all v ∈ H 1 (K). Finally, we find that |E 3 | ξ -Π hp ξ ωa , and that |E 4 | ≤ ψ a ξ -Π RT N hp (ψ a ξ) ωa .
This completes the proof of local efficiency (1.8).

Proof of the global efficiency (1.9)

For each K ∈ T , we have (

σ h + ξ + ∇u h )| K = a∈V K [σ a h + ψ a (ξ + ∇u h )]| K
by the partition of unity. Noting that each element K has (d + 1) vertices collected in the set V K since T is a simplicial mesh, the Cauchy-Schwarz inequality implies that

σ h + ξ + ∇u h 2 = K∈T σ h + ξ + ∇u h 2 K ≤ K∈T (d + 1) a∈V K σ a h + ψ a (ξ + ∇u h ) 2 K = (d + 1) a∈V σ a h + ψ a (ξ + ∇u h ) 2 ωa .
It follows from the arguments of Section 4.2 that σ a h + ψ a (ξ + ∇u h ) ωa ∇(uu h ) ωa + η a osc . Therefore,

σ h + ξ + ∇u h 2 a∈V { ∇(u -u h ) 2 ωa + [η a osc ] 2 } ∇(u -u h ) 2 + a∈V [η a osc ] 2 ,
which completes the proof of global efficiency (1.9).

Proof of Theorem 1.1

We will show here that Theorem 1.1 follows easily from Theorem 1.3, while using the finite element solution u h of (1.6) as an auxiliary ingredient of the proof. Here it is enough to take p = 1 for the polynomial degree in (1.6). Henceforth, suppose that f ∈ P p-1 (T ) and that ξ ∈ RT N p-1 (T ) and construct σ h by (3.4) using the local minimization problems (3.1).

Since

σ h ∈ H(div, Ω) ∩ RT N p (T ) satisfies ∇•σ h = Π hp f = f and σ h •n = 0 on Γ N by Lemma 3.1, we have min v h ∈H(div,Ω)∩RT Np(T ) ∇•v h =f in Ω v h •n=0 on ΓN v h + ξ ≤ σ h + ξ . (5.1)
Therefore it remains only to show that σ h +ξ is bounded by the right-hand side of (1.2). First, for each vertex a ∈ V, note that ψ a f ∈ P p (T a ) and that ψ a ξ ∈ RT N p (T a ). Thus the oscillation terms η a osc given by (4.2) vanish. Consequently, the global efficiency bound (1.9) implies that σ h + ξ + ∇u h ∇(uu h ) . Furthermore, the stability of the Galerkin method (recall that the finite element solution u h of (1.6) is an orthogonal projection of the weak solution u of (1.5) from H 1 * (Ω) to V h ⊂ H 1 * (Ω)) implies that ∇u h ≤ ∇u . Therefore we deduce that

σ h + ξ ≤ σ h + ξ + ∇u h + ∇u h ∇u . (5.2)
Since u is the weak solution of (1.5), the equivalence of primal and dual formulation of elliptic problems implies that ∇u = max

v∈H 1 * (Ω)\{0} (f, v) -(ξ, ∇v) ∇v = min v∈H(div,Ω) ∇•v=f in Ω v•n=0 on ΓN v + ξ . (5.3)
The combination of the bounds (5.1), (5.2), and (5.3) yields (1.2).

6 Proof of Theorem 1.2

As in Section 5, we merely employ the weak solution u of (1.5) and its finite element approximation u h of (1.6) as tools. For this purpose, we now set Γ D = ∂Ω and Γ N = ∅, so that H 1 * (Ω) = H 1 0 (Ω), and we choose the auxiliary polynomial degree p = 1. We will construct an equilibrated flux σ † h in the discrete minimization set of the left-hand side of (1.3a) such that σ † h + ψ † ξ is bounded by the right-hand side of (1.3b). The key idea is to write

σ † h := σ h + σ c h , σ h := a∈V w a σ a h , (6.1)
where σ h is an uncorrected high-order flux obtained from the local fluxes σ a h , for all a ∈ V (see Subsection 6.1), and where σ c h is a global, low-order, correction flux (see Subsection 6.2), and the weights w a result from

ψ † = a∈V w a ψ a in Ω, w a = ψ † (a), ∀ a ∈ V. (6.2)
The correction term σ c h is needed, since it will be found below that ∇• σ h = ψ † f -∇ψ † •(ξ + ∇u h ); thus we shall build σ c h (by posing a global low-order minimization problem) so that it satisfies ∇•σ c h = ∇ψ † •∇u h , in order to ensure that σ † h satisfies the divergence constraint required by the discrete minimization set of the left-hand side of (1.3a) (see subsection 6.3). The stability properties of σ c h then rely on Theorem 1.1 (in the low-order case p = 1), whereas the stability properties of σ a h are established by using similar ideas to those that were used to prove Theorem 1.3.

Construction of σ h

We construct locally a higher-order σ a h ∈ H(div, ω a ) ∩ RT N p (T a ) for each a ∈ V, similarly to the construction of σ a h from Section 3. We do so in the context Γ D = ∂Ω and Γ N = ∅; consequently, the sets of vertices simplify to interior and boundary ones V int = {a ∈ V; a ∈ Ω} and V ext = {a ∈ V; a ∈ ∂Ω}. For each a ∈ V, let the subset Γ a ⊂ ω a be given by those faces where the hat function ψ a vanishes, i.e. Γ a := {x ∈ ∂ω a , ψ a (x) = 0}. Equivalently, Γ a is composed of all faces F ∈ F that are contained in ∂ω a and that do not contain the vertex a; we denote this corresponding set of faces by F a . If a ∈ V int is an interior vertex, then Γ a = ∂ω a . If a ∈ V ext , then any interior face F ∈ F \ F ext such that F ⊂ ∂ω a necessarily belongs to F a . If F a only consists of such interior faces, then Γ a = ∂ω a \ Γ D = ∂ω a \ ∂Ω. However, Γ a and ∂ω a \ Γ D do not coincide for a vertex a where Γ a includes boundary faces F ⊂ Γ D opposite to the vertex a. In any case, Γ a is always a strict subset of ∂ω a for boundary vertices.

For each a ∈ V, we let the discrete space V a h be defined by

V a h := {v h ∈ H(div, ω a ) ∩ RT N p (T a ), v h •n = 0 on Γ a } ∀ a ∈ V.
For interior vertices, V a h coincides with the space V a h of Section 3. However, in general V a h = V a h for boundary vertices. In this section, V a h is used instead of V a h for a technical point appearing in the proof of Lemma 6.2 below concerning (6.8b) in situations where some elements may have all faces belonging to the boundary. For each a ∈ V, analogously to (3.1), we let σ a h ∈ V a h be defined by σ a h := arg min

v h ∈ V a h ∇•v h =g a h v h + ψ a (ξ + ∇u h ) ωa , (6.3)
where g a h is given by (3.2). Note that σ a h is well-defined for all interior vertices a ∈ V int , since (g a h , 1) ωa = 0 for a ∈ V int . In the case of boundary vertices a ∈ V ext , σ a h is also well-defined, since there are always at least some faces of ∂ω a that are not in Γ a . Finally, the extension by zero of σ a h from ω a to Ω is again H(div, Ω)-conforming. The uncorrected high-order flux function σ h ∈ RT N p (T ) ∩ H(div, Ω) is then defined by

σ h := a∈V w a σ a h . (6.4)
We have the following key result.

Lemma 6.1. Let σ a h be defined by (6.3). Then, we have

σ a h + ψ a (ξ + ∇u h ) ωa ∇(u -u h ) ωa . (6.5)
Proof. Similarly to the proof of Lemma 3.2, we apply [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF]Thm. 7] in the case of two space dimensions, and [13, Thm 2.2] in the case of three space dimensions, to deduce that

σ a h + ψ a (ξ + ∇u h ) ωa max v∈ H 1 * (ωa)\{0} (f, ψ a v) ωa -(ξ + ∇u h , ∇(ψ a v)) ωa ∇v ωa ,
where

H 1 * (ω a ) := H 1 (ω a )/R if a ∈ V int , and H 1 * (ω a ) := H 1 ∂ωa\Γa (ω a ) if a ∈ V ext . For any a ∈ V, we have ψ a v ∈ H 1 0 (ω a ) for any v ∈ H 1 * (ω a )
, where we use the fact that ψ a vanishes on Γ a by definition in the case of a ∈ V ext . We then use (1.5), recalling that Γ D = ∂Ω and that

H 1 * (Ω) = H 1 0 (Ω) here, so that (f, ψ a v) ωa -(ξ + ∇u h , ∇(ψ a v)) ωa = (∇(u -u h ), ∇(ψ a v)) ωa .
Thus we obtain (6.5) from an application of the Cauchy-Scwharz inequality and the Poincaré inequality, using ∇(ψ a v) ωa ∇v ωa for v ∈ H 1 * (ω a ), as in Section 4.2.

Construction of σ c h

We now select a global low-order correction

σ c h ∈ H(div, Ω) ∩ RT N p (T ) such that ∇•σ c h = ∇ψ † •∇u h in Ω, σ c h •n = 0 on Γ N, †
, and such that

σ c h max v∈H 1 † (Ω)\{0} (∇ψ † •∇u h , v) ∇v . ( 6
∈ V h if Γ N, † = ∂Ω), thereby yielding (∇ψ † •∇u h , 1) = (∇u h , ∇ψ † ) = (f, ψ † ) -(ξ, ∇ψ † ) = 0,
where the last identity is obtained from the hypothesis on f and ξ of Theorem 1.2 for the case Γ N, † = ∂Ω. It follows from (6.6) that

σ c h ≤ C P, † h Ω ∇ψ † ∞ ∇u h , (6.7)
where C P, † is the constant of the Poincaré inequality for H 1 † (Ω) and h Ω is the diameter of Ω.

Admissibility of σ † h

Recalling the definitions of σ h from Subsection 6.1 and of σ c h from Subsection 6.2, we define σ † h ∈ H(div, Ω) ∩ RT N p (T ) by (6.1). We now check that σ † h belongs to the minimization set of the left-hand side of (1.3b). Lemma 6.2. Let σ † h be defined by (6.1). Then

∇•σ † h = ψ † f -∇ψ † •ξ in Ω, (6.8a) σ † h •n = 0 on Γ N, † . (6.8b)
Proof. Since f ∈ P p-1 (T ) and since ξ ∈ RT N p-1 (T ) ⊂ P p (T ; R d ), we find that ∇• σ a h = ψ a f -∇ψ a •(ξ + ∇u h ) from (6.3) and (3.2). In consequence of (6.2) and of the definition of σ h in (6.4), by proceeding as in the proof of Lemma 3.1, we then obtain ∇• σ h = ψ † f -∇ψ † •(ξ + ∇u h ) in Ω. Therefore, using the definition of the flux correction σ c h , it is found that

∇•σ † h = ∇• σ h + ∇•σ c h = ψ † f -∇ψ † •ξ,
thereby showing (6.8a). Next, we show (6.8b). Recall that σ c h •n = 0 on Γ N, † by construction, and that Γ N, † = ∪ F ∈F † , where F † is the set of faces F ∈ F ext such that ψ † | F = 0. Consider a face F ∈ F † , with its corresponding element K ∈ T . Then, for any of the d vertices a ∈ F , we have w a = ψ † (a) = 0 by definition. Therefore w a σ a h •n = 0 on F trivially for all a ∈ F . Since K is a simplex, there is a unique remaining vertex a opposing the face F . Thus F ⊂ Γ a and consequently σ a h •n = 0 on F by the definition of the space V a h and of the set Γ a = {x ∈ ∂ω a , ψ a (x) = 0}. In summary, for any F ∈ F † , we have w a σ a h •n = 0 on F for all a ∈ V K , where F ⊂ K, whence the assertion (6.8b) follows.

Proof of (1.3b)

Recall that we consider the solution u of (1.5) and u h of (1.6) in the context Γ D = ∂Ω and Γ N = ∅, so that H 1 * (Ω) = H 1 0 (Ω). Thus, we obtain the bound

∇u h ≤ ∇u = max v∈H 1 0 (Ω)\{0} (f, v) -(ξ, ∇v) ∇v . (6.9)
Therefore, our last goal is to show that σ † h + ψ † ξ can be bounded in terms of ∇u h and ∇u . For each K ∈ T , we have σ

† h | K = σ c h | K + a∈V K (w a σ a h )| K from (
∇(u -u h ) 2 ωa ψ † 2 ∞ ∇(u -u h ) 2 , (6.11)
where the first inequality uses the fact that ψ † ∈ H 1 (Ω) ∩ P 1 (T ) and that the coefficients {w a } a∈V are the nodal values of ψ † at the vertices of the mesh. Therefore, the combination of (6.10) with (6.7) and (6.11) yields

σ † h + ψ † ξ C P, † h Ω ∇ψ † ∞ ∇u h + ψ † ∞ ( ∇(u -u h ) + ∇u h ).
We finally deduce (1.3b) from the above inequality and from (6.9).

A Application to a posteriori error estimates on meshes with hanging nodes

Equilibrated flux a posteriori error estimates for meshes with hanging nodes are developed in [START_REF] Dolejší | hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF] where the equilibration is performed on patches ω a corresponding to the support of hat functions ψ a associated with non-hanging nodes a of the computational mesh and forming a partition of unity of the computational domain, see [START_REF] Dolejší | hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF]Assumption 2.1]. It turns out that a slight extension of the equilibration patch ω a of [START_REF] Dolejší | hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF] enables the removal of the usual dependence of a posteriori efficiency constants on the number of levels of hanging nodes, thereby allowing for a completely arbitrary number of levels of hanging nodes. More precisely, it suffices to extend the equilibration patch ω a so that all the products h ωa ∇ψ a ∞,ωa are uniformly bounded. Then, applying Theorem 1.2 where the patch ω a of [START_REF] Dolejší | hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF] corresponds here to the domain Ω, and the hat function ψ a of [START_REF] Dolejší | hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF] corresponds here to the function ψ † , we infer that the factor h ωa max â∈ Va ∇ψ â ∞,ω â of [START_REF] Dolejší | hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF]Theorem 3.12] can be replaced by the factor h ωa ∇ψ a ∞,ωa , i.e.

h Ω ∇ψ † ∞ in the present notation. The extension of the equilibration patch is illustrated in Figure A. This extension typically entails including several layers of fine elements, so as to ensure that the factors h ωa ∇ψ a ∞,ωa are uniformly bounded. The price to pay to achieve robustness with respect to the level of hanging nodes is thus a somewhat more expensive computation of the equilibrated flux. The proof of Theorem 1.2 in Section 6 shows that this cost can be significantly reduced to the solution of two low-order systems over the extended patch, followed by local high-order corrections within the extended patch.

  by the Poincaré inequality. The upper bound (1.7) then follows from (4.1) and the Cauchy-Schwarz inequality.

Figure 1 :

 1 Figure1: Original equilibration patches of reference[START_REF] Dolejší | hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF] (left) and extended equilibration patches necessary for estimates robust with respect to an arbitrary number of levels of hanging nodes (right)

  .6) This is possible by applying Theorem 1.1. Indeed, we employ it in the setting where the Neumann part of the boundary ∂Ω is given by Γ N, † , where the Dirichlet part is ∂Ω \ Γ N, † , and where the scalar datum is given by ∇ψ † •∇u h ∈ P p -1 (T ) and the vector datum is zero. The data compatibility condition for the case where Γ N, † = ∂Ω is guaranteed, since we can then choose ψ † as the test function in (1.6) (recall that ψ †

  6.1). So, the triangle inequality and the Cauchy-Schwarz inequality imply thatK∈T a∈V K |w a | 2 σ a h + ψ a (ξ + ∇u h ) 2 K + ψ † ∇u h 2 , (6.10)where we have also used (6.2) to obtainψ † (ξ +∇u h )| K = a∈V K [w a ψ a (ξ +∇u h )]| K . Lemma6.1 and a counting argument show that K∈T a∈V K |w a | 2 σ a h + ψ a (ξ + ∇u h ) 2 K ≤ ψ †

	σ † h + ψ † ξ 2	σ c h	2 + (d + 1)
				2 ∞	a∈V	σ a h + ψ a (ξ + ∇u h ) 2 ωa
			ψ †	2 ∞
				a∈V

* This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 647134 GATIPOR).


P S Kermarrec

A S Zorko

F Bert

A S Colman

B Koteswararao

F Bouquet

Pierre Bonville

A S Hillier

A Amato

J Van Tol

E Kermarrec

R H Colman

A Ozarowski

A S Wills

P Mendels

Spin dynamics and disorder effects in the S = 1/2 kagome Heisenberg spin liquid phase of kapellasite

Keywords: number(s): 75.10.Kt, 76.60.-k, 76.30.-v, 76.75.+i

I. INTRODUCTION

The S = 1 2 Heisenberg antiferromagnet model on a kagome lattice has proven to be a fascinating and enduring problem for more than 20 years. Currently, there seems to be a broad consensus on the nature of the ground state, with recent density matrix renormalization group (DMRG) studies pointing to a resonating valence bond state that features a singlet-triplet gap, although the proposed value of = 0.13J [1,2] is still uncertain [3]. Contrastingly, fermionic approaches which suggest a gapless ground state remain debated [4][5][6][7]. Meanwhile, on the experimental side, new spin-liquid candidates have been discovered in recent years with S = 1 2 spins [8][9][10][11]. One of the most puzzling results, in view of the DMRG results, is the absence of any gap in the excitations of these experimental realizations. The reasons for such a discrepancy, and whether it is related to deviations from the ideal model, remain key open issues in the field. Several research directions have been proposed based on both theoretical developments and experiments such as Dzyaloshinskii-Moriya anisotropic interaction and next nearest-neighbors (NNN) exchange [12][13][14][15], which are believed to drive the system towards a magnetically ordered phase [16][17][18], or lattice distortions.

In this landscape of new spin-liquid candidates and drive to understand the role of perturbations to the ideal nearestneighbor (NN) Heisenberg antiferromagnetic model, kapellasite [9,19], a polymorph of herbertsmithite γ -Cu 3 Zn(OH) 6 Cl 2 , has come as a surprise. Instead of being a mere variant of herbertsmithite, it has been demonstrated to be a kagome spin-liquid candidate of a new type. μSR and inelastic neutron scattering (INS) experiments evidenced a gapless spin-liquid state down to 20 mK, with unusual dynamical short-range spin correlations. It is an experimental representative of a previously uninvestigated class of models on the kagome lattice with frustrated competing interactions, namely, ferromagnetic NN exchange J 1 and second NN (along the edges) J 2 , and antiferromagnetic NNN exchange (across the diagonal of the hexagons) J d , with J 2 J 1 ∼ J d [20,21]. The values of interactions extracted from a high-T series analysis of macroscopic susceptibility and specific heat lead to an estimate of J 1 = -12 K, J 2 = -4 K, and J d = 15.6 K [20,22], which, for classical spins, would stabilize a noncoplanar chiral state of the cuboc2 type. In kapellasite, quantum fluctuations are argued to destabilize the latter state and to lead to related dynamical correlations, as detected by INS [20].

The magnetic structure of kapellasite, α-Cu 3 Zn(OH) 6 Cl 2 , exhibits well-decoupled magnetic planes with a kagome geometry formed from the depletion of a Cu 2+ (S = 1 2 ) triangular lattice by diamagnetic Zn 2+ [Fig. 1(a)]. These planes are only weakly coupled along the third dimension (the hexagonal c axis) via interlayer Cl. . .OH hydrogen bonding, which guarantees highly bidimensional characteristics. This essential feature causes the flat shape of micron scale grains constituting our powder sample [Fig. 1 revealed a nondistorted kagome lattice, but showed a sizable shortage of Cu compared to the ideal stoichiometry [9]. The observed formula determined by ICP-AES chemical analysis is (Cu 0.73 Zn 0.27 ) 3 (Zn 0.88 Cu 0.12 )(OH) 6 Cl 2 . Importantly, the local octahedral environment of nominally Cu site is more distorted than the Zn one, and has an extension along the Cl-Cl direction [Fig. 1(c)]. The more symmetric bonding environment of the Zn site, an octahedral site coordinated by six oxygens, leads to strong preferential occupation by Zn 2+ , with 95% of the Jahn-Teller active Cu 2+ residing on the more distorted kagome sites coordinated by four oxygens and two chlorines. The major consequence of this nonstoichiometry is then the depletion of the kagome lattice by nonmagnetic Zn 2+ ions. The kagome site occupancy p = 0.73 still remains above the percolation threshold of the kagome lattice p c 0.65 [23,24].

In this paper, we investigate the magnetic properties and the effects of structural disorder in kapellasite by means of NMR, μSR, ESR, and specific-heat measurements. This paper is organized as follows. In Sec. II, we evaluate by 35 Cl-NMR the effect of disorder within the kagome planes and report evidence for local homogeneity of the magnetic properties, proving the robustness of the underlying physics against dilution. In Sec. III, we use ESR to estimate the symmetric exchange anisotropy term in the Hamiltonian and confirm that the choice of an isotropic Heisenberg Hamiltonian to describe the magnetic properties of kapellasite is reasonable in a first approach. While the main exchange interactions J 1 and J d have already been derived from magnetic susceptibility and specific-heat fits, we detail in Sec. IV an analysis of the specific heat under magnetic fields which confirms that its low-T part is magnetic and intrinsic, and thus provides a pertinent quantity for extracting the exchange integrals. Finally, in Sec. V, we investigate with NMR and μSR the spin fluctuations at low energy in the low-T frustrated antiferromagnetic phase induced by these competing interactions on the S = 1 2 kagome lattice.

II. NMR DETERMINATION OF LOCAL STRUCTURE AND SPIN SUSCEPTIBILITY

NMR experiments were performed on the most abundant Cl isotope 35 Cl (abundance of 75.8%) using an oriented powder aligned along the hexagonal c axis of the kapellasite's rhombohedral structure. The high-quality orientation leads to accurate refinements of the local structure and unambiguous determination of the local spin susceptibility and the spin dynamics. 35 Cl is a quadrupolar nucleus (I = 3 2 ), hence, one has to consider two sets of interactions resulting from the coupling of the nucleus to its magnetic and charge environments. The nuclear Hamiltonian can be written [START_REF] Slichter | Principles of Magnetic Resonance[END_REF]26]

H = H mag + H quad (1)
with

H quad H mag .
As detailed in Appendix, this generally leads for a given environment to a spectrum consisting of three lines. The central line, corresponding to the -1 2 ↔ 1 2 transition, and the satellite lines, corresponding to -3 2 ↔ -1 2 and 1 2 ↔ 3 2 transitions, which are sensitive to quadrupolar perturbations in second and first order, respectively. More details about this nuclear Hamiltonian can be found in the Appendix.

A. Local structure

In this section, we discuss the local structure of kapellasite as deduced from 35 Cl NMR experiments. For a perfect coverage of the kagome lattice with Cu 2+ ions, one would expect a unique Cl environment, where each Cl nucleus couples to the three adjacent Cu sites located at the corners of a given triangle. In disagreement with this simple scheme, we instead observe multiple lines, which are best resolved for the central - 1 2 ↔ 1 2 transition line at T = 30 K taken for θ = 0 • , where θ is the angle between the applied field and the c axis (see Fig. 2 and Sec. II B). The multiple lines can be associated with different Cl environments induced by the Cu/Zn mixing disorder, namely, four different Cu/Zn configurations on a probed triangle: Cu 3 (1), Cu 2 Zn (2), CuZn 2 (3), and Zn 3 (4). This site assignment also agrees with the high-T shifts which scale with the number of Cu to which each Cl is coupled [see Fig. 4(b)]. Since the spin-spin (T 2 ) and spin-lattice (T 1 ) corrections were found to be negligible, the integrated intensity of each line is directly proportional to the number of triangles of each configuration which is found to perfectly match with a random occupation of the Cu sites (Table I), with a filling ratio p = 0.75 (1). This agrees with the value of 0.73 extracted from neutron diffraction refinements [9].

The impact of dilution on the local structure can be traced by the satellite lines which positions depend to first order upon the parameters which govern H quad , namely, the quadrupolar FIG. 2. (Color online) 35 Cl NMR spectrum of the -1 2 ↔ 1 2 transition lines at T = 30 K, θ = 0 • (empty circles). The shift reference is indicated by the red vertical dotted line. The NMR spectrum is fitted using four Gaussian lines (black line), corresponding to various 35 Cl environments due to Cu/Zn (blue/yellow) mixing. The integrated intensity A i of each line i reveals the site probability. frequency ν q and the asymmetry η, and the orientation of the principal axis of the electric field gradient (EFG) (see Appendix). Hence, they reflect the difference in electrostatic environments at the local scale. A complete NMR spectrum measured in the θ = 0 • orientation is shown in Fig. 3(a). Three groups of NMR lines are clearly observed which correspond to the three expected transitions for a spin I = 3 2 . Only the satellite contribution from the two dominant main sites 1 and 2 can be reliably followed since the sites 3 and 4 represent a minor fraction (altogether 16%) of the total spectrum intensity. Further details on satellites 1 and 2 are displayed in Figs. 3(b) and 3(c). For θ = 0 • , if one assumes that the principal axis of the EFG remains along the symmetry axis Oz for all the Cl sites, the distance between the two satellites should be 2ν q [Eqs. (A9) and (A10)], independently of the η value. Under the latter assumption, we find two different values for sites 1 and 2, respectively, ν q,1 = 3.4(1) MHz and ν q,2 = 4.0(2) 35 Cl NMR spectrum at ν 0 = 47.950 MHz and T = 100 K. Three distinct groups of NMR lines are observed. (b), (c) NMR satellites in two different angle configurations θ = 0 • (b) and θ = 90 • (c) at T = 100 K. The total fit (black line) is the sum of the two satellites pairs of the main sites (site 1 in blue, site 2 in red). For θ = 90 • , the simulation with η 2 = 0 (green line) fails to represent the satellites of the site 2, contrary to the case where η 2 = 0.3(1) (black line).

MHz [Fig. 3(b)], which reflect two different local electrostatic environments around Cl. While for site 1 the spectrum is consistent with a η = 0 value, for site 2 the asymmetry of the spectrum obtained in the perpendicular direction (θ = 90 • ) can only be accounted for by assuming a surprisingly large value of the asymmetry parameter η = 0.3 (1). Finding such a large value might rather invalidate our initial assumption and indicate that for this site, the EFG and shift tensor principal axis do not coincide anymore. The local distortions deduced through this extensive analysis of the quadrupolar parameters can easily be attributed to the different nature of the two Cu 2+ and Zn 2+ ions. Whereas they possess the same charge and a similar radius, Cu 2+ is known to be Jahn-Teller active and to elongate the octahedron built on the nearby anions. Such local effects can then be expected when the ratio Zn/Cu is varied on a given triangle but do not affect the average threefold axis structural symmetry.

B. Spin susceptibility

In the following, we extend the concise presentation of Ref. [20] and detail the analysis of the 35 Cl NMR data taken on the central line, and further compare the susceptibilities extracted for sites 1, 2, and 3 from the line shifts taken for an applied field along the c axis, K z,1 , K z,2 from 287 K down to 1.2 K and K z,3 down to 80 K (Fig. 4). The broadening of the NMR lines observed at low temperature is discussed in Sec. VI A.

Following a typical two-step procedure, we derive the NMR spin shifts which are induced by the local susceptibilities denoted χ i : from a linear high-T fit of K z,i plotted versus the macroscopic susceptibility χ macro measured by SQUID magnetometry, K z,i = A hf,i χ macro + σ i , we extract the hyperfine coupling constants A hf,i and the chemical shifts σ i [Fig. 4(b)]. The local susceptibilities are then given by χ i = (K z,iσ i )/A hf,i . The ratios of hyperfine couplings A hf,1 /A hf,2 = 1.53(2) and A hf,3 /A hf,2 = 0.58 (2), with A hf,1 = -4.98(3) kOe/μ B , are found to be in good agreement with the ratios of the number of Cu 2+ probed, respectively, 1.5 and 0.5. This further strengthens the previous analysis of the local structure, and again indicates that at high T the S = 1 2 spins are rather uncorrelated, depicting a fluctuating paramagnetic state. This validates the overall consistency of our analysis.

The local susceptibility derived for site 1 is certainly the least perturbed by defects and is the closest to that of the ideal system. It is found to scale with χ macro over the whole temperature range, down to 1.2 K, which invalidates the possibility of a dominant contribution of defects at low T in χ macro [Fig. 4(b)]. One should note that this finding strongly contrasts with the case of many correlated oxides where a Curie-type tail in macroscopic susceptibility is commonly attributed to defects, either spin vacancies in the magnetic lattice or other species of weakly coupled spins. In these systems, the disorder contributes to the NMR linewidth only, and not to the line shift [15,[27][28][29].

This uncommon divergence as T → 0 of χ spin in comparison with the reported χ spin of NN Heisenberg kagome antiferromagnets [15,27,30] could only be explained within the framework of competing interactions on the kagome lattice, J 1 = -12 K, J 2 = -4 K, and J d = 15.6 K, evaluated by hightemperature series expansion as emphasized in Refs. [20,22]. The Curie-Weiss temperature θ CW ∼ 9 K is consequently the result of an averaging of these competing interactions.

A closer comparison of the low-T local susceptibilities χ i at sites 1 and 2 brings an additional insight into the nature of the interactions. While all shifts display a low-T upturn, a minor effect is detected below T = 50 K, where the ratio of the local susceptibilities of sites 1 and 2 χ 1 /χ 2 is not T independent anymore and shows instead a clear increase as T is lowered (Fig. 5). In order to infer a physical understanding to this deviation, χ 1 /χ 2 is fitted to the ratio of two Curie-Weiss susceptibilities C(Tθ 2 )/C(Tθ 1 ). The parameters θ 1 and θ 2 account for the interactions, respectively, probed at sites 1 and 2 and are expected to be relevant in the mean field regime, i.e., for T ≥ 2θ CW 20 K. The main result of this simple analysis is that the upturn at low T could only be reproduced when θ 1 > θ 2 . In a mean field approach and ignoring the weakest J 2 coupling, the Curie-Weiss temperature of each site i of the J 1 -J d Heisenberg model on a diluted kagome lattice can be simply expressed as θ i = -α i J 1β i J d /2. For a fully occupied kagome lattice, one has α 1,2 = β 1,2 = 1. Taking into account the random depletion of the lattice by 27%, a straightforward calculation of probabilities leads to α 1 = 0.865, β 1 = 0.73 and α 2 = 0.615, β 2 = 0.73. Therefore, one can directly infer from the relation θ 1 > θ 2 that J 1 < 0 (ferromagnetic), in agreement with Refs. [20,22].

The related scheme of interactions is summarized in the diagram of Fig. 5.

III. ELECTRON SPIN RESONANCE

ESR measurements were performed at NHMFL, Tallahassee, Florida, on a custom-made ESR spectrometer with homodyne detection working in transmission mode at the fixed frequency of ν 0 = 326.4 GHz. The same oriented sample was used as in NMR. A single ESR line was observed in the temperature range between 5 and 100 K. We notice first that this is consistent with the absence of any magnetic parasitic phase in the sample. Moreover, since the copper ions are present on both the nominally Cu and Zn sites, the observation of a single line indicates that there is a sizable exchange coupling between sites, yielding a Lorentzian line shape typical for an exchange narrowing limit.

There are three distinct T regimes as revealed in Fig. 6. Above 100 K, the ESR line could not be followed reliably due to its decreasing intensity with increasing temperature. The reduction of the signal amplitude was additionally enhanced by line broadening. Between 10 and 100 K (intermediate T ), the linewidth increases regularly with increasing T . Below 10 K, an opposite behavior is observed. In the two latter regimes, the resonance field exhibits a monotonic variation, similar to the macroscopic susceptibility.

Although the ESR line in kapellasite is similarly broad to that seen in herbertsmithite, which points to a sizable anisotropy [12], the temperature dependence of the ESR parameters, i.e., the line position B 0 (g factor; g = ν 0 μ B B 0 , where μ B is the Bohr magneton) and the linewidth B, are quite different. In herbertsmithite, a maximum in B and small temperature changes of the g factor of less than 0.02 were observed. Contrastingly, in kapellasite, a minimum is observed in the linewidth at 10 K. This behavior rather resembles that found in another Cu-based kagome compound, vesignieite, where it was assigned to a phononic contribution [13]. Lastly, we note that the measured line shifts with decreasing temperature in kapellasite are much larger than in herbertsmithite and vesignieite. This indicates that a different type of magnetic anisotropy is responsible for these shifts than the Dzyaloshinsky-Moriya (DM) interaction that had been found to be dominant in the other two cases (Fig. 6). Following, we discuss successively the information which can be extracted from the ESR linewidth and from the g factor.

A. Linewidth

In the paramagnetic regime of herbertsmithite (T J ), the ESR linewidth was constant so this asymptotic value could be used to extract the dominant magnetic anisotropy term [12]. It is notable that, in kapellasite, one cannot reach this regime even at 100 K, a temperature already as high as ∼6J d [Fig. 6(b)]. This indicates that the linewidth in the intermediate-T regime cannot be accounted for solely by the spin system and an additional line-broadening mechanism must originate from other degrees of freedom. A spin-phonon coupling, responsible for a similar increase in vesignieite [13], is quite commonly encountered in copper-based systems and is perfectly in line with the increase of the linewidth between 10 and 100 K [31][32][33][34]. This mechanism is dominant at high T so that the magnetic anisotropy cannot be extracted from the linewidth. The low-T regime (below 10 K ∼ J ) resembles more the behavior found in herbertsmithite. At such low temperatures, the phonon population vanishes and the ESR broadening becomes dominantly spin induced, reflecting the development of spin correlations.

B. g factor

Large changes of the g factor with decreasing temperature in kapellasite [Fig. 6(c)] could, in principle, be due to either a locally inhomogeneous susceptibility or a significant magnetic anisotropy of the system. The first possibility can be easily discarded. Indeed, in inhomogeneous spin systems (e.g., spin glasses, impure systems, and systems with several nonequivalent sites) nonuniform local fields develop with increasing susceptibility (decreasing T ) leading exclusively to negative B 0 shifts, irrespective of the orientation of the applied field [35,36]. In addition, line broadening that scales with the susceptibility is generally encountered. Such shifts and broadening result from staggered local fields [36]. In kapellasite, this is clearly not the case since both positive and negative line shifts are observed depending on the orientation of the applied field and the ESR width is nonmonotonic when the susceptibility is increasing steeply at low temperature. We are then led to assign these shifts to magnetic anisotropy.

When the magnetic anisotropy H is small compared to the exchange and/or Zeeman energy, Nagata's theory for g-factor shifts can be applied [37]. Then, the shift is determined by the correlation function [37,38]

g = [S -,[S + ,H ]] 2μ B B 0 S z , (2)
where [. . . , . . .] denotes a commutator, . . . canonical averaging, and S α = i S α i the α component of the total spin operator. Here, we refer to the general form of the magnetic anisotropy in the spin Hamiltonian of a spin- 1 2 system

H = 1 2 i,j =i S i • K ij • S j , ( 3
)
where the coupling tensor

K ij = d ij + D ij (4)
includes both the long-range dipolar coupling d ij and the shortrange symmetric anisotropic exchange D ij on bonds J 1 and J d . We note that the antisymmetric Dzyaloshinsky-Moriya interaction does not contribute to the ESR line shift at the level of first order [Eq. ( 2)] [38] and so we omit it. We proceed to the evaluation of both anisotropy terms from Eq. ( 4) in kapellasite, by noting that in the paramagnetic regime the Nagata's theory predicts a linear scaling of the line shift with susceptibility [37,39]

g z (T ) = S z 2μ B B 0 j =i 2K zz ij -K xx ij -K yy ij , (5)
where S z = i S z i = χ mol (T )B 0 N A gμ 0 μ B , z denotes the orientation of the applied field, and x,y the two perpendicular directions, N A is the Avogadro number, and μ 0 the vacuum permeability. A convincing scaling is observed in kapellasite for both directions in the whole investigated temperature range [Fig. 6(d)], an additional proof that the increase of the bulk susceptibility is intrinsic even at low temperatures. The scaling yields the limiting infinite-temperature values of the g factor within the kagome planes g 90 • ∞ = 2.140 and along the c axis g 0 • ∞ = 2.115. These values correspond to the values averaged over the structural triangular unit.

The contribution to the above sum related to the dipolar field arising from the finite spin S z in a polarized paramagnet can be calculated exactly,

B dd = S z μ B j =i (2d zz ij -d xx ij -d yy ij
). Following Lorentz, we write [40]

B dd = B dd + B Lor + B dem ,
where the first term is a sum of dipolar fields within a (Lorentz) sphere with a radius large enough to provide convergence, the second term B Lor = 1 3 χB 0 is the Lorentz field that originates from the polarization of the sphere's surface, and the third term B dem = -NχB 0 is the demagnetization field depending on the shape of the sample. By distributing Cu 2+ spins randomly on both Cu and Zn sites according to the sample stoichiometry, we find B 0 • dd /χ mol = -1.80 T(emu/mol Cu) -1 and B 90 • dd /χ mol = 0.90 T(emu/mol Cu) -1 . Further, based on the oblate shape of the crystalline grains (2r/d ∼ 10), we estimate N 0 • = 0.80 (5) and N 90 • = 0.10(3) [41], which yields (B Lor + B dem ) 0 • /χ mol = -1.61 T(emu/mol Cu) -1 and (B Lor + B dem ) 90 • /χ mol = 0.81 T(emu/mol Cu) -1 .

The effect of the dipolar fields (including demagnetization) is too small to explain large experimental shifts [see Fig. 6(a) for a comparison of corrected and noncorrected data]. We subtract this contribution and derive "pure" g shifts [Fig. 6(c)] arising from the symmetric anisotropic exchange. The linear scaling of the shift with macroscopic susceptibility [Fig. 6(d)] can be used to evaluate the magnitude of the exchange anisotropy. Since the Cu-Cl bonds defining the principal axis of the uniaxial local g tensor on each site are tilted from the crystallographic c direction towards the center of the triangle by 43 • , the direction of the principal axis of D ij of each ij bond is not trivial. Therefore, we only inspect the shift data for θ = 0 • , for which employing the threefold-rotational symmetry of the crystal structure, we derive

g 0 • (T ) = χ mol (T ) N A gμ 0 μ 2 B 4D cc , ( 6
)
where D cc = D cc 1 + 1 2 D cc d corresponds to an effective component of the magnetic anisotropy arising from both dominant bonds J 1 and J d .

From the slope g 0 • (T )/χ mol (T ) = -2.0 (emu/mol Cu) -1 we obtain D cc = -0.4 K, yielding the ratio |D cc /J 1 | = 3% with respect to the dominant exchange. This is somewhat larger than the conventional Moriya estimate |D/J | ∼ ( g/g) 2 ∼ 2%. However, D cc sums two contributions and, moreover, similarly enhanced values of the symmetric anisotropic exchange interaction were reported before in copper-based systems [32,34]. Our results thus suggest that kapellasite is quite close to the isotropic Heisenberg limit. We stress that the DM interaction, which is generally larger than the symmetric exchange anisotropy in copper-based systems (of the order g/g ∼ 15%) could not be assessed in this study, although it may well be present also in kapellasite. Last, we note that the large line shifts in kapellasite observed at low temperatures are due to rather high values of the intrinsic magnetic susceptibility. In contrast, in herbertsmithite and vesignieite, a suppressed susceptibility did not allow an unambiguous assignment of the symmetric exchange anisotropy based on the ESR line shift.

IV. SPECIFIC HEAT

The specific heat C(T ) was measured with a Quantum Design Physical Property Measurement System (PPMS) from 300 to 0.4 K on thin pressed pellets of kapellasite with a typical mass of 8 mg, in magnetic fields up to 9 T (Fig. 7). The analysis of C(T ) under zero magnetic field has been the subject of two recent publications, where its T dependence was used to refine the model Hamiltonian of kapellasite, based on high-temperature series calculations [20,22]. In the following, we present its field dependence and confirm that its low-T part is magnetic and intrinsic, and thus provides a pertinent quantity for extracting the exchange integrals.

In the absence of a nonmagnetic compound isostructural to kapellasite, the phonon contribution is estimated at high temperatures (50 T 160 K) with a phenomenological fit made of a linear combination of Debye and Einstein functions. The lattice contribution is found to be negligible below T = 10 K (see Fig. 7). Therefore, the broad peak singled out in C(T )/T around 2 K has clearly a magnetic origin.

The specific heat is found to be insensitive to fields smaller than 3 T, however, larger fields cause a shift of the peak towards higher T , up to 4 K for μ 0 H = 9 T, as well as a decrease of intensity. This field dependence definitely shows the magnetic character of C(T ) at low T . The peak that is observed cannot be attributed to the presence of orphan spins which would contribute as a Schottky anomaly in the specific heat as observed, for example, in herbertsmithite [42]. As well as displaying a peak, such a contribution would have a linear dependence in field, contrary to our observations (Fig. 7, inset). This is completely consistent with the fact that the upturn of the susceptibility is not driven by weakly coupled spins but is intrinsic as clearly shown by the perfect match of the measured NMR and ESR shifts with the macroscopic susceptibility. Therefore, in contrast with other frustrated systems with larger magnetic couplings, the field dependence of C(T ) observed here is due to the comparable energy scale between the applied magnetic field and the rather low magnetic exchanges in kapellasite. As a main consequence, the specific heat for μ 0 H = 0 T is characteristic of the model Hamiltonian of kapellasite and is relevant for the evaluation of the exchange interactions as performed by Bernu et al. [22].

The T dependence of the magnetic part of C(T ) could be reproduced in the range 0.4-10 K within a J 1 -J 2 -J d model, provided one assumes a ∼14% missing entropy, in this T range. In fact, the total magnetic entropy expected for a spin S = 1 2 , S = R ln 2, is fully recovered only by T = 40 K, a quite extended regime when compared to the former model. This release of the missing entropy at higher T appears in good agreement with the extended T regime where inelastic neutron scattering studies reveal the presence of magnetic correlations [20]. Two different hypotheses were proposed in Ref. [22] to account for this missing entropy based on the presence of spin vacancies: the existence of weakly coupled spins or the formation of spin singlets. The former scenario can be ruled out, as already discussed, since it would give rise to a linear field dependence of the maximum of the specific heat below 10 K. The second scenario, namely spin-singlet freezing, remains, on the contrary, plausible since the exchange along the diagonal J d ∼ 15.6 K is antiferromagnetic in nature. The identification of the involved spins is still unclear, although the singlet could result from the coupling of a a spin sitting on a kagome site and a spin defect on the center of a hexagon. Their contribution would have to show up at higher temperatures (T ∼ 10 K), and would therefore remain almost unaffected by fields μ 0 H < 9 T (∼12 K). While we are not able to single out the exact behavior of these defects from these measurements, this last scenario is consistent with our experimental results. It is not clear, however, if a spin singlet can be stabilized in the presence of other large competing interactions.

V. SPIN DYNAMICS

A. 1/T 1 NMR spin-lattice relaxation rate

The 35 Cl spin-lattice relaxation rate 1/T 1 was determined through the saturation recovery method, with a standard π/2tπ/2τπ pulse sequence (70 < τ < 100μs), from 40 to 1.3 K. The measurements were performed at the center of lines 1, 2, and 3, which gave similar average values 2 is successfully applied at high T (blue line) but no longer valid at 1.3 K where a stretched exponential curve is more appropriate (red line), due to a distribution of relaxation times τ . (c) Temperature dependence of the spin-lattice relaxation rate measured in μSR (1/T 1 , right axis, reproduced from [20]), in NMR (1/τ , left axis), and calculated from inelastic neutron scattering (INS) (see text). The color map represents the intensity of the probability density function ρ(τ ) (left axis) extracted from NMR measurements (gray hatched area is for the absence of NMR data below 1.2 K). Below 10 K, the distribution of τ is too large to allow a single relaxation time to be defined. and distributions (see following) of T 1 at low temperature. The evolution of T 1 with temperature was then extracted from measurements performed on the line with the highest intensity (1).

Within this T range, the spin dynamics can be divided into two regimes described below: (i) a high-T paramagnetic one where T 1 remains constant and (ii) a low-T fluctuating regime where a broad distribution of relaxation times is observed.

(i) High-T paramagnetic relaxation. Above 10 K, the magnetization-recovery curve M(t) is consistent with that expected for a relaxation process of magnetic origin for a nuclear spin I = 3 2 with a single spin-lattice relaxation time T 1 [43]:

M(t,T 1 ) = M sat (1 -0.1e -t/T 1 -0.9e -6t/T 1 ). ( 7
)
Fitting our recovery curves to Eq. ( 7) yields the T 1 values plotted in Fig. 8(c) with a constant T 1 = 12(3) ms value for 10 < T < 40 K.

In conventional insulating paramagnets, for T J , the spin-lattice relaxation rate 1/T 1 is indeed temperature independent. The latter was explicitly calculated by Moriya [44] to be

1/T 1 = √ 2πγ 2 N g 2 A 2 hf S(S + 1)/3z 1 ω e , ( 8 )
where z 1 = 3 is the number of Cu 2+ probed and ω e = J √ 2zS(S + 1)/3/ is the exchange frequency with z = 4 being the number of Cu 2+ neighbors on the kagome lattice.

Using experimental values A hf = -4.98 kOe/μ B for the 3-Cu triangle and the average exchange J = θ CW = 9.5 K, the formula leads to T 1 = 11 ms, which is in excellent agreement with the measured T 1 [46].

(ii) Low-T relaxation. At lower T , in the regime where correlations become noticeable, one has to take into account the q dependence of the excitations. Indeed, T 1 is given by

1 T 1 = γ 2 μ 2 B k B T q |A(q)| 2 χ ⊥ (q,ω 0 ) ω 0 , (9)
where A(q) acts as a filtering factor for some specific values of q, associated here with the triangular geometry. The functions |A(q)| 2 are plotted in Fig. 9 for the various Cl sites. It is notable that at the peak of S(q) detected in neutron experiments [20], the filtering effect varies by only ∼30% from one site to the other, with a minimal transmission of ∼60%. In the following, we give our results for the Cl site (1) situated at the center of a fully occupied triangle.

If the conventional exponential shape of relaxation is indeed observed above 10 K, the situation is remarkably different at T = 1.3 K. The saturation recovery curve there had to be fitted using the phenomenological stretched exponential variant of Eq. ( 7) [Fig. 8(b)]:

M(t) = M sat (1 -0.1e -(t/τ 0 ) β -0.9e -(6t/τ 0 ) β ). (10)
FIG. 9. (Color online) Intensity of the normalized form factor |A(q x ,q y )| 2 related to sites 1 (left) and 2 (right), shown as grayscale maps. White (black) corresponds to a full (no) transmission. The site 3 does not show any filtering effect. Blue (red) crosses indicate the wave vectors of the cuboc2 (q = 0) state. This points to a broad distribution of relaxation times τ , with a probability density function ρ(τ ), so that

M(t) = +∞ 0 M(t,τ )ρ(τ )dτ. (11)
Physically, τ 0 is related to the mean of the probability density function while β characterizes its dispersion. The probability density function itself can be derived from an inverse Laplace transform of Eq. ( 11), whose explicit expression can be found elsewhere [45], and which, following Ref. [47], can be numerically evaluated once β and τ 0 are known. In order to follow the evolution of the distribution of relaxation times with temperature, we fitted the magnetization recovery curve with Eq. ( 10) for each measured temperature. The pair of parameters (τ 0 ,β) extracted from these fits allows then to compute the ρ(τ ) profiles for 1.3 T 40 K [Fig. 8(a)]. Within the investigated T range, β is found to change smoothly from 1.0(1) at T = 40 K to 0.54(4) at T = 1.3 K. The results of these fits are presented in the color map of Fig. 8(c). Above 10 K, the distribution ρ(τ ) is well peaked around τ 0 , with ρ(τ ) ∼ δ(ττ 0 ), which legitimately defines a single relaxation time T 1 = τ 0 ∼ 12 ms, with β → 1. It then asymmetrically broadens at lower T , exhibiting a τ → 0 divergence when T → 0, while β approaches smoothly its lowest value of ∼0.5. The same qualitative behavior was observed for all Cl sites.

B. LF μSR relaxation

Complementary information on the spin dynamics can be obtained from μSR measurements in longitudinal field since it directly probes the spectral fluctuation density. As discussed in Ref. [20], the μ + polarization is described by P (t) = f e -t/T 1 + (1f ), where 1f = 0.14 is a slowly relaxing muons fraction completely decoupled for μ 0 H > 100G and thus responsible for the nonzero value observed at long time [20] (Fig. 10). The application of a small longitudinal field μ 0 H = 100 G also enables one to negate the static nuclear contribution and allows the electronic fluctuations to be probed separately. Similarly to NMR, a constant relaxation rate 1/T μ 1 is observed for T > 5 K [Fig. 8(c)]. Assuming local fluctuations, the spin-lattice relaxation times are linked through the relation 2 . At T = 20 K, the experimental scaling factor is found to be T 1 /T μ 1 ∼ 150, leading to the muon coupling constant A μ = 1.88 kOe/μ B , close to the dipolar coupling estimated to be A dip μ ∼ 1.3 kOe/μ B at the principal OH -muon site. Thus, the spin fluctuations probed in μSR above 10 K also originate from a paramagnetic regime.

T 1 /T μ 1 = (γ μ A μ /γ n A hf )
A good agreement is observed for the 3-K onset T where the relaxation rate departs from the paramagnetic limit while the evolution at lower T follows qualitatively that of the probability density function profiles extracted from 35 Cl NMR, as underlined by Fig. 8(c).

However, contrary to NMR, a broad distribution of relaxation times does not appear clearly in μSR experiments. Should the same distribution be observed, the μ + polarization for μ 0 H = 100 G would vary as P (t) = e -(t/T 1 ) β with β ∼ 0.5 at low T . Allowing a stretched exponential relaxation does not appreciably improve the quality of the fit. The β exponent is ill defined but its value is found to remain above 0.7 in the whole T range below 4 K which would suggest a moderate distribution of relaxation times as probed by the muon.

Lastly, we comment on the field dependence of the relaxation rate. In disordered spin systems with a single relaxation time and exponentially decaying correlations, the spin-lattice relaxation time T 1 exhibits a field square dependence T 1 ∝ H 2 , characteristic of a Lorentzian spectral density [48]. Such a spectral-weight function, which is centered at zero energy, is commonly used to describe conventional fluctuations in the paramagnetic regime. Here, a power-law dependence T 1 ∝ H α is observed at 60 mK, with α 1 (Fig. 10,inset). This observation is inconsistent with the existence of a single time scale, as underlined by our NMR measurements, and suggests a more exotic spectral density, such as the one at play in the spin-liquid phase of the S = 1 2 antiferromagnetic chain [49]. The contrast between the low-T relaxation extracted from the two NMR and μSR techniques is discussed in the following section.

VI. DISCUSSION

A. Effects of disorder

In this section, we discuss the information that can be extracted from the NMR spectra and the related impact of defects, that is, deviations from the ideal kagome lattice, on the ground state of kapellasite. Let us first recall that the structural refinements yield the chemical formula (Cu 0.73 Zn 0.27 ) 3 (Zn 0.88 Cu 0.12 )(OH) 6 Cl 2 which points to both the dilution of the kagome lattice (27%) and a substantial occupation of the Zn site by Cu (12%).

Spin vacancies are rather common in experimental materials, e.g., in the kagome antiferromagnet where, for a NN antiferromagnetic Hamiltonian, they are known to create a spin singlet on the sites first neighbors of the impurity and induce a staggered susceptibility at larger distances [50][51][52]. In the case of kapellasite, where the scheme of interactions is different, little is known about such an impact of the spin vacancy. From a general viewpoint, magnetic defects are expected to generate a distribution of susceptibilities. In addition, the rapid crystallization of kapellasite during the synthesis process [9] may affect the homogeneity of the defect distribution, and could further influence the disorder. The latter can yield a distribution of hyperfine constants. Although we expect the distribution of susceptibilities to dominate the linewidth of our spectra H ∼ (A hf χ ), especially in the low-T range, attributing all the width to a structural disorder-induced distribution of hyperfine constants would lead to a scaling H ∼ χ . An upper bound of A hf /A hf ∼ 33% can then be extracted from the linear slope of the plot of H versus χ in Fig. 11. In the low-T range, we note that, for site 2, the linear scaling is severely disrupted below 10 K which rather argues in favor of a distribution of susceptibilities, as expected.

In the context of a large dilution, one could wonder whether the absence of any spin freezing really reflects the spin-liquid physics of the ideal J 1 -J 2 -J d Heisenberg Hamiltonian. We first recall that the level of dilution of d = 27% is still above the percolation threshold of the kagome lattice p c 0.652 [23,24], i.e., that d < d c = 1p c 0.35. On the contrary, delafossite LaCuO 2.66 , with a 33% ( d c ) diluted kagome lattice, displays an inhomogeneous spin-frozen ground state [53]. While the extreme two-dimensional character of kapellasite could also be invoked to support the absence of freezing at finite temperature, the Mg analog of kapellasite, haydeite, has a clear transition at T ∼ 4 K [9,54]. This gives support to the intrinsic character of the spin-liquid behavior which is found in kapellasite.

To date, the existence of such a spin-liquid phase is supported within two theoretical approaches: (i) a Schwingerboson mean-field model by Messio et al., which reproduces most of the experimental observations [20], including the evolution of the μSR relaxation rate with temperature, but inherently fails to reproduce a gapless behavior as observed both in neutron and our experiments; and (ii) a pseudofermion functional renormalization group approach by Suttner et al. [55]. There, kapellasite would locate in a small nonfrozen region between the cuboc and the ferromagnetic phase, in a slightly different J 1 -J 2 Heisenberg model.

B. Inhomogeneous dynamics

We now turn to the spin dynamics in the low-temperature phase and discuss the intriguing difference observed between the relaxation probed by NMR and μSR. The form factor |A(q)| 2 in Eq. ( 9) could be a potential source of disagreement between the two measured relaxations. The muon most likely sits in a nonsymmetric site and is dominantly coupled to Cu 2+ via dipolar interaction while the 35 Cl nucleus, located at a symmetric site above a Cu 2+ triangle, probes each Cu at the corners of this triangle through hyperfine couplings.

While the distribution thought to be of hyperfine couplings mentioned previously could be responsible for the distribution of T 1 seen only in NMR, this can be ruled out as the broadening of the distribution function is clearly temperature dependent [Fig. 8(a)]; a distribution of hyperfine constants remains T independent.

The q dependence of the form factor may also filter out some modes in NMR as compared to μSR. Let us recall that the different sites are clearly resolved in the NMR spectra, and that the same stretching of the relaxation is observed at each site at low temperature. Since the intensity of the form factors differs from one site to another (see Fig. 9), it is unlikely that the q dependence of |A(q)| 2 is responsible for the T 1 distribution.

We are thus left with the firm conclusion that the distribution of T 1 observed in NMR mainly comes from the fluctuation spectrum. The difference between NMR and μSR relaxations could be ascribed to the difference of these couplings rather than being purely explained by the different time scale probed with each technique. Indeed, for both techniques, T -1 1 scales with A 2 (A is the coupling constant) so that the long-time tail of the T 1 distribution detected through NMR would yield unobservable relaxation times in μSR. Cutting off this tail by restricting the integration range of the probability density function extracted from NMR over [0,α] with α < +∞ yields a distribution of 1/T μ 1 = λ [see Eq. ( 11)] still too large to 205103-10 account for the smaller breadth of relaxation times (β 0.7) observed in μSR. It is therefore necessary to invoke very local modes that give distinct relaxations characteristics to the various sites probed by NMR and which have a distribution that is smoothed out through the long-range character of the dipolar coupling between the muon probe and the fluctuating Cu 2+ spins.

One of the characteristic features of the low-T phase of kapellasite is the well-localized intensity present in the inelastic channel measured by inelastic neutron scattering (INS) at the M point of the Brillouin zone. Both INS and NMR relaxation probe the dynamical susceptibility χ (q,ω), yet NMR has no q resolution (except via the form factor) and probes the ω → 0 region. In order to extract some information from their comparison, we convert χ (q,ω) measured by INS in an equivalent 1/T 1 by integrating χ in the ranges q = [0.2,0.8] Å and ω = [0.4,0.8] meV. The values obtained are scaled so as to match the high-T NMR and μSR ones. The result is plotted on Fig. 8(c). We note that no dramatic filtering effect is expected when applying the NMR form factor |A(q)| 2 of different sites, in agreement with the nonzero value of the form factors at the cuboc2 q positions (Fig. 9). The T evolution of T 1 generated from the INS data is in good agreement with the flat dependence observed by local probes above 10 K, but fails to reproduce any distribution at lower T . Therefore, we can conclude that at low T the NMR relaxation is dominated by extra low-energy modes that are inaccessible to INS. The low-energy fluctuation spectrum would then have spectral weight relative not only to the cuboc2 mode, but also at many different q positions, not filtered out by the form factors.

It is interesting to note that the cuboc2 phase is chiral and consequently could allow for the observation of a chiral phase transition at finite temperature, studied theoretically in details in Refs. [56,57]. This transition is combined with a proliferation of topological point defects (Z 2 vortices), which preferentially nucleate around the walls separating domains with opposite chirality. Such excitations were previously suspected to be at the origin of the unconventional spin dynamics observed in the Heisenberg triangular antiferromagnet NaCrO 2 [58]. Further work would be required to confirm that the peculiar low-energy excitations observed in kapellasite are indeed a signature of such an exotic transition.

VII. CONCLUSION

We reported the experimental investigation of the S = 1 2 Heisenberg kagome system kapellasite, where the two main competing interactions J 1 and J d conspire to stabilize a spinliquid state at low T with spin correlations of the cuboc2 type. 35 Cl NMR data have confirmed the previously established level of dilution of 27% of the kagome lattice and further evidenced its random character. Surprisingly, an identical local magnetic susceptibility is probed at high temperature for each Cu 2+ spin, irrespective of its neighborhood. This is a strong evidence of the homogeneity of the physics in the system, despite the substantial level of disorder. Additional improvements in our understanding of the disorder effect would require a set of samples with various dilution, which remains at the moment a real chemical challenge, due to the metastability of the kapellasite phase. ESR showed that the Heisenberg model is relevant to kapellasite and established a moderate magnetic anisotropy of the symmetric anisotropic exchange type of |D cc /J 1 | = 3% while no conclusion can be given about the Dzyaloshinskii-Moriya anisotropy. The intrinsic character of the specific heat and of the magnetic susceptibility are thoroughly assessed which gives a solid experimental credit to the theoretical procedure applied by Bernu et al. [22] to extract the exchange couplings of kapellasite from high-temperature series fits. Unconventional spin dynamics is consistently revealed by NMR, μSR, and inelastic neutron scattering [20] upon cooling below T = 5 K, where spin correlations effectively build up. The concomitant distribution of relaxation observed in NMR is then due to the emergence of extra low-energy modes, inaccessible to the neutron time scale. It could either be a complex signature of the effect of disorder, or of the specific excitation spectrum of the low-T cuboc2 chiral phase.

In conclusion, while a ferromagnetic kagome lattice is not frustrated, the existence of a further nearest-neighbor exchange along the diagonal J d in kapellasite, comparable to |J 1 |, restores magnetic frustration via the competition of interactions. The effect of frustration through the geometry of interactions was only relatively recently explored on the kagome lattice [55]. Up to now, most of the work on kagome magnets has been dedicated to the pure antiferromagnet Heisenberg model and the consideration of the effects of additional energy scales is still in its infancy. Second-neighbor interactions have also been argued to be necessary to explain the neutron data on herbertsmithite single crystals [55]. Kapellasite gives an order of magnitude for such interactions and more broadly opens the way towards the experimental study of the phase diagram of the J 1 -J 2 -J d Heisenberg model which displays noncoplanar chiral ordered phases and nonmagnetic regions. The discovery of its spin-liquid state certainly deserves further attention in order to reach a more comprehensive understanding, and provides new exciting perspectives to explore quantum spin-liquid physics on the kagome lattice.

  (b)]. Neutron diffraction

FIG. 1 .

 1 FIG. 1. (Color online) (a) The kagome magnetic pattern arises from the regular arrangement of Cu 2+ (S = 1 2 ) within a Cu/Zn triangular lattice (Cu in blue, Zn in yellow). The main exchange couplings J 1 , J 2 , and J d are presented on the diagram. (b) Scanning electron microscope image of kapellasite. The width of the image is 20 μm. The flat shape of the micron scale grains evidences the bidimensional character of kapellasite. (c) Octahedral coordination environments of Cu and Zn. (d) The kagome layers are only weakly coupled through Cl-OH hydrogen bonds. The 35 Cl NMR probe is dominantly coupled to three Cu. Zn atoms are not shown for clarity.

FIG. 4 .

 4 FIG. 4. (Color online) (a) 35 Cl NMR spectra of the central lines plotted versus the line shift from high (right) to low (left) temperature. The spectra are vertically shifted for clarity. The composition of the spectrum in terms of different Cl sites is shown at 80 K. The local susceptibilites can be directly extracted from the shift of the lines from the reference (dashed vertical red lines). As an example, the blue circles emphasize the shift of the NMR line corresponding to site 1 (left panel). (b) 35 Cl NMR shifts for the different Cl sites vs χ macro at high T . The hyperfine couplings are extracted from the slope of the curves [A hf,1 = -4.98(3) kOe/μ B ]. (c) Temperature dependence of the macroscopic susceptibility χ macro , the local NMR susceptibility χ spin (left) and of 1/χ macro , 1/χ spin (right). The Curie-Weiss parameter θ CW = 9.5 ± 1 K evidences a weak ferromagnetic effective coupling at high T .

FIG. 5 .

 5 FIG. 5. (Color online) Ratio of local susceptibilities χ 1 /χ 2 (black crosses) vs T . The red line is a fit using the ratio (Tθ 2 )/(Tθ 1 ) for 20 T 70 K, with θ 1 = 7.8 K and θ 2 = 5 K (see text). The expected value of 1 is found above 70 K while the steep increase at low T underlines the strengthening of antiferromagnetic correlations at the diluted site that is expected when J 1 is ferromagnetic.

FIG. 6 .

 6 FIG. 6. (Color online) Temperature dependence of (a) the ESR line position, (b) linewidth, and (c) g factor measured at 326.4 GHz. (d) Scaling of the g factor with bulk susceptibility. Open symbols include all magnetic anisotropy contributions present in kapellasite, while solid symbols represent data corrected for the dipolar (including demagnetization) fields.

FIG. 7 .

 7 FIG. 7. (Color online) Molar specific heat divided by temperature C/T per formula unit vs temperature and under magnetic field from 0 to 9 T. Below 10 K, the phonon contribution C ph is negligible. Inset: The field dependence of the peak position (T m ) is nonlinear, contrary to a Schottky behavior.

FIG. 8 .

 8 FIG. 8. (Color online) (a) Probability density functions ρ(τ ) profiles for selected temperatures. (b) Magnetization-recovery curve at T = 30 and 1.3 K. The conventional relaxation form for I = 32 is successfully applied at high T (blue line) but no longer valid at 1.3 K where a stretched exponential curve is more appropriate (red line), due to a distribution of relaxation times τ . (c) Temperature dependence of the spin-lattice relaxation rate measured in μSR (1/T 1 , right axis, reproduced from[20]), in NMR (1/τ , left axis), and calculated from inelastic neutron scattering (INS) (see text). The color map represents the intensity of the probability density function ρ(τ ) (left axis) extracted from NMR measurements (gray hatched area is for the absence of NMR data below 1.2 K). Below 10 K, the distribution of τ is too large to allow a single relaxation time to be defined.

FIG. 10 .

 10 FIG. 10. (Color online) Muon polarization under longitudinal applied magnetic fields μ 0 H = 0-2500 G, at T = 60 mK. Inset: H dependence of the spin-lattice relaxation time T 1 .

FIG. 11 .

 11 FIG. 11. (Color online) FWHM of the NMR lines corresponding to sites 1 and 2, from the spectra displayed on Fig. 4. Temperature is an implicit parameter.

TABLE I

 I

			FIG. 3. (Color online) (a)
	. Cl sites occupancies evaluated by NMR lines integrals
	A i and compared to binomial probabilities assuming a random
	disorder for a kagome site filling ratio p = 0.75.
	Sites i/config.	Int. A i	Bin. prob. (p = 0.75)
	1 Cu 3	0.43(1)	0.42
	2 Cu 2 Zn	0.42(1)	0.42
	3 CuZn 2	0.10(5)	0.14
	4 Zn 3	0.05(4)	0.02

ACKNOWLEDGMENTS E.K. acknowledges very useful discussions with B. Bernu. This work was supported in part by the European Commission under the 6th Framework Programme Contract No. RII3-CT-2003-505925, by the French Agence Nationale de la Recherche under Grants No. ANR-09-JCJC-0093-01, No. ANR-SPINLIQ-86998, No. ANR-HFM-86998, by the PHC Proteus Program No. 24322SF and by Université Paris-Sud, Grant MRM PMP. A.Z. acknowledges the financial support of the Slovenian Research Agency, Project No. BI-US/14-45-039. The NHMFL is supported by NSF Cooperative Agreement No. DMR-1157490 and by the state government of Florida.

APPENDIX: 35 Cl NMR MEASUREMENTS

Oriented sample and experiments

NMR experiments were performed on an oriented powder of kapellasite, aligned along the c axis. This sample was prepared by mixing together 200 mg of kapellasite powder, prepared following Ref. [19], with a stycast epoxy. The mixture was then placed in a cylindrical teflon mold and cured under a magnetic field of 7 T for 15 h. The level of orientation was estimated to be ∼90% from the comparison of intensities between the nonoriented and oriented NMR spectrum. Therefore, the oriented sample behaves like a single crystal in the c-axis direction, i.e., perpendicular to the kagome planes. This orientation allows the observation of narrow NMR lines instead of conventional complex powder line shape due to quadrupolar effects.

NMR experiments were performed on the most abundant Cl isotope 35 Cl (abundance of 75.8%), with a nuclear spin I = 3 2 and a gyromagnetic ratio γ n = 4.1716 MHz/T. There is only one Cl crystallographic site within the structure, which is coupled to three Cu 2+ ions of one triangle belonging to the kagome plane [Fig. 1(d)]. NMR measurements were made both in fixed and sweep field configurations, under various magnetic fields from 3 to 13 T. A 4 He flow cryostat was used to reach temperatures from 1.2 to 300 K. 35 Cl is a quadrupolar nucleus (I = 3

Nuclear Hamiltonian

2 ), hence, one has to consider two sets of interactions resulting from the coupling of the nucleus to its magnetic and charge environments. The nuclear Hamiltonian can be written [START_REF] Slichter | Principles of Magnetic Resonance[END_REF]26]

with

The first term of Eq. (A2) corresponds to the nuclear Zeeman effect under a static applied field H 0 . The second term accounts for the magnetic interaction of the nuclear spin with its local environment through the shift tensor K. Its T -dependent part directly probes the spin susceptibility of the surrounding magnetic moments. The quadrupolar Hamiltonian [Eq. (A3)] is characterized by two parameters: ν q = 3e 2 qQ/2I (2I -1)h, which results from the coupling of the nuclear quadrupole moment Q with the electrostatic field gradient (EFG) characterized by its q value along its principal z axis, and its asymmetry η in the (x ,y ) plane (0 η 1).

For an "ideal" structure of kapellasite, i.e., with perfect stoichiometry and no local distortion which might result from the Zn/Cu substitution, a threefold axial symmetry exists at the Cl site. This simplifies considerably both the hyperfine tensor and the quadrupolar Hamiltonian which then have the same principal axis along z and a rotational invariance in the (x, y) plane, hence, K x = K y and η = 0.

Under the applied fields of a few teslas, H quad is a small perturbation as compared to H mag for the experimental case considered in this work. In first order, the (2I + 1) energy levels E m write E m = m|H|m , where -I m I and |m is the eigenstate of the Zeeman Hamiltonian for an applied field H 0 . This energy splitting gives rise to the 2I resonance frequencies ν m-1,m = (E m-1 -E m )/h detected in NMR. Taking into account the previous symmetry arguments and expressing H 0 in spherical coordinates, H 0 = H 0 (sin θ cos φ, sin θ sin φ, cos θ ) x,y,z , we obtain for the "ideal" case

where ν 0 is the resonance frequency of the bare nucleus γ n H 0 corrected by the usual T -independent chemical shift. For 35 Cl-NMR (I = 3 2 ), one obtains three resonances: one with no quadrupolar shift to the first-order approximation and is only governed by magnetic effects (the central line at ν -1/2,1/2 ), and two satellites ν 1/2,3/2 and ν -3/2,-1/2 . With the applied field along the z axis (θ = 0 • ) and in the case of axial symmetry (η = 0), we simply derive

Similar formulas can be obtained for the applied field perpendicular to the z axis:

Assuming as a first approximation that the susceptibility is isotropic, we can directly extract the variation of the shift when the field is applied along the z direction K z , from the variation of the susceptibility through δK z (T ) = A hf,z δχ spin (T ) where A hf,z is the hyperfine coupling for H 0 z. As a final note to this section, we underline that any local structural distortion will introduce deviations from these simple formulas (η = 0 and different principal axis for the EFG and shift tensor not aligned along z).