A Monte Carlo simulation code applied to diffraction experiments with polarization analysis

Neutron diffraction with polarization analysis is a very useful technique to separate the coherent and incoherent contributions to the scattering cross section. This experimental technique allows the direct measurement of the spin-flip and non spin-flip contributions, which in turn are related to the coherent and incoherent contributions by simple linear combinations. Nevertheless, the correction procedures to obtain the single coherent elastic scattering are tricky. We present a method, based on Monte Carlo simulation, that allows one performing these corrections in a wide variety of coherent and incoherent samples. In this work, we have applied this method to the study of heavy water, for which we obtained the interference structure factor. We describe the experimental setup to perform the measurements in a reactor source and we detail the experimental protocol allowing the determination of the non spin-flip and spin-flip diffractograms.


Introduction
The use of neutron diffraction techniques on hydrogenated materials is limited because of the huge incoherent scattering cross section of hydrogen ( 1 H or H). For this reason, diffraction experiments on this kind of systems are performed on deuterated samples, where H is substituted with deuterium ( 2 H or D). This substitution reduces the incoherent contribution, but it has two drawbacks. First, one looses the unique information provided by a nucleus with a negative scattering length, for example, negative peaks in the radial distribution functions. Second, if isotopic effects are expected, one can never be sure whether a given feature observed in the correlation function is coming from an isotopic effect or is indeed a real feature of the correlation function.
The best way to overcome this problem is directly measure the incoherent contribution over a relevant momentum transfer. In particular, for the case of liquid systems, this relevant momentum transfer range should extend beyond 15Å −1 in order to obtain relatively good Fourier transforms of the structure factors [1]. The access to this momentum transfer range is quite limited because only few instruments worldwide have enough hot neutron flux. But how can one measure the incoherent scattering directly? The answer is that such measurement is almost impossible, at least if the aim is the determination of the total incoherent contribution. Nevertheless, a fraction of the incoherent part, the spin-incoherent scattering due to the disorder of the nuclear spins, can be determined using polarized neutrons. While coherent scattering does not change the neutron spin (non-spin-flip scattering), spin-incoherent has the neutron spin reversed (spin-flip scattering). In the case of H, over 99% of the incoherent signal is due to spin-incoherent scattering. However, this technique has been little exploited so far in studies of liquids, as the Q-range allowed by the existing instruments was too small [2,3]. Now we present a study were we have used polarized hot neutrons in order to extract the spin-incoherence on a mostly coherent liquid: the heavy water. In such a system, the incoherent contribution is small and, consequently, the corrections are small too, and the results more accurate than in a system with more incoherent scattering. This rather simple case is used to assess the method, before extending it to highly incoherent (hydrogenated) systems. This ultimate evolution is still under development and will be treated in a forthcoming publication.
Our method is based on a Monte Carlo (MC) code simulating the interaction of the neutrons with the atoms, which allows us to follow the history of the individual neutron until its detection. This method is very useful in order to evaluate the corrections (multiple scattering, attenuation, inelasticity) as a whole and not like independent corrections as usual [4]. The drawback is that it requires a model of the interactions between the neutrons and the nuclei of the particular system giving the probability of scattering at a given angle and a given energy exchange. In the case of water, we use the Synthetic Model proposed by Granada [5], which has very well shown its usefulness in the past [6].

Experimental
The experiment was performed on the Spin Polarized Hot Neutron Beam Facility (D3) at the ILL (Institut Laue Langevin, Grenoble, France) [7]. This diffractometer is the only one worldwide allowing polarization analysis at short-wavelengths and covering a wide range of momentum transfer Q (0.8 to 21Å −1 at a neutron wavelength of 0.52Å). Figure 1 shows a sketch of the D3 setup with its different components. The polarizing monochromator is a Heusler crystal, Cu 2 MnAl ((111) Bragg reflection). Between the polarizer and the sample position, a guide field preserves the spin state of the neutron and a monitor allows normalization of the eventually detected intensities to the actual neutron flux. Before and after the sample nutators and slits respectively reverses (flips) the spin state of the neutron (whether necessary) and adjust the beam size. The neutrons scattered by the sample are analyzed using a polarized 3 He spin filter [8,9], which uses the drastically different transmissions of the neutrons with a spin state parallel or antiparallel to the spin state of the helium atoms to serve as a polarization analyzer. The transmitted neutrons then reach the single 3 He detector. The detector arm rotates around the sample, covering a total angular range of 4 • to 120 • . In our experiments, we chose steps of 0.5 • up to 26 • and 1 • for the higher angles and typical data acquisition times of 1 minute/step, half the time spent on measuring the non-spin-flip (NSF), spin-flip (SF) neutron intensities respectively. Finally a beam stop is placed in the direction of the direct beam from the monochromator, to stop the transmitted neutrons and minimize the radiation background. The measurement was performed using a sample of D 2 O in liquid state at room temperature, provided by the ILL chemistry lab. The sample was loaded into a cylindrical double walled vanadium cell of 10.7 mm external diameter, 0.15 mm walls thickness and 58 mm height; thus, the sample thickness was 1.2 mm, corresponding to a volume of 0.829 cm 3 in the beam. The mass was checked before and after the experiment, confirming that the container had no leakage.
Measurements of NSF and SF diffractograms were carried out for 24 hours. Shorter times were employed to measure the instrument background, the signal from the empty vanadium sample cell and to calibrate the polarization. All diffractograms were normalized to the monitor counts. The decay of the polarization (and efficiency) of the 3 He spin filter was monitored every few hours by measuring the SF and NSF intensities from the (111) reflection of a cylindrical germanium crystal, mounted on top of the vanadium sample cell and automatically driven into the beam via a vertical translation. The measured polarization efficiencies, between 0.65 and 0.75 with a relaxation time of the order of 80 to 90 hours, were used to correct the data for the time dependent analyzer efficiency. Finally all SF and NSF diffractograms were corrected for the absorption by the empty filter.

Data treatment
The first step in data processing consists in determining the experimental NSF and SF intensities (I S exp (Q) NSF,SF ) scattered by the sample, which can be determined in a first approximation as the difference between the sample+container minus the empty container intensities. Figure 2 shows these magnitudes as obtained from our experiment on D 2 O. It is worth noting that the SF empty cell intensity is comparable with the SF total intensity, making the background subtraction a delicate issue.
Our MC code was initially developed to correct standard (non-polarized) neutron diffraction data [12] and has now been extended to include the special case of polarized neutrons. Equations (25) and (26) in Ref. [12] are extended to the case of a polarized incident beam and an unpolarized sample, thus the angular distributions of scattered neutrons (SF and NSF) after the first scattering are given by the following expressions where I 1 (k 0 , θ) represents the angular distribution of single scattering without attenuation, which is equivalent to a point-like sample with the same scattering power as the whole sample, andĨ 1 (k 0 , θ) represents the angular distribution of neutrons detected after the first scattering. In these expressions, t(E 0 ) is the transmission of the neutron beam at energy E 0 , Σ t (k 0 ) is the macroscopic total cross section, Σ s (k 0 ) NSF,SF is the macroscopic scattering cross section, σs(E 0 ) NSF,SF is the normalized microscopic differential cross section, V the sample volume, A(k 0 ) is the cross sectional area exposed to the incident beam, σ(E 0 ,E,θ) NSF,SF σs(E 0 ) NSF,SF the normalized microscopic double-differential cross section, ε(E) in the detector efficiency at energy E, and H 1 (k 0 , k) is the primary attenuation factor defined by Sears [13]. This factor takes into account the attenuation of the incident neutron beam before the first scattering and then on its way towards the detector.
These magnitudes, Eqs. (1) and (2), are not accessible by the experiment, but they can easily be simulated. Using the neutron transport theory, Sears [13] determined the analytical expressions for the multiple scattering. These expressions, involving multiple integrals, make impractical their numerical quantification, which is why we perform this evaluation through MC numerical simulation.
Finally, the corrected data are normalized to absolute units (barn/sterad/atom) using the measured differential cross section from the vanadium sample container as a reference, and taking into account the difference in volume and density of the sample and vanadium, i.e., the different number of scatterers in the beam. The normalization constant is chosen in order to obtain the good high-Q limiting value for the vanadium cross section, that is, the free scattering cross section (4.83 barns).

Monte Carlo simulation
The MC simulation allows the assessment of the different components of the scattering distributions (single scattering, multiple scattering, attenuation factors and the total scattering), which in turn will serve to determine a Q-dependent correction factor. This factor is then employed to correct the NSF and SF experimental intensity, leading to their corresponding single scattering experimental intensity, as described by Eq. (1) and (2). The MC simulation is based on the method proposed by Bischoff [14] and Copley [15], which describes all the interaction processes between the neutron and the sample since entering in the sample until it is registered by the detector at a given angle. Discrete neutron histories are tracked and averaged in a random walk directed by interaction probabilities obtained either by measurements or by models.
As input, the MC program requires information on the sample and the experimental setup. For the sample, the following information is required: the macroscopic total cross section, the absorption probabilities as a function of energy, the scattering laws, the sample geometry, and the SF interaction probability [17], which in the case of D 2 O is p SF = σ SF s /σ s = 0.1395 [16]. For the instrument, the incident neutron energy and the beam dimensions are needed.

Neutron histories
The neutron trajectory starts at a randomly chosen point in the sample and the first collision is also randomly determined taking into account the path length distribution function. At each collision the new energy and flight direction are determined using probability distributions. For the energy exchange we usually use Granada's Synthetic Model [5], which allows the estimation of the normalized double-differential cross section. In this work we employ a hybrid model which works very well for a variety of coherent and incoherent samples. This model has been previously tested to perform the data treatment in standard neutron diffraction experiments [4,12] and now has been extended for the case of diffraction experiments using polarized neutrons. For a given collision the flight direction is determined using the experimental angular distribution, whereas the energy after the collision is determined using Granada's model.

Neutron weight
The implemented MC program is based on an implicit simulation, which is more efficient and reduces the computing time as compared to analogue simulations. In the implicit simulation the neutron never gets out the sample being forced to scatter (since absorption is forbidden in this altered scheme) [14]. In this scheme, to compensate the biases in the probability distribution employed, a weight is assigned to each neutron, which decreases with the transmitted fraction along the path, the initial value being 1. The neutron history is finished when the weight drops under a predetermined cut-off value and after applying the Russian Roulette criterion [14]. For the NSF and SF scattering processes, the neutron weight of the collision i can be calculated from the weight of the preceding collision i − 1 as follows where 1−t(E, d) represents the probability of interaction within the distance d, and Σs(E,0) Σt(E,0) gives the probability that the neutron is not absorbed when it experiences a NSF or SF scattering process.

Scoring
For each collision process, the contribution of the current history to the final diffractogram is scored by the detector placed at an angle θ with respect to the incident beam. In the single scattering case, the contribution of a point r to the (NSF, SF) neutron distribution after the first scattering can be calculated by scoring the estimator where D(Q) is the normalized experimental angular distribution expressed in the elastic-Q scale (Q e = 2k 0 sin(θ/2)). Similarly (NSF, SF) neutron distribution detected after the i-th scattering can be calculated by scoring the estimator where t(E i−1 , r i , −k) is the transmission coefficient from the point r i inside the sample to the sample surface in the directionk to the detector position. It is noteworthy that for the special case i = 1, the expressions (4) and (5) averaged over a large number of histories converges to the expressions (1) and (2). Using this formalism, the multiple scattering contribution can then be calculated as:Ĩ

Correction factor
Once the different contributions are obtained by means of the numerical simulations, the NSF and SF correction factors can be calculated as where A(k 0 , θ) =Ĩ 1 (k 0 ,θ) NSF,SF I 1 (k 0 ,θ) NSF,SF is the attenuation factor. This correction factor is then used to correct the experimental intensities, thereby obtaining the microscopic differential cross section according to the expression (1), which is the sought magnitude in all diffraction experiments.

Iterative procedure
The correction process is based on an iterative method, where the outputs in each run are employed to calculate the new experimental angular distributions, in its turn serving as new input in the next run until no variation in the multiple scattering components is observed. In this case, with heavy water as the sample, the convergence is achieved after 3 iterations and each run consists in 200000 neutron histories. The iterative scheme implemented to correct the experimental distributions is where the subscript i indicates the iteration number, and D 0 (Q) NSF,SF is the measured experimental distribution. Figure 3 shows the NSF (left) and SF (right) experimental intensities for D 2 O and its different scattering components as obtained in the last iteration (single, multiple and total scattering). We found a good agreement between the experimental data and simulations and we also observe that almost all the measured signal is due to the single scattering (∼ 95 %). This is explained by the fact that the mean free path of the incident neutron (0.302 eV) is about 27 mm which is much larger than the total sample thickness (2.6 mm). The oscillations observed in the SF signal (right panel in Fig. 3) have a statistical origin due to the low incoherence of D 2 O. Once the NSF/SF correction factors are determined, the experimental intensities corrected by multiple scattering and attenuation can be expressed as

Results and discussion
A simple linear combination of the corrected NSF and SF intensities then gives the corrected coherent I corr coh (Q) and incoherent I corr inc (Q) experimental intensities: The left panel of Fig. 4 shows the corrected NSF and SF experimental intensities, while the right panel shows the coherent and incoherent experimental intensities.
The expression (10), obtained after the experimental corrections and normalization process, is still affected by coherent inelastic effects and isotope incoherent scattering (the latter is also NSF and is not eliminated by the method). In order to correct for these effects, we use a Placzek polynomial [18], as follows: Finally, subtracting this polynomial (Fig. 4, right panel) from the total coherent contribution (Eq. (10)), the normalized coherent structure factor F N (Q) is obtained. In Fig. 5 we compare our coherent structure factor with that obtained by Soper [19]. We observe a good overall agreement, specially in what concerns the intensity in absolute units and the peak positions, except for the low-Q region and the peak close to 4Å −1 . These small discrepancies are likely related to problems in the background subtraction, having more impact on the low-angle region of the diffractogram and which needs to be improved. However, they do not prevent the demonstration of the usefulness of the method.

Conclusions
We presented here a consistent method to correct data from polarized neutron diffraction experiments, which can be employed in a wide variety of coherent and incoherent samples. The proposed method corrects for multiple scattering, attenuation effects and inelasticity in a consistent way via the MC simulation of neutron trajectories, and this for both NSF and SF diffractograms. At present this method is applied to the case of D 2 O where the incoherent contribution is low, but it is currently being extended to highly incoherent samples. The results  Figure 5. The structure factor of D 2 O after applying our correction method. Our results are compared with those obtained by Soper [19].
we obtained for D 2 O are comparable with those obtained by Soper [19], showing that our procedure produces consistent data. Regarding discrepancies observed in Fig. 5 at low-angle region, a modification of the MC code is foreseen in order to take into account the attenuation of the container intensity at small angles, which will improve the outcome of our method.
Concerning the instrument D3, we show its usefulness in the precise determination of the incoherent contribution, which in turn conducts to a better determination of the interference structure factor. The main drawback is the limited counting rate for an instrument that was designed for diffraction experiments on single crystal. In order to improve the registered intensity and to reduce the counting time, a possible solution would be the use of a larger detector covering a wider range of scattering angles, as is the case of D4 [20] at ILL. Increasing the solid angle covered by the detector would increase the statistics and thus reduce the error bars.
It is important to note that a Placzek correction is still necessary, because we can not take into account the inelastic coherent scattering. The main improvement of this method is the ability of measuring a big fraction of the incoherent scattering, which makes the Placzek correction smaller, producing more accurate structure factors.
At the core of the MC method there is the knowledge of a model describing the energy exchanges between the neutron and the nuclei. We used Granada's Synthetic Model [5] for water, but it could be extended for other hydrogenated molecular liquids.
We hope that this work, and others currently in preparation or submitted [21], will help showing that the use of polarized neutrons is extremely useful for the study of hydrogenated liquids. This technique is particularly interesting for biological systems, where the isotopic substitution of hydrogen with deuterium is not always possible or desirable.