Data-driven determination of the spin Hamiltonian parameters and their uncertainties: The case of the zigzag-chain compound KCu 4 P 3 O 12
Ryo Tamura, Koji Hukushima, Akira Matsuo, Koichi Kindo, Masashi Hase
DData-driven determination of the spin Hamiltonian parameters and theiruncertainties:The case of the zigzag-chain compound KCu P O Ryo Tamura,
1, 2, 3, ∗ Koji Hukushima,
4, 5, 2, † Akira Matsuo, Koichi Kindo, and Masashi Hase ‡ International Center for Materials Nanoarchitectonics (WPI-MANA),National Institute for Materials Science, Ibaraki 305-0044, Japan Research and Services Division of Materials Data and Integrated System,National Institute for Materials Science, Ibaraki 305-0047, Japan Graduate school of Frontier Sciences, The University of Tokyo, Chiba 277-8568, Japan Komaba Institute for Science, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan Department of Basic Science, Graduate School of Arts and Sciences,The University of Tokyo, Komaba, Meguro, Tokyo 153-8902, Japan The Institute for Solid State Physics, The University of Tokyo, Kashiwa, Chiba 277-8581, Japan Research Center for Advanced Measurement and Characterization,National Institute for Materials Science, Ibaraki 305-0047, Japan (Dated: June 16, 2020)We propose a data-driven technique to estimate the spin Hamiltonian, including uncertainty, frommultiple physical quantities. Using our technique, an effective model of KCu P O is determinedfrom the experimentally observed magnetic susceptibility and magnetization curves with varioustemperatures under high magnetic fields. An effective model, which is the quantum Heisenbergmodel on a zigzag chain with eight spins having J = − . ± .
51 meV, J = − . ± .
13 meV, J = − . ± .
15 meV, and J = 6 . ± .
95 meV, describes these measured results well. Theseuncertainties are successfully determined by the noise estimation. The relations among the estimatedmagnetic interactions or physical quantities are also discussed. The obtained effective model isuseful to predict hard-to-measure properties such as spin gap, spin configuration at the groundstate, magnetic specific heat, and magnetic entropy.
I. INTRODUCTION
An effective model in materials science often explainsthe origin of physical properties in materials. Manymethods have been developed to construct an effec-tive model for a target material, and they can be di-vided into two groups. One is ab initio calculations,which determine the model parameters in an assumedeffective model by providing only basic information ofthe target material . The other is a data-driven ap-proach in which model parameters are determined so asto fit the experimentally measured data in the targetmaterial .In the latter case, many trials and errors are con-ducted to find the appropriate model parameters to de-scribe the experimental results. Data-driven analysesbased on machine learning have been extensively ex-ploited to avoid this cumbersome task. Recently, data-driven techniques are becoming indispensable in mate-rials science, because they should accelerate the dis-covery of novel materials and deepen our under-standing of materials . From the viewpoint of effec-tive model estimations, data-driven techniques are alsoefficient to accelerate automatic searches for appropri-ate model parameters and to extract relevant modelparameters .The paper estimates a spin Hamiltonian as an effec-tive model of KCu P O by a data-driven approach.Figure 1 shows the crystal structure of KCu P O ( a =7.433 ˚A, b = 7.839 ˚A, c = 9.464 ˚A, α = 108.28 ◦ , β = 112.68 ◦ , γ = 92.73 ◦ , and space group: P¯1) . Cu(II)ions have S = 1 / P O to beestimated. We determine the superexchange interactionsbetween Cu ions with uncertainty in this target modelby a data-driven approach in which the experimentallymeasured susceptibility and magnetization curves are in- OCuPK
FIG. 1. (Color online) Crystal structure of KCu P O drawn by VESTA . Black box indicates a unit cell. a r X i v : . [ c ond - m a t . m t r l - s c i ] J un putted. Once an estimated spin Hamiltonian is estab-lished, theoretical analysis of the Hamiltonian predictsvarious magnetic properties, which cannot or have notbeen measured. These properties include the magneticspecific heat, magnetic entropy, spin configuration, andspin gap. These predictions are helpful to propose a fur-ther experimental plan and design.The rest of the paper is organized as follows. SectionII shows the experimental results of the magnetic suscep-tibility and magnetization curves as functions of temper-atures in KCu P O along with the experimental meth-ods. Section III explains our data-driven approach toestimate a spin Hamiltonian. The posterior distributionof model parameters given in the data is constructed bya statistical noise model with respect to the experimen-tal observations and the prior distribution of the modelparameters. Plausible model parameters are determinedby the maximizer of the posterior distribution. Further-more, a systematic method, which allows the statisti-cal uncertainty of the model parameters to be evaluated,is presented to estimate the amplitude of noise in thenoise model. Under our formulation, multiple types ofphysical quantities can be used to estimate the effectivemodel. Section IV explains the results of the spin Hamil-tonian estimation with uncertainty by our data-drivenapproach. First, we assume the shape of the target spinHamiltonian for KCu P O from the viewpoint of thecrystal structure. Next, since the L regularization isadopted as a prior distribution to suppress an increase inthe absolute values of the magnetic interactions, determi-nation of the hyperparameter in the L regularization isdiscussed. Subsequently, by considering the observationnoise, four types of magnetic interactions are estimatedwith the error bars in a target spin Hamiltonian, andtheir relationships are discussed from the distributionsof sampling data by the Monte Carlo method. Finally,various magnetic properties of KCu P O are predicted.Section V presents the discussion and summary. II. EXPERIMENTALLY MEASUREDMAGNETIC PROPERTIES
Figure 2 (a) shows the temperature dependence ofthe magnetic susceptibility with a magnetic field of 0.01T, which was measured by a superconducting quan-tum interference device magnetometer, magnetic prop-erty measurement system (Quantum Design). Crys-talline KCu P O powder was synthesized by a solid-state reaction. Even at sufficiently low temperatures, afinite constant value of susceptibility remains. The insetof Fig. 2 (a) shows the susceptibility upon removing thisconstant term. This data is used in the data-driven ap-proach. Figure 2 (b) shows the magnetization curves attemperatures of 1.3 K, 4.2 K, 20 K, 30 K, and 50 K. Thesemagnetization curves were measured using an inductionmethod with a multilayer pulsed field magnet installedat the Institute for Solid State Physics, the University of (b) (a) FIG. 2. (Color online) (a) Magnetic susceptibility with 0.01T for KCu P O . Inset shows the results where the con-stant term in the susceptibility is removed. (b) Magnetizationcurves at various temperatures for KCu P O . Tokyo. Although a high magnetic field is imposed ( ≤ g -factor of Cu ions is determined tobe 2.08. III. DATA-DRIVEN APPROACH TO ESTIMATEAN EFFECTIVE MODEL WITH UNCERTAINTYA. Posterior distribution for effective modelestimation
Our developed effective model estimation method inRef. 39 is based on the Bayesian statistics. We con-sider that the target Hamiltonian to be estimated has K -types of model parameters: x = ( x , ..., x K ). Let y ex = { y ex ( g l ) } l =1 ,...,L be the set of experimentally mea-sured physical quantities. This set depends on the ex-ternal parameter g l , where the number of data is L . Us-ing Bayes’ theorem, the posterior distribution P ( x | y ex ),which is the conditional probability of the model param-eters given experimental data, is expressed as P ( x | y ex ) ∝ P ( y ex | x ) P ( x ) , (1)where P ( x ) is the prior distribution of the model param-eters and P ( y ex | x ) is the likelihood function of y ex given x . Assuming that the observation noise for the measure-ments follows a Gaussian distribution with a mean of zeroand a standard deviation of σ , the likelihood function isgiven by P ( y ex | x )= (cid:18) πσ (cid:19) L exp (cid:34) − σ L (cid:88) l =1 (cid:0) y ex ( g l ) − y cal ( g l , x ) (cid:1) (cid:35) . (2)Here, { y cal ( g l , x ) } l =1 , ··· ,L , which is expressed as y cal ( x ),are the physical quantities calculated from the Hamilto-nian at the external parameter g l theoretically.To express the posterior distribution briefly, we intro-duce the “energy function” as a function of x and thenoise σ is given by E ( x , σ ) = ∆( x , σ ) − log P ( x ) , (3)where an error function ∆( x , σ ) is∆( x , σ ) = 12 σ L (cid:88) l =1 (cid:0) y ex ( g l ) − y cal ( g l , x ) (cid:1) . (4)Then, the posterior distribution is expressed as P ( x | y ex ) ∝ (cid:18) πσ (cid:19) L exp [ − E ( x , σ )] . (5)From the viewpoint of the maximum a posterior (MAP)estimation, plausible model parameters to explain y ex are regarded as the maximizer of Eq. (5). Correspond-ingly, the MAP estimation is reduced to a minimiza-tion problem of the energy function. To search for themaximizer of Eq. (5) or the minimizer of Eq. (3), vari-ous optimization techniques such as the steepest-descentmethod, Markov chain Monte Carlo (MCMC) method ,and Bayesian optimization can be used.Selecting the prior distribution of the model parame-ters P ( x ) is important to obtain a physically appropriateHamiltonian . The prior distribution for the effectivemodel estimation can be regarded as the regularizationterms in the minimization problem . The most commonare L (LASSO) and L (ridge) regularization with cor-responding prior distributions P ( x ) = exp( − λ | x | ) and P ( x ) = exp( − λ (cid:107) x (cid:107) ), respectively, where hyperparam-eter λ determines the strength of regularization. If theL regularization is applied, model parameters with largecontributions can be selected based on the feature selec-tion. On the other hand, the L regularization suppressesthe increase in the absolute values of the model parame-ters. Depending on the situation and purpose, it is nec-essary to select the proper prior distribution. B. Observation noise estimation
To obtain the uncertainty of the model parameters bya MAP estimation, the width of the posterior distribu-tion around the maximizer must be estimated. The noise amplitude σ is crucial to estimate the uncertainty .In our framework, the noise amplitude σ is considered tobe a hyperparameter and a plausible value is determinedby minimizing the Bayes free-energy F ( σ ) defined as F ( σ ) := − log Z ( σ ) , (6)where Z ( σ ) is the normalization of the posterior distri-bution given by Z ( σ ) = (cid:18) πσ (cid:19) L (cid:90) Ω x d x exp [ − E ( x , σ )] , (7)where Ω x is the support of the posterior distribution de-termined by the prior distribution. To evaluate F ( σ ), it isconvenient to extend the posterior distribution P ( x | y ex )to P β ( x | y ex ) with a “finite temperature”, which is de-fined as P ( x | y ex ) := 1 Z β ( σ ) (cid:18) πσ (cid:19) L exp [ − βE ( x , σ )] , (8)where β is the inverse temperature and Z β ( σ ) is the nor-malization. By using Z β ( σ ), the Bayes free-energy is cal-culated as F ( σ ) = − (cid:90) dβ (cid:18) ddβ log Z β ( σ ) (cid:19) − log Z ( σ )= (cid:90) dβ (cid:104) E ( x , σ ) (cid:105) β + L (cid:0) πσ (cid:1) − log (cid:90) Ω x d x , (9)where the ensemble average (cid:104) E ( x , σ ) (cid:105) β with respect toEq. (8) can be obtained by the MCMC method. Here,the third term in RHS of Eq. (9) does not depend on σ and is omitted below. The noise amplitude of the exper-imental data σ ∗ is evaluated as the value where F ( σ ) isminimized. By MCMC sampling from Eq. (5) with thefixed σ ∗ , the uncertainty of the model parameters can beevaluated. C. Multiple sets of physical quantities
In our experiments of KCu P O , six different setsof physical quantities, that is, the susceptibility and themagnetization curves under different temperatures, areobtained. Then different types of physical measurementsare likely to be combined in our estimation problem.For simplicity, the different physical quantities are as-sumed to be independently obtained. Then the likelihoodfunction of a series of experimental results is defined as P ( y ex1 , y ex2 , ..., y ex N | x ) ∝ N (cid:89) n =1 P ( y ex n | x ) , (10)where the index n denotes the type of physical quantities,and N is the total number of types. Thus, the posteriordistribution is written as P ( x | y ex1 , y ex2 , ..., y ex N ) ∝ N (cid:89) n =1 (cid:18) πσ n (cid:19) Ln exp [ − E ( x , σ , ..., σ N )] , (11)where L n is the number of data points for the n -th typemeasurement and the energy function is a weighted sumgiven as E ( x , σ , ..., σ N )= N (cid:88) n =1 (cid:34) σ n L n (cid:88) l =1 (cid:0) y ex n ( g l ) − y cal n ( g l , x ) (cid:1) (cid:35) − log P ( x ) . (12)This means that the posterior distribution significantlydepends on both L n and σ n .In this paper, for simplicity, the number of data pointsfor all inputted physical quantities is arranged as L n = L for ∀ n . Furthermore, the case where σ n does not dependon the type of physical quantities is considered. Thatis, the standard deviation of the observation noise is thesame for all types of physical quantities: σ n = σ for ∀ n .To make this assumption more realistic, the contributionsof each type of physical quantity are arranged and thefollowing normalization is imposed: y ex n ( g l ) → y ex n ( g l ) − min l ( y ex n ( g l ))max l ( y ex n ( g l )) − min l ( y ex n ( g l )) , (13) y cal n ( g l ) → y cal n ( g l ) − min l ( y ex n ( g l ))max l ( y ex n ( g l )) − min l ( y ex n ( g l )) . (14)By using the normalization, the upper and lower boundsfor each type of physical quantity are 1 and 0, respec-tively. IV. SPIN HAMILTONIAN ESTIMATIONA. Target Hamiltonian
Our model estimation method assumes the shape ofthe target Hamiltonian to be estimated. KCu P O hasa complicated three-dimensional structures (Fig. 1). Onthe other hand, by only focussing on the superexchangeinteractions, that is, the Cu–O–Cu paths in the crystalstructure, the lattice structure of the Cu ions is approxi-mated by a set of independent zigzag chains. Each chainhas eight Cu ions, and the bond lengths between the Cuions are smaller than 3.2 ˚A. Note that this approximationshould be valid except for extremely low temperatures.For low temperatures, weak interactions besides these su-perexchange interactions may affect magnetism ,and this approximation is poor. In this paper, by fo-cussing on the magnetic properties except for those at lowtemperatures, the target Hamiltonian is set to a Heisen-berg model on the chain with eight spins. Figure 3 showsthe lattice structure. Cu1Cu2 Cu3 Cu4 Cu4 Cu3 Cu2Cu1
FIG. 3. (Color online) Lattice structure of the target Hamil-tonian for KCu P O . This zigzag chain is constructed byeight Cu ions. Cu ions with the same index have equivalentpositions when considering symmetry. Although symmetry considerations imply that thereare seven types of nearest neighbor interactions in eachchain, only four types of independent interactions ap-pear (i.e., J , J , J , and J in Fig. 3). Furthermore,since the magnetic ion is Cu, anisotropy should not beconsidered. Consequently, the target Hamiltonian withisotropic Heisenberg spins is written as H ( x ) = − (cid:88) i =1 J i,i +1 (ˆ s xi ˆ s xi +1 + ˆ s yi ˆ s yi +1 + ˆ s zi ˆ s zi +1 ) , (15)where J = J , = J , , J = J , = J , , J = J , = J , , and J = J , . Here, (ˆ s xi , ˆ s yi , ˆ s zi ) are the S = 1 / i th site. The magnetic properties ofthis target quantum Hamiltonian can be easily calcu-lated by the exact diagonalization method. If an effec-tive model with a large number of spins need to be es-timated, another simulation method such as the MonteCarlo method or mean-field calculation should be per-formed.Since many magnetic interactions are already trimmed,there is no need to use the L regularization for the priordistribution. On the other hand, large magnetic interac-tions are not preferable for the physical sense. To avoidan increase in the absolute values of the estimated mag-netic interactions, we adopt the L regularization as theprior distribution. Furthermore, to estimate an effectivemodel to roughly capture the magnetic properties underwide temperature and magnetic field ranges, six typesof physical quantities ( N = 6) measured in KCu P O (see Fig. 2) are used as inputted data. In addition, thenumber of data points in each inputted physical quantityis fixed as L = 100. B. Determination of hyperparameter in the L regularization To determine the hyperparameter λ in the L regu-larization, the maximizer of the posterior distribution isanalyzed. Since σ does not generally depend on the valueof the model parameters from the viewpoint of a MAPestimation, we search for the minimizer of the followingequations: E (cid:48) ( x , α ) = ∆ (cid:48) ( x ) + α (cid:107) x (cid:107) , (16)∆ (cid:48) ( x ) = N (cid:88) n =1 L (cid:88) l =1 (cid:0) y ex n ( g l ) − y cal n ( g l , x ) (cid:1) , (17)where α := 2 σ λ , and E (cid:48) ( x , α ) and ∆ (cid:48) ( x ) are the energyfunction and error function normalized by 1 / σ . Notethat determining the plausible value of α is equivalent todeciding the hyperparameter λ in the L regularizationdepending on σ .The minimizer of E (cid:48) ( x , α ) is searched by the MCMCmethod where the probability distribution is proportionalto exp[ − E (cid:48) ( x , α )]. The MCMC samplings are performedby emcee package for various α . In each MCMC sam-pling, the stretch move is used for the update schemeof states, and 8,000 states are sampled. Among the sam-pled states with fixed α , the model parameters x ∗ suchthat E (cid:48) ( x , α ) is minimized are selected. As mentioned inSec. III A, this optimization of E (cid:48) ( x , α ) can be also per-formed by various fast optimization techniques instead ofMCMC. Figure 4 shows the α dependence of E (cid:48) ( x ∗ , α ),∆ (cid:48) ( x ∗ ), and (cid:107) x ∗ (cid:107) . Here, all experimental results, thatis, the magnetic susceptibility with 0.01 T and magneti-zation curves at five temperatures (1.3 K, 4.2 K, 20 K, 30K, and 50 K), for KCu P O are inputted. The resultof E (cid:48) ( x ∗ , α ) shows an elbow curve, indicating a bound-ary between regions, where in each part, the regulariza-tion is effective or ineffective. Thus, the appropriate α is selected as the elbow position. That is, α ∗ = 10 − ,because we want to find the spin Hamiltonian not onlyto explain the experimental results but also to obtainsmall magnetic interactions as possible. This is similarto the determination technique of the cluster number inclustering analysis .In the error function ∆ (cid:48) ( x ∗ ), the elbow curve is alsoobserved, but the elbow position deviates in E (cid:48) ( x ∗ , α ).Preferably, ∆ (cid:48) ( x ∗ ) is a small value at α ∗ = 10 − , and ifthe value of α decreases, the error is almost unchanged.In addition, (cid:107) x ∗ (cid:107) monotonically decreases against α ,and it becomes sufficiently small at α ∗ = 10 − . In thisstage, the estimated magnetic interactions (i.e., x ∗ ) at α ∗ = 10 − are J = − .
65 meV, J = − .
49 meV, J = − .
96 meV, and J = 6 .
92 meV. The abso-lute values of these interactions are proper in the Cusystems . These facts indicate that the L regulariza-tion is useful to estimate a spin Hamiltonian with smallmagnetic interactions. C. Evaluation of the uncertainty and distributionof sampling data
To evaluate the uncertainty of the estimated magneticinteractions, the noise amplitude is assessed according toSec. III B. Here, the Bayes free-energy in the estimation (c)(b) -1 -2 -3 -4 -5 -6 -7 -1 -1 appropriatehyperparameter i n e ff e c t i v e r e g u l a r i z a t i o n e ff e c t i v e r e g u l a r i z a t i o n (a) FIG. 4. (Color online) (a) Hyperparameter α dependenceof E (cid:48) ( x ∗ , α ) where x ∗ is the magnetic interactions so that E (cid:48) ( x , α ) is minimized with a fixed α . Dotted lines are a visualguide. (b) α dependence of ∆ (cid:48) ( x ∗ ). (c) α dependence of (cid:107) x ∗ (cid:107) . problem for KCu P O is given by F ( σ ) = (cid:90) σ dβ (cid:104) E (cid:48) ( x , α ∗ ) (cid:105) β + N L (cid:0) πσ (cid:1) , (18)where N = 6 and L = 100. Figure 5 (a) shows the inversetemperature β dependence of (cid:104) E (cid:48) ( x , α ∗ ) (cid:105) β by MCMC cal-culations using emcee package where the Monte Carlostep is 8,000. The average values of eight independentruns are plotted, and the error bars denote the 95% con-fident intervals. It monotonically decreases as a func-tion against β . Using the results of (cid:104) E (cid:48) ( x , α ∗ ) (cid:105) β , Fig. 5(b) shows the σ dependence of F ( σ ) which is calculatedby numerical integration with the trapezoidal rule. Theminimum value of F ( σ ) is obtained at 2 σ ∗ = 10 − . Thisvalue is a plausible standard deviation for the observa-tion noise. This means that the noise amplitude is ∼ − E (cid:48) ( x , α ∗ ) / σ ∗ ) with 2 σ ∗ = 10 − and α ∗ = 10 − . Figure 6 shows the scatterplots of thesampling results between all combinations of J , J , J , J , and E ( x ). Here, the last 3,000 sampling results in theMCMC with 8,000 Mote Carlo steps are plotted. First,we determine the estimated magnetic interactions andthe error bars as 95 % confidence interval are evaluated
01 10 -1 -2 -3 -4 -5 -1 -2 -3 -4 -5 (a)(b)
23 x10 plausiblevalue -1000 0 1000 2000 -2 -1 FIG. 5. (Color online) (a) Inverse temperature β dependenceof (cid:104) E (cid:48) ( x , α ∗ ) (cid:105) β by MCMC. (b) Bayes free-energy dependenceon the noise amplitude σ . Inset is the enlarged view. by these distributions as J = − . ± .
51 meV , (19) J = − . ± .
13 meV , (20) J = − . ± .
15 meV , (21) J = 6 . ± .
95 meV . (22)Note that the inputted experimental results (Fig. 2) arenot noisy, and the evaluated error bars for the estimatedinteractions are sufficiently small. Thus, to consider morenoisy cases, the artificial noises are added to the exper-imental results of KCu P O . Figures S1 and S2 insupplemental material show the estimation results de-pending on the artificial noise. We confirm that if theartificial noise is increased, the value of σ ∗ by our noiseestimation is also increased and consequently the errorbars for the estimated interactions become large. Thus,we conclude that our estimation method can be appliedto the noisy experimental results and correctly evaluatethe uncertainty in the estimated effective model.Next, from the distributions of magnetic interactions(Fig. 6), the region of J realizing good fitting is nar-rower compared to the other interactions. That is, J is the most sensitive. In contrast, changing J has asmaller impact on the fitting error, indicating J has thehighest uncertainty. Furthermore, we can see that J and J are positively correlated, while J negatively iscorrelated with J and J . If J increases, J and J should decrease to maintain the good fitting. On theother hand, J is almost independent of the magnetic in-teractions. Consequently, J can be freely tuned withinthe error bar. We discuss these extracted correlationsfrom physical insights. Considering the lattice symme- FIG. 6. (Color online) Scatterplots of the sampling resultsbetween all combinations of J , J , J , J , and E (cid:48) ( x , α ∗ ) byMCMC around the estimated magnetic interactions. Here,the probability distribution in MCMC is proportional toexp( − E (cid:48) ( x , α ∗ ) / σ ∗ ) with 2 σ ∗ = 10 − and α ∗ = 10 − . try, J and J are placed in a similar environment, forexample, there are two places in the lattice and the in-teracted spins are not edge one. Thus, J and J wouldhave a positive correlation. Furthermore, J is an anti-ferromagnetic interaction as well as J and J , and J should become smaller with increasing J and J to keepthe energy scale of antiferromagnetic interactions in theHamiltonian, which means that J is negatively corre-lated with J and J . We note that these correlationsobtained in the estimation strongly depend on the tar-get Hamiltonian and estimated interaction parameters.Thus, these discussions are correlations of the estimatedinteraction parameters in the effective model, not cor-relations of the magnetic interactions in the real mate-rials. However, these extracted correlations can deepenan understanding of properties of the estimated effectivemodel.Figure 7 (a) compares the physical quantities betweenthe experimental and calculated results by the estimatedspin Hamiltonian with J = − .
54 meV, J = − . J = − .
90 meV, and J = 6 .
24 meV. A goodfit can be obtained except for the low temperature mag-netization curves. In particular, for 1.3 K, the fittingresult is worse and the 1/4 magnetic plateau-like behav-ior is observed in the calculation results. Since our tar-get Hamiltonian only has eight spins, the appearance ofsuch plateau cannot be avoided at low temperatures. Todescribe the low temperature magnetizations, the target estimationexperiment error ofsusceptibility error of 1.3 Kmagnetizationerror of 4.2 Kmagnetizationerror of 20 Kmagnetizationerror of 30 Kmagnetizationerror of 50 Kmagnetization e rr o r o f s u sc ep t i b ili t y e rr o r o f K m agne t i z a t i one rr o r o f K m agne t i z a t i one rr o r o f K m agne t i z a t i one rr o r o f . K m agne t i z a t i one rr o r o f . K m agne t i z a t i on (b)(a) FIG. 7. (Color online) (a) Comparison plots of the magnetic susceptibility with 0.01 T and the magnetization curves withvarious temperatures between the experimental and calculated results by the estimated spin Hamiltonian. (b) Scatterplots ofthe errors between the experimental and calculated results, which are independently evaluated in the susceptibility and fivetypes of magnetization curves.
Hamiltonian must include further magnetic interactionsbesides the four interactions.Moreover, Fig. 7 (b) shows the scatterplots of the sam-pling results by the same MCMC in Fig. 6 for all com-binations of errors between the experimental and calcu-lated results, i.e., (cid:80) Ll =1 (cid:0) y ex n ( g l ) − y cal n ( g l , x ) (cid:1) . The er-rors of susceptibility and magnetization curves with 50K, 30 K, and 20 K have strong correlations. For exam-ple, susceptibility has a positive correlation with mag-netization below 50 K, but it has a negative correlationwith magnetizations under 30 K and 20 K. This meansthat if the susceptibility is fitted, the magnetization un-der 50 K is naturally fitted, while the fittings of mag-netizations with 30 K and 20 K become worse. On theother hand, magnetization curves with 4.2 K and 1.3 Kare almost independent of the other errors around theestimated magnetic interactions. In this way, by drawingthe distributions of the magnetic interactions or physicalquantities by MCMC, their relations can be understoodand the characteristics of the estimated spin Hamiltonianare extracted. D. Prediction of the magnetic properties
The most important benefit from estimating the effec-tive model is predicting magnetic properties that cannotbe easily or have not been measured. Thus, the valueof the spin gap, which is the energy gap between theground and excited states, the spin configuration at theground state without a magnetic field, and the temper-ature dependences of magnetic specific heat and mag-netic entropy, are calculated (Fig. 8). The predictions indicate that the magnetic specific heat in KCu P O will have a peak around 30 K, but increasing the mag-netic field suppresses this peak (Fig. 8 (b)). In this mag-net, the magnetic entropy will be almost unchanged bya magnetic field over 20 K (Fig. 8 (c)). Furthermore,to predict the magnetic refrigeration property , thechange in magnetic entropy is also evaluated (Fig. 8(d)). The inverse magnetocaloric effect will be observedin KCu P O , which is a characteristic behavior ofantiferromagnets , and the magnetic entropy changeshould be quite small. In this way, using an effectivemodel determined by our data-driven approach, difficult-to-measure properties can be predicted, improving theunderstanding of the magnetic properties of KCu P O . V. DISCUSSION AND SUMMARY
We have determined the spin Hamiltonian ofKCu P O using a data-driven approach. The flow ofour prescription of data-driven approach is as follows. (i)Assume a target effective model and the posterior dis-tribution. This constitutes the difference between theexperimental and calculated results obtained by an ef-fective model and the appropriate prior distribution ofmodel parameters. (ii) Determine an appropriate hyper-parameter in the prior distribution by the elbow methodfor the MAP estimation results and estimated model pa-rameters. (iii) Obtain a plausible observation noise tominimize the Bayes free-energy. (iv) Perform MCMCsamplings using an estimated noise amplitude aroundthe estimated model parameters in (ii). From the ob-tained sampling data distributions, determine the model Spin gap: (d)(c) (b)(a)
Cu1Cu2 Cu3 Cu4 Cu4 Cu3 Cu2Cu1
FIG. 8. (Color online) (a) Spin gap value and spin configu-ration at the ground state without a magnetic field obtainedby the estimated spin Hamiltonian for KCu P O . Tem-perature dependences of (b) the magnetic specific heat, (c)magnetic entropy, and (d) magnetic entropy change when themagnetic field changes from H to 0 T are plotted. parameters with uncertainty. (v) Predict various prop-erties, which cannot be easily measured in experimentsusing the estimated effective model.Our data-driven approach found that the spin Hamil-tonian of KCu P O is the quantum Heisenberg modelon the zigzag chain with eight spins of J = − . ± .
51 meV, J = − . ± .
13 meV, J = − . ± .
15 meV, and J = 6 . ± .
95 meV. In this estimation,the magnetic susceptibility and magnetization curves atvarious temperatures are fitted. This model can describethe experimental results with very small error, exceptfor the extremely low temperature magnetization curves,demonstrating that such a data-driven approach is usefulto estimate the effective model. Our approach should beuseful in parallel with ab initio calculations.On the other hand, the results by the data-driven ap-proach should strongly depend on the inputted data. Ac-tually, we estimated different magnetic interactions whenthe inputted data is only one type of physical quantity(See Figs. S3 – S8 in supplemental material.). In thiscase, the fitting of the inputted quantities is well per-formed, but the experimental results, which are not usedin the estimation are not well fitted. This fact means thateach physical quantity can be well described by multiplesets of magnetic interactions. Consequently, preparingmultiple kinds of physical quantities as the input datawill not only estimate model parameters more uniquelybut will also produce a more reliable effective model.Recently, the parameter estimation for real materials by data-driven techniques has been frequently performed.These data driven techniques are roughly divided intotwo strategies. (i) A regression model that explains ma-terial parameters from feature variables such as materialcomposition and structure is constructed from a largeamount of collected data for known materials . Theaim of this strategy is to predict parameters in unknownmaterials using the trained regression model via data-driven approach. (ii) The material parameters are deter-mined by minimizing the discrepancy between physicalquantities of the target material and those obtained by acomputational model with the material parameters .The aim of this strategy is to understand properties ofthe target material through the estimated parameters inthe computational model. The former, however, requiresthe preparation of a large amount of data sets of materialcompositions and structures paired with target parame-ters, i.e., magnetic interactions in our problem. There-fore, for the purpose of an estimation of magnetic inter-actions, the former is not necessarily suitable, but thelatter is more suitable, and indeed our method belongsto the latter strategy. Furthermore, our proposed noiseestimation method using MCMC to evaluate the uncer-tainty of estimated parameters can be applied to variousmethods belonging to the latter strategy, and we believethat it is also useful in various data-driven techniquesto estimate materials parameters in the computationalmodel.Our data-driven approach will come into its own wheneasy-to-measure properties obtained in the laboratorysuch as SQUID are inputted. That is, our data-driven ap-proach can predict difficult-to-measure properties, whichare obtained by large-scale experimental equipment suchas neutron scattering. Thereby, our data-driven ap-proach will reduce the cost of materials development andaccelerate discoveries of novel materials.
ACKNOWLEDGEMENT
We thank Hideaki Kitazawa and Noriki Terada forthe valuable discussions and Takao Furubayashi for ESRmeasurements. This article is supported by the “Ma-terials Research by Information Integration” Initiative(MI2I) project and JST-Mirai Program Grant No. JP-MJMI18A3. M.H. was partially supported by a grant foradvanced measurement and characterization technologiesaccelerating the materials innovation at the National In-stitute for Materials Science (NIMS). The computationswere performed on Numerical Materials Simulator atNIMS and the supercomputer at Supercomputer Cen-ter, Institute for Solid State Physics, The University ofTokyo. ∗ [email protected] † [email protected] ‡ [email protected] D. Mu˜noz, I. de P. R. Moreira, and F. Illas, Phys. Rev. B , 224521 (2002). I. I. Mazin, Phys. Rev. B , 140406(R) (2007). K. Nakamura, Y. Yoshimoto, R. Arita, S. Tsuneyuki, andM. Imada, Phys. Rev. B , 195126 (2008). J. M. Pruneda, J. ´I˜niguez, E. Canadell, H. Kageyama, andM. Takano, Phys. Rev. B , 115101 (2008). M. Hirayama, T. Miyake, and M. Imada, Phys. Rev. B , 195144 (2013). M. Nuss and M. Aichhorn, Phys. Rev. B , 045125 (2014). T. Misawa and M. Imada, Nature Commun. , 5738 (2014). K. Riedl, D. Guterding, H. O. Jeschke, M. J. P. Gingras,and R. Valent´ı, Phys. Rev. B , 014410 (2016). M. Hirayama, Y. Yamaji, T. Misawa, and M. Imada, Phys.Rev. B , 134501 (2018). K. Takubo, T. Mizokawa, J.-Y. Son, Y. Nambu, S. Nakat-suji, and Y. Maeno, Phys. Rev. Lett. , 037203 (2007). R. Tamura, N. Kawashima, T. Yamamoto, C. Tassel, andH. Kageyama, Phys. Rev. B , 214408 (2011). R. A. DiStasio, J. ´E. Marcotte, R. Car, F. H. Stillinger,and S. Torquato, Phys. Rev. B , 134104 (2013). M. Hase, H. Kuroe, V. Y. Pomjakushin, L. Keller,R. Tamura, N. Terada, Y. Matsushita, A. D¨onni, andT. Sekine, Phys. Rev. B , 054425 (2015). M. Hase, M. Matsumoto, A. Matsuo, and K. Kindo, Phys.Rev. B , 174421 (2016). M. Hase, Y. Ebukuro, H. Kuroe, M. Matsumoto, A. Mat-suo, K. Kindo, J. R. Hester, T. J. Sato, and H. Yamazaki,Phys. Rev. B , 144429 (2017). K. Rajan, materialstoday , 38 (2005). A. R. Oganov and C. W. Glass, J. Chem. Phys. ,244704 (2006). G. Pilania, C. Wang, X. Jiang, S. Rajasekaran, andR. Ramprasad, Sci. Rep. , 2810 (2013). A. Seko, A. Togo, H. Hayashi, K. Tsuda, L. Chaput, andI. Tanaka, Phys. Rev. Lett. , 205901 (2015). T. Ueno, T. D. Rhone, Z. Hou, T. Mizoguchi, andK. Tsuda, Materials Discovery , 18 (2016). S. Ju, T. Shiga, L. Feng, Z. Hou, K. Tsuda, and J. Shiomi,Phys. Rev. X , 021024 (2017). R. Ramprasad, R. Batra, G. Pilania, A. Mannodi-Kanakkithodi, and C. Kim, npj Comput. Mater. , 54(2017). H. Ikebata, K. Hongo, T. Isomura, R. Maezono, andR. Yoshida, J. Comput.-Aided Mol. Des. , 379 (2017). M. Sumita, X. Yang, S. Ishihara, R. Tamura, andK. Tsuda, ACS Cent. Sci. , 1126 (2018). T. Yamashita, N. Sato, H. Kino, T. Miyake, K. Tsuda,and T. Oguchi, Phys. Rev. Materials , 013803 (2018). R. G´omez-Bombarelli, J. N. Wei, D. Duvenaud, J. M.Hern´andez-Lobato, B. S´anchez-Lengeling, D. Sheberla,J. Aguilera-Iparraguirre, T. D. Hirzel, R. P. Adams, andA. Aspuru-Guzik, ACS Cent. Sci. , 268 (2018). Z. Hou, Y. Takagiwa, Y. Shinohara, Y. Xu, and K. Tsuda,ACS Appl. Mater. , 11545 (2019). K. Terayama, R. Tamura, Y. Nose, H. Hiramatsu,H. Hosono, Y. Okuno, and K. Tsuda, Physical ReviewMaterials , 033802 (2019). J. Behler and M. Parrinello, Phys. Rev. Lett. , 146401(2007). L. M. Ghiringhelli, J. Vybiral, S. V. Levchenko, C. Draxl,and M. Scheffler, Phys. Rev. Lett. , 105503 (2015). T. L. Pham, H. Kino, K. Terakura, T. Miyake, and H. C.Dam, J. Chem. Phys. , 154103 (2016). R. Kobayashi, D. Giofr´e, T. Junge, M. Ceriotti, and W. A.Curtin, Phys. Rev. Materials , 053604 (2017). T. Suzuki, R. Tamura, and T. Miyazaki, Int. J. Quant.Chem. , 33 (2017). H. C. Dam, V. C. Nguyen, T. L. Pham, A. T. Nguyen,K. Terakura, T. Miyake, and H. Kino, J. Phys. Soc. Jpn. , 113801 (2018). K. Shiba, R. Tamura, T. Sugiyama, Y. Kameyama,K. Koda, E. Sakon, K. Minami, H. T. Ngo, G. Imamura,K. Tsuda, and G. Yoshikawa, ACS Sensors , 1592 (2018). R. Tamura, J. Lin, and T. Miyazaki, J. Phys. Soc. Jpn. , 044601 (2019). R. Tamura and K. Hukushima, PLoS ONE , e0193785(2018). H. Takenaka, K. Nagata, T. Mizokawa, and M. Okada, J.Phys. Soc. Jpn. , 124706 (2014). R. Tamura and K. Hukushima, Phys. Rev. B , 064407(2017). H. Fujita, Y. O. Nakagawa, S. Sugiura, and M. Oshikawa,Phys. Rev. B , 075114 (2018). K. Momma and F. Izumi, J. Appl. Cryst. , 1272 (2011). H. Effenberger, Zeitschrift f¨ur Kristallographie , 43(1987). C. M. Bishop,
Pattern Recognition and Machine Learning (Springer, 2006). M. Anada, Y. Nakanishi-Ohno, M. Okada, T. Kimura, andY. Wakabayashi, J. Appl. Cryst. , 1611 (2017). S. Tokuda, K. Nagata, and M. Okada, J. Phys. Soc. Jpn. , 024001 (2017). K. Obinata, S. Katakami, Y. Yue, and M. Okada, J. Phys.Soc. Jpn. , 064802 (2019). S. Katakami, H. Sakamoto, and M. Okada, J. Phys. Soc.Jpn. , 074001 (2019). M. Hase, Y. Ebukuro, H. Kuroe, M. Matsumoto, A. Mat-suo, K. Kindo, J. R. Hester, T. J. Sato, and H. Yamazaki,Phys. Rev. B , 139901 (2018). http://dfm.io/emcee, (Accessed on 1st December 2018). D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Good-man, arXiv:1202.3665 (2012). J. Goodman and J. Weare, Communications in AppliedMathematics and Computational Science , 65 (2010). R. L. Thorndike, Psychometrika , 267 (1953). R. Tibshirani, G. Walther, and T. Hastie, J. R. Statist.Soc. B , 411 (2001). K. A. Gschneidner and V. K. Pecharsky, Annu. Rev.Mater. Sci. , 387 (2000). A. M. Tishin and Y. I. Spichkin,
The Magnetocaloric Effectand its Applications (Taylor & Francis, London, 2003). K. A. Gschneidner, V. K. Pecharsky, and A. O. Tsokol,Rep. Prog. Phys. , 1479 (2005). K. Matsumoto, T. Kondo, M. Ikeda, and T. Numazawa,Cryogenics , 353 (2011). V. Franco, J. S. Bl´azquez, B. Ingale, and A. Conde, Annu.Rev. Mater. Res. , 305 (2012). V. Franco, J. S. Bl´azquez, J. J. Ipus, J. Y. Law, L. M.Moreno-Ram´ırez, and A. Conde, Prog. Mater. Sci. ,112 (2018). K. G. Sandeman, R. Daou, S. ¨Ozcan, J. H. Durrell, N. D.Mathur, and D. J. Fray, Phys. Rev. B , 224436 (2006). T. Krenke, E. Duman, M. Acet, E. F. Wassermann,X. Moya, L. Ma˜nosa, A. Planes, E. Suard, and B. Oulad- diaf, Phys. Rev. B , 104414 (2007). R. Tamura, T. Ohno, and H. Kitazawa, Appl. Phys. Lett. , 052415 (2014). R. Tamura, S. Tanaka, T. Ohno, and H. Kitazawa, J.Appl. Phys. , 053908 (2014). L. Ward, A. Agrawal, A. Choudhary, and C. Wolverton,npj Comput. Mater. , 16028 (2016). V. Stanev, C. Oses, A. G. Kusne, E. Rodriguez,J. Paglione, S. Curtarolo, and I. Takeuchi, npj Comput.Mater. Materials , 29 (2018). M. Franulovi´c, R. Basan, and I. Prebil, ComputationalMaterials Science , 505 (2009). K. M. El-Naggar, M. R. AlRashidi, M. F. AlHajri, andA. K. Al-Othman, Solar Energy , 266 (2012). Y. Tsukada, S. Takeno, M. Karasuyama, H. Fukuoka,M. Shiga, and T. Koyama, Sci. Rep.9