Rigorous derivation of hydrostatics and the separation of phases

Rigorous derivation of hydrostatics and the separation of phases

Physica 57 (1972) 267-280 o North-Holland Publishing RIGOROUS AND DERIVATION THE Co. OF HYDROSTATICS SEPARATION OF PHASES D. J. GATES Mathemati...

718KB Sizes 0 Downloads 16 Views

Physica 57 (1972) 267-280 o North-Holland Publishing

RIGOROUS AND

DERIVATION THE

Co.

OF HYDROSTATICS

SEPARATION

OF PHASES

D. J. GATES Mathematics

Defiartment,

Imperial College, London,

S. W. 7, England

t

Received 17 March 197 1

Synopsis The hydrostatic equation for the density n(y) at a point y in a fluid subject to an external potential #(y) can be written in the form n(y) =

POCU -

$4391,

where p&) is the density of a uniform system with chemical potential ,u. It is shown how this equation can be derived rigorously from the statistical mechanics of a general system of particles, with a two-body potential having a hard core. The result holds in the limit y --f 0 if the external potential is #(yx) and if m(y) is suitably defined. For the gas phase the appropriate definition is n(y) = lim,,e ni(y/y, y) where ni(x. y) is the usual one-particle distribution function. To prove the result generally, a locally spaceaveraged density is introduced.

1. latroduction. Suppose that a uniform fluid in equilibrium with density p and fixed temperature (not shown in the notation) has pressure $0(p). If an external field, with potential y(y) at the point y in the fluid, is applied, then the well-known hydrostatic eqztation. for the number density n(y) at y is ;

;j-

fio(12)= -

$ y.

(1.1)

This is essentially an experimental law. It implies that when there are two fluid phases (e.g. gas and liquid) then the field separates them spatially. In this paper we show how the equation can be derived rigorously from the principles of statistical mechanics. Kac and Thompsonl) have previously given such a derivation for the case of a one-dimensional van der Waals fluid in a uniform field. First we transform the equation into a more convenient form. If ,uo(p) denotes the chemical potential of the uniform fluid with density p, then @o@)/ap = &&)/a,. Hence (1.1) can be written t

NOW

at Mathematics Dept., Rockefeller University, New York, N.Y. 10021. 267

268

D. J. GATES

I

1

ri

*r

Fig. 1

which yields

Pd4Y)l + Y(Y)= P7

(1.2)

where .,uis a constant independent of y. If pa is the inverse function of ~0, we obtain finally the explicit formula for a(y) :

4Y) = POb - Y(Y)l.

(1.3)

For values of p where two phases are present in the uniform system, the graph of the function ,u&) has a horizontal segment, at height ,,~i say. Hence P&A) has a jump discontinuity at ,u = ~1 (fig. 1). It follows from (1.3) that if ]yI is large enough then N(Y) has jump discontinuities along surfaces in the fluid given by lu -

Y(Y) = AJ1.

(1.4)

Thus the phases are separated so that the liquid is in the potential ‘valleys” and the gas on the potential “hills”, with the gas-liquid surfaces located by (1.4). In particular if y = gyi, which represents a uniform gravitational field in the direction of the yi component of y, then the graph of the density n(yl) has the same shape as in fig. 1, with the surface at yi = (,u - pi)/g. In what sense can we expect to prove (1.3) rigorously? The work of Kac and Thompsonr) and Lebowitz and Percuss) indicates that (1.2) holds only if y(y) does not vary significantly over distances of order of the molecular size [This indeed is the experimental situation on which (1.1) is based]. Otherwise there are corrections to (1.2). This suggests we use a potential of the form y(yx) at point X. If y is small, this potential varies slowly with X. We might therefore expect (1.3) to hold exactly in the limit y + 0. To make this work we must clearly look at the density of the system on the large scale; i.e. over distances of order y-i, where y(yx) varies significantly. Hence the natural definition for s(y) in statistical mechanics is

(1.5)

HYDROSTATICS

AND

where %1(x, y) is the one-particle potential

I&X).

THE

SEPARATION

distribution

OF PHASES

269

function at the point x for the

We must also make the system infinitely

large in order to

obtain (1.3) exactly. Otherwise we will get corrections due to the presence of the container walls. Hence nr(x, y) must be a limiting function as 191+=co tJzermodynamic limit) where IQ1 is the volume of the container $2. This approach seems hopeful at first, but there are some difficult mathematical obstacles to a rigorous proof. First, no one has yet proved the existence of the limiting function ni(x, y), except at low densities (i.e. for a gas), and for the special case of a lattice gas with attractive interactions: Griffithsa) and Lebowitz and Penrosed). One would expect a general proof of the existence of ni to be very difficult for the reasons given by Fishers). (A proof of the existence, in a restricted sense, has been given recently by Ruellea).) The second obstacle is that n(y) has jump discontinuities when (1.4) holds, so it is not unique here. Hence we cannot expect the limit in (1 S) to exist when (1.4) holds. If we restrict ourselves to the gas phase then neither of these obstacles occurs and (1.3) can be proved. We do this in the next section. To obtain the result generally, we introduce a new definition in place of s(y). Let %1(x, Sz, y) be the one-particle distribution function for finite volume 101 and external potential y(yx). Let w(x) be a cube of volume 101 centred at x. Then the quantity fii(x> Q, Y, w) = (‘““‘)~j+d~‘~i’“>

Q> y),

(1.6)

represents the average density in m(x). The existence of the 191 -+ 03 limit of ffi is also too difficult to prove, but since nr 2 0 and ni is bounded above uniformly in IQ1 s), it follows that the limit infimum and the limit supremum of +ii both exist. For definiteness we choose the latter and define fir(x, 00, y, w) = lim sup #1(x, Q, y, w). IDI-=

(1.7)

By the previous argument we now consider the limit y + 0 of %(y/“~, 00, y, 0). We cannot prove the existence of this limit, but the existence of the limits infimum and supremum follows from the boundedness of Cr. For definiteness we choose the latter. Finally, to ensure that w contains a large number of particles, we let 10) --f co and define 9?(y) = lim lim sup %r(y/y, co, y, w), InI+= Y+O

(l-8)

which we call the locally space-averaged density. does exist as we subsequently show. In defining K(y) we have used the three limits first (~21+ co, then y + 0, and finally

The final limit

IwI + 00.

101 + 00

(1.9)

270

D. J. GATES

Thus we are looking at the case where (cf. Lebowitz

and Penrose7), (1.10)

where v is the number of dimensions and ~0 is the hard-core diameter. The size of the cells m is much less than the distance y-i over which y(yx) varies significantly, which in turn is much less that the size of the container. The order of these limits is important. For example, if we interchanged the limits IwI -+ 00 and y -+ 0, the result would clearly be independent of y. The main result of this paper (section 3) is to show that s(y) = p&$-y(y)] for all values of y except those on the surfaces (1.4). The result holds even for crystalline phases. We assume throughout the paper that y(y) is continuous and bounded below: (1.11)

Y(Y) 2 -% where $Gis a positive constant. core potential p(r) such that

PI 2 10,

for where C,

E

hard-

Irl -=c10,

for

v(r) = *

The particles interact via a two-body,

(1.12)

and YOare positive constants.

2. The gas phase. In this section we show that (1.3) holds for low densities using the simple definition (1.5) for n(y). Here, we define the infinitevolume, one-particle distribution function 1zi(x, y) by its fzqacity qbansiolz*) %(X>

Y) =k&3(%Y)1

(2.1)

where .zk ck _

gk(xl, y) =

in which 2

z

fjmv &+

3

*)!

s

dx2 ... dXkUk(Xl

k

... xk) exP[-~i~l

y(rW, (2.2)

(2.3)

A is the thermal wavelength, /? the reciprocal temperature and the Uk’s are the Ursell cluster functions*). The constant p has been introduced immediately. One could choose to interpret ,u as a sort of bulk chemical potential and x as the corresponding fugacity. It should be noted, however, that the usual definitions of bulk thermodynamic functions (e.g. free-energy density) in the 101 --f 00 limit, do not work here unless we impose further restrictions on y, such as making it periodic.

HYDROSTATICS

AND THE SEPARATION

OF PHASES

271

To obtain the y -+ 0 limit of mi(y/y, y) we show that its expansion converges uniformly in y for small 1~1.We have

x where

s

drz . . . drkVk(r2 . . . rk) exp[-i$4~

Vk(r2 . . . rk) = Uk(X, x +

r2,

.... x +

rk),

+ rrdll

(2.4 P-5)

which is independent of x. It follows that

(2.6) As shown by Penrose*) j drs . . . drk lvk(rs . . . rk)l < (k e2B@‘)k-sBk >

P-7)

where -2@’ is a lower bound (ip’ > 0) on the potential energy of a single particle due to v(r), and

Hence we have

where Mk = (lzl e&)k (K ezS”)k-2Bk/(R -

1) !.

(2.9)

The series & Mk has the radius of convergences)

[Bexp(1 + j$ + 2P@')l-1e Hence, by the Weierstrass M-test, the series for rti(y/y, y) converges uniformly in y for IzI < [R exp(1 + BV + 28@')1-1.

(2.10)

Small 1.~1 implies small 1~1,so that we are restricted to low densities. Now, from (1.5) and (2.1) we have (2.11) when (2.10) holds. Returning to (2.4), we note that exp[+3 & ~(y + yrt)] is bounded uniformly in y, and from (2.7) that Tr, is absolutely integrable. It follows from Lebesgue’s bounded convergence theorem that in taking y + 0 in (2.4) we can transfer the limit inside the integration. From the

272

D. J. GATES

continuity

for all t-8, giving

of y, we have y(y + yrc) + y(y)

~+W(Y/W)

(2 e-s+CY’)k @ _ 1), .s

=

The density function expansion s,

p&)

dr

2 . . . drk vk.

of a uniform

s dr

2

. . .

(2.12)

system

(y = 0) has the fugacity

drk vk.

(2.13)

Combining this with (2.11) and (2.12) gives finally the desired result (1.3), subject to the low-density condition (2.10). There is, of course, no separation of phases under this condition since only the gas phase is present. 3. The general case. In this section we prove ( 1.3) for a general state of the system, using the locally space-averaged density 2(y). We use the grand canonical ensemble. Here, the one-particle distribution function for finite volume Q is defined by

ss dx s . . .

wN(Xi . . . XN) =

z

V(Xt -

lSi
is the

XJ) +

eXp(+wN),

(3.1)

R

la

where

dXN

c

l
Y(YXi),

(3.2)

energy

-N;o$idX~ ...~d.eXp(-~~N), n

P-3)

R

is the grand partition function (GPF). Here NC(Q) is the maximum number of particles that 52 will contain without overlapping of the hard cores. For any region R in Y-dimensional space, we define the characteristic function

XR(X)=

l 0

for x E R,

(3-4

otherwise.

We consider the influence of an additional external potential Axmu which represents a uniform field in the cell 0(x’). This system has the grand partition function E(,u, 0, y, y + AX+,,). Our method is based on the observation that --

1

Sogq..., B 14 [ an

Y + kVcx9)

1 a=0

which follows from the above definitions

= ?qX’,

. ..) y),

P-5)

and (1.6). We use this to obtain

HYDROSTATICS

upper and lower bounds

AND

THE

SEPARATION

OF PHASES

on 51, and show that the bounds

coincide

273

in the

triple limit (1.9). The function log E(. . . , y + ;Z~+~,) is convex in 1, as one can show by taking its second derivative. But for any differentiable convex function f(A) we have for all A > 0 V(O) -

fc-4lln

(3.6)

I f’(O)I [f@) - f(W~.

Hence, from (3.5) we obtain 6(x',

***,y) 2

&

[log q..., y) -

1%

q..., y +

3LXO(d))l.

(3.7)

We assume 0(x’) is entirely inside 0, and denote its complement with respect to Q by 9 - 0(x’). We define the configurational partition function (3.8) Each integral over 52 can be written as the sum of two integrals, over 0(x’) and 52 - 0(x’), respectively. This gives (omitting x’ from the notation)

xfdr.+r...[

dxlv

n-o I

exp(+wN)

o

w

n-a,

eppw-n~;-npA

Qn(w y) QN-&~ - w, y),

(3.9)

where W-(m) is a lower bound, for all PZ,N and 52, on the energy of interaction due to v(r), between o and Q -

3

cc).Since (3.10)

=NF;NQ~,

it follows from (3.9) that Q,

Q, y + Axa) I emBw-=(p - 1, 0, y) =(p> LJ- w, y).

(3.11)

This gives

log q/4 G, y + &o) I log qp - A,m,y) + logw4

J2- 0, y) - BW-,

(3.12)

which provides a lower bound on the second term in (3.7). To obtain a lower bound on the first term in (3.7) we introduce a cube o’(x) centred like w(x) at x, similarly oriented but slightly smaller. The region

274

D. J. GATES

w - o’ is a shell surrounding this shell from Sz, i.e.

w’. Let 9’ be the region obtained by removing

Since 51’ C Q it follows from (3.8) that

Q&Q, Y) 2 Q&Q', 7~).

(3.13)

Each integral over 8’ in QN(SZ’, y) can be written as the sum of two integrals, over o’ and Sz - o respectively, which gives, as in (3.9),

Q@', y) 2 edBw+c"Q&J, y)QN-&J - w,y),

(3.14)

n=O

where W+(m, w’) is an upper bound for all n, N and Sz on the energy of interaction due to y(r), between o’ and Sz - o. Combining (3.10, 13 and 14) gives log qcl,

BW+.

(3.15)

co, y) + /VV--#?W+].

(3.16)

JJ, y) r log =(#K w’, y) + log qp> Q - w, y) -

Substituting

(3.12) and (3.15) in (3.7) gives

[log q/A , co’,y) -log+-1, 2 -LBnluJl

Now we take the succession of limits (1.9). The right-hand side of (3.16) is independent of D so that the limit (JJ +co does not effect it. Next we put x’ = y/y and consider the limit y + 0. We note that, by a change of integration variables,

s/-p,o’(y/y),

WI

= %4 w’(O)*Y(. +Y)lt

where co’(O) is centred at the replace y(yxg) by ~(7x6 + y) in are bounded uniformly in y, so integrals. Since y is continuous lim Y(Y~$ + Y) = w(u)

(3.17)

origin of coordinates, and y( *+y) means (3.2) and (3.3) for all i. The integrands in E that the limit y --f 0 can be taken inside the we have for all $2,

Y-4

and hence lii

q/J, O’(O)> y( * +Y)l

= So[ru -

Y(Y), w’(O)l,

where & is the GPF in zero external potential.

Treating

(3.18) the second term

HYDROSTATICS

AND

THE

SEPARATION

OF PHASES

275

in (3.16) in the same way we obtain limy s;p

+

l:;zp

%(y/~, i.4Q, Y,y)

{log Eo[p 2 --!-Bnbl

Y(Y), 0'1

log Eo[p - il - y(y), WI +

-

gw- -

BW+).

(3.19)

Next we let (WI + co. The pressure $0(p) in the grand canonical ensemble is defined by p,(p)

=

lim - ’ If+=J e-4

(3.20)

log =olp, Q),

which is known to exists). The bounds W+ and IV- on the interaction energy can be found by the method of Lebowitz and Penrose7) as shown in the appendix. The significant facts are that w-(~)/lQJl

(3.21)

as 10.11 + co,

+ 0

and w’ can be made to depend on w in such a way that

bllb’l + 1, and

Together

(3.22)

as IwI --f 00.

W+(e), W’)/l~l + 0

with (3.19) and (1.8) these imply that

fi(Y) 2 ;

{Po[cl -

W(Y)1 -

Po[cl -

for all il > 0. Since #M&U)is convexs)

Y(Y) -

it followslo)

(3.23)

U* that

POW = 3$0(/J)/% exists except at a countable in (3.23) therefore gives

G(Y) 2

POCU - Y(Y)1

number of values ,ur, ,u2, . . . of p. Letting 1 + 0

for lu -

y(Y) # ~“1,~2,

...

(3.24)

An upper bound on Fir can be obtained by the same method, using this time the left inequality in (3.6). This leads to (3.24) with the inequality reversed. Together these yield our main result fi(Y) = POCU-

Y(Y)1

for ru -

y(y) # ~1,~2,

.-..

(3.25)

4. Discussion. Again we note that (3.25) holds for all states of the system including crystalline states, so that (1.1) also holds for all states. However,

276

D. J. GATES

for a crystal po[+i(y)] no longer pressure is a tensor quantity.

represents

the pressure

at y. The true

The reason that (3.25) holds generally is that the space-averaging over o smoothes out any periodic variation in density on the molecular scale. To understand this one can consider the artificial case nr(x, y) = f(x) g(yx) where f and g are continuous and f is periodic with space-average unity. Then one can easily show that 3(y) = g(y), so that the periodicity due to f has disappeared. We note that G(y) corresponds more closely to the physically measurable density than does ni(x, y). For a physically realistic choice of p)(r) we expect (but cannot prove) that there will be just three phases: solid, liquid and gas, and hence only two points ,ur and ~2 where p&x) is discontinuous and where consequently surfaces form, as shown in section 1. The results of this paper apply also to lattice gases, the proofs requiring only slight modifications. For the lattice gas with

and

dr) I 0

for r # 0,

it is known4) that there is only a gas phase and a liquid phase, i.e. po(,u) has only one discontinuity pi =+x:(r) r

-

kTlogA”

for

v 2 2,

(4.1)

so that the positions of the surfaces are known in this case. The proofs also apply with little modification if v(r) is a positive potential, i.e. if 0 < q(r) g C Irl+-‘, where C and E are positive constants. It should be possible to extend our results to more general potentials, using the method of Fisher and Lebowit2 13). If the potential y(y) has the constant value ,u - PI for y in some region R, it follows that the region R is excluded from (3.25); i.e. we have an interphase region rather than a surface. Our result tells us nothing about what happens in R, but one might expect to get a uniform mixture of the two phases here. The proof of (3.25) works equally well if we replace either of the limit suprema in the definition of A(Y) by limit infima. Thus there are four equivalent definitions of E(Y). Our results also apply to the van der Waals system79 11). In this case we

HYDROSTATICS

AND

THE

SEPARATION

OF PHASES

277

take Y(r) = 4(r) + ?J’“KVr)* Then pe depends on y’ and its van der Waals limit PO@&0 +)

-;,;lin

PO(cL> Y’),

(4.2)

exists for all but a countable number of values of ,u and for a wide choice of the functions q and K (see Gates and Penroselr). Since we want y to vary slowly compared to 9 (i.e. y Q y’) we must take y’ -+ 0 after y + 0. Hence from (3.25) we have

for all but a countable number of values of p - v(y). This is a generalisation of the result of Kac and Thompsoni). Again (4.3) holds even for those functions K which yield crystalline states. If we choose y’ = y above, the result is quite different. The thermodynamic functions for this case were considered by Gates and Penroseli). The result one expects to get is

-Y')I,

(4.4)

P~[~(Y)I + 14~1+ !dy'fi(y') K(Y-Y') = PU,

(4.5)

3(y) = PO[~ -

Y(Y) -

S dy’fi(y’)

K(Y

or equivalently

where pa and ,ua are derived from the interaction potential q(r). The latter equation is the same as an equation of van Kampenrs) and can be derived (but not rigorously) from the variational principle of Gates and Penrose using functional differentiation. As indicated by van Kampen, this equation describes the actual structures of the gas-liquid interface, unlike (4.3). This is because here the density +i is on the same distance scale as the potential y’K(yr). Therefore, the rigorous derivation and examination of van Kampen’s equation seems an important unsolved problem. The present work can be extended by expanding nr(y/y, y) as a power series in y. The first term is just pa,$ - y(y)], while the- other terms represent corrections to the density due to the fact that the field does not vary infinitely slowly. This expansion provides an alternative to that obtained by Lebowitz and Percuss). A graphical prescription for the general term of the expansion can be given, and this will be presented in a future paper. Acknowledgement. I am grateful to 0. Penrose and E. R. Smith for helpful discussions, and’ to “The Royal Commission for the Exhibition of 185 1” for financial support.

D. J. GATES

278

APPENDIX The Bozlnds W+ and W-. We recall from (3.14) that W+(o, w’) is an upper bound, for all n, N and J2 on the energy of interaction, W(w, 0’) say, due to q(r), between o’ and D - CO,where o’ is a cube, concentric with but

smaller than o. We surrounded

w by shells of thickness

(s + 2~) as shown

in fig. 2. Let CO’have side s and w have side (s + 221), so that 10’/ = sy

/co1 = (s + 2U)“.

and

(A.1)

From (1.12) it follows that w(w, 0’) < CN@‘)

g N&d;“-‘,

(A.2)

k=l

where N(w’) is the number of particles in o.I’, Nk is the number of particles in the Kth shell around o (see fig. 2) and dk is the least distance from LO’to the Kth shell. We have from the paper of Penrose 14) that N(w’) I pc(s + 2~0)~ < pc(s + 2~ + 2ro)‘,

(A.3)

where YOis the hard-core diameter and pc is the close-packing density. Also, since the Kth shell contains (2K + 1)Y - (2K - 1)” cubes of volume w, it follows from ref. 14 that Nk I

[(2h +

1)” -

(2k -

I)‘] pc(s + 2% + 2ro)‘.

(A.4)

We also have dk = u + (k -

l)(S +

2U) 2 kU.

Second

s+ 2u

First

(A-5)

shell

shell

t-

s

43b W-W’

s+2u

w’

i;-;

Fig.12

/lo

HYDROSTATICS

AND THE SEPARATION

279

OF PHASES

With these results, (A.2) becomes W(m, w’) I

(A.6)

W+(m, w’),

where W+(c0, 0’) = C&s

+ 2% + 2re)%Y”J,,

where, in turn, Jy +-Y-E[(2k

+

< 2. 3”-l v ;

l)y -

k-l-”

(2k -

l)“]

< 00.

P-7)

k=l

Now, following

Lebowitz

v/(v + E) < ?j <

and Penrose’)

let us choose u = sn where

1.

Then we have

W+(w 4

lim

,w,

M--toJ

=;z

C&V

(s + 2ss + 2re)2’ S-V(“+e) (s + 2sv)*

= CpgJv lim s~-V@+~) = 0, S-+=0

(A-8)

which is the result quoted in (3.22). The second quantity we need to calculate is the lower bound W-(o) for all n, N and Q on the energy of interaction, W(U) say due to F(Y) between o and Q - o. We introduce the subcube m’ in w as before and write W(m) = W(fB, 0’) + W(c0, o’),

(A.9)

where W(c0, w’) is as before, and I$‘(u, LU’)is the energy of interaction between the “corridor” w - w’ and Sz - w. The argument used above for W(w, w’) also gives W(w, 0’) 2 -Cp;(s

(A.lO)

+ 2u + 2ro)%.~“-~J+

Also, by the argument leading to eqs. (4.9) and (4.11) of ref. 7 we have W(w, w’) 2 -pc[(s > -

+ 2U + 2re)D -

(s -

2re)*] @’

pcv@‘(2zL + 4ro)(s + 221 + 2ro)“-l,

(A.1 1)

where -@’ is a lower bound on the interaction of a given particle with all its neighbours. From (A.9), (A. 10) and (A. 11) we find that w(o) 2 W-(o) where W-(w) -

= -C&s vp@(2u

+ 221 + 2ro)2VzY-8 Jv + 4ro)(s + 224 + 2ra)‘-!

(A. 12)

280

HYDROSTATICS

AND

THE

SEPARATION

OF PHASES

Choosing 21= sn as before gives V(YfE)

-

p&D’2ss-1]

= 0,

(A. 13)

which is the result quoted in (3.21).

REFERENCES

1) Kac, M. and Thompson, C. J., D.K.N.V.S. Forhandlinger 42 (1969) 63 (eq. 35). 2) Lebowitz, J. L. and Percus, J. K., J. math. Phys. 4 (1963) 116. 3) 4) 5) 6) 7) 8) 9) 10) 11) 12) 13) 14)

Griffiths, R. B., J. math. Phys. 8 (1967) 484. Lebowitz, J. L. and Penrose, O., Commun. math. Phys. 11 (1968) 99. Fisher, M. E., J. math. Phys. 6 (1965) 1643. Ruelle, D., Commun. math. Phys. 18 (1970) 127. Lebowitz, J. L. and Penrose, O., J. math. Phys. 7 (1966) 98. Penrose, O., Statistical Mechanics, Foundations and Applications, IUPAP meeting, Copenhagen 1966, Editor T. Bak, Benjamin (1967) p. 101. Fisher, M. E., Arch. rat. mech. Anal. 17 (1964) 377. Bartle, R. G., Elements of Real Analysis, John Wiley and Sons Inc., (New York and London) p. 224, problem 19(f). Gates, D. J. and Penrose, O., Commun. math. Phys. 15 (1969) 255 and 17 (1970) 194. Van Kampen, N. G., Phys. Rev. 135 (1964) A362. Fisher, M. E. and Lebowitz, J. L., Commun. math. Phys. 19 (1970) 251. Penrose, O., Phys. Letters 11 (1964) 224.