Optimized signal deduction procedure for the MIEZE spectroscopy technique

A method is reported to determine the phase and amplitude of sinusoidally modulated event rates, binned into four bins per oscillation, based on data generated at the resonant neutron spin-echo spectrometer RESEDA at FRM-II.


Introduction
MIEZE (modulation of intensity with zero effort) spectroscopy is a hybrid technique combining neutron resonance spin-echo and neutron time-of-flight spectroscopy. It is routinely available at the spectrometer RESEDA at the Heinz Maier-Leibnitz Zentrum (Franz & Schrö der, 2015) and BL06 at the J-PARC Materials and Life Science Experimental Facility (Kawabata et al., 2006;Hino et al., 2013). Furthermore, MIEZE is being actively developed at the Reactor Institute Delft, the ISIS neutron source (Geerits et al., 2019) and Oak Ridge National Laboratory (Dadisman et al., 2020). In Fig. 1(a) we present a basic MIEZE setup. It uses neutron spin precession generated by two resonant (neutron) spin flippers (RSF 1 and RSF 2 ), separated by a distance L 1 and operated at individual frequencies (f 1 < f 2 ), to manipulate the spin eigenstates (Gä hler et al., 1996). The resulting interference pattern of the superposition of the spin states corresponds to a sinusoidal intensity as a function of time akin to an optical heterodyne interferometer [see Fig. 1(b)].
The modulation frequency of this intensity is given by twice the difference of the RSF frequencies, f MIEZE ¼ 2ð f 2 À f 1 Þ (Felber et al., 1998;Jochum et al., 2019). In practice, these frequencies are limited at the lower end by the neutron spin flip efficiency generated by the Bloch-Siegert shift to f min = 35 kHz (Bloch & Siegert, 1940). The limitations at the upper end are due to skin and proximity effects in the resonant flippers, as well as parasitic capacities in the resonant circuits, which currently set the maximum RSF frequency to f max = 3.6 MHz (Groitl et al., 2015;Jochum et al., 2020).
In contrast to conventional neutron spin-echo, the quantity measured in MIEZE corresponds to a sinusoidally modulated intensity in time, from which the MIEZE contrast C ¼ I 0 =I mean , with I mean the time average intensity and I 0 the amplitude of the intensity, can be extracted (Gä hler et al., 1992).
In order to increase the time resolution (the Fourier time, MIEZE ), f MIEZE has to be maximized since it is directly proportional to MIEZE via the following relationship: with the neutron mass m n , its velocity v n and the sample-todetector distance L S (cf. Fig. 1). Further details and a description of the MIEZE setup may be found in the literature Oda et al., 2020). The process of data reduction is done in two steps. In the first step, events are registered from electronic signals at the readout of the detector. These events are either accepted as neutron counts and then histogrammed on the fieldprogrammable gate array (FPGA) or rejected on the basis of event length (Schmidt et al., 2010). In the second step, contrast and phase are deduced from this 4D histogrammed data set (pixel Â pixel Â time bin Â foil). This is done by fitting the time bins, in each pixel and on every foil (see below), using a sine function. Here, we present an approach that optimizes the second step of this procedure, requiring less computing power and allowing an increase in maximum achievable Fourier time by a factor of four.
From a practical point of view the detector registers events per oscillation and histograms them according to a certain number of time bins (Kö hli et al., 2016). Thus, for a fixed number of time bins, the length of each time bin is a function of the modulation frequency. The lower limit of the time bin length is given naturally by the temporal resolution of the detector, which is limited by the electron drift time and the clock of the electronics readout of the detector (Kö hli et al., 2016). Hence, to detect signals with fast modulation, it is necessary for the number of time bins to be as low as possible. However, an insufficient number of time bins per oscillation results in a smearing of the recorded oscillation amplitude and a loss in contrast. Therefore, an optimal compromise between the two needs to be achieved.
Without loss of generality and neglecting the average background count rate, the event rate registered by a detector recording signals at discrete intervals in time is given by the integral over a harmonic oscillation with amplitude I 0 and arbitrary phase 0 : where Á ¼ 2=ðNo: of time binsÞ. In this resolution function the sinc function acts as a damping factor, the influence of which becomes smallest for Á ¼ 0, i.e. infinite time bins representing a trivial but trivially impractical solution. Moreover, an infinite number of time bins would require infinite time-stamp accuracy of every event detected. Wrongly binned events decrease the contrast. This reduction scales with the ratio between time-stamp accuracy and time-bin length. From this perspective fewer time bins are preferable as well.
The final measurement quantity extracted from a MIEZE measurement is the intermediate scattering function I (Q; ), which is determined by dividing the sample contrast by an appropriate resolution contrast: IðQ; Þ ¼ C sample =C resolution . Since all MIEZE measurements are normalized to the instrumental resolution function (which depends equally on the damping factor), the damping factor cancels out and therefore does not need to be taken into account explicitly. Nevertheless, it is important to track the damping factor, to not increase the error bars of the contrast beyond a reasonable limit.
Keeping in mind that at least three parameters I mean (the time average), I 0 (the amplitude) and 0 (the arbitrary phase) must be extracted from the signal, a minimum of three time bins is necessary for an unambiguous reconstruction. In contrast to classical neutron spin-echo, it is not possible to use a 3 He counter as a MIEZE detector. In fact, the detector requirements are quite demanding: a MIEZE detector requires high spatial and temporal resolution, while in addition the thickness of the conversion volume of the detection system in the neutron flight direction must not exceed the size of the MIEZE group which decreases with increasing MIEZE time (Schmidt et al., 2010). Currently a CASCADE detector with 16 time bins is used to detect the MIEZE signals at RESEDA (Kö hli et al., 2016). The detector consists of eight  10 B-coated detection foils, with a conversion layer thickness of 0.8-1.5 mm and a pixel size of 1.56 mm. The current CIPix ASIC preamplifier readout of the detector electronics is able to handle frequencies up to 10 MHz. The recent improvements of the RESEDA instrument (Jochum et al., 2020) have pushed the first-generation CASCADE detectors to their limits. Nevertheless the time-stamp accuracy of the detected events still has a reserve, since the internal clock and the FPGAs run at a frequency of 40 MHz, leading to a maximum binning inaccuracy of 12.5 ns. These constraints imposed by the detection system limit the maximum MIEZE frequency at RESEDA to f MIEZE = 10 MHz/16 = 625 kHz. In this regime, the damping induced by the sinc function is only 0.64%. The contrast is extracted from I mean , I 0 and 0 , which are determined through a sine fit across the 16 time channels. This fitting procedure is calculation intensive and cannot be performed in real time alongside the data acquisition.
To increase the time resolution ( f MIEZE ), a practical solution is to apply the same routine with a reduced number of time bins. Alternatively, one may find an unbiased estimate to reconstruct the parameters from the minimum necessary time bins by taking the time integration of the detector into account. In the following sections, a reconstruction procedure of the underlying parameters will be deduced using only four time bins. This relaxes the required data collection interval by a factor of four corresponding to f max MIEZE = 2.5 MHz, thereby increasing the time resolution by a factor of four. Although three time bins are the optimal choice to cover the highest frequencies, we focus here on four time bins because of their backwards compatibility with older data sets histogrammed in 16 time bins.

Reconstruction of the MIEZE signal
As a starting point for the reconstruction of the MIEZE signal, we give the mathematical description of the timedependent event rate IðtÞ as recorded by the detector. This signal may be split into a time-dependent and a time-independent contribution (I mean ). The latter describes the intrinsic background and all of the contrast reductions such as incoherent scattering, spin leakage and sample dynamics. The sinusoidal time dependence is characterized by the amplitude I 0 , the duration T ¼ 1=f MIEZE and the phase shift 0 . These combine to give IðtÞ as Since the time binning of events in the detector is equal to an integration over time of IðtÞ in the respective interval, one may write the number of detected events in the kth interval I k as IðtÞ dt; ð4Þ with k ¼ 1; 2; 3; . . . ; N for N bins. Normalizing I k by I mean corresponds to the probability of a single event occurring in the kth interval. For a subdivision into four intervals (N ¼ 4, cf. Fig. 2 gray shaded area for I 1 ) one may rewrite (4) as follows: Summing neighboring intervals and simplifying the results yields Adding the next-nearest-neighbor intervals (I k + I kþ2 ) yields only the direct component (first terms) while the phase information is lost: This is a direct consequence of the signal's harmonic periodicity.

Figure 2
A typical time-dependent sinusoidal intensity variation with phase 0 ¼ =8 and a contrast C ¼ I 0 =I mean that is defined by the amplitude I 0 and the mean value I mean . The gray shaded area I 1 normalized to I mean is the probability of a single event being detected in the first interval from the division of each oscillation of IðtÞ into four equally long time bins.
Since equations (6a)-(6d) are sums of neighboring intervals, one may use two independent but identical detector readouts to measure two separate time intervals that are =2 phase shifted relative to each other. This yields equivalent information, but allows for a doubling of f max MIEZE . Precise phase determination of harmonic signals is well established using quadrature detection in optical interferometry (Rerucha et al., 2012) or signal processing where =2 phase-shifted signals [(6a)-(6d)] are combined to reconstruct the unknown phase 0 : It is also possible to deduce the phase by subtracting equations (5a)-(5d) from each other: Equation (9) shows that in principle one interval can be neglected. However, for this approach information in the form of counts is ignored within that interval, thus reducing the overall statistics and accuracy. The average over equation (9) equals equation (8).
Using C ¼ I 0 =I mean the reconstructed (rec) contrast may be deduced as well, by combining either equations (6a) and (6c) or (6b) and (6d): Of course, the accuracy of the evaluated contrast is strongly coupled to the accuracy of the estimated phase and diverges at the singularities, i.e. when cos 0 or sin 0 tend towards zero. In order to avoid the singularities we apply equations (10a) and (10b) for the appropriate case: C rec ¼ C 1;rec ; for cos 0 ! sin 0 ; C 2;rec ; for cos 0 < sin 0:

& ð11Þ
We emphasize again that this simple reconstruction of contrast and phase, using only four time bins, allows an increase in Fourier time by a factor of four. Additionally, this reconstruction method (unlike the previously used method) does not require any computationally intensive fitting. This will speed up data reduction immensely and allow for real-time data reduction, which will make it possible to optimize measuring times and use allocated beamtime more efficiently.

Estimation of the confidence interval
Next, we discuss how many events are necessary to determine phase and contrast with a desired accuracy. For this we will compare three different attempts: (i) the 16-time-bin fitting method used so far at RESEDA (fit,16), (ii) the four-time-bin fitting method (fit,4), (iii) the four-time-bin reconstruction method (rec). The procedure does not take into account a possible phase jitter of the detector signal. The MATLAB code utilized for these calculations has been made available for reference (Jochum et al., 2021). As a first attempt, the uncertainties are estimated using Gaussian error propagation with the relative errors 1=ðI k Þ 1=2 . Deducing the partial derivatives is straightforward. Less obvious is the estimate of the total errors ÁI 1 ; . . . ; I 4 , due to their mutual dependence.  Moreover, the total errors also depend on C 0 and 0 . While a generalized Gaussian error propagation would account for the covariance between all parameters, it is not able to give a reliable answer in the limit I ! 0. Therefore, we applied simulations and executed them for various initial phases ( 0 = 0 , 15 , . . . , 120 ) and contrasts (C 0 ¼ 0:05; 0:1; . . . ; 0:95). The first ten single events with the desired sinusoidal distribution are generated using the pseudo-random generator of MATLAB and histogrammed subsequently. For a given C 0 and 0 , the probability of falling in a certain time bin is determined by equations (5a)-(5d). Subsequently, the phase and contrast are calculated according to the three different methods. Next, new events are added to this run and the evaluation is repeated recursively.
The number of added events in such a series increases logarithmically. This ensures a low computational burden over a large dynamic range of events (here over five orders of magnitude) and keeps the evaluation equally weighted in a logarithmic representation. Finally, the results are compared with each other. Figs. 3(a) and 3(d) show the phase () and contrast (C) for one run. The phases estimated for the fourpoint fitting method (green) and the reconstruction (red) are identical within error for more than 30 events.
For a low number (<30) of events, the phase and contrast values have larger deviations from the true values as a result of insufficient statistics. As expected from equation (2b), both fitting methods show biased (damped) contrast estimates. For the contrast C 0 ¼ 0:85 presented in Figs. 3(d) and 3(e), the expected damping according to equation (2b) is 0.64% Â C 0 = 0.0054 for the 16 time bins (blue) and 10% Â C 0 = 0.085 for the four time bins (green). This shows that the estimates inferred from the reconstruction method are unbiased.
To estimate the standard deviations of the phase and the contrast, the simulation was run 500 times.
From these data, the average phase ( avg ) and contrast (C avg ) as well as their corresponding standard deviations ( and C ) were calculated as a function of events I [cf. Figs. 3(b), 3(c), 3(e) and 3( f)]. While the average phase is estimated correctly, the unbiased estimate for the contrast bears the expected damping. In agreement with the experimental behavior, the estimated standard deviations and C (for the reconstruction and fitting procedures) decrease with the same asymptotic behavior as the total number of events Figs. 3(c) and 3(f )]. This proves that the applied estimator is consistent.
The relationship between standard deviation and events, for both the phase and contrast, is described by simple power laws: From a linear fit to the log-log plot of the estimated standard deviations [cf. Figs. 3(c) and 3( f)] for more than 30 events, the power-law exponents ( and C ) may be inferred: To test the generality of this power-law behavior and to determine the missing parameters ( and C ) for varying contrasts, the simulations for C 0 ¼ 0:05; 0:1; . . . ; 0:95 were repeated while keeping the initial phase fixed at 0 = 60 . In  and 4(d)] remain nearly unchanged, confirming that the use of a normal distribution to approximate a Poisson distribution is well justified. However, the coefficients and C show quantitatively distinct dependencies on the initial parameter C 0 [cf. Figs. 4(b) and 4(e)]. was observed to follow an exponential decay with increasing C 0 : with decay constants 2; and 1; which vary slightly depending on the method. The functional dependence of the contrast is less obvious, and the parabolic fits (solid lines) in Fig. 4(e) are not ideal. Since Ã and Ã (Ã ¼ or C) depend on each other, as the fits are over-parameterized, Ã was recalculated with the constraint Ã ¼ À0:5. For the sake of clarity, Ã is renamed Ã in the following if the constraint ð Ã ¼ À0:5Þ is applied. The resulting fits are plotted in Figs. 4(c) and 4( f). Compared with Fig. 4(b), the exponential dependence of is maintained. Furthermore, C can now be described well by a shifted halfnormal distribution: To show that these findings hold for the relevant range of phases, this procedure was repeated for 0 = (0 . . . 120 ) in steps of 15 . To confirm the 90 periodicity of the angular dependence, the interval was extended to 120 . This yields a set of curves comparable with the ones in Fig. 4, which are color-plotted in Fig. 5, highlighting their behavior throughout the entire parameter space. The plots confirm that the fitting parameters deduced with these techniques are practically independent of 0 . Note that, due to the periodicity of the harmonic functions, these findings are valid for all phases. Combining equations (12a), (12b), (14) and (15), we find the analytical equations for the estimate of the standard deviation: The parameters 1; , 2; , 1;C and 2;C we found for the different methods presented here are summarized in Table 1.
Equations (16a) and (16b) and Table 1 show that the deduced errors depend strongly on the initial contrast C 0 and the applied methods. The most obvious variation is observed for the parameters 1;C and 2;C . 1;C determines how quickly C drops with increasing initial contrast, whereas 2;C scales the absolute magnitude of C . We emphasize that the error bars deduced for the contrast using the four-point fitting method must be treated carefully, since the procedure of inferring the estimate is biased. Re-scaling this contrast and its error with the damping factor of 0.9 given by (2a), the same error observed for the reconstruction method is maintained. However, as long as the same procedure is used for sample and resolution measurements, these effects cancel out and can therefore be neglected.

Conclusions
We have presented an algorithm to deduce the contrast and phase of a sinusoidally modulated time series sampled at four data points per oscillation. Both contrast and phase are Color-map plots of parameters (a), (b), C (c) and C C (d) for varying initial phases 0 and contrasts C 0 using the reconstruction. All four parameters are nearly independent of the initial phase 0 . While the parameters and C remain nearly constant even through varying initial contrasts, decays exponentially and C C exhibits the same trend as shown in Fig. 3(c) for increasing C 0 . On the other hand, their accuracy is independent of the initial phase. The reconstruction trades a higher time resolution for less accurate contrast. Quantitatively, this factor is better than C0;fit;16 = C0;rec ! 0:9 compared with the fitting method, but may be compensated by increased statistics, i.e. around 20% prolonged counting time. However, using the reconstruction, there is no fitting procedure involved, which significantly reduces the required computational burden. Thus, this method may be readily applied to a large number of detector pixels as the measurements proceed in time. As mentioned above, real-time data evaluation will lead to a more efficient use of measurement time, decreasing the time needed for each experiment. Most importantly, this new method solves one of the main limitations afflicting the MIEZE resolution. Using a CASCADE-type detector (Kö hli et al., 2016) with a maximum time resolution of 100 ns (10 MHz), the maximum intensity modulation frequency for 16 time channels is 625 kHz, which, at 6 Å , with the current dimensions at RESEDA yields a MIEZE (Fourier) time of $3 ns. In contrast, the resolution limit using the four-point method is $12 ns at 6 Å or $100 ns at 12 Å .
Having extended the time resolution limit using the fourpoint reconstruction method, the next challenge for MIEZE data acquisition is the pixel size of the detector. The reason for this is that the coherence volume of the MIEZE signal is indirectly proportional to the wavelength of the incoming neutron beam, the width of the wavelength band and most importantly f MIEZE . Thus, for intensity modulation frequencies at or above 2.5 MHz, extremely flat detector surfaces are needed to minimize phase differences within a single pixel. A 10 B layer on a solid surface instead of Kapton foil could be a possible solution. Furthermore, a spherical detector foil shape would suppress the phase rings which occur on flat surfaces due to variations of path lengths (Schober et al., 2019).