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

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

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

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

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:

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.

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 _^2- .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.)

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 -

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

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

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.)

TENTATIVE COMPUTING SHEET (Cont'd.) Calculation: Instructions: Explanations:

1 r of C1 - 1,see (1) ri-1 2 rofCl,see(1) ri 3 (C3)2 S² 4 (C2)2 r 5 3 - 4 2 -r² 6 3>0.'.A f > 0'.A < C <O.-'.B<0.-.B Only for B: 7 (1)2 r² Only for B: 8 5 + 7 r² _] + - 2 Only for B: 9 8 > O.'. B"2 2 2 2 > B" >my lor ^. } 0.'.B" ^i-i ^ + s> 0 .B" 10 A or B'..2A or B' '. ri - r * B".'. 1 B".'. i- = A or B'.. +1 A or B'.'.+1 = B". .-1 B".'. -1 = 12 (10)2 r*2 13 5 + 12 r*2 + s2 _ r2 14 11(sign) x 3 s* 15 14-C3 d 16 x of C1, see (1) Xi 17 y of C1, see (1) i 18 z of C1, see (1) zi 19 > oaA of C4, see (2) EaA(v 20 16 x 19 EaA(V)Xi 21 EaT of C4, see (2) EaT(V) 22 17 x 21 EaA(V)Yi 23 20 + 22 EaA(V)Xi + EaT(V)yi 24 EaS of C4, see (2) EaS(V) 25 18 x 24 EaS(V)Zi

TENTATIVE COMPUTING SHEET (Cont'd.) Calculation: Instructions: Explanations:

26 23 + 25 fl = ,aA(V)Xi+ EaT(V)yi+ EaS(V)Zi 27 -aA of C4, see (2) ZaA(v) 28 16 x 27 f2 - fl =E,A(V)Xi 29 26 + 28 f2 30 ZsT of C4, see (2) EsT(v) 31 17 x 30 f3 - f2 = EsT(v)Yi 32 29 + 31 f3 33 Es$of C4, see (2) Ess(v) 34 18 x 33 f4 - f3 = EsS(V)Zi 35 32 + 34 f4 36 (2) of C4, see (2) E>(2) 37 16 x 36 f5 - f4 = E(A(U)Xi 38 35 + 37 f 39 :(3) of C₄ see (2) E(3) ( 40 16 x 39 f6 - f5 -= f(V)Xi 41 38 + 40 f6 42 >(4) of C4, see (2) (4)(v) 43 16 x 42 f - f6 = Z (v)xi 44 41 + 43 f 45 -¹⁰ log of R₁, see (5) -¹⁰ logA 46 45:44 d¹ 47 46>15.'.P >d..P 47 46< 15.'.Q d < d.'. Q 48 P 15 4C4 P.'.d:} Q Q.46:^ Q.'.:]d' 49 C5 + 48 t* = t+

TENTATIVE COMPUTING SHEET (Cont'd.) Calculation: Instructions: Explanations:

Only for P: 50 C1+ 11 i* = i + Only for P: 51 C' :50Ci zi* C2: 10 C2:r* C': 14 C:s* 4:C4 C:v C5:49C5 : t* From here on only Q: 52 R ₂ x44 = ,f 53 52 < 26.'.Q1 < fi'. Q1 { < 29 } * Q2 {<f2 < .>{<29 } .' Q3 {->f2 .. < 32 Q < f3 1 ' >{<325} Q4 <qf3_!.Q₄ < 35 J {< f4 > 35I .Q > i f4 Q0 < 38 <5A {>3 } . Q6 {> }A Q6 <41 - <{<f6 ' > 41.'.Q7 _ f6 'Q7 54 C3 x 46 sd¹ 55 2 x 54 2sd¹ 56 (46)2 (dl)² 57 55 + 56 2sd¹ + (d¹ )² 58 4 + 57 r² + 2sd¹ + (d¹ )² 59 58 r*

TENTATIVE COMPUTING SHEET (Cont'd.) Calculation: Instructions: Explanations:

Only for Qi: 60 C: Ci C: i C2: 59 C2:r* C3:. C3.. c4:0C4:0 C5:49C5: t* From here on only Q2,. ., Q7: Only for Q2: 61 OA of R₃, see (4) PA(v) = ( Only for Q3: 62 YT of R ₃, see (4) rPT() = V Only for Q4: 63 ps of R ₃, see (4) $s(v) =) Only for 64 C4 x (61 or vo Q2,Q3, Q4 62 or 63) 65 Q2,Q3, Q4 .64 Q2,Q3,Q4-'-V =lv Q5,Q6, Q7 (3) Q₅,Q₆,Q7.'. V 66 2x R₄2p^t 67 66- 1 2p^- 1 68 59 x 67 s' 69 C1: C1 C :i C2 59 C2 r* C3: 68 C3:' C4: 65 C4:v C5:49C5: t

From here on only Q5, Q6, Q7: 70 2x R₅2p" 71 70 - I 2p"- 1 72 59 x 71 s" 73 C: C1C:i C2': 59 C2:r* 3': 72 C' :s"C" ^: 65 C4^:v' C" :49 C':t*

TENTATIVE COMPUTING SHEET (Cont'd.) Calculation: Instructions: Explanations:

From here on only Q6, Q7: 74 2x R₆2p"' 75 74 - 1 2p"'^- 1 76 59 x 75 s'W 77 C' : C₁C1^' i C2": 59 C2": r* C3/: 76 C3": s"' C4": 65 C4": v' C': 49 C5": t^*

From here on only Q7: 78 2x R7 2p"" 79 78 - 1 2p"" - 1 80 59 x 79 s"" 81 C"tt Cl C0'": i C2"': 59 C2"' :r* C3"': 80 C3"' : s"" C4': 65 C'": v' C5"': 49 C5^" : t*

At right, sample page of John von Neumann's handwritten "tentative computing sheet" as it was used in this 1947 report.

Instruationa iExpluations t

April 2, 1947 Professor John von Neumann, The Institute for Advanced Study, School of Mathematics Princeton, New Jersey Dear Johnny:

As Stan told you, your letter has aroused a great deal of interest here. We have had a number of discussions of your method and Bengt Carlson has even set to work to test it out by hand calculation in a simple case.

It has occurred to us that there are a number of modifications which one might wish to introduce, at least for calculations for a certain type. This would be true, for example, if one wished to set up the problem for a metal system containing a 49 core in a tuballoy tamper. It seems to us that it might at present be easier to define problems of this sort than, for example, problems for hydride gadgets. It is not so much our intention to suggest that the method you are working on now should be modified as to suggest that perhaps alternative procedures should be worked out also. Perhaps one of us could do this with a little assistance from you; for example, during a visit to Princeton.

The specific points at which it seems to us modifications might be desired are as follows:

1. Of the three components A, T, S that you consider, only one is fissionable, whereas in systems of interest to us, there will be an appreciable number of fissions in the tuballoy of the tamper, as well as in the core material.

2. On the other hand, we are not likely for some time to have data enabling one to distinguish between the velocity dependence of the three functions (2) (), (3)fa (v), (v)

that you introduce so that for any particular isotope these might as well be combined into a single function of velocity with a random procedure used merely for determining the number of neutrons emerging. If there is a single such function of velocity for each of the three isotopes 25, 28, 49, the total number of function tables required would be the same as in your letter.

3. It is suggested that in the case of 25 or 28, one might wish to allow also for the possibility of one neutron emerging from fission. The dispersion of the number of neutrons per fission is not too well known but we think we could provide some guesses.

4. Because of the sensitive dependence of tamper fissions on the neutron energy spectrum, it might be advisable to feed in the measured fission spectrum at the appropriate point. This would, of course, require introduction of one or two additional random variables and would raise the nasty question of possible velocity correlation between neutrons emerging from a given fission.

5. Material S could, of course, be omitted for systems of this sort. On the other hand, when moderation really occurs, it seems to us there would have to be a correlation between velocity and direction of the scattered neutron.

6. For metal systems of the type considered, it would probably be adequate to assume just one elastically scattering component and just one inelastically scattering component. These could be mixed with the fissionable components in suitable proportions to mock up most materials of interest.

In addition, we have one general comment as follows: Suppose that the initial deck of cards represents a group of neutrons all having t = 0 as their time of origin. Then after a certain number of cycles of operations, say 100, one will have a deck of cards representing a group of neutrons having times of origin distributed from some earliest tl, to a latest, t₂ . Thus all of the multiplicative chains will have been followed until time t₁ and some of them will have been followed to various later times. Then if one wishes, for example, to find the spatial distribution of fissions, it would be natural to examine all fissions occurring in some interval At and find their spatial distribution. But unless the interval At is chosen within the interval (0,tl) one cannot be sure that he knows about all the fissions taking place on At, and the fissions that are left out of account may well have a systematically different spatial distribution than those that are taken into account. Therefore, if, as seems likely, t₁ << t₂ , it would seem to be necessary to discard most of the data obtained by the calculation. The obvious remedy for this difficulty would seem to be to follow the chains for a definite time rather than for a definite number of cycles of operation. After each cycle, all cards having t greater than some preassigned value would be discarded, and the next cycle of calculation performed with those

remaining. This would be repeated until the number of cards in the deck diminishes to zero.

These suggestions are all very tentative. Please let us know what you think of them.

Sincerely, R. D. Richtmyer. cc: S. Ulam C. Mark B. Carlson