Langmuir 1996,11, 785-792
785
Rheological Response of Surfactant Cubic Phases J. L. Jones* and T. C. B. McLeish Department of Physics, University of k e d s , k e d s LS2 9JT, U.K. Received July 19, 1994. In Final Form: December 20, 1994@ We model the response of the surfactant cubic phase to an oscillatory shear strain. The rheology of this ordered structure is determined in terms of the dynamics of the surfactant layer. Above a critical stress, “slipplanes”form giving rise to a bulk relaxation mechanism, which maintains the topology of the network while being able to capture the liquid nature of the material. The system behaves as a Maxwell liquid at low frequencies and displays Voigt cell characteristics at high frequency. We calculate the storage and loss moduli by assuming a pseudolinearresponse and by treating the periodic lattice potential perturbatively. The key results which emerge are the following: there exists a single relaxation time for the bulk, even when there is a distribution of slip plane densities and this relaxation time depends on the value of the average slip plane density. We also find that the periodicity of the lattice introduces other nonlinear effects: oscillations appear in the complex modulus components which depend on both the frequency and amplitude of the applied shear rate.
1. Introduction Within a certain compositional range, a mixture of oil, water, and surfactant forms a periodic, bicontinuous structure, the cubic phase, which exhibits quite intriguing characteristics and has been investigated by many It is well-known that such a structure displays the property of “ringing” when tapped gently, and small angle X-ray diffraction on such materialsg shows welldefined Bragg peaks, from which the size of the unit cell in the crystalline structure can be obtained. A periodic crystal-like structure emerges, in preference to say cylindrical or spherical micelles, because the cubic structure, which is a constant mean curvature surface, optimizes the curvature and is the most favorable energetically. Extensive studiedo have been made of constant mean curvature surfaces and they represent the nonEuclidean solutions to area minimization under the constraint of constant volume fraction. In the present paper we direct our attention to some recent rheological studies of surfactant cubic phases.ll These experiments have highlighted some quite novel properties of the systems. Of primary interest in this work is the observation that the systems exhibit a critical strain, above which both moduli decrease and the elasticity becomes nonlinear. Measurements of the storage modulus G(o)and the loss modulus G ( o )have been made for a range of frequencies, o,and strain amplitudes, yo. For strains greater than the critical strain, ye, both moduli fall off with strain approximately as y0-l. In addition, the frequency specAbstract published in Advance ACS Abstracts, February 15, 1995. (1)Luzzati, V.; Mariani, P.; Gulik-Krzywicki, T. Physics of AmphiphilicLayers;Meunier, J . , Langevin, D., Boccara, N., Eds.; SpringerVerlag: Berlin, 1987. (2)Fontell, K.; Ceglie, A,;Lindman, B.; Ninham, B.Acta Chem. Scand. 1986,A40,247. (3)Anderson, D.M.; Thomas, E. L. Macromolecules 1988,21,3221. (4)Radiman, S.;Toprakcioglu, C.; Faruqi, A. R.; Hjelm, R. P.; de Vallera, A. Colloq. Phys. 1990,C7, 375. (5) Barois, P.; Eidam, D.; Hyde, S. T. Colloq. Phys. 1990,C7, 25. (6)Hyde, S.T.Colloq. Phys. 1990,C7, 209. (7)Gradzielski, M.; Hoffman, H.; Oetter, G. Colloid Polym. Sci. 1990, @
268. 167. _._ I
(8)Bruinsma, R. J . Phys. IZ 1992,2,425. (9)Radiman, S.;Toprakcioglu,C.; Faruqi,A. R. J .Phys. (Paris) 1990, 57 im --, (10)Nitsche, J. C. C. Lectures on Mininmal Surfaces. Volume 1. Cambridge University Press: Cambridge, 1989,and references therein. (11)Radiman, S.;Toprakcioglu,C.;McLeish, T. C. B.Langmuir 1994, IO, 61.
trum for the loss modulus bears the hallmarks of a polymeric liquid; a broad peak and hence a relaxation process occurring with a relaxation time l/umm. The simplest interpretation of this relaxation peak is in terms of a Maxwell model for a liquid. We attempt to account for some of the observed rheological behavior for the cubic phase in terms of a “slipplane” model; i.e., large deformations occur along certain planes in the structure while maintaining relatively small deformations elsewhere. An important feature which must be encompassed by the theory is the transition from essentially a periodic solid, to liquid-like relaxation at a critical value of the stress. This must be included in such a way that the local cubic order is disrupted while the connectivity is preserved. The elasticity which emerges from our model is only truely linear before any slips appear and we have a viscoelastic solid. Above yc, slips exist and the elasticity is “pseudolinear”,in the sense that we define the complex modulus G* in the usual way except that it now contains a strain dependence and may be written as G*(w, y). The slip-plane model captures the essential physics and may be a preferred deformation mechanism to say dislocation movement, because it does not involve the breaking and forming of bonds but simply the flow of liquid in the surfactant layer. 2. A Model for the Viscoelasticity Clearly, a full description of the dynamics of the cubic phase demands an equation giving the position and shape of the surfactant layer at time t as well as the velocity field v(r,t)within the bicontinuous phases. We need to know the function R(x,y;t)which describesthe coordinates of the surface as well as v(r,t).The energy change due to the deformation of the surface is given by the Helfrich Hamiltonian12
The first term is the bending energy, where R1 and Rz are the local principal radii of curvature and K is the bending constagt. The second term is the Gauss-Bonnet term, where K is the effective chemical potential for topological changes and this term is zero if the topology remains unchanged. (12) Helfrich, W. 2.Naturforsh. 1973,28C, 693. Deuling, H. J.; Helfrich, W . J . Phys. (Paris) 1976,37, 1335.
0743-74631951241l-Q785$Q9.QQIQ0 1995 American Chemical Society
Jones and McLeish
786 Langmuir, Vol. 11, No. 3, 1995
-
u
u
L
;l /
I l l
(
l
1
1
I I
l
I (
l 1
l l I
l 1
!~ I
~ 1
!
l
1
Figure 3. Forces acting on layer n of the bulk sample. Figure 1. (a) Unstrained surfactant skeleton. The skeleton is a simple cubic lattice formed from a system of “pipes” of surfactant. (b)A strain y is applied producing a displacement x, - xn-1 between the layers. Sliding in the surfactant layer causes the formation of two junctions one, say (A), with functionality 5 and the other (B) with functionality 3.
3ooc
innr
J U U L
Xn
fast, and the time constant for this is given in ref 11.More importantly though, if this were the only relaxation mechanism, it provides no means of relieving the residual stress from the bulk sheared cubic symmetry. The introduction of slip planes does provide a mechanism for this and allows longer time relaxations. The central idea which we use to model the elastic response ofthe cubic phase is that the relaxation processes are governed by the motion of the surfactant layer. We therefore consider the forces acting on the surfactant and assume that the surfactant skeleton can be described as above by the parameterx,(t). We consider a microdomain of N layers with a surfactant layer which forms a Schwarz-P surface, i.e., a surfactant skeleton which is a simple cubic lattice. We apply a shear strain to a single crystal along the (100) direction. The forces which act on layer n ofthe skeleton are shown in Figure 3. There is a local dissipation term from the viscous forces which act on the layer from the bulk below layer n which is n(%,- %,-Ju. 7 is the viscosity of the bulk material, a the cell size of the lattice, and xn is the displacement of the nth layer from its original position. The corresponding elastic force is a x , - xn-Ja where G is the bulk modulus. For the case where no slip planes exist, the equilibrium condition on the forces gives the following equation of motion
7-nFigure 2. (a) Deformation occurs by sliding along the planes in the (100) direction rather than as shown in Figure 2b. This is because the local curvature of the surfactant layer is altered less by this mechanism. (b)This mode of deformation changes the local curvature of the surfactant skeleton considerably.
Let us now consider the deformation of a single crystal. The picture for the system simplifies a little ifwe maintain that the connectivity of the phase does not alter during the deformation process and that strain is accommodated by layers in the lattice structure sliding over each other. The sliding takes place in the surfactant wall. Our coordinate describing the position of the surfactant layer in this case only needs to contain the position of the “pipe ends” as shown in Figure 1. Furthermore, if we assume that all the lattice cells in one layer remain the same size and all slide the same distance relative to their adjacent layers, then we only require a single variable to describe the surfactant, x,(t) which is the displacement of the nth layer from its unstrained position at time t. There is another relaxation mode which we will mention briefly but which is dealt with in detail in a previous paper.ll Figure 2a shows that the surfactant tubes are all either parallel or perpendicular to the strain. The perpendicular ones have undergone a relaxation process from being initially strained as shown in Figure 2b. The relaxation takes place by flow within a tube, it is relatively
-
and for small local deformations we can take the continuum limit as n s and this gives
(3) The above equation is for the bulk material and we can define a relaxation time TI= v/G. The steady-state solution of this equation has the following form
x(n,t) = uP(t)n+ D(t>
(4)
where D(t> is zero if we impose the fured boundary conditions x(0,t) = 0 and x(N,t) = Nay, but will become a necessary part of the solution when slip planes are introduced. In that case, eq 4 is true in a piecewise fashion for elements of crystal between slip-planes (see below), and D ( t )gives information on the distance slipped. The quantity P(t) is defined as P(t) = a-%(n,t)/&z and is the gradient of the displacement in the bulk sample. For the no slip case this is equal to the strain y(t).
Rheological Response of Surfactant Cubic Phases a) Strain is uniform across all layers
i=n
b) Slip has occurred and the strain in one layer is large but in all others is very small.
slip-p
X8
.....
Na
Langmuir, Vol. 11, No. 3, 1995 787
..,, .., . 1
i=l ,
I
l
(
l
I
,
l
l
I
l
I
Figure 4. Slip plane formation. Shown here is the skeleton when a fixed strain y is applied to the system. xn is the displacement of the nth layer from its unstrained position. In (a) all the layers have the same strain y , but in (b) a slip has occurred most of the strain is accommodated by the slip plane. The strain in all the other layers is very small. This can lead to a lowering of the energy provided there is sufficient energy t o create a slip.
Let us now consider applying a fixed strain yo to the sample. When no slipping occurs, P(t) = constant = yo. We can write a general expression for the stress exerted at the top surface of the sample and this is
The first term gives the elastic stress and the second is the viscous component. For our simple case, this gives ~7 = Gyo and here G is a constant. Correspondingly, the stress is a constant with time and no stress relaxation occurs. The system behaves as a solid at long times and is effectively a Voigt cell. 3. Viscoelasticity with Slip Planes It is clear that the model outlined in the previous section does not permit bulk relaxation and essentially describes a solid even at long times. However, if the structure forms slip planes we are able to maintain the same crystalline structure but allow relaxation processes. Energetically there are good reasons why this could occur. Figure 4 shows the formation of a slip plane. In order for a slip plane to arise the energy must be lower, when most of the strain is taken up by a large value of xn-xn-l a t the slip plane than when the strain is equally distributed among all the planes. These are the energetic considerations of slip-plane formation. Moreover, it is highly likely that there exist weaknesses in the crystal lattice which yield to the formation of slip planes at relatively low stresses. The density of the slip planes in the crystal is, under most conditions, probably kinetically determined and a function of the irregularities in the lattice made at formation. However, if the activation energy to create slip planes is not large, then an argument for the equilibrium density after a bulk shear may be made as follows. The elastic energy u is periodic in the relative displacement of the layers i.e. u = -KOcos Ax27da where KO is an elastic force constant. Consider applying a relatively large strain yo 1and that all the force constants KO are the same (i.e., the structure is a uniform perfect crystal). We can calculate the relative number of slip planes nslipto unslipped nunslippedusing an argument which minimizes the total energy at coexistence of the slipped and unslipped planes. The total energy of the system is ET = Eslipnslip + Eunslippednunslipped. When this energy is minimized it is found that the slip plane density N-' =
-
4.
We can then investigate the case where weaknesses exist in the structure, i.e., there exists a Gz at some layer
such that Gz *E G. In this case slip can occur a t a much lower strain. This is explained in the followingway: the stress across the sample is continuous so u = Gybulk = Gayslip. Therefore yslipis much larger than ybulk, the strain in a nonslipping plane. In fact if we assume that slip occurs when the local strain is unity, then the critical strain for slip is ycrit = G2/Gu1k.Both thermodynamic and kinetic arguments therefore support the conclusion that under finite shear, the soft crystal will minimize the stored elastic energy by forming slip planes of locally very high deformation. We examine their dynamical consequences in the following sections.
4. One Slip Plane Linear Response Once a slip has occurred, the dynamics are dramatically changed. Because deformation in the slip planes is so large, the response in the surfactant layer is due to viscous dissipation. The effective viscosity will have a large contribution from the surfactant flow along the surface, so we invoke a slip-plane viscosity V Z which differs from the bulk viscosity 8. The dissipation in the slip plane is then ~ Z per W cycle. This allows us to treat the dynamics in the slip plane as primarily viscous interactions with the lattice potential added as a perturbation. We look at the case where the periodic potential is added in a later section. We now consider the dynamics when the periodic potential is zero in the slip plane. The equation for the bulk is the same as for the nonslip case, but there is a new boundary condition which occurs at the slip plane. Balancing the forces in this plane requires GP
82 . + 8P - -(x,, a
- k,,-l) = 0
(6)
wherexn -x,-1 is the distance slippedD(t1. (We note that D(t)can represent a slip over many lattice spacings, which will become important in section 7 below). Applying a fured strain to the system gives the boundary conditions xo = 0 and X N = Nya and the distance slipped is then D(t) = Nya - ( N - 1)Pa. Equation 6 may be solved by substituting for x,, - x n - l and we find, independent of the position of the slip plane, an exponential relaxation with a time constant t = (17 Nr2)/G (forN >> 1). We therefore have
+
P(t) = yo exp(-t/t)
(7)
and the stress is given by
a=GP+qP
(8)
and so the relaxation modulus is
G(t)= (G
+ r/t)exp(-tlt)
(9)
In this case, stress relaxation can occur and the behavior has some of the characteristics of a Maxwell cell. It can be seen that t depends on the domain size N , which in this example is fixed. This has been discussed in the context of general viscoelastic response from fluids with domains by Doi and Ohta.13 At this stage, we note an important point that it is not applicable to calculate G* from the transformation of a t ) , i.e., from
G*(w)= iw
G(t)eiWtdt
(10)
There are two important reasons why this is the case. Firstly, the stress is actually a superposition ofresponses, and although P(t) is exponential, there are some Voigt (13) Doi, M.; Ohta, T. J. Chem. Phys. 1991, 95 (2).
Jones and McLeish
788 Langmuir, Vol. 11, No. 3, 1995 elements in the system which give solid-like behavior in the high-frequency range. Moreover, when we add back the periodic potential of the lattice, the nonlinearities here give rise to a nonsinusoidal stress. We must therefore calculate the storage and loss moduli directly from the stress-strain relationship when a periodic potential is applied to the system. Consider applying a strain y = yo sin wt. Then for a linear response, the resulting stress will be given by (T
= y o ( G ( o ) sin ot
+ G ( o )cos wt)
(11)
where G(o)is the storage modulus and represents the stress contribution which is in phase with the strain. G ( w ) is the loss modulus and gives the contribution from the stress which is n12 out of phase with the strain. To proceed with the calculation of the modulus components, we write the total displacement at the top layer XN in terms of D(t) and P(t) which for large N is x N = Nay, sin ot = NP(t)a
+D(t)
Storage and Loss Moduli versus Frequency G=IO .ya=O.O1. r,=I.O
Storage and L O ~Moduli S versus Frequency. G=10 ,yo=O.Ol,r,=O.I
I- -
6.0
(12)
G',(w) G",(U)
which allows b(t)to be expressed in terms of P(t). As for the fixed strain case this is then substituted into the boundary condition (6)to give a time-dependent equation for P(t). We have the following simple equation for the dynamics for the case where the lattice potential is zero
G +Po= -Po r12Nayow cos ot B B
(13)
+
where B = a(q Nqz)and t =BIGa. Solving this equation leads to a solution for P&) of the form
! I -2 0
0.0
/
'I
/ /
00
40
20
1IT
log,,0
60
l/(TTd In
Figure 5. The zeroth order storageand loss moduli for different values of the slip plane density and viscosity ratio (a, top) r,, = 1.0, (b, bottom) r,, = 0.1.
+
where A = q2N/(q Nq2). The total stress can be found from (5) and in this zeroth order case the moduli can be extracted from eq 11 and are given by G,'(w) = AG
, where
tl
(022- 0 2 2 , 2 ) (1 + w 2 2 )
L I U
(1
(15)
+ 022)
= q/G as in section 2.
5. Discussion of Zeroth Order Results We analyze the results in terms of two important parameters, one of which is N-l, the slip plane density and the other is the viscosityratio rv = qIq2. The important timescales in the system are defined in Figure 5: t is the bulk relaxation time and is determined from the peak of the relaxation spectrum; tl is the local relaxation time; tmin = (ttl)l'z which is obtained from the minimum in the spectra. In this notation the ratio of the timescales for large N is tltl = N/rv. The plots of G ( w )and G ( w ) are shown in Figure 5 and these show how the position ofthe peak llt ofthe spectrum 115 varies with the slip plane density N-l. Around the peak, the behavior is that of a Maxwell fluid. At the lowfrequency end of the spectrum there appears to be a terminal zone and €or t >> tl then G w2t2and G = 072. At high frequency the dynamics are all dominated by the
bulk viscosity q since slipping does not occur. The behavior is therefore solid-like at high frequency, with a small viscous dissipation coming purely from the bulk viscosity. We can consider two regions of the frequency spectrum: w 2 t t 1 > 1, for which the phase angle 6 is given by tan 6 = G"/G = l / w ( t - tl)and tan 6 = w t t l / ( t - tl), respectively. For the case when t >> tl (i.e., when N 1and r? 1) we get a phase angle given by tan 6 llot. When w t > 1, the phase difference is tan 6 w t l which for high frequencies gives 6 n12. In this case the system has a Voigt cell response. If we examine the case where t t1 (which occurs when rv and t are of the same order of magnitude), then we find the high- and low-frequency phase lags are both still nJ2, but in this case the value of the prefactor A is less than unity. In the case where the two relaxation times are widely separated A 1. In summary, varying the slip plane density N-l and the viscosity ratio rv slides the spectrum along the frequency axis. A large ratio t I z 1 stretches the curves along the frequency axis as well as changing the height of the peaks by a factor of up to l/2. At both the high- and low-frequency limits the value of 6 is n/2 for all cases.
-
-
-
-
-
-
-
6. Polydispersity in the Slip Plane Density
All of the above assumes that the slip plane density is a constant. Let us now consider the effect of introducing
Langmuir, Vol. 11, No. 3, 1995 789
Rheological Response of Surfactant Cubic Phases polydispersity into our system. Let n be the number of slip planes and Ni the number of layers between two slip planes with the labels i and i 1. N is the total number of layers in the sample. The ith slip plane carries a slip distance Di and the ith domain has a deformation yi and modulus gi. The stress exerted by the top of the ith layer is ai and because the stress is a constant throughout the sample
+
0. c = g .c yc.
+ q.p. c c = g . . + r.y. jYj
j
j
vij
vary from one domain to another one. In this case we have
vi
VzJ’i
=-
a
(27)
Using the fact that u is constant throughout the sample, we may write
(17)
ua D.=-
The total displacement at the top layer for N >> n is
(28)
72i and then using the constancy of y , eq 20 becomes n
i
(29)
and so differentiating
(19) The stress in the slip planes themselves is uniquely equal to v&/a for all i, and so for ~2 = constant, this condition determines a universal slip rate among the slip planes. This allows the first sum in eq 19 to be simplified and so we obtain
By writing the stress using the forms of both (21) and (291, we have the following equation of motion for y .
(30) Averaging the viscosities over all planes gives a time constant of
(20)
For the present let us assume that g, and vi are the same for all layers and label them G and 7 , respectively. Then if the stress exists, it is given by
~ i ( t=>~ p i ( t + ) GyJt)
(21)
and the uniqueness gives
pi(t) = S, exp
-G(t-t’)/?j
a(‘’)dtt
(22)
In other words all the y(t)’s are equal Vt. In a slip plane we can write the stress in terms of a purely viscous contribution (23)
However, we also require that this is equal to the stress in the rest of the sample and so we may write using (201, (211, and (23)
This may then be rewritten to obtain
which when solved for yi gives a relaxation time t=
V
+ (Ni)V2 G
Thus we arrive at the initially surprising result that we get a single time constant when there is polydispersity in the system rather than a spectrum of different ones. We can now consider the effect of allowing the values of q-2 to
(31) where N/n(l/~z)= (Ni)/(1/~2).Essentially it can be seen that there is only ever one time constant for the system, which is found by averaging the important parameters. We defer experimental comparison of these results to the discussion. 7. Nonlinear Theory The linear model which we have described so far misses features which are due to the connectivity of the structure. Therefore we introduce a weak elastic force at the slip plane which reflects the translational symmetry of the background structure. The periodic potential perturbs the purely viscous behavior. The response is now nonlinear and when a sinusoidal strain is applied, the stress contains higher harmonics than just the first and we must redefine G*. The choice of G* is determined by the experimental method used to measure it. This involves the measurement of the fundamental frequency and the amplitude. G* is then assumed to have a sinusoidal form with an appropriate amplitude and frequency. In order to extract the storage and loss moduli in the nonlinear case, we are going to use a “pseudolinear))theory defined as follows. We put all the nonlinear terms into the calculation of the stress for a sinusoidal strain and then only consider the lowest harmonic terms in the stress expression. This enables us to extract the two moduli G ( w ) and G ( w ) from an equation of the form
G(t)= G ( w ) y , sin ot
+ G ( w ) y , cos w t + higher harmonics (32)
This approach will become clearer in later sections when we actually calculate the moduli. Let us examine the situation when the nonlinear lattice potential is small and treat it as a perturbation. We begin by describing, in general terms, the equations which are appropriate when a strain y = yo sin wt is applied to the
Jones and McLeish
790 Langmuir, Vol. 11, No. 3, 1995
system. The new boundary condition a t the slip plane when the lattice potential is included is
GP + y P - :(An - An-,)
Storage and Loss Moduli versus Frequency. E=O.Ol, G=1O5,y,=O.OI, r,=0.1, N=1000
60r--
I----
’
n a
G’,(w)
I
- G,sin 2 4 x n - xn-,) = 0
A-
(33)
where GZis the elastic modulus at the slip plane due to the lattice potential. The displacement at the top layer is given by (12). Expanding in the variable E , which is our perturbation parameter, we write the variable P(t) as
P(t)= POW + €Pl(t)+ E2P,(t)+ ...
(34)
and Gz = EG. Examination of the terms which are first order in E in the perturbation expansion produces the following equation from the boundary condition
+ BP, - G sin(-)2nD =0 a
GP,
(35)
’
0.0 00
I I
I 20
Storage and Logs Moduli versus Frequency. ~=0.001,G=IO ,yo=O.Ol,r,=O.l, N=IOOO 60
---- G’,b)
The solution of this equation can be expressed in integral form using the Green function Pl(t) and admitting only terms up to Po@) in the source term. The expression for the first-order correction to the value of P(t)is therefore
--------
- G”,(w) G’,(u) G”,(o)
-
J, exp(t’/t) sin 2 n ( - ~ ~ ~ (+t ’ )
~ , ( t=) exp(-t/z)$
60
40
l%,W
I
Nyo sin ut’) dt’ (36) Ifthe form ofPo(t)from eq 14 is substituted, then this may be rewritten as
P,(t) = exp(-t/t)-
sb
1 t
exp(t’/z) sin(a sin(#
+ wt’) dt’ (37)
+
+
1 wz
yoexp(s’/wz) sin(a sin s’) ds’ (38)
+
where tan q.5 = (1 02zzl)/w(t - 71) and a = 2jtNy((l (r,/(rv N ) ) 2 w 2 t z ) / ( 1 Pl(t) can be evaluated relatively simply by expressing the sine in the integrand as a sum of Bessel functions as follows
+
+
m
J,,,(a) sin((2Z + 1)s) (39)
sin(a sin s ) = 2 l=o
The value of a determines the first nonlinear correction to the response. The full correction term is given by
wz(21
+ 1)cos((2Z + 1)(++ 0z))l (40)
We choose only the first harmonic terms for the definitions of G and G“ and so dropping the higher harmonics the relevant part of the solution becomes
P, =
2J,(a) [(cos # (1 w 2 P )
+
I
2.0
4.0
1 6.0
‘O&P
and then making the change of variables s = q.5 wt’ we have the following expression for the first-order correction term
P,(t) = exp[-(s/wz)l-
I
0.0 0.0
+ ut sin $1 sin wz + (sin + - cut cos #) cos or1 (41)
Using the above defining equation, we are then able to calculate the two components of the complex modulus,
Figure6. First-order correctionsto the storage and loss moduli. GT’and GT”are the total moduli with the first-order correction added. Go’and G{ refer to the zeroth order moduli. (a, top) E = 0.01; (b, bottom) E = 0.001. taking into account the first-order corrections from the harmonic potential of the lattice. The results are shown in Figure 6 . It is clear that the asymptotic value of a determines the maximum value of the correction. For consistency within our perturbation theory it is important that this term is small compared to the zeroth order one. Taking the asymptotic limit of a, then a 2zNy as w 0 leads to the following criteria for the perturbation expansion to be valid: E
