Force on two touching dielectric spheres in a parallel field

.s ;

Journal of

ELECTROSTATICS ELSEVIER

Journal of Electrostatics 35 (1995) 297-308

,

Force on two touching dielectric spheres in parallel field

a

Richard D. Stoy* Department of Electrical Engineerin9, Widener University, Chester, PA 19013, USA Received 16 November 1994; accepted after revision 23 January 1995

Abstract

An interactive linear multipole model is used to derive equations for the electric force between two uncharged dielectric spheres aligned in a parallel field. Up to 2800 multipoles are used to compute the force between identical touching spheres. Where necessary, extrapolation to an infinite number of multipoles is performed by an empirical equation. The model and equations apply equally well to a magnetizable sphere in a parallel magnetic field. The data presented allow the calculation of the force for values of electric (or magnetic) susceptibilities up to 21 500.

1. Introduction

There has been a great interest in recent years concerning interparticle forces in electric and magnetic fields. Several modern technologies utilize the chaining of magnetic particles in magnetic fields [1 ]. The understanding of many-particle systems begins with the systematic study of simpler systems. The two-sphere system has been elevated to the status of the basic building block in modeling the strong interactions of more complex particulate systems [ 1], The early analyses of systems with two uncharged dielectric spheres were performed in bispherical coordinates [ 2 - 4 ] . That coordinate system has the limitation that no two spherical surfaces can ever touch [4]. The present work is essentially a coordinate-free multipolar approach. Here we start from the known response of a dielectric sphere to a point charge to find the response to an arbitrary axial multipole. This response is combined with the dipole moment of a sphere which is its response to a uniform external electric field. With two spheres aligned in the field direction, matrix equations for the axial multipole moments in both spheres are derived. Then the formulas for the electric field and attractive force experienced by each sphere are calculated from the multipole moments. The case in which the external field is perpendicular to the line joining the centers of the two spheres is not treated in this paper because the multipoles * E-mail:[email protected] 0304-3886/95/$09.50 (~) 1995 Elsevier Science B.V. All rights reserved SSDI 0 3 0 4 - 3 8 8 6 ( 9 5 ) 0 0 0 0 3 - 8

R.D. Stoy/Journal of Electrostatics 35 (1995) 297-308

298

are not axial and the force is repulsive. With the maximum number of multipoles set at 2800, the force is numerically computed for evenly spaced numbers o f multipoles at relative susceptibilities up to 21 500 for touching spheres. These results apply equally to magnetizable spheres in a magnetic field. The actual force is obtainable only with an infinite number o f multipoles. This force is found by extrapolation using an empirical equation from numerical data for finite numbers (up to 2800) o f multipoles when the susceptibility exceeds 200.

2. Interactive multipole model for two dielectric spheres In this work the rationalized MKSA system o f units is used exclusively. Stratton [5] defines axial multipoles of the mth order in terms o f the potential as (~m

~---

1 (m) Pm(cOs 0) ~-~eP r-g;5 ,

m -- 0, 1,2 ....

(1)

where the Pro(COS0) are the Legendre polynomials in cos 0. Stratton [5] also defines an axial multipole o f order m in terms o f the axial multipole o f order m - 1 as

p(m) = mp(m-l)lm_l,

m = 1,2,3...

(2)

lm-I ~ 0 and p(m-1) ~ oc in such a manner that p(m-l)lm_ 1 remains finite. The conceptual construction of the axial multipole of order m and positive moment p(m) directed along the positive z axis occurs by placing a multipole o f order m - 1 and negative moment - m p (m-l) at the origin and an other multipole o f order m - 1 and positive moment m p (m-') a small distance lm 1 along the positive z axis. This is equivalent to a positive multipole o f order m - 1 at the origin and a negative multipole o f order m - 1 a small distance lm-1 along the negative z axis. Consider a dielectric sphere of radius Ra and permittivity s] whose center is located at the origin o f a spherical coordinate system. The medium outside the sphere is characterized by its permittivity /33. We locate a point charge q = p(0) at z = D > R 1 on. the positive z axis. According to Stratton [5] the potential at any point outside the sphere is when

1

q

q~+ = ~b0 ÷ ~bli - 47z~3 r2

q

oo ~~

-/'/(/31-~-~3)-

R~n+lrn(c°sO)

(3)

4~z~;3nz==0n~:l ÷ (rt ÷ 1 )/33 D n+l n +r~ '

where r2 is the distance from the point o f observation to the location o f the point charge. Let us model a point dipole as two displaced point charges. The dipole +p(1) is located on the z axis at z = D: the charge - q = _p(0) is located at z = D and the charge + q = +p(O) is located at z = D + d. The dipole-induced potential outside the sphere is now

qb~(1) _ p(l) ~ 4roe3

n(n + 1)(el - s 3 ) R ~ n+l Pn(cosO) r n+l n=0 nel + (n + 1)-~3 ~

(4)

299

R.D. Stoy/Journal of Electrostatics 35 (1995) 297 308

N o w we place an axial quadrupole p(2) on the z axis at z = D: a dipole - 2 p {~) at z = D and a dipole + 2 p (|) at z = D + d. The quadrupole-induced potential outside the sphere is /?2n+l n(n + 1)(n + 2)(el - e3)~-i

(~+(2) __ _p(2) 2(4rc~3 )

n ~ ~ ( n ~ i)-e3

n=0

D "+3

P,(cos 0)

(5)

r "+'

An axial octupole p(3) is now positioned on the z axis at z = D: an axial quadrupole - 3 p (2) at z = D and an other axial quadrupole 3p (2) at z = D + d. The resulting octupole-induced potential is (D; ( 3 ) -

pz,+| Pn(cOs 0) n ( n + 1)(n + 2)(n + 3)(el -- ~33)*'1

__P(3)

6(4r~83)

ne,| + (n + 1)e3

n=0

D "+4

(6)

r n+|

An axial multipole p(rn) is positioned at z = D. For the multipole-induced potential outside the sphere we can write o+(m ) _

(--1)rap (m)

n(m + n)!

4roe3

m! n!

n=0

v2.+| P , ( c o s O) el -- e,3 "'l n81 + (n + 1 )e3 D m+"+l r "+|

(7)

An axial multipole p(m) located at z = - D can be generated from an axial multipole ( - l ) m p (m) located at z = D coupled with the coordinate transformation 0 --~ 0 + g, and noting that cos(0 + g) = - cos 0, and P r o ( - cos 0) = ( - 1)rePro(cos 0). Thus, the potential outside o f sphere 2 o f radius R2 and permittivity e2 caused by an axial multipole p(m) located at z = - D is oc

,h +(m) --

p(m) Z ( _ l ) n n ( m + n )

•e2

4roe3

n=0

!

m! n!

/?2n+l ez--e3 "'2 nc 2 + (n + 1 )e3 D m+'+l

P.(cos O) r "+|

(8)

If we are concerned only with the potential and field in the region o f space which is outside o f a sphere, we can say that the sphere responds to an axial multipole o f order m at z = + D with an infinite series o f axial multipoles whose order ranges from 1 to oc which are all located at the center o f the sphere. N o w let us place sphere 1 on the z axis with sphere 2 also on the z axis but displaced a positive distance D from sphere 1. Now let both spheres be immersed in a uniform electric field o f strength E0 in the positive z direction. The response o f each sphere, initially neglecting interactions, is that o f a dipole also oriented in the positive z direction. The non-interacting moments o f sphere 1 and sphere 2, p(|~) and /020, -(|) are

p(ll0) = 4:gg,3~ R ~ E 0 ,

81 -l- Z£ 3

p(|) = 4ge 3 ~ R 3 E 0 20

~2 ÷ 293

.

(9)

300

R.D. Stoy/Journal of Electrostatics 35 (1995) 29~308

When D is finite the spheres see each other and each sphere has located at its center an infinite number of multipoles. The multipole moments of sphere 1 caused by the external field and the multipole moments of sphere 2 are £

_(1)

~

=oJ.v,0 +

m=O

,gl -- go3

(m)n( m + n)!

xm

- ( - 1 ~ p:

mi;~i

R 2n+l

nel + (n + 1)~3 D m+n+l '

n=0,1,2...

(10)

The Kronecker delta finn has the value 1 for m = n, and the value 0 for m ¢ n. The multipoles of sphere 2 caused by the external field and the multipoles of sphere 1 are

x _(1)

~

~2 -- ~3

.n (m)n(m+n) !

----~qnY2o +

-(-1)

Pl

mini

m=O

R i n+ I

n/32 q- (/l -Jr- 1 )g3 D m + ' + l '

n =0,1,2...

(11)

Initially, the only non-zero multipoles are P(l1) = ~'10"(')and p~l) = p(l~). Subsequently, the operations described by Eqs. (10) and (11) are repeated until the values of the multipoles have converged sufficiently. Matrix inversion methods can also be used; these will be discussed later. Eqs. (10) and (11) bear a significant resemblance to the integral Eqs. (18) and (19) of Lindell et al. [6]. Those workers applied the electrostatic image method to two dielectric spheres in a homogeneous electric field. The z component of the electric field in the region of sphere 2 caused by the extemal field and by all the multipoles of sphere 1 is directly calculated to be

1 ~

(m + 1)p]")

Ez,2 = Eo + ~

rm+2

(12)

m=O

The nth partial derivative of Ez,2 with respect to z is then ~3nEz,2 _ ( - 1 ) n ~

4g~;3

(~zn

m=0

(m + n + 1)!P(lm) m i r ~-+~~'~ "

(13)

Similarly the z component of the electric field at the center of sphere 1 caused by the external field and by all the multipoles at the center of sphere 2 is

-1 ~

(-1)m(m+l)p~ m)

Ez, I = Eo + 4-~e3 L m

rm+2

(14)

0

The nth partial derivative of Ez, 1 with respect to z is OnEz, n

_

4Its3--1 ~ m=O

(--1)m(m + n + 1)! P2-(m)

,,.2-,r ~ 2 -

( 15 )

R.D. StoylJournal of Electrostatics 35 (1995) 297-308

301

The electric force on axial multipoles using the definitions presented here is given in terms of the derivatives of the electric field as given by Jones [7]: ~

F~=

p(n) c~Ez n! ?z"

(16)

n=0

Combining results, we find that the electric force on sphere 1 due to the field of sphere 2 is FI = ~

~

m=On=O

- ( - 1 ) m (m + n + 1 ) , p~n)p~m) 4rc/~3 m! n! D re+n+2 "

(17)

While the electric force on sphere 2 due to the field of sphere 1 is F2 = ~ m=0n=0

( - 1)" (m + n + 1)! Pl(m)P2(n) 47ze3 m! n! O m+n+2 "

(18)

All product terms in Eqs. (17) and (18) represent an attractive force. These two equations resemble Eq. (16) of Fowlkes and Robinson [8] for the solution of a somewhat different problem: the electrostatic force on a dielectric sphere resting on a conducting substrate. In principle there are no restrictions on ~1 or e2, except that no pole can be introduced into Eqs. (10) or (11). Letting el, /~2 ~ 0(3 produces the case of conducting spheres. For the equivalent magnetic case (el -~ #~, e2 ---' #z, e3 ~ /~3, and Eo ~ Ho) the superconducting limit can be obtained as p~, /~2 --~ 0.

3. Matrix representation of multipoles and solution The infinity symbols in Eqs. (10), (11), (17) and (18) must be taken literally. That is, the true value of the force is obtained only with the two infinite sets of multipoles which satisfy Eqs. (14) and (15). However, any realistic calculation of the multipoles and the force can involve only a finite number of multipoles. Let us define N as the maximum numbers of multipoles in any particular calculation. This means that all multipoles of order > N are identically zero. In the following, the matrices A, A1 and A2 are N by N; the column vectors P, PI and Pz are of length N; while the column vectors P0, P10 and /'20 have non-zero first elements with all other elements set to zero. The identity matrix is denoted by 1 and the inverse of any square matrix A is designated A -1 so that AA -1 = I. Eqs. (10) and (11) can be cast in the form: Pl = PIo + A1P2,

P2 = P20 + A2P1.

(19)

The closed-form solutions are PI = ( h i 1 -- A2)-I(P2o + A~-IPIo),

(20)

P2 = (A2 ! - ,41 ) - I ( P l o + A21P20)-

(21)

R.D. Stoy/Journal of Electrostatics 35 (1995) 297-308

302

The iterative solutions are Pl,new = PI0 + A1P2, old,

(22)

P2, new ----P20 + A2PJ,ol~.

For identical spheres [el = e2, Rj = R2, p~m) = ( - 1 ) m p l m)] we find that Eqs. (10) and (11) can be written as (23)

P=Po+AP. The solutions, closed-form and iterative, are

P = (I - A ) - ' Po,

(24)

Pnew = Po + APold.

4. Qualitative numerical results for identical spheres Let N again be the number o f multipoles used for values of/3~, R1 and D.

0 p(lm) ~3N

> 0

for m < N

and

OF1

~-

> 0.

(25)

As N is increased each multipole asymptotically approaches its true value, as does the force. It has also been found that (m), (re+l)]

PZ /Pl c~N

[

> 0

for

(m + 1)~
(26)

c3F1 < 0, ~3D

~pl '~) - < 0. c3D

(27)

~3F~

c~p~m) - > 0.

(28)

> 0,

Any decrease in the distance between the spheres or any increase in the permittivity o f the spheres increases both the magnitude o f each multipole and the force. The observations made here were crucial for the developments in the next section.

5. Force between identical touching spheres This a case o f great interest which cannot be solved in bispherical coordinates. The numerical solution o f Eq. (24) to yield data on the force has been carried out over a wide range o f electric susceptibility Zc (=/31/e3- 1): from 10 -10 to 21 500. These data are presented in Fig. 1. The methods presented in this work have been used to produce numerical data only for positive values o f Ze. Negative values o f Xe could be used provided that no singularity is introduced into Eqs. (10) or (I 1 ). Normalization was carried out as follows: R1 = 1 m, E0 = 1 V/m, D = 2Rl, and 133 1 F/m. =

R.D. StoylJournal of Electrostatics 35 (1995) 297-308

303

10

lO

10 s

Lr-

10 ° (9

.o o 1 0 .5 E "~-

gl0

-1(~

0 _J

10 "15

Y

1 0 2c

10 "~ 10 "1°

10 -s

10 0

10 s

Logarithm of susceptibility Fig. l, Log-log graph of the attractive force between two touching spheres as a function of the susceptibility (10-1°~
5.1. Values of susceptibility less than 1 In this range of ge only a few multipoles are necessary to find the normalized force. The initial slope of 2 in the log-log presentation is due to the Pl0-P20 (1) (~) interaction. The normalized force in Newtons is given simply by F ~ ( z e ) = 0.5236Z~

(N).

(29)

5.2. Values of susceptibility less than 100 In this range of ge it was always possible to achieve numerical convergence of the force (values of the force which are close enough to their asymptote) by using a sufficient number of multipoles, Equipment limitations kept the maximum number of multipoles at 2800. These data are shown in Fig. 2. In the author's attempts to find an equation for the force as an asymptotic function of the number of multipoles, success was achieved using a constant plus a sum of power (or exponential) functions. Eq. (30), which follows, fits the data exactly in the numerical sense. There is no apparent physical reason why Eq. (30) works with

R.D. Stoy/Jow:nal of Electrostatics 35 (1995) 297 308

304 10 3

: i i i iiJ-:,

I

I

I I ;A-~tl

10 2 t

t

8

J

I

/f'f .....

10 l

,,9

//

~5

~1o

Itl

o

, I

I I ] I

,I t"

IP"

t. tll

8*

IX

10"

I

I}1 Itl

I I IJl I I III II III

.tr

11" I

/I 10 .2

,e'~

10

0

10 Logarithm of susceptibility

II

Ill

11

I11

I I II]11 I i I I Ill I I [1111

I

10 3 10 -1

tJI

1 1 t I~1 I 1 1[1

I.,,4"I

.d

. ]!,!!!

.

1

I II I l l l

ld

Fig. 2. Log-log graph of the attractive force between two touching spheres as a function of the susceptibility ( lO-I ~Ze ~ 102) .

such precision, but it does. In the present range o f Z~, Eq. (30) predicts the values o f subsequent data points with ten-digit accuracy. The normalized force was computed in increments of 50 multipoles from N = 0 up to 2800 multipoles in each sphere. Here we define M -- N/50 where M goes from Ml (>~0) to M2 (~<56). The resulting data was used to determine the constants on the right-hand side o f the following equation: k

F,(ze,MI + j ) = Foo(Ze) + Z fi/~/,

k = ½(M2 - M~ ),

j = 0, 1,2 .... 2k. (30)

i=l

The term Fo~(Ze) is the asymptotic force for N ~ oo and also for j ~ oo with IBi] < 1. The method for determining the 2k + 1 coefficients on the right-hand side of Eq. (30) from an equal number o f evenly spaced data points on the left-hand side is given by Hamming [9]. An example is given in Appendix B. The researcher does not know beforehand the minimal value of k. Selecting too small a value for k produces scattered results, while a larger value for k (i.e. k --~ k + 2, k + 4, etc.) produces additional complex conjugate pairs of both fi and Bi without any significant corrupting effects on Foo(Z¢). In fact, the effect is more to stabilize the value of F ~

R.D. Stoy l Journal of Electrostatics 35 (1995) 297308

probably due to the inclusion of more in Ft. The numerical stability of Fo~(Z~) for data points solidified the author's belief force to an exceptionally high degree of

305

data points with a greater spread of values different sets of 2k + 1 out of a maximum 57 that Eq. (30) is a valid representation of the precision for Z~ < 100.

5.3. Values of susceptibility 9reater than 100 In the previous section we saw that Eq. (30) represented force data for ge < 100 and N ~<2800. There is no evident reason why this equation should suddenly become invalid for ~(e values > 100. In this range of Ze the use of Eq. (30) to compute Fo~(Ze) involves extrapolation. In general extrapolation is a leap into the unknown. Here the data fits the model so precisely that we can extrapolate with confidence. It is only in the decade between Ze - - 104 and X~ -- 105 that the limited range of force data precludes extrapolation to consistent and stable values of F ~ . Only stable values of Fo~ are reported in this work. The greatest value of Z~, 21 500, produces an extrapolated value of force which is about 12.2 times the value of force for 2800 multipoles. Fig. 3 shows the relevant data. 7

110 2

10 3

10 4

10 s

Logarithm of susceptibility Fig. 3. Log-log graph of the attractive force between two touching spheres as a function of the susceptibility (102~
306

R.D. Stoy/Journal of Electrostatics 35 (1995) 297-308

5.4. Discussion of numerical data The normalized force data can be used to find the force for arbitrary values of ~33, E0 and Rl. Noting Eqs. (9)-(10), we see that each multipole p(ln) is proportional to e3EoR~l+2 and that each term of the force in Eqs. (17) and (18) is proportional to p(m) _(n)/. Dm+n+2or e3E~R~ and that the force in Newtons is given by

1 //2 /g3~o,1

F(Ze'e3'E°'RI)=F~(Ze)\

(83, actual'~/E0, actual'~2/R1 actual)2 1F/m } \ l V / m J

\

lm

(N).

(31)

The case of two magnetic spheres aligned in a uniform external magnetic field can be treated as the dielectric case has. Here el ~ #l, 83 ~ ,//3 and Eo ~ Ho. The magnetic susceptibility ;(m ( = # 1 / / 2 3 - - l ) exactly replaces the electric susceptibility ~e and

F~(Zm ) = F~(Xe):

F(Zm'#3'H°'RI)=F~(Zm)

(]/3, actual) (H_0,actual'~2(RI actual)2 1H/m

\ 1A/mJ

\

lm-

(N).

(32)

We should remember that no effects of dielectric or magnetic saturation can be included here since we are always assuming that e~ and /tl are constants and are not functions of E0 and H0 respectively. The greatest value of Ze or Zm, 21 500, seems sufficient to encompass most linear dielectric polarization or magnetization phenomena. Three graphs of F ~ as a function of Zc are given. The normalized force F ~ is basically a power function of the susceptibility. The entire data set seems to indicate that the attractive force between two touching spheres is unbounded as Xe or ~m ----+ ~ . The extension of this argument to the case of touching conducting spheres is accompanied by difficulties with the boundary conditions at the point of contact as discussed by O'Meara and Saville [10].

6. Conclusion The interactive multipole model has shown its power by allowing us to calculate the force between two touching spheres in a uniform electric or magnetic field for values of susceptibility up to 21 500, and for all values of field strength and radius. The data show that the force is a power function of the susceptibility with the exponent ranging between 2.00 and 1.54.

Appendix A. Computational aspects The numerical solutions of Eqs. (18) and (24) were performed on the an IBMcompatible personal computer (PC) which contains a Datatech mainboard (PLM-33001 ) with an Intel i80486/DX2 66 central processor unit (CPU) and a Weitek 4167-33 coprocessor. The mainboard holds 65.536 Megabytes (MB) of dynamic random-access memory (DRAM) to accomodate the 2800 x 2800 array of double precision values

R.D. StoylJournal of Electrostatics 35 (1995) 297 308

307

(62.72 MB). The program to solve Eq. (18) was written in the Fortran 90 language using a Microway NDP486 compiler. Lower-upper triangular (LU)-decomposition was found to be the best method for Eq. (24) when Ze exceeds 1000. Double precision and the i80486's own coprocessor were used throughout. The production o f the 57 values o f force for use in Eq. (30) requires about 36h. The solution o f Eq. (30) was developed on the same PC and transferred to a Cyber 930 mainframe computer which has the equivalent o f quadruple precision on a PC.

Appendix B. Example of curve-fitting procedure Consider that we have a set o f seven data points which are evenly spaced in the independent variable. These seven data points could be a subset o f a larger set o f points. Using Eq. (30), we get

FI(ze,Mj + j ) = FI0, F l l , F12, FI3, F14, FIS, FI6; k =3,j

= 0, 1 , 2 , 3 , 4 , 5 , 6 = 2k.

(B.1)

This is equivalent to writing the following: ~" B 0, Fio = F .~ + f i B ° + f e B ° + j3 3

F1, = F ~ + f , B I + UzB~ +.f3B~, F,2 = F ~ + f,B~ + f 2 e 2 + f3B~, F,3 = F ~ + f lB~ + f zB 3 + J'3B~,

(B.2)

F14 = F ~ + U,B 4 + f e B 4 + f 3B4, F,5 = F ~ + U,B~ + f 2B59+ f 3B~, Fi6 = F ~ + f l B

6 + f 2 B 6 + f 3 B6.

The extension to any odd number o f data points is obvious. The number o f data points is always equal to the number o f constants to be determined. Since the IBil are all less than 1, F ~ is the asymptotic value as j ~ ~x~. Without a theoretical basis for Eq. (30) it is not possible to assign a physical significance to the f i and the B i.

References [1] C. Tan and T.B. Jones, Interparticle force measurements on ferromagnetic steel balls, J. Appl. Phys., 73 (1993) 3593-3598. 12] A. Goyette and A. Navon, Two dielectric spheres in an electric field, Phys. Rev. B, 13 (1976) 43204327. [3] R. Ruppin, Surface modes of two dielectric spheres, Phys. Rev. B, 26 (1982) 3440-3444. [4] R.D. Stoy, Solution procedure for the Laplace equation in bispherical coordinates for two spheres in a uniform external field: Parallel orientation, J. Appl. Phys., 65 (1989) 2611-2615. [5] J.A. Stratton, Electromagnetic Theory, Vol. 174, McGraw-Hill, New York, 1941, p. 204.

308

R.D. Stoy/Journal of Electrostatics 35 (1995) 297-308

[6] I.V. Lindell, J.C.-E. Sten and K.I. Nikoskinen, Electrostatic image method for the interaction of two dielectric spheres, Radio Sci., 28 (1993) 319-329. [7] T.B. Jones, Dielectrophoretic force in axisymmetric fields, J. Electrostatics, 18 (1986) 55~52 [8] W.Y. Fowlkes and K.S. Robinson, The electrostatic force on a dielectric sphere resting on a conducting substrate, in: K.L. Mittal (Ed.), Particles on Surfaces: Detection, Adhesion and Removal, Plenum Press, New York, 1988, pp. 143-155 . [9] R.W. Hamming, Exponential Approximation - Sums of Exponentials, in: Numerical Methods for Scientists and Engineers, Dover, New York, 2nd ed., 1986, pp. 617-627. [10] D.J. O'Meara and D.A. Saville, The electrical forces on two touching spheres in a uniform field, Q. J. Mech. Appl. Math. 34 (1981) 9-26.