jp076174l

2782

J. Phys. Chem. B 2008, 112, 2782-2794

Diffusion NMR for Determining the Homogeneous Length-Scale in Lamellar Phases Ingrid Åslund,* Celia Cabaleiro-Lago, Olle So1 derman, and Daniel Topgaard Physical Chemistry 1, Lund UniVersity, P.O.B. 124, SE-221 00 Lund, Sweden ReceiVed: August 2, 2007; In Final Form: December 10, 2007

The size of the anisotropic domains in a lyotropic liquid crystal is estimated using a new protocol for diffusion NMR. Echo attenuation decays are recorded for different durations of the displacement-encoding gradient pulses, while keeping the effective diffusion time and the range of the wave vectors constant. Deviations between the sets of data appear if there are non-Gaussian diffusion processes occurring on the time-scale defined by the gradient pulse duration and the length-scale defined by the wave vector. The homogeneous length-scale is defined as the minimum length-scale for which the diffusion appears to be Gaussian. Simulations are performed to show that spatial variation of the director orientation in an otherwise homogeneous system is sufficient to induce non-Gaussian diffusion. The method is demonstrated by numerical solutions of the Bloch-Torrey equation and experiments on a range of lamellar liquid crystals with different domain sizes.

Introduction Molecular motion in many natural and synthetic materials, such as wood, brain tissue, and mesoporous silica, is anisotropic as a result of the ordered arrangement of the anisotropic structural elements. The experimentally measured diffusion anisotropy is influenced by the spatial heterogeneity of the preferred direction of the structure. This problem is of special relevance for diffusion tensor imaging (DTI) of the human brain,1,2 where the measured anisotropy depends on the pixel size in the image. With low spatial resolution it is more likely that the signal from one pixel contains contributions from nerve fibers with many different directions, resulting in an apparent low anisotropy. In the limit of a highly disordered structure, such as brain gray matter,3 the diffusion appears to be isotropic. In order to study the structure of such materials, we require methods that are sensitive to the anisotropy and where the length-scale defined by the experiment can be varied. Lamellar liquid crystals are not only technically important for a range of applications, e.g., as soaps or templates for mesoporous inorganic materials,4 but also convenient model systems when developing methods for investigating diffusion anisotropy. A schematic representation of the different lengthscales of a multidomain lamellar phase is shown in Figure 1. On the nanometer scale, the structure consists of alternating layers of surfactants and water. Unless special measures have been taken to prepare a uniform domain orientation, the sample is composed of an ensemble of randomly oriented domains. The domain size depends on the details of the sample preparation,5 but is usually in the range of 1-100 µm. On a macroscopic scale, the sample appears to be isotropic. When studying a system with a hierarchy of structural lengthscales and dynamic time-scales, it is necessary to consider the length- and time-scales defined by the measurement technique. As an example, quadrupolar splittings appear in the 2H NMR spectrum of 2H2O if the molecules are located in an anisotropic domain of uniform orientation for a time longer than the inverse of the value of the splitting (in Hz). A typical lamellar phase gives rise to a 1 kHz splitting corresponding to a 1 ms time* Corresponding author. E-mail: [email protected].

scale and 2 µm length-scale, assuming a diffusion coefficient of 2 × 10-9 m2/s. If the domains are much smaller than 1 µm, the sample appears to be isotropic. 2H2O exchange6,7 and diffusion-diffusion exchange spectroscopy (DEXSY)8 are two methods that have been devised in order to study the time-scale for exchange between domains of different orientation. The distribution of domain orientations gives rise to a distribution of the NMR observables: 2H splittings and diffusion coefficients for the 2H2O exchange and DEXSY experiments, respectively. In 2D correlation experiments, where the two observation periods are separated by an adjustable mixing time, off-diagonal peaks appear if significant exchange between domains of different orientation has occurred. The characteristic time for exchange can be estimated by varying the mixing time and observing when the off-diagonal features appear. The estimated time can then be converted to a length if the diffusion coefficient is known. For technical reasons, the diffusion NMR experiment invariably operates on time-scales longer than about 1 ms and lengthscales larger than 1 µm. When analyzing or modeling diffusion NMR experiments, it is thus sufficient to quantify the material properties, for example the anisotropy, on a coarse-grained length-scale (see Figure 1).9 With further coarse-graining, the material will start to appear homogeneous and isotropic. In analogy with Spowart et al.,10 we denote the length for this latter transition as the homogeneous length-scale λhom. Here, we describe an experimental protocol for determining λhom using the pulsed-field-gradient spin echo (PFG SE) or stimulated echo (PFG STE) experiment.11-14 Two time-scales and one length-scale, all of which can be varied over wide ranges, are defined by the diffusion NMR experiment. In the case of isotropic Gaussian diffusion, the detected NMR signal is given by the Stejskal-Tanner equation,15 incorporating the effects of all of the experimentally defined scales. Deviations from this expression can be observed if any of the structural or dynamical scales of the system coincide with the ones defined by the diffusion NMR experiment. By a proper choice of experimental parameters, one can determine the minimum length-scale at which the system appears homogeneous. The new protocol is based on previous work investigating the effect

10.1021/jp076174l CCC: $40.75 © 2008 American Chemical Society Published on Web 02/14/2008

Homogeneous Length Scale in Lamellar Phases

J. Phys. Chem. B, Vol. 112, No. 10, 2008 2783

Figure 2. PFG STE NMR pulse sequence. The 90° RF pulses produce a stimulated spin echo a time 2τ1+τ2 after the first 90° pulse. The echo amplitude E is encoded for translational motion using pairs of magnetic field gradient pulses with strength and direction given by the vector G. The effective diffusion time td is given by ∆ - δ/3, where ∆ is the time between the onset of the gradient pulses, and δ is the duration of one gradient pulse. For experimental series with different δ, the values of ∆ and G are varied simultaneously in order to keep td and q ) γGδ/2π constant. f(t) is the time modulation of the effective gradient.

Figure 1. Illustration of the different length-scales of a multidomain lamellar phase. Top: isotropic on a macroscopic scale. Middle: anisotropic on the micrometer length-scale. The colors symbolize different orientations of the anisotropic domains. Bottom: nanometer length-scale lamellar structure with alternating layers of surfactant (yellow slabs/red balls symbolizing oily interior/hydrophilic headgroups) and water (not shown). The homogeneous length-scale λhom denotes the length-scale of coarse graining required to make the material appear homogeneous.

of the durations of the displacement-encoding gradient pulses.16-20 In addition to experiments, an extensive set of simulations is performed to study the diffusion properties of idealized geometric models of, or formal analogs to, polydomain lamellar phases.

where γ is the magnetogyric ratio of the observed nucleus, and D(r) is the local diffusion tensor. The second term of eq 2 is written in a form that makes m(r,t) continuous, while all the information of the transport properties of the system is given by the not necessarily continuous functions F(r) and D(r).17 This form is particularly useful when numerically solving the partial differential equation for the unknown function m(r,t). Approximations concerning either the system or the gradient modulation are necessary in order to solve eq 2. Assuming a homogeneous system, i.e., that D and F are independent of r, M(t) can be shown to be

M(t) ) exp[-γ2G‚D‚G

M(t) )

∫F(r)m(r,t)dr

(1)

The signal is encoded for motion using magnetic field gradients where the amplitude and direction is given by the vector G and the time dependence is contained in the function f(t), as shown in Figure 2. Neglecting nuclear relaxation processes, m(r,t) evolves according to the Bloch-Torrey equation:21,22

∇‚[F(r)D(r)‚∇m(r,t)] ∂m(r,t) ) - iγG‚rf(t)m(r,t) + ∂t F(r)

(3)

at the time of echo formation. For the Stejskal-Tanner version of the diffusion NMR experiment,15 f(t) is

f(t) ) -H(t) + H(t - δ) + H(t - ∆) - H(t - ∆ - δ)

(4)

where H is the Heaviside step function. The echo attenuation E is the ratio between the signals with and without applied gradients. For this particular experiment, E is defined as

Theoretical Basis The detected NMR signal at the time t is proportional to the macroscopic complex transverse magnetization M(t), which can be obtained from a volume integral of the local complex transverse magnetization m(r,t) weighted by the normalized local spin density F(r):

∫0t (∫0t′ f(t′′)dt′′)2dt′]

E)

M(∆ + δ) MG)0(∆ + δ)

(5)

Combining and evaluating eqs 3, 4, and 5 yields

E ) exp[-γ2G‚D‚Gδ2(∆ - δ/3)]

(6)

This equation is more compactly expressed as

E ) exp(-b:D)

(7)

In eq 7, “:” implies a generalized scalar product evaluated with

b:D )

∑R ∑β bRβDRβ; R,β ) x,y,z

(8)

bRβ ) γ2GRGβδ2td

(9)

in which

(2)

2784 J. Phys. Chem. B, Vol. 112, No. 10, 2008

Åslund et al.

and td is the effective diffusion time:

td ) ∆ - δ/3

DRβ(t) ) (10)

The standard way of evaluating diffusion in isotropic Gaussian systems is by fitting the isotropic version of eq 7, E ) exp(bD), to the experimental data. A plot of ln E vs b is often referred to as a Stejskal-Tanner plot. For heterogeneous systems, where F and D are functions of r, it is convenient to quantify the motion using the propagator P(r0,r,t), which is the probability density that a molecule starting at position r0 moves to position r during the time t. In order to facilitate numerical calculations, we choose to normalize P(r0,r,t) with F(r). On account of this normalization, our equations are slightly different from the ones used in standard texts.12 In the so-called short gradient pulse (SGP) approximation, i.e., when δ f 0 while δG ) constant, the gradient pulses define the start and end of the diffusion time, and t ) ∆. In this case, the echo attenuation is written as23

E)

∫∫F(r0)F(r)P(r0,r,∆) exp[-i2πq‚(r - r0)]drdr0

(11)

where the wave vector q is

q ) γGδ/2π

(12)

The evolution of P(r0,r,t) is calculated by solving21

∂P(r0,r,t) ∇‚[F(r)D(r)‚∇P(r0,r,t)] ) ∂t F(r)

(13)

1 1 D(t) ) trace{D(t)} ) [Dxx(t) + Dyy(t) + Dzz(t)] 3 3

Within the SGP approximation the diffusion time t equals ∆, implying that 〈R2〉 and D(t) can be determined from the initial, low-q, decay of E irrespective of the functional form of P(R,∆). For the purpose of this paper, it is convenient to define an isotropic apparent diffusion coefficient D(td,δ) depending on both td and δ:

D(td,δ) ) -

∂ ln E(G,δ,td) 1 lim Gf0 γ δ td ∂G2 2 2

(21)

D(td,δ) as defined in eq 21 equals D(t ) td) if the SGP approximation is valid or if the propagator is Gaussian. A dependence of D(td,δ) on δ implies that the diffusion is nonGaussian on the time-scale of δ, which is often an order of magnitude smaller than td. Evaluating experimental data using eq 21 is therefore a simple way of extending the effective timescale for detecting non-Gaussian, restricted, diffusion. Three length-scales are defined by the diffusion NMR experiment: the diffusion length during the gradient pulse lδ, the diffusion length during the diffusion time lt, and the inverse of the wave vector q-1. The diffusion lengths are given by

(14)

In eq 14, δ denotes the Dirac delta function (not to be confused with the gradient pulse length). Just as in eq 2, the chosen normalization makes the unknown function P(r0,r,t) continuous throughout the system. Solutions to eq 13 exist for a series of simple geometries such as cylinders and spheres.24-26 Equation 11 can be expressed as23

∫P(R,∆) exp(-i2πq‚R)dR

(15)

where R ) r - r0 is the displacement, and P(R,∆) is the average propagator

∫F(r0)F(R + r0)P(r0,R + r0,∆)dr0

(16)

Insertion of an anisotropic Gaussian average propagator into eq 15 reproduces eq 6, the only difference being that ∆ - δ/3 is replaced with ∆ as required by the SGP approximation. A series expansion of eq 15 shows that

lim E(q,∆) ) exp(-2π2q‚〈R2〉‚q) qf0

(17)

where 〈R2〉 denotes the mean-square displacement tensor with the elements

〈RRRβ〉 )

(20)

lδ ) x2D(t ) δ)δ

P(r0,r,0) ) δ(r - r0)

P(R,∆) )

(19)

The isotropic average D(t) is calculated from

with the initial condition

E(q,∆) )

〈RRRβ〉 2t

∫RRRβ P(R,t) dR

lt ) x2D(t ) td)td

while q was defined in eq 12. Deviations from Gaussian behavior will be observed if any length-scale of the heterogeneities of the system coincides with any of the NMR scales. For a rough estimate of lδ and lt, it is sufficient to replace D(t ) δ) and D(t ) td) with any value between the short- or long-t limits of D(t). Modeling Diffusion in Multidomain Liquid Crystals For liquid crystals, D(r) can be quantified with the local director n(r), using spherical coordinates expressed with the angles θ(r) and φ(r) against the direction of the magnetic gradient, and the diffusion coefficients parallel and perpendicular to the director, D|(r) and D⊥(r). Each domain orientation gives rise to a specific diffusion coefficient Dθ according to

Dθ ) D| cos2 θ + D⊥ sin2 θ

(23)

In an idealized lyotropic liquid crystal, we may assume that the density and the principal components of the diffusion tensor are constant throughout space. In this case, it is sufficient to specify the spatially constant values of D|| and D⊥ and the local values of θ(r) and φ(r). The values of D| and D⊥ are reported using the anisotropy η and the isotropic average 〈D〉 given by

(18) η)

〈R2〉

quantifies the width of P(R,∆). The elements of the timedependent apparent diffusion coefficient27 tensor D(t) are defined as

(22)

and

D| D⊥

(24)


D ⊥ + D| 2

(25)

2D⊥ + D| 3

(26)

〈D〉 ) 〈D〉 )


Equations 25 and 26 are valid for 2D and 3D systems, respectively.28 For typical multidomain liquid crystals, the direction of n(r) varies on the micrometer length-scale, thus giving rise to variations of the observed Dθ on the same length-scale. As will be demonstrated with numerical calculations below, such variations lead to non-Gaussian diffusion. The length-scale of the spatial variation of the domain orientations can be quantified by the director correlation function C(R):

C(R) )

∫n(r + Rn)‚n(r)dr

(27)

The functional form of C(R) for a liquid crystal composed of well-defined “crystallites” is markedly different from a system where the director reorientation as a function of position occurs smoothly. In the simple case of an exponential correlation

C(R) ) exp(-R/λp)

(28)

the decay constant λp, the persistence length, is a natural choice for the domain size.6 For the more realistic case where there is no simple functional form for C(R), the definition of the domain size is more problematic. One solution is to define the domain size as the maximum value of R where C(R) exceeds a certain threshold value. In analogy with the exponential correlation, this threshold value could be chosen as e-1, but it is important to realize that this particular choice is not more “natural” than any other value such as, say, 0.5 or 0.01. Having these considerations in mind, the homogeneous length-scale λhom can be defined as the value of R where C(R) has dropped to a threshold value close to zero. If the properties of the material are averaged on a length-scale larger than λhom, the material appears homogeneous and isotropic. We can separate several regimes for the analysis of the diffusion NMR experiment. In the limit of large anisotropic domains, i.e., λhom is much larger than all the NMR lengthscales lδ, lt, and q-1, multi-Gaussian diffusion behavior is observed. Integrating the contributions from an ensemble of randomly oriented domains gives rise to a multiexponential decay of the diffusion NMR signal:29

E(b) )

∫0π/2 sin θ exp[-b(D| cos2 θ + D⊥ sin2 θ)]dθ

(29)

with the initial slope given by

lim E(b) ) exp(-b〈D〉) bf0

(30)

In the limit of small domains, where λhom is much smaller than all of lδ, lt, and q-1, the signal decay is single exponential with a decay constant D∞:

E(b) ) exp(-bD∞)

(31)

The value of D∞ merits special attention. Naively, one would assume that D∞ ) 〈D〉, and this would also be the result of the Ka¨rger model30 describing exchange between multiple domains with Gaussian diffusion. Now consider a finite size domain with fast diffusion surrounded by domains with slow diffusion. The observed diffusion coefficient for the fast domain will start to

Figure 3. 2D models of a polydomain lamellar liquid crystal composed of square domains arranged in a checkerboard pattern with (a) sharp or (b) smooth boundaries between adjacent domains. The spin density F and the trace of the diffusion tensor 〈D〉 are constant throughout the x,y-plane, while the diffusion in the x- and y-directions, Dx and Dy, vary according to the figures. The grayscale corresponds to (Dx - Dy)/ 〈D〉 (i.e., white: fast x; gray: isotropic; black: fast y). The blue arrows indicate the direction of maximum diffusivity. The lengths of the arrows are scaled by |(Dx - Dy)/〈D〉| (i.e., no arrow: isotropic). d is the size of the unit cell in the periodic pattern. The anisotropy η as defined in eq 24 corresponds to the anisotropy in the center of each crystallite.

decrease when the diffusion displacement approaches the size of the domain, leading to the general result D∞ < 〈D〉. Below, this effect is demonstrated by numerical solution of eq 13 for 2D models of polydomain liquid crystals. The presence of restrictions in the boundaries between the differently oriented domains would lead to further reduction of D∞.31 The intermediate length-scale, i.e., λhom, is about the same as any of lδ, lt, and q-1 and is more difficult to analyze, and we do not expect to find a general expression for the echo intensity. However, some features are quite independent of the exact geometry, and we will utilize this fact to experimentally determine λhom using a method that first will be demonstrated by numerical calculations for a 1D system designed to have diffusion properties qualitatively similar to the ones for a 3D lamellar phase. If we for a moment forget about the technical limitations of the diffusion NMR experiment, we would expect the following variations of D(t) as the diffusion time is increased. In the limit where the diffusion displacement is much smaller than the width of the water domains in the lamellar phase, D(t) would be close to the value for bulk water D0. Note that this regime would require diffusion times smaller than 1 µs, which is way out of the range of accessible time-scales for diffusion NMR. Increasing t to the regime where the displacements are much larger than the lamellar spacing, but still smaller than the domain size, leads to D(t) ) 〈D〉. For an idealized lamellar phase, D| ) D0 and D⊥ ) 0, resulting in 〈D〉 ) 2D0/3 according to eq 25. Further


Åslund et al.

increase of t leads to a transition to the long-time limit D(t) ) D∞. The ratio 〈D〉/D∞ depends on the details of the spatial arrangement of the domains, the presence of barriers at the grain boundaries, and the nature of these potential barriers. Diffusion Propagator in a 2D Model of a Liquid Crystal. Simple 2D models of polydomain liquid crystals were constructed by arranging square domains in the x,y-plane as shown in Figure 3, and the diffusion propagator was estimated using a finite difference method to numerically solve eq 13. The systems are discretized in a 2D rectangular grid of nodes labeled with index m,n. The distances between adjacent nodes in the x- and y-directions are ∆x and ∆y. In order to estimate an average propagator according to eq 16, the calculation is repeated for a series of starting positions (x0,y0) labeled with a common index i. The (x,y) interval is defined with (x0,y0) in the center. The index m,n thus corresponds to a certain value of the displacement (x - x0,y - y0). The displacements (x x0)m,n and (y - y0)m,n, the diffusion coefficients in the x- and x y and Dm,n,i , and the equilibrium spin density y-directions Dm,n,i Fm,n,i are specified for each node. The evolution of the system is calculated for a series of time steps tj labeled with index j and with the duration ∆t. For the first time step, j ) 1, the propagator equals Pm,n,i,j ) 1/(∆x∆y) for the node with position (x - x0 ) 0, y - y0 ) 0), while being zero for all others. In each time step, a new propagator Pm,n,i,j+1 is calculated by

Pm,n,i,j+1 ) x x y y - Jm,n,i,j + Jm,n-1,i,j - Jm,n,i,j (32) Pm,n,i,j+1 + Jm-1,n,i,j

where

Pm,n,j∆x∆y ) 1 ∑n ∑ m

In order to investigate the shape of the propagator in the short-, intermediate-, and approaching the long-time limits, a first series of simulations were performed with a fixed starting position, (x0 ) d/6, y0 ) d/8). The system was discretized in 64 × 64 square array of nodes, and the evolution of the system was followed in 500 time steps. The starting position was located in the center of the system. The diffusion coefficient was chosen to give (2〈D〉tmax)1/2/(xmax - xmin) ) 0.13. Three values of the repeating unit d were used: 512∆x, 32∆x, and 8∆x. A second series of simulations were performed to investigate the apparent diffusion coefficient as a function of time and anisotropy. In this case, the system was specified in a 128 × 128 grid and followed over 2200 time steps. The repeating unit was covering 16 × 16 grid points, and the average propagator was obtained by averaging the propagators for all these 256 starting positions according to eq 35. The diffusion coefficient was chosen to give (2〈D〉tmax)1/2/d ) 1.0. The two geometries shown in Figure 3 were studied for three values of the anisotropy: η ) 1, 3, and 10 (as defined in eq 24). The simulations were implemented in Matlab running on a Mac PowerPC G5 with dual 2.5 GHz processors and 1 GB RAM. The first series of simulations required less than 15 min of computer time, while the second series took about 20 h. 1D Analog of a Liquid Crystal. As a suitable analog to a liquid crystal, we choose a system where the diffusion coefficient D(z) varies with position z according to

D(z) ) Dmin + (Dmax - Dmin )

x ) Jm,n,i,j

(

)

x x + Fm+1,n,iDm+1,n,i ) Pm+1,n,i,j Pm,n,i,j ∆t(Fm,n,iDm,n,i (33) 2 Fm+1,n,i Fm,n,i 2(∆x∆y)

is the flow from node m,n to m+1,n, and y Jm,n,i,j )

(

)

y y + Fm,n+1,iDm,n+1,i ) Pm,n+1,i,j Pm,n,i,j ∆t(Fm,n,iDm,n,i (34) 2 Fm,n+1,i Fm,n,i 2(∆x∆y)

is the flow from node m,n to m,n+1 at time step j. The average propagator Pm,n,j is evaluated by summing

Pm,n,j )

∑i F0i Fm,n,iPm,n,i,j ∑i

(35)

where F0i ) F(x0,y0) at the starting position (x0,y0) with index i. For each time step, the apparent diffusion coefficient D(tj) is calculated with

D(tj) )

∑ ∑ (x - x0)m,nPm,n,j∆y∆x m n 2tj

(36)

The denominator in eq 35 and the choice of initial value of the propagator ensure that it is normalized according to

cos(2πz/d) + 1 2

(38)

where d is the repeating distance of the periodic structure. The minimum and maximum diffusion coefficient Dmax and Dmin correspond to the diffusion along and across the lamellar planes, respectively. Note that this system should not be seen as a geometric model of a lamellar phase, but rather a useful analog32,33 for which the long- and short-scale limits can be calculated analytically as shown in the Appendix. In the shorttime limit, E(b) is multiexponential:

E(b) )

1 d

d/2 exp[-bD(z)]dz ∫-d/2

(39)

where the initial, low-b, decay is given by

〈D〉 )

F0i Fm,n,i

(37)

1 d

d/2 D(z)dz ) ∫-d/2

Dmax + Dmin 2

(40)

In the long-time limit, E(b) is single exponential with the decay constant

D∞ )

[∫ 1 d

d/2

-d/2

1 dz D(z)

]

-1

) xDmaxDmin

(41)

The similarity between the expected outcome of the diffusion NMR experiment for a true 3D liquid crystal and the 1D analog is quite evident. For the simple 1D model, numerical solutions are sufficiently fast to permit careful investigations of the multidimensional parameter space, which is required to properly study the intermediate time- and length-scales. The Bloch-Torrey equation can be solved analytically only for a small set of special cases, thus making the use of numerical calculations more or less inevitable when theoretically examin-



Figure 4. 1D analog of a polydomain lamellar phase. The diffusion coefficient D(z) varies with position z according to eq 38 using the parameters maximum diffusion coefficient Dmax ) 2 × 10-9 m2/s, minimum diffusion coefficient Dmin ) 2 × 10-10 m2/s, and repeating distance d ) 50 µm. Edge effects are minimized by calculating the evolution of the entire system, but only recording the NMR signal from the central unit cell (shown with crosses).21

ing diffusion NMR. Here we use a finite difference method based on the works of Hwang et al.34 and Novikov et al.35 The simulated system consists of a 1D grid with distance ∆z between adjacent nodes labeled with index i. The position zi, the diffusion coefficient Di, and equilibrium spin density Fi are specified for each node. The evolution of the system is calculated for a series of time steps with duration ∆t. The NMR pulse sequence, as shown in Figure 2, is defined by the value of the gradient modulation fj for each time step labeled with index j. For the first time step, j ) 1, the complex transverse magnetization equals mi,j ) 1 for all i. For each successive time step, mi,j is calculated by first rotating the magnetization under the action of the gradient

m′i,j ) mi,j exp(iγGfjzi∆t)

(42)

and then applying the effect of diffusion

mi,j+1 ) m′i,j + Ji-1,j - Ji,j

(43)

where

Ji,j )

(

)

∆t(FiDi + Fi+1Di+1) m′i+1,j m′i,j Fi+1 Fi 2(∆z)2

(44)

is the flow of magnetization from node i to i+1 at time step j. The echo attenuation E is calculated for the final time step j ) N using the following formula:

E)

∑ Fimi,j)N ∑ Fi

(45)

where the summation is performed over all the nodes belonging to the central unit cell of the system (see Figure 4). For all simulated systems, the NMR pulse sequence parameters were the same: 64 values of q logarithmically spaced between 2.5 × 103 and 2.5 × 105 m-1, 32 values of td logarithmically spaced between 4 and 400 ms, and 32 values of δ logarithmically spaced between 1 ms and 1.5 times the value of td (corresponding to ∆ ) δ; cf. eq 10). This gives in total 65536 different parameter combinations for each system. The first series of simulations was performed for systems with Dmax ) 2 × 10-9 m2/s, Dmin ) 2 × 10-10 m2/s, and three values of d (10, 20, and 50 µm). In the second series of systems, Dmin was varied (2 × 10-10, 4 × 10-10, and 1 × 10-9 m2/s), while keeping Dmax and d constant at 2 × 10-9 m2/s and 20 µm. The values of ∆t and ∆z were set to 0.2 ms and 1 µm, and the total system size was 400 µm. These parameters were

determined from simulations with Dmax ) Dmin ) 2 × 10-9 m2/s and d ) 50 µm, successively decreasing ∆t and ∆z and increasing the size21 until the average diffusion coefficient estimated from the simulated data deviated less than 0.5% from the input value for all combinations of δ and td, and no significant deviation from single-exponential decay could be observed in graphs of the types shown in Figure 7 and Figure 9. An illustration of the model is shown in Figure 4. A total computational time of about 50 h was needed to evaluate the echo attenuation for all parameter combinations and systems. Experimental Materials. 1,4-Bis(2-ethylhexyl) sodium sulfosuccinate (AOT) was obtained from Aldrich. Didodecyldimethylammonium bromide (DDAB) was obtained from TCI. 2H2O was purchased from Dr. Glasser, Basel. Sample Preparation. The DDAB sample (25%) was prepared by mixing the appropriate amounts of surfactant and H2O/ D2O 1:1 (in weight) mixture into a 5 mm disposable o.d. NMR tube. The tube was then flame sealed and centrifuged in alternating directions up to 3 h to ensure homogeneity, and left to equilibrate at 25 °C. AOT samples were prepared using two different procedures. Appropriate amounts of AOT and H2O/ D2O 1:1 mixture were weighed into NMR tubes and then flame sealed in both cases. In the first preparation procedure, method A, the sample was centrifuged in alternating directions for more than 4 h. After mixing, they were left for equilibration at 25 °C. In the second procedure, method B, the sample was centrifuged until it was homogeneous. Then it was placed in a thermostatic bath at 110 °C for 24 h. After that, the sample was allowed to equilibrate while the temperature was gradually decreased to 25 °C during a time period of 3 h. 2H NMR. The internal structure of a lamellar phase depends on the method of preparation and equilibration time.5 On the basis of published phase diagrams of AOT36,37 and DDAB,38 a series of different samples were prepared and examined with crossed polarizers to confirm the existence of a lamellar phase. 2H NMR was subsequently used to get a rough and quick estimate of the domain size and select a series of samples for further studies with diffusion NMR. The 2H spectrum gives information about the interaction between the quadrupole moment of the 2H nucleus and the electric field gradients.39,40 In isotropic solutions, this interaction is completely averaged by molecular reorientations, thus giving a single resonance line. For an anisotropic environment, e.g., for water in a lamellar liquid crystal, the averaging of the quadrupole interaction is not complete, and a quadrupolar splitting is observed. The value of the splitting depends on the orientation of the lamellar director with respect to the main magnetic field, giving rise to the typical Pake doublet for a sample consisting of an ensemble of randomly oriented domains. The Pake doublet is only present if the water molecules reside in a domain with uniform orientation during the 2H NMR time-scale (∼1 ms). If molecular exchange between domains of different orientation occurs on a much shorter time-scale, no quadrupolar splitting will be observed. A smeared Pake doublet is the result of exchange on an intermediate time-scale. 2H experiments were performed at 25 °C using a Bruker DMX-100 spectrometer operating at 15.35 MHz 2H resonance frequency. The spectra were recorded after a 10.5 µs 90°-pulse using a spectral width of 3 kHz. Diffusion NMR. Diffusion NMR experiments were performed at 25 °C using a Bruker DMX-200 spectrometer


Åslund et al.

Figure 5. 2D propagators P(x,,x0,y,y0,t) for the 2D models of a polydomain lamellar phase. The color scale is linear between zero (black) and the maximum value (white). The data in panels a-c/d-f refer to the geometry in Figure 3, with panels a and b having the maximum anisotropy η ) 10, defined in eq 24. The reduced diffusion length l/t , defined in eq 46, is increasing from left to right in the sequence (a,d): l/t ) 0.012, (b,e): l/t ) 0.25, and (c,f): l/t ) 1.0. The arrows indicate the direction of maximum diffusivity and have the same meanings as in Figure 3. The traces to the right of and on top of the surface plots are the projections of the propagator in the x- and y-directions, respectively.

operating at 200.13 MHz 1H resonance frequency with a Bruker DIF-25 gradient probe driven by a BAFPA-40 unit. Diffusion was measured using the PFG STE sequence.41 Four logarithmically spaced values of δ in the range of 2-20 ms were used for each value of td (20, 60, 240, and 640 ms). Only results from experiments with the minimum and maximum values of δ are reported below since the observed trends were always monotonous with increasing δ. For all experiments, the range of logarithmically spaced q values (from 9 × 103 to 2 × 105 m-1) was kept constant by adjusting the gradient strength, which was set to the maximum value of 9.6 T/m for the shortest pulse length. Results and Discussion This section is organized as follows: First, we present results from the simulations starting with the shape of the propagator and the apparent diffusion coefficients for the 2D models shown in Figure 3, and verify our conjecture that the apparent diffusion coefficient decreases with time, even when the spin density and the trace of the diffusion tensor are constant throughout space. We continue with the 1D analog shown in Figure 4, and verify that the evolution of the average propagator is in line with the expected outcome for a true 3D lamellar phase and the simulations of the 2D model. The 1D analog is further investigated by reproducing the analytical expressions for the echo attenuation in the short- and long-time limits, and evaluating the apparent diffusion coefficient as a function of diffusion time and pulse length. Then we proceed to q-space analysis of the data, showing how the pulse length and diffusion time can be adjusted to separate the diffraction features resulting from either diffusion within a single pore or pore-to-pore hopping. From the q-space analysis, the concept of the homogeneous length-scale λhom emerges naturally as the border between the low-q range, being independent of the pulse length, and the high-q range, where the features originating from restricted diffusion are affected by the value of the pulse length. The value of λhom is investigated for the series of simulated

systems and for the influence of experimental parameters. Finally, we determine λhom experimentally for a series of lamellar phases prepared to cover a range of domain sizes. Propagator and D(t) for the 2D Models. In order to compare the results from the different simulated systems, it is convenient to report the time-scale in terms of the reduced diffusion length l/t

l/t ) x2D(t)t/d

(46)

which gives the fraction of the repeating distance that the molecules have traveled during the diffusion time. The results from the simulations to study the shape of the propagator at various time-scales are displayed in Figure 5. At short times, when l/t , 1, the propagator is Gaussian, but with a width reflecting the diffusivity in the various directions. For a polydomain sample with many crystallite orientations, the average propagator is a sum of the individual Gaussian propagators, thus giving rise to multi-Gaussian diffusion as studied by Callaghan and So¨derman.29 At longer times, when l/t is approaching 1, the propagators acquire a shape mirroring the shape of the individual domains. By comparing (b) and (d), one can note that the geometry with sharp grain boundaries gives rise to sharp edges of the propagator. Although the average propagator obtained by averaging all the individual nonGaussian propagators may be rather featureless, the effect of changing the parameters in the diffusion NMR experiment will not follow the simple predictions of the Stejskal-Tanner equation. At even longer times, the propagator again approaches a Gaussian shape, but this time it is isotropic since the materials are constructed to give equal preference for domain orientations in the x- and y-directions. A similar crossover from short-time anisotropic to long-time isotropic diffusion for elliptical particles was recently investigated by Han et al.42 In this study it was concluded that diffusion is Gaussian in both the long- and shorttime limits, while non-Gaussian features are maximized for the intermediate time-scale.



Figure 6. Apparent diffusion coefficients D(t) for the 2D models of a polydomain lamellar phase. The data in panels a and b correspond to the geometry shown in Figure 3a,b. The results are normalized with the isotropic part of the 2D diffusion tensor 〈D〉, and are presented as a function of the reduced diffusion length l/t defined in eq 46. The value of the anisotropy η is indicated in the figure. Note that the scale of the y-axis is different in panels a and b.

The results from the simulations to study the apparent diffusion coefficient as a function of diffusion time and anisotropy are shown in Figure 6. For both simulated systems, there is a decrease of D(t) with both increasing t and η. This decrease is more pronounced for the system with sharp grain boundaries. The control simulations performed with η ) 1 yield D(t)/〈D〉 over the full range of t, thereby proving that the size of the simulation box is large enough for edge effects to be negligible. Although both 〈D〉 and F are constant throughout space, the spatial variation of the director is sufficient to induce the features characteristic of restricted diffusion (viz. D(t) decreasing with t and non-Gaussian propagators) normally observed in porous systems. This fact is in contrast to the common assumption that the long-t diffusion coefficient is the same as the isotropic average of the diffusion tensor used by many authors (including the present ones).6,43 Echo Attenuation for the 1D Analog. In order to compare the results from the different simulated systems, we report the parameters of the diffusion NMR experiment (δ, td, q, and b) using the reduced variables

l/δ ) x2D∞δ/d l/t ) x2D∞td/d q* ) qd b* ) bD∞

(47)

The definition of l/t is slightly different from the one in eq 46. The values of l/t and l/δ represent how far the molecules have traveled with respect to the length of the repeating unit during

Figure 7. Stejskal-Tanner plots with echo attenuation E vs the reduced Stejskal-Tanner variable b* for simulations with (a) l/t ) 0.06 , 1, (b) l/t ) 0.35 ≈ 1, and (c) l/t ) 2.3 . 1. The red and blue solid lines represent the short- and long-time limits given by eq 39 and eq 31 with eq 41. The initial slope for the short-time limit, eq 30 with eq 40, is indicated with the dashed red line. The black circles show the simulated data points for the chosen value of l/t and all values of l/δ. Data points for the same value of l/δ are connected with a solid line to guide the eye. The data points for different l/δ values are overlapping in panels a and c, but not in b.

their respective duration of time. For the present systems with Dmax ) 2 × 10-9 m2/s and Dmin ) 2 × 10-10 m2/s, application of eqs 40 and 41 yields 〈D〉 ) 1.1 × 10-9 m2/s and D∞ ) 6.9 × 10-10 m2/s. As shown in Figure 7, the signal is independent of δ in the limits l/t , 1 and l/t . 1. For the case of l/t , 1, a multiexponential decay with an initial slope given by 〈D〉 is observed, in agreement with eq 29. A single-exponential decay with slope D∞ is observed for the long-time limit l/t . 1. When l/t ≈ 1, both the initial slope and the general shape of the curves depend on the value of δ, indicating non-Gaussian diffusion. These observations agree with the results from the simulations of the 2D models. Hence, it is possible to mimic the diffusion behavior of a polydomain liquid crystal using a simple 1D analog permitting rapid calculations. The apparent diffusion coefficient D(td,δ) was evaluated by fitting a single exponential to the initial decay (1 > E > 0.99) of data of the type shown in Figure 7, and the results for the first series of systems with Dmax ) 2 × 10-9 m2/s and Dmin )


Åslund et al.

Figure 8. Apparent diffusion coefficient D(td,δ) evaluated with eq 21 for the simulated systems with Dmax ) 2 × 10-9 m2/s and Dmin ) 2 × 10-10 m2/s. Data for the three systems with d ) 10, 20, and 50 µm are plotted vs the reduced variables l/t and l/δ, defined in eq 47, in order to compare results from the different systems and cover the whole range of time-scales. The color scale is linear with maximum and minimum values 〈D〉 ) 1.1 × 10-9 m2/s and D∞ ) 6.9 × 10-10 m2/s.

2 × 10-10 m2/s are shown in Figure 8. As expected, l/t , l/δ , 1 yields 〈D〉, while l/t , l/δ . 1 gives D∞. In the intermediate regime Dapp(td,δ) decreases with both increasing td at constant δ and increasing δ at constant td. As noted previously,44 the long-time limit is reached already at l/t ≈ 1. Although this latter regime is uninteresting from the point of view of D(td,δ), useful features appear using a q-space analysis, as will be shown next. Simulated data for the systems with Dmax ) 2 × 10-9 m2/s and Dmin ) 2 × 10-10 m2/s is shown as q-space plots in Figure 9. Since the 1D model is periodic, the data displays pronounced coherence features, which have been observed experimentally for a number of locally ordered systems.45-48 For more disordered structures, such as the lamellar phases studied here, the coherence minima and maxima are smeared leaving a featureless, but nevertheless non-Gaussian, tail at the high-q end of the data.20,49-51 Although we do not expect to observe diffusion diffraction in the experimentally studied systems, let us digress for a moment and comment on this part of the simulation results. In analogy with conventional scattering theory, the total echo attenuation E(q) for a system of interconnected pores (or unit cells) contains contributions from the form factor P(q) and the structure factor S(q):12,52,53

E(q) ) P(q)S(q)

(48)

P(q) describes the contribution from diffusion within a single pore, and S(q) relates to the migration between the units in the interconnected network. Several papers address the dependence of P(q) and S(q) on the system geometry and the experimental variables.16,18,46,49,50,54,55 For the present simulated case, where the size of the unit cell is the same as the distance between the unit cells, P(q) displays minima, and S(q) features maxima at the same positions: q* ) 1, 2, 3, and so forth. The minima of P(q) become deeper as l/t is increasing and move toward higher values of q* as l/δ increases. The maxima of S(q) are visible for the condition l/t g 1, and the amplitude increases with increasing l/δ.53 As seen in Figure 9, optimal conditions for observing the pore-hopping maxima occur at l/δ ≈ 0.6. Interestingly, this observation is general for all the simulated systems with different Dmin. The data in Figure 9 can be divided into two regions: one at high q where E depends on the value of δ, and one at low q where E is independent of δ. With increasing l/t the border qb

Figure 9. Simulated echo attenuation E vs the reduced diffusion distance during the gradient pulse l/δ and the reduced wave vector q* at three chosen values of the reduced diffusion distance during the effective diffusion time l/t : 0.45 (a), 1.0 (b), and 1.6 (c). The color scale is logarithmic between E ) 10-4 (blue) and E ) 1 (red). Values of E below 10-4 are shown in black. Data were obtained for the systems with Dmax ) 2 × 10-9 m2/s, Dmin ) 2 × 10-10 m2/s, and d ) 50 (a), 20 (b), and 10 (c) µm. The arrows indicate the smallest value of q* where E values obtained at different l/δ deviate for the given threshold levels (1.02, 1.05, and 1.15 from left to right).

between the two regions moves toward lower q until a plateau value is reached for approximately l/t > 0.4. The long-td limit of qb defines the homogeneous length-scale λhom. At q values smaller than 1/λhom, the sample appears homogeneous, corresponding to Gaussian diffusion. The exact value of λhom obviously depends on the criteria used to evaluate where the signal starts to depend on δ. Here, qb(td) was defined as the minimum value of q where the ratio between the maximum and minimum E at this particular td was larger than a chosen threshold level. qb for three threshold levels (1.02, 1.05, and 1.15) are plotted vs l/t in Figure 10. λhom is then defined as the value of qb-1 at the plateau reached at long td. On the basis of these observations, we suggest that λhom can be estimated using the following simple protocol for the diffusion NMR experiment: record E(q) for an array of δ and


Figure 10. Border between Gaussian and non-Gaussian behavior qb vs reduced diffusion distance l/t for the threshold levels 1.02 (circles), 1.05 (squares), and 1.15 (triangles). Data were obtained for the simulated system with Dmax ) 2 × 10-9 m2/s, Dmin ) 2 × 10-10 m2/s, and d ) 20 µm.

Figure 11. Homogeneous length-scale λhom for the two series of simulated systems using the threshold levels 1.02 (circles), 1.05 (squares), and 1.15 (triangles) (a) λhom vs the repeat distance d for Dmax ) 2 × 10-9 m2/s and Dmin ) 2 × 10-10 m2/s. The lines are the results of linear regressions. (b) λhom vs Dmin/Dmax for d ) 20 µm.

td values, determine qb for each td, and estimate λhom from the long-td plateau of qb. If the order of magnitude of λhom is already known, it is sufficient to choose a set of parameters consisting of a single value of td (giving l/t ≈ 1), two values of δ (one with l/δ , 1 and one with l/δ ≈ 1), and, say, 10 values of q (covering the range around q* ≈ 1). Thus, λhom can be estimated from as little as 20 data points. Considering technical issues such as the maximum available values of G, and system-related parameters such as relaxation times and diffusion coefficients, we estimate that λhom in the range of 1-100 µm can be determined for typical colloidal systems and biological tissues. In order to test our protocol, λhom was determined for the two series of simulated systems as shown in Figure 11. As expected, there is a linear relation between λhom and d. Regressing the equation λhom ) Ad onto the data yields A ) 4.9, 3.2, and 2.1 for the three threshold levels, 1.02, 1.05, and 1.15, respectively. A decrease of λhom with increasing Dmin/Dmax, corresponding to smaller amplitude of the heterogeneity of the system, is shown in Figure 11b. This decrease is less pronounced for the higher threshold levels. Experiments. Experimental 2H spectra are shown in Figure 12. A well-defined powder pattern is exhibited by the AOT


Figure 12. 2H spectra for lamellar liquid crystals of 2H2O/H2O and AOT method B (top), AOT method A (middle) and DDAB (bottom). All spectra were recorded after a 90°-pulse of 10.5 µs.

sample prepared by method B, indicating the presence of lamellar domains that are large with respect to the length/timescales defined by the 2H NMR spectrum as explained above. Small deviations from the ideal Pake pattern can be observed at the 90° peaks, possibly arising from slight preferential domain alignment at the glass walls of the tube. The Pake pattern is smeared for the sample prepared by method A, pointing to intermediately sized domains. In the case of DDAB, the spectrum displays a predominant smeared Pake doublet combined with a small-amplitude, well-defined doublet. The two different overlapping patterns originate from a majority lamellar phase having small domains, and a minority lamellar phase with larger domains. On the time-scale of weeks, the DDAB sample evolved to produce a single well-defined Pake doublet (data not shown) indicating domain growth. On the basis of the width and shape of the 2H spectra, we can grade the samples according to their domain size: DDAB with small domains, AOT (method A) with intermediate size domains, and AOT (method B) with large domains. A random distribution of the lamellar orientations can be assumed because of the method of preparation and the observed powder patterns (even if some were smeared). More quantitative information about the internal structure of the system can be obtained from the diffusion NMR experiments. The echo attenuations for the H2O in the three samples are shown in Figure 13. Just like in the simulations, the low-q part of the data is independent of δ, indicating Gaussian behavior. For the samples with large and intermediate size domains, as determined with 2H NMR, the signal depends on the value of δ in the high-q region of the data. For the experimental data, qb was determined with the same algorithm used for the simulations. However, an extra step with linear interpolation between the data points was added in order to compensate for the limited sampling of q-space. A threshold level of 1.15 was used as a compromise between the influence of experimental noise and systematic deviations. qb is only weakly dependent on the value of td as shown in Figure 14. We chose to estimate λhom from qb at the longest value of td, yielding λhom ) 39, 15, and

jp076174l

Recommend Documents