Phase diagrams and free-energy landscapes for model spin-crossover materials with antiferromagnetic-like nearest-neighbor and ferromagnetic-like long-range interactions

TL;DR

This paper uses advanced Monte Carlo simulations to explore phase diagrams and free-energy landscapes of a spin-crossover model with competing interactions, revealing complex critical phenomena and hysteresis behaviors relevant to experimental materials.

Contribution

It introduces an efficient simulation method to analyze thermodynamics of a spin-crossover model with both short- and long-range interactions, uncovering new critical phenomena.

Findings

Identification of tricritical points and critical end points in phase diagrams

Observation of horn-shaped metastability regions

Insights into hysteresis loops in spin-crossover materials

Abstract

We present phase diagrams, free-energy landscapes, and order-parameter distributions for a model spin-crossover material with a two-step transition between the high-spin and low-spin states (a square-lattice Ising model with antiferromagnetic-like nearest-neighbor and ferromagnetic-like long-range interactions) [P. A. Rikvold et al., Phys. Rev. B 93, 064109 (2016)]. The results are obtained by a recently introduced, macroscopically constrained Wang-Landau Monte Carlo simulation method [C. H. Chan, G. Brown, and P. A. Rikvold, Phys. Rev. E 95, 053302 (2017)]. The method's computational efficiency enables calculation of thermodynamic quantities for a wide range of temperatures, applied fields, and long-range interaction strengths. For long-range interactions of intermediate strength, tricritical points in the phase diagrams are replaced by pairs of critical end points and mean-field critical points that give rise to horn-shaped regions of metastability. The corresponding free-energy landscapes offer insights into the nature of asymmetric, multiple hysteresis loops that have been experimentally observed in spin-crossover materials characterized by competing short-range interactions and long-range elastic interactions.

Equations

$$\mathcal{H}=J\sum _{⟨i,j⟩}{s}_i{s}_j-HM-\frac{A}{2L^2}{M}^2,$$

$$H=\frac{1}{2}(k_{\mathrm{B}}T\ln r-D),$$

$$m$$

$$m_s$$

$$g(E(H,A),M,M_s)=g(E(0,0),M,M_s)$$

$$E(H,A)=E(0,0)-HM-\frac{A{M}^2}{2L^2}.$$

$$Z_{m,m_s}=\sum _{E}g(E,m,m_s)e^{-E/T}.$$

$$Z_{\mathrm{all}}=\sum _{m,m_s}{Z}_{m,m_s}.$$

$$P(m,m_s)\Delta m\Delta m_s=\frac{Z_{m,m_s}}{Z_{\mathrm{all}}},$$

$$F(m,m_s)=-T\ln Z_{m,m_s}.$$

$$P(m)\Delta m$$

$$P(m_s)\Delta m_s$$

$$F(m)$$

$$F(m_s)$$

$$⟨m⟩$$

$$⟨m_s⟩$$

$$u_{m_s}=1-\frac{⟨m_s^4⟩}{3⟨m_s^2{⟩}^2}.$$

$$\text{var}(m)={\chi }_{m}T=L^2(⟨m^2⟩-⟨m{⟩}^2),$$

$$(\Delta H{)}^{2/3}\propto T_x-T,$$

Full text

I introduction

Spin-crossover (SC) materials are molecular crystals in which the individual molecules contain transition metal ions that can exist in two different spin states: a low-spin ground state (LS) and a high-spin excited state (HS). Molecules in the HS state have larger volume and higher effective degeneracy than those in the LS state Miyashita et al. (2008); Nakada et al. (2011, 2012); Enachescu et al. (2005); Nishino et al. (2007); Konishi et al. (2008); Bousseksou et al. (2011); Halcrow (editor). Due to its higher degeneracy, crystals of such molecules can be brought into a majority excited HS state by increasing temperature, changing pressure or magnetic field, electrochemical stimuli, or exposure to light Miyashita et al. (2009); Konishi et al. (2008); Nishino et al. (2010); Enachescu et al. (2005); Gütlich et al. (1994); Chong et al. (2010); Ohkoshi et al. (2011); Asahara et al. (2012); Chakraborty et al. (2014); Miyamoto et al. (2016); Mathonie et al. (2017). The size difference between the HS and LS molecules causes local elastic distortions that lead to effective long-range elastic interactions mediated by the macroscopic strain field Nishino et al. (2007); Teodosiu (1982); Zhu et al. (2005). In addition to such long-range interactions, the materials will also typically have local interactions caused by, e.g., quantum-mechanical exchange or geometric restrictions. These intermolecular interactions may cause first-order phase transitions that can render the HS state metastable and lead to hysteresis when exposed to time-varying fields Gütlich et al. (1994); Miyashita et al. (2005). In the case of optical excitation into the metastable phase, this phenomenon is known as light-induced excited spin trapping (LIESST) Gütlich et al. (1994); Miyashita et al. (2009); Mathonie et al. (2017). The metastable properties in combination with the SC materials’ sensitivity to a wide range of external stimuli make them promising candidates for applications such as switches, displays, memory devices, sensors, and actuators Ohkoshi et al. (2011); Chakraborty et al. (2014); Kahn and Jay Martinez (1998); Linares et al. (2012).

In the SC literature, the phase transitions caused by the short-range and long-range interactions are often discussed using an Ising-like pseudospin formulation, in which the HS state is represented as $s=+1$ and the LS state as $s=-1$. This is the representation we will use in this paper. It has the advantage of a high degree of symmetry, and it enables easy reference to studies of other Ising-like models. To minimize the strain energy, the elastic long-range interaction favors different molecules being in the same state (LS-LS or HS-HS). In this pseudospin language it is therefore called ferromagnetic-like, or simply ferromagnetic. The short-range interactions depend on the particular material and may either be ferromagnetic-like, or they may favor neighboring molecules in opposite states (LS-HS), which is analogously called antiferromagnetic-like, or simply antiferromagnetic. We emphasize that this nomenclature only represents an analogy and does not imply a magnetic origin of the interactions. In the remainder of this paper, we will use the simplified terms, ferromagnetic and antiferromagnetic, for interactions that favor uniform and checkerboard spin-state arrangements, respectively.

If the short-range interaction is ferromagnetic, it has been found that adding even a very weak long-range interaction causes the universality class of the critical point to change from the Ising class to the mean-field class Nakada et al. (2011, 2012). On the other hand, if the short-range interaction is antiferromagnetic, the critical line will terminate at a certain point, with the appearance of metastable regions in the phase diagram, bounded by sharp spinodal lines Brown et al. (2014); Rikvold et al. (2016). Then, with sufficiently strong long-range interaction, new mean-field critical points emerge in the phase diagrams – a phenomenon which is not predicted by simple Bragg-Williams mean-field theory Rikvold et al. (2016). These new mean-field critical points also become the end points for the spinodal lines bounding the metastable regions.

In some SC materials, the transition between the LS and HS phases proceeds as a two-step transition via an intermediate phase Köppen et al. (1982); Zelentsov et al. (1985); Petrouleas and Tuchagues (1987); Bousseksou et al. (1992); Jakobi et al. (1992); REAL92 ; Boinnard et al. (1994); Bolvin (1996); Chernyshov et al. (2004); Bonnet et al. (2008); Pillet et al. (2012); Buron-Le Cointe et al. (2010); Lin et al. (2012); Klein et al. (2014), giving rise to complex, asymmetrical hysteresis loops. In the case of Fe(II)[2-picolylamine]3Cl${}_{2}\cdot$Ethanol Köppen et al. (1982), x-ray diffraction has revealed an intermediate phase, characterized by long-range order on two interpenetrating sublattices with nearest-neighbor molecules in different states (HS-LS) Chernyshov et al. (2003); Huby et al. (2004). Several of these experimental results were recently reviewed Shatruk et al. (2015); Brooker (2015). This situation can be modeled by an Ising-like model with antiferromagnetic nearest-neighbor interactions. Various mean-field approximations to this model have been considered, both without Bolvin (1996) and with Zelentsov et al. (1985); Bousseksou et al. (1992); Chernyshov et al. (2004) a long-range ferromagnetic term.

Recently, Rikvold et al. used standard importance-sampling Monte Carlo (MC) simulations to obtain phase diagrams and hysteresis curves for such an Ising model with nearest-neighbor antiferromagnetic interactions and ferromagnetic long-range interactions approximated by a mean-field equivalent-neighbor (Husimi-Temperley) term Rikvold et al. (2016). (See Hamiltonian in Sec. II.) To locate the various transition lines in the phase diagram, this method requires separate simulations for different values of temperature, field, and long-range interaction strength. This procedure is very computationally intensive, and phase diagrams could therefore only be drawn for three different interaction strengths.

In the present paper, we provide detailed phase diagrams for this system with a range of different long-range interaction strengths, from quite weak to quite strong. In addition to phase diagrams, we also obtain free-energy landscapes and order-parameter probability densities in terms of the model’s two order parameters, magnetization ($M$) and staggered magnetization ($M_{s}$). To obtain these results with a reasonably modest computational effort, we use a recently proposed, macroscopically constrained Wang-Landau (WL) MC algorithm Ren et al. (2016); Chan et al. (2017). With this method, a simple analytic transformation of the system energy $E$ enables us to extract results for any combination of temperature, applied field, and long-range interaction strength from one single, high-precision simulation of the joint density of states (DOS), $g(E,M,M_{s})$, for a simple square-lattice Ising antiferromagnet in zero field. The details of how to use the algorithm to calculate the joint DOS, and how to extract from it free-energy landscapes, order-parameter probability densities, and phase diagrams are given in our recent papers, Ref. Chan et al. (2017, 2017). Here, we concentrate on the physical aspects of this model SC material and, in particular, their dependence on the long-range interaction strength. In the process, we also obtain improved estimates for the positions and shapes of the first-order coexistence lines in the phase diagrams.

Studies of Ising models with long-range interactions have a long history. Some notable examples are work on Ising models with weak long-range interactions by Penrose, Lebowitz, and Hemmer Penrose and Lebowitz (1971); Hemmer and Lebowitz (1976); Penrose and Lebowitz (1979), and with long-range lattice coupling by Oitmaa and Barber Oitmaa and Barber (1975). Herrero studied small-world networks with both ferromagnetic Herrero (2002) and antiferromagnetic interactions Herrero (2008). Hasnaoui and Piekarewicz Hasnaoui and Piekarewicz (2013) recently used an Ising model with Coulomb long-range interaction to simulate nuclear pasta in neutron stars. It should also be mentioned that the Ising model with long-range interactions decaying as $r^{-(d+\sigma)}$ with $d=1,2,3$ and $0<\sigma<d/2$ was studied by Luijten and Blöte Luijten and Blöte (1997), and the effect of long-range interactions on phase transitions in short-range interacting systems were studied by Capel et al. Capel et al. (1979).

The remainder of this paper is organized as follows. In Sec. II we present the Ising-like model Hamiltonian and its interpretation as a model for SC materials. In Sec. III we briefly discuss the macroscopically constrained WL algorithm and present the analytic energy transformation that enables us to extract data for arbitrary model parameters from a single simulated joint DOS. We also show how constrained partition functions are obtained from the joint densities of states, and how the partition functions lead to free-energy landscapes and order-parameter probability densities. Sec. IV contains our main results: phase diagrams, as well as probability densities and free-energy landscapes at selected phase points. All these are obtained for several values of the long-range interaction strength, ranging from quite weak to quite strong, and producing a number of topologically different phase diagrams. Section V contains a brief summary and conclusions. Details of our estimates of finite-size and statistical errors are given in Appendix A.

II 2D Ising-ASFL model

To approximate a SC material with antiferromagnetic-like nearest-neighbor interactions and ferromagnetic-like elastic long-range interactions, we here employ the model introduced by S. Miyashita and first used in Refs. Brown et al. (2014); Rikvold et al. (2016). This is a $L\times L$ square-lattice nearest-neighbor Ising antiferromagnet with ferromagnetic equivalent-neighbor (aka Husimi-Temperley) interactions. It is defined by the Hamiltonian,

[TABLE]

with $J>0$. We name it the two dimensional Ising Antiferromagnetic Short-range and Ferromagnetic Long-range (2D Ising-ASFL) model. The first two terms constitute the Wajnflasz-Pick Ising-like model Wajnflasz and Pick (1971), in which the pseudo-spin variable $s_{i}$ denotes the two spin states at site $i$ ($-1$ for LS and $+1$ for HS), and $M=\sum_{i}s_{i}$ is the pseudomagnetization. The effective field term,

[TABLE]

contains $D>0$, which is the energy difference between the HS and LS states, and $r$, which is the ratio between the HS and LS degeneracies. $T$ is the absolute temperature, and $k_{\rm B}$ is Boltzmann’s constant. (Changing the temperature in the physical SC system therefore corresponds to a combined change in temperature and effective field in this pseudospin model. See Figs. 5(a) and 8 of Ref. Rikvold et al. (2016).)

The last term in Eq. (1) approximates the elastic long-range interactions of the SC material as in Refs. Mori et al. (2010); Nakada et al. (2011); Rikvold et al. (2016). Since it lowers the energy of more uniform spin-state configurations (mostly +1 or mostly $-1$) in a quadratic fashion, it is a ferromagnetic term. Throughout the paper, temperature ($T$), energy ($E$), magnetic field ($H$), and long-range interaction strength ($A$), will be expressed in dimensionless units ($|J|=k_{B}=1$).

The order parameters of this model are magnetization ($M$) and staggered magnetization ($M_{s}$). They can be normalized as $m=M/L^{2}$ and $m_{s}=M_{s}/L^{2}$. If we break the two-dimensional square lattice into two sublattices (A and B), like the black and white squares on a chessboard, $m$ and $m_{s}$ can be expressed in terms of the normalized magnetizations ($m_{A}$, $m_{B}$) of these two sublattices as

[TABLE]

The usual order parameter for SC materials is the proportion of HS molecules, $n_{\rm HS}$, which is related to the pseudospin variables as $n_{\rm HS}=\left(m+1\right)/2$.

The equilibrium (stable) and metastable phases at zero temperature were obtained from the Hamiltonian by simple ground-state calculations in Rikvold et al. (2016). We briefly repeat the results here for convenient reference, also introducing the following short-hand notation for the low-temperature ordered phases:

antiferromagnetic (which is doubly degenerate), is called AFM;

ferromagnetic with majority of $s_{i}=+1$, is called FM$+$;

and ferromagnetic with majority of $s_{i}=-1$, is called FM$-$.

$A<8$: AFM is stable for $-4+A/2<H<4-A/2$, metastable against transition to FM+ for $4-A/2<H<4$, and metastable against transition to FM$-$ for $-4<H<-4+A/2$. FM+ is stable for $H>4-A/2$, and metastable for transition to AFM or FM$-$ for $4-A<H<4-A/2$. FM$-$ is stable for $H<-4+A/2$, and metastable for transition to AFM or FM$+$ for $-4+A/2<H<-4+A$.

$A>8$: AFM is never the stable ground state, but it is metastable for $-4<H<4$. FM+ is stable for $H>0$ and metastable for $4-A<H<0$. FM$-$ is stable for $H<0$ and metastable for $0<H<-4+A$.

III Method

III.1 Obtaining joint density of states

The results presented in this paper are all based on the joint DOS, $g(E,M,M_{s})$, determined once for $H=A=0$, which corresponds to a simple square-lattice Ising antiferromagnet. Using this, the joint DOS for any arbitrary value of $(H,A)$ can be obtained by

[TABLE]

where

[TABLE]

Note that this is an alternative, but equivalent way to express the content of Eq. (10) in Ref. Chan et al. (2017). This result is based on the fact that all the microstates are equally shifted in energy when a field-like parameter couples to a function of the global property $M$, as shown in Eq. (1). With the joint DOS, all thermodynamic quantities can be calculated, as demonstrated in Chan et al. (2017). From $g(E,M,M_{s})$ at different $(H,A)$, we can obtain $g(E,M)$ and $g(E,M_{s})$, as shown in Ref. Chan et al. (2017).

To obtain an accurate $g(E,M,M_{s})$ at $H=A=0$, the macroscopically constrained WL method is used Chan et al. (2017); Ren et al. (2016). With the help of simple combinatorial calculations in the $(M,M_{s})$ space, the method converts what would otherwise be a time-consuming multi-dimensional random walk in the $(E,M,M_{s})$ space into many independent, one-dimensional random walks in $E$, each constrained to a fixed value of $(M,M_{s})$. Through further, symmetry based simplifications Chan et al. (2017), the method can obtain an accurate estimate of $g(E,M,M_{s})$ in a relatively short time.

As the details of how to arrive at these results have already been presented in Chan et al. (2017), here we simply focus on the physics of the model SC material as $A$ is changed. All the phase diagrams, free-energy landscapes, and probability densities shown in Sec. IV are obtained with $L=32$.

III.2 From joint density of states to thermodynamic quantities

We define the constrained partition function of any macrostate $(m,m_{s})$ as

[TABLE]

The overall partition function of the system is then

[TABLE]

The joint probability of finding the system in a macrostate $(m,m_{s})$ is

[TABLE]

where $\Delta m$, $\Delta m_{s}$ are the order-parameter step sizes, both chosen to be the same value, around $0.03$. The free energy of macrostate $(m,m_{s})$ is

[TABLE]

We will plot these quantities in terms of $(m_{A},m_{B})$ which have a one-to-one relation with $(m,m_{s})$ (see Eqs. (3) and (4)).

Summing over the contributions of the joint probability (Eq. (9)) in one direction, we obtain the marginal probability densities as

[TABLE]

With these densities, we can calculate the expectation values of the order parameters and other quantities. We can express the free energy in terms of one order parameter as

[TABLE]

The presence of the long-range interaction induces metastable phase regions in the phase diagrams. A very important point is that when we consider values of $(T,H,A)$ lying in those regions, the stable phase will be the phase that has larger total area in the marginal probability density, rather than the phase that shows the higher peak. Systems lying on the coexistence line between two phases will have equal areas in the marginal probability density.

In a free-energy contour plot or joint probability density plot, against $m$ and $m_{s}$ (or against $m_{A}$ and $m_{B}$), the region around $(m,m_{s})=(1,0)$ [or $(m_{A},m_{B})=(1,1)$] corresponds to the FM+ phase. Similarly, the region around $(m,m_{s})=(-1,0)$ [or $(m_{A},m_{B})=(-1,-1)$] corresponds to the FM$-$ phase. The region around $(m,m_{s})=(0,1)$ [or $(m_{A},m_{B})=(1,-1)$] corresponds to the AFM+ phase, and the region around $(m,m_{s})=(0,-1)$ [or $(m_{A},m_{B})=(-1,1)$] corresponds to the AFM$-$ phase. Finally, the region around $(m,m_{s})=(0,0)$ [or $(m_{A},m_{B})=(0,0)$] corresponds to the disordered phase. However, these are just the most extreme cases. Some AFM phases have significant ferromagnetic properties, and some FM phases may be quite disordered.

In our model, for a particular $(T,H,A)$ triple, if the system can exist as a disordered phase, it cannot exist as an AFM phase, and vice versa. However it may happen that a disordered phase shows strong AFM properties. Changing $(T,H,A)$ may let the system change from one phase to another through a continuous phase transition, as it crosses the critical line between the two phases. In the Ising-ASFL model, a critical line only exists between the disordered phase and the AFM phase. The phase boundary between the ferromagnetic phase and the disordered phase is a coexistence line, and it ends with a mean-field critical point for sufficiently strong long-range interaction $A$. This critical point is located where the two spinodal lines meet.

The expectation values of the two order parameters can be obtained easily as

[TABLE]

As the two AFM phases always exist in pairs and the probability of finding the system in both are the same, $\langle m_{s}\rangle=0$.

As $\langle m_{s}\rangle=0$, we define the corresponding fourth-order Binder cumulant as Landau and Binder (2015); Binder (1981); Binder and Landau (1984); Challa et al. (1986),

[TABLE]

Here we only define the cumulant for the order parameter $m_{s}$, as only the critical line will be located by the cumulant. When we take the ensemble average, we have to exclude all the phase points that belong to the metastable FM+ or FM$-$ phase. That is, when we look at $F(m)$, if we find more than one minimum (i.e. more than one phase are found), we neglect the states that have values of $|m|$ greater than the separating value of $m$. The critical line in this model is commonly accepted to be in the Ising universality class Lourenço and Dickman (2016), which (assuming isotropy, periodic boundary conditions, and a square shape as in the present study) has a cumulant fixed-point value of $0.6106924(16)$ Kamieniarz and Blöte (1993); Chen and Dohm (2004); Selke and Shchur (2005); Salas and Sokal (2000). We therefore locate the critical line by finding the phase point within a temperature range where the cumulant is close to $0.61$, and does not deviate from $0.61$ by more than $0.01$. The resulting critical line for $A=0$ is included in Fig. 1 together with the analytically approximated critical line for the pure square-lattice Ising antiferromagnet in the thermodynamic limit from Ref. Wu and Wu (1990). Within the resolution of this figure, our $L=32$ data coincide with this highly accurate approximation.

The variance of the order-parameter $m$, which is proportional to the susceptibility times the temperature,

[TABLE]

is considered as we use its maxima to separate the FM$\pm$ phases from the disordered and AFM phases. All the coexistence lines that we show are located by using this quantity. Note that this quantity is very difficult to measure through importance-sampling MC, while our approach can directly calculate it using $g(E,M,M_{s})$. Further details on the method are given in Ref. Chan et al. (2017).

In next section, we consider the phase diagrams for different values of $A$ and study selected phase points. These are the main results of the present paper. Notice that all the phase diagrams are symmetric about the $T$ axis, with an exchange between FM$+$ and FM$-$. For $A=0$, the model reduces to the standard square-lattice antiferromagnetic Ising model Lourenço and Dickman (2016); Chan et al. (2017).

IV Phase diagrams

IV.1 Weak long-range interaction, $A=1,4$

It is reasonable to assume that adding a ferromagnetic long-range interaction $A$ to the pure antiferromagnet must favor the appearance of the ferromagnetic phases, and thus push the critical line towards lower values of $|H|$. Figure 1 supports this assumption. Moreover, the critical lines also terminate at lower $|H|$ and higher $T$ for larger $A$. The phase diagrams in Fig. 2 show that the critical lines end with the appearance of a metastable region in the phase diagram, and that the metastable region grows as $A$ increases. All phase diagrams shown in this paper are symmetric under simultaneous reversal of $H$ and $m$. Error bars including statistical and finite-size errors are included with every data point in this and all subsequent phase diagrams. With the exception of Fig. 6, they are everywhere smaller than the symbol size. A discussion of how the errors were estimated is found in Appendix A.

Introducing the long-range interaction $A$ with the $M^{2}$ term makes it much weaker than the $HM$ term for small $M$, so that the long-range interaction effect is negligible when $H$ and $A$ are small, and so it does not significantly affect the critical temperature near $H=0$. On the other hand, when we increase $H$, the $M^{2}$ term will eventually be larger than the $M$ term, and finally causes a local free-energy minimum to show up in the FM$+$ region, corresponding to a metastable FM$+$ phase region in the phase diagram (Figs. 3 (a) and (b)). A new FM$+$ peak also appears in the joint probability density ($P(m_{A},m_{B})$) and marginal probability densities ($P(m)$ and $P(m_{s})$). One peak may be much smaller than the other, such that it may not be easy to discover the presence of metastability through looking at the probability density (Figs. 3 (b) and (d)). Notice that although one phase may have much smaller probability density than the other, the lifetimes for these metastable phases increase exponentially with system volume, $e^{cL^{2}}$ for a two-dimensional system, so that they are still macroscopic, and thus cannot be neglected Mori et al. (2010); Rikvold and Gorman (1994); Binder and Virnau (2016).

The AFM and FM+ phases are separated by the coexistence line in the metastable region, and we observe that when $T$ is low, the coexistence line is a practically straight line at constant $H$ in the phase diagram. Note that this result is different from Rikvold et al.’s former result Rikvold et al. (2016) for $A=4$, which indicates a reentrant behavior of the coexistence line at low $T$. This discrepancy is probably due to incomplete ergodicity in the importance-sampling MC with mixed initial conditions used in Ref. Rikvold et al. (2016).

For any point lying on that straight vertical segment of the coexistence line, as in Fig. 4 (a)-(c), the coexisting AFM phases and the FM+ phase are located at their extreme locations, i.e., $m=+1,m_{s}=\pm 1$. Increasing $T$ bends the coexistence line toward lower $|H|$ values. Simultaneously, the AFM phases and the FM$+$ phase move away from the extreme positions and towards each other, as shown in Fig. 4 (d)-(f). The coexistence line finally joins the critical line at the tricritical point, where the two AFM phases and the FM+ phase become indistinguishable at the continuous phase transition point. Figure  4 (g)-(i) represent a point lying on the coexistence line, below the tricritical point. We see from the joint probability density in Fig. 4 (g) that the ferromagnetic phase and the AFM phases are coalescing. However, the marginal probability along the $m$ axis in Fig. 4 (h) still has two peaks. We therefore regard the system as in AFM/FM$+$ coexistence, with this small system fluctuating easily between the two phases. Extrapolation of the end points of the two spinodal lines gives the merging temperature, which corresponds to the tricritical point. When the two spinodal lines merge, the distance between them ($\Delta H$) varies against temperature as Newman and Schulman (1980)

[TABLE]

where $T_{x}$ represents the tricritical or critical temperature, where the coexistence line ends. After obtaining the tricritical temperature, we can estimate the tricritical field as the average of the extrapolation points of the two spinodal lines. Figure 4 (j)-(l) show data at the tricritical point for $A=4$, where the AFM phases and the FM$+$ peak finally join together into one phase.

IV.2 Medium long-range interaction, $A=6,7,8$

As mentioned above for small $A$, moving along the coexistence line toward the critical line, one approaches a tricritical point, where the two AFM phases and the FM phase become indistinguishable. Below the tricritical temperature, the three phases are distinct. Then it is reasonable to expect that, if $A$ is big enough, the two AFM phases may combine into one disordered phase at a lower $T$ than the one where they further combine with the FM phase. In this scenario, we will find that the critical line, which represents the AFM/disordered phase transition, intersects the coexistence line at a critical end-point, and new metastable regions (horn regions) emerge in the phase diagram as shown for $A=6,7$, and $8$ in Fig. 5.

Figurte 6 is a closer look at the horn region for $A=7$. The coexistence line separates the FM phase from the AFM phases at low $T$. After passing through the critical end-point, it separates the FM phase from the disordered phase. At a higher $T$, it ends in a mean-field critical point, where the disordered and FM phases become indistinguishable.

Figure 7 shows the case near the critical end-point. As this point is the intersection of the critical line and the coexistence line, it has properties of both lines. Since it is on the coexistence line, the combined AFM/disordered phase is equally probable as the FM$+$ phase, as shown in (c) and (e). Since it is on the critical line, the AFM peaks are connected through the middle disordered region as it corresponds to a continuous phase transition between the AFM phases and the disordered phase (shown in (b)). For the marginal probability density function $P(m_{s})$, if we remove the contribution from the FM$+$ phase as shown in the inset, the height ratio between a AFM peak to the central point in the middle between the two peaks is around $26/1$, which is close to the established value of about $22/1$ Nicolaides and Bruce (1988). Figure 8 shows a point close to the mean-field critical point at $(H,T)=(0.566,2.63)$ for $A=7$, where we see that the two peaks in $P(m)$ have coalesced into one single peak. We note that the position of the critical point found here is consistent with the one found in Ref. Rikvold et al. (2016) by importance-sampling MC with system sizes up to $L=1024$, $H=0.561(1)$ and $T=2.61(1)$.

Figure 9 shows results as we move along the coexistence line to a point near the mean-field critical point. The disordered phase peak gradually contracts to $m_{s}=0$ as the AFM fluctuations weaken (refer to the first row of the figure), and the FM$+$ peak slowly merges with the disordered phase peak until only one peak is left in the marginal probability along the $m$ direction (refer to the second row of the figure). We see that the two peaks in the marginal probability density, $P(m)$, along the FM axis, which correspond to two different phases, become less sharp and merge. Note that the joint probability density in (g) seems to show only one peak, but after summing up all the contributions from different $m_{s}$, the marginal probability density in (h) shows two peaks, and we still regard them as two phases even though they are strongly connected by fluctuations.

Figure 10 shows the results observed at four points that are equidistant from the coexistence line, but lie in four different phase regions, with the critical end-point nearly at the center, as shown by the four red dots in Fig. 6. Parts (a)-(b) show a point lying in region III, which has stable disordered phase and metastable FM$+$ phase; (c)-(d) show a point lying in region IV, which has stable AFM phase and metastable FM$+$ phase; (e)-(f) show a point lying in region V, which has stable FM$+$ phase and metastable disordered phase; and (g)-(h) show a point lying in region VI, which has stable FM$+$ phase and metastable AFM phase.

The phase diagram for $A=7$ is well suited for comparison with a number of experimental results for SC materials that show asymmetric, two-step thermal hysteresis loops Köppen et al. (1982); Zelentsov et al. (1985); Petrouleas and Tuchagues (1987); Bousseksou et al. (1992); Jakobi et al. (1992); REAL92 ; Boinnard et al. (1994); Bolvin (1996); Chernyshov et al. (2004); Bonnet et al. (2008); Pillet et al. (2012); Buron-Le Cointe et al. (2010); Lin et al. (2012); Klein et al. (2014); Chernyshov et al. (2003); Huby et al. (2004); Shatruk et al. (2015); Brooker (2015). Such a two-step loop, obtained directly from the joint probability density, $P(m,m_{s})$, along a path between $(H,T)=(1.0,2.5)$ and $(-1.5,1.75)$ in Fig. 6 is shown in Fig. 11. This path corresponds to the parameters $\ln r=20/3$ and $D=44/3$ in Eq. (2). The narrow high-temperature loop corresponds to the crossings of the spinodal lines in the horn region, while the wide low-temperature loop corresponds to crossings of the spinodals in the negative-$H$ region. In order to calculate these hysteresis loops, at each point along the hysteresis path we first located the local maximum in $F(m)$ that separates the two phases. Then, $\langle m\rangle$ and $\langle|m_{s}|\rangle$ were obtained by summing over $P(m,m_{s})$ as described in Sec. III.2. Although we do not show other examples of hysteresis loops here, we emphasize that our macroscopically constrained WL method enables the calculation of such loops for any value of $A$ and any choice of hysteresis path, solely based on the DOS data for the pure Ising antiferromagnet, without any further MC simulations. The hysteresis loop shown here is fully consistent with the one obtained by importance-sampling MC simulations for the same parameters in Ref. Rikvold et al. (2016) DD . The only significant differences are the slopes of the $\langle|m_{s}|\rangle$ curve where the path crosses the critical line, which in both cases are due to finite-size effects. On the other hand, finite-size effects in the positions of the spinodals are negligible, as discussed in Appendix A.

The phase diagrams for $A=6,7,8$ in Fig. 5 show several additional, noteworthy features. First, the phase diagrams shown are symmetric about the $T$ axis, with an exchange between FM$+$ and FM$-$. This is because the FM$+$ spinodal line is just touching the $T$ axis at $T=0$ for $A=4$ (Fig. 2) (c). Further increases of $A$ beyond $4$ will make a FM$-$ spinodal line show up in the positive $H$ region. Thus, the strong $AM^{2}/(2L^{2})$ causes a FM$-$ metastable region to appear in the positive $H$ field region. Figure 12(c)-(d) illustrate the case of a point lying on the coexistence line between the FM$+$ phase and AFM phases, inside the FM$-$ metastable region. The small drop in the free energy in (d) near $m=-1$ indicates the metastable FM$-$ phase. Figure 12(a)-(b) illustrate a point at $H=0$ and at a low $T$, where both FM phases are metastable.

Second, observe that the coexistence lines turn toward stronger $|H|$ when approaching the mean-field critical points (Fig. 5). This is because the disordered phase is more favorable than the FM phases at high $T$, so a stronger $|H|$ field is required to balance this effect.

Third, when $A$ increases, the mean-field critical temperature also increases, which makes the area of the horn region increase. This is because the ferromagnetic effects increase with $A$ according to the Hamiltonian (1), so a stronger disordering effect (higher temperature) is required to balance it.

Fourth, the coexistence line moves toward lower $|H|$ as $A$ increases, which makes the stable AFM region shrink and the stable FM regions grow. This is because strong $-AM^{2}/(2L^{2})$ stabilizes the ferromagnetic phases at lower $|H|$. The coexistence line for $A=8$ is at $H=0$ for low $T$ (Fig. 5(c)). In that case, the two AFM phases and the two FM phases are equally probable as shown in Fig. 13(a). When $T$ increases to a high enough value, disorder effects start to show up, making $|m_{s}|$ decrease from 1 (Fig. 13(c)). At low $|H|$ and high $T$, the disordered phase is preferred over the ferromagnetic phase. This effect starts to show up before reaching the critical temperature, making the coexistence line turn away from $H=0$ before it crosses the critical line, as shown in Fig. 14(a).

Fifth, the FM$-$ spinodal line continues moving toward higher $T$ when $A$ increases as the ferromagnetic phase is getting stronger. At $A=8$, the FM$-$ spinodal line has moved above the critical line. This produces a region (Fig. 14(a)) that is stable in the FM$+$ phase, and metastable in both the FM$-$ and disordered phases (region VIII), and another region that is stable in the disordered phase, and metastable in both FM$\pm$ phases (region VII).

Observe from Figs. 6 and 14(a) that the coexistence line makes a relatively large bend at the critical end-point. This is because after passing through this point, the AFM phase changes to the disordered phase, which is favored at high temperature, making the coexistence line have a smaller slope. Therefore, a relatively large bend in the coexistence line is found at the critical end-point. This agrees with the previously observed result that $d^{2}H/dT^{2}$ along the coexistence line reaches a maximum at the critical end-point Wilding (1997a, b); Tsai et al. (2007). Note that the location of the coexistence lines given by Rikvold et al.’s Ref. Rikvold et al. (2016) is different from the current result for $A=7$. The former result may be due to incomplete ergodic sampling by the mixed-start importance-sampling MC method used in that work to locate the coexistence lines. This might also affect experimental attempts to accurately detect phase coexistence. Further analysis of the discrepancy between the importance-sampling MC using the mixed-start method and the present method in locating coexistence lines is in progress Chan et al. (b).

IV.3 Transitional long-range interaction strength $A=8.1,9$

From the ground-state analysis in Ref. Rikvold et al. (2016), $A=8$ is the dividing line for the stable phase at $T=0$. For $A>8$, the stable phase at $T=0$, $H>0$ can only be the FM$+$ phase. Figures 13(a)-(d) show that increasing $A$ from $8$ to $8.1$ makes the FM phases overtake the AFM phases and become the stable phases below the critical line. The Bragg-Williams mean-field approximation Bragg and Williams (1934, 1935) also suggests that phase diagrams having $A>8$ belong to the same group (large long-range interaction group) and possess the same nature Rikvold et al. (2016). While Ref. Rikvold et al. (2016) has already pointed out that the Bragg-Williams mean-field approximation fails in predicting the existence of the horn regions (Figs. 6 and 14), here we find that the existence of the horn region induces a range of transitional long-range interaction strengths, between the medium long-range interaction and the strong long-range interaction. $A=8.1$ (Fig. 14(b)) and $A=9$ (Fig. 14(c)) belong to this range.

In this transitional range of $A$, we notice several things. First, the coexistence lines still exist, but the FM phases have pushed them to meet the $T$ axis at high temperatures, and this intercept temperature increases with $A$ (Fig. 14). Second, while for $A=8$ and when $T$ is low, the AFM phases and the FM$\pm$ phases are equally stable along the $T$ axis (Fig. 13(a)). Increasing $A$ makes the FM$\pm$ phases overtake the AFM phases along the $T$ axis. Figure 13 demonstrates this by comparing four points on the $T$ axis for $A=8$ and $8.1$. Third, the FM phases push the two spinodal lines originating from the mean-field critical point toward $H=0$. As a result, at around $A=9$ (Figs. 14(c) and 16(a)), the disordered spinodal lines nearly touch the $T$ axis before the two mean-field critical points from the $\pm H$ side of the phase diagram coalesce at even higher $A$.

While Fig. 13(h) shows a point close to the coexistence line for $A=8.1$, which has the disordered phase spread to the two AFM corners without connecting to the two FM$\pm$ peaks, Fig. 15 shows a point close to the coexistence line for $A=9$, which has the disordered phase connected to the two FM$\pm$ peaks. The connecting bridge should disappear and the three peaks should become sharper, as the system size is increased.

IV.4 Strong long-range interaction, $A=9.5,11$

When the long-range interaction is sufficiently strong, the two mean-field critical points in the $\pm H$ horn regions will coalesce into one critical point as shown in the phase diagrams for $A=9.5$ and $11$ in Fig. 16(b)-(c). Above this mean-field critical temperature, the system is in a disordered phase. If we increase $H$, the system undergoes a continuous crossover from the disordered phase to the FM$+$ phase, but there is no sharp transition point. The combined mean-field critical point is also the end point of the FM$\pm$ spinodal lines. When $H>0$ and $T$ is below the FM$-$ spinodal line (region IX in Fig. 16), the marginal probability density $P(m)$ has two peaks, and the system has a stable disordered/FM$+$ phase and a metastable disordered/FM$-$ phase (Fig. 17). As there is a continuous crossover between the disordered phase and the FM$+$ phase above $T_{c}$, it is natural that near the mean-field critical point, the marginal probability density has a large peak at a value of $|m|$ that is smaller than 0.5. Moreover, the metastable phase can show very strong disordered properties, so we consider the metastable phase below the FM$-$ spinodal line to be a disordered/FM$-$ phase. The topology of the phase diagrams for $A=9.5$ and $11$ is the same as found for $A=10$ in Ref. Rikvold et al. (2016).

Fig. 18(a)-(b) show probability densities at the coalesced mean-field critical point. It is found by extrapolation of the FM$-$ spinodal line and Eq. (19). Note that, as the critical point is in the mean-field universality class, at $T=6$, which is below the critical point for $A=11$ as shown in Fig. 18(c)-(d), we regard it as having stable FM$\pm$ phases, connected by fluctuations resembling the disordered phase. However, we do not regard the system as having a metastable disordered phase. The fluctuation connection has disappeared at around $T=5.5$. At $T=2.5$ as shown in Fig. 19(a)-(b), the system is close to the AFM and disordered spinodal line, the free energy in (b) shows a flat maximum around $m=0$. Figure 19(c)-(d) shows the case at $T=2.35$ for $A=11$, which is a point in region X in the phase diagram of Fig. 16(c). The free-energy contour and the free-energy drop near $m=0$ indicate the existence of the metastable disordered phase. Further reduction in $T$ below the critical line brings the system to the stable FM$+$ phase with metastable AFM phases, i.e. region XI in the phase diagram of Fig. 16(c), as shown in Fig. 19(e)-(f) for $T=2.2$. The free-energy contour, and the drop in free energy near $m=0$, indicate the existence of the metastable AFM phases. As $A$ increases (Fig. 16), the disordered/AFM spinodal line merges with the critical line. We expect region X, the disordered metastable phase region, to disappear when $A$ becomes very large.

V Summary and Conclusion

In this paper we have presented detailed phase diagrams, free-energy landscapes, and order-parameter distributions for a model SC material with antiferromagnetic-like nearest-neighbor and ferromagnetic-like long-range interactions Rikvold et al. (2016), covering a wide range of temperatures $T$, fields $H$, and long-range interaction strengths $A$. This was accomplished with a relatively modest computational effort by a recently developed, Macroscopically Constrained WL method for systems with multiple order parameters Chan et al. (2017). The method produces DOS for given values of the system energy $E$, magnetization $m$, and staggered magnetization $m_{s}$ for a square-lattice Ising antiferromagnet (i.e., $A=0$) in zero field. The DOS for arbitrary values of $H$ and $A$ are then found by a simple transformation of $E$ [Eq. (5)], without the need for additional simulations. From the transformed DOS, we obtain free-energy landscapes and $(H,T)$ phase diagrams, including metastable regions important to applications of SC materials Ohkoshi et al. (2011); Chakraborty et al. (2014); Kahn and Jay Martinez (1998); Linares et al. (2012). Topologically different phase diagrams are obtained, depending on the strength of $A$. For $A=0$, the numerically well-known phase diagram for the square-lattice antiferromagnet is recovered (Fig. 1).

For weak long-range interactions, $0<A\lesssim 4$, the high-temperature critical line terminates in a tricritical point at a nonzero temperature, from which sharp spinodal lines marking the extent of metastable phase regions extend to $T=0$ (Fig. 2). In this parameter range, the phase diagram is topologically identical to what is predicted by a simple Bragg-Williams mean-field approximation as discussed in Ref. Chan et al. (2017).

At a value of $A$ between 4 and 6 (which we have not attempted to determine accurately), the tricritical point decomposes into a critical end-point and a mean-field critical point at a higher temperature. The resulting horn structure of the phase diagram, which is not seen in simple Bragg-Williams mean-field calculations, is illustrated in Fig. 5 for the intermediate interaction strengths, $A=6$, 7, and 8. The phase diagram obtained for $A=7$ (Fig. 6) is in excellent agreement with that obtained by computationally intensive importance-sampling MC simulations in Ref. Rikvold et al. (2016). The only clear difference is the shape of the AFM/FM coexistence lines. A detailed investigation of this issue is in progress Chan et al. (b). (Very recently, horn regions and asymmetric, two-step hysteresis loops, analogous to those seen in the model studied here, have also been observed for a model with antiferromagnetic-like nearest-neighbor interactions and genuine elastic interactions NISH17 .) The horn structure gives rise to asymmetric, two-step hysteresis loops (see example in Fig. 11) that are similar to experimental observations in several different SC materials Köppen et al. (1982); Zelentsov et al. (1985); Petrouleas and Tuchagues (1987); Bousseksou et al. (1992); Jakobi et al. (1992); REAL92 ; Boinnard et al. (1994); Bolvin (1996); Chernyshov et al. (2004); Bonnet et al. (2008); Pillet et al. (2012); Buron-Le Cointe et al. (2010); Lin et al. (2012); Klein et al. (2014); Chernyshov et al. (2003); Huby et al. (2004); Shatruk et al. (2015); Brooker (2015).

For $A>8$, the AFM phase is no longer a possible ground state of the model. In the transitional region, $8<A\lesssim 9$, the horn region shrinks as shown in Fig. 14, until the two mean-field critical points coalesce into a single critical point at $H=0$ for a value of $A$ somewhere between 9 and 9.5. (This value we also have not attempted to determine accurately.) To our knowledge, this regime of transitional interaction strengths has not been investigated before. Phase diagrams for the strong-interaction case, represented by $A=9.5$ and 11, are shown in Fig. 16. These are topologically identical to the one shown for $A=10$ in Ref. Rikvold et al. (2016). We believe our results can contribute to the interpretation of the fascinating phase diagrams and hysteresis loops observed in many SC materials and other systems with competing short- and long-range interactions.

ACKNOWLEDGMENTS

The Ising-ASFL model was first proposed by Seiji Miyashita, and we thank him for useful discussions. The simulations were performed at the Florida State University High Performance Computing Center. This work was supported in part by U.S. National Science Foundation grant No. DMR-1104829.

Appendix A Finite-size effects and error estimates

The questions of finite system sizes and error estimates are intimately connected, and it is reasonable to ask whether the system size of $L=32$ that we use here is sufficient to ensure reliable results. The fourth-order Binder cumulant presumably leads to cancellation of leading corrections to scaling Binder (1981) and is a remarkably accurate method to locate critical points. The most general way to utilize the method is to look for the crossings between plots of cumulant vs temperature or field for different system sizes. However, the model studied here fulfills all the symmetry requirements to yield a fixed-point value of $0.6106924(16)$ Kamieniarz and Blöte (1993); Chen and Dohm (2004); Selke and Shchur (2005); Salas and Sokal (2000). As a consequence, it is possible to obtain good estimates of critical points as the phase points where the cumulant is near this value for a single system size, as we have done here. This is demonstrated in Fig. 20, where we compare critical lines obtained here using the macroscopically constrained WL method with $L=32$, with those obtained in Ref. Rikvold et al. (2016) by importance-sampling MC using the standard method of cumulant crossings for $L\leq 1024$. The differences are indeed very small, and although they are included as error bars in all the phase diagrams shown in this paper, they only exceed the symbol size in the lower right quadrant of the enlarged view of the horn region for $A=7$, shown in Fig. 6. The finite-size effects are even smaller for the spinodal lines (not shown here), and again the error bars obtained from the differences with the results of Ref. Rikvold et al. (2016) are only visible in Fig. 6. Statistical errors were reduced below the level of the finite-size effects by averaging the DOS over ten independent macroscopically constrained WL simulations as described in Appendix C of Ref. Chan et al. (2017).