Principal component analysis of diffuse magnetic neutron scattering: a theoretical study

The study of quantum matter has emerged in recent decades as a major field of scientific endeavour. The behaviour of many-body systems is quite well understood at relatively high temperatures where it is dominated by classical forces and entropy and where it can be simulated efficiently using classical computers. However quantum effects including entanglement and particle indistinguishability make the equivalent, low-temperature problem much harder in principle. Even so, a good understanding has emerged for various ordered ground states, including the Landau Fermi liquid and states showing magnetic, superconducting or topological order [1, 2]. Strongly-correlated quantum matter [3], on the other hand, shows quantum correlations persisting at intermediate energy scales and is less well understood with many outstanding questions. For example, the precise relationship between the intermediate-temperature, 'liquid' states and the various ground states in its phase diagrams remains unknown [4, 5].

In recent years the arsenal available to tackle such challenging problems has been enlarged by the application of machine learning (ML). For instance, artificial neural networks have been used to efficiently encode the wave function of a many-body Hamiltonian, searching for the ground state by reinforcement learning [6]; to predict the properties of one material from those of other substances, without involving a model Hamiltonian [7]; and to detect phase transitions from piezoelectric relaxation measurements [8] and spectroscopic imaging scanning tunnelling microscopy [9].

Here we propose an application of ML to magnetic neutron scattering (NS). The neutron's intrinsic magnetic moment and availability of neutron beams with wavelengths of the order of an angstrom make NS one of the most powerful probes of magnetism in materials [10]. NS has provided, for instance: a thorough characterisation of the magnetic excitation spectrum of the cuprates [11–15]; strong evidence of magnetic monopole excitations in 'spin ice' frustrated magnets [16, 17]; and a quantitative understanding of quantum phase transitions in magnetic insulators [18–21].

Our approach to magnetic NS is based on principal component analysis (PCA), a well-established technique for dimensionality reduction [22, 23] that can be regarded as a form of unsupervised ML. Formally, PCA is equivalent to a linear autoencoder [24] and is a suitable initial step for a wide range of classification problems. In recent years PCA and auto-encoders have been applied to data obtained through numerical simulation of many-body systems [25, 26]. It has been shown that, when provided with detailed information on large representative samples of microstates of such systems, such algorithms are capable of 'discovering' important features in their phase diagrams, including order parameters and phase transitions [25]. On the other hand experimentally-accessible information is normally limited and does not provide access to individual microstates, consisting instead of thermal averages. The question emerges: can PCA still identify important features from such averages?

Recently an autoencoder-based approach to magnetic diffuse NS data on the 'spin-ice' material Dy₂Ti₂O₇ has been demonstrated [26]. The autoencoder is trained on a set of simulated neutron-scattering images. The simulations correspond to a class of candidate model Hamiltonians. The trained autoencoder is then used to describe real experimental data. It was found that the autoencoder provides a compact description of the experimental data, facilitating the identification of optimal model Hamiltonians, and that it can also recognise distinct physical regimes in the simulations. The latter suggests the question highlighted above can be answered in the affirmative.

Here we ask whether PCA can be used to infer relevant features from the data even in the absence of prior knowledge of a class of applicable Hamiltonians (or in cases where the Hamiltonians might not be tractable). This would require training the algorithm directly on the experimental data. Given the scarcity of neutron flux, the training set would need to consist of a limited number of scattering images, each with limited resolution. The trained algorithm would then be used to filter additional (but similarly limited) data sets and it would have to produce qualitative signatures of any relevant features. Of particular interest is the ability to infer ground-state phase boundaries from finite-temperature data.

We perform a theoretical study to address the above questions, focusing on a class of models describing spin-1/2, anti-ferromagnetic, ring-shaped molecular magnets. The field-dependent phase diagrams of all instances of our model feature one or more level crossings (LC) where the nature of the ground state changes discontinuously. One of these LCs is a so-called 'entanglement transition' (ET). At the ET the ground state factorises exactly. These changes in the system's ground state have clear signatures in simulated low-temperature, high-resolution diffuse magnetic NS cross-sections [27]. The main question we tackle here is whether a PCA can detect these features using a more limited number of lower-resolution images. We will see that this is indeed possible. Moreover the ground-state phase boundaries can be accurately determined from finite-temperature data. In order to achieve all this we introduce the notion of a 'score bifurcation plot' for principal components (PCs). This tool is only a small variation of the usual 'score plots' in general use with PCA algorithms but it shows in a particularly transparent way how the different quantum ground states in the model emerge from the high-temperature phase as the temperature is lowered. We argue that this can be a useful tool for the identification of quantum ground states from experimental data.

For the purpose of our study we consider a spin-1/2, anisotropic Heisenberg ring in an applied magnetic field perpendicular to the plane of the ring (for an illustration, see figure 1). Assuming nearest-neighbour interactions only, the system has the Hamiltonian

$$\begin{align}\hat{H}& =\sum\limits _{j=1}^{N}\left\{-J\left[\left(1+\gamma \right){\hat{S}}_{j}^{x}{\hat{S}}_{j+1}^{x}+\left(1-\gamma \right){\hat{S}}_{j}^{y}{\hat{S}}_{j+1}^{y}\right.\right.\\ & \quad \left.\left.+{\Delta}{\hat{S}}_{j}^{z}{\hat{S}}_{j+1}^{z}\right]-h{\hat{S}}_{j}^{z}\right\}.\end{align} \tag{ 1 }$$

Here N is the number of magnetic ions in the ring, which we assume to be even, J and h are, respectively, the exchange and field energies, and γ and Δ are two dimensionless parameters describing the anisotropy of the spin–spin interaction. The operator ${\hat{S}}_{j}^{\alpha }$ represents the ${\alpha }^{\underline{\text{th}}}$ component of the spin at the ${j}^{\underline{\text{th}}}$ magnetic site and the labels x, y, z refer to the local magnetic axes at that site. The x and y axes rotate from site to site so as to preserve the C_N rotational symmetry of the molecule around the z axis, which is fixed. Note that we have assumed that the interaction is diagonal in this basis. An illustration of the geometry of the model can be found in reference [27, figure 1]. We take all three components of the interaction to be anti-ferromagnetic, corresponding to J < 0, and assume without loss of generality 0 ⩽ γ ⩽ 1 and Δ > 0. ¹ The boundary conditions are enforced by setting N + 1 ≡ 1.

Zoom In Zoom Out Reset image size

The behaviour of the model defined by equation (1) has been studied extensively [27–32]. For fixed J, γ, Δ it has N/2 ground state degeneracies at h = h₁, h₂, ..., h_N/2 = h_f, where the last degeneracy takes place at the N-independent factorisation field ${h}_{\mathrm{f}}=J\sqrt{{\left(1+{\Delta}\right)}^{2}-{\gamma }^{2}}$. The simulated diffuse magnetic NS function $S\left(\mathbf{q}\right)$ for the geometry under consideration and with the scattering vector q within the plane of the molecule [27] shows qualitative changes from anti-ferromagnetic correlations for h < h_f to ferromagnetic ones for h > h_f. This is consistent with an ET from anti-parallel bell states to parallel bell states, respectively, known to take place at h_f Amico et al [33]. Less striking, but well-defined changes also occur at the other LC. Specifically, numerical evidence for a jump of $S\left(0\right)$ in the ground state taking place at all N/2 LC has been obtained for N = 4, 6 [27], 8 and 10 [34]. At finite temperatures the jumps become crossovers which get smoother as the temperature is raised further.

The codes we used for this study are freely available as open source from reference [35] (NS simulations) and [36] (PCA). They require only the Octave computer language [37]. Further details are given in the appendix A.

For any given set of values of the parameters of our model, the scattering function $S\left(\mathbf{q}\right)$ mentioned above can be interpreted as an image ² . Giving the parameters different values allows us to generate different images which can be subject to PCA. Quite generally, the result of a PCA of any set of images is a complete, orthogonal basis set that can be used to reconstruct exactly, through linear superposition, any of the images in the original training set. The advantage of this new basis is that the PCs are ordered, with the first basis element capturing the largest amount of variance within the original data set, the second capturing the second largest amount, and so forth. For images comprising solely random pixel values this would offer no advantage but if the images are strongly-correlated then a very good approximation to all the images in the training set can be obtained using just the first few PCs. PCA can thus be regarded as a technique for dimensionality-reduction and this forms the basis of its application to problems such as face recognition [23]. Its effectiveness relies on correlations within the training set: if correlation is high (e.g. all images represent human faces) then a small number M of PCs can capture most of the variance in the data set. The details of our PCA procedure are given in the appendix A.

In our problem we expect to achieve significant reduction because all images have been derived from instances of the same class of Hamiltonians. Let us fix the number of magnetic moments in the molecule N and vary the parameters γ, Δ, h and temperature T (all four energies are measured in units of J). For each set of values we can use the method in references [27, 34] to compute $S\left(\mathbf{q}\right)$ for a fixed set of wave vectors q. This results in a set of images which can then be classified by a standard PCA algorithm. Our expectation is borne out by the scree plots in figure 2(a). Specifically, we find that M = 4 PCs suffice to capture 99% of the variance in the data set for the range of values of N shown in the graph.

Zoom In Zoom Out Reset image size

In an experimental situation, we expect the parameters defining the strength and anisotropy of magnetic interactions to be fixed for a given material, while the strength of the externally-applied magnetic field and temperature can vary. Scree plots for a representative case (γ = 0.6, Δ = 0) are shown in figure 2(b). ³ We find that now M = 2 captures 99% of the variance in the training set.

Our results indicate that the number M of PCs necessary to reproduce to very high accuracy all the images in the training set is bound by the number M_p of free parameters used to generate the training images. This might suggest an unbiased (model-independent) way to constrain experimentally the number of independent parameters describing a class of related materials—an important step in the derivation of a model Hamiltonian. We note, however, that M is more likely related to the number of distinct states in the phase diagram which, although related to M_p, can potentially be larger than it. Elucidating this will require the analysis of a broader range of models and is beyond the scope of the present work.

We highlight that unlike an autoencoder, where the number M_n of neurons in the hidden layer has to be fixed a priori, our approach enables us to find out a posteriori the number M of PCs needed to describe the data accurately. One would in principle expect M_n ∼ M and thus for our model M_n ∼ M_p. In contrast, for the model of reference Samarakoon et al [26] with M_p = 4 dimensionless parameters M_n = 30 was found to strike a good balance between overfitting and underfitting. This would suggest that in our case we achieve greater dimensionality reduction.

Although PCA is equivalent to an autoencoder in the linear limit, we must note that the two studies are rather different. In Samarakoon et al [26] the aim is to automatically interpret structure factor data, from noisy measurements and to disambiguate among many possible solutions of the inverse scattering problem, whereas our approach is to achieve an embedding of the data that is amenable to human interpretation. Our work does not explicitly address experimental noise. However, we note that both autoencoder and PCA are known to have noise suppressing properties and therefore we anticipate that our method would also be applicable to experimental data.

Once the value of M has been set, other images not included in the original training set can be accurately reconstructed using the first M PCs. For this to be possible, the following two conditions must be met: (i) the new images must be of the same type as those in the training set; (ii) the training set must be sufficiently representative. The latter criterion will only be met if the number of images in the training set M_t is sufficiently large. This limits the feasibility of obtaining training sets experimentally.

Figure 3 shows three versions of the first 2 PCs for a particular instance of our model (N = 4, γ = 0.6, Δ = 0). Each version has been obtained using a different training set, shown graphically in the first row of panels: a large training set obtained from 500 randomly-chosen values of $\left(h,T\right)$; a much smaller training set generated from 9 randomly-chosen values of $\left(h,T\right)$; and a minimal set formed by 3 values of $\left(h,T\right)$, chosen strategically (one with high T one with low T and low h, and one with low T and high h). ⁴ Appendix B gives further details about the smaller of the two random sets, including a display of the individual simulated NS images that compose it. Our results indicate that M_t can indeed be very small. This is consistent with the small number of PCs necessary to describe the data accurately (see previous section) and is confirmed by figure 4 which shows the reconstruction of a particular instance of $S\left(\mathbf{q}\right)$ using the three sets of PCs. We note that the reconstructed image is not present in any of the training sets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 2–4, taken together, imply that a relatively small number of initial observations can be used to determine a few PCs in terms of which data obtained subsequently can be accurately described—in other words, by projecting new observations onto the low-dimensional space spanned by the PCs, we obtain a low-dimensional representation of the data in terms of the PC scores. The dependencies of such scores on system parameters such as magnetic field h or temperature T can then be used to identify features in the phase diagram. Such approach, when applied to microscopic data on classical states, has been shown capable of detecting, for example, symmetry-breaking phase transitions [38]. Here we ask whether the same benefit can be obtained when working with observable averages such as $S\left(\mathbf{q}\right)$ for our quantum magnet model. We will answer in the affirmative and moreover present a useful analytical tool based on this idea, which we call a 'score bifurcation plot'. Our simulations indicate that this technique may facilitate the detection of qualitative changes in the ground state of a quantum system, even from data taken at relatively high temperatures.

To illustrate our method we consider first the simplest case of our model where N = 2. Such S = 1/2 dimer has two possible ground states: a low-field anti-ferromagnetic state |↑↓⟩ − |↓↑⟩ with anti-parallel entanglement and a high-field ferromagnetic state |↑↑⟩ + δ|↓↓⟩ with parallel entanglement (as h → ∞ the amplitude δ → 0 resulting in the classical state |↑↑⟩). The two states are degenerate at the factorisation field h_f. Panels (a) and (b) of figure 5 show the weights, or scores, of the first two PCs, wPC1,wPC2, obtained by projecting the PCs onto the simulated NS function $S\left(\mathbf{q}\right)$, for two different training sets. In panel (a) the weights have been obtained for the scattering functions in the original training set. We can appreciate a marked difference between the high-temperature data, concentrated in a small region of wPC1–wPC2 space, and the low-temperature data which cover a much wider area. This is reminiscent of results obtained for microstates of Ising-type models—cf figure 3(b) in reference [38], where we also see that, at high temperatures, the data are clustered in a small, high-entropy region of PC space while, at lower temperatures, the data fan out as distinct ground states are selected by their energies. In our case, however, we are examining the statistical average $S\left(\mathbf{q}\right)$ which is a function of parameters h, T that can, in an experimental situation, be controlled externally. It is therefore possible to explore the space of scattering functions systematically by varying h and T and projecting the new measurements onto the PCs deduced from the training set. A simulation of that approach is presented in panel (b), which shows the evolution of wPC1 and wPC2 with temperature for different fixed values of the magnetic field (fixed-field temperature scans). We observe a marked difference between the curves corresponding to fields lower than the factorisation field h_f and those greater than that field. Each constant-field temperature scan is represented by a single curve in the space defined by the two PCs. As the temperature is lowered, the curve starts to bend in a direction that indicates the nature of the ground state (ferromagnetic if h > h_f and anti-ferromagnetic if h < h_f). This manifests as a marked bifurcation, with the curve corresponding to h = h_f (highlighted in cyan) marking the bifurcation where the directions of this bending changes sign. At this field, below some finite temperature T* the system gets 'stuck' at a particular point (wPC1^*, wPC2^*) and does not evolve further. This suggests that such 'score bifurcation plots' can be used to elucidate systematically the ground-state phase diagram from finite-temperature data, even in systems such as the one we model where there are no finite-temperature phase transitions. Appendix D describes a geometric construction that can be used to determine h_f accurately from limited finite-temperature data.

Further evidence of the above hypothesis is provided in figure 5(c) where similar data are presented for N = 4. As we reviewed in section 2, two special values of the field h₁, h_f are expected to emerge at low temperatures in this case. Indeed we find two bifurcations occurring at those fields, within the resolution given by our field increment Δh = 0.05J. We note however that the bifurcation at h_f, where the ground state factorises leading to the vanishing of entanglement measures, is detectable at a higher temperature than that at the level crossing field, where entanglement is suppressed but does not vanish [27]. We have verified that the factorisation field is also seen in similar score bifurcation plots obtained for N = 6.

Zoom In Zoom Out Reset image size

In section 3 we showed that PCA could achieve good dimensionality reduction for NS data simulated using our model. In section 4 we further showed that this could be achieved using surprisingly small training sets, in the sense of containing very few observations. We will now address the question of experimental resolution—specifically, how many pixels each individual observation needs to have for the score bifurcation plots introduced in the last section to yield accurate values of the critical fields. We note that modern NS instruments allow a trade off between neutron flux and resolution at the time when the measurement is being made [39]. In addition, one can always group pixels together, post-measurement, to form a more coarse-grained, but less noisy image. We note that the authors of reference [26] have addressed the question of noise in a different way, namely by applying their methodology directly to noisy images, and their conclusions are similar to ours.

To address our question we have repeated some of our previous calculations using lower-resolution images both at the training and analysis stages. Specifically, we replace our previous 24 × 24 pixel matrices with 8 × 8 matrices which corresponds to nearly an order of magnitude reduction in the amount of data in each observation. The results, for the N = 4 case, are shown in figure 6. Clearly, the score bifurcation plot obtained from these lower-resolution images is as useful as that obtained before, and in particular it allows us to pinpoint the critical fields h₁ and h_f to the same values, within our field-scan accuracy Δh = ±0.025J.

Zoom In Zoom Out Reset image size

In this work we have presented a method to obtain quantitative information about the phase diagram of a quantum system from experimentally observable data, namely the diffuse magnetic NS function $S\left(\mathbf{q}\right)$. Our method is based on a simple form of unsupervised ML, PCA, and provides a visual representation of the data that facilitates a greater understanding of the underlying physics. We addressed our research question theoretically by analysing simulated scattering functions for a simple model of a cluster quantum magnet.

Our method is based on using a small training set of observations to determine PCs describing the data. We then analyse subsequent observations by projection onto those PCs. We found that this procedure can achieve a large degree of dimensionality reduction i.e. very few PCs suffice for an accurate description of subsequent observations. Consequently, we found that effective training requires only a small number of observations. Moreover, each observation may be a very low-resolution image, which should facilitate greatly the implementation of our methodology in a real experimental setting. Finally, we devised a way to use the trained PCA algorithm to characterise the system's phase diagram.

Our method to investigate phase diagrams relies on temperature scans at fixed values of another control parameter (in our case, magnetic field h—however the method can be trivially generalised to other control parameters such as pressure). By plotting the evolution of the system in PC space we find paths with bifurcations occurring at the values of the field that are known to correspond to changes in the system's ground state.

It is interesting to speculate why a linear technique, PCA, can provide such a compact and enlightening description of a system dominated by strong correlations. In this respect we note that our calculation of the scattering function S(q) was carried out in the linear-response regime [27, 35]. This is standard in the theory of magnetic NS and is justified by the fact that the neutron acts as a weak perturbation [10]. Whether this can be used as the starting point for a justification of the applicability of PCA is an interesting question but lies outside the scope of the present work.

We hasten to add that the linear dependence of S(q) on the PCs does not imply that S(q) is a linear function of the model parameters γ, Δ, T, h. The scattering function is given by

$$\begin{equation}S\left(\mathbf{q}\right)=\sum\limits _{n}\enspace {w}_{n}\left(\gamma ,{\Delta},h,T\right){S}_{n}\left(\mathbf{q}\right)\end{equation} \tag{ 2 }$$

where the score w_n (γ, Δ, h, T) of the ${n}^{\underline{\text{th}}}$ PC is in general a highly non-linear function of γ, Δ, h and T. Thus in our dimensionally-compressed representation of the data S(q) is a linear function of just a few S_n (q)'s but the non-linearity of the model is still reflected in the dependence on system parameters of the expansion coefficients. It is nevertheless remarkable that the non-linear features of the model can be accurately captured in this way. An explicit example of such non-linear behaviour of the PC scores is provided in appendix C.

Our method provides an efficient way to extrapolate the ground-state phase diagram from finite-temperature data, even in a system that is effectively of finite size and therefore lacks well-defined, finite-temperature phase boundaries. One may speculate that applying such methodology to poorly-understood real systems such as copper-based high-temperature superconductors might offer a fresh perspective on an old conundrum in the theory of strongly-correlated electron systems: are the correlated quantum 'liquid' phases found at finite temperature best thought of as manifestations of the quantum critical points separating the different quantum-ordered phases (as proposed by Laughlin and co-workers [4])? Or are they best regarded, instead, as condensations of the 'gaseous' phase existing at higher temperatures, which become susceptible to different quantum ordering transitions as the temperature is lowered further (as put forward by Anderson [5])? Addressing this question will require applying our method to experimental data on real systems. For instance, a score bifurcation plot of measurements taken in the 'strange metal' phase of the cuprate phase diagram may provide evidence for proposals that puzzling crossovers observed in that region are due to a Widom line associated with a hidden, low-temperature critical endpoint [40].

JQ wishes to acknowledge useful discussions with S T Carr, G Möller, S Ramos, and T Tula and financial support from the Engineering and Physical Sciences Research Council (Grant No. EP/P00749X/1).

No new data were created or analysed in this study.

The PCA was implemented in the Octave (v3.4.3) programming language [37] using the singular value decomposition function which returns singular vectors normalized to unit length. A is a matrix, in which each column is a scattering function image represented in column-wise concatenated form. The data is centred by subtracting from each column its mean pixel intensity value, thereby forming matrix X. This is followed singular value decomposition from which PCs, PC scores, and scree plots are straight-forward to obtain. We reproduce the key part of our code here:

• (a)  
  : % Data centering:
• (b)  
  : X = A-ones(size(A)(1), 1)^*mean(A)
• (c)  
  : % Singular value decomposition:
• (d)  
  : [V, lambda, junk] = svd(X'^*X);
• (e)  
  : % Principal components:
• (f)  
  : U = X^*V^*lambdaˆ(−.5);
• (g)  
  : % Principal component scores:
• (h)  
  : S = U'^*X;
• (i)  
  : % Data for scree plot:
• (j)  
  : scree_data = diag(lambda);

Adapt the code for Matlab by replacing line 2 with X = A-ones(size(A, 1), 1)'^*mean(A);

The complete code is available in reference [36] as a Github repository. The codes used to generate the matrix A can be found in reference [35].

We note that our re-centering procedure does not involve subtracting an average over observations, as is common in other PCA implementations. This essentially means that the first PC describes that average, and the corresponding score quantifies how much a specific observation deviates from it. Our tests indicate that, for the type of data studied here, this yields a clearer and more complete description of the underlying correlations.

Here we give details of one of the training data sets, namely the nine-sample data set corresponding to N = 4, γ = 0.6, Δ = 0 and the values of temperature T and magnetic field h shown in figure 3(b). Table 1 shows the numerical values of h and T for this set, along with an index a, b, c, ... used to identify each particular (h, T) combination. Figure 7 shows the corresponding simulated NS functions. Each image constitutes one element of the training set (an 'observation'). This nine-image training set was used to obtain the PCs shown in figures 3(e) and (h). Projecting the 9 images shown in figure 7 onto those PCs we obtain the coordinates of the 9 filled circles shown in figure 5(b). The lines shown in the latter figure are obtained by projecting other images (not forming part of the training set) onto the same PCs. A similar procedure was followed when using the three-image training set corresponding to the (h, T) values in figure 3(a) and the 500-image training set corresponding to figure 3(c).

Zoom In Zoom Out Reset image size

Here we discuss briefly the field- and temperature-dependences of the individual PC scores wPC1, wPC2. These are shown in figure 8 for the same calculation used to obtain figure 5(b). As we can see the temperature-dependence of the PC scores changes qualitatively at the factorisation field—for instance at low temperatures the derivative of wPC1 with respect to T is positive for h > h_f and negative for h < h_f. At exactly h = h_f (highlighted in cyan) this derivative vanishes in the limit T → 0. We note, however, that in a general situation we will not have a constant-field scan taken at exactly h = h_f (see footnote 3 on page 4). Because of this, in general locating the value of the field where the PC scores become temperature-independent involves an extrapolation. This must be contrasted with the score bifurcation plot where wPC1 and wPC2 are shown in combination and the factorisation field can be determined precisely from finite-temperature data.

Zoom In Zoom Out Reset image size

Here we describe a simple geometric construction for determining precisely the location of the factorisation field from finite-temperature data. We illustrate this by focusing on the N = 2 data in figure 5. The same data is reproduced in figure 9, but with additional symbols showing the points along the constant-field curves where the temperature reaches two particular values, namely T = 0.10 and T = 0.15. By inspection of that plot it is quite clear that the two curves will cross at the point where the bifurcation takes place, whose precise coincidence with the factorisation field we have determined numerically. It is evident that this construction will work, and give the same value, irrespective of which two temperatures we choose, as long as they are low enough to show the inflection of the curve to the right of the bifurcation. A similar construction can be made for the N = 4 plots in figures 9(c) and 6.

Zoom In Zoom Out Reset image size