A new mathematical approach to finding global solutions of the magnetic structure determination problem.

Determination of magnetic structure is an important analytical procedure utilized in various fields ranging from fundamental condensed-matter physics and chemistry to advanced manufacturing. It is typically performed using a neutron diffraction technique; however, finding global solutions of the magnetic structure optimization problem represents a significant challenge. Generally, it is not possible to mathematically prove that the obtained magnetic structure is a truly global solution and that no solution exists when no acceptable structure is found. In this study, the global optimization technique called semidefinite relaxation of quadratic optimization, which has attracted much interest in the field of applied mathematics, is proposed to use as a new analytical method for the determination of magnetic structure, followed by the application of polarized neutron diffraction data. This mathematical approach allows avoiding spurious local solutions, decreasing the amount of time required to find a tentative solution and finding multiple solutions when they exist.

Introduction

Development of method for determining crystal and magnetic structure is an important issue[1-3]. The magnetic structure of a compound is directly related to the microscopic origins of various intriguing magnetic phenomena observed in physics, chemistry, biology and geology[4-11]. Accurate information on the experimentally determined magnetic structures can be potentially used to effectively design functional materials and opens multiple opportunities for utilizing advanced first-principle calculations and informatics approaches. Furthermore, the responses of the analysed magnetic structure to external factors such as temperature, pressure, electromagnetic field, light and environment are usually related to the corresponding changes in the macroscopic properties of minerals and biomolecules.

Since the pioneering work of Shull et al. published in 1951[1], many studies on the determination of the magnetic structures of various materials by neutron diffraction have been performed. Thus, Izymov, Kovalev and Bertaut developed the irreducible representation (IR) theory, which successfully classified magnetic structures by their magnetic point symmetries obtained through the analysis of neutron diffraction data[12-14]. Moreover, the modern neutron diffractometers exhibit relatively high statistical accuracy and resolution[15,16]. However, analysing neutron diffraction data to determine the correct magnetic structure represents a challenging task. The major issue here is the existence of so-called local solutions to the optimization problem (Fig. 1A). The analysis procedure is mathematically categorized as the nonlinear optimization, in which the obtained diffraction data points are fitted with the function of neutron scattering cross-section characterized by high-dimensional nonlinear parameters such as magnetic moment vectors[17]. After performing nonlinear optimization, it is very difficult to prove that the best solution obtained via various numerical models is the global one. Furthermore, multiple or many solutions are hidden in some cases, and no solutions exist in other cases. To mitigate these issues, global optimization techniques including a simulated annealing method, the Monte Carlo method and genetic algorithms can be used[18-23]. However, none of these approaches can determine the global solution(s) with 100% probability.

Figure 1

Schematic diagrams illustrating the non-linear optimization process. (A) Before relaxation. (B) SDR technique. (C) Duality theorem combined with SDR. The objective function f is globally minimized to the global optimum value . When the duality gap is close to zero (corresponding to the square root of the machine epsilon), the global optimization process is complete according to the duality theorem. In panels (B) and (C), the solid balls represent the paths toward feasible solutions chosen by the interior point method.

In this study, a new analytical method for the determination of magnetic structure is proposed, which allows (i) judging whether the obtained solution is a truly global one; (ii) evaluating the uniqueness of the solution obtained for a given set of experimental data; and (iii) completing the global optimization procedure in a very short time (typically, less than several seconds). Furthermore, the applicability of the proposed method is demonstrated using polarized-neutron powder diffraction data.

Mathematical Formulation for Applications of Semidefinite Relaxation Method

Relaxation techniques have been previously developed in the field of mathematical programming. In particular, semidefinite relaxation (SDR) combined with semidefinite programming (SDP) has found many applications in applied mathematics and engineering (for example, see Chapter 2.2 in ref.[24]). SDR is an efficient technique for solving nonlinear optimization problems, such as quadratic programs (QPs), namely the minimization of multivariate quadratic polynomials under constraints. It should be noted that SDR could provide both fast convergence and a numerical proof of the global optimality of the obtained solutions using the global convergence property and duality theorem of the convex programming methods, as schematically shown in Fig. 1(B,C).

In the field of optical imaging, Candès et al. published a pioneering work (called the PhaseLift method), in which a sparse modelling approach was adopted for general phase retrieval[25]. Furthermore, one of the authors (ROT) developed an SDR-based mathematical approach to investigate whether a crystal structure could be uniquely identified using only diffraction data and independent atomic model[26].

In this study, this method was applied to magnetic structure analysis by replacing the optimized parameters with magnetic moments. The problem of magnetic structure determination can be mathematically described by the following set of N quadratic equations:where and denote the observed and calculated integrated magnetic diffraction intensities at the magnetic reflection , respectively;  = (x1, …, x)T and represents the coefficient matrix numerically obtained using the scattering cross-section formula for each , which consists of the absolute intensity scale factor, magnetic form factor, Lorentz factor, multiplicity and structural weight factor (the details of this formula are summarized in Supplementary Information).

The symbol indicates that each value includes an experimental error. In order to incorporate all errors, the following optimization problem is solved:

The inequality can be also removed if the permitted range of (, for example, the upper limit of magnetic moment) is uncertain while maintaining the validity of the subsequent discussion of the SDR and SDP techniques.

Equation (2) is classified as a so-called l1-norm minimization problem, which is equivalent to the following form of the quadratic programming problem:

The problems described by equations (2) and (3) are equivalent since they have the same set of solutions and the minimized values of the objective functions (more details are provided in Supplementary Information). The basic strategy utilized in this work is to apply SDR to equation (3) in order to solve equation (2) and then provide a computational proof on the global convergence property of the solution. The least-squares minimization of the function is avoided because a small increase in the polynomial degree considerably magnifies the size of the SDR problem (it increases proportionally to the power of N). Thus, it is possible to locate the global optimum by performing l1-norm minimization, which can be subsequently used as the initial parameter of the normal least-squares method to calculate the refined parameters and then compare them with literature data.

Application to Experimental Data Analysis

Experimental results

The developed SDR method was verified by applying it to the experimental data obtained for pyrochlore Nd2Ir2O7, which served as a proximate material for a three-dimensional Weyl semimetal and a component of spintronic devices on the basis of its geometrically frustrated magnetism and 5d-electron configuration[27-32]. This state theoretically corresponds to the all-in all-out type of magnetic structure described by the magnetic propagation vector mag = (0, 0, 0), in which all the magnetic moments are oriented either towards the centre of the participating tetrahedron or in the opposite direction[30,32]. This prediction was experimentally confirmed in previous unpolarized neutron diffraction studies[33,34]. However, the superposition of the nuclear and magnetic reflections observed for the magnetic structures with mag = (0, 0, 0) produces ambiguous results during their separation. Furthermore, it is not possible to mathematically prove with 100% certainty the absence of other acceptable solutions to the problem of magnetic structure determination. Thus, in this work, polarized neutron diffraction studies were performed for this material, and the obtained magnetic structure was verified mathematically.

Figure 2 shows the representative neutron diffraction data obtained for Nd2Ir2O7 at the minimum temperature T = 1.4 K. The non-spin-flip and spin-flip parameters IOFF and ION roughly correspond to the nuclear and magnetic reflection components Inuc and Imag, respectively. The magnitude of IOFF is systematically larger than ION at all temperatures, thus confirming the necessity of conducting polarized neutron diffraction experiments. For ION, the intensity of the 113 reflection increases with decreasing temperature, indicating the existence of strong temperature dependence for this reflection, which is not observed for the 222 reflection. Moreover, as the flipping ratio of the neutron beam in the actual experiments is finite, the magnitudes of ION and IOFF can be expressed by the following formulas:where rmix = N−/(N+ + N−) is the mixing rate, and N+ (N−) denotes the number of majority-spin (minority-spin) neutrons. The value of rmix is selected to satisfy the condition Imag = 0; hence, ION = rmix·Inuc at the paramagnetic T = 40 K. As a result, the pure magnetic intensity Imag can be determined by combining IOFF and ION. A more detailed explanation and data obtained for other reflections are presented in Supplementary Information, while the calculated Imag,obs values are listed in Table 1.

Figure 2

Typical polarized neutron diffraction data obtained for pyrochlore Nd2Ir2O7. ON/OFF denotes the spin-flip/non-spin-flip scattering channel. The lower panel displays both the overall OFF and ON data, whereas the upper panel depicts the low-intensity ON data at a higher magnification. The solid and dotted curves show the results of Gaussian fitting.

Table 1

Experimentally observed and calculated integrated magnetic scattering intensities. The symbol ± denotes the experimental errors.

h mag k mag l mag	Imag,obs (1.4 K)	I mag,cal,Γ3	I mag,cal,Γ5	I mag,cal,Γ7	I mag,cal,Γ9
111	0 ± 53	0	7	0	457
200	0 ± 36	0	0	0	1
220	932 ± 54	932	932	0	1
113	1796 ± 68	1795	4	0	860
222	0 ± 63	0	0	0	0
400	0 ± 37	0	0	0	0
331	349 ± 54	580	5	0	349
420	884 ± 62	588	541	0	2

Magnetic structure analysis results

The arbitrary magnetic structure with mag = (0, 0, 0) in the space group is described as the function of arb = (m(Nd)1, …, m(Nd)12, m(Ir)1, …, m(Ir)12) containing narb = 24 variables generated by the eight magnetic atoms of the unit cell and three Cartesian coordinates. Here, the magnetic symmetry was classified by performing IR analysis (this procedure is generally not required for the proposed method; however, it may be potentially useful because of the reduction of the number of variables if the number of observed reflection points is limited). The possible magnetic structures are described by the representations ΓNd = Γ3 + 2Γ5 + 3Γ7 + 6Γ9 and ΓIr = Γ3 + 2Γ5 + 3Γ7 + 6Γ9, where Γ3 corresponds to the all-in all-out type of the magnetic structure, and the coefficients denote the numbers of basis vectors summarized in Table 2[18,35]. Hence, the total numbers of variables are reduced to nΓ3 = 2, nΓ5 = 4, nΓ7 = 6 and nΓ9 = 12. Furthermore, representation Γ9 includes six ferromagnetic basis vectors. However, the Nd2Ir2O7 structure exhibits only extremely weak spontaneous magnetization (around 10−4μB/formula unit)[29], which is significantly below the detection limit of the neutron diffraction technique. Therefore, the ferromagnetic basis vectors are precluded and the value of nΓ9 is further reduced from 12 to 6 during the analysis of the neutron diffraction data.

Table 2

Basis vectors of the magnetic structures represented by the Γ3, Γ5, Γ7 and Γ9 IRs[18,35]. The symbol – in the Γ9 IR indicates that the ferromagnetic basis vectors are precluded.

IRs	Variables x	Nd1, Ir1	Nd2, Ir2	Nd3, Ir3	Nd4, Ir4
		(0,0,0), (1/2,0,0)	(1/4,1/4,0), (3/4,1/4,0)	(0,1/4,1/4), (1/2,1/4,1/4)	(1/4,0,1/4), (3/4,0,1/4)
Γ3	m(Γ3,Nd)1, m(Γ3,Ir)1	(1,1,1)/3	(−1,−1,1)/3	(1,−1,−1)/3	(−1,1,−1)/3
Γ5	m(Γ5,Nd)1, m(Γ5,Ir)1	(1,−1,0)/2	(−1,1,0)/2	(1,1,0)/2	(−1,−1,0)/2
	m(Γ5,Nd)2, m(Γ5,Ir)2	(1,1,−2)/6	(−1,−1,−2)/6	(1,−1,2)/6	(−1,1,2)/6
Γ7	m(Γ7,Nd)1, m(Γ7,Ir)1	(0,−1,1)/2	(0,−1,−1)/2	(0,1,−1)/2	(0,1,1)/2
	m(Γ7,Nd)2, m(Γ7,Ir)2	(1,0,−1)/2	(1,0,1)/2	(−1,0,−1)/2	(−1,0,1)/2
	m(Γ7,Nd)3, m(Γ7,Ir)3	(−1,1,0)/2	(1,−1,0)/2	(1,1,0)/2	(−1,−1,0)/2
Γ9	m(Γ9,Nd)1, m(Γ9,Ir)1	(1,1,0)/2	(−1,−1,0)/2	(−1,1,0)/2	(1,−1,0)/2
	m(Γ9,Nd)2, m(Γ9,Ir)2	(0,1,1)/2	(0,1,−1)/2	(0,−1,−1)/2	(0,−1,1)/2
	m(Γ9,Nd)3, m(Γ9,Ir)3	(1,0,1)/2	(1,0,−1)/2	(−1,0,1)/2	(−1,0,−1)/2
	—	(1,0,0)	(1,0,0)	(1,0,0)	(1,0,0)
	—	(0,1,0)	(0,1,0)	(0,1,0)	(0,1,0)
	—	(0,0,1)	(0,0,1)	(0,0,1)	(0,0,1)

Thus, the variables used for the magnetic structure analysis in this work are defined as follows: Γ3 = (m(Γ3,Nd)1, m(Γ3,Ir)1), Γ5 = (m(Γ5,Nd)1, m(Γ5,Nd)2, m(Γ5,Ir)1, m(Γ5,Ir)2), Γ7 = (m(Γ7,Nd)1, m(Γ7,Nd)2, m(Γ7,Nd)3, m(Γ7,Ir)1, m(Γ7,Ir)2, m(Γ7,Ir)3) and Γ9 = (m(Γ9,Nd)1, m(Γ9,Nd)2, m(Γ9,Nd)3, m(Γ9,Ir)1, m(Γ9,Ir)2, m(Γ9,Ir)3). The goal is to determine the globally optimal solutions for representations Γ3, Γ5, Γ7 and Γ9 using the Imag,obs magnitudes listed in Table 1. After that, the solution characterized by the best fit can be selected.

The output values produced by the SDP solver are as follows.

The convergence procedure results in the following coefficients:

andwhere the values of are expressed the μB units, and the negative sign indicates the nearest ferromagnetic Nd–Ir correlation. While a single optimum solution is determined for representation Γ3, multiple solutions are obtained for Γ5 and Γ9, respectively, in accordance with the rank values of XΓ5(opt) and XΓ9(opt), which indicate the existence of an infinite number of solutions with the same objective functions and , where and are arbitrary. Without SDR, it is difficult to find the optimum solutions for this type of problems and prove that better solutions do not exist. Furthermore, the zero magnitude of Γ7(opt) is obtained for representation Γ7. Indeed, when Γ7 values are finite (non-zero), the values for indices increase more rapidly than those for indices. Thus, the SDR method automatically overcomes these issues.

The obtained indicator values for representations Γ3, Γ5, Γ7 and Γ9 are as follows.

The corresponding duality gaps determined by the SDP solver are equal to , , and , respectively. Their magnitudes are close to zero, indicating that the convergence to the global optimums is achieved.

The minimized objective function  = 4.5, 19, 32 and 27, and the more familiar R-factor,  = 0.13, 0.63, 1.0 and 0.81, respectively. These results show that the smallest values are obtained for representation Γ3.

The goodness of the fit  = 2.6, 14, 25 and 20, and the related statistical parameter , 4.8, 16 and 14, respectively. The value closest to the expected value (1 for ; for the latter parameter) is obtained for Γ3.

All these facts suggest that Γ3 represents the best solution, whereas representations Γ5, Γ7 and Γ9 are excluded for the first time. For clarity, the values of Imag,cal are listed in Table 2 and also shown in Fig. 3. The magnitude of Imag,cal-Γ3 matches Imag,obs very well, whereas the globally optimal Imag,cal-Γ5, Imag,cal-Γ7 and Imag,cal-Γ9 parameters substantially differ from the corresponding Imag,obs values.

Figure 3

Comparison of the Imag,obs (1.4 K) and Imag,cal values. The diagonal straight line represents the condition Imag,obs = Imag,cal.

After proving mathematically the uniqueness of representation Γ3, the refined values mNd = −1.22(5)μB and mIr = 0.14(5)μB are obtained using the Γ3(opt) values as the initial parameters in the normal least-squares method. These results are comparable with the reported values of mNd = −1.27(1)μB and mIr = 0.34(1)μB[34].

Discussion

Representation Γ9 (the only IR containing ferromagnetic basis vectors[35]) is rejected as the optimal solution despite the existence of weak ferromagnetism[29]. The observed inconsistency suggests that this ferromagnetism is symmetrically decoupled from the bulk magnetic structure (for example, as a surface or interface-protected property). In fact, by performing careful measurements of the macroscopic magnetization and electrical resistivity of isomorphic Cd2Os2O7, it was found that the surface ferromagnetism coupled with novel spin-polarized conductivity emerged on the walls between the all-in all-out and all-out all-in antiferromagnetic domains[36]. The macroscopic observation is consistent with the findings of this study verified both microscopically and mathematically, indicating their potential applicability in domain wall spin electronics.

We discuss the expected application scope of the SDR method. After determining the list of (hmag, kmag, lmag, Imag,obs) values, the described SDR method can be used to find the corresponding magnetic structures as the global solutions. First, the configurations with mag = (0, 0, 0) are considered, indicating their high potential applicability in various fields (including magnet materials). Second, this method is not restricted to polarized neutron diffraction experiments, but can be also used in studies involving unpolarized neutrons when the magnetic structure is characterized by mag ≠ (0, 0, 0) or mag = (0, 0, 0) with detectable magnetic moment. Both the magnetic and crystallographic structures can be simultaneously refined by the normal least-squares method, in which the obtained SDR solutions are utilized as the initial values of the magnetic structural parameters. Third, the SDR technique is able to easily process thousands of independent variables. Therefore, magnetic structures of arbitrary types can be theoretically determined using the advanced diffractometers that provide a relatively large number of reflection points (even for the target materials with complex compositions).

The diffraction intensities are not represented by quadratic functions of the atomic positions ; hence, the SDR technique seems to be inapplicable for determining the values of (crystal structure). However, the neutron and X-ray diffraction intensities are represented by those of the nuclear and electron densities (generalized crystal structure), respectively. Likewise, a magnetic structure is also generalized to the magnetic moment density. Thus, the SDR technique is expected to enable analysing the global solutions of the generalized structures together with the aforementioned high-volume processing ability, such as protonic/ionic distributions in the conductors and electronic spin-orbital distributions.

Materials and Methods

Calculations

The system of quadratic equations was solved using the SDP solver SDPA[37]. To find the global optima of SDP problems, interior point methods were efficiently used[38].

Experiments

The polarized neutron elastic scattering experiments were performed using the HB1 thermal neutron three-axis spectrometer located at the High Flux Isotope Reactor (HFIR) of the Oak Ridge National Laboratory (ORNL). Heusler alloy 111 reflection crystals were utilized as the monochromator and analyser. The flipping ratio R = 10 was obtained using the nuclear 222 reflection at a paramagnetic temperature of 40 K corresponding to the beam polarization P = 0.82. The polarization vector was set parallel to the scattering vector (//). The incident energy of neutrons was Ei = 13.5 meV. The horizontal collimator sequence was 48′(open)–80′–80′–240′. A pyrolytic graphite Bragg-reflection filter was used to efficiently eliminate the contamination caused by higher-order wavelengths.

A powder Nd2Ir2O7 sample was synthesized by a solid-state reaction method inside a quartz tube. About 4.8 g of the sample was wrapped in thin aluminium foil and shaped to a hollow cylinder with a thickness of 0.8 mm and diameter of 20 mm to mitigate the effect of the strong neutron absorption of Ir nuclei. The cylinder was stored in an aluminium container filled with He gas and placed under the cold heads of a He-closed-cycle (Displex) or liquid-He-type (Orange) cryostat.