Bayesian Inference as a Tool to Optimize Spectral Acquisition in Scattering Experiments

Abstract

Nowadays, an increasing number of scattering measurements rely on the use of large-scale research facilities, which is usually granted after highly competitive peer-reviewing and typically for short-time lapses. The optimal use of the allocated time requires rigorous estimates on the reliability of the data analysis, as inferred from the limited statistical accuracy of the measurement. Bayesian inference approaches can significantly help this endeavor by providing investigators with much-needed guidance under challenging decisions on experimental time management. We propose here a method based on the real-time data analysis of running experiments, which fully exploits the core strengths of Bayes theorem. The procedure is implemented in sequential steps in which the spectral measurement is adjourned by summing to it successive acquisition runs, and the spectral modeling is upgraded accordingly. At each stage, the statistical accuracy of the measurement improves, and a more grounded joint posterior distribution is drawn and used as a prior in the subsequent data acquisition stage. The gradual reduction in the model parameters’ uncertainty down to the targets set a priori by experimenters provides a quantitative “success criterion,” which helps prevent oversampling during acquisition. A similar “on the fly” data modeling, might substantially change the way large-scale facilities operate.

1. Introduction

Nowadays, fundamental and applied research in Condensed Matter Physics relies heavily on the use of large research infrastructures. These include continuous or spallation neutron sources or X-rays synchrotron facilities present today worldwide. Often, sources of neutrons and X-rays are indeed found in the same geographical place, given the recognized complementarity of these two powerful spectroscopic techniques for the study of matter. This is the case of the European Photon and Neutron (EPN) campus in Grenoble, France, which hosts the Institut Laue-Langevin (ILL) [1], and the European Synchrotron Research Facility (ESRF) [2]. Similarly, both neutron and X-ray facilities are hosted by the Rutherford Appleton Laboratory in Oxfordshire, UK, (ISIS and Diamond, respectively) [3, 4], by the Paul Sherrer Institut in Villigen, Swiss (SINQ and SLS) [5, 6], and soon by the city of Lund, Sweden (ESS and Max IV Laboratory, respectively) [7, 8]. Large-scale facilities are accessible to scientists for beam-time allocation through a highly competitive proposal selection carried out by expert panels through peer-review processes. Based on this peer-review outcome, the number of days (or even hours) assigned to an experiment is thoroughly pondered. It readily appears how critical is to establish an optimal experimental strategy enabling to gather the most informative and precise data out of an approved measurement. For this purpose, one needs to evaluate not only the ideal number of samples and related physical and chemical conditions, but frequently (if not always, in neutron scattering experiments) the time needed for a certain number of ancillary measurements that are mandatory to achieve a clean set of data. These include accurate measurements of the resolution function, the background signal, and spurious intensity effects, in which the raw measurement needs to be precisely corrected for. Therefore, an optimal use of the beam time assigned to an experiment would greatly benefit from a quantitative criterion to take sensible decisions during the measurement. Here, we propose a simple method to achieve such a criterion based on Bayesian statistics and its inferential capabilities [9, 10, 11]. In Section 2, we briefly describe an inelastic neutron or X-ray measurement and the main concerns rising when deciding its duration. In Section 3, we focus on the output of a Brillouin Neutron scattering (BNS) experiment: the spectrum of density fluctuations of a system; in particular, we show how one can use a Bayesian approach to model this observable. In the same section, we recall a fundamental property of the Bayes theorem that makes it suited to a recursive use for data analysis purposes. To demonstrate the potentialities of this approach, we reproduce the results of a typical BNS measurement by generating simulated experimental spectra. We then summarize the results of an on-the-fly data modeling of these spectra, which enables us to draw a joint posterior distribution for the adopted model parameters eventually guiding the decision on when conveniently stop a spectral acquisition.

Such a running analysis should establish the premises for developing a Measurement Integration Time Optimizer (MITO), a computational tool to assist scattering experiments in large-scale research facilities. In Section 4, we will shortly mention aspects of the approach described which deserves attention or caution; finally, in Section 5, conclusions and possible perspectives are outlined.

2. Neutron and X-ray scattering measurements

The main outcome of neutron and X-ray scattering measurements is the rate $\dot{N}$ of neutrons or photons scattered at an angle $2\mathrm{\Theta }$ with energy changed by an amount $\mathrm{\hslash }\omega$, and ultimately captured by an array of detectors intercepting a finite solid angle. Aside of instrumental factors such as flux, detector efficiency, or detector sensitive area, the intensities recorded by the detectors depend on the physicochemical properties of the sample via its spectrum of density fluctuation, $S\left(Q,E\right)$ [12], which conveys insights on the structure (positions) and the dynamics (movements) of the atoms in the sample.

Obviously, the longer the detector counting, the more accurate the spectral shape determination. In fact, the number of neutron (x-ray) counted within an integration time $t$, $\dot{N}t$, obeys to a Poisson distribution, its standard deviation thus being $\sqrt{N}$. As the integration time $t$ increases, the counting statistics improves as the relative experimental errors ($\sim 1/\sqrt{N}\equiv 1/\sqrt{\dot{N}t}$) decreases. Hence, the chance to detect interesting details of the spectral shape ultimately depends, of course, on the sample properties, but also on the accuracy of the intensity measurement. A difficulty to be faced in typical INS spectral acquisitions is that the measurement might be terminated prematurely, that is, before providing the information sought for. This possibility appears especially penalizing if the counting statistics achieved is not accurate enough, the spectral features not well-defined, or the signal sought for very weak. Conversely, data can also be integrated longer than strictly needed to capture the effect under scrutiny. In this case, further prolonging the counting would not complement the insight of the measurement significantly, and, beyond some time lapse, would not even improve its quality. Even worse, it could jeopardize the ultimate success of the experiment due to the time waste, which could prevent the accomplishment of the full experimental plan.

Without digging into computational details, here we outline a strategy to support experiments with a Bayesian protocol providing useful assistance in the measurement’s planning and decision making. This will help investigators to determine when the integration time of a spectral acquisition can be safely stopped, either because all useful information was gathered, or because the predetermined target established for relative uncertainties was reached. For the sake of simplicity, we will focus on an exemplary inelastic neutrons scattering (INS) case, with the implicit assumption that the method can be safely extended to X-rays scattering (IXS), in fact being generally valid for any spectroscopy measurement. We stress that the case we are considering is very likely also the most demanding in terms of computational effort for reasons that will be briefly illustrated later in this chapter.

3. Inelastic neutron scattering

As mentioned, the general aim of an INS measurement is to measure the spectrum of density fluctuations $S\left(Q,E\right)$, which conveys insights on positions and movements of the atoms in a sample. Oversimplifying, depending on the spectrometer we use to determine it, we can have access to different aspects of $S\left(Q,E\right)$, either relating to collective movements of atoms (e.g., acoustic waves, structural relaxation processes), single-particle ones (e.g., translational diffusion, rotations, librations …), or both. To measure $S\left(Q,E\right)$ of a given system one can use two different types of neutron spectrometers: triple axis spectrometers (TAS) and time-of-flight (TOF) ones.

Here, we assume to execute measurements with a TOF spectrometer [13] where $S\left(Q,E\right)$ surfaces are sampled, ideally, for each $Q$ and $E$ values simultaneously [14]. The rate of neutrons scattered at the different scattering angles $2\mathrm{\Theta }$ (see Appendix A) and impinging on the sensitive area of the detector after a time (of flight) $t$ defines the time-dependent intensity function $I\left(2\mathrm{\Theta },t\right)$. The latter is converted into an intensity $I\left(Q,E\right)$, which is a function of the momentum, $Q$, and the energy, $E=\mathrm{\hslash }\omega$, exchanged between sample and probe particles, with $\mathrm{\hslash }$ and $\omega$ being the reduced Planck constant and the exchanged frequency, respectively. In Appendix A, a sketch of the BRISP spectrometer and few hints about the principles of the TOF technique are shortly recalled.

Aside from instrumental effects such as energy resolution and signal background, $I\left(Q,E\right)$ is proportional to $S\left(Q,E\right)$, which is the physical variable, INS (and IXS) investigators usually seek for. To sample and gather this intensity function with adequate counting statistics, providing us the needed information, a certain acquisition time is required, which depends on the characteristics of the instruments (incident neutron flux, detector efficiency, resolution …) and on the scattering properties of the sample. These are embodied in its double differential cross section $d^2\sigma /d\mathrm{\Omega }\mathit{dE}$ [12, 14] defined as the number of neutrons deviated in a second into the small solid angle $\mathrm{\Delta \Omega }$ subtended by a detector along the $2\mathrm{\Theta }$ direction, with final energy included in the interval between $E$ and $E+\mathrm{\Delta }E$ [14]. More explicitly:

$$\dot{N}=\frac{d^2\sigma }{d\mathrm{\Omega }\mathit{dE}}J\mathrm{\Delta \Omega \Delta }E.$$E1

where $J$ represents the incident flux of neutrons.

To summarize, the image of the $S\left(Q,E\right)$ is being built up as the measurement runs, and the larger the acquisition time, the more precise is the $S\left(Q,E\right)$ rendering. To avoid data loss of the entire measurement in case of instrument failure, the measurement is usually split into different sub-runs which, for the sake of simplicity we will hereafter assume to have the same acquisition time.

3.1 A Brillouin neutron scattering measurement

Once the $S\left(Q,E\right)$ surface is measured, different constant $Q$ cuts of it are determined by interpolation. As an example of a typical INS measurement outcome, in Figure 1 we show the spectrum of liquid silver at $Q=6{\mathit{nm}}^{-1}$ measured with the Brillouin spectrometer BRISP [15] at the High Flux Reactor of the Institut Laue Langevin (Grenoble, France) [16, 17]. The spectral intensity has the typical shape of the spectrum from a disordered samples, which consists of a central peak broadened around the elastic energy $E=0$ sided by a pair of inelastic peaks shifted from the elastic energy by an amount of $\pm \mathrm{\hslash }{\omega }_s$, which defines the energy of the excitation. The peaks sitting at positive and negative energy are respectively the well-known Stokes and Anti-Stokes lines of a Brillouin spectrum [18].

Figure 1.

Dynamic structure factor of liquid silver measured on the Brillouin spectrometer BRISP at a momentum transfer $Q$ = 6 nm${}^{-1}$ and an incident energy $E_0$ = 83.9 meV. The experimental data points are affected by resolution broadening. The dashed dot line is a fit obtained from an oversimplified model: A central Lorentzian line is summed to a damping harmonic oscillators (DHO) function to describe the inelastic excitations [16].

3.2 Modeling through a Bayesian approach

Let us assume that measured data can be described by a chosen model specified by the vector $\mathrm{\Theta }=\left({\theta }_1,{\theta }_2,\cdots ,{\theta }_m\right)$, whose generic component ${\theta }_m$ is a model parameter. For the sake of generality, we include the possibility that some $\mathrm{\Theta }$ components, instead of being the parameter of a well-identified model, designate instead one model option among several competitive ones whose reliability is to be concurrently tested [10, 19]. The vector $y=\left(y_1,y_2,\cdots ,y_n\right)$ indicates instead the measured dataset with $n$ being the sample size that is the number of data points. With this notation in mind, we can express the Bayes theorem [20] as follows:

$$P\left(\mathrm{\Theta }|y\right)=\frac{P\left(y|\mathrm{\Theta }\right)P\left(\mathrm{\Theta }\right)}{P\left(y\right)},$$E2

where $P\left(\mathrm{\Theta }|y\right)$ is the posterior distribution of the parameters built according to the experimental outcome, $P\left(\mathrm{\Theta }\right)$ is the prior distribution (or simply prior) of the parameters, that is, that in one’s hands before any data measurement, $P\left(y|\mathrm{\Theta }\right)$ is the likelihood of the data, that is, the probability of observing the data conditional on a certain parameter vector, and $P\left(y\right)$ is the marginal probability of the data, which plays the role of normalizing constant, so that Eq. (2) has a unit integral over the variable $\mathrm{\Theta }$. We stress that the prior probability includes all our initial knowledge (or ignorance) and can be more or less informative depending on the preliminary insight we got on the problem at hand. Bayes’ theorem is therefore a prescription on how to learn from experience, insofar as it gives a golden rule to update one’s beliefs in light of the accrued data.

Now let us imagine to have achieved a portrayal of $S\left(Q,E\right)$ of a given sample from a measurement run of a certain duration $t_1$. We can ideally try to fit this first rough $S\left(Q,E\right)$. This will provide a first joint multi-dimensional posterior distribution of the parameter vector $\mathrm{\Theta }$, which likewise improves our knowledge of model parameters with respect to the prior we started with. It is meaningful to think about this posterior as an updated prior to feed back into the Bayes theorem, as we keep gathering new data in the experiment.

Unfortunately, the posterior distribution has no explicit analytical expression, thus being of hardly any use to feed back in the Bayes theorem again. We then measure a second run, which, with no loss of generality, can be assumed for simplicity of the same duration $t_2=t_1$ of the first one. New data can be certainly add to the old ones to get a new, more accurate, dataset.

Upon indicating data gathered during the first and second run, respectively, as $y=\left(y_1,y_2,\cdots ,y_n\right)$ and $y^{\prime }=\left(y_1^{\prime },y_2^{\prime },\cdots ,y_n^{\prime }\right)$, we can formally express the posterior distribution of the parameters vector, conditionally on the complete collection of data, as

$$P\left(\mathrm{\Theta }|y,y^{\prime }\right)=\frac{P\left(y^{\prime }|\mathrm{\Theta },y\right)P\left(\mathrm{\Theta }|y\right)}{P\left(y^{\prime }|y\right)}$$E3

which is a mere formulation of the Bayes theorem. We observe how the prior we have now is just the posterior distribution for the vector parameter $\mathrm{\Theta }$ having already observed the dataset $y$.

The datasets $y^{\prime }$ and $y$ being independent, we have that $P\left(y^{\prime }|\mathrm{\Theta },y\right)=P\left(y^{\prime }|\mathrm{\Theta }\right)$ and $P\left(y^{\prime }|y\right)=P\left(y^{\prime }\right)$. On the other hand, we can apply again the Bayes theorem to get $P\left(\mathrm{\Theta }|y\right)=P\left(y|\mathrm{\Theta }\right)P\left(\mathrm{\Theta }\right)/P\left(y\right)$. Doing the substitutions, Eq. 3 becomes:

$$P\left(\mathrm{\Theta }|y,y^{\prime }\right)=\frac{P\left(y^{\prime }|\mathrm{\Theta }\right)P\left(y|\mathrm{\Theta }\right)P\left(\mathrm{\Theta }\right)}{P\left(y^{\prime }\right)P\left(y\right)}$$E4

Once again, $y^{\prime }$ and $y$ being independent, we have $P\left(y^{\prime }|\mathrm{\Theta }\right)P\left(y|\mathrm{\Theta }\right)=P\left(y^{\prime },y|\mathrm{\Theta }\right)$ and $P\left(y^{\prime }\right)P\left(y\right)=P\left(y^{\prime },y\right)$ for the property of the joint probability of independent variables. Eq. (4) becomes:

$$P\left(\mathrm{\Theta }|y,y^{\prime }\right)=\frac{P\left(y,y^{\prime }|\mathrm{\Theta }\right)P\left(\mathrm{\Theta }\right)}{P\left(y,y^{\prime }\right)}$$E5

We finally observe that the posterior probability for the vector $\mathrm{\Theta }$ given the datasets $y^{\prime }$ and $y$ can be obtained viaEq. (3) provided the posterior we derived after the first measurement is used as a new prior. This is equivalent to using as a prior the one we started with, yet multiplied for the likelihood pertinent to (inclusive of) all data collected thus far.

We can thus apply Eq. (5) in a recursive fashion to analyze on the fly neutron scattering data as we collect them. We would like to determine the most appropriate total acquisition time based on solid statistical arguments and with the prospect of inferring something about the quality of the data collected. The ultimate goal would be to have the possibility of ending the acquisition when the maximum level of information that can be obtained from the measurement has already been reached. Further prolonging the acquisition would not bring any extra relevant information. Certainly deciding when the measurement can reasonably be interrupted is at the discretion of the experimenter who may still want a specific precision from the measurement.

3.3 Simulation of a Brillouin neutron scattering experiment

Let us imagine to perform a BNS measurement. The instrument will acquire scattering intensity for a certain time, splitting such an acquisition in separate runs of the same duration. We can focus now on the spectrum corresponding to a constant $Q$ cut of the $S\left(Q,E\right)$ we are measuring. To visualize this, we can generate simulated experimental data as they were actually measured. Table 1 provides the parameters setting used in the simulation study. These otherwise arbitrary parameter values are chosen to reproduce a typical spectrum collected in a neutron or X-ray scattering measurement on an amorphous system. More specifically, we have tuned the parameters so to obtain barely resolved inelastic contributions to the spectrum, which is a typical problem faced in routine measurements. Indeed, because of the limited instrument resolution and finite counting statistics, whenever the excitation features are not sufficiently sharp, that is, adequately separated and broadened by small damping, they can blend to each other and partially merge into the dominant central peak, which makes their detection challenging. In these cases, the limited acquisition time of a single measurement become an even more constraining factor. Therefore, the chosen parameter set is suitable to mimic a typical scenario encountered in scattering studies on disordered systems. In Figure 2, we show one of these datasets randomly generated using the following simple model:

Parameter	Value
$A_e$	33
${\mathrm{\Omega }}_1$	4
${\mathrm{\Omega }}_2$	8
${\mathrm{\Gamma }}_1$	0.8
${\mathrm{\Gamma }}_2$	1
$A_1$	10
$A_2$	5

Table 1.

Absolutes parameter values for the model from which the simulated datasets were drawn.

Figure 2.

Generated spectrum as drawn from the model in Eqs. (6) and (7) at a $Q$ value of 5 nm${}^{-1}$. This spectrum simulates the data as they could appear after a very short acquisition run. The plotted quantity is in fact the scattered intensity, to which the dynamic structure factor is proportional.

$$\tilde{S}\left(Q,E\right)=R\left(E\right)\otimes S\left(Q,E\right)$$E6

where $\otimes$ indicates the convolution product, $R\left(E\right)=\frac{1}{\sqrt{2\pi }\zeta }exp\left(-\frac{E^2}{2{\zeta }^2}\right)$ is the instrument resolution function, assumed to have a zero-centered Gaussian profile with ${\zeta }^2$ variance, which in the present case gives a $\mathit{\text{FWHM}}$ = 3.1 meV. Once again, we suppose to use the BRISP spectrometer with an incident energy $E_0$ = 80 meV, as achieved by using the (004) reflection from a Pyrolytic Graphite monochromator [15]. The dynamic structure factor is here approximated as:

$$S\left(Q,E\right)=A_e\left(Q\right)\delta \left(E\right)+\left[n\left(E\right)+1\right]\frac{E}{k_{B}T}\left\{\sum _{k=1}^2\frac{2}{\pi }{A}_k\left(Q\right){\mathit{DHO}}_k(Q,E)\right\}$$E7

where $\delta \left(E\right)$ is the Dirac Delta function describing the elastic response of the system modulated by an intensity factor $A_e\left(Q\right)$, ${\mathit{DHO}}_k$ are $k$ inelastic contributions to the spectrum described by Damped Harmonic Oscillator (DHO) functions [21], $n\left(E\right)={\left(e^{E/k_{B}T}-1\right)}^{-1}$ is the Bose factor expressing the detailed balance condition, $k_B$ is the Boltzmann constant, and we have chosen the temperature $T$=1337 $K$. Finally, the simulated experimental data points are corrupted by an additive random fluctuation $\epsilon \left(Q,E\right)$:

$$y\left(Q,E\right)=\tilde{S}\left(Q,E\right)+\epsilon \left(Q,E\right),$$E8

with $\epsilon \left(Q,E\right)\sim \mathcal{N}\left(0,{\sigma }^2\tilde{S}\left(Q,E\right)\right),$ for any $Q$ and $E$, where the symblol $\sim$ means “distributed according to,” $\mathcal{N}$ denotes the Gaussian distribution and ${\sigma }^2$ is a constant factor to be estimated on the experimental data.

We can then generate as many experimental runs as we wish and sum them as usually done. When a scattering measurement is actually performed in the large-scale facilities above mentioned, the beam time granted to researchers is probably the most critical requisite for a successful experiment. The number of equal duration measurement runs on a given sample provides the total time allotted for that sample. In Figure 3, we show data as they result from the sum of 20 runs of identical integration time to qualitatively visualize the improvement in data precision and spectral shape definition that can be achieved by enhancing the counting statistics through a factor 20 increase of the acquisition time.

Figure 3.

Sum of 20 generated spectra as drawn from the model in Eqs. (6) and (7).

We will try now to fit our experimental data using a Bayesian Markov Chain Monte Carlo (MCMC) [22] algorithm equipped with a Reversible Jump option (RJ) [23], as explained in detail in Refs. [10, 11]. This algorithm allows to draw values from a distribution which is only known up to a normalization constant and thus to simulate the joint posterior distribution of the parameter vector of the model, $\mathrm{\Theta }$, as defined in Eq. 7. The analytical evaluation of the normalization constant is in fact usually really hard if not impossible at all. We again stress that in this simplified model the number $k$ of inelastic components contributing to the spectrum is in itself a free model parameter to be estimated conditionally on available data. Notice that the RJ option allows the MCMC algorithm to explore models with different numbers $k$ of inelastic components with $k=1\dots k_{\mathit{max}}$, $k_{\mathit{max}}$ being the maximum number of excitations allowed. As a first step, the first measurement run is best fitted by the model and the first-level posterior distribution of the model parameters vector is thus obtained. Once the corresponding marginal distribution $P\left(k|y\right)$ has been computed, the best-fit model conditional to the first experiment run can be determined. After a second run of data $y^{\prime }$ is available, this is added to the previous one and the MCMC algorithm is used on the complete data to obtain a new joint posterior $P\left(\mathrm{\Theta }|y,y^{\prime }\right)$. Parameter estimates from the first run can be used as starting values to speed up convergence. The process repeats itself as new runs become available.

Here below (Figures 4–6) we show the posterior distribution for some of the model parameters when the data are analyzed considering 5, 10, 20, 40 and finally 60 runs. We emphasize again that at each step we are applying the Bayes theorem, feeding back the posterior we obtained at a previous step as a prior (new knowledge about the data) for the following step. The likelihood is enriching itself more and more as long as we acquire new data. In the present example, the algorithm finds as best model the one with two inelastic modes as it should be desirable since the generation model is the one in Eq. (7). In fact because of the random error added to simulate a real experimental dataset, especially if the inelastic modes are chosen to have frequencies close to each other and/or large damping, it is not straightforward to find the number of modes predicted by the model whatever is the amount of data considered.

Figure 4.

Posterior distribution of the two excitation frequencies as estimated from the Bayesian analysis after 5 (green), 10 (purple), 20 (yellow), 40 (brown), and 60 (blue) experimental runs.

Figure 5.

Posterior distribution of the dampings of the two excitations as estimated from the Bayesian analysis after 5, 10, 20, 40, and 60 experimental runs.

Figure 6.

Posterior distribution of the amplitudes of the two excitations estimated from the Bayesian analysis after 5, 10, 20, 40, and 60 experimental runs.

Figure 4 clearly demonstrates that the posterior distribution for the inelastic shifts ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ of the spectrum get sharper upon increasing the number of runs considered in the analysis. The distribution becomes better shaped, peaked, and symmetric providing, of course, a better estimate of the model parameters. Still, if we take the mean of the posterior distribution for ${\mathrm{\Omega }}_1$ (or for ${\mathrm{\Omega }}_2$) after 40 runs and we compare it with the one obtained after 60 runs, the difference between these two means is about 1$\%$; when comparing instead the mean after 5 runs with the mean after 60 runs, the difference amounts to less than 3$\%$ which very likely is already smaller than experimental uncertainties typically reported in dispersion curves displayed by scientific papers. Similar considerations hold for the other two parameters defining the damped harmonic oscillators, namely the peak amplitudes $A_1$ and $A_2$ and the dampings ${\mathrm{\Gamma }}_1$ and ${\mathrm{\Gamma }}_2$ (Figures 5 and 6) [24].

In Figures 7 and 8, the best fit after 5 runs and 60 runs sum spectra is shown along with the estimated DHO components and the central elastic contribution.

Figure 7.

Simulated spectra (blue dots) at a $Q$ value of 5 nm${}^{-1}$ as obtained summing 5 runs drawn by the model in Eq. (7). The best-fit model (red curve) to the drawn data and the estimated DHO components (black and green line) are also shown. The dash dot magenta line is the estimated elastic central component.

Figure 8.

Simulated spectra (blue dots) at a $Q$ value of 5 nm${}^{-1}$ as obtained summing 60 runs drawn by the model in Eq. (7). The best-fit model (red curve) to the drawn data and the estimated DHO components (black and green line) are also shown. The dash dot magenta line is the estimated elastic central component.

In Figure 9, we report the posterior distribution function for the number of the detected inelastic modes as a function of the number of runs considered in the analysis. It is evident that as the experimental evidence becomes more precise, the probabilities associated with a higher number of modes progressively vanish.

Figure 9.

Top panel: Posterior distribution function for the number $k$ of inelastic modes detected in the simulated Brillouin spectrum as a function of acquisition time. From top to bottom, the results after 5, 10, 20, 40, and 60 experimental runs. Bottom panel: As in the top panel but at the very top of the figure the posterior of $k$ after only 1 run also is shown. In the insets, two different priors $P\left(k\right)$ are shown. In the top panel a uniform prior for $k$ is plotted. In the bottom panel, a modified (see text) binomial prior distribution has been chosen for comparison.

To conclude we briefly draw the reader’s attention on the necessity to assess convergence of the MCMC algorithm before using its output for inferential purposes. In literature, a great deal of effort has been spent in developing convergence diagnostic tools for MCMC. Some of these tools are specifically intended to check convergence of the Markov chain to the stationary distribution, or to check for convergence of summary statistics, such as sample means, to the corresponding theoretical quantities. For a recent review of the subject, see Ref. [25]. Although many convergence criteria and stopping rules with sound theoretical foundation have been proposed, in practice MCMC users often decide convergence by applying empirical diagnostic tools, in particular graphical methods. The most common graphical convergence diagnostic method is the trace plot, which is a time series plot showing the values of the model parameters at each sweep against the sweep numbers. The trace plot enables to visualize the capability of the Markov chain in exploring the parameter space. For example, the presence of flat bits reveals that the MCMC algorithm gets stuck in some part of the parameter space and is a symptom of slow convergence. This happens when too many proposals are rejected consecutively. On the other hand, when proposals are too easily accepted, the algorithm may move slowly not exploring the parameter space in an efficient way. In this case, the trace plots would show visible trends or changes in spread, implying that stationarity has not been reached yet. Often Bayesian statisticians refer to a “hairy caterpillar” when describing trace plots and what they should look like. In Figure 10, we report trace plots for the excitation frequencies and for the number of inelastic modes in the spectrum. Another helpful graphical method is the running mean plot, which shows parameters’ time-average estimates against the iterations. This line plot should stabilize to a fixed value as iteration increases (Figure 11).

Figure 10.

Left panel: Trace plot for the excitation frequencies in the spectrum obtained summing 20 experimental runs. Right panel: Trace plot for the number $k$ of DHOs for the same data.

Figure 11.

Left panel: ${\mathrm{\Omega }}_{1,2}$ time-average estimates as a function of algorithm sweep. Right panel: Cumulative occupancy fraction for the most visited models.

4. A few caveats and additional remarks

For some neutron scattering techniques, the raw data are not immediately available for a reliable lineshape analysis since collected intensities are affected by many spurious contributions. Indeed, the intensity ultimately detected contains, beside the genuine signal from sample, the unwanted ones coming from the empty cell, multiple scattering, sample environment, background; furthermore, it partly results from misguiding effects such as sample auto-shielding and detector efficiency. The importance or even the presence of such spurious effects changes for different neutron techniques; in this perspective, dedicated considerations and suitable adjustments to the recursive method proposed might be required. Data in the test discussed here are assumed to be already corrected for all these effects, that is, already cleaned from any unwanted contribution. Therefore, this is an ideal case, which might be not always straightforward applicable. Differently a reducing data routine has to be performed in advance. Nevertheless, depending on the technique an effort might be done to recognize some quality parameter to draw out the same conclusions that we can get from the model parameters we have seen here above. It is also true that in other neutron techniques, the possibility to reduce the data rapidly letting them available for an on the run analysis is preventing the problems here aforesaid and this is even more true for IXS.

Overall, the results of the test discussed are not surprising, as an improvement of the statistical accuracy expectedly enhances the precision of the parameters’ determination. Nonetheless, this simple analysis shows how informed decisions about ending or continuing a measurement can be taken based on quantitative grounds. The knowledge of an entire multi-dimensional joint posterior distribution, the evolution of its mode, and overall shape upon increase of the acquisition time could help us to establish not only if data collected are sufficiently integrated but if further counting can enhance the measurement’s insight. To this purpose, we can illustrate briefly an example slightly different from the one proposed before. Let us assume that the dynamic structure factor features two pairs of inelastic peaks (shoulders), for example, not only the one observed in the previous example, but, as normally the case for liquids [16, 26], an additional one which, for some reason, does not emerge clearly from the spectral shape. For instance, at low $Q$‘s, the first pair can be completely submerged by the resolution tails, while at larger $Q$‘s, the paired modes can become highly damped barely standing out from the background, while, at even larger $Q$‘s, they can move out of the energy window covered. Furthermore, if the counting statistics is poor, the marginal posterior distribution for the number of modes can convey ambiguous information and, in some instances, the presence of the second pair of inelastic modes can be overlooked. With an on-line analysis of spectra under collection, one can likely appreciate the possible evolution of such a distribution upon increasing the integration time. For instance, at the early stages of the experiment such a distribution may lead to infer a single pair of inelastic modes, while two (or more) pairs can be inferred as the measurement progresses. However, the incorporation of the Occam’s razor principle [11, 27] in the Bayes theorem represents a safe antidote against the risk of overparametrization, especially when the counting statistics is still poor. This can be sufficient to keep the value of $k$ from exceeding its true value, which in this example is known to be 2. Figure 9 (top panel) illustrates how the posterior distribution of the number of modes $k$ evolves as a function of the integration time. This trend has a straightforward explanation if one considers the gradually improving of statistical accuracy. At the beginning of the measurement, the algorithm could struggle to establish the true value of $k$ assigning not negligible probability to models with a redundant number of modes. As the data are further harvested, the posterior distribution becomes more accurate and the number of modes converges to the most plausible one (i.e., $k=2$). In Figure 9 (bottom panel), we show the evolution of the posterior distribution of $k$ as a function of the measurement acquisition time when a different prior for $k$ is chosen. Since we have simulated experimental runs from a model with two DHOs, in this specific case, we know that the best model to fit the data must have two inelastic modes. Therefore, the chosen prior is a modified binomial distribution which privileges a solution with two inelastic modes. In this case, the prior was:

$$P\left(k\right)=\left(\begin{array}{c}{k}_{\mathit{max}}-1\\ k-1\end{array}\right){\pi }^{k-1}{\left(1-\pi \right)}^{k_{\mathit{max}}-\left(k-1\right)}$$E9

where $k_{\mathit{max}}$ is the maximum number of modes contemplated by the model and we set $\pi =0.3$. With this $\pi$ value, the variable $\left(k-1\right)\sim \mathcal{B}\mathit{\text{in}}\left(k_{\mathit{max}}-1;\pi \right)$ and the different values of the prior are reported in Figure 9 (inset of bottom panel). In this figure, we have included also the results obtained by considering only a single experimental run. It appears that, when we have a firm prior knowledge about the system at hand, the posterior distribution converges to its asymptotic value even faster. In fact after five experimental runs we obtain a probability $P\left(k=2\right)$ already close to 90$\%$. In Figure 12, we show instead the values of the posterior distribution for $k$ attained after a single experimental run. When, as in this case, the counting statistics is really poor the probabilities associated to values of $k$ different from the expected value, that is, the one of the generating model ($k=2$) is not negligible. Finally, the evolution of the posterior for the low-frequency excitation shift strikingly emerges from Figure 13 in which we compare the results obtained either after a single run or after 60 runs.

Figure 12.

Posterior distribution function for the number $k$ of inelastic modes detected in the spectrum after only one experimental run.

Figure 13.

Posterior distribution for the lower excitation frequency after one experimental run and 60 runs.

5. Conclusions and perspectives

This chapter deals with a topic of pivotal interest for scattering experiments at large-scale research infrastructures, as the optimal use of the usually short beamtime allocated for the measurement. The analysis of simulated measurements presented here demonstrates that the assistance of a Bayesian inference protocol can provide a decisive advantage in the decision making and time optimization processes of routine inelastic scattering experiments, and, more in general, of any scattering or diffraction measurement. Specifically, we considered a prototypical neutron scattering study split into shorter acquisition runs; Bayesian inference is used to analyze partial acquisitions obtained by summing an increasing number of individual runs to ultimately guide the investigator in his/her difficult decision on when to stop the beam counting. Such a decision is based upon previously established success criteria, as the achieved evidence for a physical phenomenon affecting the spectral shape, or the met targets in the experimental uncertainties associated with a given lineshape modeling. In this perspective, the development of a dedicated Measurement Integration Time Optimizer protocol could be especially beneficial, as it would provide conventional neutron or X-ray investigations with real-time Bayesian inference assistance. We believe that the availability of a similar on-the-fly data analysis tool can drastically minimize the time wasted in beamtime measurement, also holding the potential for a drastic revision of the beamtime allocation process. In fact, with this novel data analysis tool, decisions on beamtime assignment can be taken on the ground of spectral simulations in which the spectra to be successively measured can be analyzed as obtained with different integration times. We anticipate that these novel inference tool can mark a discontinuity in the workflow of typical scattering experiments at large-scale research facilities.

Conflict of interest

The authors declare no conflict of interest.

Time-of-Flight neutrons instruments are a class of spectrometers which allows measuring inelastic neutron scattering, providing insights into the dynamics of matter. In Figure 14, we show, as an example, BRISP, a direct geometry Brillouin spectrometer once installed in the reactor hall of the High Flux Reactor of ILL. The neutrons scattered by the sample are collected by a highly pixeled detector covering a certain angular range. The scattering angle $2\mathrm{\Theta }$ is defined as the angle between the direct beam axis and the direction of the scattered neutrons (Figure 15). A Fermi chopper device splits the continuous beam coming from the monochromator into 10 $\mathit{\mu s}$ bursts of neutrons and fixes for each burst an initial reference time. The wavevector ${\mathbf{k}}_{\mathbf{0}}$, hence the energy $E_0$ of the neutrons impinging on the sample, are known and so is the time $t_0$ such incident neutrons take to fly from the reference initial time to the sample position. The detector electronics suitably synchronized with the chopper provides us a measure of the total time-of-flight $t_{\mathit{tof}}$ from the reference time to a specific pixel and of course the position and hence the distance traveled by each scattered neutron. If we call $t_1$ the time the scattered neutron takes to fly from an interacting atom in the sample to a specific detector pixel we have that:

Figure 14.

Sketch of the Brillouin spectrometer BRISP. A monochromatized continuous beam, severely collimated is converted in a pulsed beam by a Fermi chopper which labels each pulse with an initial reference time as the starter in a running race. Once the neutrons interact with the sample they are scattered and finally collected by 2D-pixeled detector. Reproduced from ref. [28].

Figure 15.

Sketch of the kinematic scattering triangle. Incident neutrons characterized by a wavevector ${\mathbf{k}}_{\mathbf{0}}$ are scattered by a sample with a wavevector ${\mathbf{k}}_{\mathbf{1}}$. The angle $2\mathrm{\Theta }$ between the incident and the scattered radiation is the scattering angle. The vector $\mathbf{Q}$=${\mathbf{k}}_{\mathbf{0}}$-${\mathbf{k}}_{\mathbf{1}}$ is the transferred momentum in the scattering process.

$$t_{\mathit{tof}}=t_0+t_1=\frac{L_0}{v_0}+\frac{L_1}{v_1}$$E10

where $L_0$, $L_1$ are the distances between the chopper and the sample and between the sample and the detector, respectively, and $v_0$ and $v_1$ are the initial and final velocities of the incident and scattered neutron.

From Eq. (10) it is straightforward to obtain the energy $E_1$ with which the neutron reaches the detector and then the energy $E=\mathrm{\hslash }\omega$ transferred from the probe particle to the sample. It is, in fact:

$$E_1=\frac{1}{2}{\mathit{mv}}_1^2=\frac{1}{2}m{\left(\frac{L_1}{t_{\mathit{tof}}-\frac{L_0}{t_0}}\right)}^2=\frac{1}{2}m{\left(\frac{L_1}{t_{\mathit{tof}}-L_0\sqrt{\frac{m}{2E_0}}}\right)}^2$$E11

and

$$E=E_1-E_0=\mathrm{\hslash }\omega$$E12

Thanks

We would like to thank U. Bafile, E. Guarini, F. Formisano, R. Magli, and M. Maccarini for the very stimulating discussions. The open access fee was covered by Institut Laue Langevin (ILL), Grenoble France.

Abbreviations

BNS	Brillouin neutron scattering
MITO	Measurement Integration Time Optimizer
INS	inelastic neutron scattering
IXS	inelastic X-ray scattering
BRISP	BRillouin spectrometer
DHO	damped harmonic oscillator
MCMC	Markov Chain Monte Carlo
RJ	reversible jump
TOF	Time Of Flight