# Fileset

[ydc4-cfrm.pdf](https://mdr.nims.go.jp/filesets/1a1aa933-0eab-4d58-9fa4-0c511f5a9d03/download)

## Creator

[Hiroshi Oike](https://orcid.org/0000-0001-6866-7774), [Hidemaro Suwa](https://orcid.org/0000-0002-6026-2630), Yasunori Takahashi, [Fumitaka Kagawa](https://orcid.org/0000-0002-1763-6799)

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Thermally quenched metastable phase in the Ising model with competing interactions](https://mdr.nims.go.jp/datasets/47270ca5-7f50-481d-84c2-7d92de403d1a)

## Fulltext

Thermally quenched metastable phase in the Ising model with competing interactionsPHYSICAL REVIEW B 112, 064409 (2025)Thermally quenched metastable phase in the Ising model with competing interactionsHiroshi Oike ,1,2,* Hidemaro Suwa ,3 Yasunori Takahashi,4 and Fumitaka Kagawa 5,61PRESTO, Japan Science and Technology Agency (JST), Saitama 332-0012, Japan2Research Center for Materials Nanoarchitechtonics (MANA), National Institute for Materials Science (NIMS), Ibaraki 305-0047, Japan3Department of Physics, The University of Tokyo, Tokyo 113-0033, Japan4Department of Applied Physics and Quantum-Phase Electronics Centre (QPEC), The University of Tokyo, Tokyo 113-8656, Japan5Department of Physics, Institute of Science Tokyo, Meguro 152-8551, Japan6RIKEN Center for Emergent Matter Science (CEMS), Wako 351-0198, Japan(Received 16 April 2025; revised 2 July 2025; accepted 9 July 2025; published 6 August 2025)Thermal quenching has been employed to discover various metastable materials such as hard steels andmetallic glasses. More recently, quenching-based phase control has been applied to correlated electron systemsthat exhibit metal-insulator, magnetic, or superconducting transitions. Despite the discovery of metastableelectronic phases, however, the microscopic origin of metastability remains elusive, as the system’s degreesof freedom—such as electrons—can vary even at low temperatures. Here, using Monte Carlo simulations,we demonstrate a thermally quenched metastable phase in the Ising model that does not conserve the totalmagnetization. When multiple types of interactions that stabilize different long-range orders are introduced,the ordering kinetics exhibits significant slowing down, and the characteristic timescale for ordering divergesas the temperature decreases. As a result, the system can reach low temperatures without undergoing orderingif the cooling rate is sufficiently high. Quantitative analysis of the divergent behavior suggests that the energybarrier for eliminating the local structure of competing orders is the origin of this metastability. Thus, the presentsimulation demonstrates that competing interactions play an essential role in achieving metastability.DOI: 10.1103/ydc4-cfrmI. INTRODUCTIONThe thermodynamically most stable phase of a material isuniquely determined when thermodynamic parameters suchas the temperature, pressure, and chemical composition arefixed. In contrast, multiple metastable phases may exist fora given thermodynamic parameter. For example, although anaggregate of carbon atoms has graphite as the most stablephase under ambient conditions, several allotropes, such asdiamond and lonsdaleite, exist as metastable phases [1]. Usingmetastable phases as an exploration space, solid-state chem-istry and metallurgy have been developed [2–6]. A metastablephase is often achieved via thermal quenching, as exemplifiedby the glassy state that occurs when a liquid is cooled sorapidly that crystallization is kinetically avoided [7–10]. Morerecently, quenching-based phase control has been applied tocorrelated electron materials, realizing metastable electronicphases such as charge glasses [11–13], charge density waves[14], magnetic skyrmions [15–18], orbital disordered phases[19], ferromagnetic phases [20,21], and superconductivity[22]. Furthermore, the metastable phase of correlated elec-*Contact author: oike.hiroshi@nims.go.jpPublished by the American Physical Society under the terms of theCreative Commons Attribution 4.0 International license. Furtherdistribution of this work must maintain attribution to the author(s)and the published article’s title, journal citation, and DOI.trons is also realized through intense excitation of electronsusing pulsed lasers [23,24]. Thus, metastability is becoming awidely researched topic in materials science.Along with the development of methods for makingglasses, their design principles have also been clarified. Su-percooled liquids tend to have a high glass forming abilitynear a eutectic point [Fig. 1(a)], as exemplified by moltenalloys [8], pressurized silicon [25], and salt-added water [26].The microscopic mechanism underlying this behavior hasbeen investigated [27–29]. A similar guiding principle maybe applied to metastable electronic phases. However, thisprinciple is not straightforward because crystallization differsfrom metal-insulator transitions or magnetic phase transitionsin terms of the dynamics of the microscopic components.Atomic diffusion often divergently slows at low tempera-tures according to the Vogel-Fischer law [30], resulting inmetastabilization of the supercooled liquid. In contrast, inmetal-insulator transitions and magnetic phase transitions,the internal degrees of freedom of the electrons can vary tominimize the free energy of the electronic system even atlow temperatures. Therefore, a thermally quenched electronicphase may not be metastable.The present study aims to generalize the dynamics andmetastable states near a eutectic-like triple point. The Isingmodel is suitable for this purpose because it is a typicalstatistical mechanics model for describing phase transitionsinvolving the rearrangement of atoms [31,32] and the inter-nal degrees of freedom of electrons [33,34]. When eutecticcrystallization is modeled with the Ising model, one element2469-9950/2025/112(6)/064409(11) 064409-1 Published by the American Physical Societyhttps://orcid.org/0000-0001-6866-7774https://orcid.org/0000-0002-6026-2630https://orcid.org/0000-0002-1763-6799https://ror.org/00097mb19https://ror.org/026v1ze26https://ror.org/057zh3y96https://ror.org/057zh3y96https://ror.org/05dqf9946https://ror.org/03gv2xk61https://crossmark.crossref.org/dialog/?doi=10.1103/ydc4-cfrm&domain=pdf&date_stamp=2025-08-06https://doi.org/10.1103/ydc4-cfrmhttps://creativecommons.org/licenses/by/4.0/OIKE, SUWA, TAKAHASHI, AND KAGAWA PHYSICAL REVIEW B 112, 064409 (2025)(a) (b)Chemical potential, Pressure...TemperatureLiquid phaseCrystalline phase ACrystalline phase BEutectic pointJ1J3JPJ22-by-1 order 4-by-4 order2-by-1 order 4-by-4 orderphase phaseDisordered phaseTJ30.0 0.5 1.0 1.50246(c) J1 = 1.0J2 = 1.5JP = −1.0first orderfirst orderfirst orderfirst orderFIG. 1. (a) Typical eutectic phase diagram. (b) Schematic diagram of the interactions between spins in the Ising model and the orderedphases of interest. (c) T -J3 phase diagram where J1, J2, and Jp are fixed to 1.0, 1.5, and −1.0, respectively. Since almost all crystallizationprocesses are first-order phase transitions, the parameters were set so that the ordered phase would also be formed by a first-order phasetransition (see Appendix A). A eutecticlike triple point exists at parameters of (T, J3) = (2.0, 0.75).is mapped to the up spin and another to the down spin. Theconservation of the number of atoms corresponds to the con-servation of magnetization, which is defined by the differencein the numbers of up and down spins. This conservation lawimposes constraints on changes in the spin configuration,resulting in energy barriers in the process corresponding tothe exchange of atoms and the movement of an atom intoa vacancy [35,36]. In contrast to eutectic crystallization, theconservation laws are sometimes broken in correlated electronsystems, such as in the phase transition from paramagneticto ferromagnetic phases. Therefore, in the present study, ourfocus is on the nonconserved Ising model, in which the spinof each site can change independently. The pioneering worksclarified the metastability in the phase transition between or-dered phases [37,38] and the metastability of a supercooleddisordered phase in the Ising model with plaquette interac-tions [39–42], but the metastability of the disordered phasenear a eutectic-like triple point has not been clarified. Never-theless, we find that a thermally quenched disordered phaseexhibits metastability near the eutecticlike triple point in theIsing model with competing interactions that favor differ-ent ordered structures. This result suggests that the designprinciples for glass are also applicable to correlated electronsystems.The remainder of this paper is organized as follows.Section II describes the simulation method. In Sec. III, wepresent the emergence of a metastable phase, followed bythe temperature-dependent ordering kinetics of the systemin Sec. IV. In Sec. V, we discuss the ordering kinetics andstructural factors to explore the relationship between the en-ergy competition of two ordered phases and the metastabilityof a disordered phase. Section VI discusses the implicationsof our results for metastable materials and strongly corre-lated electron systems. Finally, Sec. VII provides concludingremarks.II. SIMULATION METHODTo implement a eutecticlike triple point in the phase di-agram of the Ising model, we consider a two-dimensionalsquare lattice and introduce the Hamiltonian with severaltypes of interactions following a previous study [43]:H = J1∑n.nσiσ j +J2∑2ndσiσ j +J3∑3rdσiσ j +Jp∑ringσiσ jσkσl ,(1)where J1, J2, J3, and Jp are the nearest-neighbor, next-nearestneighbor, third-nearest neighbor, and four-body ring interac-tions, respectively, as shown in Fig. 1(b). σi represents a spin,which indicates the local configuration at lattice site i andtakes a value of +1 or −1. J1 favors a checkerboard antiferro-magnetic order, J2 favors a 2-by-1 order, and J3 favors a 4-by-4order. Jp represents the many-body spin interactions arisingfrom higher-order effects of electron itinerancy [44,45] andcontributes to the ordering temperature and discontinuity ofthe phase transition. When the interaction parameters (J1, J2,Jp) are fixed to (1.0, 1.5, −1.0), the most stable ordered phasechanges from a 2-by-1 order to a 4-by-4 order with variationof J3, as shown in Fig. 1(c). The disordered, 2-by-1, and4-by-4 phases meet at the parameters of (T, J3) = (2.0, 0.75),which is a eutectic-like triple point. Thus, our simple Hamilto-nian successfully realizes a phase diagram with a eutecticliketriple point.The isothermal phase transition kinetics are then simulatedat each point of the T -J3 phase diagram. The σi values in aninitial state are randomly set to 1 or −1 to simulate a stateimmediately after thermal quenching from an infinite temper-ature with an infinite cooling rate. The time evolution of σiis simulated by the single spin-flip and heat-bath methods, inwhich the values of σi are updated after a Monte Carlo stepwith the probabilitiesP(σi = 1) = exp[−H (σi = 1)/T ]Z, (2)P(σi = −1) = exp[−H (σi = −1)/T ]Z, (3)Z = exp[−H (σi = 1)/T ] + exp[−H (σi = −1)/T ]. (4)During one Monte Carlo step, all the sites are updated in arandom order. In this study, a Monte Carlo step was con-sidered a unit of time, as in previous studies [38,46,47]. Weimposed periodic boundaries and performed simulations andanalyses. In the main text, we present the results for a system064409-2THERMALLY QUENCHED METASTABLE PHASE IN THE … PHYSICAL REVIEW B 112, 064409 (2025)(c)(a)t = 0 10 102 103 104 105T = 1.6t = 0 10 102 103 104 105T = 0.1J3 = 0.7J3 = 0.7(b)t = 0 10 102 103 104 105T = 1.6J3 = 0.7FIG. 2. (a) Isothermal time evolution of the spin configuration at parameters of (T, J3) = (1.6, 0.7). The system size N is 642. (b) Three-colored representation of the time evolution shown in panel (a). (c) Three-colored representation of the time evolution at parameters of (T, J3) =(0.1, 0.7). The red and blue areas represent 2-by-1 and 4-by-4 orders, respectively, and the white area represents all other configurations. Themethod of converting the spin configurations to three colors is described in the main text.size N of 642, which is large enough to study the fundamentalproperties, as discussed in Appendix B.III. EMERGENCE OF METASTABILITYTo examine metastability, the isothermal time evolutionis studied for the parameters of (T, J3) = (1.6, 0.7) and(0.1, 0.7), where J3 is slightly decreased from the criticalvalue of the eutectic point, 0.75, so that the 2-by-1 orderbecomes the most stable phase. At a temperature of T = 1.6,stripe-structured domains representing the 2-by-1 order ap-pear, and the domains gradually increase in size, becoming auniform 2-by-1 order as time evolves [Fig. 2(a)]. To make thephase transition process with time easier to understand, thethree types of local structures (2-by-1 order, 4-by-4 order, anddisordered) are colored differently according to the followingprocedure. First, a site is selected, and 3-by-3 sites containingthe selected site at the center are extracted. The selected site iscolored red if the 3-by-3 sites match the 2-by-1 order, coloredblue if they match the 4-by-4 order, or colored white if they donot match either of these orders. Thus, many small 2-by-1 and4-by-4-order domains appear in the initial stage, and the 2-by-1 domains coarsen over time until they become comparableto the system size [Fig. 2(b)]. The ordering kinetics from thesupercooled disordered phase presented below were obtainedthrough this procedure.In contrast to that at T = 1.6, the real-space configurationat T = 0.1 is approximately independent of time from t = 10to t = 105 [Fig. 2(c)]. This time-independent behavior meansthat the disordered phase is stable on the timescale of thesimulation. This disordered phase is considered metastablebecause the free-energetically most stable phase is a uniform2-by-1 order. The comparison of the configurations at t = 0and t = 105 indicates that the metastable phase is not a ran-dom state of the initial state but a fine domain structure of2-by-1 and 4-by-4 orders. The formation of a fine structuresuggests that local energy optimization occurs in the earlystages but that global optimization has not yet occurred.IV. TIME-TEMPERATURE-TRANSFORMATION(TTT) DIAGRAMTo obtain an overview of the temperature dependence, weexamine the TTT diagram at J3 = 0.7, as thermally quenchedmetastable phases are often discussed using the TTT dia-gram [7,13,48], which summarizes the temperature-dependentisothermal time evolution of an order parameter. To character-ize the extent to which the ordering progresses, we define theorder parameters for the 2-by-1 and 4-by-4 phases as follows:φ(2by1) = 1N√√√√√∑q=qx,qy∣∣∣∣∣N∑j=1σ jeiq·r j∣∣∣∣∣2, (5)φ(4by4) = 1N√√√√√∑q=qp,qm,−qp,−qm∣∣∣∣∣N∑j=1σ jeiq·r j∣∣∣∣∣2, (6)where N is the number of sites and qx, qy, qp, and qm are(π, 0), (0, π ), (π/2, π/2), and (π/2,−π/2), respectively.These definitions lead to φ(2by1) and φ(4by4) being equalto 1 in the uniformly ordered states and approximately 0 inan initial state with a random configuration. The isothermaltime evolution for several temperatures is obtained by aver-aging 64 individual runs, as shown in Fig. 3(a). By making acontour map of the time evolution of φ(2by1) up to t = 105for T = 0.001, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2.0, and2.2, a TTT diagram is obtained, as shown in Fig. 3(b).As shown in Fig. 3(b), φ(2by1) rapidly develops over timewithin a temperature range of 1.5 < T < 2.0, whereas itsevolution dynamics slow at higher and lower temperatures.064409-3OIKE, SUWA, TAKAHASHI, AND KAGAWA PHYSICAL REVIEW B 112, 064409 (2025)FIG. 3. (a) Isothermal time evolution of order parameters φ(2by1) and φ(4by4) at parameters of (T, J3) = (0.7, 0.7). τst and τcp representthe times at which the order parameter of the stable order exceeds 0.1 and that the order parameter of the competing order decreases to halfits value at t = 1, respectively. The stable order is a 2-by-1 order for J3 less than 0.75 and a 4-by-4 order for J3 greater than 0.75, with thecompeting order being the opposite. (b) TTT diagram of the Ising model at a J3 of 0.7. (c) Arrhenius plot of the growth time of the stable 2-by-1order τst and the destruction time of the competing 4-by-4 order τcp. �st and �cp were extracted from the slopes of the Arrhenius plots of τstand τcp, respectively. (d) Schematic representation of the time evolution. From the initial random arrangement, fine 2-by-1 and 4-by-4 domainsrapidly form, and when the process of overcoming the energy barrier �st , which is comparable to �cp, is activated by thermal fluctuations, auniform 2-by-1 order is formed.In particular, within a temperature range of T < 0.4, theordering hardly progresses, and the system does not reachthe most stable 2-by-1 phase, at least within the timescaleof the simulation, indicating the emergence of metastability.This nonmonotonic temperature dependence is commonly ob-served in supercooled liquids and arises from two competingtemperature dependencies: the energy barriers in the phasetransition pathway are less likely to be overcome at low tem-peratures, whereas the free-energy driving force for orderingbecomes greater at low temperatures. The TTT diagram forthe present Ising model suggests the existence of an energybarrier.To characterize the timescale over which a stable orderedphase is formed, we define τst as the time at which the orderparameter of the stable order exceeds 0.1. The examinationof τst at various temperatures reveals that τst follows anArrhenius-type behavior [Fig. 3(c)], indicating the presenceof a thermal activation barrier �st [Fig. 3(d)]. For supercooledliquids, the diffusion of atoms is widely believed to contributeto the energy barrier [35,36]. However, note that a finite �stappears, although individual degrees of freedom are vari-able without a diffusion process in the present nonconservedsystem.To gain insight into the energy barrier, the dynamics ofthe disappearing competing order are characterized by thetime τcp at which φ(4by4) is half of the value at t = 1. Theexamination of τcp at various temperatures reveals that τcpalso shows an Arrhenius-type temperature dependence witha thermal activation barrier �cp [Fig. 3(c)]. A possible reasonfor the close values of �st and �cp is that the rate-limitingprocess of the phase transition is breaking competing orders.V. ENERGETIC COMPETITION AND METASTABILITYBecause J3 controls the relative stability of the 2-by-1 and4-by-4 orders, we then examine the J3 dependence of thephase transition kinetics and metastability. To overview theJ3 dependence, the values of φ(2by1) − φ(4by4) are plottedin Fig. 4(a) as contour maps in the J3-T plane at t = 1,10, 102, 103, 104, 105. As the 2-by-1 (4-by-4) order grows,the red (blue) color becomes clearer in the contour maps[Fig. 4(a)]. At a J3 close to 0 or 1.5, the ordering appreciablyprogresses at t = 10 and is almost complete at t = 103. AsJ3 becomes closer to that at the eutecticlike point (J3 = 0.75),both φ(2by1) and φ(4by4) remain at values less than 0.1 evenat t = 105 within a temperature range of T < 0.5.To characterize the region in which the metastable phasesexist in the J3-T phase diagram, we define the upper tempera-ture limit of metastability Tm as the temperature below whichφ(2by1) and φ(4by4) are less than 0.1 at t = 105, and the J3dependence of Tm is also plotted in the J3-T phase diagram att = 105 [Fig. 4(a)]. Tm is found to take a maximum value nearJ3 = 0.75, indicating that the system tends to be metastablenear the phase boundary between the 2-by-1 and 4-by-4phases. The energy barriers �st and �cp also take a maximumnear J3 = 0.75 [Fig. 4(b)] and are positively correlated with064409-4THERMALLY QUENCHED METASTABLE PHASE IN THE … PHYSICAL REVIEW B 112, 064409 (2025)FIG. 4. (a) Contour plot of the time evolution of the order parameter in the J3 range of 0–1.5. φ(2by1) − φ(4by4) is plotted, as the formationof the 2-by-1 order is indicated by a positive value (red) and the formation of the 4-by-4 order is indicated by a negative value (blue). Tm is thetemperature below which φ(2by1) or φ(4by4) is less than 0.1 at t = 105. (b) J3 dependence of the energy barriers �st and �cp estimated withthe Arrhenius plot shown in Fig. 3(c). (c) Correlations between Tm and the energy barriers �st and �cp. A guide to the eye is used to show theproportional relationship.Tm [Fig. 4(c)]. Thus, the J3-dependent energy barrier is foundto play a key role in controlling the metastability in our modelsystem with a eutecticlike triple point.The good correspondence between the evolution of theenergy barriers and that of the metastability suggests thatthe local configuration is correlated with the metastability.To address this issue in more detail, below we examine thedisordered metastable state in terms of the structure factor.At a temperature of T = 0.1 and J3 values of 0.4, 0.7, 0.8,and 1.1, the metastable disordered phases contain both 2-by-1and 4-by-4 domains, the fractions of which depend on J3[Figs. 5(a)–5(d)]. To quantify the domain fractions, the real-space configurations are converted into the structure factorS(q) asS(q) = 1N∣∣∣∣∣∣N∑j=1σ jeiq·r j∣∣∣∣∣∣2, (7)as shown in Figs. 5(e)–5(h). The intensities of S(q) near(±π, 0) and (0,±π ) correspond to the fraction of the 2-by-1domains [Fig. 5(i)], and those near (±π/2,±π/2) correspondto the fraction of the 4-by-4 domains [Fig. 5(j)]. The do-main fractions of the 2-by-1 and 4-by-4 orders, f (2by1) andf (4by4), are estimated by integrating S(q) over the circularareas near these q values, as shown in Figs. 5(i) and 5(j).Whereas φ(2by1) and φ(4by4) reflect S(q) only for the qvalues corresponding to the uniformly ordered structure,f (2by1) and f (4by4) integrate S(q) around the q values. Be-cause the Fourier transformation of a small domain includesS(q) around the q values, f (2by1) and f (4by4) are consideredmore suitable for evaluating the fractions of 2-by-1 and 4-by-4domains in the fine domain structure. f (2by1) and f (4by4)are found to vary monotonically as a function of J3 [Fig. 5(k)],reflecting the J3 dependence of the domain structures, asshown in Figs. 5(a)– 5(d).Having clarified the microscopic structure of themetastable phase, we discuss how the energy barrier foundin the simulation can be evaluated on the basis of themicroscopic structure. We consider a single spin-flip to breakthe uniform competing order, which corresponds to Fig. 5(l)for J3 = 0.75 − 1.5 and to Fig. 5(m) for J3 = 0 − 0.75. Theenergy barrier �Eo required for these processes is calculatedand plotted in Fig. 5(n) and is found to be much largerthan �st . This discrepancy in the energy barrier reflectsthe fact that the single spin-flips in the uniformly orderedphase are not a direct factor determining the robustnessof the metastability. Because the metastable phases thatwe are discussing are highly disordered and the locallyformed order is often surrounded by disordered states, thediscrepancy between �Eo and �st is rather reasonable. Toaccount for the probability that a locally formed competingorder is surrounded by disordered states, the energy barrieris scaled by the domain fraction as �Es = f (4by4)�Eo for064409-5OIKE, SUWA, TAKAHASHI, AND KAGAWA PHYSICAL REVIEW B 112, 064409 (2025)FIG. 5. (a)–(d) Three-colored representation of the spin configurations at T = 0.1 and t = 104 for J3 = 0.4 (a), 0.7 (b), 0.8 (c), and 1.1(d). (e)–(h) Structure factor S(q) of the metastable phase shown in panels (a)–(d). (i), (j) Structure factor S(q) of the uniform 2-by-1 order (i)and 4-by-4 order (j). (k) J3 dependence of the estimated domain fractions of the 2-by-1 and 4-by-4 orders, f (2by1) and f (4by4). See the maintext and Appendix C for details of the estimation method. (l), (m) Single spin-flip process for the uniform 2-by-1 order (l) and 4-by-4 order(m). (n) J3 dependence of the energy barriers �Eo and �Es. �Eo is the energy increase due to the single spin-flip processes shown in panel (l)for 0.75 < J3 < 1.5 and panel (m) for 0 < J3 < 0.75. The scaled energy barrier �Es equals f (4by4)�Eo for 0 < J3 < 0.75 and f (2by1)�Eofor 0.75 < J3 < 1.5.0 < J3 < 0.75 and f (2by1)�Eo for 0.75 < J3 < 1.5. We thenfind that values of �Es are close to those of �st [Fig. 5(n)].The correspondence between the scaled energy barrier �Esand that obtained by the Arrhenius law �st suggests a methodfor quantitatively estimating energy barriers from microscopicstructures.VI. DISCUSSIONSIn the mean-field approximation, the metastable solution ofa disordered phase is lost when the system is supercooled tobelow the spinodal point [49]. Similarly, in the present Isingmodel at J = 0.7, the local minimum close to φ(2by1) = 0 islost at T = 2.18, which is only a few percent below Tc [seeFig. 6(c) in Appendix A]. However, the disappearance of ametastable solution does not necessarily indicate the absenceof a practically stable phase. Indeed, the present study showsthat a practically metastable phase exists even when the sys-tem is deeply supercooled. This type of metastability is oftenobserved in solid solutions produced by rapid cooling. Forexample, in the TiO2-VO2 system, spinodal decompositionoccurs at temperatures below 830 K, but under rapid cooling,a stable solid solution with a uniform distribution of Ti and Vions can be obtained at room temperature [50], meaning thatspinodal decomposition is kinetically arrested. Such kineticarrest of spinodal decomposition is widely accepted in thecontext of conserved systems, as the diffusion processes ac-companying decomposition reasonably account for the energybarriers responsible for the arrest. In nonconserved systems,however, whether kinetic arrest occurs is nontrivial, sincethe spinodal growth does not require decomposition, and theenergy barriers depend strongly on the underlying interactionsof many-body systems. The appearance of a metastable phasein the present Ising model demonstrates that kinetic arrest ofspinodal growth can also occur in a nonconserved system,implying the universality of kinetic arrest beyond obviousfree-energy minima.Finally, we discuss the possible relevance of the presentsimulations for correlated electron materials. Transition metaloxides are typical examples of materials with eutectic-type064409-6THERMALLY QUENCHED METASTABLE PHASE IN THE … PHYSICAL REVIEW B 112, 064409 (2025)geometries in electronic phase diagrams [51–53]. In man-ganites with a composition in which two electronic phasescompete, a glassy electronic state appears under the coolingprocess with a standard cooling rate due to the effect of ran-domness [52]. This sensitivity to randomness may indicate atendency for a metastable disordered phase to appear underthermal quenching. In terms of the energetic competition be-tween multiple ordered phases, geometrical frustration, as inthe case of triangular lattices, can also be effective in realizingmetastable phases. For example, in quarter-filled triangularlattice systems θ -(BEDT-TTF)2X , several charge-ordered pat-terns are energetically close to each other, indicating energeticcompetition [54–57]. Metastable charge glass phases are lesslikely to appear under a cooling process with a certain coolingrate when the triangular lattice is distorted [12,58], indicatingthat frustration is effective for inducing metastable phases.Thus, frustrated systems, which have attracted much attentionfor the exploration of exotic quantum phases, are promisingfor the realization of metastable electronic phases.VII. CONCLUDING REMARKSThe present study aims to generalize the principle ofmetastability established in metallurgy with the goal of apply-ing it to correlated electron systems. We have demonstratedthat a thermally quenched disordered phase exhibits metasta-bility near the eutecticlike triple point in the Ising model withcompeting interactions. This simple model of metastabilityenables us to discuss how metastability varies with the inter-action strength and how it correlates with local structures. Thenature of metastability revealed by the model is likely relevantto experimentally observed metastable electronic phases, sug-gesting a link to correlated electron systems.Microscopic simulations address phenomena of a limitedtime range due to computational constraints. Therefore, theappearance of a metastable phase may simply indicate that thetime range is insufficient to reach the most stable phase. In ad-dition, the physics research of correlated electron systems hasfocused mainly on their most stable phases. Thus, numericallyfound metastable phases have not received much attention.However, many metastable materials have been created byexploring states before the most stable state is reached. Ifwe focus on the physics of metastability, then a numericallyobtained metastable phase would become a valuable researchtarget.ACKNOWLEDGMENTSH.O. thanks H. Ohtani, M. Enoki, K. Nakanishi, and H.Kageyama for fruitful discussions. Calculations were per-formed using computational resources of the SupercomputerCenter at the Institute for Solid State Physics, the Univer-sity of Tokyo. This work was supported by JST PRESTO(Grant No. JPMJPR21Q2), the Center of Innovation for Sus-tainable Quantum AI (SQAI), JST (Grant No. JPMJPF2221),JST CREST (Grant No. JPMJCR24I1), and JSPS KAKENHI(Grants No. 22H01164, No. 23K22435, No. 23H04861, No.22K03508, and No. 24H01609). MANA is supported byWorld Premier International Research Center Initiative (WPI),MEXT, Japan.FIG. 6. (a), (b) Temperature dependence of the order param-eters φ(2by1) for J3 = 0.0 − 0.7 (a) and φ(4by4) for J3 = 0.8 −1.5 (b). (c) Free-energy landscape for J3 = 0.7 near the orderingtemperature. The arrows indicate the minimum free energy nearφ(2by1) = 0.DATA AVAILABILITYThe data that support the findings of this article are notpublicly available. The data are available from the authorsupon reasonable request.APPENDIX AThe phase transition temperature of the Ising model wasdetermined from the temperature dependence of the orderparameter when the free energy reaches a minimum value(Fig. 6). When J3 is less than 0.75, a 2-by-1 order forms atthe phase transition temperature, whereas when it is greaterthan 0.75, a 4-by-4 order forms. The order parameter showsa discontinuous change at the phase transition temperature,064409-7OIKE, SUWA, TAKAHASHI, AND KAGAWA PHYSICAL REVIEW B 112, 064409 (2025)T = 0.6J3 = 0.7 t = 103L = 64 L = 128FIG. 7. Spin configurations at parameters of (T, J3) = (0.6, 0.7),t = 103, and L = 64 and 128.indicating a first-order phase transition. The free energy wascalculated using the Wang-Landau algorithm with the systemsize N = L2 = 322 [59]. The behavior discussed in Fig. 6(c),in which the local minimum near φ(2by1) = 0 disappears be-tween T = 2.18 and T = 2.2, was also observed for L = 24.Thus, L = 32 is large enough to reliably capture the temper-ature dependence of the free-energy minimum. Furthermore,using L = 64, we confirmed the order-parameter jump at thefirst-order transition temperature, shown in Figs. 6(a) and6(b). The temperatures at which the order parameter jumpsand at which the minimum point of the free energy shifts inFig. 6(c) are almost the same, indicating that the size depen-dence is sufficiently small.APPENDIX BTo check the system-size dependence, a comparison be-tween system sizes N (= L2) of 642 and 1282 is shown inFig. 7. Although the typical domain size is not dependent onthe system size, the order parameter is by definition depen-dent on the system size for the following reasons. The orderparameter φ(2by1) can be rewritten asφ(2by1) =√√√√√ 1N⎛⎝ ∑q=qx,qyS(q)⎞⎠, (B1)and the structure factor S(q) can be rewritten asS(q) = 1NN∑j,k=1〈σ jσk〉eiq·(r j−rk ), (B2)where 〈σ jσk〉 is the correlation function of spins at sites j andk. σ jσkeiq·(r j−rk ) is 1 if sites i and j are in the same domain ofthe ordered structure with wave number q and can otherwisetake a random value whose absolute value is 1. Therefore, thevalue obtained by taking the sum of this quantity over k isapproximately equal to N jq , which is the domain size of theordered structure with wave number q containing site j. Thestructural factor can then be expressed approximately asS(q) ∼ 1NN∑j=1〈N jq〉. (B3)FIG. 8. (a), (b) Isothermal time evolution of φ(2by1) for(T, J3) = (0.6, 0.7) (a) and (T, J3) = (0.6, 0.7) (b). (c), (d) Isother-mal time evolution of Lφ(2by1) for (T, J3) = (0.6, 0.7) (c) and(T, J3) = (0.6, 0.7) (d). The size dependence was examined atL = 64, 128, 256, and 512.The domain size and the number of domains are approxi-mately proportional to ξ 2 and N/ξ 2, respectively, where ξis the correlation length. The structure factor can then beapproximated by taking the sum over j asS(q) ∼ ξ 2. (B4)The order parameter is then expressed asφ(2by1) ∼ ξL. (B5)Thus, φ(2by1) depends on the system size when ξ doesnot. We observed that Lφ(2by1) remains independent of thesystem size up to the timescale at which Lφ(2by1) ∼ ξ ap-proaches L (Fig. 8). When the typical domain size is equalto the system size, the order estimate described above isno longer valid. Therefore, this scaling does not hold whenφ(2by1) is close to 1. In the present study, we focus on thedynamics when φ(2by1) is less than 0.5 to discuss the prop-erties that are independent of the system size. In this φ(2by1)range, Lφ(2by1) shows similar behavior for L = 64, 128, 256,and 512, as seen in the TTT diagrams (Fig. 9), so we examinedthe T -J3 dependence in detail for L = 64.As mentioned above, the present simulations investigatethe growth of ξ . According to previous studies [60,61], scalinganalysis of the free energy shows that the bulk term scalesas Ld , while the surface term scales as Ld−1, where d isthe spatial dimension. Drawing an analogy from this scalinganalysis, by substituting L with ξ and setting d = 2, the bulkfree energy scales as NDξ 2, and the interfacial energy betweendomains scales as NDξ , where ND is the number of stable064409-8THERMALLY QUENCHED METASTABLE PHASE IN THE … PHYSICAL REVIEW B 112, 064409 (2025)FIG. 9. (a)–(d) TTT diagram plotted based on Lφ(2by1) forJ3 = 0.7 and L = 64 (a), 128 (b), 256 (c), and 512 (d).ordered domains. Since ND ∼ (L/ξ )2, the bulk and interfacialenergy contributions scale as L2 and L2/ξ , respectively.This scaling analysis suggests that the interfacial free en-ergy decreases monotonically as ξ increases, while the bulkcontribution remains constant. The monotonic decrease offree energy with increasing ξ aligns with the monotonicgrowth of ξ over time observed in our simulations. However, itdoes not account for the thermal activation barriers associatedwith quantities such as �st , �cp, �Eo, and �Es. Microscopicsimulations—such as those analyzing the internal structureof interfaces in our study—can provide valuable insight intothese thermal activation barriers.APPENDIX CThe fractions of domains of 2-by-1 and 4-by-4 orders inthe disordered phase were estimated from the structure factoras follows. In a uniform 2-by-1 ordered state, S(q) is equalto N if q is either (π, 0) or (0, π ), which are the pointsindicated in Fig. 4(i), and 0 if q is any other value. In a uniform4-by-4 ordered state, S(q) is equal to N if q is (π/2, π/2),(π/2,−π/2), (−π/2, π/2), or (−π/2,−π/2), which are thepoints indicated in Fig. 4(j), and 0 if q is any other value.When the 2-by-1 and 4-by-4 orders form a domain structure,S(q) takes a finite value around the values of q correspondingto the 2-by-1 and 4-by-4 orders. The values obtained by inte-grating the S(q) of these diffuse spots over the q range shownin Figs. 4(i) and 4(j) indicates how much area the 2-by-1and 4-by-4 domains occupy. Therefore, by normalizing thesevalues by N , the domain fractions of the 2-by-1 and 4-by-4orders are estimated asf (2by1) = 1N∑q∼qx,qyS(q), (C1)f (4by4) = 1N∑q∼qp,qm,−qp,−qmS(q). (C2)[1] B. Derjaguin, D. Fedoseev, V. Varnin, and S. Vnukov, Thenature of metastable phases of carbon, Nature (London) 269,398 (1977).[2] W. Sun, S. T. Dacek, S. P. Ong, G. Hautier, A. Jain, W. D.Richards, A. C. Gamst, K. A. Persson, and G. Ceder, Thethermodynamic scale of inorganic crystalline metastability,Sci. Adv. 2, e1600225 (2016).[3] W. Sun, A. Holder, B. Orvañanos, E. Arca, A. Zakutayev, S.Lany, and G. Ceder, Thermodynamic routes to novel metastablenitrogen-rich nitrides, Chem. Mater. 29, 6936 (2017).[4] A. J. Martinolich and J. R. Neilson, Toward reaction-by-design:achieving kinetic control of solid state chemistry with metathe-sis, Chem. Mater. 29, 479 (2017).[5] M. Bykov, S. Chariton, H. Fei, T. Fedotenko, G. Aprilis,A. V. Ponomareva, F. Tasnádi, I. A. Abrikosov, B. Merle, P.Feldner et al., High-pressure synthesis of ultraincompressiblehard rhenium nitride pernitride Re2(N2)(N)2 stable at ambientconditions, Nat. Commun. 10, 2994 (2019).[6] H. Ito, Y. Nakahira, N. Ishimatsu, Y. Goto, A. Yamashita, Y.Mizuguchi, C. Moriyoshi, T. Toyao, K.-I. Shimizu, H. Oikeet al., Stability and metastability of Li3YCl6 and Li3HoCl6,Bull. Chem. Soc. Jpn. 96, 1262 (2023).[7] D. Uhlmann, A kinetic treatment of glass formation, J. Non-Cryst. Solids 7, 337 (1972).[8] A. L. Greer, Metallic glasses, Science 267, 1947 (1995).[9] C. A. Angell, Formation of glasses from liquids and biopoly-mers, Science 267, 1924 (1995).[10] P. G. Debenedetti and F. H. Stillinger, Supercooled liquids andthe glass transition, Nature (London) 410, 259 (2001).[11] F. Kagawa, T. Sato, K. Miyagawa, K. Kanoda, Y. Tokura, K.Kobayashi, R. Kumai, and Y. Murakami, Charge-cluster glassin an organic conductor, Nat. Phys. 9, 419 (2013).[12] H. Oike, F. Kagawa, N. Ogawa, A. Ueda, H. Mori, M.Kawasaki, and Y. Tokura, Phase-change memory function ofcorrelated electrons in organic conductors, Phys. Rev. B 91,041101(R) (2015).[13] S. Sasaki, K. Hashimoto, R. Kobayashi, K. Itoh, S. Iguchi, Y.Nishio, Y. Ikemoto, T. Moriwaki, N. Yoneyama, M. Watanabeet al., Crystallization and vitrification of electrons in a glass-forming charge liquid, Science 357, 1381 (2017).[14] M. Yoshida, Y. Zhang, J. Ye, R. Suzuki, Y. Imai, S. Kimura,A. Fujiwara, and Y. Iwasa, Controlling charge-density-wavestates in nano-thick crystals of 1T-TiSe2, Sci. Rep. 4, 7302(2014).064409-9https://doi.org/10.1038/269398a0https://doi.org/10.1126/sciadv.1600225https://doi.org/10.1021/acs.chemmater.7b02399https://doi.org/10.1021/acs.chemmater.6b04861https://doi.org/10.1038/s41467-019-10995-3https://doi.org/10.1246/bcsj.20230132https://doi.org/10.1016/0022-3093(72)90269-4https://doi.org/10.1126/science.267.5206.1947https://doi.org/10.1126/science.267.5206.1924https://doi.org/10.1038/35065704https://doi.org/10.1038/nphys2642https://doi.org/10.1103/PhysRevB.91.041101https://doi.org/10.1126/science.aal3120https://doi.org/10.1038/srep07302OIKE, SUWA, TAKAHASHI, AND KAGAWA PHYSICAL REVIEW B 112, 064409 (2025)[15] H. Oike, A. Kikkawa, N. Kanazawa, Y. Taguchi, M. Kawasaki,Y. Tokura, and F. Kagawa, Interplay between topological andthermodynamic stability in a metastable magnetic skyrmionlattice, Nat. Phys. 12, 62 (2016).[16] G. Berruto, I. Madan, Y. Murooka, G. M. Vanacore, E.Pomarico, J. Rajeswari, R. Lamb, P. Huang, A. J. Kruchkov,Y. Togawa, T. LaGrange, D. McGrouther, H. M. Rønnow,and F. Carbone, Laser-induced skyrmion writing and erasingin an ultrafast cryo-Lorentz transmission electron microscope,Phys. Rev. Lett. 120, 117201 (2018).[17] X. Yu, D. Morikawa, T. Yokouchi, K. Shibata, N. Kanazawa, F.Kagawa, T.-H. Arima, and Y. Tokura, Aggregation and collapseof skyrmions in a non-equilibrium state, Nat. Phys. 14, 832(2018).[18] M. T. Birch, R. Takagi, S. Seki, M. N. Wilson, F. Kagawa, A.Štefančič, G. Balakrishnan, R. Fan, P. Steadman, C. J. Ottley,M. Crisanti, R. Cubitt, T. Lancaster, Y. Tokura, and P. D. Hatton,Increased lifetime of metastable skyrmions by controlled dop-ing, Phys. Rev. B 100, 014425 (2019).[19] T. Katsufuji, T. Kajita, S. Yano, Y. Katayama, and K. Ueno,Nucleation and growth of orbital ordering, Nat. Commun. 11,2324 (2020).[20] K. Matsuura, H. Oike, V. Kocsis, T. Sato, Y. Tomioka, Y.Kaneko, M. Nakamura, Y. Taguchi, M. Kawasaki, Y. Tokuraet al., Kinetic pathway facilitated by a phase competitionto achieve a metastable electronic phase, Phys. Rev. B 103,L041106 (2021).[21] K. Matsuura, Y. Nishizawa, M. Kriener, T. Kurumaji, H. Oike,Y. Tokura, and F. Kagawa, Thermodynamic determination ofthe equilibrium first-order phase-transition line hidden by hys-teresis in a phase diagram, Sci. Rep. 13, 6876 (2023).[22] H. Oike, M. Kamitani, Y. Tokura, and F. Kagawa, Kinetic ap-proach to superconductivity hidden behind a competing order,Sci. Adv. 4, eaau3489 (2018).[23] L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D. Svetin, S.Brazovskii, and D. Mihailovic, Ultrafast switching to a stablehidden quantum state in an electronic crystal, Science 344, 177(2014).[24] I. Vaskivskyi, J. Gospodaric, S. Brazovskii, D. Svetin, P. Sutar,E. Goreshnik, I. A. Mihailovic, T. Mertelj, and D. Mihailovic,Controlling the metal-to-insulator relaxation of the metastablehidden quantum state in 1T-TaS2, Sci. Adv. 1, e1500168(2015).[25] V. Molinero, S. Sastry, and C. A. Angell, Tuning of tetra-hedrality in a silicon potential yields a series of monatomic(metal-like) glass formers of very high fragility, Phys. Rev. Lett.97, 075701 (2006).[26] M. Kobayashi and H. Tanaka, Possible link of the V-shapedphase diagram to the glass-forming ability and fragility in awater-salt mixture, Phys. Rev. Lett. 106, 125703 (2011).[27] Y. Waseda and K. Suzuki, Structure of molten silicon and ger-manium by X-ray diffraction, Z. Phys. B 20, 339 (1975).[28] T. A. Weber and F. H. Stillinger, Local order and structuraltransitions in amorphous metal-metalloid alloys, Phys. Rev. B31, 1954 (1985).[29] H. Tanaka, Relationship among glass-forming ability, fragility,and short-range bond ordering of liquids, J. Non-Cryst. Solids351, 678 (2005).[30] H. Tanaka, Possible resolution of the Kauzmann paradox insupercooled liquids, Phys. Rev. E 68, 011505 (2003).[31] T.-D. Lee and C.-N. Yang, Statistical theory of equations ofstate and phase transitions. II. Lattice gas and Ising model,Phys. Rev. 87, 410 (1952).[32] C. Van Baal, Order-disorder transformations in a generalizedIsing alloy, Physica 64, 571 (1973).[33] G. F. Newell and E. W. Montroll, On the theory of the Isingmodel of ferromagnetism, Rev. Mod. Phys. 25, 353 (1953).[34] C. Castellani, C. Di Castro, D. Feinberg, and J. Ranninger, Newmodel Hamiltonian for the metal-insulator transition, Phys. Rev.Lett. 43, 1957 (1979).[35] P. Fratzl and O. Penrose, Kinetics of spinodal decomposition inthe Ising model with vacancy diffusion, Phys. Rev. B 50, 3477(1994).[36] P. Fratzl, O. Penrose, R. Weinkamer, and I. Žižak, Coarsening inthe Ising model with vacancy dynamics, Phys. A (Amsterdam,Neth.) 279, 100 (2000).[37] P. A. Rikvold, H. Tomita, S. Miyashita, and S. W. Sides,Metastable lifetimes in a kinetic Ising model: Dependence onfield and system size, Phys. Rev. E 49, 5080 (1994).[38] P. A. Rikvold, G. Brown, S. Miyashita, C. Omand,and M. Nishino, Equilibrium, metastability, and hysteresisin a model spin-crossover material with nearest-neighborantiferromagnetic-like and long-range ferromagnetic-like inter-actions, Phys. Rev. B 93, 064109 (2016).[39] A. Lipowski and D. Johnston, Metastability in a four-spin Isingmodel, J. Phys. A: Math. Gen. 33, 4451 (2000).[40] A. Lipowski and D. Johnston, Crystallization of a supercooledliquid and of a glass: Ising model approach, Phys. Rev. E 64,041605 (2001).[41] P. Dimopoulos, D. Espriu, E. Jané, and A. Prats, Slow dynamicsin the three-dimensional gonihedric model, Phys. Rev. E 66,056112 (2002).[42] G. Savvidy, The gonihedric paradigm extension of the Isingmodel, Mod. Phys. Lett. B 29, 1550203 (2015).[43] D. Landau and K. Binder, Phase diagrams and critical be-havior of Ising square lattices with nearest-, next-nearest-,and third-nearest-neighbor couplings, Phys. Rev. B 31, 5946(1985).[44] D. Thouless, Exchange in solid 3He and the Heisenberg Hamil-tonian, Proc. Phys. Soc. 86, 893 (1965).[45] G. Misguich, C. Lhuillier, B. Bernu, and C. Waldtmann, Spin-liquid phase of the multiple-spin exchange Hamiltonian on thetriangular lattice, Phys. Rev. B 60, 1064 (1999).[46] A. J. Bray, Theory of phase-ordering kinetics, Adv. Phys. 43,357 (1994).[47] M. Naskar, M. Acharyya, E. Vatansever, and N. G. Fytas,Metastable behavior of the spin-s Ising and Blume-Capel fer-romagnets: A Monte Carlo study, Phys. Rev. E 104, 014107(2021).[48] T. Sato, K. Miyagawa, and K. Kanoda, Electronic crystalgrowth, Science 357, 1378 (2017).[49] S. Miyashita, Collapse of Metastability (Springer, Berlin, 2022),pp. 28–30.[50] Z. Hiroi, H. Hayamizu, T. Yoshida, Y. Muraoka, Y. Okamoto,J.-I. Yamaura, and Y. Ueda, Spinodal decomposition in theTiO2-VO2 system, Chem. Mater. 25, 2202 (2013).[51] Y. Tokura and Y. Tomioka, Colossal magnetoresistive mangan-ites, J. Magn. Magn. Mater. 200, 1 (1999).[52] Y. Tokura, Critical features of colossal magnetoresistive man-ganites, Rep. Prog. Phys. 69, 797 (2006).064409-10https://doi.org/10.1038/nphys3506https://doi.org/10.1103/PhysRevLett.120.117201https://doi.org/10.1038/s41567-018-0155-3https://doi.org/10.1103/PhysRevB.100.014425https://doi.org/10.1038/s41467-020-16004-2https://doi.org/10.1103/PhysRevB.103.L041106https://doi.org/10.1038/s41598-023-33816-6https://doi.org/10.1126/sciadv.aau3489https://doi.org/10.1126/science.1241591https://doi.org/10.1126/sciadv.1500168https://doi.org/10.1103/PhysRevLett.97.075701https://doi.org/10.1103/PhysRevLett.106.125703https://doi.org/10.1007/BF01313204https://doi.org/10.1103/PhysRevB.31.1954https://doi.org/10.1016/j.jnoncrysol.2005.01.070https://doi.org/10.1103/PhysRevE.68.011505https://doi.org/10.1103/PhysRev.87.410https://doi.org/10.1016/0031-8914(73)90010-4https://doi.org/10.1103/RevModPhys.25.353https://doi.org/10.1103/PhysRevLett.43.1957https://doi.org/10.1103/PhysRevB.50.3477https://doi.org/10.1016/S0378-4371(99)00527-0https://doi.org/10.1103/PhysRevE.49.5080https://doi.org/10.1103/PhysRevB.93.064109https://doi.org/10.1088/0305-4470/33/24/304https://doi.org/10.1103/PhysRevE.64.041605https://doi.org/10.1103/PhysRevE.66.056112https://doi.org/10.1142/S0217984915502036https://doi.org/10.1103/PhysRevB.31.5946https://doi.org/10.1088/0370-1328/86/5/301https://doi.org/10.1103/PhysRevB.60.1064https://doi.org/10.1080/00018739400101505https://doi.org/10.1103/PhysRevE.104.014107https://doi.org/10.1126/science.aal2426https://doi.org/10.1021/cm400236phttps://doi.org/10.1016/S0304-8853(99)00352-2https://doi.org/10.1088/0034-4885/69/3/R06THERMALLY QUENCHED METASTABLE PHASE IN THE … PHYSICAL REVIEW B 112, 064409 (2025)[53] K. Shibuya, M. Kawasaki, and Y. Tokura, Metal-insulator tran-sition in V1−xWxO2 (� x � 0.33) epitaxial thin films) epitaxialthin films, Appl. Phys. Lett. 96, 022102 (2010).[54] H. Seo, Charge ordering in organic ET compounds, J. Phys.Soc. Jpn. 69, 805 (2000).[55] Y. Pramudya, H. Terletska, S. Pankov, E. Manousakis, and V.Dobrosavljević, Nearly frozen Coulomb liquids, Phys. Rev. B84, 125120 (2011).[56] S. Mahmoudian, L. Rademaker, A. Ralko, S. Fratini, and V.Dobrosavljević, Glassy dynamics in geometrically frustratedCoulomb liquids without disorder, Phys. Rev. Lett. 115, 025701(2015).[57] L. Rademaker, Z. Nussinov, L. Balents, and V. Dobrosavljević,Suppressed density of states in self-generated Coulomb glasses,New J. Phys. 20, 043026 (2018).[58] T. Sato, F. Kagawa, K. Kobayashi, A. Ueda, H. Mori, K.Miyagawa, K. Kanoda, R. Kumai, Y. Murakami, and Y. Tokura,Systematic variations in the charge-glass-forming ability ofgeometrically frustrated θ -(BEDT-TTF)2X organic conductors,J. Phys. Soc. Jpn. 83, 083602 (2014).[59] F. Wang and D. P. Landau, Determining the density ofstates for classical statistical models: A random walk algo-rithm to produce a flat histogram, Phys. Rev. E 64, 056101(2001).[60] J. Lee and J. M. Kosterlitz, New numerical methodto study phase transitions, Phys. Rev. Lett. 65, 137(1990).[61] J. Lee and J. M. Kosterlitz, Finite-size scaling and Monte Carlosimulations of first-order phase transitions, Phys. Rev. B 43,3265 (1991).064409-11https://doi.org/10.1063/1.3291053https://doi.org/10.1143/JPSJ.69.805https://doi.org/10.1103/PhysRevB.84.125120https://doi.org/10.1103/PhysRevLett.115.025701https://doi.org/10.1088/1367-2630/aab8cehttps://doi.org/10.7566/JPSJ.83.083602https://doi.org/10.1103/PhysRevE.64.056101https://doi.org/10.1103/PhysRevLett.65.137https://doi.org/10.1103/PhysRevB.43.3265