Analogies between Analogies

Preferred Citation: Ulam, S. M. Analogies Between Analogies: The Mathematical Reports of S.M. Ulam and his Los Alamos Collaborators. Berkeley: University of California Press, c1990 1990. http://ark.cdlib.org/ark:/13030/ft9g50091s/

Analogies between Analogies

The Mathematical Reports of S.M. Ulam and His Los Alamos Collaborators

S.M. Ulam

UNIVERSITY OF CALIFORNIA PRESS

Berkeley · Los Angeles · Oxford

© 1990 The Regents of the University of California

Foreword

"Good mathematicians see analogies between theorems or theories. the very best ones see analogies between analogies." Stefan Banach

Slanislaw Ulam's affiliation with Los Alamos National Laboratory spanned over two-thirds of his professional life. There was no aspect of its mathematical activity during this period in which he was not involved. either centrally or tangentially.

His catholic view of the role of mathematics vis-à-vis other sciences extended far beyond into the mathematical and scientific community at large, as did his genius for problem formulation and for applying the most abstract ideas from the foundations of mathematics to computing. physics. and biology. In addition he possessed the ability to excite others many of them not trained as mathematicians -and involve them in his researches. The impact of his work is still felt both at the Laboratory and in those larger communities, for he liked to disseminate his ideas orally in an ever wildening round of lectures and seminars from where they took on a life of their own. The Monte Carlo method which he originated with von Neumann in order to study neutron scattering and other nuclear problems at Los Alamos is one such example. Its offshoots are now so universal that they are even applied to regulate traffic lights!

His influence, along with that of John von Neumann. the brilliant lHungarian nlathemrlatlician. contributed to the establislhment at the Laboratory of an atmosphere and a tradition that fostered and supported an exceptional- if not unique--interaction between mathematics and science. Extensive testimony and documentation concerning the integral role that Staln Uilam played in this interaction can be found in "From Cardinals to (haos. Reflexions on the Life and Legacy of Stanislaw Ulam" published by (ambridge University Press in 1989.

From 1944 until his death in 1984. while connected with Los Alamos in a variety of ways. from staff member, to group leader. to research advisor. to 'no-fee consultant'-one of his favorite expressions-he wrote Laboratory reports (mnany with the help of trusted and talented collaborators) that show a breadth of scientific interests unusual for a mathematician. They cover pioneering work. in the horse-and-buggy days

― x ―

of computing. on mathematical modeling of physical processes, nuclear rocketry, space travel. and biomathematics. (Another eleven. weaponsrelated reports, with Evans. Everett. Fermi. Metropolis, von Neumann. Richtmyer. Teller. Tuck. and others are still classified and unavailable for publication.)

Mathematically speaking, three motifs run through Ulam's theoretical and applied work: the iteration or composition of functions. or relations; the use of evolving computer capability in the exploration of analytically intractable problems: and the introduction of probabilistic approaches-while knowing that most practical applications are made in the presence of uncertainty. The fusion of these themes is characteristic of the central contributions of this collection.

As to the quotation from which the title of this book is derived, one must remember that Ulam held Banach. along with von Neumann and Fermi. "as one of the three great men whose intellects impressed me the most." Stefan Banach was an outstandingly original Polish mathematician and one of the founders of the now famous Lw6w school of mathematics (see Chapter 20. Preface to "The Scottish Book.") Banach was Ulam's friend and mentor in Poland before World War II. His influence on Ulam was profound and Ulam liked to quote his comment on the ability of some mathematicians to see "analogies between analogies." There is no question that in the practice of his craft and his art. Ulam was guided by this principle, and that he. in turn. epitomized its application. In addition, to Ulam the idea of analogy was itself amenable to mathematical discussion.

In 1983. when D. Sharp and M. Simmons. editors of this Los Alamos Science series, asked Ulam to gather his unclassified-and declassified-reports for publication in one volume, they intended to omit a few that had appeared elsewhere. After his death in 1984. it was decided to publish them all as many represent preliminary studies of subjects that were subsequently expanded elsewhere, leading, in several instances. to the development of new and extensive theories.

Ulam had dictated brief introductory notes and a sketch for a preface which he intended to develop. Rather than put words in his mouth. it was thought more appropriate to reproduce his notes in their short and unpolished form. It was also decided to leave the style and substance of the reports untouched. as evidence that scientific advances do not usually arise in their final, definite form. More often than not they are the product of sequences of tentative. sometimes repetitive. and even at times inaccurate steps. Two appendices complete the volunie: a list of Ulam's publications and a brief biographical chronology. (More detailed biographical material can be found in his autobiography. "Adventures of a Mathematician." as well as in "From Cardinals to Chaos.")

― xi ―

This collection represents an important complement to the selection of papers and problems, mostly in pure mathematics, published by MIT Press in 1974 as "Sets, Numbers, Universes," and to a volume of essays, "Science, Computers, and People," published by Birkhiuser in 1985. And, whereas these two books are composed of papers readily available albeit scattered in the scientific literature, the Los Alamos Reports have been for the most part difficult of access and little known. Their historic value is therefore very real.

As mentioned earlier, many of these reports and much of Ulam's work was done in collaboration. He liked to stress the importance of working with collaborators, with whom, he said, the nature of their "shared ideas and techniques" depended "on the personality and experience of the individuals." We add a few words about these colleagues as well as about the work they and Ulam engaged in.

As early as 1944, when he joined the Manhattan Project, Ulam and David Hawkins (Chapter 1) were "playing" with the notion of branching processes, or multiplicative systems, as they called them, motivated by their application to atomic physics.

David Hawkins, a philosopher of science by profession and mathematician by inclination, was, in Ulam's words "the best amateur mathematician I know," and they became fast friends. Hawkins is presently professor emeritus at the University of Colorado.

The work in the Ulam-Hawkins report was subsequently developed with Everett in the three extensive reports grouped in Chapter 3, which are reproduced here for the first time. While clearly motivated by the need to understand neutron multiplication in fission processes, the reports lay the foundations for-in their own words-a "formalism general enough to include as special cases the multiplication of bacteria, radioactive decay, cosmic ray showers, diffusion theory and the theory of trajectories in mechanical systems."

C. J. Everett, who died in 1987, was a mathematician with whom Ulam worked on a conceptual as well as technical level in Wisconsin and at Los Alamos, where he became a member of Ulam's group. An eccentric, shy, and witty man, he was quite probably the only person who ever opted for bus transportation to come to Los Alamos for a hiring interview, and he was known for having turned in a monthly progress report-in which staff members were supposed to describe their research-which said tersely "progress was made on last month's progress report."

The first written proposals for the Monte Carlo method put together in a 1946 "report" called "Statistical Methods in Neutron Diffusion" appear in Chapter 2. This method of approaching precise but intractable problems through the introduction of random processes and

― xii ―

probabilistic experimentation, has found wide application not only in areas close to those motivating its origin but others more removed, such as operations research, and combinatorics. In fact, the "report" --of which only eight copies were made-consists of two letters and handwritten calculations photographed and stapled together. Its cover specifies that the "work" was "done" by Ulam and von Neumann and "written" by von Neumann and Robert Richtmyer-then head of the Laboratory's Theoretical Division. Its informality attests to the casual manner in which information was disseminated through the Laboratory at the time.

The long term professional and personal rapports between von Neumann and Ulam need not be recounted here--references to them can be found in the books already mentioned. Suffice it to say that though there exist few papers and abstracts under their joint names, von Neumann's extensive correspondence with Ulam attests to their interacting interests in pure mathematics, in pioneering computer technology and techniques, and in cellular automata and the brain. (The correspondence is now stored in the archives of the Philosophical Society in Philadelphia.)

Ulam's collaboration with Enrico Fermi initiated the computer simulations of nonlinear dynamical systems that lead to the evolution of a major field of research popularly labeled "chaos theory." Fermi called this work, which was developed with the programming assistance of John Pasta, "a minor discovery," a modest understatement given the seminal character of this investigation (Chapter 5.) Chapters 10 and 11, with P. R. Stein, address the subject of nonlinear transformations in greater detail.

Fermi, with whom Ulam became acquainted in Los Alamos during the war, was a man of simple tastes and life style. The Ulams had an opportunity to sample this while motoring together across France one summer. Feeling ill at ease during a lunch in a recommended temple of gastronomy, Fermi decreed he would select the night's lodgings. Meandering through a picturesque valley he chose a modest inn by a babbling brook where, after dinner, sitting under the stars they discussed physics and new mathematical problems to experiment with after the vibrating string calculations. However the night's encounter with fleas, bedbugs, and mosquitoes made him admit the next morning that the higher-class hostelry next door that Ulam had eyed, might perhaps have provided a more restful night.

In the area of space technology, Ulam investigated schemes for nuclear rocketry with Everett and with Conrad Longmire, a physicist from the mountains of Tennessee who played a mean banjo. Chapter 7 describes a way to propel very large space vehicles by a series of small

― xiii ―

Foreword external nuclear explosions which later developed into Project Orion. Chapter 9 deals with the propulsion of space vehicles by extraction of gravitational energy from planets. Schemes based on a similar idea are now used in "flyby" missions to the outer planets and to provide part of the energy for spacecrafts going beyond the planets.

A study of patterns of growth, with Robert Schrandt (Chapter 12) investigates how simple recursively defined codes can give rise to complex objects. Such studies have become a growth industry of their own in the improved computer graphics world of today.

Several other reports are devoted to biomathematical questions. Their findings have opened new fields of biomathematical research. Abstract schemata of mathematics are applied to pattern recognition with the help of computers investigating, for example, the way visual pictures are recognized. Using metrics in molecular biology shows how a new mathematical concept of distance between finite sequences or objects can be applied to reconstruct the evolutionary history of biological organisms.

Closely involved with Ulam in this work was Paul Stein, a physicist turned mathematician under Ulam's influence who became an invaluable collaborator able to implement and develop the gist of Ulam's directions. William Beyer, a gifted fellow mathematician, also collaborated on a conceptual and technical level in the biological investigations.

John Pasta, Mary Tsingou-Menzel, Robert Schrandt, and Myron Stein lent their talents to creative programming, at a time when the art was in its infancy and pre-microchip-era machines with names like "Eniac," "Maniac," "Johnniac," presented storage and timing constraints. Pasta, who died in 1980, was a self-made man of Italian descent. He had furthered his education and became a physicist while working on the New York city police force.

The mathematician Al Bednarek, one of the editors of this volume and coauthor of this foreword, also collaborated with Ulam on problems of parallel computation (Chapter 18). He is a former chairman of the mathematics department at the University of Florida.

Shortly after his arrival at Los Alamos in 1944, Ulam was asked by a colleague what it was that he was doing. Since at the time he was a very pure mathematician and had not yet familiarized himself with the nature of the work, his Socratic answer was "I supply the necessary don't know how!" Stan Ulam's "necessary don't know how" as well as his modestly unenunciated "know how" are sorely missed by all who were privileged to have known him or worked with him.

Last but not least, the editors wish especially to thank Peggy Atencio, Ben Atencio, Janet Holmes, Debi Erpenbeck, Gary Benson, Chuck

― xiv ―

Calef and Gloria Sharp, among other members of Los Alamos Information Services Division, and above all Chris West and Pat Byrnes, for their herculean efforts in transposing into print and formatting these extremely difficult reports, and also Patricia Metropolis for her invaluable informal advice and help. The editors assume full responsibility for any existing discrepancies or inaccuracies. They also gratefully acknowledge the permission granted by Rozprawy Matematyczne to reprint their edition of the report "Non linear transformations studies on electronic computers." The preparation of this book was done under the generous auspices of the Los Alamos National Laboratory.

A. R. Bednarek, Gainesville, Florida Francoise Ulam, Santa Fe, New Mexico

― xv ―

Preface

The collection of these reports, which appeared over the considerable span of years that I spent at Los Alamos, concerns a great variety of topics. Its very heterogeneous nature illustrates the diversity of the programs and of the areas of research that interested the laboratory.

Before World War II it was almost exclusively in the universities, in the graduate schools of the larger institutions that scientific research used to take place. The Bureau of Standards and a very few large industrial companies such as Bell Telephone, General Electric, and some pharmaceutical firms were the exception to the rule.

This little collection may bear witness, in a very modest way, to the wide-ranging changes, which are still going on in the organization and practice of research in this country and abroad. Because of the novel problems which confronted its scientists during the wartime establishment of Los Alamos, the need arose for research and ideas in domains contiguous to its central purpose. This trend continues unabated to the present.

Problems of a complexity surpassing anything that had ever existed in technology rendered imperative the development of electronic computing machines and the invention of new theoretical computing methods. There, consultants like von Neumann played an important role in helping enlarge the horizon of the innovations, which required the most abstract ideas derived from the foundations of mathematics as well as from theoretical physics. They were and still are invested in new, fruitful ways.

An enormous number of technological and theoretical innovations were initiated at this laboratory during these forty years. To mention but a few, besides the advances in computing, one can name research on nuclear propulsion of rockets and space vehicles, in molecular biology, and on the technology of separating cells.

The growing importance of research laboratories such as this one became a not exclusively American phenomenon. For instance the aspect of academic research has changed almost beyond recognition in France. What used to be, before World War II, almost exclusively the province of universities, has now shifted to the French National Center of Research (Centre National de Recherche.)

― xvi ―

This collection of Los Alamos Reports ranges over almost four decades and may illustrate, I hope, how a mathematical turn of mind, a mathematical habit of thinking, a way of looking at problems in different subfields of physics, astronomy, or biology can suggest general insights and not just offer the mere use of techniques. Ideas derived from even very pure mathematical fields can provide more than mere "service work," they may help provide true conceptual contributions from the very beginning.

The period in question has seen the origin and development of the art of computing on a scale which vastly surpasses the breadth and depth of the numerical work of the past. In at least two different and separate ways the availability of computing machines has enlarged the scope of mathematical research. It has enabled us to attempt to gather, through heuristic experiments, impressions of the morphological nature of various mathematical concepts such as the behavior of solutions of certain nonlinear transformations, the properties of some combinatorial systems, and some topological curiosities of seemingly general behaviors. It has also enabled us to throw light on the behavior of solutions of many problems concerning complicated systems, by allowing numrerical computations of very elaborate special physical problems, using both Monte Carlo type experiments and extensive but "intelligently chosen" brute force approaches, in hydrodynamics for example.

A number of such experiments have revealed, surprisingly, a nonclassical ergodic behavior of several dynamical systems. They have showed unexpected regularities in certain flows of dynamical systems, in the mechanics of many-body problems, and in continuum mechanics. Recently they have been applied to the study of elementary particle physics set-ups and interactions.

And now there appear some most exciting vistas in the applications of mathematics to biology that deal with both the construction and the evolution of living systems, including problems of the codes, which seem to define the basic properties of organisms and ultimately may provide us with a partial understanding of the working and evolution of the nervous system and some of the powers of the brain itself.

In addition these reports show the varying involvments of my collaborators and myself. I particularly want to stress the importance of

― xvii ―

the role of collaborators. An ever increasing number of publications of mathematical research is proof of the advantages derived when two or more authors share ideas and techniques. The nature of this exchange varies from case to case, depending on the personality and experience of the individuals.

These few very sketchy remarks are merely intended to emphasize how the necessity of defense work at the frontiers of science has continued to this day to stinmlate research in a multitude of directions.

S. M. Ulan Santa Fe, February 1984

― xviii ―

Faded cover of the original 1944 report after it was released from classification in 1956. Note its low number.

― 1 ―

1—
Theory of Multiplicative Processes:
With David Hawkins (LA-171, November 14, 1944)

This report treats branching processes, of neutron proliferation for instance, through a mathematical theory involving compositions of generating functions.

It is a precursor of the studies of multiplicative systems in several variables written with C. J. Everett in 1948 (LA-683, -690 and -707) which develop an elaborate and basic theory of "multiplicative" (branching) processes. See for example a book by T. E. Harris: The Theory of Branching Processes, published by Springer in 1963. (Author's note.)

Abstract

Generalpropertiesofstatistics ofmultiplicativesystemsare discussedtogetherwiththestudyoffluctu ations in the number of particlesin such systems. A general methodisindicatedthroughwhichonemay study the fluctuations in the casewhereonetakesinto accountthefactorsofgeometryand time-dependence ofconstants.

The statistical theory of multiplicative chain processes does not compare in completeness to date with the corresponding theory of additive processes. The present paper is intended primarily as an exposition of a simple theory of the statistics of multiplication, permitting application to a variety of special problems.

― 2 ―

Analogies between Analogies The simplest (the "Bernoullian") case may be described as follows: A particle can produce, with probabilities po, pi, P2. , Pn,.. . a number 0, 1, 2, 3, ..., n,... of similar particles in one generation. We assume that each particle produced has again the same probabilities of producing n offspring. We also assume that each particle dies at procreation. Required is the probability law pk(n) for any generation k.

We remark parenthetically that this formulation makes the multiplicative process essentially discrete and finite. The statistics of neutron multiplication involves a continuous process as well, namely a random distribution in energy, space, and time. We disregard this aspect initially. Later we shall show that the admission of such continuity leads to a generalization of the methods described below. There are, in the meantime, two physically accurate interpretations of a discrete series: (1) one can represent the chain process as a graph; the n particles in the kth generation are the n lines connecting the kth and the k + 1st - branch points in a chain or set of chains; (2) the n particles are those in existence at the kth unit of time, where the probability law pl(n) is the distribution one unit of time after the introduction of a single particle. If a time unit be chosen equal to the average time between fissions, the distinction is in many cases not crucial. Frankel, and later Feynman, studied the continuous process. We shall show later that their differential equations of the random process correspond to the infinitesimal transformations of the group in which our iteration (see Theorem I) may be imbedded.

1. The first problem to consider is this: we are given an amount and arrangement of active material. In this system a neutron produces on the average n neutrons with probability p(n); C=)p(n) = 1.p(O) is the average probability of leakage or absorption, without subsequent production of neutrons. p(n) normalized for n > 0 is a nuclear constant, so far purely empirical, known as to its first moment and less accurately as to its second mnoment. Required is the probability of having n neutrons after k generations (or units of time). This problem is solved, in principle, by:

Theorem I. Let f(x) be the generating function of the distribution of the lnmber of offspring, i.e., f(x) =n =op(n)xn. Then the generating function for the kth generation fk(x) = fk(x), the kth iterate of f(x). [The kth iterate is defined as follows: f¹ (x) = f(x), fk(x) = f(fk-l(x)). The theorem asserts that the probability pk(n) is given as the coefficient of xn in the ascending polynomial or power series expression of fk(x). The physical multiplication of

― 3 ―

the random variable is reflected in the iterated substitution by which x- f(x)].

Proof: Starting with one neutron in the 0th generation we obtain, with probability pk(n), n neutrons in the kth generation. Beginning with r neutrons, denote the corresponding probability by p() (n). Now assume that a chain is started by one neutron. We have Pk(n) = EPk-l(r)p) (n). r=O

Now if f(x) is the generating function of the distribution pl (n), the generating function of the distribution p(r(n) is [f(x)]r. This follows from the assumption that contemporary neutrons are independent in procreative powers, and from the theorem (of Laplace) that the generating function of a sum of independent random variables is the product of their generating functions. The above proposition may also be verified for r = 0, since p(l)(n) = 0 for all n > O. Substituting generating function for probability in the above equation, we have:

Two remarks may be made at this point. (a) The simple proof above sustains a more general theorem if the distribution generated by f(x) is not constant, but time- or generation-dependent. Instead of the iterate ffff...(f(x)), we will have some fgh...(q(x)). By the mode of argument established, the chain process may be analyzed one step further. Let g(y) = ay +b be the generating function for the probabilities b of loss or absorption of a single neutron and a of producing fission, with a + b = 1. Let h(x) = cl x + c₃x³ + ... be the generating function of distribution of neutrons per fission. Then if the two are combined by the transformation y -, h(x), we have that the distribution of neutrons per neutron is generated by f(x) = g[h(x)]. (b) If on the other hand we start from a single fission, and wish to know the distribution for the number of first-generation fissions, this is given by F(y) = h(g(y)). The iterates of f(x) and F(y) are connected by simple and evident relations.

There remains the practical problem of determining coefficients and other properties of fk(x), given f(x). To this end we first shall establish some general properties of iteration.

― 4 ―

2. Let f(x) be a monotone function. Assume, for example, f(x) increasing, i.e., if x< y, f(x) < f(y). A fixed point for f(x) is a value xo such that f(xo) = xo. The set of fixed points for a continuous function is closed, i.e., the points which are not fixed form a collection of disjoint intervals, whose endpoints are fixed points. If we form the sequence fk(x) for a given x we obtain a sequence of points converging to a fixed point xo which forms the endpoint of the interval in which x is situated. In fact, there are two cases possible, either f(x) < x or f(x) > x. From the monotone character of f(x) it follows that we shall have correspondingly either fk(x) < fk-l(x) or fk(x) > fk-1(x) for all k. Unless these sequences tend to -oo or +oo, they will have limit points. If now limk-, fk(x) = xo, we must have f(xo) = xo. In fact limk-, f(fk(xo)) = f(xo) = xo. In addition, it is easy to see that xo is the next fixed point to x (on the left or right depending on whether f(x) < x or f(x) > x. This follows from the fact that if f(x) is monotone and f(xo) = xo, f(xl) = xl, then for all x such that xo < < < x, we have f(xo) = xo< f(x) < f(x₁) = xi.

In our case f(x) is a power series with all coefficients non-negative, f(O)>0,f(l) = 1. This function is certainly monotone and increasing for all non-negative x. Let xo be the first (non-negative) fixed point, xo certainly exists, the set of fixed points being closed. From these conditions it follows that limko,fk(O) = x₀ . But if the variable in a generating function is set = 0, the value of the function is the probability that the random variable takes the value 0. Hence xo gives us the limit of the probability of mortality in the system. The probability of immortality is therefore simply 1 - x₀ , where xo is the smallest nonnegative root of the equation f(x) = x. It is easy to see that if, as in our case, all the coefficients in the expansion of f(x) are non-negative and f(1) = 1, then from f'(1) > 1 it follows that there is a root, and only one root x₀ , which is non-negative and < 1. If f'(1) < 1,xo = 1 is the smallest positive root. We obtain immediately therefore the familiar fact that neutrons in a subcritical gadget without source will, with probability 1, die out in a finite time. For the supercritical gadget the probability of indefinite production can be obtained by solving the equation f(x) = x.

The kth iterate of a function can be obtained by a simple graphical or mechanical method which is based on the fact that along the diagonal, f(x) = x. Thus we may for given x replace this argument by f(x), getting f² (x) graphically, then repeating f⁴ (x) and so forth. In the case of the generating function under discussion this shows that fk(x) very rapidly approaches its asymptotic form: for the critical or subcritical case the asymptote in the interval 0 <x < 1 is limkoo fk(x)-xl; for the supercritical case in the interval 0 <x < I the asymptote is

― 5 ―

limk,, fk(x) _ xo. This implies that for all positive powers of x in fk(x) the coefficients approach 0 uniformly, i.e., the mass of probability is either absorbed altogether into the zero region (subcritical case), or is spread out in an infinitely long tail (supercritical case). In the region of criticality the distribution has an infinitely long tail with mass approaching zero as the probability of mortality approaches one.

3. One of the important properties of generating functions is that they permit the calculation of moments. Thus if p, is the distribution itself, f(x) = '°_=opnXⁿ its generating function, we have, because obviously f(1) = 1, the first moment or expected value of the random variable = _=onpn = f'(1) = the first derivative of f(x) at x = 1. Similarly the second moment of the number of neutrons can be found if we know the second derivative. In fact

Similarly the rth moment can be found easily from the values of the first r derivatives of f(x) at x = 1. (The rth derivative at x = 1 is sometimes called the rth combinatorial moment.)

Our generating function is the kth iterate fk(x). It turns out that its first m derivatives depend only on the first m derivatives off(x) itself in a rather simple way. We have, in fact:

Theorem II. (a) [fk(x)]= = f'(1) = n (i.e., the proof of the intuitively obvious result that the expected number of neutrons after n generations is n.) (b) [fk(x)]=" =f"(x)=l₁ . [(f(l))k + (yf(l))k+1 + ..(f (1))² k- ]).

The proof is immediate by induction: [f(fk-l(x))] = f'(fk-l()))[fk-l()]'.

But for x = 1, fk-l(x) = x = 1; therefore since by assumption [fk-l(x)] =₁= [f'(l)]k-l we obtain our formula (a). By differentiating twice we obtain (b). Somewhat more complicated formulae hold for higher derivatives:

Their derivation is through recursive relations as follows: by differentiating the identical equation fk() = f(fk-l(x)) repeatedly, and in all places substituting xo for fk-l(xo), we obtain a

― 6 ―

sequence of linear first-order difference equations. Representing dr/dx^' fk (x) = Mk,r(Ml,r = ir) we obtain Mk,l =M lMk-l,1 Mk,2 =l M₂M2₁+ M₁Alk-l,2 Mk,3 = M3M₁3_ ,₁ + 3M2Ik_₁,l Mk-1,2 + MlMk-l,₃ each is of the form Xk = Ak-1 + Mlxk-1 whose general solution is k Xk = E M As-₁ + M=¹ x1 s=2 Solutions for the first three derivatives are

Since in the function under discussion xo = 1 is a fixed point, these derivatives are the combinatorial moments of the distribution. We may now consider the three cases where M₁(== v of common use) is > 1, < 1, or = 1.

a) In the supercritical case where M1> 1, it is clear from the method of deriving these factorial moments that if the random variable n is measured in units of M₁, all moments approach a finite asymptotic form. Computation of moments for this asymptotic distribution may be greatly simplified as follows: Let us define a function O(X) = xl/Ml, the kth iterate being _k(x) = x(1/M^l) k. The generating function fk(Qk(x)) if expanded in powers of x(1/MI)k has the

― 7 ―

same coefficients as fk(x) but these are now probabilities associated with the number of particles measured as fractions of the expected number. This is to say that the distribution is scaled in units of Mk = -Uk, and its first moment = 1. Since for the supercritical case all moments approach a constant value as k -- oo when scaled in this way, and since the generating function is monotonic in the region (0, oc), there exists a common limiting value, g(x) of both f(k[k(x)] and fk-l[k-¹ (x)]. Since fk[ck(x)] = f[Ok-1[jk(O(x))]], we may write in the limit: g(x) = f[g[((x)]], Q(x) = xl/Ml, f(x) given, and from this functional equation for g, its moments may be obtained from the second, third, etc., derivatives of g by solving only linear algebraic equations.

b) In the exactly critical case, M₁= 1, the moments are Mk,1 = 1 Mk,2 = k • M2 Mk,3 = kM₃+ 3 2

This is a distribution in which Pk,o - 1 - 1/kM₂, and such that if the system has not died in the kth generation, the expected number of neutrons is _kM₂.

c) In the subcritical case all moments converge to zero, but are approximately proportional to the first moment.

4. We may consider here briefly a simple special case, in which the iteration problem may be solved exactly.

Let f(x) = (ax+b)/(cx+d); we have here a three-parameter family of functions (one of the four constants a,b, c, d, is immaterial). We can adjust them so that f(1) = 1, and f'(1) = 7. We can then impose another condition, either on f"(l), or so that f(xo) = xo, where xo is the "true" probability of mortality. Functions of the above sort form a group under substitution. This can be verified directly by substituting. (They form the so-called projective group of the line.) A fortiori the iterated function fk(x₎ = (akX + bk)(ckx + dk)

By expanding fk(x) in a power series in x, we obtain the exact solution of our problem in this fairly general case. We determine the constants by the following three relations:

― 8 ―

(1) Because f(1) = 1, we have for every k fk (1) = 1 which gives ak + bk = Ck + dk.

(2) Similarly, for the second fixed point xo of f(x), i.e., the root xo= 1 of f(x) = x, we have f(xo) = xo, and therefore for all k fk(xo) = xo or akxo + bk = Ckxo + dk and xo = -b/c, from axo + b = xo(cxo + d).

(3) From the results of section 3, we know that [fk(x)]l=₁= [f,(1)]k =-k This gives a(cx + d) - (ax + b) (cx + d)² or taking account of (1) (a-c)(c + d) and therefore for all k(ak-ck)k (Ck - dk)

From the above three relations it is easy to calculate the constants ak, bk, Ck, dk in terms of v and one arbitrary parameter. By eliminating ak, bk and developing into a power series, we get, noting that ck/dk /₁Vk - 1, assuming, e.g., v > 1, the result in the form fk(x) = [(ax + b)/vk+l]{l + (1- l/vk)x + (1 - l/vk)² x2+ ... (1 - 1/vk)nxn +...

This constitutes a complete solution of our problem. It is interesting to note that the probability of having n neutrons decreases geometrically with n; the ratio of the successive terms is in the case v> 1, k large extremely close to 1. The distribution has the form of an exponential, decreasing very slowly. Asymptotically the probability of having exactly n neutrons is independent of n. This result shows also the possibility of enormous fluctuations in multiplicative systems. The "law of large numbers" in its ordinary formulation is not true for multiplicative processes. In fact the probability of having more (or less) than f times the expected value of neutrons tends to a positive constant (dependent on £). The following form of the law of large numbers is valid, as the examination of the distribution shows at once:

― 9 ―

Theorem III. Given an e > O, there exists an N such that for all k > N, the probability of the number n of neutrons in the kth generation being such that (_ e)k <n < (+-te)k is greater than 1 -e: P{(v-_)k < n < ( + E)k⁾ > 1-.

It remains to discuss the most general form of the distribution. We hope to do this later through two methods, one consisting of the consideration of functions of the form hfh^- (x), where f is of the projective linear form discussed above, and h(x) is an arbitrary monotonic function. The kth iterate then is simply hfkh-l(x). The function h(x) will give us more arbitrary parameters for our real distribution. The second method consists in developing f(x) into a series of functions whose terms have the "projective" form.

Finally it may be remarked that the limiting distribution obtained above is formally identical to those obtained by Frankell and Feynman who used a continuous time parameter instead of our discretegenerations model. Their physical model is somewhat different and leads to the finding of the infinitesimal transformation of the continuous, abelian, one-parameter group into which the group of iterates of a function can be imbedded.

5. There are many other problems besides the question of the probable number of neutrons after k generations which can be solved by operational methods. The first we shall consider is that of a subcritical system (i7 < 1) with a source. We suppose that the distribution of neutrons entering the system in a given generation has the generating function X(x),f(x) being the generating function of the system itself as before, we shall have

Theorem IV. The generating functions in the zero, first, second generations are the functions: O(X), (x). * [f(x)], O(X) - ) [f(x)]- [f2(x)]. Proof is completely analogous to that of Theorem I. In general, letting Fk(x) represent the distribution in the kth generation a) Fk(x) = ¢(x) . Fk-l[f(x)] If the system is subcritical, but sustained at a definite level by the source, we shall have the limiting distribution-or its limiting generating function-as a nonsingular function of x : limk-, Fk(x) = F(x), F(1) = 1. Passing to the limit on both sides of our equation a) we get

― 10 ―

b) F(x) = F(x) [ F[f(x)] where O(x), f(x) are given. One has to determine F(x) from this functional equation. Even without doing it one can obtain at once useful statistical information, for example the moments of F(x), by differentiating b). Thus: F'(1) = (1- f') F"(1) =_ + 20f' + f"+/' (1 - f,² ) (1 - f'² )(1 - f')

giving us a way to compute standard deviations, and similarly, more complicated expressions for the higher derivatives and moments. The first derivative--the expected value--being inversely proportional to the degree of subcriticality becomes infinite if f'(1) approaches 1.

6. We come now to the probability distribution of the sum of all neutrons in the system from the first to the kth generation. We have established previously that if f(x) = E'=o0 pnXn is the generating function for the probabilities of n particles in the first generation then the generating function of the kth generation is given by the kth iterate fk(x).

If we want the generating function for probabilities of having the total of n particles from the first to the kth generation, we shall proceed as follows.

The total of n particles can be obtained by any one of the following mutually exclusive cases: we can have I in the first generation and n-1 in the remaining k - 1, or 2 in the first generation and n - 2 in the remaining k - 1; in general we can have r in the first and n - r in the remaining k - 1 generations. The required probability is therefore the sum of q(n) =-Pr-P (n) . r

Here pn-r(n) denotes the pro the probability that, starting from r in the first generation, we shall attain from these r a total of n - r in k - 1 generations. But the r particles are independent of each other. The probability of getting the total of n- r from them is therefore the probability of n - r in the sum of these r variables. The generating function for the sum of the independent variables is the product of the generating functions corresponding to each of them. In our case it is the rth power of f(x). We are looking for the coefficient of x^n- r in [fk-1(x)]r. Our required probability qk equals therefore the sum with

― 11 ―

respect to r of coefficients of xn-r in [fk-l(x)]r, or the sum of the coefficients of xⁿ in r PnX [f k-(x)].l

But the coefficient of x" in ZrPnXr[fk-l(X)]r is the same as this coefficient in f(x fk-l(x)). This is true for all n. Therefore the generating function for q, is f(xfk-l(x)). Since n here is arbitrary we get:

Theorem V. The generating function for the time sum is: uk(x) = f[xukl(x)].

If we 'count" the original particle, this multiplies the generating function by x; expressing this slightly modified form recursively, we obtain the more convenient expression: uk(x) = xf[uk-(x)]

As we know we have, in general, a relation between moments of the nth order of a distribution function and the nth derivative of the generating function. We shall now show how one can compute the derivatives of uk(x) for any k in an explicit manner.

Since, as was shown above, uk(x) = xf[k-1(x)],

we may obtain the desired results by repeated differentiations, and by solving the resulting finite difference equations. But if k is allowed to approach infinity, and if the system is subcritical, li ux)limuk( k-i1(x)= u(x) k-—oo k-oo

Hence for the distribution of the total number produced, we have u(x) = x f[u(x)], differentiating, we obtain: u/'() = (1 - f (1)) u( -If" '(-f')] (1 - f')³

― 12 ―

These examples show how moments of the distributions can be computed for various problems in our discrete model. Otto Frisch has shown how, from the form of these moments, one can write their correct form for the continuous model, without having to solve the partial differential equations of the problem. This correspondence between the two models will be taken up later. It may be said that a generality of method has been established by the foregoing results, which demonstrate that the iteration of suitable operators corresponds to various physical observables connected with chain processes. For example it may be mentioned that the transformation x -(1/x)f(x) gives us the probability-distribution for differences between the number of neutrons in a generation and the number in the next generation. Thus fk-l[(1/x) .f(x)] generates probabilities of this kind. The mathematical description of a multiplicative chain process is seen to involve the iteration of a functional operator U. These operators U act on the domain of all monotone functions g(x), g(l) = 1. To summarize again just a few examples:

(1) U(g) = f(g), f here is a given monotone function, g represents any function of the domain on which U operates, i.e., g(x) is monotonic, g(1) = 1. This operator U is the only one that has been studied extensively in the literature. Its iteration leads to the simple iteration process: g(x), f(g(x)), f[f(g(x))]...fk(g(x)) ...

(2) U(g) = f(x-g), f a given function. The domain of the operator, i.e., the admissible g are the same, but there seems to be very little known about the iterates of this operator. This operator is tied to the probability law of the total number of particles produced.

(3) U(g) = ¢(x) g(f(x)); ((x), f(x) are given. The iterates of this operator give us the distribution of the number of particles produced when a source with given distribution 5(x) is acting constantly.

(4) U(g) = f[(/x) . g] . This operator relates to the probability distribution of the difference of the number of particles in successive generations. The study of conjugates, fixed points, etc., for such operators seems to be important. We hope to undertake this study later.

We turn now to a more complex version of the problem. Hitherto it has been assumed that the generating function was independent of temporal and geometrical factors. However, our methods are extensible beyond these limitations.

― 13 ―

7. The calculation of the probability distributions in the general case of heterogeneous particles will now be considered. So far we have assumed that the probability of generating n neutrons is the same independently of the parent neutron. If one takes the real situation where the system of the active material is of finite extent, then obviously the probability of leakage and absorption is a function of position of the parent nucleus. It is obvious that in general chemical or nuclear chain-reaction processes one has to deal with several kinds or even a continuous variety of the elementary generating functions.

In order to explain our methods of iteration of functional operators for this general case we shall take the simplest case of two kinds of particles. If we divide, for the first approximation, the sphere of the active material into two parts, an inner sphere and the outer shell, we shall characterize the neutrons generated in the one part by the subscript x, the others by subscript y. An x-particle can generate either x-particles again or penetrating to the outer shell y-particles, or, of course, leak out or be absorbed; the same, though with different probabilities, applies to the y-particles. In reality we should consider a one-dimensional variety of kinds of particles corresponding to all values of their distance r from the center of the sphere or even a twodimensional one if we want to take into account different velocities. To simplify the presentation we shall limit ourselves here to just two kinds (x and y).

We assume that the following elementary probabilities are given by the nuclear constants and by the integrals of the geometry involved. An x-particle can produce n(> 0) x-particles with the probabilities Pn and n(> 0) y-particles with probabilities qn. The probability of dying out absorption or leakage will be denoted by po. For the y-particles the corresponding probabilities will be denoted by Pn,n, and po. It is because of the geometry of the system that po and po are certainly different.

We now write the two functions of two variables each: f(x, y) = po + plx + . pn^xn + ... + qly + . qnyⁿ + . g(x, y) = po⁺ pl^x⁺...pn +.' + qlY + ... qnY + ...

The coefficients of f(x, y) give the probabilities of having in the first generation a given number of x- or y-particles starting with one xneutron. Those of g(x, y), if we start with a y-neutron.

Required are the probabilities of finding in the next generation a given number of x- and y-particles. Let us form the function f2(x, y) = f[f(x, y) g(x, )] .

― 14 ―

By reasoning exactly as in the proof of Theorem I (or Theorem III) we see that the probability of having n x-particles and m y-particles is given by the coefficient of xrly^M in f₂(x, y). If we started in 0th generation with a y-particle we will get these probabilities as the coefficients of xny^t in .g[f(x,y),g(x,y)]. By an obvious induction we obtain:

Theorem VI. The probabilities of having n x-particles and m y particles in the kth generation are given by the coefficient of xny^t in f[Tk-l(r)] or g[Tk-l(r)] (depending on whether we started from an x- or from a y-particle). Tk(p) is a transformation of the plane (x, y) into itself defined as follows: if p =(x, y) then T'(p) = T(p) [f(x, y), g(x, y)]; Tk(r) = T(Tk-l(p)]. Without going into the details of the proof or actual computations of moments we wish to conclude by the following remarks:

(1) In the case of 3 or any finite number r of different kinds of particles, the formalism necessary to obtain the generating function for the kth generation is the same. It consists of iterating a given set of r functions or a transformation in r dimensions (variables xl, x2... xr).

(2) One fairly general case where the coefficients of the mixed powers of the variables xa xa2 ... xar can be computed explicitly in a closed form for any number k of generations is when the given transformation is the r dimensional generalization of our projective transformations on the line, i.e., p = (xl,x2...xr); p' = T(r) = (x4x2...x4) where x: = fl(xl ...Xr) =(alIxi + . alrr + bl)/(cllx +...clrXr + dl) .................................................................. _r= fr(xi . ..r) = (ariXl + . . arrxr + br)/(CrlXl + .. Crr + dr).

(3) The computation of moments of the distribution in the most general case does not involve the explicit knowledge of Tk(r), but can be obtained through the knowledge of the moments of the r given elementary functions fl(Xl ...Xr) ...fr( ...Xr)

The role of the numerical multiplication of moments is here taken over by matrix multiplication.

(4) The other operators corresponding, e.g., to U(g) = f(x g), etc., have not been so far investigated in the r-dimensional case.

― 15 ―

Conclusions Regarding Applications

The expected value of the number of neutrons per fission "v" is known with fair accuracy. The critical mass and the expected number of neutrons in a gadget depend on this constant alone. Very little seems to be known, however, about the distribution function of the number of neutrons or even only about its second moment. The great fluctuations in multiplicative systems discussed above are of some practical interest for the following reasons:

1. The correct timing of the initiation of the gadget is vital for high efficiency. Even with good sources there will be an uncertainty of several generations time -due to fluctuations in multiplication.

2. The fluctuations of multiplication are of interest in all "integral" experiments.

3. For gadgets large in comparison with the mean free path for fission, the spatial fluctuations may destroy the initial spherical symmetry.

In dealing with such problems it is useful to develop a uniform technique for describing the statistics of multiplicative phenomena. This paper constitutes a first step consisting essentially in the observation that the iterated substitution (of a function, or more generally of a functional operation) represents exactly the statistical laws of multiplicative processes. In the sequel, it is hoped to apply this technique to the study of the problems of geometrical- and time-dependence of the process.

Reference

1. Stanley P. Frankel, The Statistics of the Hypercritical Gadget, LAMS-36, January 8, 1944.

― 16 ―

Title page of the original 1947 report which was declassified in 1959. Los Alimos documents used to list the persons who did the work as well as those who wrote the reports. This practice was abandoned soon after the end of the war.

― 17 ―

2—
Statistical Methods in Neutron Diffusion:
With J. von Neumann and R. D. Richtmyer (LAMS-551, April 9, 1947)

This report, written in 1947 by J. von Neumann and R. D. Richtmyer, is about work done by J. von Neumann and myself-as the title page indicates. It gives the first published ideas and proposals for the Monte Carlo Method. (Author's note.)

Abstract

There is reproduced here some correspondence on a method of solving neutron diffusion problems in which data are chosen at random to represent a number of neutrons in a chain-reacting system. The history of these neutrons and their progeny is determined by detailed calculations of the motions and collisions of these neutrons, randomly chosen variables being introduced at certain points in such a way as to represent the occurrence of various processes with the correct probabilities. If the history is followed far enough, the chain reaction thus represented may be regarded as a representative sample of a chain reaction in the system in question. The results may be analyzed statistically to obtain various average quantities of interest for comparison with experiments or for design problems.

This method is designed to deal with problems of a more complicated nature than conventional methods

― 18 ―

based, for example, on the Boltzmann equation. For example, it is not necessary to restrict neutron energies to a single value or even to a finite number of values, and one can study the distribution of neutrons or of collisions of any specified type not only with respect to space variables but with respect to other variables, such as neutron velocity, direction of motion, time. Furthermore, the data can be used for the study of fluctuations and other statistical phenomena.

THE INSTITUTE FOR ADVANCED STUDY Founded by Louis Bamberger and Mrs. Felix Fuld Princeton, New Jersey School of Mathematics March 11, 1947

VIA AIRMAIL: REGISTERED Mr. R. Richtmyer Post Office Box 1663 Santa Fe, New Mexico

Dear Bob: This is the letter I promised you in the course of our telephone conversation on Friday, March 7th.

I have been thinking a good deal about the possibility of using statistical methods to solve neutron diffusion and multiplication problems, in accordance with the principle suggested by Stan Ulam. The more I think about this, the more I become convinced that the idea has great merit. My present conclusions and expectations can be summarized as follows:

(1) The statistical approach is very well suited to a digital treatment. I worked out the details of a criticality discussion under the following conditions:

(a) Spherically symmetric geometry. (b) Variable (if desired, continuously variable) composition along the radius of active material (25 or 49), tamper material (28 or Be or WC), and slower-down material (H in some form). (c) Isotropic generation of neutrons by all processes of (b). (d) Appropriate velocity spectrum of neutrons emerging from the collision processes of (b), and appropriate description of the

― 19 ―

cross-sections of all processes of (b) as functions of the neutron velocity; i.e., an infinitely many (continuously distributed) neutron velocity group treatment. (e) Appropriate account of the statistical character of fissions, as being able to produce (with specified probabilities), say 2 or 3 or 4 neutrons.

This is still a treatment of "inert" criticality: It does not allow for the hydrodynamics caused by the energy and momentum exchanges and production of the processes of (b), and for the displacements, and hence changes of material distribution, caused by hydrodynamics; i.e., it is not a theory of efficiency. I do know, however, how to expand it into such a theory (cf. (5) below).

The details enumerated in (a)-(e) were chosen by me somewhat at will. It seems to me that they represent a reasonable model, but it would be easy to make them either more or less elaborate, as desired. If you have any definite desiderata in this respect, please let me know, so that we may analyze their effects on the set-up.

(2) I am fairly certain that the problem of (1), in its digital form, is well suited for the ENIAC. I will have a more specific estimate on this subject shortly. My present (preliminary) estimate is this: Assume that one criticality problem requires following 100 primary neutrons through 100 collisions (of the primary neutron or its descendants) per primary neutron. Then solving one criticality problem should take about 5 hours. It may be, however, that these figures (100 x 100) are unnecessarily high. A statistical study of the first solutions obtained will clear this up. If they can be lowered, the time will be shortened proportionately.

A common set-up of the ENIAC will do for all criticality problems. In changing over from one problem of this category to another one, only a few numerical constants will have to be set anew on one of the "function table" organs of the ENIAC.

(3) Certain preliminary explorations of the statistical-digital method could be and should be carried out manually. I will say somewhat more subsequently. (4) It is not quite impossible that a manual-graphical approach (with a small amount of low-precision digital work interspersed) is feasible. It would require a not inconsiderable number of computers for several days per criticality problem, but it may be possible, and

― 20 ―

it may perhaps deserve consideration until and unless the ENIAC becomes available. This manual-graphical procedure has actually some similarity with a statistical-graphical procedure with which solutions of a bombing problem were obtained during the war, by a group working under S. Wilks (Princeton University and Applied Mathematics Panel, NDRC). I will look into this matter further, and possibly get Wilks' opinion on the mathematical aspects.

(5) If and when the problem of (1) will have been satisfactorily handled in a reasonable number of special cases, it will be time to investigate the more general case, where hydrodynamics also comes into play; i.e., efficiency calculations, as suggested at the end of (1). I think that I know how to set up this problem, too: One has to follow, say, 100 neutrons through a short time interval At; get their momentum and energy transfer and generation in the ambient matter; calculate from this the displacement of matter; recalculate the history of the 100 neutrons by assuming that matter is in the middle position between its original (unperturbed) state and the above displaced (perturbed) state; recalculate the displacement of matter due to this (corrected) neutron history; recalculate the neutron history due to this (corrected) displacement of matter, etc., etc., iterating in this manner until a "self-consistent" system of neutron history and displacement of matter is reached. This is the treatment of the first time interval At. When it is completed, it will serve as a basis for a similar treatment of the second time interval At; this, in turn, similarly for the third time interval At; etc., etc.

In this set-up there will be no serious difficulty in allowing for the role of light, too. If a discrimination according to wavelength is not necessary, i.e., if the radiation can be treated at every point as isotropic and black, and its mean free path is relatively short, then light can be treated by the usual "diffusion" methods, and this is clearly only a very minor complication. If it turns out that the above idealizations are improper, then the photons, too, may have to be treated "individually" and statistically, on the same footing as the neutrons. This is, of course, a non-trivial complication, but it can hardly consume much more time and instructions than the corresponding neutronic part. It seems to me, therefore, that this approach will gradually lead to a completely satisfactory theory of efficiency, and ultimately permit prediction of the behavior of all possible arrangements, the simple ones as well as the sophisticated ones.

(6) The program of (5) will, of course, require the ENIAC at least, if not more. I have no doubt whatever that it will be perfectly

― 21 ―

tractable with the post-ENIAC device which we are building. After a certain amount of exploring (1), say with the ENIAC, will have taken place, it will be possible to judge how serious the complexities of (5) are likely to be.

Regarding the actual, physical state of the ENIAC my information is this: It is in Aberdeen, and it is being put together there. The official date for its completion is still April 1st. Various people give various subjective estimates as to the actual date of completion, ranging from mid-April to late May. It seems as if the late May estimate were rather safe.

I will inquire more into this matter, and also into the possibility of getting some of its time subsequently. The indications that I have had so far on the latter score are encouraging.

In what follows, I will give a more precise description of the approach outlined in (1); i.e., of the simplest way I can now see to handle this group of problems.

Consider a spherically symmetric geometry. Let r be the distance from the origin. Describe the inhomogeneity of this system by assuming N concentric, homogeneous (spherical shell) zones, enumerated by an index i = ,..., N. Zone No. i is defined by ri-_<r<ri, the ro, rl, r2,...,rN-1, rN being given: 0 = ro < rl < r₂< ... < rN-1 < rN = R, where R is the outer radius of the entire system.

Let the system consist of the three components discussed in (1), (b), to be denoted A, T, S, respectively. Describe the composition of each zone in terms of its content of each of A, T, S. Specify these for each zone in relative volume fractions. Let these be in zones Nos. ixi, Yi, Zi, respectively.

Introduce the cross sections per cm³ of pure material, multiplied by l⁰ loge = .43..., and as functions of the neutron velocity v, as follows: Absorption in A, T, S:EaA(),EaT(V), EaS(o) Scattering in A, T, S: EAA(v),ErT(V), EsS(v) Fission in A, with production of 2, 3, 4 neutrons: fA(V), (.fA(V), (4) .V

Scattering as well as fission are assumed to produce isotropically distributed neutrons, with the following velocity distributions:

― 22 ―

If the incident neutron has the velocity v, then the scattered neutrons' velocity statistics are described for A, T, S, by the relations V = V=(A(V), V = VT(V), V' = Vs(V) Here v' is the velocity of the scattered neutron, (OA(v), ((T(V), (PS(V) are known functions, characteristic of the three substances A, T, S (they vary all from 1 to 0), and v is a random variable, statistically equidistributed in the interval 0, 1.

Every fission neutron has the velocity v_o. I suppose that this picture either gives a model or at least provides a prototype for essentially all those phenomena about which we have relevant observational information at present, and actually for somewhat more. It may be expected to provide a reasonable vehicle for the additional relevant observational material that is likely to arise in the near future. Do you agree with this?

In this model the state of a neutron is characterized by its position r, its velocity v, and the angle 0 between its direction of motion and the radius. It is more convenient to replace 0 by s =r cos 0, so that v/r s² is the "perihelion distance" of its (linearly extrapolated) path. Note that if a neutron is produced isotropically, i.e., if its direction "at birth" is equidistributed, then (because space is three-dimensional) cos will be equidistributed in the interval -1, 1, i.e., s in the interval -r, r.

It is convenient to add to the characterization of a neutron explicitly the No. i of the zone in which it is found, i.e., with ri_-< r < ri. It is furthermore advisable to keep track of the time t to which the specifications refer.

Consequently, a neutron is characterized by these data: i, r, s, v,t .

Now consider the subsequent history of such a neutron. Unless it

suffers a collision in zone No. i, it will leave this zone along its straight path, and pass into zones Nos. i + 1 or i - 1. It is desirable to start a "new" neutron whenever the neutron under consideration has suffered a collision (absorbing, scattering, or fissioning-in the last-mentioned case several "new" neutrons will, of course, have to be started), or whenever it passes into another zone (without having collided).

Consider first, whether the neutron's linearly extrapolated path goes forward from zone No. i into zone No. i + 1 or i - 1. Denote these two possibilities by I and II.

― 23 ―

If the neutron moves outward, i.e., if s > 0, then we have certainly I. If the neutron moves inward, i.e., if s < 0, then we have either I or II, the latter if, and only if, the path penetrates at all into the sphere ri- 1. It is easily seen that the latter is equivalent to s2 >>r² r2_1^. So we have: s > O.-. A s< {r' 1 + S2 _ r2 < O B'} . I . s < 0.'- B{ri2_^qs2 _²^- .B ^O . r?+s² _r² > '0B"} 11

The exit from zone No. i will therefore occur at *= ri for I . = ri-1 for II. It is easy to calculate that the distance from the neutron's original position to the exit position is d = s*- s, where *--2 + for II

The probability that the neutron will travel a distance d¹ without suffering a collision is ₁₀^-fd, where f =(a ) + EA(v) + ) + ) + ) (v)) X + (aT() ZT(v)) Yi+ (aS⁽ V) + Ea (V)) Zi

It is at this point that the statistical character of the method comes into evidence. In order to determine the actual fate of the neutron, one has to provide now the calculation with a value A, belonging to a random variable, statistically equidistributed in the interval 0, 1; i.e., A is to be picked at random from a population that is statistically equidistributed in the interval 0, 1. Then it is decreed that ₁₀^- fd has turned out to be A; i.e., d'10 log Af -logA *

From here on, the further procedure is clear. * Today this older notation would read d¹ = -log₁₀ A/f. (Eds.)

― 24 ―

If d'>d, then the neutron is ruled to have reached the neighboring zone No. i ± 1 ( + fr II) without having suffered a collision. The "new" neutron (i.e., the original one, but viewed at the interzone boundary, and heading into the new zone), has characteristics which are easily determined: i is replaced by i ± 1, r by r*,s is easily seen to go over into s*, v is unchanged, t goes over into t* = t+d/v. Hence, the "new" characteristics are i ± 1,r*,s*,v,t*.

If, on the other hand, d¹ < d, then the neutron is ruled to have suffered a collision while still within zone No. i, after a travel d¹ . The position at this stage is now r* = /r² + 2sd + (dl⁾² , and the time t* =t+.v

The characteristic contains, accordingly, at any rate i,r*, and t* in place of i, r, and t. It remains to determine what becomes of s and v.

As pointed out before, the "new" s will be equidistributed in the interval -r*, r*. It is therefore only necessary to provide the calculation with a further value p', belonging to a random variable, statistically equidistributed in the interval 0, 1. Then one can rule that s has the value = r*(2p' - 1)

As to the "new" v, it is necessary to determine first the character of the collection: Absorption (by any one of A, T, S); scattering by A, or by T, or by S; fission (by A) producing 2, or 3, or 4 neutrons. These seven alternatives have the relative probabilities fi =EaA(V)Zi + EaT⁽ v)yⁱ⁺ EaS(v)Z f2 - fl= E A(V)Xi f3- f2 =YsT(1V)i4-f3sS(v)zi f-f4= f((V)xi f6_f5 = fA (v)xii f - Z 4 (V) f-f26= ()x -

― 25 ―

We can therefore now determine the character of the collision by a statistical procedure like the preceding ones: Provide the calculation with a value p belonging to a random variable, statistically equidistributed in the interval 0, 1. Form fi = -f; this is then equidistributed in the interval 0, f. Let the 7 above cases correspond to the 7 intervals 0, fl; f, f2; f2, f3; f3, f4; f4, f5; f5, f6; f6, f, respectively. Rule, that one of those 7 cases holds in whose interval f actually turns out to lie.

Now the value of v can be specified. Let us consider the 7 available cases in succession.

Absorption: The neutron has disappeared. It is simplest to characterize this situation by replacing v by 0.

Scattering by A: Provide the calculation with a value v belonging to a random variable, statistically equidistributed in the interval 0, 1. Replace v by V = VpA(V) Scattering by T: Same as above, but v = VPT(V) Scattering by S: Same as above, but v^' = v(S(v) ·

Fission: In this case replace v by vo. According to whether the case in question is that one corresponding to the production of 2, 3, or 4 neutrons, repeat this 2, 3, or 4 times, respectively. This means that, in addition to the p',s' discussed above, the further p", s";p"', s"'; p"", s"" may be needed.

This completes the mathematical description of the procedure. The computational execution would be something like this: Each neutron is represented by a card C which carries its characteristics i, r, s, v, t, and also the necessary random values A, p,v,p pl?' ",,p, ,,

I can see no point in giving more than, say, 7 places for each one of the 5 characteristics, or more than, say, 5 places for each of the 7 random variables. In fact, I would judge that these numbers of places

― 26 ―

are already higher than necessary. At any rate, even in this way only 70 entries are consumed, and so the ordinary 80-entry punchcard will have 10 entries left over for any additional indexings, etc., that one may desire.

The computational process should then be so arranged as to produce the card C' of the "new" neutron, or rather 1 to 4 such cards C', C", C"', C"" (depending on the neutrons' actual history, cf. above). Each card, however, need only be provided with the 5 characteristics of its neutron. The 7 random variables can be inserted in a subsequent operation, and the cards with v = 0 (i.e., corresponding to neutrons that were absorbed within the assembly) as well as those with i = N +1 (i.e., corresponding to neutrons that escaped from the assembly) may be sorted out.

The manner in which this material can then be used for all kinds of neutron statistic investigations is obvious.

I append a tentative "computing sheet" for the calculation above. It is, of course, neither an actual "computing sheet" for a (human) computer group, nor a set-up for the ENIAC, but I think that it is well suited to serve as a basis for either. It should give a reasonably immediate idea of the amount of work that is involved in the procedure in question.

I cannot assert this with certainty yet, but it seems to me very likely that the instructions given on this "computing sheet" do not exceed the "logical" capacity of the ENIAC. I doubt that the processing of 100 "neutrons" will take much longer than the reading punching, and (once) sorting time of 100 cards, i.e., about 3 minutes. Hence, taking 100 "neutrons" through 100 of these stages should take about 300 minutes, i.e., 5 hours.

Please let me know what you and Stan think of these things. Does the approach and the formulation and generality of the criticality problem seem reasonable to you, or would you prefer some other variant? Would you consider coming East some time to discuss matters further? When could this be?

With best regards; Very truly yours, John von Neumann

― 27 ―

Tentative Computing Sheet

Data: (1)ri, Xi, Yi, Zi as functions of i = 1,..., N¹ . (ro = 0.) (2) EaA(V),aT(),EaS(),EsA(), EsT(V), Ess(vu), w(f f(3() ,fAfA(v as functions of v > 0, < vo^.2 (3) vo. (4) (A(vV), (T(V), (PS() as functions of v> 0, < 1.2 (5) -⁰ logA as function of A > 0, < 1.2 Card C: C1 i C₂r C3 s C4 V C5 t Random Variables: R₁ A R₂U R₃v R₄p' R₅pl R₆P"' R₇p"" 1 Tabulated. (Discrete domain.) 2 Tabulated, to be interpolated, or approximated by polynomials. (Continuous domain.)

― 28 ―