1. Introduction
The study of disordered magnetic systems has attracted much attention in statistical physics due to the rich physics associated with frustration that disrupts an otherwise ordered state. Among various models, the ±J random-bond Ising model (RBIM) stands out as a classic example. Its Hamiltonian is given by
$\begin{eqnarray}H(\{{S}_{i}\})=-\displaystyle \sum _{\left\langle ij\right\rangle }{J}_{ij}{S}_{i}{S}_{j},\end{eqnarray}$
where Si = ± 1, i = 1, …, N are the Ising spins associated with the lattice sites, $\left\langle ij\right\rangle $ denotes all nearest neighbor pairs, and Jij are quenched random variables chosen independently from a bimodal distribution $\begin{eqnarray}P({J}_{ij})=p\delta ({J}_{ij}-J)+(1-p)\delta ({J}_{ij}+J).\end{eqnarray}$
In this work, we shall focus exclusively on the two-dimensional (2D) square lattice with J = 1 and periodic boundary conditions.Extensive work has been done to map the phase diagram of the 2D ±J RBIM on the p − T plane as shown in figure 1 (see [1–4]). There exists a special line, known as the Nishimori line (NL) [5–7], along which the energy and specific heat as a function of temperature T and ferromagnetic (FM) bond percentage p can be computed exactly. The NL is given by
$\begin{eqnarray}{{\rm{e}}}^{2{J}/{T}}=\frac{{p}}{1-{p}}.\end{eqnarray}$
It crosses the ferromagnetic phase boundary at the multicritical point (MCP), separating it into upper (Ising-like) and lower parts. While transition on the high temperature side is governed by the 2D Ising fixed point, with disorder playing only a minor role [8–10], its nature on the low-temperature side remains controversial [2, 11–13].
Figure 1. Phase diagram of the 2D ±J RBIM. The ferromagnetic phase boundary at high temperatures (black line) starts from the pure case at TIS/J = 2.269. . . and ends at the multicritical point on the Nishimori line (red dashed line). Below the multicritical point, the phase boundary curves slightly to the left and reaches the zero-temperature ferromagnetic to spin-glass transition point at pc0 ≃ 0.9. Thus, between pN and pc0, the paramagnetic phase ‘re-enters' at sufficiently low temperatures. These phase boundaries were drawn based on data from [1–3]. |
Conventional Monte Carlo methods are inefficient in equilibrating frustrated systems at low temperatures [14, 15]. Enhanced sampling algorithms, such as the replica exchange Monte Carlo method [16–18], show better performance, but still have difficulty accessing equilibrium and true ground states. In contrast, exact matching algorithms can identify individual ground-state configurations of the 2D ±J RBIM on a planar graph [2]. Finite-size scaling (FSS) of the resulting ground-state properties yields a zero-temperature critical point at pc0 ≃ 0.897. However, since the degeneracy of the ground-state grows exponentially with the size of the system, a more complete sampling of ground state configurations is required to determine relevant statistics, especially around the FM phase boundary at zero and low temperatures.
Entropic sampling methods [19–22], which iteratively build the density of states (DOS) of the system by biasing Monte Carlo moves to flatten the energy histogram, offer a direct route to ground-state entropy. However, the original Berg's scheme [21], which has been successful in treating first-order phase transitions, converges slowly towards the ground state. Wang and Landau [23, 24] (WL) introduced a continuous iterative scheme that gained wide acceptance. An alternative algorithm to accelerate DOS convergence through extrapolation was implemented by Tang [25]. The latter has been shown to be effective in treating frustrated systems [26, 27] and is adopted in this study. Enhanced with noise filtering and DOS extrapolation, this method delivers high-precision estimates of the ground-state entropy across a wide range of p values. By careful examination of how the ground-state DOS changes as p decreases, we found compelling evidence for a mixed-order transition from the FM to the spin-glass (SG) phase. Our results also shed light on the shape of the ferromagnetic phase boundary at low temperatures.
The paper is organized as follows. In section 2 , we present the details of the entropy sampling algorithm used in our work. In section 3 , we apply the FSS analysis to the Binder cumulant and the disorder-averaged magnetization to map out the phase boundary below MCP. Section 4 presents the analyses of zero-temperature entropy and magnetization, revealing first-order–like behavior. Conclusions are drawn in section 5 . The appendix provides the simulation details.
2. Entropic sampling
The Metropolis Monte Carlo algorithm samples spin configurations {Si} according to the Boltzmann weight P(E) ∼ e−E/T, where E is the total energy of the system. Each proposed spin flip that changes the system energy from E1 to E2 is accepted with the probability 4 ). This procedure yields a flat energy histogram when running the simulation for a long enough time.
$\begin{eqnarray}P({E}_{1}\to {E}_{2})=\min \left(\frac{P({E}_{2})}{P({E}_{1})},1\right).\end{eqnarray}$
In entropic sampling, one adopts instead $\begin{eqnarray}P(E)\sim 1/g(E)={{\rm{e}}}^{-S(E)},\end{eqnarray}$
where g(E) is the density of state at energy E and $S(E)={\mathrm{ln}}\,g(E)$ is the system entropy. The spin flips are performed again following equation (As g(E) is not known a priori, Berg introduced an iterative procedure that starts with a trial sampling weight P0(E). In each round of simulation, the energy histogram hk(E) is collected over a sufficiently long period, from which the sampling weight for the next round is constructed according to
$\begin{eqnarray}{P}_{k+1}(E)=\frac{{P}_{k}(E)}{{h}_{k}(E)+a},\qquad \qquad k=0,1,\ldots \end{eqnarray}$
where a ‘floor value' a ∼ 1 is introduced for E values with insufficient statistics. The above equation can be rewritten in logarithmic form, $\begin{eqnarray}{R}_{k+1}(E)={R}_{k}(E)+{\mathrm{ln}}\,[{h}_{k}(E)+a],\end{eqnarray}$
where ${R}_{k}(E)=-{\mathrm{ln}}\,{P}_{k}(E)$ is the estimated entropy function up to a constant.Figure 2(a) shows Rk(E) obtained from simulations of the 2D ±J RBIM on a 48 × 48 square lattice at p = 0.9, starting from P0(E) = const. that corresponds to the Boltzmann distribution at T = ∞. The total number of trial spin flips increases by a factor of 1.2 in successive iterations, with measured histograms hk(E) are shown in the insert. Upon sufficient sampling, we have 7 ) becomes increasingly inefficient as one enters the low temperature regime.
$\begin{eqnarray}{h}_{k}(E)\sim g(E){P}_{k}(E)={{\rm{e}}}^{S(E)-{R}_{k}(E)}.\end{eqnarray}$
Let ${E}_{k,{\rm{\min }}}$ be the value of E where Rk(E) switches from S(E) to a constant value, as shown in figure 2(a). The condition ${h}_{k}({E}_{k,{\rm{\min }}}-{\rm{\Delta }}{E}_{k,{\rm{\min }}})=a$ yields $\begin{eqnarray}{\rm{\Delta }}{E}_{k,{\rm{\min }}}\equiv {E}_{k,{\rm{\min }}}-{E}_{k+1,{\rm{\min }}}\simeq {T}_{k}{\mathrm{ln}}\,(\bar{h}/a),\end{eqnarray}$
where ${T}_{k}^{-1}=\partial S({E}_{k,{\rm{\min }}})/\partial E$ and $\bar{h}={h}_{k}({E}_{k,{\rm{\min }}})$. Thus, the iterative scheme (
Figure 2. Comparison of two iteration schemes for the 2D ±J RBIM with N = 48 × 48 and p = 0.9. (a) Berg's iterative scheme, (b) the polynomial fitting and extrapolation scheme. Both panels show the evolution of relative entropy and histogram as a function of the number of iterations. The dashed lines indicate regions where hk = 0, while the solid lines correspond to regions where hk > 0. The polynomial fitting and extrapolation scheme exhibits significantly faster convergence compared to Berg's original method. |
To overcome this difficulty, we perform a least-squares fit of the right-hand-side of equation (7 ) to a polynomial over a range of E values around ${E}_{k,{\rm{\min }}}$ where sufficient statistics is warranted. This polynomial function is used to construct a trial Rk+1(E) up to a new ${E}_{k+1,{\rm{\min }}}\lt {E}_{k,{\rm{\min }}}$, beyond which Rk+1(E) is smoothly connected to a linear function of E at temperature Tk+1 < Tk. This procedure produces results shown in figure 2(b), with the length of simulation in each round identical to the previous case. It is evident that this modified scheme is much more efficient in reaching the low energy states of the model.
It should be noted that, as one approaches the ground state under a given bond configuration, the entropy function S(E) may develop structures not well approximated by a polynomial of low degree. Therefore, once we have reached energies quite close to the ground state, we switch back to equation (7 ) for further iteration and convergence. Furthermore, to examine in greater detail the degenerate ground states with different magnetization M = ∑iSi, we extended the flat-histogram scheme to the 2D (E, M) plane [23, 28]), building the joint DOS g(E, M) using a 2D form of equation (7 ).
The above procedure allows us to sample the entire energy range for systems up to 64 × 64. Since equation (8 ) only yields the entropy function $S(E)={\mathrm{ln}}\,g(E)$ up to a constant, we use the condition
$\begin{eqnarray}\displaystyle \sum _{E}g(E)={2}^{N}\end{eqnarray}$
to compute the absolute value of g(E) and S(E). The 2D density of states g(E, M) satisfies ∑Mg(E, M) = g(E). We have also imposed the symmetry hk(E, − M) = hk(E, M) in constructing the 2D histogram.3. Ferromagnetic phase boundary and critical exponents
In this section, we examine the critical properties of the magnetization along the phase boundary below the MCP. We apply the FSS ansatz to the high-precision data generated under the entropic sampling scheme. Despite the relatively small range of system sizes, accurate estimates of the transition line and the critical exponents are obtained, in excellent agreement with previous studies using a variety of numerical methods.
Table 1. Summary of critical point and critical exponent estimates below the MCP for the 2D ±J RBIM on the square lattice. In the table, ${L}_{{\rm{\max }}}$ denotes the largest system size examined in each study, pc is the estimated critical point, ν is the correlation length exponent, and β is the magnetization exponent. |
| Method | Year | Temperature | ${L}_{{\rm{\max }}}$ | pc | ν | β |
|---|---|---|---|---|---|---|
| Series expansion [31] | 1979 | 0 | — | ∼0.901 | — | — |
| Matching algorithm [32] | 1989 | 0 | 210 | 0.895(10) | — | — |
| Matching algorithm [33] | 1994 | 0 | 300 | 0.892∼0.905 | — | — |
| Exact Ground States [34] | 1997 | 0 | 32 | 0.896(1) | 1.30(2) | — |
| Real-space renormalization group [35] | 2001 | 0 | — | 0.8951(3) | — | — |
| 0.5 | — | 0.8919(4) | — | — | ||
| Exact Ground States [36] | 2003 | 0 | 42 | 0.8967(1) | 1.49(2) | — |
| Matching algorithm [2] | 2004 | 0 | 700 | 0.897(1) | 1.55(1) | 0.09(1) |
| Monte Carlo [3] | 2009 | 0.5 | 64 | 0.8925(1) | 1.50(4) | 0.092(2) |
| 0.645 | 64 | 0.8915(2) | 1.50(4) | 0.099(3) | ||
| Pfaffian technique [4] | 2011 | below MCP | 512 | — | 1.52(5) | — |
In entropic sampling, a single simulation is sufficient to obtain the temperature dependence of any physical quantity X for a given bond configuration. One simply keeps track of the mean value X(E) of X over configurations at energy E, and use the reweighting formula
$\begin{eqnarray}\langle X\rangle (T)=\frac{\displaystyle {\sum }_{E}X(E)g(E){{\rm{e}}}^{-E/T}}{\displaystyle {\sum }_{E}g(E){{\rm{e}}}^{-E/T}},\end{eqnarray}$
at the end of the simulation. The result is further averaged over different disorder realizations, denoted as [⟨X⟩](T), where [ ⋯ ] denotes the disorder average.In the case of total magnetization M, an alternative way is to first collect the 2D histogram g(E, M) during simulation. Its moments can then be computed according to Appendix .
$\begin{eqnarray}\langle |{M}^{n}|\rangle (T)=\frac{\displaystyle {\sum }_{E,M}|{M}^{n}|g(E,M){{\rm{e}}}^{-E/T}}{\displaystyle {\sum }_{E,M}g(E,M){{\rm{e}}}^{-E/T}},\end{eqnarray}$
followed by average over disorder realizations. The Binder cumulant $\begin{eqnarray}B(T)=\frac{[\langle {M}^{4}\rangle ]}{{[\langle {M}^{2}\rangle ]}^{2}}\end{eqnarray}$
can be calculated this way. Additional details of our implementation can be found in 
Figure 3. Binder cumulant B as a function of p below the MCP. Panels (a)–(d) correspond to temperatures T = 0, 0.5, 0.8, and 0.95, respectively. |
A more systematic approach to estimate pc(T) and the finite-size scaling exponent ν(T) is to fit the Binder cumulant data at a given T, but different sizes L and FM bond concentration p to the scaling form,
$\begin{eqnarray}B(L,p)=\tilde{B}\left((p-{p}_{c}){L}^{1/\nu }\right).\end{eqnarray}$
In [3, 29, 30], the quality of the fit is assessed using the reduced chi-square static, χ2/DOF, where a good fit corresponds to χ2/DOF ≈ 1. Here, DOF stands for the number of degrees of freedom in the fit. In our implementation, the scaling function $\tilde{B}(x)$ is approximated by a polynomial of degree ${n}_{{\rm{\max }}}=3$, $\begin{eqnarray}\tilde{B}={B}_{0}+\displaystyle \sum _{n=1}^{{n}_{\max }}{B}_{n}{x}^{n},\end{eqnarray}$
and we perform a least-squares fit to the Binder cumulant data, treating B0, {Bn}, pc, and ν as fitting parameters.Including data from very small system sizes degrades the quality of the fit, often resulting in χ2/DOF ≫ 1. As the minimum system size Lmin increases, χ2/DOF decreases and eventually approaches 1. We select Lmin as the smallest system size beyond which further increases do not significantly reduce χ2/DOF. To estimate statistical uncertainties of the fitted parameters, we employ the bootstrap method. Throughout this study, we report uncertainties at the 3σ (99.7%) confidence level, unless otherwise stated.
Table 2. Fitting results for the FSS of the Binder cumulant B at different temperatures using equation ( |
| Temperature | χ2/DOF | pc | ν | B0 |
|---|---|---|---|---|
| 0 | 47.7/40 | 0.89456(13) | 1.50(8) | 1.0755(10) |
| 0.1 | 47.7/40 | 0.89456(13) | 1.50(8) | 1.0755(10) |
| 0.2 | 47.8/40 | 0.89456(13) | 1.50(8) | 1.0755(10) |
| 0.3 | 47.7/40 | 0.89426(12) | 1.50(8) | 1.0774(9) |
| 0.4 | 29.3/40 | 0.89327(7) | 1.54(7) | 1.0831(8) |
| 0.5 | 33.5/40 | 0.89233(9) | 1.55(8) | 1.0882(9) |
| 0.6 | 36.0/40 | 0.89155(10) | 1.55(7) | 1.0937(9) |
| 0.65 | 36.4/40 | 0.89122(10) | 1.56(6) | 1.0969(9) |
| 0.7 | 36.4/40 | 0.89094(10) | 1.56(6) | 1.1006(10) |
| 0.8 | 35.7/40 | 0.89055(10) | 1.55(6) | 1.1096(10) |
| 0.9 | 35.2/40 | 0.89054(9) | 1.53(5) | 1.1207(10) |
| 0.95 | 36.2/40 | 0.89075(8) | 1.51(5) | 1.1269(9) |
| TNL | 36.3/40 | 0.89078(8) | 1.58(5) | 1.1273(10) |
With pc and ν determined above, we analyze the mean magnetization m = ⟨M⟩/N across the phase boundary. Assuming the FSS ansatz,
$\begin{eqnarray}m(L,p)={L}^{-\beta /\nu }{\bar{u}}_{h}(p)\tilde{m}\left((p-{p}_{c}){L}^{1/\nu }\right),\end{eqnarray}$
where β (not to be confused with the inverse temperature) is the magnetization exponent at a given T, ${\bar{u}}_{h}(p)$ is a function of p related to the magnetic nonlinear scaling field [3], and $\tilde{m}$ is the scaling function. According to this FSS form, we fit the numerical data using $\begin{eqnarray}\begin{array}{rcl}{\mathrm{ln}}\,m & = & -\frac{\beta }{\nu }{\mathrm{ln}}\,L+\displaystyle \sum _{n=0}^{{n}_{\max }}{a}_{n}{(p-{p}_{c})}^{n}{L}^{n/\nu }\\ & & +\displaystyle \sum _{m=0}^{{m}_{\max }}{b}_{m}{(p-{p}_{c})}^{m}.\end{array}\end{eqnarray}$
The last term represents the analytic corrections [3], and its inclusion significantly reduces χ2/DOF, bringing it closer to 1. In the fitting, we set ${n}_{\max }=3$ and ${m}_{\max }=1$, and fix pc using the value listed in table 2. The fitting results are summarized in table 3. Our estimates for β are in close agreement with the results reported in [3]: β = 0.092(2) at T = 0.5, and β = 0.099(3) at T = 0.645; as well as with those in [2]: β = 0.09(1) at T = 0.Table 3. Fitting results for the FSS of the magnetization m at different temperature, using the equation ( |
| Temperature | χ2/DOF | β/ν | β |
|---|---|---|---|
| 0 | 45.0/40 | 0.0552(7) | 0.083(5) |
| 0.1 | 45.0/40 | 0.0552(7) | 0.083(5) |
| 0.2 | 45.0/40 | 0.0552(7) | 0.083(5) |
| 0.3 | 43.8/40 | 0.0563(6) | 0.085(5) |
| 0.4 | 34.4/40 | 0.0598(5) | 0.091(5) |
| 0.5 | 40.4/40 | 0.0631(7) | 0.097(5) |
| 0.6 | 44.7/40 | 0.0666(8) | 0.103(5) |
| 0.65 | 46.4/40 | 0.0688(8) | 0.107(5) |
| 0.7 | 47.8/40 | 0.0712(8) | 0.109(5) |
| 0.8 | 49.4/40 | 0.0771(9) | 0.120(5) |
| 0.9 | 48.8/40 | 0.0843(9) | 0.130(5) |
| 0.95 | 48.0/40 | 0.0883(8) | 0.134(5) |
| TNL | 50.8/40 | 0.0886(9) | 0.141(5) |
4. Ground state and mixed-order transition
In section 3 , we analyzed the behavior of the system's magnetization near the FM phase boundary within the framework of FSS, under the implicit assumption that the transition is continuous. However, the small value of β/ν obtained could also be consistent with a first-order-like transition, where the total magnetization exhibits a discontinuous jump rather than vanishing continuously as the fraction of antiferromagnetic (AFM) bonds increases. In this section, we explore evidence supporting this alternative scenario.
4.1. Magnetization-resolved degeneracy of low-lying states
We begin by examining how the magnetization profile of the ground state and low-lying excited states evolve as AFM bonds are gradually introduced into the system. Let Na denote the number of AFM bonds, which is initially set to 0. At each step, an FM bond is randomly selected and converted into an AFM one, thereby increasing Na by one. We then compute the joint DOS g(E, M) of the system, from which the mean magnetization at a given E is obtained.
Figure 4 shows the DOS g(E, M) of the ground state and the first, second, and third excited states of a 32 × 32 system near the transition point pc0, following a particular sequence of bond flips. Spin configurations at a given E break up into two well-separated branches: the first branch on the right with only a small fraction of spins in the opposite direction, and the second branch on the left where a substantial fraction of spins are in the opposite direction. For both branches, the spacing between curves at successive En = E0 + 4n, n = 0, 1, …, is largely uniform on logarithmic scale, so that peaks at different En align at similar magnetization values.
Figure 4. Density of states g(E, M) from the ground state E0 to the third excited state E3, shown for four different values of Na in a system of size L = 32. Each subplot corresponds to a specific Na. A significant change in the ground-state magnetization distribution is observed between Na = 188 and Na = 189, where a bimodal structure emerges. Another abrupt transition occurs between Na = 193 and Na = 194, where the right peak in the ground-state distribution disappears, leaving only the left peak. Both transitions are induced by the addition of a single AFM bond. |
When the AFM bonds on the lattice are isolated from each other, the ground state is unique, with M = N and an energy Eg=−2(N − Na). For the case shown in figure 4(a), one obtains Eg = −1672 which differs from the actual E0 by only six bond energies. This indicates that the majority of AFM bonds remain unconnected from each other even as the system approaches the phase boundary. Flipping either of the two spins connected by an isolated AFM bond increases the system energy by 4. This observation leads to an estimate for the multiplicity of the nth excited state
$\begin{eqnarray}g({E}_{n},M-2)\simeq 2{N}_{a}g({E}_{n-1},M),\end{eqnarray}$
which is in agreement with the observation made above.Low energy configurations with M much smaller than N cannot be formed by flipping spins connected to isolated AFM bonds. Instead, they arise from flipping an entire cluster of spins, with minimal cost in domain-wall energy. As Na increases, this process first occurs in the excited states and gradually descends to the ground state. Figures 4(a) and (b) illustrate how this happens upon the flipping of a single bond: a subset of configurations within the left branch of the first excited state reduce their energy by 2, creating the left branch of the new ground state; on the other hand, ground-state configurations in 4(a) all increase their energy by 2, making up the right branch of the new ground state. In contrast, the situation depicted in figures 4(c) and (d) is slightly different. Upon a single bond flip, the new ground state descends from a subset of configurations in the ground state in 4(c), whereas the rest of the previous ground state configurations increase their energy by 2. As a result, the high magnetization branch of the new ground state is eliminated.
More generally, for a given spin configuration, flipping a single bond either increases or decreases the system energy by 2, depending on whether the bond was previously satisfied (i.e., with energy −1) or unsatisfied (with energy 1). Consequently, the evolution of g(E, M) can be understood as an exchange of spin configurations between neighboring energy states. Since configurations with smaller m tend to have more nearest-neighbor pairs of opposite spins, their energies are more likely to decrease upon a bond flip compared to those with larger M. As a result, for the ground state and low-lying excited states, the magnetization profile g(E, M) at a given E shifts towards smaller M as Na increases.
Figure 5(a) shows the evolution of average magnetization for low-lying states as Na increases. Up until Na < 194, these states remain FM with m(En) > 0.8. At Na = 194, an abrupt drop in ground-state magnetization (Δm ≈ 0.7) is observed, due to loss of the high magnetization branch (figure 4(d)). In contrast, the magnetization of the excited states show only a moderate decrease (Δm ≈ 0.2), without losing their high magnetization branch at this point.
Figure 5. (a): Average magnetization m as a function of the number of AFM bonds Na from the ground state (E0) to the third excited state (E3). For this sequence of bond configurations, the ground-state magnetization exhibits a first-order transition. (b): Entropy gap ΔS between the first excited state and the ground state versus Na. The dashed line, ${\mathrm{ln}}\,(2{N}_{a})$, represents the expected contribution from completely independent (single spin) excitations. A sudden increase in ΔS at Na = 194 suggests a change in the excitation mechanism at this point. |
Figure 5(b) shows the entropy gap ${\rm{\Delta }}S={S}_{1}-{S}_{0}\,={\mathrm{ln}}\,[g({E}_{1})/g({E}_{0})]$ between the first excited and ground states against Na in the neighborhood of the transition. Based on equation (18 ), we expect ${\rm{\Delta }}S\simeq {\mathrm{ln}}\,(2{N}_{a})$ at sufficiently small Na, which is borne out by data in the inset up until very close to the transition. This confirms that dominant excitations above the FM ground state are single spin flips. At Na = 194, however, ΔS increases significantly, indicating that more complex excitations emerge. The somewhat chaotic behavior of ΔS to the right is attributed to large fluctuations in g(E0) around the transition.
To investigate the microscopic origin of this abrupt change, we compare ground-state spin configurations at the largest g(E0, M) before and after the transition. Figure 6(a) shows the most probable configuration at Na = 193: a nearly uniform ferromagnet configuration with only a few isolated spins in the opposite direction. By contrast, panel (b) displays the ground state at Na = 194, which is no longer a single domain but instead splits into two oppositely magnetized regions of comparable size. The sudden emergence of multi-domain configurations underlies the observed collapse of the net magnetization and the concurrent jump in entropy. We have confirmed that the observations described above hold generally for different sequences of bond flips, although the precise value of Na to destroy the FM ground state differs from sample to sample.
Figure 6. Spin and bond configurations in the ground state for Na = 193 (a) and Na = 194 (b). Blue and yellow cubes denote spins Si = 1 and Si = −1, respectively, while red bonds indicate frustrated bonds with JijSiSj = −1. For Na = 193, both (a) and (b) represent degenerate ground states, whereas for Na = 194, only configuration (b) remains as the ground state. The primary difference between the two configurations is the presence of a large spin cluster. |
4.2. Finite-size scaling analysis of ground-state magnetization
In section 3 , we perform a standard FSS analysis of the zero-temperature Binder cumulant and magnetization, yielding critical exponents ν = 1.50(8) and β = 0.083(5). Based on these values, the data collapse of the ground-state magnetization, shown in figure 7(c), exhibits excellent quality. However, the resulting scaling behavior does not align with our anticipated physical picture.
Figure 7. Average ground-state magnetization m and its FSS under different parameters. (a) Ground-state magnetization m as a function of 1 − p. (b) and (c) FSS analysis of m using (pc0, β) = (0.898, 0) in (b) and (pc0, β) = (0.8945, 0.083) in (c). |
We propose our FSS analysis based on our physical picture. Before the transition, the system behaves exactly as a ferromagnet, a sudden decay of the ground-state magnetization is expected only after the transition. Thus we consider two scaling regions for the magnetization near the critical point19 ) or (20 ). Nevertheless, the sudden drop of ground state magnetization in a given disorder configuration as seen in figure 4 gives credence to a jump of system magnetization at the transition.
$\begin{eqnarray}m(p,L)=\left\{\begin{array}{ll}m(p),\quad & p\gt {p}_{c0},\\ \tilde{m}((p-{p}_{c0}){L}^{1/\nu }),\quad & p\lt {p}_{c0}.\end{array}\right.\end{eqnarray}$
We perform a data collapse of the segmented scaling function using pc0 = 0.898 and ν = 1.5, as shown in figures 7(a) and (b). In figure 7(c), we present a best fit of the same data using the conventional finite-size scaling form at ν = 3/2, $\begin{eqnarray}m(p,L)={L}^{-{\rm{\beta }}/{\rm{\nu }}}\tilde{m}\left((p-{p}_{c0}){L}^{1/{\nu }}\right),\end{eqnarray}$
with β = 0.083 and pc0 = 0.8945. Although the quality of data collapse is comparable to those of figure 7(a) (negative side) and figure 7(b) (positive side), the transition point is quite a bit further away from the best estimate pc0 = 0.897(1) given in [2]. Given the small range of sizes accessible in our simulations, we cannot completely rule out either equation (To further examine the scaling behavior, we fix ν = 1.5 and fit the ground-state magnetization data in the range 0.102 < 1 − p < 0.112 (as indicated by the dashed region in figure 7(a)) using the following form
$\begin{eqnarray}m=\displaystyle \sum _{n=0}^{{n}_{\max }}{a}_{n}{(p-{p}_{c0})}^{n}{L}^{n/\nu }.\end{eqnarray}$
We gradually increase the minimum system size Lmin to assess the quality of the fit. The corresponding results are summarized in table 4. As Lmin increases, χ2/DOF approaches 1, indicating that β = 0 is reasonable in the FSS analysis for larger system sizes. Based on the results in table 4, we conservatively estimate the zero-temperature critical point to be pc0 = 0.898(2). This result is very close to the large-system simulation value pc = 0.897(1) reported in [2].Table 4. Fitting results for the FSS of the magnetization m at different ${L}_{{\rm{\min }}}$, using the equation ( |
| Lmin | χ2/DOF | pc0 |
|---|---|---|
| 24 | 3345.9/69 | 0.9008(4) |
| 32 | 974.2/53 | 0.9001(4) |
| 40 | 217.6/37 | 0.8991(7) |
| 48 | 99.6/31 | 0.8988(6) |
| 56 | 16.5/13 | 0.8982(12) |
In addition to the average ground-state magnetization, we analyze the quartiles of the magnetization as a characterization of its distribution. Near the zero-temperature critical point, the quartile magnetization mq exhibits scaling behavior analogous to that of the average magnetization.
We define the lower quartile m25, the median quartile m50, and the upper quartile m75 of the magnetization based on the cumulative distribution function of the magnetization distribution P(m, p, T), which is normalized such that ∑mP(m, p, T) = 1. The cumulative distribution function is given by $Q(m,p,T)={\sum }_{| {m}^{{\prime} }| \leqslant | m| }[P({m}^{{\prime} },p,T)]$. Then, the quartiles are defined as the magnetization values satisfying
$\begin{eqnarray}Q({m}_{25},p,T)=0.25,\end{eqnarray}$
$\begin{eqnarray}Q({m}_{50},p,T)=0.50,\end{eqnarray}$
$\begin{eqnarray}Q({m}_{75},p,T)=0.75.\end{eqnarray}$
The interquartile range (IQR) of the magnetization is defined as $\begin{eqnarray}{m}_{{\rm{I}}{\rm{Q}}{\rm{R}}}={m}_{75}-{m}_{25}.\end{eqnarray}$
The same scaling analysis procedure as described in the analysis of ground-state magnetization in this section is employed, with the critical exponents fixed at β = 0 and ν = 1.5. To properly account for the FSS behavior on the SG side of the transition, only data satisfying mq < 0.7 are considered. Data with mq > 0.7 are predominantly influenced by contributions from the FM order, especially in the case of m75. As shown in figure 8, the data exhibit a good collapse in the SG phase for the lower, median, and upper quartiles, denoted m25, m50, and m75, respectively. The corresponding estimates of the critical point are pc0 = 0.8969(6), pc0 = 0.8980(7), and pc0 = 0.898(2), respectively.
Figure 8. FSS analysis of ground-state magnetization quartiles, performed with fixed critical exponents β = 0 and ν = 1.5. Panels (a)–(c) display scaling collapses for the lower quartile (m25), median (m50), and upper quartile (m75), respectively. (a) m25: Best-fit critical point pc0 = 0.8969(6), minimum system size ${L}_{{\rm{\min }}}=24$, and goodness-of-fit χ2/DOF = 81.1/60.(b) m50: pc0 = 0.8980(7), ${L}_{{\rm{\min }}}=32$, χ2/DOF = 47.9/36. (c) m75: pc0 = 0.898(2), ${L}_{{\rm{\min }}}=40$, χ2/DOF = 25.1/15. |
The mixed-order nature of the zero-temperature transition suggests the existence of a bimodal structure in the ground-state magnetization distribution in the thermodynamic limit, indicative of phase coexistence. This feature is also observed in figure 4. The upper and lower quartiles, m75 and m25, are identified as the right and left peaks of the distribution, respectively. The IQR of the magnetization, mIQR, serves as a measure of the separation between the two peaks. A finite value of mIQR in the thermodynamic limit constitutes evidence for a bimodal magnetization distribution.
We compare the behavior of mIQR at two critical points of the 2D ±J RBIM. At the pure Ising critical point, which corresponds to a well-established second-order phase transition, the IQR of magnetization exhibit excellent FSS collapse at the critical point when rescaled using the exact critical exponents, as shown in figure 9(a). The data of mIQR scales as L−1/8 and vanishes in the thermodynamic limit, indicating the absence of a bimodal distribution at the critical point. In contrast, at the zero-temperature critical point, the data collapse well with β = 0 and ν = 1.5, and a critical point pc0 = 0.893(5), as shown in figure 9(b). The critical exponent β = 0 implies that mIQR does not vanish with increasing system size, but instead converges to a finite value. This strongly suggests the presence of a bimodal structure in the ground-state magnetization distribution in the thermodynamic limit, providing direct evidence that the transition is mixed-order in nature.
Figure 9. FSS analysis of the IQR of magnetization at the Ising critical point and zero-temperature critical point. (a) FSS at the Ising critical point using TIs = 2.269, βIs = 1/8, and νIs = 1. (b) FSS at zero-temperature critical point, with best-fit pc0 = 0.893(5) and goodness-of-fit χ2/DOF = 88.5/77, using fixed exponents β = 0, ν = 1.5, and minimum system size ${L}_{{\rm{\min }}}=24$. |
5. Discussion
In this work, we investigated the 2D ±J RBIM using the entropic sampling method with an efficient iterative algorithm that we implement. Our algorithm has been successfully applied to the 2D ±J RBIM, reaching system sizes up to L = 64, allowing for detailed FSS analysis of the Binder cumulant and magnetization.
From these analyses, we determined the reentrant phase boundary below the MCP and extracted the associated critical exponents ν and β. Detailed FSS results are summarized in tables 2 and 3. The exponent ν varies between 1.50(8) and 1.56(6) below the MCP, suggesting that the zero-temperature critical point and the MCP share the same correlation length exponent. The magnetization exponent β decreases from 0.134(5) at T = 0.95 to 0.083(5) at zero temperature. The critical point pc shifts from 0.89075(8) at T = 0.95 to 0.89456(13) at T = 0.
By analyzing the magnetization in the ground state and low-energy excited states, we observed clear signatures of a mixed-order phase transition in finite-size systems, marked by a sudden drop in magnetization upon introducing a single AFM bond. Further analysis of the DOS, entropy gap and spin configurations of this sample reveals that this drop corresponds to the disappearance of the ground-state high magnetization branch and the emergence of multi-domain structures. Based on these results, we argue that the FM-SG transition at zero temperature is a mixed-order phase transition. By considering two distinct scaling regions, we find that setting β = 0 in the FSS analysis of both the average magnetization and its quartiles is justified. This result supports the interpretation that the zero-temperature transition is mixed-order in nature, in the sense that the ground-state magnetization experiences a sudden decay only after the transition in the thermodynamic limit, where a critical exponent ν = 1.5 is observed on the SG side. The presence of a bimodal magnetization structure near the zero-temperature critical point was confirmed via FSS of the IQR of magnetization. Based on the FSS analysis of the ground-state magnetization m and its quartiles m25, m50, m75, we conservatively estimate the zero-temperature critical point to be pc = 0.898(2).
Our findings at the zero-temperature critical behavior deviate slightly from previous studies [2], where exact matching algorithms were used to simulate systems up to L = 700. Although their method accurately determines the ground state, magnetization was computed using zero-temperature single-spin-flip dynamics, which may not fully capture the bimodal nature of the magnetization distribution near the critical point. This limitation could lead to inaccuracies in magnetization-related observables.
