Liquid scintillators neutron response function: a tutorial

This tutorial is devoted to the understanding of the different components that are present in the neutron light output pulse height distribution of liquid scintillators in fusion relevant energy ranges. The basic mechanisms for the generation of the scintillation light are briefly discussed. The different elastic collision processed between the incident neutrons and the hydrogen and carbon atoms are described in terms of probability density functions and the overall response function as their convolution. The results from this analytical approach is then compared with those obtained from simplified and full Monte Carlo simulations. Edge effect, finite energy resolution, light output and transport and competing physical processes between neutron and carbon and hydrogen atoms and their impact on the response functions are discussed. Although the analytical treatment here presented allows only for a qualitative comparison with full Monte Carlo simulations it enables an understanding of the main features present in the response function and therefore provides the ground for the interpretation of more complex response functions such those measured in fusion plasmas. Although the main part of this tutorial is focussed on the response function to mono-energetic 2.45 MeV neutrons a brief discussion is presented in case of broad neutron energy spectra and how these can be used to infer the underlying properties of fusion plasmas via the application of a forward modelling method.


Introduction
The measurement of the neutron emission from deuterium-deuterium (DD) and deuterium-tritium (DT) fusion reactions is one of the most important methods of assessing the performance of present and future fusion reactors. Since a neutron is released for each fusion reaction occurring in the plasma, the measurement of neutron flux emitted from the plasma is directly correlated to the fusion power. The emitted neutrons' energy spectrum is characterized by two main components at 2.45 and 14.1 MeV from DD and DT reactions respectively 1 . The neutron energy spectrum is affected, among other things 2 , by the fusion plasma operating conditions. For example, the broadening of the Gaussian energy spectrum for 2.45 and 14.1 MeV neutrons is proportional to the square root of the plasma fuel ion temperature due to the relative velocity distribution of the reactants [1]. In addition, the neutron energy spectrum is affected by the different additional heating schemes that are normally employed in fusion devices such as neutral beam injection and radio-frequency heating. For example, in the case of neutral beam heating the emitted neutron energy spectrum is the combination of a thermal component due to the plasma fuel ions, a beamthermal component from the reactions between beam and fuels ions and a beam-beam component. The relative intensity of these components is in turn affected by the plasma conditions. In present day devices, the beam-thermal component is dominating while in future fusion devices, the thermal component will be dominant. In the case of radio-frequency heating, neutrons with energies above the DD and DT energies are observed due to fusion reaction from ions accelerated to energies of a few MeV. The measurement of the spatial and temporal evolution of the neutron emission from fusion plasmas in terms of its flux and energy spectrum can therefore provide important information on the plasma itself which can be used to optimize fusion power production.
Different diagnostics are used to measure the neutron emission from fusion plasmas [1,2] but all rely on the conversion of the neutron into a charged particle that can then be detected. This conversion takes places either thanks to neutron induced nuclear reactions in which heavy charged fission fragments are produced 1 or via elastic scattering of neutrons with light atoms, typically hydrogen. In their simplest form, neutron diagnostics can just be used as counters where the measured counts are proportional to the number of neutrons emitted by the plasma and therefore to the fusion power. The proportionality constant is determined via the absolute calibration of such neutron counters and depends on several parameters such as the counter's efficiency and its position with respect to the neutron source (the plasma) and the fusion reactor [3,4]. Neutron spectrometers are more sophisticated diagnostics in which the measured neutron energy spectrum is linked in a non-trivial way to the neutron source. Since the energy of fission fragments does not reflect the energy spectrum of the incident neutron, fission chambers can not be used as spectrometers, if one excludes the very crude spectroscopic capability offered by threshold reactions. For this reason, neutron spectrometers are all based on the conversion of the neutron into a light recoil particle via elastic scattering. Neutron spectrometers can be distinguished by the way in which the scattered neutron and the recoil particle are processed. In compact spectrometers, the recoil particle deposits its energy into the scattering medium which, depending on the material, can emit a pulse of scintillation light that is detected: in this case the scatterer itself acts as the detector. The light emission is then converted into a voltage signal and the voltage pulse height spectrum generated by the recoil particles is measured 2 . Since the voltage signal is proportional to the energy deposited by the recoil particle and this depends on the incident neutron energy, a detector based on scintillation material can, in principle, be used as a neutron spectrometer. Large spectrometers, instead, are based on the measurement of the scattered neutron or recoil particle in a detector other than the scatterer. For example, the recoil protons ejected by neutron scattering on a thin, hydrogen-rich foil can be detected in an array of detectors after they have been momentum and energy separated [5]. Alternatively, the time difference between the two scintillation events generated by the same incident neutron in two spatially separated scintillators provides the neutron time of flight and therefore its energy [6]. Recently, time-of-flight spectrometers are also taking advantage of the information on the amount of energy deposited by the recoil particles in the scatterer and in the detector to suppress the contribution from unwanted random coincidences that, especially in fusion devices, typically affect such instruments [7].
The often complex relation between the incident neutron energy spectrum and the output signal from the detector is referred to as the detector response function. Since the neutron energy spectrum at the detector's location is not mono-energetic and since the response function is usually dependent on the incident neutron energy, it is necessary to determine the response function for all the neutron energies of interest. The resulting set of response functions (one for each incident neutron energy) is referred to as the detector response function matrix. 1 It is the knowledge of this response function that allows to infer the characteristics of the neutron energy spectrum and ultimately of the plasma itself. Two different approaches can be used to relate the measured neutron energy spectrum to the fusion plasma source: i) forward modelling and ii) inversion algorithms. Forward modelling relies on the accurate modelling of all the processes that in the plasma affect the neutron emission, of the neutron transport from the source to the detector, of the conversion of the incident neutron field into recoil particles and eventually into the detector output signal [8]. Inversion algorithms make no assumption or modelling on the neutron source which instead is obtained by different least-squares minimization methods based on the knowledge of the detector response function [9]. Both methods have advantages and disadvantages but the one thing they have in common is the requirement for a very well characterized detector response function. Response functions for scintillator detectors are usually measured experimentally with well characterized neutron sources [10,11] and interpreted with the help of dedicated Monte Carlo codes such as NRESP [12] or more general radiation transport codes such as MCNP [13] and MCNP-PoliMi [14] and GEANT4 [15] (which has the additional feature of simulating the scintillation light transport too). Very good agreement is found between measured and simulated neutron response functions for all the cited codes. For the purpose of this tutorial, NRESP will be used. In NRESP, the relevant cross-sections (and differential crosssections) for neutron interactions with hydrogen and carbon in the energy range 0.02 to 20 MeV are included. The detector material composition and geometry including the liquid scintillator housing and in the optical window connecting to the photo-cathode are modelled. The deposited energy is then converted into a light pulse height distribution taking into account the finite detector energy resolution and experimentally measured light output functions for all the generated particles in the medium (recoils and α-s from nuclear reactions on carbon).
The interpretation of the measured or simulated response functions for liquid scintillator can be quite difficult if only the light output pulse height distribution is given. For a detailed understanding of the origin of the different features present in the response function it is necessary to analyse the contributions from the individual processes occurring in the liquid scintillator such as single and multiple elastic scattering, nuclear reactions and so on. This information can be obtained from the Monte Carlo codes discussed above but this is usually not trivial. In addition, although these codes can provide the expected light output pulse height spectrum, for example, for neutrons that have collided three times with protons, they do no provide an explanation for why it has that particular shape. The aim of this tutorial is to provide such an explanation by means of an analytical derivation of the expected light output pulse height distribution combined with a very simple Monte Carlo code used for its verification. The analytical approach is limited to a few simple cases but it provides nevertheless the basic understanding of how the real response function arises from a combination of multiple individual neutron interactions. This tutorial focuses in particular on the detailed explanation of the different contributions to the response function of a liquid scintillator to incident mono-energetic neutrons. As it will be shown, the interpretation of the response function even in this simple situation is far from trivial. The results obtained in this particular case are easily generalized for a broad neutron energy spectrum. The tutorial is therefore structured as follows. Section 2 is dedicated to a brief overview of the scintillation process, of the resulting recoil particle light output response function and of a simplified Monte Carlo code used for the study of multiple elastic scattering. Section 3 is devoted to the study of the light response function in the case of single elastic scattering processes. Section 4 is dedicated to the determination of the light response function in the case of multiple elastic scattering of a neutron with particles of the same species (for example, only with protons or only with carbon atoms). Section 5 discusses the mixed double scattering case in which the neutron scatters elastically with two different particles (for example first with a carbon atom and then with a hydrogen atom). In section 6, the light response functions calculated with full Monte Carlo simulations are qualitatively interpreted with the help of the response functions derived in the previous two sections for each elastic scattering process. Section 7 briefly addresses some aspects affecting the light response function that have been neglected in the previous sections and elucidate the use of the response function matrix for the interpretation of the liquid scintillator response function when non mono-energetic neutron sources are present such as in the case of fusion reactions. In addition, the forward modelling method is also discussed. Final comments and remarks are presented in section 8.

Neutron response function of liquid scintillators
Scintillators are among the most common type of detectors both for γ-rays and neutron radiation detection and operate on the principle of induced fluorescent light emission upon the interaction of the radiation within the material. A detailed description of the physical principles, material composition and application of scintillators can be found in [16]. For some scintillation materials, light emission depends on the type of incident radiation and therefore it is possible to discriminate between γ-rays and neutrons, a feature that is essential in fusion application as neutron diagnostics often operates in mixed fields. Liquid scintillators based on organic materials such as benzene (C 6 H 6 ), tuolene (C 6 H 5 ·CH 3 ) and xylene (C 6 H 4 ·(CH 3 ) 2 ), that is hydrogen-rich materials, belong to this category. The main interaction mechanism between γ-rays and liquid scintillators is Compton scattering with the electrons while neutrons interacts by elastic scattering with the hydrogen and carbon nuclei 1 . In both cases a recoil particle is generated (an electron for γ-rays and a proton or carbon nucleus for neutrons). Coulomb interactions between recoil particles and the organic scintillator molecules result in the conversion of the recoil kinetic energy into molecular excitation energy. Part of this excitation energy is then dissipated via thermal quenching and part via fluorescent light emission in the UV region of the visible spectrum. Typical light pulses have a very fast rise time (few nanoseconds) and decay constants between 20 and 200 ns. The fluorescent light is then transported (either directly of via reflection with the liquid scintillator housing walls) to a photocathode: depending on the volume of the scintillator, significant attenuation of the fluorescent light in the scintillator medium itself can occur. The fluorescent light is then converted into electrons at the photocathode via the photoelectric effect and the initial few ejected photoelectrons are subsequently accelerated, focussed and multiplied via secondary electron emission by multi-stage dynodes photomultiplier resulting in large gains (10 5 -10 9 ). The conversion of photons into electrons and their multiplication is a process that, under the correct experimental conditions, is highly linear. Non-linearities occur for example if the photomultiplier is operated at high gains, at high counting rates and in presence of even weak magnetic fields. The electron current at the anode of the photomultiplier is converted into a voltage via a load resistor and the detector voltage output is fed into the acquisition system by co-axial cables (usually tens of meters in fusion experiments). Attenuation and distortion of the voltage signal in the cables will occur but these processes are linear and easily modelled. In present days fusion neutron diagnostics the detector voltage signal is digitized at very high sampling frequencies ranging from 250 MHz up to 4 GHz with high resolution from 10 to 14 bits. The ADC process can be considered linear if the integral an differential non-linearities are negligible which is often the case. Several conversion mechanism contribute therefore to the final signal that is measured for a single neutron interacting with the Figure 1 -Liquid scintillator light output as a function the recoil particle energy. Data adapted from [17].
detector. If the liquid scintillator is to be used as a spectrometer it is then important that a good linearity exists between the recorded signal and in the incident neutron energy. Under the correct experimental conditions all the above processes are linear but one: the fluorescent light output for recoil particles from neutron elastic collision is inherently non-linear as shown in figure 1. This nonlinearity complicates significantly the neutron response function as discussed in detail in sections 3, 4 and 5.
Examples of the response function calculated with NRESP of two liquid scintillators with the same diameter (6 cm) and two different thickness (0.5 and 3.0 cm) for a beam of incident monoenergetic neutrons with an energy of 2.45 MeV are shown in figure 2. As can be seen, the response function exhibits a rich series of features some of which depends on the detector thickness. For the purpose of this tutorial a Simplified monTe cArlo neutron Response funcTion Simulator (STARTS) has been written to calculate the probability density function of the energy and light output distributions of multi-scattered neutrons and of all recoil protons and carbon nuclei in any possible combination. Contrary to the full Monte Carlo codes discussed above, STARTS is specifically aimed at the calculation of the pulse height spectra for specific types of elastic collisions. The following simplifications have been made in STARTS: i) no cross-section dependence is included as the particles involved in the elastic scattering are selected by the user, ii) all neutrons undergo a fixed number of elastic scattering defined by the user, iii) the detector is considered infinite in size, iv ) recoil particles deposit all their energy in the detector, v ) the incident neutrons are mono-energetic and vi) no energy resolution broadening is included in the calculation of the response function and vii) the light output function used is taken from [12]. The first simplification implies that the relative intensity of the contribution to the response function from neutron-proton (np) and from neutron-carbon (nC) elastic collisions is neglected: this does not affect the shape of the response function. Simplifications iii) and iv ) do not affect in any qualitative way the response function while simplification vi) has an important effect which is however well known and easily understood. The impact of all these simplifications is briefly discussed in section 3.

Response function for single neutron scattering
Consider an incident neutron with initial energy E n,0 making an elastic collision with a proton, assumed to be at rest, as depicted schematically in panel (a) of figure 3. After the elastic collision the energies of the neutron and of the recoil protons will be E n,1 and E p,1 where the index "1" indicates that these quantities refer to the energy of the particles after the first collision. From classical mechanics, invoking the conservation of energy and linear momentum, it can be shown that in the case of a neutron colliding with a generic target particle these energies are given by: where the index t identifies the recoil target, θ is the scattering angle in the centre of mass reference system, α = (A − 1) 2 /(A + 1) 2 and A = m t /m n that is the ratio between the target and neutron masses. It is useful, especially for the discussion in the following sections, to observe that E n,1 /α represent the maximum possible energy from which a neutron with energy E n,1 could have originated in a single elastic scattering. Under the assumption of an isotropic cross-section for the elastic scattering in the centre of mass, the scattering angle can assume any value between 0 and π with equal probability 1 . It can then be shown that the probability of the scattered neutron to have an energy in the interval [E n,1 , E n,1 + dE n,1 ] is given by: where p(E n,1 ) is the Probability Density Function (PDF) of the continuous random variable E n,1 . Note that p(E n,1 ) is a properly normalized PDF 1 . In fact, the total probability of the neutron having an energy in the interval [αE n,0 , E n,0 ] after one collision is the sum (integral) of all the probabilities of the neutron having an energy in the interval [E n,1 , E n,1 + dE n,1 ] over all possible energies, which is equal to: as it can be easily verified by inserting equation (3) into the above expression. note that equation (4) is the total probability law for mutually exclusive events for continuous random variables. . Similarly, the probability of the recoil target of having an energy in the interval [E t,1 , E t,1 +dE t,1 ] is given by: According to equations (3) and (5), the probability of observing a scattered neutron and recoil target with energies in the range [E, E + dE] is non-zero and constant within the appropriate energy interval as shown in panel (b) of figure 3. In the case of np elastic scattering, α = 0 which implies that the recoil proton energy E p,1 can assume any value between 0 MeV ("grazing" collision) and E n,0 ("head-on" collision) with equal probability. In the particular case depicted in panel (b) of figure 3, the energy of the scattered neutron is E n,1 and the energy of the recoil proton is E p,1 = E n,0 − E n,1 . Since for mono-energetic neutrons with energy E n,0 the recoil protons can equally assume energies in [0, E n,0 ] it follows that, in this simplified scenario, the response function of the detector is the one depicted in panel (b) of figure 3, that is a "box-like" function. As can be seen from figure 1, a recoil proton of energy E p,1 = E n,0 = 2.45 MeV would generate a light output yield of approximately L ≈ 0.8 which is exactly the upper edge of the response functions shown in figure 2.
This simple "box-like" response function is modified by the non-linear dependence of the light output function on the recoil particle's energy. Following [16], if one assumes for the light output function the relation: then: and therefore: where k ≈ 0.21 MeV −β and β ≈ 3/2 give a good approximation to the light output function for recoil protons shown in figure 1 up to 3 MeV. Equations (5) (with α = 0) and (8) can be combined to obtain the light response function for the recoil protons: which implies that the light output response function increases for low light yields. Figure 4 shows the light response function for the 3 cm thick detector calculated by NRESP, by equation (9) and by STARTS. As can be seen, the overall features at low light output yields (L 0.2) and for L ≈ 0.8 are well described qualitatively 1 but not in the range in between.
The "bump" seen in figures 2 and 4 at intermediate L values can not be explained as the result of the detector finite resolution, of edge effects for detectors of finite size nor as a consequence of the contribution from nC scattering. The effect of the detector finite energy resolution is mainly to "smear out" the sharp edge at the maximum light output as shown by the red curves in figure  2. For a detector of finite size it is possible for the recoil protons generated near the outer surface of the scintillator to escape after having deposited in the medium only a fraction of their energy. The mean free path λ of 2.45 MeV protons in liquid scintillators is of the order of λ 0.1 mm. Even assuming that all the recoil protons generated within a distance λ from the outer surfaces escape, the overall light response function would be reduced in its amplitude by approximately 1 % but its shape would not be modified. The effect would be even smaller for recoil protons of lower energies. According to equation (2), the maximum energy for a recoil carbon, for which α ≈ 0.716, is E c,1 = (1 − α)E n,0 ≈ 0.7 MeV. The corresponding light output is L ≈ 0.01 (see figure 1) and therefore the contribution to the light response function from recoil carbon nuclei is confined to the very low end of the response function. On the scale used in figure 4 this is hardly visible.
It is possible to conclude therefore that the "box-like" response function for single np elastic scattering combined with the light output non-linearity together with the detector finite energy resolution describe quite well the simulated response function especially for the thin detector. However, for thick detectors this is not the case: in [16] this feature is described as the result of neutron double elastic scattering with protons in the scintillator. Section 4 is dedicated to the understanding of the origin and shape of this feature.

Response function for multiple neutron scattering
In the case of an incident neutron making two elastic collisions four different outcomes are possible as shown in figure 7: double elastic collision on two protons (track "c") or on two carbon nuclei (track "f") or mixed scattering first on a carbon nucleus and then on a proton (track "d") or the other way around (track "e"). As discussed in section 3, the contribution to the total light output from recoil carbon nuclei is negligible and therefore the double scattering on carbon nuclei can also be neglected and is not further discussed here. In a similar fashion, the total light output in the case in which the neutron first collides with a proton and then with a carbon nucleus is almost equivalent to a single elastic scattering with a proton. Tracks of type "e" therefore contributes to the response function as single np events described in section 3. The reponse function from double scattering on protons (track "c") is discussed in this section while track "d" (collision on carbon nucleus followed by collision on proton) is discussed in section 5. The discussion of the neutron double elastic scattering on proton is divided in four parts. The probability density function for the energy of the neutron after two elastic collisions is derived in section 4.1 while the probability density function for the energy of the second recoil target is derived in section 4.2. The probability density function for the total deposited energy in two elastic collisions by the neutron is derived in section 4.3 and finally the PDF for the corresponding total light output is obtained in section 4.4.

Doubly scattered neutron energy probability density function
The probability of a neutron with initial energy E n,0 having energy in the range [E n,2 , E n,2 + dE n,2 ] after two elastic collisions given that after the first collision it had an energy in the range [E n,1 , E n,1 + dE n,1 ] is given by the conditional probability: where p(E n,2 ∩E n,1 ) is the PDF of the joint events resulting in the energies E n,1 and E n,2 . Recalling equation (3), the probability p(E n,2 | E n,1 )dE n,2 is given by: which is equal to equation (3) with E n,0 and E n,1 replaced by E n,1 and E n,2 respectively. The probability p(E n,2 )dE n,2 of a neutron with initial energy E n,0 having energy in the range [E n,2 , E n,2 +dE n,2 ] after two elastic collisions regardless of which energy it had after the first collision can be obtained, using the law of the total probability, as the sum of the probabilities of all the possible combination of two collisions that could have resulted in E n,2 being in such range. Each joint event has a probability given by p(E n,2 ∩ E n,1 )dE n,2 which, using equation (10) and replacing p(E n,1 ) with equation (3) in combination with equation (11), can be written as: The total probability p(E n,2 )dE n,2 is then obtained by integrating over all possible energies E n,1 : Note that E n,2 ranges between the neutron initial energy E n,0 , corresponding to two subsequent "grazing" collisions, and the minimum energy α 2 E n,0 corresponding to two subsequent "head-on" collisions (θ = π). If E n,2 ∈ [αE n,0 , E n,0 ], then E n,1 can take any value between E n,2 and E n,0 (see panel (b) of figure 6) so that equation (13) results in: If E n,2 ∈ [α 2 E n,0 , αE n,0 ], then E n,1 can range only between αE n,0 and E n,2 /α since in a single collision it is not possible for a neutron with initial energy E n,1 ∈ [En, 2/α, E n,0 ] to reach a final energy in the range [α 2 E n,0 , E n,2 ] (see panel (c) of figure 6). In this case then, equation (13) gives: To summarize, the probability p(E n,2 )dE n,2 of a neutron with initial energy E n,0 having energy in the range [E n,2 , E n,2 + dE n,2 ] after two elastic collisions is given by: 0 otherwise (18) and the corresponding PDF is shown in panel (a) of figure 7. In the particular case of the two targets being protons, α = 0 and therefore E n,1 and E n,2 can both take any value in the interval [0, E n,0 ] which implies that equation (13) should be integrated in the interval: p(E n,2 )dE n,2 = dE n,2 En,0 resulting in:

Second recoil target energy probability density function
Consider now the two recoils particles generated from the first and second elastic scattering with a single neutron with energy E n,0 . The probability for the first recoil to have an energy in the range [E t,1 , E t,1 + dE t,1 ] is given by equation (5) while the probability of observing the second recoil particle with an energy in the range [E t,2 , E t,2 + dE t,2 ] given that after the first collision the neutron has an energy E n,1 is given by the conditional probability: The probability of observing the joint events in which the first scattered neutron has energy E n,1 and collides with a second target transferring to it the energy E t,2 is: Substitution of the corresponding expressions for the PDFs results in: The probability p(E t,2 )dE t,2 of observing a second recoil target with an energy in the range [E t,2 , E t,2 +dE t,2 ] is then obtained using the law of total probability, that is, by integrating equation (23) over all possible energies E n,1 that could result in second recoil particle to be in such energy range. Note that the second recoil particle can have an energy in the interval E t,2 ∈ [0, (1−α)αE n,0 ] regardless of E n,1 , so the corresponding probability p(E t,2 )dE t,2 is: If E t,2 ∈ [(1 − α)αE n,0 , (1 − α)E n,0 ] then the incident neutron can only have energies in the interval [E t,2 /(1 − α), E n,0 ] so the corresponding probability p(E t,2 )dE t,2 is: To summarize, the probability of the second recoil particle to have an energy in the interval [E t,2 , E t,2 + dE t,2 ] is given by: The PDF p(E t,2 ) is shown in panel (b) of figure 7. In the particular case of the two targets being protons, equation (28) reduces to: (29) Figure 8 shows the probability density functions for the energy of two recoil protons and two recoil carbon nuclei calculated according to equations (5) and (28)

Total deposited energy for doubly scattered neutrons
Most liquid scintillators are small in size compared to the distance travelled by neutrons with energies in the MeV range in the time required for the ADC to record a few samples (a few nanoseconds). If a neutron was to scatter twice within this time interval then the two recoil particles would be generated on a time scale so short that even present day fast data acquisition system would see the two collisions as a single event with an energy equal to the sum of the energy deposited by each individual recoil particle. As discussed at the beginning of section 4, carbon contributes to the light response function only for L ≪ 0.1 and therefore, even if the energy deposited can be a substantial fraction of the initial neutron energy, it is not discussed further.
Consider instead an incident neutron scattering elastically with two protons: the total energy deposited, in this case, is E d = E p,1 + E p,2 . The probability density function p(E d ) can be obtained by observing that if the first recoil particle deposits an energy E p,1 , then the probability of observing a deposited energy E d is equal to the probability of the second recoil proton to have energy E p,2 , that is: The probability of observing E p,1 and E d is given by: Integration over all possible energies E p,1 , recalling that p(E p,1 ) = 1/E n,0 , gives then the probability of observing E d 1 : Figure 9 shows the PDF for the total energy deposited by the two recoil protons from a neutron with E n,0 = 2.45 MeV calculated according to equation (33) and with the STARTS code.

Total light output for doubly scattered neutrons
The non-linear relation between the proton recoil energy deposited in the scintillator and the emitted light output transforms non-linearly the two continuous random variables E p,1 and E p,2 into the two continuous random variables L 1 and L 2 thereby transforming non-linearly their outcome space as well. In particular, the outcome space for the total deposited energy by the two recoil protons given by: becomes where if: (37) Figure 10 shows the effect of this non-linear transformation of the probability density function p(E p,2 | E p,1 ) into p(L 2 | L 1 ) for E n,0 = 2.45 MeV calculated by STARTS. For "head-on" collisions, the neutron energy is transferred only to one proton and the light output is L(E p ) ≈ 0.8. For θ < π, resulting for example in E p,1 = 0.425 MeV, then E p,2 ∈ [0, 2.025] MeV (see left panel of figure 10). The recoil proton with energy E p,1 = 0.425 MeV will result in a light pulse L 1 = 0.05 and therefore L 2 ∈ [0, 0.60] (see right panel of figure 10). The maximum light output in this case is L = 0.65 for a total deposited energy of 2.45 MeV. The insert on the right panel of figure 10 shows p(L 2 | L 1 = 0.05) calculated by STARTS. Note however that the observable quantity is not the total deposited energy E d but the light output L = L 1 + L 2 where L 1 and L 2 must satisfy the condition: (38) The probability density function for the total light output p L is then calculated as the convolution p L1 ⊗ p L2 1 . The first term, p L1 , can be written in terms of p(E p,1 ) observing that the probability p(E p,1 )dE p of a recoil proton to have an energy in the range [E p,1 , E p,1 + dE p,1 ] is equal to the probability p(L 1 )dL that the corresponding light output is in the range [L 1 , L 1 + dL 1 ] from which follows 2 : Replacing p(E p,1 ) with its corresponding expression (see equation (5)), equation (39) becomes: In a similar fashion, the PDF for the conditional probability p(L 2 | L 1 )dL 2 is then: and therefore, the probability density function p(L) of observing a total light output L resulting from two recoil protons with light pulses L 1 and L 2 is: The PDFs p(L 1 ) and p(L 2 | L 1 ) are non-zero only if L 1 and L 2 satisfy the condition given in equation (38). This implies that for a given E n,0 and L 1 there is a maximum value L 2,max for which this condition is satisfied and is given by: The dependence of L 2,max on L 1 is shown by the dashed black line in the left panel of figure 11:  figure 11). The corresponding deposited energy is shown by the blue line in the right panel of figure 11. Since in this case E d < E n,0 , then p(L 1 ) and p(L 2 ) are non-zero for all L 1 ∈ [0, L obs ]. Conversely, consider now the case in which the observed total light output is L obs = 0.65 (see red line in the left panel of figure 11). In this case, p(L 1 ) and p(L 2 ) are non-zero only if L 1 ∈ [0, L 1,a ] ∪ [L 1,b , L obs ] where L 1,a and L 1,b are the light outputs for which L obs = L 2,max (L 1 ). For L 1 ∈ [L 1,a , L 1,b ] it turns out that E p (L 1 ) + E p (L 2 ) > E n,0 (as shown by the red line on the right panel of figure 11) which is of course not physically possible. It is clear then, that as L increases so does the integral p(L) until L = L B,1 with: where L * 1 is the point of tangency, i.e.: The PDF p(L) reaches its maximum for L = L B,1 and goes to zero as L → 0 since the regions where the PDF p(L 1 ) and p(L 2 ) are non-zero become vanishing small 1 . The way in which the convolution  integral in equation (42) is calculated as a function of L is elucidated in figure 12 which shows how the integrand depends on L 1 for three different values of the total light output L. Figure 13 shows instead the PDF p(L) given by equation (42) for all possible values of the total light output, that is for L ∈ [0, L M ], compared with the results from STARTS. The predicted contribution to the total light output response function due to a neutron of a given initial energy making two elastic scattering with protons shown in figure 13 can be individuated in the NRESP response function shown in figure 4: the sharp knee for L ≈ 0.539 can now be understood in terms of the maximum energy repartition between the two recoil protons. A closer comparison between figures 4 and 13 reveals however that, for L < L B,1 , p(L) does not drop as predicted by equation (42) which indicates that double scattering on protons is not sufficient to reproduce NRESP response function. For this to happen, it is necessary to consider the contribution from triple np scattering. In a fashion similar to what has been done for the double scattering case, it is possible to write the PDF for the light output for the third recoil proton as: The probability density function of the sums L = L 1 + L 2 + L 3 is then calculated as the convolution and is shown in figure 14 together with the one calculated by STARTS. As can be seen, p(L) does not drop much for L ∈ [L B,2 , L B,1 ] where L B,2 corresponds to situation in which E p (L 1 ) + E p (L 2 ) + E p (L 3 ) = E n,0 . It is clear from the results shown in figures 4, 13 and 14 that a linear combination of the response functions from one, two and three recoil protons can reproduce NRESP output. This is however postponed until section 6 as the last important component to the total light response function, that is the one arising from recoil protons from a neutron that has undergone a prior elastic collision with a carbon atom (track "d" of figure 7), is discussed in the next section.

Response function for neutron scattering with different targets
From the discussion in sections 3 and 4.2 it is clear that, even if the light output from the recoil carbon makes a negligible contribution to the response function, the energy repartition between the recoil and scattered particles will affect the light output of the scattered proton in the second elastic collision. The corresponding energy and light output PDFs can be derived by generalizing to the case where an incident neutron with initial energy E n,0 makes two elastic collisions with atoms characterized by A 1 and A 2 such that A 1 > A 2 (and therefore α 1 > α 2 ). The probability of observing a neutron with energy between [E n,1 , E n,1 + dE n,1 ] after the scattering with A 1 is: The probability of the twice scattered neutron to have an energy in the range [E n,2 , E n,2 + dE n,2 ] after the scattering with A 2 is obtained by applying the law of total probability to the joint event p(E n,2 ∩ E n,1 ) = p(E n,2 | E n,1 )p(E n,1 ) where the integral is carried out over all possible E n,1 : The integration limits above depends on E n,2 possible ranges. In particular then: In the case α 1 > 0 and α 2 = 0, the PDF p(E n,2 ) then becomes: The probability density function for the energy of the recoil target A 1 is given by (5) with α replaced by α 1 if E A1,1 ∈ [0, (1 − α 1 )E n,0 ] and zero everywhere else. The PDF for the energy of recoil target A 2 is given by equation (49) but the integration limits are now given by: In this case, (1 − α 2 )E n,0 corresponds to A 2 's maximum possible energy if the neutron has lost no energy in the collisions with A 1 (grazing collision) while (1 − α 2 )α 1 E n,0 is A 2 's maximum possible energy if the neutron has lost the maximum energy possible in a "head-on" collision with A 1 . Integration of equation (49) with E n,1 in the intervals specified in (53) gives the probability for the second recoil target to have an energy in [E A2,1 , E A2,1 + dE A2,1 ]: In the case α 1 > 0 and α 2 = 0, the probability of observing the recoil proton with energy in the range [E p,1 , E p,1 + dE p,1 ] is therefore: Panel (a) of figure 15 shows the PDF for the energy of the recoil proton as a function of E p,1 calculated by equation (55) together with STARTS results. The corresponding PDF for the light output is given by applying equation (9) to equation (55) and is shown in panel (b).

Comparison with full Monte Carlo simulations
It is now possible to proceed to the qualitative comparison between a combination of the theoretical PDFs derived in the previous sections with the light output response function obtained by full Monte Carlo simulation such as NRESP. In this context, a full Monte Carlo simulation refers to a simulation in which the proper cross-sections of the relevant processes are taken into account, the the density of hydrogen and carbon atoms in the scintillation material are specified and a realistic geometrical model of the detector is used. Without this level of details, neither STARTS nor the analytical approach can properly estimate the relative importance (amplitude) of the different components constituting the light output response function although one can reasonably assume that the role of multiple np scattering becomes more important as the size of the detector increases and that collisions with carbon nuclei followed by collisions on hydrogen are always present given that the cross-sections for these processes are of the same order of magnitude. The weights of the light output response functions corresponding to the different elastic scattering processes are obtained by a simple multiple linear regression where the dependent variable is the overall response function calculated by NRESP. The relative amplitude of the different components so obtained has just an indicative (qualitative) nature and for a proper estimate of these weights one should exclusively use full Monte Carlo codes. With these words of caution well present in mind, the comparison between the NRESP light output response function and the one from the linear regression based on STARTS calculations is shown in figure 16 for both the "thin" and "thick" detectors. As can be seen, the light response function for the "thin" detector can be well matched neglecting the contribution from triple np scattering while this component is essential to match the response function for the "thick detector".
As an example of how the analytical approach can be used to make predictions regarding the expected response function of a liquid scintillator, consider the case of a scintillator with a very large active volume. It is obvious that multiple elastic scattering of the incident neutron on protons should make a large contribution to the total light output response functions. Observing that the response function for multiple np scattering can be calculated by multiple convolutions of the same probability density function, although each weighted by a different factor which depends on the scattered neutron energy before each collision, and using the central limit theorem 1 then it is possible to predict that expected response function should approximate a normal probability density function. This is indeed the case as shown in figure 17 where the detector response function of a 12.7-by-12.7 cm cylindrical EJ-309 liquid scintillator measured experimentally [11] is compared with a normal distribution.
Having now described, albeit qualitatively, the building blocks of the light response function of a liquid scintillator in the simplest case, it is now possible to tackle its interpretation in presence of those effects neglected so far and in more realistic scenarios such as those encountered in fusion devices.  [11]. The experimental data have been acquired with an acquisition threshold which explain the reason why the response function goes to zero at low light output.

Response functions in realistic scenarios
In this section additional aspects affecting the response function of liquid scintillators that have been neglected in the previous sections are discussed. These can be divided in two categoties: those aspects that are intrinsic to the scintillator properties which are discussed in sections 7.1 and 7.2 and those which depend on the incident neutron energy spectrum, discussed in section 7.3.

Detector properties affecting the light output response function
The response function discussed so far has been limited to the case of 2.45 MeV mono-energetic incident neutrons for which the predominant interaction mechanism is elastic scattering on hydrogen and carbon nuclei. For neutrons of higher energies the following competing processes occur: inelastic nC scattering (for E n > 4.8 MeV) and the nuclear reactions 12 C(n, α) 9 Be (for E n > 7.4 MeV) and 12 C(n, n ′ ) 9 Be * (for E n > 10 MeV). At these high energies, nuclear reactions induced by the neutron interacting with the detector material surrounding the liquid scintillator cell are also possible. Inelastic scattering results in the scattered neutron and recoil particles to have a lower energy than in the case of elastic scattering. The light output from these recoils particles have smaller amplitude and will be distributed continuously from zero to the maximum light output possible. The light output from α-s of few MeV is much higher than that of recoil carbons of similar energy but still lower by a factor of about 10 compared to the one for recoils protons (see figure 1). The contribution from α-s appear as an additional "edge" at the low end part of the light output response function superimposed to the "box-like" response for protons. Indicating with L p the light output resulting from a "head-on" collision of a high energy neutron with a proton and with L α the light output produced by an α particle then L α /L p ≈ 0.1 for E n,0 = 15 MeV [10]. Due to the complexity of these competing events as well as the strong angular dependence of their cross section, full Monte Carlo simulations are indispensable for a correct modelling of experimental response functions at neutron energies above few MeVs.
The transport of the scintillation light from its point of generation inside the active detector volume to the PMT becomes important for large detectors or detectors with one dimension much larger than the other. This is due to the fact that the scintillator itself is a self-absorbing medium and therefore the same deposited energy will result in two light pulses of different amplitude depending on the location of the recoil particle generation with respect to the PMT. The effect of light transport and self-attenuation is usually negligible for small volume liquid scintillators (a few centimetres). For larger scintillators, self-attenuation results in light pulses being distributed at lower light output than expected. For example, in the case of incident mono-energetic neutrons onto a long cylindrical scintillator with a PMT at one end of the cylinder, self-absorption will result in the "smearing" of the "head-on" collision edge towards lower amplitudes thus worsening the detector energy resolution. If liquid scintillators are to be used for neutron spectroscopy it is therefore essential that they are of limited dimensions. Codes such GEANT4 [15] and MCNP-PoliMI [14] are capable of treating properly the transport of photons and should be used for scintillators of shape more complex than the standard cylindrical one and for large volumes.
Two additional quantities that affect the light response function are the light output as a function of the deposited energy by different recoil particles and the energy resolution function. The light output function depends by the size and geometry of the scintillator and affects the relative position of the different features in the light output pulse height spectrum and in particular the position of the "edge" corresponding to the deposition of the full neutron energy in a single elastic scattering with a proton ("head-on" collision). The energy resolution dependency on the neutron energy results in the smearing of the response function and it is mostly predominant at the location of the "head-on edge". The energy resolution is affected in turn by the geometry and size of the scintillator and typically ranges between 5 and 30 % in the neutron energy range 0 -3 MeV. For a discussion on the different methods for measuring the light output and energy resolution functions the reader is referred to [10,11] the references therein. Since the "head-on edge" region is the one generated by neutrons that have not undergone any collision, the spectral information on the neutron source is mainly contained in this region and therefore the determination of the light output and energy resolution functions is very important especially at the energy of interest (for example at 2.45 and 14.1 MeV for fusion plasmas).

Detector efficiency
The efficiency of a liquid scintillator to a mono-energetic neutron with energy E n,0 for single elastic scattering on hydrogen is given by [16]: where d is the detector thickness, n H and n C are the number density of hydrogen and carbon atoms respectively, and σ H and σ C are the neutron scattering cross-section on hydrogen and carbon. The exponential terms represents the fraction of neutrons of energy E n,0 passing through a scintillator of thickness d (i.e. the uncollided neutrons) and therfore the terms within curly brackets represents the fraction of neutrons interacting in the scintillation material. The detector efficiency can be calculated using full Monte Carlo codes as the integral of the pulse height spectrum: the efficiency of a 1.5 cm thick scintillator as a function of the incident neutron energy is shown in figure 18. Since in reality, pulse height spectra are measured with an acquisition threshold, as shown in figure 17, the lower limit of the integration must be accurately determined. This is typically done by calibrating the detector light output against the energy deposited by recoil electron using standard γ-rays calibration sources (such as 22 Na, 137 Cs and 207 Bi) and then converting the experimental light output threshold into a deposited energy threshold using the corresponding light output function [18,19]. The light output so calibrated is given in units of MeV electron-equivalent (MeVee) as shown in figure 17. In this case, the light output threshold is approximately 0.1 MeVee resulting in a deposited energy threshold of approximately 0.8 MeVee for recoil protons (see figure 1). In particular, the efficiency and the scintillator cross-sectional area determine the number of counts per seconds that are measured for a given incident neutron flux. As shown in figure 18, the probability of a neutron to interact with the scintillator increases as its energy decreases. As a results, neutrons that have undergone (multiple) scattering before reaching the detector can make a significant contribution to the pulse height spectrum. Neutrons that have collided first with the scintillation material will have an increased probability of making multiple collisions especially for thick detectors as the detection efficiency increases. It is clear then that the dependency of the efficiency on the neutron energy has a non-trivial effect on the light response function both in shape and amplitude. Full Monte Carlo codes ensure that all these effects are included in the response function.

Neutron energy spectra in fusion devices and forward modelling
In fusion devices, liquid scintillators are typically located inside thick neutron shield and view the neutron source (the plasma) via long collimators and the energy spectrum of the neutrons at the detector are anything but mono-energetic. A typical fusion neutron energy spectrum contains spectral components from i) DD and DT thermal emission (approximately gaussian with a width which depends on the plasma ion temperature), ii) fusion reactions involving ions with energies much larger than the thermal ions produced by additional heating and iii) scattered neutrons, that is, neutrons that before their first interaction with the scintillator have already suffered one or multiple collisions with the material surrounding it (from the fusion device itself, to the neutron diagnostic in which the detector is installed). In such a situation, the interpretation of the experimental response function requires the accurate estimate of the neutron flux and energy spectrum at the detector location. In turns this requires accurate plasma physics models for the calculation of the DD and DT neutron sources and of the neutron emission along the line of sight by which the neutron source is seen by the detector to provide the direct (not collided) neutron flux and energy spectrum at the detector. In addition, Monte Carlo neutron transport simulations are required to estimate the scattered neutron flux and its energy spectrum given a specific neutron source and line of sight. The resulting neutron energy spectrum at the detector will contains all these spectral components at different levels.
An example of the different components of the neutron spectrum at a detector location in a DD plasma with neutral beam heating is shown in panel (a) of figure 19 which can be thought of as the neutron emissivity multiplied by the detector area and the solid angle. The response function for each incident neutron energy is then calculated resulting in the response function matrix in which, as shown in panel (b), the columns contain the light output response function at for a particular neutron incident energy. The response functions for incident neutron with energies 2.45 and 4.00 MeV are shown in panels (c) and (d): each response function is then folded with the energy resolution function. The neutron energy spectrum can then be thought as a set of weights that multiply the response function matrix for each incident neutron energy resulting in the weighted response function matrix shown in panel (e). Adding together the contribution from all neutron energies present in the neutron energy spectrum is then equivalent to the column-wise sum of the elements of the weighted response function matrix: the result in shown in panel (f) together with experimentally measured light output pulse height spectrum integrated over a given time interval and multiplied by the detector efficiency. The comparison shown in panel (f) highlights one final aspect of the forward modelling, that is, the conversion of the matrix of normalized response functions multiplied by the absolute neutron energy spectrum into absolute counts. The geometry of the line of sight is already included in the neutron energy spectrum shown in panel (a) so this is achieved by multiplying the neutron energy spectrum by the integration time used to calculate the measured pulse height spectrum in panel (f) and by the detector efficiency.
In this case, the agreement is not perfect. Several reasons could be responsible for this. One possibility is that the original neutron spectrum might not have been properly estimated because one or more plasma parameters that are used in the modelling codes are measured or estimated incorrectly. Among such quantities are the fuel ion and the fast ion (from the additional heating) distribution functions in space and energy, the plasma temperature and density. If one assumes that the detector response function is well known, then such discrepancy can be used to improve the models used for the description of the plasma. This method of proceeding in the modelling of the measured light pulse height spectra is therefore called "forward modelling" since, as the name implies, it starts from the source of the neutrons and then propagates the resulting neutron field taking into account all the aspects affecting the neutron transport from the source to the detector (including realistic geometries and material compositions of the fusion device and of the diagnostic) and then converts the neutron field at the detector into the light output pulse height spectrum via the response function matrix. It is clear therefore how crucial it is that the liquid scintillator response function is determined as accurately as possible if information on the plasma itself is required from the neutron measurements from liquid scintillators.

Summary
The interpretation of the neutron response function of liquid scintillators can be rather subtle as discussed in this tutorial even in the simplest case imaginable of a mono-energetic neutron source. This case, which can be approximately obtained in neutron calibration facilities, however is very useful to understand the different processes that contributes to the neutron light response function which can then easily generalized to the case of fusion neutron sources. Response functions can be easily obtained with Monte Carlo codes as discussed here but their use and the interpretation of the calculated results can be quite complex and as shown in this tutorial not strictly necessary for a qualitative understanding of the response function. Monte Carlo codes are however indispensable whenever a quantitative comparison with experimental measurements is crucial. This is particularly true in the case where one wants to infer absolute quantities from measurements based on the response of liquid scintillators from fusion plasmas such for example the absolute neutron yield (which requires an accurate determination of the detector efficiency) or the underlying fuel ions distribution functions from the measured neutron energy spectrum. In such cases, forward modelling is a robust although quite complex method that enables an understanding of the underlying plasma properties and therefore can provide useful feedback information on the theory and tools used to model the physics of fusion plasmas. This, in turns, requires a very good overall characterization of the response function of liquid scintillators. It is hoped that this tutorial provides an useful introduction to the understanding of the response functions obtained by full Monte Carlo codes.