Machine learning approach to muon spectroscopy analysis

Machine learning (ML) methods are now widely used in many areas of physics, usually as a tool to analyse large amounts of data [1–3]. These techniques are particularly useful in regression, classification and dimensionality reduction tasks which are often required in processing scientific data. Specifically in condensed matter physics, ML is well suited for many tasks ranging from predicting materials properties based on existing databases and pattern recognition in specific experimental data to analysing theoretical models of quantum materials. Prominent examples include the prediction of novel materials [4–6], identification of phase transitions in models of magnetic materials starting from Ising models [7–12], reaching complex spin liquids in Heisenberg systems [13] and the detection of entanglement transitions from simulated neutron scattering data [14]. ML algorithms were also proven to be state of the art techniques in simulations of wave functions [15] or density matrices [16–19] for many-body quantum systems and the tomographic reconstruction of many-body wave functions from experimental data [20].

Much of the research in this area so far is concerned with simulation or analysing simulated data, however it has also been shown that such techniques can detect phase transitions from piezoelectric relaxation measurements [21] or discovering existence of translational symmetry-breaking states from real, electronic quantum matter images [22]. Here we want to apply a simple dimensionality reduction algorithm to real data from muon spin rotation (μSR) experiments [23] to see if we can detect phase transitions for a range of different materials. We decided to use the data from this type of experiment since models used in μSR data analysis require previous understanding of the local environment of muons inside probed sample, which is not always easily available. Therefore, as an alternative, we propose the use of linear principal component analysis (PCA), a simple unsupervised ML technique which does not make any prior assumption, yet is known to reveal correlations within the data. By demonstrating that this approach works, we propose that it may serve as a more unbiased way of detecting phase transitions observed in μSR experiments. In this paper we apply PCA to μSR data from a small number of superconducting and magnetic materials whose physics are known to differ widely from each other. In particular we explore the technique for data from time reversal symmetry breaking (TRSB) superconductors, which are among the most difficult to analyse, since changes in experimental data are very subtle. Other materials that we have tested are a symmetry breaking antiferromagnet (BaFe₂Se₂O) and a spin liquid (LuCuGaO₄). We find some evidence that PCA can detect important features such as phase transitions. We also find that when the system is trained on all the materials, taken together, the results improve—even though the materials chosen have different underlying physics.

The paper is organised as follows. In section 2, we briefly present the set up of the muon spectroscopy experiment and the current method of analysing the data from it. In section 3, we present the principal component (PC) analysis in general and how we used it in practice. Then, in section 4, we move on to results of applying PCA to data from different materials and discuss in detail how the method performs. We summarise the results in section 5.

The general setup of a μSR experiment design to measure the local magnetic environment consists of spin-polarised muons being implanted into a sample, which is surrounded by multiple positron detectors. Once they enter the sample, muons will interact with the atoms causing them (muons) to thermalise and eventually implant themselves at some sites of the system. The spin of the muons will start to precess due to the local magnetic field and the muons will eventually decay into positrons and neutrinos with a mean life time of 2.2 μs. The positron velocity direction is directly connected to the muon spin orientation at the time of decay [24–27] and therefore the intrinsic magnetic field of the sample will affect the final distribution of positron detection events.

A commonly used setup is to have symmetrical detectors in front of (F) and behind (B) the sample (with respect to the muon beam). The quantity that we are interested in is the difference in number of counting events between the two detectors as a function of time N_i (t), i ∈ {F, B}, called the asymmetry function

$$\begin{equation}A\left(t\right)=\frac{{N}_{\text{B}}\left(t\right)-{N}_{\text{F}}\left(t\right)}{{N}_{\text{B}}\left(t\right)+{N}_{\text{F}}\left(t\right)}.\end{equation} \tag{ 1 }$$

The analysis of the data involves fitting specific asymmetry curves to the experimentally-obtained curve. Given some knowledge of the underlying physics for a particular material and/or some justified assumptions, a model can be formulated, and the asymmetry curve can be derived from it. In some simple cases appropriate closed-form expressions can be derived [28, 29], though more generally ad hoc calculations are necessary [30]. For some systems, our understanding is still not sufficiently developed for such predictions—for instance, the theory of zero-field muon spin relaxation (ZF-μSR) in superconductors with broken time-reversal symmetry (TRS) is still in its infancy [31].

In practice, for complex systems it is customary to use a phenomenological expression featuring several adjustable parameters. Electronic order can then manifest as a temperature-dependence of those parameters. For instance, in ZF-μSR investigations of superconductors [32] one often fits:

$$\begin{equation}{A}_{\text{phen.}}\left(t\right)={A}_{0}{G}_{\text{KT}}\left(\sigma ,t\right)\mathrm{exp}\left(-\lambda t\right)+{A}_{\text{bckg}},\end{equation} \tag{ 2 }$$

where G_KT(σ, t) is the Kubo–Toyabe function describing coupling to static, randomly-oriented magnetic moments [28, 29, 33, 34] with relaxation rate λ and Gaussian magnetic field strength distribution with standard deviation σ. The parameters σ, λ, A₀, A_bckg are then interpreted to describe distinct relaxation mechanisms. In conventional superconductors these parameters tend to evolve smoothly through the superconducting critical temperature, T_c. In other systems, marked changes in some of these parameters occur at T_c [35]. These are often interpreted as evidence of broken TRS and in some systems this has been confirmed by Kerr effect or SQUID magnetometry. Quite frequently, it is found that only one of the fitting parameters in equation (2) depends on temperature. This is usually either σ or λ, which naturally leads to a classification of TRS-breaking superconductors. We note, however, that the relaxation rates involved are very small, meaning that only a small portion of the curve described by equation (2) is represented in the experimental data sets (due to the finite lifetime of the muon). As a result, this classification may not always be as robust as would be desirable. For instance, some superconductors that are expected to have very similar underlying physics can fall in different classes. Such is the case of the proposed nonunitary triplet superconductors LaNiC₂ [36] and LaNiGa₂ [37], whose asymmetry functions are best described by a temperature-dependent σ and λ, respectively, in spite of experimental [38, 39] and theoretical [32] evidence of very similar underlying physical mechanisms. Likewise, the muon spin relaxation rate in spin glasses can often be described by a stretched exponential function (with temperature-dependent exponent), reflecting the variation in local spin fluctuation rates as well as non-exponential decay at muon sites [40–42]. However, fitting experimental data can give parameter values that are not expected from standard models/numerical analysis [43]. In conclusion, it would be highly desirable to have a way of analysing the temperature-dependence of μSR spectra that can detect electronic ordering transitions without the need to assume any a priori fitting functions.

To analyse the data from a muon spectroscopy experiment without making assumptions about the physical nature of the materials, we decided to use an unsupervised ML technique called PCA [7, 44, 45]. The concept behind it—in the context of muon spectroscopy experiment and asymmetry functions—is presented in figure 1. We can think about different experimental measurements as points in some data space with N dimensions. In the case of muon spectroscopy, each dimension i = 1, 2, ..., N represents a time window t_i , within which the positron detections are measured. If the measurements are not random but correspond, for example, to the same material at different temperatures, we expect correlations between those points. PCA can detect these correlations by first removing the average of all experimental curves, then measuring the covariance for each dimension and linearly transforming the coordinates so that the new basis of the data space consists of only few directions that capture most of the covariance. The vectors of this new basis are called PCs and can be thought of as the most common deviations from the average curve. We can reconstruct all of the measurements used in the analysis by adding to the average a linear combination of PCs. We can also represent each curve by specifying its projections onto the PCs, which are often called PC 'scores'. Thus, PCA provides us with a more compact description of the experimental data and additionally we can recover information about linear correlations from their magnitudes (or PC scores) and shapes (the PCs, or PC vectors).

Zoom In Zoom Out Reset image size

In the example shown in figure 1(c), most of the data lies in two-dimensional space π₁ × π₂. PCA finds new orthogonal directions (π'₁, π'₂), because there exist linear correlation between the π₁ and π₂ coordinates of data points. We can now specify each asymmetry curve by its projection onto π'₁, whereas before we would have to state both π₁ and π₂ coordinates. We do lose some information about the individual data points in this way, but we gain in the more compact representation of asymmetry curves. Usually, more than one PC is needed to represent the data well. The number of important PCs varies with different data sets and can be decided by looking at how much covariance each PC holds.

We now present a more specific description of the PCA method. Each measurement can be represented as a vector ${\mathbf{a}}_{j}={\left({A}_{j}\left({t}_{1}\right),{A}_{j}\left({t}_{2}\right),\dots ,{A}_{j}\left({t}_{N}\right)\right)}^{\mathrm{T}}$, ⁶ with its values equal to the values of asymmetry function at specific times and the index j = 1, 2, ..., M taken to label the distinct measured asymmetry curves that we want to analyse by the algorithm. We further assume that all measurements were recorded for the same set of N measurement times t_i , taken relative to the time for implanting the muon into the material. We combine the vectors a_j in column form to construct a matrix A

$$\begin{equation}\mathbf{A}=\left[\begin{matrix}\hfill {A}_{1}\left({t}_{1}\right)\hfill & \hfill {A}_{2}\left({t}_{1}\right)\hfill & \hfill \dots \hfill & \hfill {A}_{M}\left({t}_{1}\right)\hfill \\ \hfill {A}_{1}\left({t}_{2}\right)\hfill & \hfill {A}_{2}\left({t}_{2}\right)\hfill & \hfill \dots \hfill & \hfill {A}_{M}\left({t}_{2}\right)\hfill \\ \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill {A}_{1}\left({t}_{N}\right)\hfill & \hfill {A}_{2}\left({t}_{N}\right)\hfill & \hfill \dots \hfill & \hfill {A}_{M}\left({t}_{N}\right).\hfill \end{matrix}\right]\end{equation} \tag{ 3 }$$

In the next step we remove the mean of each vector dimension (i.e., averaging over the column index) so that the whole data is centred around the coordinate origin, as shown in figure 1. We end up with a matrix X with elements given by

$$\begin{equation}{\left[\mathbf{X}\right]}_{ij}={A}_{j}\left({t}_{i}\right)-\frac{1}{M}\sum\limits _{k=1}^{M}{A}_{k}\left({t}_{i}\right).\end{equation} \tag{ 4 }$$

The most common way for obtaining PCs is to perform a singular value decomposition of X. To this end, we evaluate the covariance matrix

$$\begin{equation}\mathbf{S}=\frac{1}{M-1}\mathbf{X}{\mathbf{X}}^{\mathrm{T}},\end{equation} \tag{ 5 }$$

such that the eigenvectors of S are the PCs and the corresponding eigenvalues indicate the amount of covariance captured by the given PC. If we write the eigenvectors into a matrix U, then a table of scores C for each measurement can be obtained by the matrix product

$$\begin{equation}\mathbf{C}={\mathbf{U}}^{\mathrm{T}}\mathbf{X},\end{equation} \tag{ 6 }$$

and the full reconstruction of the initial experimental data is expressed as

$$\begin{equation}\mathbf{R}=\mathbf{U}{\mathbf{C}}^{\mathrm{T}}.\end{equation} \tag{ 7 }$$

The previously discussed usefulness of the method derives from the fact that we can choose only the few PCs that capture most of the covariance in order to accurately reconstruct the initial data. Naturally, a large reduction in the number of relevant PCs does not have to arise for all possible data sets, as singular value decomposition only performs a linear transformation—in particular, if the data has non-linear correlations the method will not perform well. Fortunately, looking at the eigenvalues of S, one can decide if the linear PCA is sufficient, based on the decrease of PC scores which is often illustrated in a so-called scree plot of the PC scores against their index.

In order to account for the experimental noise in the data, we have re-binned raw data into new time windows according to the measurement error. Since the error increases with time, wider time windows are required at larger times to get comparable errors. Hence, available measurement points are more widely spaced at later times, as can be seen in figures 2(a)–(c). It is important to re-bin all of the measurements simultaneously because all time-windows t₁, t₂,...,t_N in our matrix A have to be the same for all columns for the PCA to be well defined. Note that this specification mirrors the treatment in regression methods, where less weight is attributed to data at long times to account for the larger measurement errors. Our binning procedure is discussed in detail in appendix A.

Zoom In Zoom Out Reset image size

3.1. Philosophy of our PCA approach

4.1. PCA for simulated data

To illustrate characteristic results of performing PCA on asymmetry functions, we first consider an example application to synthetic data generated from model Kubo–Toyabe functions G_KT(σ, t) with added error E(t). Each such simulated asymmetry function was taken from the general form given by

$$\begin{equation}{A}_{\text{sim}}\left(t;T\right)={A}_{0}{G}_{\text{KT}}\left(\sigma \left(T\right),t\right)\mathrm{exp}\left(-{\Lambda}t\right)+{A}_{\text{bckg}}+E\left(t\right),\end{equation} \tag{ 8 }$$

where we have further encoded a dependency of σ(T) associated with a symmetry-breaking phase transition in many superconductors with TRSB, such that σ(T) ∼ constant for T > T_c, while varying linearly below T_c. The error values were generated from a Gaussian distribution N(μ = 0, Σ_sim) centred on zero and with a standard deviation ⁷ Σ_sim(t) depending on time after muon implantation as:

$$\begin{equation}{{\Sigma}}_{\text{sim}}\left(t\right)=R\left({a}^{t}+b\right).\end{equation} \tag{ 9 }$$

Errors observed in real measurements increase with time t due to the overall smaller number of events detected at later times. The parameters A₀, σ(T), Λ, A_bckg, R, b, and a were chosen to match experimental data of one of the TRSB superconductors we studied (LaNiGa₂). In addition to the parameters reflecting experimental conditions, we studied the effect of different error amplitudes R (which in experiments would correspond to experiments undertaken with different amounts of time allocated for integrating the signal) in order to verify robustness of the PCA approach. Our results from the application of PCA to these simulated data is displayed in figure 3. We have included four possible cases of 'noise' amplitudes ranging between no error and twice the error we expect from our measurements. The PCA on clean data clearly captures the transition temperature T_c assumed in the simulated data, which separates regions of temperature with or without variation of the PC scores with T. The first principal score dependency is found to be very robust to added noise, even for the cases where the error is much larger than expected experimentally. By contrast, the second PC does not seem to hold any useful information for realistic noise level. Note also the small overall scale of the second PC score. Nevertheless, the phase transition is always clearly visible in the 1st PC, which motivates using PCA for experimental data.

Zoom In Zoom Out Reset image size

4.2. PCA for experimental data

We applied PC analysis to data from zero-field μSR experiments for a range of different materials. Among them are TRS breaking superconductors ⁸ (LaNiGa₂, LaNiC₂, LaNi_1−x Cu_x C₂), spin liquid (LuCuGaO₄) and an antiferromagnet (iron oxyselenide BaFe₂Se₂O). We first performed the analysis for each material separately. The shape of the two most important PCs and the dependence of the scores on temperature are presented on figure 4.

Zoom In Zoom Out Reset image size

Our technique worked best for the antiferromagnetic material (first column on figure 4), for which both expected phase transitions are clearly visible. Although the magnetic behaviour of the antiferromagnet BaFe₂Se₂O is relatively simple [54], this understanding has been challenging to arrive at: (a) T_N ∼ 240 K is clear from neutron powder diffraction experiments but is more subtle in magnetic susceptibility measurements [54–56] due to the layered nature of the material; (b) magnetic susceptibility data collected on several samples suggest a magnetic phase transition at ∼115 K [54, 55] which is now thought to be due to Fe₃O₄-related impurities and is not intrinsic to the main phase [54]; (c) there's no evidence for the low-temperature ∼40 K phase transition from neutron powder diffraction [54] or heat capacity data [56] and this phase transition is thought to involve freezing of spin fluctuations. It is striking that this unsupervised ML analysis correctly identified the two phase transitions intrinsic to BaFe₂Se₂O without the need for complementary data. We think that this reflects the strength of both the PCA analysis and the μSR technique.

The changes in the asymmetry function are more subtle for the superconducting materials (second-fifth column on figure 4), but the behaviours of PC scores still change at the critical points. We stress that in conventional superconductors we would not expect any change of zero-field muon-spin relaxation at the superconducting transition temperature T_c. By contrast, LaNiC₂ and LaNiGa₂ are known to exhibit such changes, and this is believed to be a manifestation of their internally-antisymmetric, non-unitary triplet (INT) pairing states with TRSB [32]. In the case of LaNiC₂ (third column), we only have one point above the phase transition and therefore we do not expect visible change. That is confirmed in PC score plots. Worth mentioning is also the LaNi_0.9Cu_0.1C₂ case, in which there seem to be more than one critical point, at least in the behaviour of the 1st PC score. That might be caused by some other phase transition but more probably it is caused by limitations of the method. One solution to that problem would be to look also at the 2nd PC score, where only one transition point is prominent. Overall, linear PCA seems to be performing better for the spin liquid and antiferromagnetic materials than for the TRS breaking superconductors analysed in this paper, as is evidenced in our scree plots (the last row on figure 4). For the first four materials, even the last few PCs hold a significant amount of covariance ⁹ . That may imply that the data has non-linear correlations or that we did not have not enough data available for these types of materials, since most ML algorithms perform better the more data is provided. It is important to note that we can still resolve the changes in the scores of 1st and 2nd PCs, at least for LaNiGa₂ and LaNi_0.9Cu_0.1C₂.

The last studied material is a proposed spin liquid—LuCuGaO₄. Muons have been used as a proof of a spin liquid state, as it can be argued that the resultant dynamics could show a plateau in the relaxation rate with reducing temperature where no long range order is detected [57]. In our case the PCA shows no evidence for a phase transition, even though a plateau is observed, likely indicating there is no phase transition as the proposed liquid state is entered.

Because the most significant PCs for the TRS breaking superconductors look similar for all the cases studied, in hope of improving results for TRSB systems, we proceed to apply PCA to all of the experimental data simultaneously. The results for PC scores are presented on figure 5 and the improvement of PCs are shown on figure 6. The PCs are now much smoother functions and additionally, only three PCs are sufficient to capture 80% of the observed covariance. The scores of the first PC did not change much for all materials, despite their different physical properties. This is probably connected to the fact that all data come from the same type of experiment and all asymmetry functions are similar in general.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that the data for both LaNi_1−x Cu_x C₂ materials was previously unpublished. Using PCA, we were able to see signs of the superconducting phase transition in zero-field muSR experiments, which is the first evidence that the TRSB of LaNiC2₂ also exists in these Cu-doped materials. The increased onset temperature is consistent with the known enhancement of T_c with Cu doping [53].

We have proposed the use of PC analysis to process muon spin spectroscopy data, and in particular to aid with the identification of features relating to phase transitions in the probed materials. Our results demonstrate that the representation of the observed asymmetry functions in the space of PC vectors is sensitive to changes in the physics of the observed system. In particular, the evolution of PC scores as a function of tuning parameters provides insights into the location of possible phase transitions. Comparing this analysis to a more conventional approach, based on regression analysis using standard fitting functions, we find that PCA is typically at least as sensitive, if not more. More importantly, the PCA approach is free from any underlying assumptions about the physics of the observed material: rather than assuming a specific form of a fitting function (e.g. Kubo–Toyabe or stretched exponential), PCA discovers the PCs that describe a given system, without human intervention. This is the salient feature of the method we put forward and it means that the same, universal analysis can be applied to any material.

In addition to the universality of our method, we have found that the quality of the results is enhanced when data for multiple materials are analysed as a joint dataset, even when the underlying physics of each system being considered are quite different. The ability to thus enhance understanding gained from a new experiment based on existing data goes beyond the possibilities of preexisting approaches, where data for each material is necessarily analysed and fitted in isolation, and overarching commonalities are anticipated in advance by the formulation of a suitable fitting function. We anticipate this could offer great advantage when deployed in large-throughput user facilities. In particular, given the advantages gained from combining multiple data sets, our results suggest a new way to leverage recently-developed open-data tools and policies [58, 59].

We hope that our unsupervised ML approach to muon spectroscopy data analysis could become one of the standard tools used in that field. In addition to its virtue, noted above, of providing a unified way of treating all muons data, we believe our approach can also accelerate future experiments, as the treatment within this framework will require less data to be collected before signatures of the physics can emerge—especially if data from previous experiments is used to enhance the analysis of new materials as outlined above. In addition, the simplicity of the analysis means that it could easily be performed immediately while experimental measurements are being taken, thus opening the possibility to inform the conduct of the experiment in real time. At the other end of the spectrum, it is also possible to conduct experiments where much larger data sets are gathered [60]. Our simulations suggest that our method applied to such data might yield valuable new insights into phase transitions. They also would be ideal additions to such past-experiment data bank. Given the advantages gained from combining multiple data sets, our results should encourage the community to gather historic and future measurements in a common database in order to harvest the benefits of this approach.

A possible future extension of our work would be to deploy an additional unsupervised ML technique to analyse the output of our analyses as presented here. The principal score dependencies (as shown in figure 5) in the method presented here still need to be processed by human eye to establish phase transition temperatures. There exist ML tools that could categorise different phases from the PC scores, that have been shown to work well with model data [61]. It would be interesting to apply them to our problem.

T Tula implemented the PCA algorithm and performed the analysis of the simulated and experimental data presented in this paper under supervision from J T Quintanilla and G Möller. S R Giblin, A D Hillier, E E McCabe and S Ramos provided, formatted and commented on the experimental data. D S Barker performed a preliminary study of PCA applied to experimental and simulated μSR data under supervision of J Quintanilla, with further input from S Gibson. T Tula wrote the manuscript in close consultation with G Möller and J Quintanilla and with input from all co-authors.

We would like to thank Stephen Blundell, Tom Lancaster and Roberto De Renzi for helpful discussions about the content of the paper.

TT is supported by the EPSRC via a DTA studentship under Grant No. EP/R513246/1 and by the School of Physical Sciences, University of Kent. SR and EEM are grateful to Mr Ben Coles (for BaFe₂Se₂O synthesis) and to Dr Fiona Coomer (experimental support) for the μSR data from reference [54]. JQ acknowledges support from the EPSRC under the project 'Unconventional superconductors: new paradigms for new materials' (Grant No. EP/P00749X/1). GM gratefully acknowledges support by the Royal Society under University Research Fellowship URF∖R∖180004.

The data that support the findings of this study are available upon reasonable request from the authors.

The error in muon spectroscopy measurements increases at later times due to the smaller number of overall positron detection events. Since PCA treats each dimension (time window) equally, one needs to pre-process the raw data by re-binning the time windows. While this differes from the exact way in which errors are treated in a standard fitting procedure, where a weight function is applied to give less weight to data with larger errors, our approach is broadly equivalent in that it makes sure that the standard errors of the rebinned time points are roughly the same for every measurement. Specifically, we have set up an algorithm for re-binning the data, such that each new time bin holds the same magnitude of error averaged over all measurements. To illustrate how the algorithm proceeds, let us consider the matrix E, holding the raw values of standard deviations at each time point and for each asymmetry curve (similarly to the matrix A from equation (3)):

$$\begin{equation}\mathbf{E}=\left[\begin{matrix}\hfill {E}_{1}\left({t}_{1}^{\text{raw}}\right)\hfill & \hfill {E}_{2}\left({t}_{1}^{\text{raw}}\right)\hfill & \hfill \dots \hfill & \hfill {E}_{M}\left({t}_{1}^{\text{raw}}\right)\hfill \\ \hfill {E}_{1}\left({t}_{2}^{\text{raw}}\right)\hfill & \hfill {E}_{2}\left({t}_{2}^{\text{raw}}\right)\hfill & \hfill \dots \hfill & \hfill {E}_{M}\left({t}_{2}^{\text{raw}}\right)\hfill \\ \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill {E}_{1}\left({t}_{N}^{\text{raw}}\right)\hfill & \hfill {E}_{2}\left({t}_{N}^{\text{raw}}\right)\hfill & \hfill \dots \hfill & \hfill {E}_{M}\left({t}_{N}^{\text{raw}}\right)\hfill \end{matrix}\right],\end{equation} \tag{ A.1 }$$

where M is the number of asymmetry functions that we consider in the analysis and N corresponds to the number of time windows in raw data. We set the first time window t₁ to be equal to ${t}_{1}^{\text{raw}}$, which holds the average error of

$$\begin{equation}{\bar{E}}_{{t}_{1}}=\frac{1}{M}\sum\limits _{i=1}^{M}{E}_{i}\left({t}_{1}^{\text{raw}}\right).\end{equation} \tag{ A.2 }$$

Then we iterate over ${t}_{j}^{\text{raw}}$ to create new bins in the following way: suppose we created a new bin t_k−1 by including raw data up to the original bin at time ${t}_{j-1}^{\text{raw}}$. We then evaluate

$$\begin{equation}\frac{{\bar{E}}_{{t}_{1}}}{{\bar{E}}_{{t}_{j}^{\text{raw}}}},\end{equation} \tag{ A.3 }$$

with ${\bar{E}}_{{t}_{j}^{\text{raw}}}=\frac{1}{M}{\sum }_{i=1}^{M}{E}_{i}\left({t}_{j}^{\text{raw}}\right)$. We know that (A.3) is smaller than 1, since the standard deviation is increasing with time due to decreasing muon counts. If (A.3) is close to one, then the amount of averaged error is similar to the first bin and we can leave the time window ${t}_{k}={t}_{j}^{\text{raw}}$. However, if it reaches certain threshold, we add another time window ${t}_{j+1}^{\text{raw}}$ and evaluate:

$$\begin{equation}\frac{{\bar{E}}_{{t}_{1}}}{\sqrt{{\left({\bar{E}}_{{t}_{j}^{\text{raw}}}\right)}^{2}+{\left({\bar{E}}_{{t}_{j+1}^{\text{raw}}}\right)}^{2}}}.\end{equation} \tag{ A.4 }$$

We repeat this procedure until

$$\begin{equation}\frac{{\bar{E}}_{{t}_{1}}}{\sqrt{\sum\limits _{l=0}^{{L}_{k}-1}{\left({\bar{E}}_{{t}_{j+l}^{\text{raw}}}\right)}^{2}}}\approx 1,\end{equation} \tag{ A.5 }$$

and we set a new time bin ${t}_{k}=\frac{1}{{L}_{k}}{\sum }_{l=0}^{{L}_{k}-1}{t}_{j+l}^{\text{raw}}$.

When the different asymmetry functions come from the same material, we expect our method of re-binning to work well. One might expect that problems could arise if we consider sets of measurements for different materials of strongly different amount of statistics. However, one can prove that as long as the time dependency of the error follows the same functional behaviour for those sets (i.e., the errors differ just in a scale factor), the re-binning will not be affected. Generically, we expect that the overall envelope of the number of counts is set by the exponential decay of muons, which is set by the universal muon lifetime, and material-specific details will provide sub-dominant changes to this overarching behaviour.

We wrote implementation of PCA for data from experiment as a package for python. Current version of the code, finished at the time of publishing this article, can be found in [62].