Felipe García-Sánchez

Chuanzeng

Chuanzeng Zhang
email: c.zhang@uni-siegen.de

A comparative study

Keywords: Hypersingular traction BEM, Frequency-domain BEM, Laplace-domain BEM, Time-domain BEM

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

Introduction

Dynamic crack analysis has many important applications in engineering sciences such as in fracture and damage mechanics, quantitative non-destructive material testing, geophysics and geomechanics. Since analytical solutions for dynamic crack problems can be obtained only for very simple cases, numerical methods have to be applied in general to solve the arising initial boundary-value problems. Among many available numerical methods, the boundary element method (BEM) is an accurate and efficient numerical tool for dynamic crack analysis, at least for homogeneous and linear elastic solids.

From the mathematical points of view, a crack in twodimensional (2-D) elastic solids is a line with two coincident faces, which leads to a degeneration of the classical displacement BEM formulation over both crack-faces. This problem can be avoided by using the dual BEM, where the displacement boundary integral equations (DBIEs) are used over one of the crack-faces while the traction boundary integral equations (TBIEs) are applied to other crack-face. TBIEs can be obtained by the partial differentiation of the DBIEs and the subsequent application of the Hooke's law. Another remedy to overcome the degeneration of the DBIEs for crack analysis is the use of the hypersingular TBIEs only on one of the crack-faces, where the crack-openingdisplacements (CODs) are the fundamental unknown quantities. For cracked solids with finite boundaries, a combined use of both the DBIEs and the TBIEs is promising. In this case, the TBIEs are applied to one of the crack-faces while the DBIEs are used on the external boundaries of the cracked solids. An overview of different possibilities for crack analysis in cracked elastic solids with bounded domains by BEM can be found in [START_REF] Ch | A 2-D hypersingular time-domain traction BEM for transient elastodynamic crack analysis[END_REF]. The numerically computed CODs can be used in a data post-processing to obtain the stress intensity factors (SIFs), which are the most important crack-tip characterizing parameters in linear elastic fracture mechanics.

Over the past years, three different BEM formulations, namely, the frequency-domain [START_REF] Dineva | BIEM for 2D steadystate problems in cracked anisotropic materials[END_REF][START_REF] García-Sánchez | Two-dimensional time-harmonic BEM for cracked anisotropic solids[END_REF], the Laplacedomain [START_REF] Albuquerque | Dual reciprocity boundary element method in Laplace domain applied to anisotropic dynamic problems[END_REF] and the time-domain [START_REF] Albuquerque | The boundary element method applied to time dependent problems in anisotropic materials[END_REF][START_REF] Tan | A time-domain BEM for transient wave scattering analysis by a crack in anisotropic solids[END_REF][START_REF] Tan | A time-domain collocation-Galerkin BEM for transient dynamic crack analysis in anisotropic solids[END_REF][START_REF] Ch | Transient elastodynamic antiplane crack analysis of anisotropic solids[END_REF][START_REF] Ch | Transient dynamic crack analysis in anisotropic solids[END_REF] BEM, are often applied to transient elastodynamic crack analysis. To analyze their accuracy and efficiency, a comparative study is performed in this paper. Hypersingular TBIEs are applied for this purpose. A collocation method for the spatial discretization of the hypersingular BIEs is adopted. Hypersingular integrals are dealt with by a regularization technique based on a suitable variable change [START_REF] García | Traction boundary elements for cracks in anisotropic solids[END_REF][START_REF] García-Sánchez | Two-dimensional time-harmonic BEM for cracked anisotropic solids[END_REF]. In the frequency-domain and the Laplace-domain BEM, hypersingular BIEs in the transformed domain are first solved numerically for discrete frequency and Laplace-transform parameters. To obtain the time-dependent solutions, fast Fourier inverse transform and Stehfest's Laplace-inversion algorithm [START_REF] Stehfest | Comun. Algorithm 368: numerical inversion of Laplace transform; an overview and recent developments[END_REF] are applied in the frequency-domain and the Laplacedomain BEM, respectively. The time-domain BEM uses a convolution quadrature formula of Lubich [START_REF] Lubich | Convolution quadrature and discretized operational calculus, Part I[END_REF][START_REF] Lubich | Convolution quadrature and discretized operational calculus, Part II[END_REF] for approximating the arising convolution integrals and it leads directly to time-dependent solutions. Numerical examples for computing transient elastodynamic SIFs in homogeneous and linear elastic solids of general anisotropy are presented to compare the accuracy and the efficiency of the three different BEM formulations.

Problem statement and elastodynamic BIEs

We consider a finite crack in a 2-D, homogeneous, linear elastic and anisotropic solid. The deformation of the cracked solid is assumed to be in a state of either generalized plane strain or generalized plane stress. In the absence of body forces, the cracked anisotropic solid satisfies the equations of motion [START_REF] Achenbach | Wave propagation in elastic solids[END_REF]

σ ij,j = ρ üi , ( 1 )
Hooke's law

σ ij = C ijkl u k,l , ( 2 )
the initial conditions

u i (x, t) = ui (x, t) = 0, for t = 0, (3)
and the traction boundary conditions on the crack-faces

p i (x, t) = σ ij (x, t)n j (x) = p * i (x, t), x ∈ c . ( 4 )
In Eqs. ( 1)-( 4), σ ij and u i denote the stress and the displacement components, C ijkl is the fourth order elasticity tensor, p i (x, t) is the traction vector, n j is the outward unit normal vector, c = + c + - c represents the upper and the lower crack-faces, ρ is the mass density, and p * i (x, t) is the prescribed crack-face loading, respectively. Also, a comma after a quantity represents partial derivatives with respect to spatial variables, and superscript dots stand for the time differentiations of the quantity. Unless otherwise stated, the conventional summation rule over double indices is implied, and the indices i and j take the values 1 and 2.

The displacements can be represented by the following boundary integral

u i (ξ , t) = - + c p * ij (ξ , x, t) * u j (x, t)d , ( 5 )
where x = (x 1 , x 2 ) and ξ = (ξ 1 , ξ 2 ) represent the field and the source points, p * ij are the traction fundamental solutions, u j (x, t) are the CODs defined by

u i (x, t) = u i (x ∈ + c , t) -u i (x ∈ - c , t), ( 6 )
and an asterisk * denotes Riemann convolution which is defined by

g(x, t) * h(x, t) = t 0 g(x, t -τ )h(x, τ )dτ . ( 7 )
The traction fundamental solution p * ij is related to the displacement fundamental solution u * ij by

p * ij = C jkln ∂u * il ∂x n n k (x). ( 8
)
Substituting Eq. ( 5) into Hooke's law [START_REF] Albuquerque | The boundary element method applied to time dependent problems in anisotropic materials[END_REF] we obtain an integral representation formula for the traction components as

p j (ξ ) = - + c s * ij (ξ , x, t) * u j (x, t)d , ( 9 )
where

s * ij = C ikln ∂p * lj ∂ξ n n k (ξ ), ( 10
)
and n k (ξ ) represents the outward unit normal vector to the boundary at the collocation point. By taking the limiting process ξ → + c , the following time-domain traction BIEs are obtained

p j (ξ ) = -= + c s * ij (ξ , x, t) * u j (x, t)d , (11)
where = stands for the finite-part integral of Hadamard.

Frequency-domain and Laplace-domain TBIEs can be obtained directly by applying the Fourier-transform and the Laplace-transform to Eq. [START_REF] Lubich | Convolution quadrature and discretized operational calculus, Part I[END_REF].

Elastodynamic fundamental solutions

Unlike in the isotropic case, elastodynamic fundamental solutions for homogeneous, anisotropic and linear elastic solids cannot be given in closed-forms and they are much more complicated. For 2-D case, fundamental solutions can be represented by a line-integral over a unit circle, while for 3-D case they can be given as a surface-integral over a unit sphere. In the present analysis, two different elastodynamic fundamental solutions are adopted, namely, the frequency-domain and the Laplace-domain dynamic fundamental solutions.

Two-dimensional frequency-domain displacement fundamental solution has the following expression [START_REF] Wang | Elastodynamic fundamental solutions for anisotropic solids[END_REF]

u * ij (ξ , x, ω) = 1 8π 2 |η|=1 2 m=1 γ m ij ρ c 2 m (ω|η • (x -ξ )|/c m )dS η , (12)
where η = (η 1 , η 2 ) is the wave propagation vector, ω is the circular frequency, γ m ij = V im V jm (no sum over m) is the projection operator, with V im and ρc 2 m being the eigenvectors and the eigenvalues of the Christoffel ten-

sor ij ij (η 1 , η 2 ) -ρc 2
m δ ij V jm = 0, (no sum over m), [START_REF] Melenk | The partition of unity finite element method. Basic theory and applications[END_REF] in which δ ij denotes the Kronecker delta, c m represents the phase velocities of elastic waves, and the Christoffel tensor is defined by

ij (η 1 , η 2 ) = C ikjl η k η l . ( 14
)
The function (ζ ) is given by

(ζ ) = iπ e iζ -2 [cos(ζ ) ci(ζ ) + sin(ζ ) si(ζ )] , (15)
where ci and si denote the cosine and the sine integrals which are defined by

ci(ζ ) = -- ∞ ζ cos t t dt; si(ζ ) = - ∞ ζ sin t t dt, (16)
with ζ being real andindicating the Cauchy principal value integral. Two-dimensional Laplace-domain displacement fundamental solution has been recently published in [START_REF] Wang | 3-D and 2-D dynamic Green's functions and time-domain BIEs for piezoelectric solids[END_REF] for piezoelectric materials. From [START_REF] Wang | 3-D and 2-D dynamic Green's functions and time-domain BIEs for piezoelectric solids[END_REF], 2-D Laplace-domain displacement fundamental solution for homogeneous, anisotropic and linear elastic solids can be obtained as

u * ij (ξ , x, s) = 1 8π 2 |η|=1 2 m=1 γ m ij ρ c 2 m (s|η • (x -ξ )|/c m )dS η , ( 17
)
where s is the Laplace-transform parameter,

(z) = -e -z Ei(z) + e z Ei(-z) , (18)
with z being complex and Ei being the exponential integral defined by

Ei(z) = -- ∞ -1 e -zt t dt; Ei(-z) = - ∞ 1 e -zt t dt. ( 19
)
It can be easily shown that both functions and have a logarithmic singularity. To deal with this singularity it is advantageous to split the fundamental solutions into a singular static part and a regular dynamic part as [START_REF] Wang | Elastodynamic fundamental solutions for anisotropic solids[END_REF][START_REF] Wang | 3-D and 2-D dynamic Green's functions and time-domain BIEs for piezoelectric solids[END_REF]

u * ij (ξ , x, ω) = u R * ij (ξ , x, ω) + u S * ij (ξ , x), ( 20
)
u * ij (ξ , x, s) = u R * ij (ξ , x, s) + u S * ij (ξ , x), ( 21
)
where the superscripts R and S denote the regular dynamic part and the singular static part, respectively. The singular static part is independent of ω and s and has the following expression [START_REF] Sih | On cracks in rectilinearly anisotropic bodies[END_REF] while the regular dynamic part can be written as

u S * ij (ξ , x) = - 1 4π 2 |η|=1 2 m=1 γ m ij ρc 2 m log |η • (x -ξ )|dS η ,
u R * ij (ξ , x, ω) = 1 8π 2 |η|=1 2 m=1 γ m ij ρ c 2 m R (ω/c m , |η•(x-ξ )|)dS η , ( 23
) u R * ij (ξ , x, s) = 1 8π 2 |η|=1 2 m=1 γ m ij ρ c 2 m R (s/c m , |η•(x-ξ )|)dS η , (24)
with the regular continuous functions

R (x, y) = (x y) + 2 log y, (25) R (x, y) = (x y) + 2 log y. ( 26
)
The static displacement fundamental solution ( 22) can be reduced to the following explicit expression [START_REF] Eshelby | Anisotropic elasticity with applications to dislocation theory[END_REF][START_REF] García | Traction boundary elements for cracks in anisotropic solids[END_REF][START_REF] García-Sánchez | Two-dimensional time-harmonic BEM for cracked anisotropic solids[END_REF]

u S * ij (ξ , x) = - 1 π Re 2 m=1 P jm Q mi ln z m -z 0 m -ln(i + µ m ) , (27)
where

z m = x 1 + µ m x 2 , z 0 m = ξ 1 + µ m ξ 2 (28)
are the complex counterpart of the integration and the collocation points, the matrices P jm , Q mi and the complex constants µ m are given in the Appendix. By substituting Eqs. ( 20)-( 24) into Eqs. ( 8) and ( 10), the corresponding traction and higher-order traction fundamental solutions p * ij and s * ij can be obtained as

p * ij (ξ , x, ω) = p R * ij (ξ , x, ω) + p S * ij (ξ , x), ( 29
)
p * ij (ξ , x, s) = p R * ij (ξ , x, s) + p S * ij (ξ , x), ( 30
)
s * ij (ξ , x, ω) = s R * ij (ξ , x, ω) + s S * ij (ξ , x), ( 31
)
s * ij (ξ , x, s) = s R * ij (ξ , x, s) + s S * ij (ξ , x), (32)
where

p S * ij (ξ , x) = 1 π Re 2 m=1 L jm Q mi µ m n 1 (x) -n 2 (x) z m -z 0 m , ( 33
) p R * ij (ξ , x, ω) = 1 8π 2 |η|=1 2 m=1 S m ij ρ c 2 m • ω c m • (ω|η•(x-ξ )|/c m ) ×sign η • (x -ξ ) dS η , ( 34
) p R * ij (ξ , x, s) = 1 8π 2 |η|=1 2 m=1 S m ij ρ c 2 m • s c m • (s|η•(x-ξ )|/c m ) ×sign η • (x -ξ ) dS η , ( 35
)
s S * ij (ξ , x) = 1 π Re 2 m=1 T m ij µ m n 1 (x) -n 2 (x) z m -z 0 m 2 , ( 36
) s R * ij (ξ , x, ω) = 1 8π 2 |η|=1 2 m=1 U m ij ρ c 2 m ω c m 2 × (ω|η • (x -ξ )|/c m )dS η , ( 37
) s R * ij (ξ , x, s) = 1 8π 2 |η|=1 2 m=1 U m ij ρ c 2 m s c m 2 × (s|η • (x -ξ )|/c m )dS η . ( 38
)
In Eqs. ( 33)-( 38),

(ζ ) = d (ζ )/dζ , (z) = d (z)/dz,
and the auxiliary functions L jm , S m ij , T m ij and U m ij are given in the Appendix.

Numerical implementation of the BEM

To solve the hypersingular TBIEs [START_REF] Lubich | Convolution quadrature and discretized operational calculus, Part I[END_REF] as well as their counterparts in the frequency-domain and the Laplacetransformed domain, a collocation method for the spatial discretization by using quadratic elements is developed. Discontinuous quadratic elements are adopted over the crack-face. It should be mentioned here that the use of discontinuous elements is necessary in order to fulfill the C 1 -continuity requirement of the CODS in the hypersingular TBIEs.

Treatment of hypersingular integrals

As mentioned previously, the TBIEs [START_REF] Lubich | Convolution quadrature and discretized operational calculus, Part I[END_REF] and their counterparts in the transformed domain involve hypersingular integrals of the type 1/(z mz 0 m ) 2 when the integration point coincides with the collocation point. After the discretization, the hypersingular integrals to be computed have the following expression

I m == e µ m n 1 (x) -n 2 (x) (z m -z 0 m ) 2 φ d , ( 39
)
where e is the boundary element under consideration, φ denotes the quadratic shape-function and n(x) is the outward unit normal vector to the boundary. By means of a suitable variable change [START_REF] García | Traction boundary elements for cracks in anisotropic solids[END_REF][START_REF] García-Sánchez | Two-dimensional time-harmonic BEM for cracked anisotropic solids[END_REF], the hypersingular integral (39) can be regularized as described in the following.

Introducing a complex distance between the collocation and the integration points as a new variable

χ m = z m -z 0 m = (x 1 -ξ 1 ) + µ m (x 2 -ξ 2 ), ( 40
)
the Jacobian of this transformation is then given by

dχ m d = dχ m dx 1 dx 1 d + dχ m dx 2 dx 2 d = µ m n 1 (x) -n 2 (x), ( 41
)
where the following relations

dx 1 d = cos(θ ) = -n 2 (x), dx 2 d = sin(θ ) = n 1 (x) (42)
have been used (see Fig. 1). By using Eqs. ( 40) and (41), the hypersingular integral (39) can be rewritten as Considering the shape function φ as a function of the complex variable χ m , and using the first two terms of its Taylor-series expansion at χ m = 0, i.e.,

I m == e 1 χ 2 m φ dχ m . ( 43
)
φ(χ m ) = φ(χ m = 0) + dφ dχ m χ m =0 χ m + O χ 2 m = φ 0 + φ 0 χ m + O χ 2 m , (44)
the integral I m can be recast into

I m = e φ -(φ 0 + φ 0 χ m ) χ 2 m dχ m + φ 0 = e 1 χ 2 m dχ m +φ 0 - e 1 χ m dχ m . ( 45
)
The first integral in Eq. ( 45) is regular and can be computed numerically by using standard Gaussian quadrature formula. The second and the third integral are hypersingular and strongly singular, but they can be evaluated analytically as

= e 1 χ 2 m dχ m = - 1 χ m e , - e 1 χ m dχ m = log (χ m )| e . ( 46
)
4.2 Frequency-domain and Laplace-domain BEM

In the frequency-domain and the Laplace-domain BEM, the boundary value problem is first solved numerically for discrete values of the frequency and the Laplacetransform parameters. Subsequently, the corresponding time-domain solutions are obtained by using the fast Fourier inverse transform in the frequency-domain BEM and the inverse Laplace-transform algorithm of Stehfest [START_REF] Stehfest | Comun. Algorithm 368: numerical inversion of Laplace transform; an overview and recent developments[END_REF] in the Laplace-domain BEM.

According to the Stehfest's inversion algorithm [START_REF] Stehfest | Comun. Algorithm 368: numerical inversion of Laplace transform; an overview and recent developments[END_REF] a time-dependent function f (t) can be approximated by

f (t) = ln 2 t N n=1 ν n f (n ln 2/t), ( 47
)
where f (s) is the Laplace-transform of f (t), and

ν n = (-1) n+N/2 × min{n,N/2} k=(n+1)/2 (2k)! k N/2 (N/2 -k)!k!(k -1)!(n-k)!(2k-n)! . ( 48
)
Stehfest suggested to use a single precision arithmetic and N = 10 for the truncation limit in order to obtain accurate results. In this work, N =12 and 24 with a double precision arithmetic are used.

In the course of this study, we have also tested the Durbin's method for the inverse Laplace-transform, which has been preferred and suggested in some previous investigations on dynamic problems [START_REF] Gaul | A comparative study of three boundary element approaches to calculate the transient response of viscoeleastic solids with unbounded domains[END_REF][START_REF] Narayanan | Numerical operational methods for time-dependent linear problems[END_REF]. Our own experiences have confirmed that the Durbin's method is generally more accurate but also more complicated and computationally much more expensive than the Stehfest's algorithm. The Durbin's method uses a complex Laplace-transform parameter, while the Stehfest's algorithm applies a real Laplace-transform parameter and is easier to implement. For the dynamic crack problems considered in the present study, no significant changes in the numerical results have been noted by using these two different inversion methods. For this reason, the Stehfest's inversion algorithm is applied in our numerical examples, which will be presented in Sect. 6.

Time-domain BEM

In the time-domain BEM, the Riemann convolution integral

f (t) = g(t) * h(t) = t 0 g(t -τ )h(τ )dτ (49)
is approximated by the quadrature formula of Lubich [START_REF] Lubich | Convolution quadrature and discretized operational calculus, Part I[END_REF][START_REF] Lubich | Convolution quadrature and discretized operational calculus, Part II[END_REF], which is given by

f (n t) = n j=0 ω n-j ( t)h(j t), n = 0, 1, 2, . . . , N, ( 50
)
where the time-interval t is divided into N equal time-steps t, and the weights ω n ( t) are defined by

ω n ( t) = r -n N N-1 m=0 g δ(ζ m ) t e -2πinm/N . ( 51
)
In Eq. ( 51), -g(•) is the Laplace-transform of g(t),

-δ(ζ m ) = 2 j=1 (1 -ζ m ) j /j,
ζ m = re 2πim/N , and -r = 1/2N , with being the numerical error in computing the Laplace-transform g(•).

It should remarked here that a backward differential formula (BDF) of the second order is used for δ(ζ m ), and the error in computing ω n ( t) is of the order O( √ ) [START_REF] Lubich | Convolution quadrature and discretized operational calculus, Part I[END_REF][START_REF] Lubich | Convolution quadrature and discretized operational calculus, Part II[END_REF]. In this analysis, different -values between 10 -6 and 10 -12 have been tested to verify the effects of on the numerical accuracy. Our numerical tests have shown that in this -range there are no notable influences of on the numerical results. In contrast to the conventional time-domain BEM [START_REF] Tan | A time-domain BEM for transient wave scattering analysis by a crack in anisotropic solids[END_REF][START_REF] Tan | A time-domain collocation-Galerkin BEM for transient dynamic crack analysis in anisotropic solids[END_REF][START_REF] Ch | A 2-D hypersingular time-domain traction BEM for transient elastodynamic crack analysis[END_REF], the present method uses the Laplace-domain instead of the time-domain elastodynamic fundamental solutions. This is advantageous for cases where timedomain dynamic fundamental solutions are not available but their Laplace-transforms can be obtained. As representative examples we just mention the viscoelastic and the dynamic poroelastic problems, which have been investigated in details by Schanz [START_REF] Schanz | A boundary element formulation in time domain for viscoelastic solids[END_REF], Gaul and Schanz [START_REF] Gaul | A comparative study of three boundary element approaches to calculate the transient response of viscoeleastic solids with unbounded domains[END_REF] and Schanz [START_REF] Schanz | Wave propagation in viscoelastic and poroelastic continua[END_REF] using the quadrature formula of Lubich [START_REF] Lubich | Convolution quadrature and discretized operational calculus, Part I[END_REF][START_REF] Lubich | Convolution quadrature and discretized operational calculus, Part II[END_REF].

For the present crack problem in an infinite anisotropic domain, the Laplace-domain system matrix, denoted by A(s m ), can be computed by using the following equation

A(s m ) = E e=1 e s * ij ξ , x, ζ m t φ(x) d (52)
According to Eq. (51), the system matrix at the nth timestep A n can be obtained as

A n = r -n N N-1 m=0 A(s m ) e -2π inm/N , ( 53
)
where

s m = δ(ζ m )/ t.
Finally, a system of linear algebraic equations for the discrete CODs can be obtained as

n j=0 A n-j u j = p j . ( 54
)
By invoking the zero initial conditions (3), Eq. (54) leads to the following explicit time-stepping scheme [START_REF] Ch | Transient elastodynamic antiplane crack analysis of anisotropic solids[END_REF][START_REF] Ch | A 2-D time-domain BIEM for dynamic analysis of cracked orthotropic solids[END_REF][START_REF] Ch | Transient dynamic crack analysis in anisotropic solids[END_REF]]

u n = A 0 -1 ⎛ ⎝ p n - n-1 j=1 A n-j u j ⎞ ⎠ (55)
for computing the unknown CODs at the nth time-step.

In Eq. ( 55), A 0 -1 is the inverse of the system matrix A 0 at the time-step n = 0.

Computation of SIFs

In all three BEM as presented in previous sections, straight quarter-point elements are adopted at the crack-tips in order to capture the local √ r-behavior of the CODs near the crack-tips. This allows us to compute the dynamic SIFs very efficiently and accurately.

In the vicinity of the crack-tip, the displacement field has the following asymptotic expressions [START_REF] Sih | On cracks in rectilinearly anisotropic bodies[END_REF]

u 1 = 2r π K I Re 1 µ 1 -µ 2 µ 1 p 2 cos θ + µ 2 sin θ -µ 2 p 1 cos θ + µ 1 sin θ + K II Re 1 µ 1 -µ 2 p 2 cos θ + µ 2 sin θ -p 1 cos θ + µ 1 sin θ , ( 56
)
u 2 = 2r π K I Re 1 µ 1 -µ 2 µ 1 q 2 cos θ + µ 2 sin θ -µ 2 q 1 cos θ + µ 1 sin θ + K II Re 1 µ 1 -µ 2 q 2 cos θ + µ 2 sin θ -q 1 cos θ + µ 1 sin θ , ( 57
)
where r and θ are polar coordinates with the origin at the crack-tip, K I and K II are the mode-I and the mode-II SIFs, and

p m = b 11 µ 2 m + b 12 -b 16 µ m , q m = (b 12 µ 2 m + b 22 -b 26 µ m )/µ m . ( 58
)
In Eq. (58), b ij (i, j = 1, 2, 6) are the materials compliance coefficients. By substituting Eqs. ( 56) and (57) into Eq. ( 6), we obtain the following relationship between the dynamic SFIs and the CODs

K I K II = π 8d • 1 H 11 H 12 H 21 H 22 u 1 u 2 , ( 59
)
where

H ij = ⎡ ⎣ Im q 1 -q 2 µ 1 -µ 2 Im p 2 -p 1 µ 1 -µ 2 Im µ 1 q 2 -µ 2 q 1 µ 1 -µ 2 Im µ 2 p 1 -µ 1 p 2 µ 1 -µ 2 ⎤ ⎦ , (60) = H 11 H 22 -H 12 H 21 . (61)
In Eq. ( 59), u 1 and u 2 are the CODs at the collocation point closest to the crack-tip either in the frequency-domain, or in the Laplace-transformed domain or in the time-domain depending on the applied BEM, and d is the distance between this collocation point and the crack-tip. Using Eq. ( 59), dynamic SIFs can be obtained directly once the CODs have been computed numerically by the BEM. In the frequency-domain and the Laplace-domain BEM, an inverse transform is subsequently applied to obtain time-dependent SIFs. In contrast, time-domain BEM yields directly time-dependent SIFs.

Numerical results

To compare the three different BEM, both straight and curved cracks are analyzed. The straight crack has a finite length 2a, and is considered to be either in an isotropic or a fully anisotropic unbounded plane subjected to a tensile or a shear impact crack-face loading (see Fig. 2). As an example for curved cracks, we consider a circular arc-shaped crack with a radius r and opening angle 2α embedded in a fully anisotropic plane, which is subjected to a radial impact crack-face loading as shown in Fig. 3. In all cases, the accuracy and the efficiency of the three BEM are compared. Furthermore the stability of the time-domain BEM is studied by using different time-steps.

Numerical calculations have been carried out by using ten elements as depicted in Fig. 4. For the straight crack, the size-ratio between two consecutive elements is constant, which is equal to 2 between the central and the crack-tip element. For the curved crack, the straight crack-tip elements have a length of rα/30, while the first three elements from the crack-center have an angle 0.50α, 0.30α and 0.15α, respectively. The size of the fourth element from the crack-center is determined by the remaining size of the crack. For convenience of presentation and discussion, the following abbreviations are introduced:

-FD = frequency-domain BEM, -LD = Laplace-domain BEM, -TD = Time-domain BEM, -FFT = Fast Fourier Transform, -LP = Laplace-transform parameter. Fig. 4 BEM meshes used for straight and arc-shaped (opening angle = 2α) cracks

Isotropic solids

In order to compare our numerical results with available analytical solutions we first consider a straight finite crack of length 2a in an unbounded, homogeneous, isotropic and linearly elastic solid. The crack is subjected to a tensile impact crack-face loading σ 22 (t) = σ 0 H(t) t c L / a t c L / a Fig. 5 Normalized dynamic SIFs for a crack in an isotropic solid subjected to a tensile and a shear impact crack-face loading or a shear impact crack-face loading σ 12 (t) = τ 0 H(t), where σ 0 and τ 0 are the loading amplitudes and H(t) is the Heaviside step function. Poisson's ratio is taken as 1/4.

Numerical results are presented and compared with the analytical solutions [START_REF] Thau | Transient stress intensity factors for a finite crack in an elastic solid caused by a dilatational wave[END_REF] in Figs. 5 and6. Here, the dynamic SIFs are normalized by their corresponding static values K st I = σ 0 √ πa and K st II = τ 0 √ πa, and the dimensionless time tc L /a is used for convenience, where c L is the velocity of the longitudinal wave. In the timedomain BEM, a time-step c T t = a/20 is selected. As we can see in Figs. 5 and6, the numerical results obtained by all three BEM show very good agreement with the analytical solutions, though some small discrepancies occur near the peak values of the normalized dynamic SIFs.

To study the stability of the time-domain BEM, additional calculations for three different time-steps have been performed. The corresponding numerical results are presented in Fig. 6, which imply that the present time-domain BEM is quite insensitive to the selected time-steps. This is an important advantage over the classical time-domain BEM [START_REF] Tan | A time-domain BEM for transient wave scattering analysis by a crack in anisotropic solids[END_REF][START_REF] Tan | A time-domain collocation-Galerkin BEM for transient dynamic crack analysis in anisotropic solids[END_REF][START_REF] Ch | A 2-D hypersingular time-domain traction BEM for transient elastodynamic crack analysis[END_REF], which suffers from the stability problem.

General anisotropic solids

Straight crack

Next, we consider a straight finite crack of length 2a in an unbounded, homogeneous and linear elastic solid with a general anisotropy. The following elastic constants have The selected material constants correspond to a Graphite-epoxy composite [START_REF] Ch | Transient dynamic crack analysis in anisotropic solids[END_REF]. The same BEM mesh as in the isotropic case has been applied. For a tensile impact crack-face loading, the normalized mode-I and mode-II dynamic SIFs versus the dimensionless time t c T /a are presented in Fig. 7, while the corresponding numerical results for a shear impact crack-face loading are given in Fig. 8. Here, the shear wave velocity is defined by c T = √ C 66 /ρ, and a timestep c T t = a/20 is used in the time-domain BEM. The results are compared with that obtained by Zhang [START_REF] Ch | Transient dynamic crack analysis in anisotropic solids[END_REF], who used a Galerkin-method for the spatial discretization.

The good agreement between our numerical results and that of Zhang [START_REF] Ch | Transient dynamic crack analysis in anisotropic solids[END_REF] as shown in Figs. 7 and8 verifies the accuracy of the three implemented BEM. As in the isotropic case, our time-domain BEM is quite stable.

Circular arc-shaped crack

Now we consider a circular arc-shaped crack embedded in an unbounded anisotropic domain with the same material constants as used in the preceding case. The crack is defined by the radius r of the circle and the opening angle 2α as shown in Fig. 3. A tensile impact crack-face loading in the radial direction in the form of σ r = σ 0 H(t) is considered.

The numerical results are presented as normalized mode-I and mode-II dynamic SIFs versus dimensionless time, which is defined as in the preceding case for a straight crack. The dynamic SIFs are normalized by σ 0 √ π a, where σ 0 is the loading amplitude and 2a is the arc-chord length of the crack. A comparison of the numerical results obtained by the three BEM is shown in Figs. 9, 10, and 11 for the two crack-tips. Also here, a time-step c T t = a/20 is used in the time-domain BEM. As expected, the differences between the dynamic SIFs at both crack-tips increase as the curvature of the arc-shaped crack increases. The numerical results show that, except for the mode-II SIF at the crack-tip B, the time at which the peaks of the SIFs are achieved does not vary substantially with the curvature in the cases of α = 15 • and α = 30 • . In the case of α = 60 • , the peak K I -factor at the crack-tip A is shifted to a smaller time instant. The results obtained by the time-domain BEM and the frequency-domain BEM agree very well, while the results obtained by the Laplace-domain BEM using 12 Laplace-transform parameters show some differences. The main difference lies in the tendency that the Laplace-domain BEM In order to investigate the influence of the used number of the Laplace-transform parameters on the accuracy of the Laplace-domain BEM, numerical calculations have also been carried out for 24 Laplace-transform parameters. The results show that the agreement with other two BEM are improved now in the first part of the plot, i.e., for a dimensionless time below about 5.5. However, in the large time-range, oscillation of the numerical results appears which indicates an instability of the method.

The stability of the time-domain BEM is tested here again for three opening angles of the arc-shaped crack. The numerical results for the normalized dynamic SIFs are presented in Figs. 12, 13, and 14 versus the dimensionless time tc T /a. For three quite different time-steps selected, no essential differences in the numerical results are observed, which confirms again that the present time-domain BEM is quite insensitive to the used timesteps. Moreover, the curvature of the crack-faces does not affect the stability behavior of the present timedomain BEM, at least for the investigated cases in this analysis. Mesh-sensitivity of the dynamic BEM is a critical issue for dynamic problems. For a straight crack in an anisotropic solid subjected to a tensile impact loading, a mesh-sensitivity study of the three presented BEM is carried out for different numbers and distributions of the boundary elements. The used material constants are given in Eq. (62). The corresponding numerical results for the normalized dynamic SIFs are presented in Figs. 15, 16, 17, 18, 19, 20, and 21. From this study the following conclusions can be drawn:

-The three BEM presented in this paper are quite insensitive to the used element-number.

-Ten quadratic elements are sufficient to obtain accurate numerical results. -A uniform mesh is also adequate for the dynamic crack problem under consideration, although a nonuniform mesh is applied in most part of the analysis.

In addition, the following comments to the mesh-sensitivity of the dynamic BEM should be made:

-In the classical time-domain BEM, the spatial discretization and the time discretization are not independent of each other. This implies that the mesh-sensitivity and the time-step sensitivity cannot A too small element-size (or too small time-step) may lead to an instability of the numerical scheme, while a too large element-size (or too large time-step) may cause a physically unrealistic large numerical damping of the results. A brief review on the subject can be found in reference [START_REF] Ch | 3-D transient dynamic crack analysis by a novel time-domain BEM[END_REF]. To the best knowledge of the authors, yet there are no rigorous mathematical stability proofs for the classical time-domain BEM based on collocation methods both for spatial and time discretizations. However, many previous including our own numerical experiences show that for isotropic elastic solids stable and accurate numerical results can be obtained by using c L t ≤ l ≤ 2c L t [START_REF] Dominguez | Boundary elements in dynamics[END_REF], where c L is the longitudinal wave velocity, l is the element-size, and t is the used time-step. The situation for anisotropic elastic solids becomes even more tangled due to the directional dependence of the wave velocities. The present time-domain BEM by using Lubich's quadrature is in any case less sensitive to the used mesh-size in comparison to the classical time-domain BEM, since the present timedomain BEM is less sensitive to the used time-steps, which implies implicitly a less sensitivity of the method to the spatial mesh-size. -In the frequency-domain BEM for time-harmonic analysis, the mesh-sensitivity is governed by the ratio of the nodal spacing and the wave-length. To reliably capture the sinusoidal waveform by using polynomial shape functions in the frequency-domain BEM, sufficient nodes within the wave-length are required.

A commonly quoted and accepted rule of thumb recommends 8-10 nodes per wave-length (i.e., four or five quadratic elements) are needed. For large scale problems and high frequencies (short wave-lengths) this may cause some serious difficulties due to the substantially large memory, storage and computing time, which may exceed the available computer resources. To overcome this difficulty, several advanced methods have been proposed in recent years. Here, we just mention the fast multipole method (FMM) [START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF][START_REF] Rokhlin | Rapid solution of integral equations of classical potential theory[END_REF][START_REF] Rokhlin | Rapid solution of integral equations of scattering theory in two dimensions[END_REF], the wave boundary element method (WBEM) or the wave basis functions method [START_REF] Bettess | Short wave scattering: problems and techniques[END_REF][START_REF] Perrey-Debain | P-wave and S-wave decomposition in boundary integral equation for plane elastodynamic problems[END_REF][START_REF] Perrey-Debain | Wave boundary elements: a theoretical overview presenting applications in scattering of short waves[END_REF] based on the partition of unity method (PUM) [START_REF] Melenk | The partition of unity finite element method. Basic theory and applications[END_REF]. A review of some advanced computational methods for wave simulation in high frequency range has been presented by Bettess [START_REF] Bettess | Short wave scattering: problems and techniques[END_REF].

For the present transient dynamic crack analysis, numerical calculations for different frequencies are are computed via a regularization technique based on a suitable variable change. This allows us to recast a hypersingular integral into a regular integral plus a strongly singular integral and another new hypersingular integral. Regular integrals are computed numerically by using standard Gaussian quadrature formula, while the strongly singular and the new hypersingular integrals are evaluated analytically.

In the frequency-domain and the Laplace-domain BEM, the boundary value problem is first solved in the transformed domain for discrete frequency and Laplacetransform parameters. Then, time-domain results are obtained by inverse Fourier-transform and inverse Laplace-transform. In the time-domain BEM, the quadrature formula of Lubich is applied for approximating the Riemman convolution integrals, which requires only Laplace-domain fundamental solutions instead of time-domain fundamental solutions. At crack-tips, quarter-point elements are adopted in all three BEM to describe the local behavior of the CODs at the cracktips properly. Dynamic SIFs are obtained directly from the numerically computed CODs.

The implemented BEM presented in this paper are general and can be used for straight and curved cracks. Numerical examples for computing dynamic SIFs show that the implemented BEM are very accurate and robust, and they are valid even for isotropic material properties as a special case of the anisotropic materials. It is well known that the isotropic case is a mathematically degenerated case of the general anisotropic BEM formulation, which may cause some numerical difficulties.

Regarding the computational efficiency or cost, the time-domain BEM is the fastest method while the Laplace-domain BEM with the Stehfest's inversion algorithm is the most expensive one, see Fig. 21. In order to obtain the numerical results for N time-steps, the system matrices have to be computed for -N/2 + 1 Laplace-transform parameters in the timedomain BEM using the Lubich's quadrature, -approximately 2N frequencies in the frequencydomain BEM using the fast Fourier inverse transform, and -n × N Laplace-transform parameters in the Laplacedomain BEM using the Stehfest's inversion algorithm, with n being the truncation limit (n =12 or 24 in this analysis).

On the other hand, the Laplace-domain BEM with the Stehfest's inversion algorithm have two important advantages: it is easy to implement and it uses only a real Laplace-transform parameter.

Though the present crack analysis is shown in this paper for unbounded anisotropic solids, the extension of the three BEM to transient dynamic crack analysis in bounded anisotropic solids is straight-forward and the corresponding results will be reported in future.

Fig. 1

 1 Fig. 1 Differential boundary element geometry around the integration point with the outward unit normal vector

Fig. 2 AFig. 3 A

 23 Fig. 2 A straight crack in an unbounded domain subjected to a tensile impact loading (left) and a shear impact loading (right)

Fig. 6 ⎣ C 11 C 12 C 16 C 22 ρ = 1 ,

 6111216221 Fig.[START_REF] Dominguez | Boundary elements in dynamics[END_REF] Influences of time-steps on the normalized dynamic SIF for a crack in an isotropic solid subjected to a tensile impact crackface loading

Fig. 7 Fig. 8

 78 Fig. 7 Normalized mode-I and mode-II dynamic SIFs for a crack in an anisotropic solid subjected to a tensile impact loading

Fig. 9 Fig. 10

 910 Fig. 9 Normalized mode-I and mode-II dynamic SIFs for an arc-shaped crack with α = 15 •

Fig. 11 Fig. 12

 1112 Fig. 11 Normalized mode-I and mode-II dynamic SIFs for an arc-shaped crack with α = 60 •

Fig. 13

 13 Fig. 13 Normalized mode-I and mode-II dynamic SIFs for an arc-shaped crack with α = 30 • by using three different time-steps

Fig. 15 Fig. 16

 1516 Fig. 15 Mesh-sensitivity of the frequency-domain BEM (uniform mesh)

Fig. 17

 17 Fig. 17 Mesh-sensitivity of the Laplace-domain BEM (uniform mesh)

Fig. 19 Fig. 20

 1920 Fig. 19 Mesh-sensitivity of the time-domain BEM (uniform mesh)

7 ConclusionsFig. 21

 721 Fig.[START_REF] Schanz | Wave propagation in viscoelastic and poroelastic continua[END_REF] Comparison of the computational costs of the three BEM (t L =19 minutes 32 seconds)

  (C irk1 + µ m C irk2 ) n r (ξ )L jm Q mk , (no sum over m),

	The auxiliary functions S m ij , T m ij and U m ij in Eqs. (34)-(38)
	are given by
	S m ij = C jklr γ m ir n k (x)η l ,	(66)
	T m ij = (67)
	U m ij = -C irkl S m kj n r (ξ )η l .	(68)

Acknowledgements This work is supported by the German Research Foundation (DFG) under the project number ZH 15/5-1 and ZH 15/5-2, which is gratefully acknowledged.

Appendix: Auxiliary functions

In Eq. ( 27), µ m are the complex roots of the following characteristic equation

The matrices P and Q are determined by

where

(no sum over m).

(65)


Vilius Kurauskas

Emmanuelle Weber

Audrey Hessel

Isabel Ayala

Dominique Marion
email: dominique.marion@ibs.fr

Paul Schanda
email: paul.schanda@ibs.fr

Cross-Correlated Relaxation of Dipolar Coupling and Chemical-Shift Anisotropy in Magic-Angle Spinning R 1ρ NMR Measurements: Application to Protein Backbone Dynamics Measurements

Keywords: protein, rotating frame relaxation, numerical spin simulations, spin-diffusion, J-doublet

come

Introduction

The structure of a protein is determined by a large number of individually weak interactions, which constantly re-shuffle at ambient temperatures. This dynamic exchange between different conformers is often important for functions such as binding or enzymatic turnover. NMR spectroscopy is very well suited to provide direct atomic-level access to the amplitudes and time scales of such motions. Whereas solution-state NMR is a well-established tool for such studies, magic-angle spinning (MAS) solid-state NMR (ssNMR) is rapidly establishing in studies of biomolecular dynamics, complementing its solution-state counterpart in cases where limited solubility or large size hamper the application of solution-state methods. In addition to opening new fields of applications, ssNMR also has theoretical advantages for determining dynamics, compared to solution-state NMR. In particular, the absence of overall tumbling motion enables the study of motions without "blind windows" of time scales, i.e. over all time scales from picoseconds to seconds, in contrast to its solution-state counterpart, which is severely challenged on the nanosecond to several-microseconds time scale. [START_REF] Schanda | Studying Dynamics by Magic-Angle Spinning Solid-State NMR Spectroscopy: Principles and Applications to Biomolecules[END_REF] In MAS ssNMR, motions can be studied through (i) their averaging effect on anisotropic interactions (in particular, the averaging of dipolar couplings and chemical-shift anisotropies), which report on the amplitude of motion seen by these interactions, and (ii) through the investigation of nuclear spin relaxation rate constants, which depend on amplitudes and time scales of motions. Many biological processes, such as binding, gating or enzymatic turnover, occur on time scales of hundreds of nanoseconds to milliseconds (ns-ms). Studying these time scales is, therefore, of particular relevance when attempting to decipher the molecular basis of these processes. In MAS ssNMR, motion on these time scales can be accessed by measurements of transverse relaxation rate constants, such as the decay of spin coherence, R 2 . However, quantitative measurements of transverse relaxation rates in ssNMR are generally hampered by the fact that spin coherences in solid samples decay through two mechanisms, (i) due to random field fluctuations (i.e., dynamics), which is what we are interested in herein, and (ii) due to evolution under the incompletely averaged anisotropic interactions, in particular due to [START_REF] Schanda | Studying Dynamics by Magic-Angle Spinning Solid-State NMR Spectroscopy: Principles and Applications to Biomolecules[END_REF] H-1 H dipolar couplings. The latter, termed dipolar dephasing, strongly depends on the presence of (multiple) dipolar couplings, and on the employed MAS frequency. Accessing the dynamics-induced component requires that coherent contributions to signal decay (in particular dipolar dephasing) are suppressed to negligible levels. Several reported studies have shown that heteronuclear ( 13 C, [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N) R 1ρ rate constants, i.e. the decay of C x (N x ) in the presence of a spin-lock RF field, can be measured accurately without significant interfering effects from dipolar dephasing if (i) high MAS frequencies and/or (ii) extensive deuteration are applied. [START_REF] Krushelnitsky | Microsecond Time Scale Mobility in a Solid Protein as Studied by the 15N R 1ρ Site-Specific NMR Relaxation Rates[END_REF][START_REF] Lewandowski | Site-Specific Measurement of Slow Motions in Proteins[END_REF][START_REF] Ma | Probing Transient Conformational States of Proteins by Solid-State R1rho Relaxation-Dispersion NMR Spectroscopy[END_REF][START_REF] Lamley | Unraveling the Complexity of Protein Backbone Dynamics with Combined 13 C and 15 N Solid-state NMR Relaxation Measurements[END_REF] R 1ρ experiments under such conditions have been reported so far for [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N and [START_REF] Korzhnev | An NMR Experiment for the Accurate Measurement of Heteronuclear Spin-Lock Relaxation Rates[END_REF] C' sites in small model proteins (SH3, GB1, ubiquitin). [START_REF] Krushelnitsky | Microsecond Time Scale Mobility in a Solid Protein as Studied by the 15N R 1ρ Site-Specific NMR Relaxation Rates[END_REF][START_REF] Lewandowski | Site-Specific Measurement of Slow Motions in Proteins[END_REF][START_REF] Ma | Probing Transient Conformational States of Proteins by Solid-State R1rho Relaxation-Dispersion NMR Spectroscopy[END_REF][START_REF] Lamley | Unraveling the Complexity of Protein Backbone Dynamics with Combined 13 C and 15 N Solid-state NMR Relaxation Measurements[END_REF][START_REF] Lamley | Intermolecular Interactions and Protein Dynamics by Solid-State NMR Spectroscopy[END_REF][START_REF] Ma | Observing the Overall Rocking Motion of a Protein in a Crystal[END_REF] R 1ρ experiments in solids are particularly informative also when the RF field strength approaches the rotary resonance conditions, as has been shown for small molecules [START_REF] Quinn | Quantifying Conformational Dynamics Using Solid-State R₁ρ Experiments[END_REF][START_REF] Quinn | Monitoring Conformational Dynamics with Solid-State R₁ρ Experiments[END_REF] and proteins. [START_REF] Ma | Probing Transient Conformational States of Proteins by Solid-State R1rho Relaxation-Dispersion NMR Spectroscopy[END_REF][START_REF] Lamley | Intermolecular Interactions and Protein Dynamics by Solid-State NMR Spectroscopy[END_REF] As an alternative method to probe ns-ms motions, it has been proposed to measure cross-correlated relaxation (CCR) between 1 H-15 N dipolar/ 15 N chemical-shift anisotropy (CSA) interactions in R 2 'type experiments (i.e. in the absence of a spin-lock RF field). These latter experiments measure the difference of the apparent transverse decay rate constants of the two [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N scalar-coupling doublet components, i.e. the decay of N x -[H α ] and N x -[H β ] components in a spin-echo experiment. [START_REF] Chevelkov | Quantitative Measurement of Differential 15-N-H Α/β T2 Relaxation Rates in a Perdeuterated Protein by MAS Solid-State NMR Spectroscopy[END_REF] Even though each individual rate constant itself contains contributions from dipolar dephasing (and can, therefore, not be interpreted quantitatively), the respective difference can be quantitatively interpreted in terms of dynamics on ns-ms time scales. [START_REF] Skrynnikov | Asymmetric Doublets in MAS NMR: Coherent and Incoherent Mechanisms[END_REF] This experiment suffers, however, from the fact that the R 2 ' decay rates in MAS ssNMR are large, because of the dipolar-dephasing contribution. As a consequence, detection of a small differential relaxation rate from two individually large relaxation rate constants corresponding to the two doublet components (R 2 '+Γ CSA/D and R 2 '-Γ CSA/D ) is challenging. Furthermore, the rapid signal loss throughout the experiment, based on scalar-coupling transfers, [START_REF] Chevelkov | Quantitative Measurement of Differential 15-N-H Α/β T2 Relaxation Rates in a Perdeuterated Protein by MAS Solid-State NMR Spectroscopy[END_REF] results in relatively low detection sensitivity.

Here we investigate the effects of dipolar/CSA CCR in heteronuclear ( 13 C, 15 N) R 1ρ experiments in MAS ssNMR. In solution-state NMR, such cross-correlated relaxation effects have been reported [START_REF] Massi | NMR R 1ρ Rotating-Frame Relaxation with Weak Radio Frequency Fields[END_REF][START_REF] Korzhnev | An NMR Experiment for the Accurate Measurement of Heteronuclear Spin-Lock Relaxation Rates[END_REF] , but have so far been neglected in MAS ssNMR. We use explicit numerical spin dynamics simulations, as well as experimental measurements, and show that CCR rate constants under spin-lock RF fields are measurable, and that they report on molecular motion. The CCR rate constant has a similar dependence on dynamics as the in-phase R 1ρ rate constant, as we show with numerical simulations. We assume here that the CCR rate constant in R 1ρ measurements is essentially free from coherent contributions, for several reasons. First, it has been convincingly shown that CCR rate constants, derived from the above-mentioned experiments with free evolution (R 2 ') are primarily reflecting dynamics, even though each of the two doublet components itself is not free from coherent contributions. [START_REF] Chevelkov | TROSY Effects in MAS Solid-State NMR[END_REF] In R 1ρ -based experiments the situation can be assumed to be even more favorable, because the decay of each of the two components is much less prone to coherent decay. [START_REF] Lewandowski | Site-Specific Measurement of Slow Motions in Proteins[END_REF] Furthermore, we show that experimentally determined R 1ρ -derived CCR rate constants (referred to here as Γ 1 ρ CSA / D ) reveal the expected patterns of motions in a well-studied protein, and correlate with R 1ρ rate constants, which have been shown to contain only minimal coherent contributions. However, both CCR and R 1ρ contain contributions from stochastic proton flips: we provide evidence that proton spin-diffusion reduces the CCR rate constants due to averaging of the doublet components. Furthermore, proton spin diffusion also leads to an increase of the auto-relaxation (R 1ρ ) rate constants, in a manner that depends on the RF field strength. Our experimental data and simulations thus suggest that knowledge of the proton spin-diffusion rate is required for accurately determining dynamics, either through in-phase R 1ρ measurements or CCR measurements under spin-lock. Finally, we also discuss pulse schemes to suppress the effects of CCR, and thus eliminate the bi-exponentiality of decay in in-phase R 1ρ measurements.

Materials and Methods

Samples and NMR spectroscopy

In this study we investigate the cross-correlated relaxation between a heteronuclear CSA tensor ( 15 N) and the 1 H- [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N dipolar coupling in amide moieties; our conclusion apply also to other 1 H-X spin pairs, such as 1 H- [START_REF] Korzhnev | An NMR Experiment for the Accurate Measurement of Heteronuclear Spin-Lock Relaxation Rates[END_REF] C in backbone or side chain sites. We use u-[ 2 H, [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] H]-labelled ubiquitin, crystallized with 2-methyl-2,4-pentane-diol (MPD) [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] [START_REF] Korzhnev | An NMR Experiment for the Accurate Measurement of Heteronuclear Spin-Lock Relaxation Rates[END_REF] C, [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N frequencies, spinning at 38 kHz MAS frequency. The effective sample temperature, calibrated with internal standard using the bulk water frequency, was set to 300 K . Pulse sequences used in this study are shown in Figure S1 in the Supporting Information, and use 1 H-15 N out-andback cross-polarization, and 1 H detection. In particular, the sequence in Figure S1B was used for the CCR measurement, and (in-phase) R 1ρ rate constants with and without 1 H π pulse decoupling are shown in Figure S1A and C, respectively. Finally, the proton spin-diffusion experiment was performed with the sequence shown in Figure S1D. [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N spin-lock RF field frequencies were calibrated with a nutation experiment: after initial 1 H-15 N cross-polarization, a 5 ms spin-lock period was applied (to take into account possible detuning during a spin-lock, as this is relevant for a R 1ρ experiment), followed by a 90 degree phase-shifted nutation pulse. If this pulse is a π/2 pulse, then the signal is zero. We determined the amplitude of this pulse by parametric incrementation. We find that two alternative calibration methods result in very similar calibrations: (i) a simple nutation pulse, applied at the beginning of the solvent-saturation period (see Figure S1), i.e. disregarding the spin-lock period, and (ii) the observation of the rotary resonance condition. For the latter, the signal intensity is monitored as a function of the 15 N spin-lock RF field strength in an experiment that applies a 5 ms spin-lock just after the initial 1 H-15 N CP transfer (Figure S1A); a minimum is observed when the RF field matches the n=1 rotary-resonance condition ν RF = ν MAS . We find that the three methods provide essentially identical calibration results. Processing and analysis was done with nmrPipe [START_REF] Delaglio | NMRPIPE -a Multidimensional Spectral Processing System Based on Unix Pipes[END_REF] , and NMRView (OneMoon Scientific) as well as in-house written scripts in python language. Error estimates of rate constants were based on spectral noise estimates (3 times standard deviation of the signal in empty regions of the spectrum).

Numerical simulations

Simulations of the evolution of spin coherences were performed using in-house written programs using the GAMMA simulation library. [START_REF] Smith | Computer Simulations in Magnetic Resonance. An Object-Oriented Programming Approach[END_REF] The simulation approach includes two spins ( 1 H , 15 N), existing in two exchanging sub-systems, spanning a composite Liouville space. [START_REF] Schanda | Studying Dynamics by Magic-Angle Spinning Solid-State NMR Spectroscopy: Principles and Applications to Biomolecules[END_REF][START_REF] Skrynnikov | Asymmetric Doublets in MAS NMR: Coherent and Incoherent Mechanisms[END_REF][START_REF] Abergel | On the Use of the Stochastic Liouville Equation in Nuclear Magnetic Resonance: Application to R 1ρ Relaxation in the Presence of Exchange[END_REF][START_REF] Kubo | Stochastic Liouville Equations[END_REF] The simulation approach, as well as numerous results of similar kinds of simulations were reported recently. [START_REF] Schanda | Studying Dynamics by Magic-Angle Spinning Solid-State NMR Spectroscopy: Principles and Applications to Biomolecules[END_REF] Briefly, the evolution of the appropriate coherences (N x , 2N x H α , 2N x H β , H z ) is followed over the course of a [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N spin-lock. The two exchanging sub-states differ in the orientations of the 1 H-15 N dipolar coupling tensor and 15 N CSA tensor. The interaction tensor parameters were assumed as follows: 1 H- [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N dipolar coupling 22954 Hz (corresponding to 1.02 Å bond length), axially symmetric [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N tensor with Δσ=-170 ppm, i.e. σ zz =-113 ppm (6780 Hz at 14.1 T B 0 field strength). The 1 H CSA parameters have negligible effect in this context. The time-domain data were fitted with python scripts in order to obtain relaxation rate constants, using mono-or bi-exponential functions as reported in the text.

Results and Discussion

We first aim at providing a quantitative understanding of the effects of CSA/D CCR in R 1ρ experiments by numerical simulations. Although a very well-established theoretical framework for nuclear spin relaxation exists (Redfield theory [START_REF] Redfield | On the Theory of Relaxation Processes[END_REF] ), this theory is not strictly valid in cases of slow motion, which are the focus of this paper. Recent systematic analyses have established that Redfield results may disagree with numerical simulations, [START_REF] Schanda | Studying Dynamics by Magic-Angle Spinning Solid-State NMR Spectroscopy: Principles and Applications to Biomolecules[END_REF][START_REF] Ma | Observing the Overall Rocking Motion of a Protein in a Crystal[END_REF] even in cases where the actual relaxation is more than one order of magnitude slower than the dynamic event. In order to circumvent these problems we resort to numerical simulations for investigating the properties of R 1ρ relaxation in rotating solids. The simulation approach is based on an explicit jump model of a 1 H- [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N moiety exchanging between two conformations. The employed composite Liouville-space approach has been used similarly before [START_REF] Skrynnikov | Asymmetric Doublets in MAS NMR: Coherent and Incoherent Mechanisms[END_REF][START_REF] Abergel | On the Use of the Stochastic Liouville Equation in Nuclear Magnetic Resonance: Application to R 1ρ Relaxation in the Presence of Exchange[END_REF][START_REF] Kubo | Stochastic Liouville Equations[END_REF] .

Figure 1 shows an example of a simulated decay of N x coherence for a two-site exchange, using powder averaging and MAS. The observed decay (red curve in Figure 1A) is clearly nonexponential, and the exponential fit (black) unsatisfactorily describes the decay. The nonexponentiality may be ascribed to two effects. The first is that the NMR interactions that lead to relaxation (dipolar coupling, CSA) are orientation-dependent, and therefore the relaxation rate constants are orientation-dependent, too. When summing all contributions from different crystallite orientations across the sample, as required when considering a sample with randomly oriented crystallites, the resulting relaxation is always multi-exponential, as has been described and reviewed earlier. [START_REF] Schanda | Studying Dynamics by Magic-Angle Spinning Solid-State NMR Spectroscopy: Principles and Applications to Biomolecules[END_REF][START_REF] Torchia | Spin-Lattice Relaxation in Solids[END_REF][START_REF] Giraud | Quantitative Analysis of Backbone Dynamics in a Crystalline Protein from Nitrogen-15 Spin-Lattice Relaxation[END_REF] If this is the only reason for non-monoexponential behavior, then individual crystallites (i.e. with a single orientation within the rotor-fixed frame) should show mono-exponential decays. However, we find that also when using single crystallites rather than a powder distribution the decay is bi-exponential (data not shown), pointing to an additional mechanism.

This second possible origin of non-exponential decay is the presence of cross-correlated relaxation due to [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N-CSA/ 1 H-15 N dipolar-coupling interference. When observing the two components of the initial N x coherence, namely the coherences N x -[H α ] and N x -[H β ], the decays are nearly exponential when integrated over a powder (Figure 1B). If one considers individual crystallites then the decay of the two coherences are indeed perfectly exponential (Figure S2C,D). One can readily show that the bi-exponential decay of N x is due to interference of CSA and dipolar mechanisms: removing either of the two mechanisms (i.e. setting the jump angle for one of the two interactions to zero, or, equivalently, setting the interaction strength to zero) leads to identical relaxation of the two components, and mono-exponential decay of N x (data not shown). Furthermore, if a 1 H π pulse is inserted, thus inverting the sign of the dipolar interaction, the near-monoexponential decay is restored, as shown in Figure 1A (grey line). This is expected, as the effect of the cross-correlated relaxation is eliminated.

Numerical simulations also allow determining the difference of the rate constants of the two components, i.e. the CSA/D CCR rate constant in the R 1ρ experiment, which we term here Γ 1 ρ CSA / D .

Figure 2 establishes that the differential relaxation rate constant depends on the amplitude and time scale of the motion, as expected for a cross-correlated relaxation, and on the difference of RF field strength and MAS frequency, again as expected. [START_REF] Schanda | Studying Dynamics by Magic-Angle Spinning Solid-State NMR Spectroscopy: Principles and Applications to Biomolecules[END_REF][START_REF] Kurbanov | The Nuclear Magnetic Resonance Relaxation Data Analysis in Solids: General R1∕R1ρ Equations and the Model-Free Approach[END_REF] It exhibits a maximum when the motion is on the time scale of microseconds to tens of microseconds, in a manner that depends on the spin-lock RF field strength (Figure 2 E,F).

For completeness and illustration purposes, we also provide the analytical expression for the CCR rate constant under spin-lock, Γ 1 ρ CSA / D , in the framework of Redfield theory (which is only approximately valid 1 ):

Γ 1ρ CSA/ D =δ NH ⋅δ CSA ⋅P 2 (cos(θ)) [ 3 8 J (ω N )+ 1 12 ⋅J (ω RF -2 ω MAS )+ 2 12 ⋅J (ω RF -ω MAS )+ 1 12 ⋅J (ω RF +2 ω MAS )+ 2 12 ⋅J (ω RF +ω MAS )] (1)
Here, δ NH is the anisotropy of the dipolar coupling tensor, δ CSA is the chemical-shift anisotropy and P 2 (cos(θ)) is the second Legendre polynomial involving the angle between the principal axes of the dipolar and CSA interactions, θ . It is interesting to confront this equation to the expression of the R 1ρ relaxation rate constant due to the dipolar interaction and CSA, which is given as: 1,23

R 1ρ =R 1 + sin 2 θ eff [R 1 Δ - 1 2 R 1 ]
where

R 1 Δ ={ δ NH 4 } 2 [3⋅J (ω N )+ 1 3 ⋅J (ω RF -2 ω MAS )+ 2 3 ⋅J (ω RF -ω MAS )+ 1 3 ⋅J (ω RF + 2ω MAS )+ 2 3 ⋅J (ω RF +ω MAS )] + δ CSA 2 [ 1 3 ⋅J (ω RF -2ω MAS )+ 2 3 ⋅J (ω RF -ω MAS )+ 1 3 ⋅J (ω RF +2ω MAS )+ 2 3 ⋅J (ω RF +ω MAS )] (2)
Comparison of ( 1) and (2) reveals that the two relaxation rate constants provide similar information.

In this work we are interested primarily in motions slower than hundreds of nanoseconds (because in this range Γ 1 ρ CSA / D and R 1ρ are largest, as shown below); in this regime one can neglect the J(ω N ) term; furthermore, R 1 is very small compared to R 1ρ . Neglecting thus the J(ω N ) terms one can see that the two rate constants depend on the same spectral density functions (see Equations 1 and 2). We thus expect that the two rate constants are correlated, an observation that is indeed confirmed both experimentally and through simulations (discussed further below, Figures 2G and3C). In the remainder of this manuscript we do not use analytical calculations but use exclusively numerical simulations, because this allows to take into account additional factors, such as relaxation of the proton spin due to external mechanisms, as described below. By using numerical simulations we also circumvent the fact that Redfield theory is outside its range of validity for slow motions (microseconds). [START_REF] Schanda | Studying Dynamics by Magic-Angle Spinning Solid-State NMR Spectroscopy: Principles and Applications to Biomolecules[END_REF] Experimental detection of CSA/D cross-correlated relaxation

In order to experimentally verify the effects of cross-correlated relaxation in R 1ρ experiments, we prepared a sample of 2 H, [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N labelled microcrystalline ubiquitin. The protein was crystallized in a buffer that results in a re-protonation of exchangeable sites, such as the amide sites that we are interested in here, at a level of 60%. In such highly deuterated samples, narrow 1 H line widths can be achieved, and the simplicity of the spin system, that can be approximated as a two-spin system as far as [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N relaxation is concerned, also allows for straightforward analysis of relaxation data. [START_REF] Reif | Ultra-High Resolution in MAS Solid-State NMR of Perdeuterated Proteins: Implications for Structure and Dynamics[END_REF][START_REF] Schanda | Solid-State NMR Measurements of Asymmetric Dipolar Couplings Provide Insight into Protein Side-Chain Motion[END_REF][START_REF] Knight | Fast Resonance Assignment and Fold Determination of Human Superoxide Dismutase by High-Resolution Proton-Detected Solid-State MAS NMR Spectroscopy[END_REF] Of particular relevance here, narrow [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N lines are obtained even without 1 H decoupling, and therefore one can separately detect the decay of the doublet components. Figure 3 shows such [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N R 1ρ data for the [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N doublet components of amide sites in ubiquitin, obtained at a MAS frequency of 38 kHz, and a [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N RF field strength of 12 kHz. For the majority of residues the decay rate constants of the two components are detectably different. The differential relaxation rate constant of the two components,

Γ 1 ρ CSA / D ={R 1ρ ([N x H α ])-R 1ρ ([N x H β ]}/2
, varies across the protein sequence, ranging from close to zero to about 20 s -1 for different residues (Figure 3A). Residues located in secondary structure elements have small values, below 1 s -1 in most cases. Largest values of Γ 1 ρ CSA / D , exceeding 10 s -1 , are observed for residues G10 and K11 in the β1-β2 loop, which has been demonstrated before to undergo nanosecond motion in the microcrystals we are studying. [START_REF] Schanda | Quantitative Analysis of Protein Backbone Dynamics in Microcrystalline Ubiquitin by Solid-State NMR Spectroscopy[END_REF][START_REF] Haller | Amplitudes and Time Scales of Picosecond-to-Microsecond Motion in Proteins Studied by Solid-State NMR: A Critical Evaluation of Experimental Approaches and Application to Crystalline Ubiquitin[END_REF] Crosscorrelated relaxation rates are furthermore enhanced for residues located in the loop connecting β2 and ubiquitin's α-helix, the α-β3 loop and the loop connecting β4 to the short 3 10 helix. These loop regions have also elevated R 1ρ rate constants (Figure 3B). Figure 3C shows the correlation of residue-wise Γ 1 ρ CSA / D and R 1ρ rate constants. Because the majority of residues are located in secondary structure elements, residue-wise variations are generally small. Basically, only residues located in loop regions have significantly enhanced Γ 1 ρ CSA / D and R 1ρ rate constants. This limits the precision with which the correlation can be determined. The slope of the correlation is 1.48 when taking into account all residues (red line) and 1.96 when neglecting the two data points with largest rate constants. These experimental data unequivocally establish that the correlation has a non-zero offset: while R 1ρ rate constants exceed 2 s -1 for all residues, many residues have a close-to-zero CSA/D cross-correlated relaxation rate constant. This finding cannot be understood based on the simulations in Figure 2. To understand this offset, one needs to take into account partial proton selfdecoupling, as shown in the following section.

For completeness, we have also shown previously determined CSA/D CCR rate constants, determined with a spin-echo element instead of a spin-lock period. The drawback of the former experiment is that the decay with free evolution (R 2 ') is significantly faster than decay under spinlock (in ubiquitin by up to one order of magnitude). Consequently, it is more difficult to determine small CCR rates from the difference of the two comparably large decay rate constants corresponding to the two doublet components. Qualitatively, the three observables highlight the same residues undergoing slow motions.

It shall The correlation of the Γ 1 ρ CSA / D rate constant

Proton spin diffusion and its effect on heteronuclear transverse relaxation

Longitudinal relaxation of the proton spin, or spin diffusion (i.e. magnetization exchange between protons) leads to interconversion of the spin states, α and β. Accordingly, the states N x [H α ] and N x [H β ] inter-convert, and this would result in (partial) self-decoupling of these two [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N doublet components. In the context of cross-correlated relaxation, such self-decoupling is an important parameter, because it is expected to reduce the CCR rate constant. [START_REF] Hu | Interference between Transverse Cross-Correlated Relaxation and Longitudinal Relaxation Affects Apparent J-Coupling and Transverse Cross-Correlated Relaxation[END_REF] In this section we study the consequences of proton longitudinal relaxation or spin diffusion on heteronuclear relaxation and cross-correlated relaxation, using numerical simulations. At difference to the simulations shown in Figure 2, we added an ad hoc random field contribution acting on the proton spin. The effect of this random field is to relax the proton spin, and therefore interchange spin states H α and H β . In principle one might also attempt to explicitly simulate the proton spin diffusion or relaxation by using a large spin system with multiple protons. In practice, however, such simulations are prohibitively long if the spin system is larger than about 4-5 spins, and the accuracy of the proton evolution likely is not better than in our ad hoc external relaxation approach. Figure 4 also shows that the 1 H spin flips have an effect on the R 1ρ rate constants. For the range of R 1ρ rate constants that are observed experimentally, i.e. in the range of R 1ρ of a few s -1 to tens of s -1 at most, the proton relaxation leads to an increase in R 1ρ . The increase of R 1ρ depends on the RF field strength and the MAS frequency (Figure S4). The origin of this increased R 1ρ is the fact that stochastic proton spin flips counteract the MAS-induced recoupling of the dipolar interaction, as discussed in more detail in the Supporting Information. Taken together, the simulations of Figure 4 establish that the effect of the additional 1 H spin relaxation is (i) a reduction of the CCR rate constant, and (ii) an increase of the R 1ρ rate constant. The decreased CCR and increased R 1ρ rate constants change the slope of the correlation of the two rate constants, and result in an offset of the two rate constants (Figure 4D). It is remarkable that this situation is qualitatively found in the experimental data, shown in Figure 3C, which reproduce an offset of R 1ρ rate constants. It shall be noted that Figure 4 assumes a single amplitude and varies the time scale of the motion, a situation that is not likely to apply to the range of amide sites across a protein. Nonetheless, qualitatively the behavior observed in Figure 4 applies irrespectively of the precise motional model.

The amount by which the CCR rate is reduced depends on the rate at which the 1 H spin flip (Figure 4). We therefore need to estimate which spin diffusion rate constants are realistic in our experimental setting, i.e. at a MAS frequencies of 38 kHz and with the deuteration level we employed (60% back-exchange). Figure 5 shows experimental measurements of spin-diffusion between amide protons. By comparison of the spectra obtained with very short mixing time (100 μs, red) and 15 ms mixing (black), one can observe that there are additional peaks pointing to 1 H-1 H spin diffusion. Figure 5B highlights a few cases of well-resolved cross-peaks, corresponding to spin diffusion in a β-turn region (residues 45-47); for the cases of these cross-peaks, the corresponding amide sites are within ~4.2 Å at most. Figure 5C shows the time evolution of representative diagonal-and cross-peaks. The diagonal peaks show non-monoexponential decay, which we ascribe to contributions of spin-diffusion/cross-relaxation (fast component) and relaxation of the magnetization, presumably at the most dynamic sites, or through exchange with solvent. When fitting these curves with bi-exponential functions, respective rate constants are several tens per second for the fast rate (about 30-40% amplitude) and about 5-10 s -1 (60-70% amplitude), respectively. If attempting to describe the decay with a single rate constant, approximated by the time point where the signal decay to 1/e ~37%, apparent rate constants of 20 s -1 are obtained. The observable cross-peaks show a buildup within tens of milliseconds, similar to the fast rate observed in the diagonal-peak decay curves. Taken together, these data allow one to estimate that proton relaxation occurs at a rate of approximately 10-20 s -1 . (We note that the spin-diffusion rate is expected to decrease at higher MAS frequency and higher deuteration levels. A very recent report, which became available during the revision of the present manuscript, 30 used higher MAS frequencies (60 kHz), but lower deuteration (100% back-protonation of exchangeable sites), and values of about 10-20 s -1 were reported for the majority of the residues. While this is in rather good agreement with the present estimates, no apparent multi-exponential decay has been reported there, nor did the authors report on spin-diffusion cross-peaks.)

Using this estimate of 1 H relaxation, one can simulate its effect on measured Γ 1 ρ

CSA / D
and R 1ρ rate constants, and possibly correct for them. It has been proposed in solution-state NMR that 1 H-1 H NOE (i.e. relaxation-induced spin-diffusion) can be taken into account when analyzing crosscorrelated relaxation, [START_REF] Reif | Cross-Correlated Relaxation for the Measurement of Angles between Tensorial Interactions[END_REF] using an analytical approach. In principle, such an approach could also be applied in solids. However, the effect of the proton spin flips on the [START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N relaxation depends on the RF field and the MAS frequency, and an appropriate formalism is currently not developed. However, one may resort to numerical simulations, such as the ones we use here, to evaluate the effects of proton spin flips on heteronuclear R 1ρ or CCR rate constants. Experimentally determined proton decay rate constants, or at least global estimates, may be used to this end.

We illustrate this with the data shown in Figure 4. For example, assume that the experimentally determined R 1ρ rate constant is 3 s -1 , i.e. including the effects of proton relaxation. The actual R 1ρ rate constant, corrected for a 1 H relaxation rate constant of 20 s -1 (10 s -1 ), is 1.8 s -1 (2.6 s -1 ), i.e. the excess decay induced by the 1 H relaxation is of the order of 0.4 to 1.2 s -1 . In principle, one can take into account these effects, and use the corrected values when fitting motional models to transverse relaxation data.

A similar estimation for Γ 1 ρ CSA / D suggests, for example, that in the presence of a 1 H relaxation with a rate constant of 10 s -1 , a measured Γ 1 ρ CSA / D of 1 s -1 points to an actual Γ 1 ρ CSA / D rate constant of 3.6 s -1 . Arguably, it is more difficult to correct Γ 1 ρ CSA / D for 1 H relaxation effects, because the values of this rate constant are smaller than those for R 1ρ , and 1 H relaxation has a stronger effect on

Γ 1 ρ CSA / D than on R 1ρ .
It is noteworthy that such a correction appears more difficult for the CCR rate constant, which may become very small and thus error-prone, than for R 1ρ . Of course, one may attempt to suppress the effects of proton spin flips even more, by using higher MAS frequencies and/or higher deuteration levels. Increasing MAS frequency will, however, not suppress proton magnetization decay which is actually due to dynamics or chemical exchange, such that it appears unlikely that one may suppress its effects completely.

Suppression of the effect of CSA/D CCR in R 1ρ experiments

Finally, we note that the presence of CSA/D CCR might hamper accurate extraction of R 1ρ rate constants. In solution-state NMR, schemes for suppression of the effect of CSA/D CCR have been reported. [START_REF] Korzhnev | An NMR Experiment for the Accurate Measurement of Heteronuclear Spin-Lock Relaxation Rates[END_REF][START_REF] Massi | NMR R1rho Rotating-Frame Relaxation with Weak Radio Frequency Fields[END_REF] They generally require irradiation of 1 H spins. Suppression of CCR in R 1ρ experiments, even in solution-state NMR, is complicated by the fact that 1 H RF fields applied during heteronuclear spin-lock may partially re-introduce the scalar coupling Hamiltonian in solution and solids, [START_REF] Massi | NMR R 1ρ Rotating-Frame Relaxation with Weak Radio Frequency Fields[END_REF][START_REF] Korzhnev | An NMR Experiment for the Accurate Measurement of Heteronuclear Spin-Lock Relaxation Rates[END_REF] and it may also reintroduce the dipolar coupling in the solid state. In MAS ssNMR particular care has to be taken to avoid rotary resonance conditions, i.e. [START_REF] Schanda | Studying Dynamics by Magic-Angle Spinning Solid-State NMR Spectroscopy: Principles and Applications to Biomolecules[END_REF] H multi-pulse or CW schemes which counteract the averaging brought about by MAS. To investigate these effects, we have tested two different suppression schemes, namely (i) a single 1 H π pulse, applied in the center of the relaxation period and (ii) a train of 1 H π pulses, spaced by 5 ms delays. Both schemes, or variants thereof, are used in solution-state NMR. [START_REF] Massi | NMR R 1ρ Rotating-Frame Relaxation with Weak Radio Frequency Fields[END_REF][START_REF] Korzhnev | An NMR Experiment for the Accurate Measurement of Heteronuclear Spin-Lock Relaxation Rates[END_REF][START_REF] Palmer | Characterization of the Dynamics of Biomacromolecules Using Rotating-Frame Spin Relaxation NMR Spectroscopy[END_REF] Interestingly, we find that the multi-pulse scheme, even if the 1 H pulse spacing (5 ms) is very long compared to the rotor period (~20-30 μs), the heteronuclear coherences decay significantly faster than in the absence of 1 H pulses or with a single 1 H π pulse (Figure S5). The use of a single 1 H π pulse appears as a straightforward method to suppress the effects of CCR. We note, however, that even with a suppression of CCR, one does not expect mono-exponential decay, due to the inherent multi-exponential nature of decays in solids, outlined above.

Conclusions

We have shown here that CSA/D cross-correlated relaxation can be measured in MAS ssNMR of deuterated proteins. This relaxation rate constant reports on dynamics on nanosecond-tomicrosecond time scales, as we illustrate with an application to microcrystalline ubiquitin. In many practical applications, however, the more easily accessible R 1ρ rate constant of the decoupled line has advantages over the measurement of CSA/D CCR under spin-lock, in particular due to higher sensitivity and resolution in a decoupled vs a non-decoupled spectrum. We also provided evidence, by numerical simulations, that rapid 1 H spin flips, as mediated by 1 H-1 H spin diffusion, reduces the CCR rate constant, and also contributes to the observed R 1ρ rate constant. This finding also shows that for accurate quantitative analyses of R 1ρ rate constants one should correct the measured rate constants for the contribution arising from proton relaxation. The present report of cross-correlated relaxation adds to the increasing arsenal of techniques to study protein dynamics on long time scales by MAS ssNMR, and stresses the importance of considering hitherto ignored effects when interpreting relaxation rate constants in ssNMR.

Figures showing the pulse sequences used in this study; figures with additional simulations, and experimental data; extended discussion about the increase of R 1ρ rate constants due to proton spin relaxation. This material is available free of charge via the Internet at http://pubs.acs.org.

Figure 4

 4 Figure4shows the results of numerical simulations of the heteronuclear R 1ρ experiment in the presence of this additional proton relaxation, simulated here for the case of small-amplitude motion of the HN bond. Simulations with larger-amplitude motions are shown in FigureS3. These simulations establish that proton longitudinal relaxation induced by external relaxation mechanisms indeed reduces the CCR rate constant. The reduction of the Γ 1 ρ CSA / D rate constant depends on the rate constant of the proton relaxation, as can be seen by comparison the red curves shown in panels (A) to (C). This reduction of Γ 1 ρ CSA / D can be understood by the averaging of the two components.

Figure 1 :

 1 Figure 1: Simulated evolution of 15 N spin coherences under a spin-lock RF field. (A) Time evolution of in-phase N x coherence (red) and a single-exponential fit, which clearly results in poor agreement with the decay. For the simulation shown in grey, a 1 H π pulse was inserted in the center of the spin-lock period. An exponential fit is shown in black. (B) Simulation of the doublet components, N x -[H α ], N x -[H β ] (red), and exponential fits (black). The jump angle was 20 degrees and the jump rate constant 30000 s -1 .

Figure 2 :

 2 Figure 2: Simulated relaxation rate constants of a 15 N spin due to its CSA tensor and 1 H-15 N dipolar coupling, within the two-site jump model. Shown are rate constants for the individual doublet components (A, B), the average rate constant R 1ρ (C) and the CSA/D CCR rate constant (D), as a function of the jump angle and jump rate constant. (E) shows R 1ρ and CCR rate constants for a constant jump angle (20º) as a function of the jump rate constant, assuming a MAS frequency of 40 kHz and a RF field strength of 12 kHz. (F) uses parameters as in (E), but with an RF field strength of 32 kHz. (G) Correlation of R 1ρ and CCR rate constants, assuming a jump angle of 5º and a MAS frequency of 40 kHz and a RF field strength of 12 kHz. The jump angle can be related to the often-used order parameter S 2 as S 2 =1/4(3 cos 2 (θ)+1) . 34 According to this equation, jump angles of 5, 10, 15, 20, 25 degrees, for example, correspond to S 2 = 0.994, 0.977, 0.95, 0.91, 0.866, respectively.

Figure 3 :

 3 Figure 3: Experimental observation of CSA/D CCR in[START_REF] Igumenova | Assignments of Carbon NMR Resonances for Microcrystalline Ubiquitin[END_REF] N spin-lock relaxation experiments in deuterated, 60% back-exchanged ubiquitin, obtained at a MAS frequency of 38 kHz. (A) Differential relaxation rate constants as a function of the residue number, determined from the relaxation rate constant of the two doublet components. (B) R 1ρ relaxation rate, i.e. obtained from the average of the two rate constants. In (C), a correlation of the two rate constants is shown. The red line corresponds to a correlation obtained over all available residues, with a slope of 1.48 and an intercept of 2.74 s -1 . The insert focuses on the data points at low rate constants. The blue line has a slope of 1.96 s -1 . Panel (D) shows previously published[START_REF] Schanda | Quantitative Analysis of Protein Backbone Dynamics in Microcrystalline Ubiquitin by Solid-State NMR Spectroscopy[END_REF] CSA/D CCR relaxation rate constants, obtained with a spin-echo (delay -π pulse -delay) sequence instead of a spin-lock, at the same B 0 field strength. The insert shows the correlation of these data with the R 1ρ -based CCR rate constants of panel (A).

Figure 4 :

 4 Figure 4: The effect of proton relaxation on 15 N transverse relaxation parameters, from numerical simulations with an ad hoc added random field fluctuation leading to proton relaxation. Shown are the simulated relaxation rate constants of N x (i.e. the relaxation with a 1 H π pulse in the center of the spin-lock relaxation period), R 1ρ , and the cross-correlated relaxation rate constant, Γ 1 ρ CSA / D , obtained from the differences of relaxation of the two doublet components. Different rate constants of the 1 H relaxation are assumed: 0 s -1 (A), 5 s -1 (B) and 20 s -1 (C); R 1ρ and Γ 1 ρ CSA / D rate constants are shown in black and red, respectively. Panel (D) correlates the R 1ρ and Γ 1 ρ CSA / D relaxation rate constants for different 1 H relaxation rate constants. The jump angle was set to 5 degrees, corresponding to an order parameter of S 2 =0.994. More simulations with a large motional amplitude are shown in Figure S3.

Figure 5 :

 5 Figure 5: Experimental investigation of proton relaxation and spin diffusion in the deuterated, 60% back-exchanged ubiquitin sample at 38 kHz MAS (14.1 T). (A) Pulse sequence used to detect 1 H relaxation/spin diffusion. (B) Overlay of two spectra obtained with different mixing times as indicated. Selected amide peaks, and the corresponding spin-diffusion cross-peaks are indicated, and the distances between these sites are indicated on the right. (C) Representative time traces observed at peak positions corresponding to amide sites (diagonal peaks), and two spin-diffusion cross-peaks, obtained at spectral positions at which intensity is observed only with a mixing time. The grey lines correspond to bi-exponential fits, and the fitted rate constants are indicated in brackets. The identity of the peaks for which the graphs are shown in C is indicated in the spectrum with (i)-(iii).

Acknowledgements

This work was financially supported by the European Research Council (ERC-StG-311318 ProtDyn2Function) and the Commissariat a l'énergie atomique et aux énergies alternatives (CEA), the Centre National de la Recherche Scientifique (CNRS) and the Université Grenoble Alpes. This work used the platforms of the Grenoble Instruct Center (ISBG; UMS 3518 CNRS-CEA-UJF-EMBL) with support from FRISBI (ANR-10-INSB-05-02) and GRAL (ANR-10-LABX-49-01) within the Grenoble Partnership for Structural Biology (PSB). We thank Matthias Ernst for stimulating discussions and for assistance with implementing the proton random field fluctuations in the GAMMA simulations.

Supporting Information Description

TOC Figure