Origin of mean-field behavior in an elastic Ising model

Simple elastic models of spin-crossover compounds are known empirically to exhibit classical critical behavior. We demonstrate how the long-range interactions responsible for this behavior arise naturally upon integrating out mechanical fluctuations of such a model. A mean-field theory applied to the resulting effective Hamiltonian quantitatively accounts for both thermodynamics and kinetics observed in computer simulations, including a barrier to magnetization reversal that grows extensively with system size. For nanocrystals, which break translational symmetry, a straightforward extension of mean-field theory yields similarly accurate results.

Figure 1

Fourier-space effective potential for the triangular lattice, Eq. (13). Note that ${\tilde{V}}_{\mathbf{q}}$ is smooth everywhere except $\mathbf{q}=0$ since ${\tilde{V}}_0=0$ but ${lim}_{\mathbf{q}\to 0}{\tilde{V}}_{\mathbf{q}}=8$.

Figure 2

Comparison between MC and MFT results for the free energy as a function of magnetization and for the critical temperature. (a) Free-energy profile for a periodic triangular lattice with $N=168$ at $T=6$, $p=0$. The MC curve was computed with umbrella sampling for the model of Eq. (2) (see [15] for details), while the MFT curve was obtained from Eq. (22). (b) MFT estimates for the triangular lattice critical temperature $T_c^{\text{MF}}=2\overline{V}\approx 7.31$ and the corresponding value of the Binder cumulant $U^{\text{MF}}\left(T_c^{\text{MF}}\right)$ [12] closely predict the intersection point of MC Binder cumulants $U$ for different system sizes. Specifically, MC indicates that $T_c^{\text{MC}}\approx 7.2$, so the MFT result is accurate within $\approx 2\%$. MC results for the Binder cumulants were computed by sampling ${10}^6$ configurations at each temperature. These configurations were generated using the effective energy function of Eq. (12) rather than Eq. (2) in order to avoid statistical errors associated with insufficient sampling of mechanical fluctuations.

Figure 3

Magnetization dynamics after a quench from temperature $T_0=8$ to $T=4$ on the triangular lattice. (a) Representative configurations from a single quench trajectory with $N=168$. The time $t$ following the quench is measured in MC sweeps. (b) Time evolution of the $q=\left|\mathbf{q}\right|=0$ (longest-wavelength) Fourier component of the spin structure factor $\mathcal{M}(\mathbf{q},t)$ following the quench. This mode grows rapidly at short times and saturates at the equilibrium value of $m^2$. Solving the mean-field master equation, Eq. (26), for a system size $N=168$ and initial condition $P(m,0)=exp[-F_{\text{MF}}(m;T=T_0)/k_B{T}_0]$ yields a prediction $\mathcal{M}(0,t)={\sum }_{m=-1}^1{m}^{2}P(m,t)$ (labeled MFT in the plot) which closely agrees with the MC result. (c) Time evolution of a short-wavelength Fourier component of $\mathcal{M}$ with $q=4\pi /3$ (corresponding to a corner of the first Brillouin zone of the triangular lattice) computed from MC simulations. This mode decays rapidly, consistent with the apparent lack of short-wavelength structure in the configurations. MC curves in both (b) and (c) were obtained by averaging over ${10}^3$ independent trajectories initialized from equilibrium configurations sampled at $T_0=8$ and propagated with Metropolis spin-flip dynamics at $T=4$. All MC simulations here were performed using ${\mathcal{H}}_{\text{eff}}$, Eq. (12).

Figure 4

Hysteresis loop at $T=6$ from MFT and MC. The MFT curve was obtained by numerical solution of Eq. (21). Solutions to this equation which are also local free-energy minima (of which there are at least one and at most two) comprise the mean-field hysteresis loop. MC Dyn. (dynamic) results were obtained by sweeping the field from $h=-1.0$ to $h=0.9$ and back again (direction indicated by the black arrows) for a simulated system with $N=2688$ using ${\mathcal{H}}_{\text{eff}}$ [Eq. (12)]. For each field value, there were ten MC sweeps of equilibration and ten MC sweeps of production. MC Eq. (equilibrium) results were obtained by locating the local minima of the free energy as a function of magnetization [computed with umbrella sampling using $\mathcal{H}$, Eq. (2)] for different values of $h$ and a system size of $N=168$. Inset schematics illustrate the fact that, in the thermodynamic limit, barrier crossing does not occur; the system can escape from a metastable well only once it has reached the limits of metastability. The yellow star indicates the values of $m$ and $h$ used as a starting point for dynamics in Fig. 5.

Figure 5

Average magnetization versus time for phase change dynamics at the end of the hysteresis loop (left) with $T=6$ and the corresponding free-energy profile (right) with $T=6$, $h=0.5$. Gray dashed lines indicate the positions of the formerly metastable well and the single stable well. The yellow star indicates the initial state, marked with the same symbol in Fig. 4. The MFT result for $\langle m\left(t\right)\rangle$ was calculated via Eq. (34) for a system size $N=168$, and the mean-field free energy is given by Eq. (22). The MC result for $\langle m\left(t\right)\rangle$ was computed by averaging over ${10}^4$ independent trajectories of a system with $N=168$. These trajectories were propagated by Metropolis MC according to the effective energy ${\mathcal{H}}_{\text{eff}}$ [Eq. (12)]. Their initial configurations were sampled from an equilibrium trajectory whose magnetization $m=-0.7$ was fixed by performing Kawasaki dynamics [41]. The MC result for the free energy was computed via umbrella sampling of a system with $N=168$ using $\mathcal{H}$ [Eq. (2)].

Figure 6

Probability distribution of the first-passage time for phase change at the end of the hysteresis loop with $T=6$, plotted on (a) linear and (b) logarithmic scales. For each MC trajectory, the first-passage time was defined as the number of MC steps taken en route from the formerly metastable initial state ($m=-0.7$) to the bottom of the stable well ($m=0.76$). Results are shown for a system of size $N=168$. MFT predictions were computed via Eq. (35) for the same system size. The long-time exponential tail is characteristic of diffusion on a bounded interval; its corresponding decay rate is set by the least negative eigenvalue of $\mathit{\Omega }$ [42]. The entire eigenvalue spectrum of $\mathit{\Omega }$ is sensitive to changes in $N$, highlighting the size dependence of first-passage time statistics that appears to be well captured by MFT.

Figure 7

Pair interaction function $V_{\mathbf{R},{\mathbf{R}}^{\prime }}$ for different locations $\mathbf{R}$ of a tagged atom (outlined in black). The value of $V_{\mathbf{R},{\mathbf{R}}^{\prime }}$ for interaction with another atom at ${\mathbf{R}}^{\prime }$ is indicated by color according to the scale shown.

Figure 8

Numerical solutions to Eq. (38) for the position-dependent mean-field magnetization of a hexagonal nanocrystal of size $N=271$ with triangular lattice structure. Curves with different shades of red represent the magnetization of different sites in the nanocrystal. Sites near the perimeter of the crystal have smaller magnetization than sites well within the interior; all sites transition from zero to nonzero magnetization at a temperature $T_c\approx 6.2$. The vertical dashed line marks the bulk value for $T_c$; open boundary conditions thus suppress the nanocrystal $T_c$ compared to its bulk value.

Figure 9

Nanocrystal magnetization as a function of temperature for different system sizes. Bulk magnetization versus temperature (for a system with $N=168$ subject to periodic boundary conditions) is included for comparison. Nanocrystal MFT curves were obtained as the numerical solutions of Eqs. (38) and (39). Nanocrystal MC points were obtained as the minima of free-energy profiles computed via umbrella sampling of the effective Hamiltonian, Eq. (37), for each system size. The bulk MFT curve was computed using Eq. (21), and bulk MC points were obtained from bulk free-energy minima computed via umbrella sampling using Eq. (2).

Figure 10

Nanocrystal free energies per atom for different system sizes at $T=3$, $h=0$. Curves were computed via umbrella sampling of the effective Hamiltonian, Eq. (37), for each system size. These profiles strongly suggest a free-energy barrier which grows linearly with $N$.