A Monte Carlo solution of multiple scattering in cold neutron scattering experiments

A Monte Carlo solution of multiple scattering in cold neutron scattering experiments

NUCLEAR INSTRUMENTS AND METHODS I0I (1972) 263-266; © N O R T H - H O L L A N D P U B L I S H I N G CO. A M O N T E CARLO S O L U T I O N O F M ...

272KB Sizes 0 Downloads 30 Views

NUCLEAR

INSTRUMENTS

AND

METHODS

I0I

(1972) 263-266; © N O R T H - H O L L A N D P U B L I S H I N G CO.

A M O N T E CARLO S O L U T I O N O F M U L T I P L E SCATTERING IN COLD N E U T R O N SCATTERING E X P E R I M E N T S H. KALLI

Department of Technical Physics, Helsinki University of Technology, Otaniemi, Finland Received 17 January 1972

This paper discusses a Monte Carlo solution of multiple scattering in cold neutron scattering experiments. In the calculations performed free gas cross-sections with the scatterer mass A = I, 2 and 3 are used to investigate the order of magnitude of multiple scattering in hydrogenous samples. The results seem to strengthen the

assumption that multiple scattering cannot be ignored. The expected leakage probability (ELP) technique is used in the Monte Carlo solution. Some other means of reducing the variance are studied.

1. Introduction

2. Discussion of the calculations 2.1. FREE GAS CROSS-SECTIONS

The primary starting point of the present work was the discrepancy found between the theoretical and experimental results in the inelastic scattering of cold neutrons by methane1'2). In ref. 2 the authors assume that "it is possible that one of the causes for this discrepancy is multiple scattering". In a previous work 3) conducted in our reactor laboratory the amount of multiple scattering was calculated theoretically using the analysis of ref. 4. The correction for multiple scattering is found to afford better agreement between the theoretical scattering model and experimental results. Similar results are also obtained in ref. 5 concerning the results of ref. 2. In this paper the order of magnitude of the multiple scattering of cold neutrons by thin slab samples of hydrogeneous materials was studied using the free gas model cross-sections. The bindings of hydrogen atoms in different molecules were imitated by performing the calculations for three different mass numbers A = 1,2 and 3. The second starting point of the work was an interest in Monte Carlo techniques. In an earlier work 6) the author discussed various variance reduction techniques in Monte Carlo solution. Some of these techniques have been applied in the present paper. (For Monte Carlo papers on the multiple scattering of neutrons see refs. 7 and 8.) In the following sections we shall first give the formula for the cross-sections utilized, then we shall report the method and the variance reduction techniques used in the computer programs and finally we shall list and discuss the results.

263

The following expressions used in our Monte Carlo programs for the differential scattering cross-section as(Ei ~ El; p) are, for example, derived in ref. 9: 4~\~/

~

x

where t¢z = 2 [ - E f + E i + 2 / . t 4 ( E f E i ) ' ] ,

= Ei-Ef,

Ei Ef p A

is the initial energy of the scattering neutron, is the final energy of the scattering neutron, is the cosine of the scattering angle, is the ratio of the scatterer atom mass to the neutron mass, a b is the microscopic bound atom cross-section, T is the temperature of the scatterer. (The units are chosen such that k = 1, h = I and the neutron mass = 1.) The total scattering cross-section as(Ei) has the expression

o"s

=

a----L-f .I(2 A8 r + 1)erf[~/(A~r) ] +

2 AcT (

264

H. KALLI

where af is the microscopic free atom cross-section: ~rf=20.3b,

au =

--

af,

eT = EilT.

I/v-absorption was assumed, % = 0.332 b. The scatterer temperature T = 0.01 eV corresponding to approximately 120 K was used in all the calculations. The macroscopic cross-sections were chosen such that the total cross-section for incoming 0.005 eV neutrons was 1.0 c m - 1 2.2. MONTE CARLO PROGRAMS The geometrical configuration used in the calculations is an infinite slab. The coordinate system is the three-dimensional xyz-system. The slab is situated perpendicularly to the z-axis, its surface planes are z = 0 and z = d, where d is the thickness of the slab. The incoming source neutron beam meets the plane z - - 0 perpendicularly at origin with a discrete energy Ein. (The value E i n = 0.005 eV was used in all the calculations.) The number I,(12o, Enut) of the nth-order scattered neutrons streaming out of the slab through the plane z = d in a certain direction ~2o with a certain energy Eout [(1/s • solid angle • eV)] can be written as a multiple integral iX)

In(~'~O'Eout)=f~f

~/n(r' ~'~' E) x

collisions and the respective I,($20, Eout)'s were calculated at 40 energy points in three directions. Samples of 1000 histories were used; the c.p.u, time for one fixed scatterer mass value and slab thickness was about three minutes. The estimated standard error is normally smaller than 5% and always below 10% in all the results obtained. 2.3. REDUCTION OF VARIANCE

Since all the calculations were made for relatively thin slabs (d < 0.2 m.f.p.) it is evident that the efficiency of the so-called analog Monte Carlo method (AMC) would be rather poor. (Most histories would leak out of the slab without any collisions.) The method of expected leakage probability (ELP) makes it possible to study limit values numerically when the slab thickness approaches zero. In this case A M C becomes completely ineffective. In A M C the distance s between two collisions is sampled from the well-known exponential distribution p (s) = ~ t (E) exp [ - 2~t(E)s];

If the distance in the direction of the flight to the outer boundary is a, then the leakage probability

PL = 1 --

/o 0

s,(e) where 7',(r,S2,E)

is the nth-order collision density,

Z s ( E , ~ E ' , f ~ ' ) is the macroscopic differential scattering cross-section, is the macroscopic total cross-section, p(r,12o) is the distance from the point r to the plane z = d in the direction I] o, V is the slab volume. This multiple integral is a standard type in Monte Carlo work. The Monte Carlo solution of this integral is discussed in detail in ref. 6 which also gives the original references. Here we are only concerned with the special features of the programs used. Three different program versions (with different variance reduction techniques) were written for a U N I V A C 1108 computer. Multiple scattering was studied by varying the scatterer mass and the slab thickness d. The histories were generated up to four z~t ( E )

p(s)ds = e x p [ - - I ; t ( E ) a ] .

The non-leakage probability

VO 4n × 2~s(E, ff'~---~Eout, Do) e-$,(Eout)a(,, ao) d r d E d l 2

s > 0, (homog. medium).

Pno, = 1 - PL. In ELP it is assumed that from every collision or source point a weight PLW leaks out of the slab, while a part P.o~W is retained. (w is the weight of a history.) Then the distance s between two "collision" points is sampled from the following normalized distribution

p*(s) = Xt(E) exp[-J~t(E)s];

0 < s < = a.

Pnon As a result of the application of ELP a history never leaks out of the slab. In a thin slab the weight very quickly decreases toward zero. In our basic version, which was used as a standard in efficiency comparisons, ELP was applied in sampling the path lengths. The energy after a collision was sampled from a multigroup distribution obtained by integrating the actual free gas scattering cross-section. The cosine/~ of the scattering angle was sampled from the distribution

p(l~)=A~exp(bip);

-l__
A M O N T E C A R L O S O L U T I O N OF M U L T I P L E

TABLE 2

TABLE 1 Ratio

Ratio Imulti//1 for A = 1. Eout (eW)

265

SCATTERING

Scattering angle 30 °

Eout '. (eV)

for A = 2.

15 °

Scattering angle 30 °

60 °

0.01 0.05 0.1 0.2

0.0035 0.0145 0.027 0.051

0.0067 0.028 0.052 0.097

0.0135 0.056 0.104 0.195

0.010

0.01 0.05 0.1 0.2

0.013 0.053 0.099 0.185

0.010 0.048 0.090 0.170

0.015 0.064 0.120 0.226

0.0155 0.067 0.126 0.24

0.015

0.01 0.05 0.1 0.2

0.027 0.112 0.210 0.394

0.020 0.083 0.155 0.290

0.019 0.079 0.148 0.279

0.0175 0.075 0.141 0.268

0.020

0.01 0.05 0.10 0.2

0.043 0.175 0.327 0.610

0.031 0.126 0.235 0.440

0.024 0.099 0.185 0.347

d (cm)

15 °

0.005

0.01 0.05 O. 1 0.2

0.0033 0.014 0.028 0.051

0.0064 0.027 0.050 0.095

0.0135 0.056 O. 105 O. 199

0.005

0.010

0.01 0.05 0.1 0.2

0.0075 0.032 0.06 0.114

0.0087 0.037 0.07 0.138

0.0140 0.060 0.113 0.215

0.015

0.01 0.05 0.1 0.2

0.012 0.052 0.099 0.191

0.0116 0.050 0.094 0.181

0.020

0.01 0.05 0.1 0.2

0.015 0.064 0.120 0.228

0.014 0.061 0.115 0.218

60 °

where At is a normalization factor and parameter b~ was chosen to give the correct average cosine of the scattering angle for every energy group. During the solution the deviations from the actual free gas probability distributions were compensated by the reduction of the history weight as described in ref. 6. In the second version importance sampling of the cosine Pl of the first scattering angle was combined with the basic version. The incoming source neutrons meet the slab perpendicularly. The cosine P l of the first scattering angle was sampled from distributions strongly favouring values near zero (i.e. paths parallel to the slab surface planes). Our aim was to reduce the standard error in the results of the second collision, but this was a minor success. A little efficiency was gained in higher energies but nothing in the lower energies compared with the calculations with the basic version. In the third version the antithetic variates technique was applied by performing the calculations in pairs of two negatively correlated neutron histories (cf. refs. 6 and l 0). This correlation was generated in the selection of the first path length. The final estimator was the average of the pair. This technique was found to be rather effective in calculating the results for the first collision. However, in the following collisions the correlation becomes weaker.

d (cm)

Imulti/ll

3. Results First let us repeat the values of constants in the calculations performed: scatterer temperature T = = 0 . 0 1 e V , energy of incoming source neutrons Ein=0.005 eV, macroscopic total cross-section of incoming n e u t r o n s ~ ' t ( E i n ) = 1.0cm -1. The three tables list the computed values of the ratio Imuiti/I1 4

where ImuJtl = ~ In" The scatterer mass A is used as the n=2

main classification factor. Inside each table the energy Eout of outgoing neutrons, the scattering angle and the slab thickness d are varied. The estimated standard error is normally below 5% and always smaller than 10% of the results. In fig. 1, 11, /multi and the s u m /tot = 11-[-/multi a r e shown as functions of the neutron wavelength for A -- 3, d - - 0.2 cm and the scattering angle = 15 °. 4. Discussion of the results The results of the calculations performed seem to strengthen the assumption that multiple scattering cannot be ignored in cold neutron scattering experiments. Our general result is thus parallel to the results in refs. 3 and 5. The proportion of multiply scattered neutrons seems to be relatively high, even for thin slabs, and it is therefore difficult to eliminate the effect completely by making the scatterer samples thinner. Fig. 1 shows that in spite of the peaked behaviour of function

266

H. K A L L I TABLE 3

Eout (eV)

d (cm)

A=3

50

Ratio Imulti/ll for A = 3.

15 °

Scattering angle 30 °

60*

0.0035 0.0146 0.027 0.051

0.0068 0.028 0.052 0.098

0.0135 0.056 0.103 0.194

kT : 0.01 eV Ein= 0.005 eV d= 0.2cm sccttering o

"~ 40 == .A 8 30 E v

0.005

0.01 0.05 0.1 0.2

t ~gLe:

15

2,

._,~°10 0.010

0.01 0.05 0.1 0.2

0.025 0.083 0.154 0.290

0.014 0.059 0.I 1 0.206

0.016 0.066 0.126 0.239

0.015

0.01 0.05 0.1 0.2

0.057 0.228 0.424 0.793

0.032 0.130 0.242 0.435

0.022 0.091 0.171 0.323

0.020

0.01 0.05 0.1 0.2

0.11 0.44 0.815 1.52

0.06 0.24 0.447 0.848

0.032 0.128 0.238 0.445

11 the multiple scattered component Imutti seems to be a relatively fiat function. The calculations performed clearly show that by using the ELP technique the Monte Carlo method can be effectively used in multiple scattering studies. Likewise the use of the antithetic variates technique seems promising. This work was supported by the National Research Council for Technical Sciences. The author would also like to thank Prof. E. Tunkelo for discussions in the course of the work.

/y, 2

3

4

5

WAVELENGTH (£)

Fig. 1. Dependence of /tot, 11 and /multi on the neutron wavelength.

References 1) G. Solt, Physica 37, no. 2 (1967) 253. 2) B. A. Dasannacharya and G. Venkataraman, Phys. Rev. 156, no. 1 (1967) 196. a) A. Tiitta and E. Tunkelo, Physica 54, no. 3 (1971) 393. 4) E. L. Slaggie, Nucl. Sci. Eng. 30 (1967) 199. 5) A . K . Agrawal and S.G. Das, Trans. Am. Nucl. Soc. 14 (1971) 351. 6) H. Kalli, On the variance reduction techniques in Monte Carlo solutions of neutron transport problems, Report TKK-F-A157 (Helsinki University of Technology, Department of Technical Physics, Otaniemi, Finland, 1971). 7) A. Choudry, Nucl. Instr. and Meth. 71 (1969) 221. 8) F. G. Bischoff, W. E. Moore and M. L. Yeater, Monte Carlo evaluation of multiple scattering and resolution effects in douple-differential neutron scattering cross section measurements, RPI-328-211 (1971). 9) M. M. R. Williams, The slowing down and thermalization o f neutrons (North-Holland Publ. Co., Amsterdam, 1966). 10) j . M . Hammersley and D . C . Handscomb, Monte Carlo methods (Methuen, London, 1965).