Reconstruction and uncertainty quantification of lattice Hamiltonian model parameters from observations of microscopic degrees of freedom

Phase transitions between dissimilar phases underpin multiple areas of condensed matter physics and materials science. From an applications viewpoint, phase transitions are intricately involved in the process of material formation in virtually all technologies from ceramics to metals to semiconductors. Electrical, gating, and chemical control of phase transitions have emerged as a promising pathway for low-voltage electronics in tunneling barriers and field effect devices.^1,2 Beyond condensed matter physics, phase transitions in Bose–Einstein condensates allow new classes of quantum information technology devices.³ 

The ubiquity of the phase transitions in virtually all areas of applied and fundamental science has made them one of the central areas for theoretical exploration. Among the existing paradigms for exploring phase transitions, a special place is held by the lattice models such as Ising,⁴ Heisenberg,⁴ Kitaev,⁵ etc. Here, the system is represented as a set of interacting spins on a spatial lattice. The spins represent the individual degrees of freedom of the material and can be magnetic spins, electrical dipoles, atomic species, etc. The lattice defines the geometry of the interactions. Depending on the conservation laws, character of the spins, interactions, and geometry, these lattice models can represent an extremely broad range of phenomena from ordering in metal alloys to ferroic transitions to surface adsorbate dynamics.

The universality of the lattice model approach has spurred intensive analytic and numerical investigations of these models.^6,7 These studies mainly seek to explore the phase diagrams of the lattice models, i.e., correspondence between the global variables such as temperature and field and preponderant behaviors and patterns for spin arrangement. Similarly of interest are global thermodynamics parameters such as magnetization, heat capacity, and susceptibility. The spatial arrangements of the spins are traditionally explored in the form of the ground-state pattern as well as correlation functions representing the distributions of the defects and disorder, i.e., over which length scales the system is ordered.⁸ Notably, these behaviors allow for straightforward exploration only for systems with well-defined ground states, whereas the presence of competing interactions leads to frustrated and poorly defined ground states. Nonetheless, the combination of efficient sampling strategies and high-performance computing has allowed exploration of even complex multidimensional models.

Over recent years, the emergence of machine learning (ML) tools and ready availability of information-compression algorithms, such as principal and independent component analysis, support vector machines and deep learning tools such as variational autoencoders, have provided new impetus within the field.^9–12 Even though the applications of ML tools to lattice models are relatively recent, the field has experienced an exponential growth in the last three years, opening pathways to exploration of exotic quantum phases, off-lattice models, etc.^13–16 

However, experimental selection of the appropriate lattice model corresponding to a specific phase transition and determination of its parameters, termed the “inverse Ising problem,” presents a more difficult problem;^17–20 an in-depth view of that field is explored in an article by Nguyen et al.²¹ Traditionally, the symmetry and type of order parameter are determined based on macroscopic measurements and scattering studies, whereas numerical values of exchange integrals are obtained based on the matching between theoretical predictions and macroscopic measurements of thermodynamic parameters. However, these studies typically require high-quality samples for which the corresponding critical behaviors can be reliably determined. Similarly, in the materials with high disorder, non-ergodic or non-thermalized ground states determination of appropriate lattice models can be highly non-trivial.

The development of high-resolution imaging tools such as scanning transmission electron microscopy and scanning probe microscopy has allowed direct insight into the atomic-scale structure of materials.^22–24 Beyond atomically resolved and element specific imaging, these techniques provide information on the minute symmetry-breaking distortions, allowing for direct mapping of order parameter fields such as polarization,^25–29 octahedra tilts,^30–33 strains,³⁴ and more complex emergent phenomena.³⁵ Correspondingly, of interest is the extraction of the parameters of the underpinning physics model from the atomically resolved observations. As one approach, this can be achieved using the direct matching of the mesoscopic order parameter field to Ginzburg–Landau type models.^36,37 Alternatively, the mesoscopic descriptors such as hysteresis loops or relaxation curves can be used as proxy identifiers for local behavior.³⁸ Finally, it has been proposed that the direct observation of the mesoscopic degrees of freedom can be directly compared to the lattice model via statistical distance minimization.^39–43 Indeed, naively expecting to minimize loss between theoretical and experimental datasets without the appropriate metric for thermodynamic systems is incorrect and will lead to incorrect reconstructions.⁴⁴ Previously, we have reported the principles of the statistical distance minimization^41,42 and applications for the establishing of the parameters of the non-ideal solid solutions in layered superconductors,⁴³ information fusion between the surface and depth-resolved chemical information in manganites,^41,43 and pair- and triple-atom model comparison for segregation in layered chalcogenides.⁴⁰ Here, we explore the reconstruction of the Ising model parameters in the presence of competing interactions with statistical distance minimization and explore the sensitivity of the reconstruction to the incomplete knowledge of the model as well as deriving uncertainty bounds on the reconstructions. The guidelines for experimentally driven model selection are proposed.

We explore the behavior of the Ising model with the nearest neighbor and next nearest neighbor interactions. Depending on the sign and magnitude of the exchange integrals, this model can give rise to a rich phase diagram containing ferromagnetic, antiferromagnetic, paramagnetic, and frustrated phases. Here, first we apply the ML tools for rapid simulation of the associate phase diagram based on the local statistics.

The classical Ising Hamiltonian model realized on N²-evenly spaced lattice sites of a square lattice with nearest neighbor (NN) and next nearest neighbor (NNN) interactions is the chosen model for this study. The Hamiltonian of a given configuration σ in the absence of an external magnetic field is given by

where J₁ and J₂ are the exchange integrals corresponding to the nearest neighbor interactions and next nearest neighbor interactions, respectively; ΔJ₁ and ΔJ₂ are the corresponding disorders in exchange integrals; and ⟨i,j⟩ and ⟨i,k⟩ are the sums over all nearest neighbors and next nearest neighbors, respectively.

A Monte Carlo (MC) simulation was performed on a square lattice with N = 40 using the Metropolis algorithm. The system was equilibrated for $1000×N2$ Monte Carlo steps, and the data are acquired on the next $1000×N2$ Monte Carlo steps. A spin configuration reversal of each lattice site was attempted at every step. The flip was accepted when the energy of the resultant configuration was lower than the previous step. Else, the probability of the spin flip was estimated using the Boltzmann distribution given by Eq. (2), where $β=1/kBT$ is the inverse temperature, with T being the absolute temperature and k_B being the Boltzmann constant, and the denominator is the partition function,

Macroscopic properties as a function of temperature and a function of exchange integrals are explored. Setting J₂ and ΔJ₂ to zero gives the traditional 2D Ising model with only the nearest neighbor (NN) interactions. Macroscopic properties of such model as a function of reduced temperature are shown in Figs. 1(a)–1(d). Clearly, all four properties delineate the phase transition from the ferromagnetism to the paramagnetic region. Phase transitions as a function of exchange integrals at a constant temperature are shown in Figs. 1(e) and 1(f). Energy per site is illustrated in Fig. 1(g). The corresponding magnetization is shown in Fig. 1(e), clearly highlighting the ferromagnetic region and transition between the ferromagnetic and the paramagnetic states (diagonal line). The specific heat is shown in Fig. 1(h), delineating the primary phase transitions between ferromagnetic, antiferromagnetic, paramagnetic, and frustrated states. Finally, susceptibility is shown in Fig. 1(f), showing the divergence across ferromagnetic–paramagnetic phase transition. Overall, Fig. 1 provides a traditional approach for mapping the phase diagram of lattice models based on the macroscopic thermodynamic parameter analysis.

An alternative approach for mapping phase diagrams in simulated lattice models is based on the machine learning and statistical analysis of the local spin configurations and correlation between spins. For example, multiple works showed that using a supervised learning techniques such as convolutional neural networks, which are traditionally used for image recognition tasks, one can identify phase transition and map phase boundaries in the classical and quantum lattice models using two-dimensional and three-dimensional “images” of the spin configurations around the critical points as a training set.^9–11 At the same time, the more conventional multivariate statistics and probabilistic learning tools such as the kernelized principal component analysis and the variational autoencoder allowed learning phase transition and the associated order parameter in the classical Ising model from local spin configurations in an unsupervised manner.^45,46

Traditional descriptors of such microstates are spin configurations of all the lattice points.⁴⁷ In this scenario, the number of rows (R) in the feature matrix is R = N_sim* N_mc, where N_sim is the number of distinct sets of simulation parameters and N_mc is the number of Monte Carlo steps considered for the analysis. The columns (C) are the number of features used for the analysis. Spin configurations of the entire lattice at each simulation are used as attributes. The first few principal components corresponding to such matrix closely resemble the macroscopic properties’ (magnetization, energy, specific heat, and susceptibility) trends as a function of reduced temperature.⁴⁷ 

Here, we propose the approach where phase classification is performed based on the relative frequencies of local configurations at a given set of simulation parameters. These frequencies of local configurations describe the microstates of the system exhaustively. For a given case, these are averaged over all the microstates generated by the Monte Carlo simulations. The local configurations considered for further analysis are described in Fig. 2(e). The relative frequencies of local configurations of ferromagnetic, anti-ferromagnetic, paramagnetic, and the frustrated systems are shown in Figs. 2(a)–2(d). In this case, the number of elements of the feature matrix is 5400 as opposed to 10⁹ in the study referred. Monte Carlo simulations are performed at a reduced temperature of T_r = 0.96 and at 900 evenly spaced points on a grid space of exchange integrals J₁ = [−1.25, 1.25] and J₂ = [−0.75, 0.75] in the units of k_BT. From here on, all the values of exchange integrals are reported in the units of k_BT. Phase diagrams of complex systems can be mapped using the statistical analysis of the local structure values.

To establish the phase diagram of the system over the parameter space, we use a simple k-means clustering approach. A similar analysis was proposed by Canabarro et al. in which the pairwise correlations among all spins are used as the descriptors.⁴⁸ This algorithm segregates the input data into k-clusters based on minimizing the variance within each cluster.⁴⁹ The number of clusters is apriori unknown but can be estimated based on the quality of separation and the target output. Here, when the feature matrix is segregated into four clusters using the k-means clustering technique, the resulting clusters [Fig. 3(a)] exactly resembled the energy [Fig. 1(a)] and specific heat [Fig. 1(c)] phase diagram. In Fig. 3(a), phases 1, 2, 3, and 4 represent the ferromagnetic, anti-ferromagnetic, frustrated, and paramagnetic regions, respectively. Representative configurations for the four phases are plotted in Fig. 3(d). Ferromagnetic and anti-ferromagnetic transitions to paramagnetic phase can clearly be observed in Fig. 3(a). It is also observed that for k = 2, the clustering plot reproduced the magnetization and susceptibility phase diagrams. The clusters are also plotted in two-dimensional space of first two principal components⁴⁹ [Fig. 3(b)] and as a dendrogram [Fig. 3(c)] to show the hierarchy between the clusters. Figure 3(b) illustrates the separability of the phases and their general configurations in reduced (principal component) space. This provides information on which phases are separated by continuous and discontinuous transitions. Dendrogram [Fig. 3(b)] shows the intra-cluster distances between different clusters in the “histogram space.” The subsequent increase in the clusters’ number leads to the emergence of additional regions concentrated at the boundary between frustrated and paramagnetic states; however, primary ferromagnetic, antiferromagnetic, and frustrated regions remain invariant.

We note that classically the exploration of the lattice models relies on the macroscopic thermodynamic descriptors such as those shown in Fig. 1 as well as structure factors and correlation functions. Here, we rapidly map the relevant phase diagram using the comparison of the local motifs. In addition to the obvious and previously reported applications for the exploration of the phase diagrams of the lattice models, these observations suggest that the statistics of local spin configurations contain the information on the microscopic interactions in Eq. (1).

Here, we extend the statistical distance minimization method to reconstruct the interaction in the Ising model within the full parameter space. During reconstruction, we have access to the spin configurations (synthetic experimental data) but not the values of Ising model parameters (exchange integrals). Correspondingly, the goal of the reconstruction is to determine the exchange integrals (parameters of the generative physics model) based on the observations of local configurations and establish the associated uncertainties.

To determine the unknown model parameters, we run the Monte Carlo simulations at a given set of simulation parameters to generate the spin configurations of the base (pseudo-experimental or synthetic experimental data) case. Here, 20 snapshots from the last 400 steps of the MC simulations are used as experimental observables for the base case. To reconstruct the exchange integrals of a given base case, simulations were carried out on the entire grid space of the exchange integrals. The local configurations of the base case are then compared to the simulated cases using the statistical distance metric.⁵⁰ This metric is employed to quantify the distinction between pair thermodynamics systems and is given by

where p_i and q_i are the probabilities of finding a local configuration i in the base case and model simulations, respectively, and the summation runs over all possible (six in this case) local configurations. The probabilities of the local configurations are determined as the maximum likelihood estimates from the experimental (p) and simulated (q) histograms. Statistical hypothesis testing indicates that as the measured distance metric reduces to zero, the target and model system measurements become indistinguishable. A unique characteristic of this metric, which sets it apart from popular measures such as Kullback–Leibler divergence, is that it automatically incorporates the statistical weights of the collected data, thus allowing for separation of weak signals (viz., thermal fluctuations) from the statistical noise inherent in samples of a thermodynamic system. Extracting these fluctuations from the data is critical for improving the predictive capabilities of the optimized model, as they encode the system's response to changes in external conditions.

Since sampling in any real experiment is limited, the key question of exactly how many samples are required to achieve a certain tolerance on the inferred parameters is an important one. Toward this aim, we consider an averaging method utilizing the samples at hand. In this case, we divided the parent dataset of 400 images (last 400 MC steps of synthetic experimental data) into 20 sets of uncorrelated images. Interaction parameters are reconstructed for each subset at a given reduced temperature, and mean, standard deviation, and the probability density function of the reconstructed values are reported.

First, we explore the reconstruction of the model parameters for the full model. The reconstruction results for all four quadrants of the phase diagram are shown in Figs. 4 and 5. For all values of J₁ and J₂, the reconstruction shows that the statistical distance minimization yields an unbiased reconstruction of the model parameters. Here, two base cases (J₁= + 0.75, J₂ = + 0.35 and J₁ = − 0.75, J₂ = + 0.35) were reconstructed that have a well-defined ground state where the nearest and next nearest neighbor interactions are not frustrated. The uncertainty is quantified using the experimental averaging technique discussed above.

The reconstruction fails below the phase transition region as there is no information in the spin configurations. Exchange integrals were reconstructed with high certainty at temperatures just above the phase transitions, where the configurations are strictly defined by the interaction parameters. The uncertainty of reconstruction increases with increasing temperature due to random spin configurations. Still, the reconstruction is possible for temperatures that are an order greater than the transition temperature before the uncertainty impedes the interpretation.

For the next two cases (J₁ = − 0.75, J₂ = − 0.35 and J₁ = + 0.75, J₂ = − 0.35), the nearest and next nearest neighbor interactions are incompatible and lead to the frustrated ground states. The ground states of these systems are not well defined. Proposed reconstruction technique works even in this scenario (of no well-defined ground states). Since there are no apparent phase-transitions in these regions, the reconstruction works even at low temperatures (Fig. 5). High temperature behavior for these cases shows similar behavior to that described for the first two cases.

To complement the analysis of the uncertainties via the corresponding point estimates of mean and dispersion and the marginalized posterior distribution for individual exchange integral values, we also analyze corresponding joint distributions to perform uncertainty quantification (UQ) in the 2D parameter space. In the limit of large sample, where each sample corresponds to an individual local configuration, the likelihood function attains the shape of the normal distribution centered around the limiting distribution with variance equal to 1/4 as a result of the variance equalizing transformation built into the probability space with metric s. The log-likelihood of a sample at a distance s from the limiting probability distribution is then proportional to −2ns², where n is the number of samples (individual local configurations).^50,51 As the distance s increases, the log-likelihood decreases, and it does so with a linear dependence on the number of samples n. The likelihood over the exchange integral space is shown in Figs. 5(c) and 5(d) at different points of reduced temperature. The increase in uncertainty with temperature can be visualized as the increase in spread of likelihood.

In this section, we explore the cases where the pseudo-experimental (base) case is guided only by the interactions between the nearest neighbors (J₁ = 0.75) and is reconstructed by simulations incorporating both nearest and next nearest neighbors (J₁ and J₂). This case is equivalent to the traditional 2D Ising model with NN interactions only. The reconstructions follow a similar trend as discussed in the previous cases [Fig. 5(a)]. It works with very high certainty post the transition temperature and the uncertainty increases with the temperature. Interestingly, for this case, the zero value of J₂ can be reliably determined.

We further explore the reconstruction of an underdetermined model [Fig. 6(b)] where the experimental image is guided by both nearest and next nearest neighbors (J₁ = + 0.75, J₂ = + 0.35) while it is modeled using the simulations with only nearest neighbor interactions (J₁). The model attempts to fit the data by overestimating the exchange integral, but the estimated value is different at different temperatures. It can be inferred from the reconstruction curves that a model with only the nearest neighbor interactions will not be able to reproduce the experimental results that are governed by both nearest and next-nearest neighbors. The likelihood for both these cases in the exchange integral space is shown in Figs. 7(a) and 7(b), respectively.

Reconstruction until this section has been performed on a square lattice in which the values of exchange integrals in the pseudo-experimental dataset were constant over the entire lattice. In this section, we explore the reconstruction when a bond-disorder is introduced. To execute this, every lattice site is associated with a different value of exchange integral J + ΔJ, where ΔJ is normally distributed with zero mean and a specified standard deviation. The disorder is added to both exchange integrals J₁ and J₂. The standard deviation is referred to as disorder for the rest of the section. The bond disorder was only introduced in the pseudo-experiment case while the simulation cases were run at ΔJ = 0.

Reconstructions are performed for disorder levels of 10%, 100%, and 1000% of the values of exchange integrals, and the results are shown in Fig. 8. The reconstruction works well at 10%. With 100% disorder, the reconstruction starts to fail at low temperatures since at these temperatures the values of exchange integrals have a very strong effect on configuration. The configurations get more random with the increase in temperature, and this is represented by the increase in uncertainty in reconstructions. At 1000% disorder, the reconstruction completely fails, and the proposed methods cannot be applied at these high levels of disorder.

The likelihood over the exchange integral space is shown in Figs. 9(a)–9(c) at different points of reduced temperature. Note that for temperatures close to transition temperature and low disorder, the recovered exchange integrals are strongly correlated, suggesting that the system behavior can be well defined by single values of exchange integrals. For higher disorder level, reconstruction is impossible.

Above the transition, the reconstruction yields unbiased estimates of exchange integrals for low and intermediate disorder, whereas for high disorder the reconstructed integrals are below the true value. This trend continues for higher temperatures. Finally, in all cases, the temperature increase yields broader distribution of posterior values, rendering reconstruction less reliable. Furthermore, the effects of lattice size on reconstruction and uncertainty quantification are discussed in the supplementary material.

To summarize, here we have developed an approach for the reconstruction of the microscopic parameters of the lattice model with two competing interactions from the observations of the microscopic degrees of freedom. We established the uncertainties and reliability of such analysis in a broad parameter-temperature space. Below the phase transition, the reconstruction is impossible and can lead to the biased estimates of the parameters. However, these situations can be readily identified from the experiment as the presence of large single-phase regions and correspondingly extremely poor sampling of the configuration space. At the same time, the reconstruction is possible well above the phase transition (1–2 orders of magnitude) and in the regions of the parameter space when the macroscopic ground state of the system is poorly defined due to frustrated interactions. Similarly, the reconstruction is robust with respect to the frozen disorder. Interestingly, bond disorder tends to affect reconstruction close to transition temperatures leading to biased estimates of exchange integrals, whereas at high temperatures unbiased estimates with large uncertainty emerge.

As an ancillary task, we have developed a machine learning approach based on histogram clustering to predict phase diagrams efficiently using a reduced descriptor space. We illustrated that the phase evolution in the vicinity of regions comprising frustrated phase gives rise to the hierarchy of local configurations.

The explored approach can be applied to traditionally complex problems in condensed matter physics such as ferroelectric relaxors and morphotropic phase boundary systems, spin and cluster glasses, and quantum systems once the local descriptors linked to the relevant physical behaviors are known. Correspondingly, of interest becomes the role of precision in the determination of local descriptors on parameter reconstruction, necessitating the development of appropriate Bayesian frameworks to link the experimental data and reconstructions.

See the supplementary material for discussion of the effect of lattice size on reconstruction and uncertainty quantification.

This work was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), Materials Sciences and Engineering Division (M.V., L.V., S.V.K., and R.K.V.). A portion of this research was performed and partially supported (M.Z.) at the Center for Nanophase Materials Sciences (CNMS), a U.S. Department of Energy Office of Science User Facility.

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