A variational approach to non-linear dielectrics: Application to polyelectrolytes

Journal of Electrostatics, 20 (1987) 219-232

219

Elsevier Science Publishers B.V., Amsterdam - - Printed in The Netherlands

A V A R I A T I O N A L A P P R O A C H TO N O N - L I N E A R D I E L E C T R I C S : A P P L I C A T I O N TO P O L Y E L E C T R O L Y T E S

ANGELO MORRO

Biophysical and Electronic Engineering Department, Viale Causa 13, 16145 Genova (Italy) and MAURO PARODI

Institute of Electrical Engineering, Piazza d'Armi, 09100 Cagliari (Italy) (Received November 10, 1986; accepted in revised form May 3, 1987)

Summary The paper deals with the determination of the electric field and the ionic distribution in a dielectric with an electric-field dependent dielectric permittivity. First it is shown that the associated nonlinear problem can be cast in a variational formulation. Then the functional is specialized to the case of rod-like polyelectrolytes in water solutions and the effects of the permittivity on the electric field and the ionic distribution are investigated. The numerical results are compared with the corresponding ones for a constant permittivity. Besides other features, the comparison shows that the nonlinear effects on the electric field are usually remarkable near the macroion surface.

Introduction The dependence of the dielectric permittivity on the field intensity is at the basis of outstanding electric phenomena. Its theoretical characterization, however, is far from being established. Meanwhile the experimental data are difficult to obtain for a number of substances (e.g. gaseous dielectrics) for which the variation of dielectric permittivity e is negligibly small, compared with itself, also for extremely large values of the field, say 108 V/m. In liquid dielectrics measurable variations of e can occur for field intensities of 107 V/m or less [ 1 ]. In these situations the non-linear behaviour of e has to be taken into account whenever the values of the field are expected to be so high in a significant part of the domain under consideration. A typical example of this situation is provided by the electric field near a polyelectrolyte macroion in a water solution. For instance, the relative dielectric permittivity of water near a DNA macroion is estimated to be a few units (say 4-8) while its bulk value is 80 [ 2 ]. This variation is likely to play a decisive role in the behaviour of the ionic

0304-3886/87/$03.50

© 1987 Elsevier Science Publishers B.V.

220

macroio-~__~-n -'~ :!1 ]:

-

R Fig. 1. A cylindrical cell with a negatively charged maeroion.

atmosphere surrounding the macroion. Of course, from a mathematical viewpoint the dependence of e on the electric field results in additional difficulties in finding the solution to the pertinent system of differential equations. The polyelectrolyte behaviour is becoming more and more interesting in several contexts such as classical electrochemistry, the modelling of electric phenomena in biology [ 3 ] and molecular electronics [ 4 ]. Motivated by this extent of applicability, this paper is devoted to polyelectrolyte solutions viewed as non-linear dielectrics. The aim of the paper is twofold. First, to show that, for any dependence of e on the electric field, the dielectric (and, specifically, the polyelectrolyte solution) can be described by a genuine variational formulation. Second, on exploiting the variational formulation, to obtain an accurate solution and to examine in detail the behaviour of the electric field and the mobile ions. To do this we adhere to the current literature [ 5 ] and regard the polyelectrolyte solution as constituted by cylindrical elementary cells (Fig. 1 ). In any cell the mobile ions (counter-ions) are allowed to obey a Boltzmann distribution law. This scheme results in a system of non-linear equations in the electric potential and the ionic concentration. After having set up the variational formulation the authors prove that the minimality of the functional is ensured by the positivity of the differential dielectric permittivity. T h e n the electric potential in the cell is determined by means of a proper application of Euler's method.

221 The model for a polyelectrolyte solution

The polyelectrolyte solution is viewed as a lattice of cylindrical elementary cells [ 5 ]. Each macroion is modelled as a cylindrical rod whose axis lies along that of the cell. On the basis of the arguments given by Fuoss et al. and the related experimental data [ 6 ], the rod-like macroion is regarded as isolated. The charged groups of the macroion are taken to be distributed uniformly. The solvent carries the mobile ions and occupies the cylindrical corona a ~
Although our purpose here is to set up a variational formulation for a cylindrical elementary cell of a polyelectrolyte solution, we develop an approach which applies equally well to any charge distribution. On adopting the Boltzmann distribution for the mobile ions, Poisson's law for the solution can be written as

( z0)

V" (eJzq~) +qzco exp --~-~ = 0

(1)

Here ~ is the unknown electric potential while q is the proton charge, z the valence of the mobile ions, k the Boltzmann constant, and T the absolute temperature. The non-linearity of the solution is described by the dependence of the dielectric permittivity e of the solvent on the electric field in the form

For the sake of convenience the potential ~ is taken to vanish at the macroion

222 surface. Hence the parameter Cotakes the physical meaning of volume concentration of the mobile ions near the macroion surface. The volume concentration

C=Coexp

( z0) -~-~

qN

and then 0 is subject to a constraint. Precisely, since is the overall charge of mobile ions in V, the electric potential must satisfy the condition

'z'c° f exp( -qzO~

(2)

V

The first step toward a variational formulation for the problem under consideration consists in the determination of a functional [ ( 0 ) , if it exists, whose Euler-Lagrange equation coincides with eqn. (1) to within a possible inessential common factor [ 7 ]. This question is answered through the general theory pertaining to the inverse problem of the calculus of variations [ 8 ]. Equation (1) is a second-order non-linear partial differential equation in the unknown 0. Setting aside unnecessary mathematical details, we recall that a second-order differential equation

L(u) =0

(3)

in the unknown function u, admits a variational formulation relative to the L2 (V) inner product if and only if the potentialness condition [9 ]

(~U,j

,h

is satisfied; here a comma followed by a letter j, h means differentiation with respect to the (Cartesian) coordinates xj, xh with j,h = 1,2,3. Moreover the theory shows that if eqn. (4) is satisfied then eqn. (3) is the Euler-Lagrange equation of the functional 1

V

0

w h e r e / * depends on an arbitrarily chosen reference value u* of u and, possibly, on the parameter Co. In the present case u = 0 and L is the left-hand side of eqn. (1). To show that eqn. (4) holds in connection with eqn. (1) is quite straightforward, though highly non-trivial because of the dependence of e on Once this is proved we apply eqn. ( 5 ) for the determination of the functional [ ( 0 ) . Letting V0* = 0 it follows that

IV04.

223 Lv¢l2

qz¢ [(¢)=f*+I [ ½ I e(~)d~+kTcoexp(--~--~)]dz, v

(6)

o

Accordingly, apart from the additional term f*, the functional (6) is the volume integral of the Lagrangian density I!g¢l 2

Z(¢,VO)=½ ~ e(~)d~+kTcoexpl---~/ o

Account for constraints, such as eqn. (2), is usually accomplished through the multiplier rule. In the present case this rule results in the addition of the quantity

v

to the functional (6). However the Lagrange multiplier/1 is to be viewed as an additional unknown for the variational problem and this makes the procedure scarcely applicable. A more profitable way is based on the fact that f(@) has to be extremum with respect to the parameter Co. Specifically

Of Of* exp(_qz¢~ = c)--Co- O-~Co-t-kT J \ k T ; d~, Using eqn. ( 2 ) we can write

Of*

NkT

OCo- fzlco Hence

f* = NkT ln co+ Jzl

Physically, a proper expression for the constant $ allows f to be viewed as the free energy of the cell [ 9]. The value of ~, though, is completely irrelevant to the description of the cell. That is why, for formal simplicity, we let ~ vanish. In conclusion

[(~) = ! [ ½,v!,~~(~) d~+kTco exp ( - ~T)I d~'-NkTlnc° Izl

(7)

Henceforth we consider the specific case of the cylindrical cell of the polyelectrolyte solution. T h e n V is again the region

224

V= { (r,O,x3): a<~r<~R, 0~<0 ~<2n, 0~x~ ~
VO(R) = 0

(8)

On taking advantage of the arbitrariness of an additive constant in the definition of 0, for the sake of mathematical convenience we assume that ~(a) = 0

(9)

Moreover, because ¢ is independent of x3, letting n be the unit outward normal to OV we have

Vo'n=O

atx3=0,l

(10)

Then we may state the variational problem as follows: Among the functions

belonging to C2( V ) and satisfying the boundary conditions (8-10) find the field O which makes the functional (7) stationary. Of course, from an operative viewpoint, the search for the solution 0 is greatly simplified if f is not only stationary but also extremum at the solution. The next section clarifies how the functional can be extremum at the solution. Extremality of the functional and differential permittivity We examine now whether, given a value of the concentration Co, the functional ( 7 ) has an extremum at the solution of the variational problem stated above. Let h (r) • C2 [ a,R ] with

h(a) =0 Then, letting 0 stand for the solution to the variational problem, we represent the class of admissible functions as 0 + ah, ~ • [0,~] for some positive &. The functional f has a minimum at ¢ if df(~+ah)

i, =o=O

(11)

and d2f(O+ah) da

>0,

a•[0,&]

(12)

for any admissible h. Upon evaluating df/da and applying the divergence theorem we can write

-d--d

-- -

v. (evo) +qzco exp

kT]] ov

225 The volume integral vanishes because ~ is the solution to the Euler-Lagrange equation (1). Meanwhile the surface integral vanishes because of the boundary conditions. T h e n eqn. (11) holds. Let e' denote the derivative of e with respect to the argument, i.e. e ' = d e / d e 2. The evaluation of the second derivative of f with respect to a yields

d2f(~+ah) da2 - f { 2e'(IV~)+aVhl2)[(V~+a['h).Vh] 2 V

+e(]V¢+aVh]2)]Vh]e +c°

kT

exp

dp

It is apparent that eqn. (12) holds if and only if 2 e ' ( ] V ~ + a V h ] 2 ) [ ( V ~ + a V h ) . V h ] 2 + e ( ] V ¢ + o t V h ] 2) ]Vh] 2 >10

(13)

for any admissible function h. This leads us to investigate the validity of eqn. (13) for non-linear dielectrics. In this regard, letting

E=-(V~+aVh),

v=Vh

we write eqn. (13) as 2e'(E 2) ( E - v ) 2 + e (E2) v 2 >/0

(14)

Quite remarkably, the inequality (14) expresses the positiveness of the differential permittivity. To see that it is so observe that from

D=e(E2)E we have

OD 2 -O-E=e(E )•+2 e'(E 2) E ® E where I is the identity tensor and ® denotes the dyadic product. Then, letting v be any vector we have

D ) = e ( E 2) v 2 + 2 e' ( E 2 ) ( E - v ) 2 v- ( O ~-~v This provides the sought result and then allows us to view inequality (14) as a generally valid property for the permittivity. Observe now that if e'>/0 then inequality (14) is obviously satisfied. Usually, though, e' < 0 and then not every function e (E 2) guarantees the minimum of/. Indeed, since

(E.v)2 <~E2v2

226

only those functions ferential inequality

e(E 2) guarantee the minimum of f which satisfy the dif-

2~'E2+~>~0

(15)

For the sake of definiteness we show that there is a family of functions satisfying (15). Because m <

2E 2

E~o+ 2E 2

for any non-vanishing Eo, the solution e to E02+ 2 E 2 satisfies (15). In terms of ¢(0) and ¢(oc) such a solution is

e = e ( ~ ) Jr [ l + 2E2/E211/2

(16)

The permittivity functions (!6), parametrized by e ( 0 ), e ( ~ ), and E 2, satisfy (15) and are likely to constitute a physically sound model of non-linear dielectrics. In conclusion, whenever the permittivity satisfies (15), as for example in the case of eqn. (16), the functional (7) with the conditions (8)-(10) has a minimum at the solution. Numerical solution by Euler's method

For the sake of generality the implementation of the variational approach is performed by employing non-dimensional quantities. They are defined as follows q~

-kT Co-

2na2 co 2

r

R

~=a

F =-a

~'=aV

q22 fly - 2 n e v k T

evbeing the vacuum dielectric permittivity. Hence the electric field can be written as qa

being regarded as a function of ~. Accordingly the Poisson-Boltzmann equation (1) becomes

227

Letting 0 = k T/qalEol, we can write the permittivity function (16) as

~(0) -E(oo) [1 + 202 (d~/d~) 2] 1/2 The boundary conditions for {# are d~

e#(1):0,

~--~-(1)=ft.,

d~

~-(F) :0

The quantity fla represents the normalized electric field intensity at the macroion surface and is determined, via Gauss' theorem, by the equation

i{

~(o)-E(oo) ],#~: #v -2 .q;72j

In terms of the potential ~ and the parameter Cothe functional (7), to within inessential constants, can be given the form r

1

+flvgo exp(-zq~)

(d~o'~2 + 1

v

[1

]]

(d(~2 1/2

1 ~d~-~ln~o

Unlike Ref. [5], the adoption of the Ritz method [7,10] for the determination of the solution is quite unpractical here mainly because the needed set of trial functions turns out to be much too large. The main reason for this is that, as is shown in the next section, the resulting electric field versus ~ is a markedly decreasing function. Roughly speaking, it can drop by more than one order of magnitude within the first few percent of the whole interval [ 1,F]. A proper formatulation of Euler's finite difference method [11] leads to accurate results while taking a much simpler work for the computer implementation. Since the ideas underlying the general formulation of Euler's method are well-known, emphasis is given here only to the parts which are specific to our problem. The normalized electric field, namely d~/d~, is represented through a polygonal curve. Since the field is (expected to be) a rapidly decreasing function near the rod ( ~ -- 1 ), it is convenient to employ a denser set of segments within the first part of the domain, [ 1,U] say. Let m be the number of segments in the interval [ 1, U] and n the number of segments in [ U,I"]. The abscissae of the vertices of the polygonal curve are

228

( 1+ k/ll, l U+(k-m)A2,

k=O,l,..,m k=m+l,m+2,..,m+n

with/11 = ( U - 1 )/m, A2= ( F - U ) / n . Letting Yhdenote the corresponding values of the field we have Yo =fla,

Ym+~--0

Correspondingly the values of the potential are approximated via the trapezoidal rule as (

A1

1

(tim+/12

| ~k

+ Y_____~.fk ~ Y i ,

2

i=1

Ym-bYk F ~ Y i 2 i=rn+l

k=O,1,..,m, ,

Such approximation makes the functional [ into a function F of the unknown ordinates Yl,Y2,..,Ym+ n- 1. These ordinates are determined by requiring that the function F be extremum, namely OF

Oy,

=0,

k=l,2,..,m+n-1

(17)

while OF

0~o

=o

(18)

The computer algorithm, giving the solution of the system (17)- (18) is based on the Newton-Raphson method [ 12 ]. Since the matrix of the partial derivatives of the system of equations is symmetrical, a considerable saving of memory space is possible because only the upper triangular part of the matrix has to be evaluated and stored during the iteration procedure. The whole algorithm has been implemented in FORTRAN language and double precision on a HP 1000 computer by using the VIS routines for the vector processing. The resulting computational efficiency is rather good: for instance, when m + n = 90 each iteration step takes less than 15 s. The numbers m, n, and U are adjusted depending on the particular situation. Some indications about their values are given in the next section. As a general criterion, a good choice can be validated by testing the resulting solution properly. In view of Gauss' theorem the exact solution must fulfil the equation ck ~oevflvZ j e x p ( - - z ~ ) ¢ d~--e(fla)fla--~(Yh)Y~h 1

(19)

229 ~ke [ 1,F]. T h e n the accuracy of an approximate solution may be estimated through the difference between the two sides of eqn. (19) for any k. Results

As an application, detailed results have been derived in the outstanding case of DNA macromolecules. According to the literature the value of the inner radius is taken to be 10 A and the spacing between two electronic charges of the chain 1.7 P,, which implies that ~ = 0.588 ~, -1. The outer radius R depends on the molar concentration M of the DNA solution. By taking M in m o l / d m 3 and R in A we have ( 1027~ ~1/2

R=\~]

where NA is Avogadro's number. Owing to the high value of the electric field the water dipoles close to the DNA macromolecule are markedly structured. This makes the dielectric permittivity e to be approximately equal to 4 ev, which value is usually assumed for e (or) [ 2 ]. Following these indications, a value of 6 ev for e has been assumed in correspondence of the rod surface. Then, because the permittivity of the bulk water is equal to 80 ~, we get 0 -- 0.0247. The mobile ions have been regarded as monovalent and the temperature equal to 300 K. The typical features of the main quantities of the cell near the rod of a DNA macromolecule are shown in Fig. 2. The value of R is 72 ,~ (i.e. F = 7.2 ), which corresponds to M=30 m M / d m 3. Because of the large variation of the quantities considered, the behaviour versus ~ is represented in the region around the rod, 1 < ~ ~<2. In order to emphasize the role of the non-linearity a comparison is made with the case ~ -- 80 ev (dashed line ). The normalized charge density ~ drops by three orders of magnitude in the interval 1 ~<~ ~<2 (which corresponds to a distance of one molecular radius). The relative permittivity e/ev is also plotted; it indicates that the non-linear behaviour is actually confined to 1 ~<~ ~<2. As ~ > 2 the permittivity is approximately constant and equal to the bulk value. This may be of interest in the modelling of polyelectrolytes. Generally speaking, the behaviour of the field and the charge density in Fig. 2 stress the value of the algorithm adopted. An accurate solution requires a large number of sampling points in the region 1 ~<~ ~
230

~-/~

102

d~O

10

4

-71 I0

2

.01

i

I

i

i

i

1,5 ~;= r/a

I

i

i

t

2

1 ' ''" .05 .1

'

'

I'

~ '~" 1

), [~,-']

Fig. 2. Plot of the electric field, the charge density, and the relative dielectric permittivity near a D N A macroion (continuous curves). The dashed line represents the electric field in the case of the constant permittivity ~ = 80 ev. The outer radius R of the cylindrical cell corresponds to a D N A concentration M = 30 raM/din a. Fig. 3. Plot of the normalized average distance of the ionic cloud versus the charge parameter (continuous line). The dashed line plots the same quantity in the case of the constant permittivity ~:= 80~,,.

varying 2 over more t h a n one decade and keeping all other parameters fixed. A suggestive account is given by the ratio cr

rg a

-/

dv

1 Jv

a

c dv

v

which yields the normalized average distance of the ionic cloud from the rod. Its behaviour versus ~ is plotted by the continuous line in Fig. 3. The dashed line plots the same quantity in the case when the permittivity is taken equal to the bulk value. The curves are practically coincident for low values of)~ but differ more and more as ~ increases. The case ~ = 0.59 A - 1, representing a DNA macromolecule, lies within the region where the non-linear effects cannot be ignored. These features are emphasized by Fig. 4 where the charge density Go and the electric field fla are plotted versus )~; as before, the dashed lines repre-

231

102

°o ~rOo 10

leo

yY .1

=

.05

.1

f

i

i

i

i i i

1

Fig. 4. Plot of the electric field, the charge density, and the relative dielectric perrnittivity versus the charge parameter 2. Continuous lines deal with a field-dependent dielectric permittivity; dashed lines represent the corresponding quantities in the case e = 80 ev.

sent the same quantities in the case when the permittivity is taken equal to the bulk value. Both Coand fla show a large range of variation when the non-linearity is considered. Again non-linearity effects become more and more remarkable as ~t increases. Both the curves of Fig. 4 and those of Fig. 3 show that the non-linearity of the medium can be ignored only when ~ is sufficiently low (less than 0.1, say), while the non-linear effects play an important role for highly charged rods. Conclusions

In this paper a variational approach has been elaborated for the determination of the electricfield in a medium with a dielectricpermittivity which depends on the electric field itself.Although the variational formulation is rather general, a detailed analysis has been performed for the particular,but prominent, model of rod-like polyelectrolytesin water solution described by the Poisson-Boltzmann equation. The electricfieldand the ionic distribution have been determined in various situations through a suitable application of Euler's method; for the sake of definiteness the physical parameters of the D N A macromolecule have been considered and a physically well-grounded function for the dielectricpermittivityhas been adopted. The electricfieldand the ionic distribution have been contrasted with the analogs in the case of

232 c o n s t a n t permittivity. T h e comparison shows the reinforcing effect of the nonlinearity on th e counter-ion condensation and the shielding effect on the electric field near the macroion rod. Finally, t he dependence of t he electric field on the p a r a m e t e r 2 (charge density) has been analyzed. T h e analysis indicates that, for values of ~ close to t h a t for D N A or higher, the non-linear effects c a n n o t be ignored, t hough t h e y occur mainly in a r a t h e r small region around the macroion rod. To our mind the results obtained in this paper might be viewed as a basic tool for setting up appropriate models of polyelectrolytes in various situations. Meanwhile, because of its generality, t he variational formulation is likely to be applied successfully to several problems where the investigation of nonlinear effects of th e p e r m i t t i v i t y is in order.

Acknowledgments T h e authors wish to t h a n k Prof. A. Chiabrera for having drawn their attention to the possible role of t he differential permittivity. T h is research was supported by the Italian C N R t h r o u g h the Research Project "Materiali e Dispositivi per l'Elettronica a Stato Solido".

References 1 C.J.F.BSttcher, Theory of Electric Polarisation, Elsevier, Amsterdam, 1952. 2 E. Clementi, Structure of water and counterions for nucleicacids in solutions, in: E. Clementi and R.H. Sarma (Eds.), Structure and Dynamics: Nucleic Acids and Proteins, Academic Press, New York, NY, 1983,p. 321. 3 F. Thoma, Th. Koller and A. Klug, Involvement of histone H1 in the organization of the nucleosomeand of the salt-dependentsuperstructuresof the chromatin, J. CellBiol.,83 (1979) 403. 4 F.L.Carter (Ed.), MolecularElectronic Devices, Marcel Dekker, New York, NY, 1983. 5 M. Parodi, Behaviour of mobile ions near a charged cylindrical surface: application to linear polyelectrolytes,J. Electrostat., 17 (1985) 255. 6 R. Fuoss, A. Katchalsky and R. Lifson, The potential of an infinite rod-like moleculeand the distribution of counter-ions, Proc. Natl. Acad. Sci., 37 (1951) 579. 7 I.M. Gelfand and S.V. Fomin, Calculus of Variations, Prentice-Hall, EnglewoodCliffs, NJ, 1963. 8 F. Bampi and A. Morro, The inverse problem of the calculus of variations applied to continuum physics, J. Math. Phys., 23 (1982) 2312. 9 A. Morro and M. Parodi, Extensive force in polyelectrotytesolutions, J. Mol. Liq., 32 (1986) 279. 10 W.F. Ames, Nonlinear Partial Differential Equations in Engineering, Academic Press, New York, NY, 1965, Chap. 5. 11 M.L.Krasnov, G.I. Makarenko and A.I. Kiselev, Problems and Exercises in the Calculus of Variations, MIR Publ., Moscow, 1975. 12 N.V.Kopchenovaand I.A. Maron, Computational Mathematics, MIR Publ., Moscow, 1975.