Phase Diagram for a model of Spin-Crossover in Molecular Crystals

TL;DR

This paper introduces a simple theoretical model for spin-crossover in molecular crystals, analyzing its phase diagram through mean field and Monte Carlo methods to understand the transition behaviors and hysteresis phenomena.

Contribution

It provides a combined mean field and Monte Carlo analysis of a new model for spin-crossover, elucidating the phase diagram and transition properties.

Findings

Good agreement between analytical and simulation results in certain parameter regions

Mean field approach breaks down where correlations and cooperativity are significant

Model captures key features like abrupt, gradual, and incomplete transitions

Abstract

Spin-crossover has a wide range of applications from memory devices to sensors. This has to do mainly with the nature of the transition, which may be abrupt, gradual or incomplete and may also present hysteresis. This transition alters the properties of a given sample, such as magnetic moment, color and electric resistance to name some. Yet, a thorough understanding of the phenomenon is still lacking. In this work a simple model is provided to mimic some of the properties known to occur in spin-crossover. A detailed study of the model parameters is presented using a mean field approach and exhaustive Monte Carlo simulations. A good agreement is found between the analytical results and the simulations for certain regions in the parameter-space. This mean field approach breaks down in parameter regions where the correlations and cooperativity may no longer be averaged over.

Figures

Tables

Table 1. Table 1: Mean spin at low temperature.

$\alpha$	$\lambda$	$n$
$\sqrt{5/2}$	$0.25$	$2$
$\sqrt{5/2}+1/2$	$0.25$	$2$
$\sqrt{5/2}-1/2$	$0.25$	$0$ or $2$
$\sqrt{5/2}$	$0.5$	$2$
$\sqrt{5/2}+1/2$	$0.5$	$2$
$\sqrt{5/2}-1/2$	$0.5$	$2$

Equations

$$\hat{H}=\sum _{i=1}^N{\hat{H}}_i+\sum _{⟨i,j⟩}{\hat{V}}_{ij},$$

$${\hat{H}}_i=ϵn_i+\left({\hat{a}}_i^{†}{\hat{a}}_i+\frac{1}{2}\right)-\alpha n_i\left({\hat{a}}_i^{†}+{\hat{a}}_i\right),$$

$${\hat{V}}_{ij}=-\frac{\lambda }{4}\left({\hat{a}}_i^{†}+{\hat{a}}_i\right)\left({\hat{a}}_j^{†}+{\hat{a}}_j\right)$$

$$\hat{H}={\hat{H}}_e+{\hat{H}}_{ph},$$

$${\hat{H}}_e$$

$${\hat{H}}_{ph}$$

$$U_{ij}$$

$$s(\mathbf{q})$$

$$H_e^{mf}$$

$$Z_{mf}$$

$$\sum _{i,j}{U}_{ij}\approx -N\frac{{\alpha }^2\lambda z}{2}\frac{1}{1-\frac{2}{\lambda z}}.$$

$$Z_{mf}=\exp \left[-\beta \frac{n^2{\alpha }^{2}N}{1-\frac{\lambda z}{2}}\right]\times$$

$${\left(g_0+g_1{e}^{-\beta \left(ϵ-\frac{2n{\alpha }^2}{1-\frac{\lambda z}{2}}\right)}+g_2{e}^{-2\beta \left(ϵ-\frac{2n{\alpha }^2}{1-\frac{\lambda z}{2}}\right)}\right)}^N.$$

$$f_{mf}(n,T)$$

$$f_{mf}^{SR}(n,T)$$

$$f_{mf}\approx ϵ\eta -\frac{{\alpha }^2}{1-\frac{\lambda z}{2}}\left({\eta }^2-{\left(n-\eta \right)}^2\right)-T\ln \sum _{i=0}^2{g}_i,$$

$$f_{mf}^{SR}\approx \{\begin{array}{l}{n}^2{\alpha }^2\left(1+\frac{\lambda z}{2}\right),\text{for }n<ϵ/2{\alpha }^2\left(1+\lambda z/2\right)\\ (n-2{)}^2{\alpha }^2\left(1+\frac{\lambda z}{2}\right)+2ϵ-T\ln g_2-4{\alpha }^2\left(1+\frac{\lambda z}{2}\right),\\ \text{for }n>ϵ/2{\alpha }^2\left(1+\lambda z/2\right)\end{array}$$

$$n_{int}(T)=\frac{2ϵ-T\ln g_2}{4{\alpha }^2\left(1+\frac{\lambda z}{2}\right)}.$$

$$T=f_{mf}^{SR}(n_{int}(T),T),$$

$$s(q\to 0)=\frac{z}{2}.$$

$$U_{ij}=-\frac{{\alpha }^2\lambda z}{2N}\frac{1}{1-\frac{2}{\lambda z}}.$$

$$s(\mathbf{q})=\sum _{k=1}^d\cos (q_k)\approx \frac{z-a^2{q}^2}{2}.$$

$$U_{ij}\approx \frac{{\alpha }^2\lambda z^2}{8\pi |{\mathbf{r}}_i-{\mathbf{r}}_j|\rho a^2}f\left(|{\mathbf{r}}_i-{\mathbf{r}}_j|,\kappa \right),$$

$$f(r,\kappa )=\{\begin{array}{l}{e}^{-\kappa |{\mathbf{r}}_i-{\mathbf{r}}_j|/a}{\kappa }^2>0\\ \cos \left(\kappa |{\mathbf{r}}_i-{\mathbf{r}}_j|/a\right){\kappa }^2<0\end{array}.$$

$${\kappa }^2=\frac{2}{\lambda }-z$$

$$V(r_{ij})=\{\begin{array}{l}{U}_{ij},r_{ij}\ne 0;\\ 0,r_{ij}=0\end{array}.$$

$$\sum _{i,j}V(r_{ij})=\iint d^3\mathbf{r}{d}^3{\mathbf{r}}^{\prime }\rho (\mathbf{r})\rho ({\mathbf{r}}^{\prime })V(|\mathbf{r}-{\mathbf{r}}^{\prime }|),$$

$$\rho (\mathbf{r})=\frac{1}{V}\sum _{i=1}^N\delta \left(\mathbf{r}-{\mathbf{r}}_i\right).$$

$$\{\begin{array}{l}\mathbf{r}=\frac{{\mathbf{r}}_0}{2}+\mathbf{R}\\ {\mathbf{r}}^{\prime }=-\frac{{\mathbf{r}}_0}{2}+\mathbf{R}\end{array}.$$

$$\sum _{i,j}V(r_{ij})\approx \rho V4\pi {\int }_0^{\infty }d{r}_0{r}_0^{2}V(r_0).$$

$$\sum _{i,j}{U}_{ij}\approx N\frac{{\alpha }^2\lambda z^2}{2a^2}{\int }_0^{\infty }d{r}_0{r}_{0}f(r_0,\kappa ).$$

$${\alpha }^2\left(1+\frac{\lambda z}{2}\right)N+\sum _{i,j}{U}_{ij}$$

Full text

I Introduction

Research around several phenomena in the overlap between solid state and condensed matter have the peculiar behavior of having periods of time where it is quiet and some other periods of time where a lot of research is being done. Such is the case of spin-crossover (SCO) phenomena, which spans broadly nine decades (see Ref. Halcrow (2013) for a compilation of research over the years together with Refs. Gütlich and Goodwin (2004); Halcrow (2011)). The SCO phenomenon is the the transition between a low spin (LS) and a high spin (HS) state on a metal ion with $d^{4}-d^{7}$ electronic configuration. Experimental evidence suggest that the cause of this transition has to do with the competition between the strength of the field of ligands and the spin-coupling energy between electrons Shriver, Atkins, and Langford (1994). For instance, octahedral compounds Transition Metals Series $3d^{4-7}$ present this type of transition, in which they can be in a high or low spin state, depending on whether the ligand field is stronger or weaker than the matching energy.

Furthermore, in the case of thermally induced SCO transition, the free energy difference between both states should be of the order of thermal energy, i.e., $k_{B}T$ (henceforth we consider $k_{B}=1$) Gütlich, Hauser, and Spiering (1994). In this sense, high temperature favors HS whereas low temperature favors LS. Interestingly, the phenomenon is rather ubiquitous in nature. For instance, in the case of the transition in $\text{Fe}^{II}$ complexes between $t_{2g}^{6}$ ($S=0$) and $t_{2g}^{6}$ ($S=2$) configurations, SCO is responsible for oxygen transport in hemoglobin and probably for the change under pressure of ferropericlase in the Earth’s mantle Yang et al. (2015). In solids, SCO can be found in many transition metal oxides, organometallic complexes, inorganic salts or organic radicals and has a cooperative nature, frequently leading to abrupt changes of macroscopic physical properties and hysteresis. Perhaps even more appealing are the applications of SCO which range from display and memory devices and electroluminescent devices to MRI contrast agent Létard, Guionneau, and Goux-Capes (2004); Halcrow (2013). Additionally, the use of SCO combined with the properties of nanoporous metal-organic frameworks may also be used in molecular sensing Halder et al. (2002).

The HS and LS phases have different properties that depend on the electronic distribution in 3d orbitals. In this sense, the HS to LS transition has a large impact on the physical properties of a material, such as the magnetic moment, the color, the dielectric constant as well as the electrical resistance, among other properties. In other words, features such as optical, vibrational, magnetic and structural differ between one phase and the other. Hence, measuring these properties serve as a proxy to measure and monitor the SCO induced by external perturbation, such as light, pressure or temperature, for instance Gütlich and Hauser (1990); Chernyshov et al. (2007); Konishi et al. (2008); Ohkoshi et al. (2011); Pinkowicz et al. (2015). Several techniques such as magnetic susceptibility measurements, as well as optical and vibrational spectroscopy of the kind of UV, IR, Raman and Mossbauer spectroscopy are used for this end. However, the holy grail is predicting the spin curve for a given material under cooling and heating together with the critical temperature and the hysteresis loop. There has been various efforts in this directions Chernyshov et al. (2004, 2007), yet there is still a lack of a theory, which should not come as a surprise given the great variety of materials having SCO. A complete microscopic description requires the consideration of three basic ingredients: i) the spin and vibrational states of individual octahedral complexes, ii) the interaction between them leading to cooperative effects and iii) the coupling with external factors such as temperature, pressure as well as external electric or magnetic field . In this regard, for instance, in oxides cooperativity is attributed to electronic exchange while in molecular crystals to electron-phonon coupling Nesterov et al. (2017).

In a previous paper by two of the authors Rodríguez-Castellanos, Plasencia-Montesinos, and Reguera-Ruiz (2018), a theoretical approach to SCO in mononuclear molecular crystals containing $Fe^{II}$ ions was presented, where a simple effective interaction between neighboring local breathing modes was considered. Only electrons in eg states are linearly coupled to vibrations, being that the case there is no need for two breathing modes with different vibrational frequencies and coupling parameters. Furthermore, decoupling breathing modes with successive canonical transformations leads to a lattice model where short range and long range ferromagnetic and antiferromagnetic interactions arise in a strightforward manner. In the present paper we further study the phase diagram by means of Monte Carlo simulations and analytical derivations. The structure of the paper is as follows: In the next section we present the model and its features, in section III we solve the model in the thermodynamical limit, in section IV we compare the analytics with the simulations and discuss the results, finally section V is for conclusions.

II Model

Consider a $d$- dimensional periodic array of $i\;(1,...,N)$ mono-nuclear metal complexes, each containing an ion $\text{Fe}^{\text{II}}$ in an octahedral site surrounded by non-magnetic ligands. The number of electrons occupying $e_{g}$ states at the $i$th ion will be denoted with $n_{i}$ which represent the spin on the lattice site $i$. We further consider the degeneracy for each state, $g_{n}$, such that $g_{0}=1$, $g_{1}=9$ and $g_{2}=15$. Since the charge density distribution of the antibonding $e_{g}$ states is more localized near their octahedral neighbors, these electrons are strongly coupled with a local breathing vibration mode described by the operators $\hat{a}_{i}^{\dagger}$ and $\hat{a}_{i}$. Breathing modes at neighboring sites interact via acoustic phonons. Thus, the effective Hamiltonian for this electron-local vibrations system is:

[TABLE]

where

[TABLE]

is the Hamiltonian for the $i$th ion, a harmonic oscillator and a term coupling both. Whereas,

[TABLE]

takes into account the coupling between breathing modes localized in neighbor sites. The parameter $\epsilon>0$ is the excitation energy per $e_{g}$ electron, which is obtained by subtracting the splitting energy between $t_{2g}$ and $e_{g}$ states and the pairing energy $P$ in $t_{2g}$ states, while $\alpha$ and $\gamma$ are coupling parameters. We will assume $\epsilon=10$ unless stated otherwise explicitly. Moreover, we consider the phonon energy of the local breathing mode when the ion is in HS-state $\hbar\omega=1$ as well as the Boltzmann constant $k_{B}=1$.

Coupling to local modes induces virtual transitions that renormalize electron energies, give rise to an effective electron-electron interaction and shift atomic positions (see below). The interaction Hamiltonian $\hat{V}_{ij}$ is the simplest approximation to an effective inter-site coupling which could result from averaging over degrees of freedom (acoustical phonons) connecting local breathing modes at neighboring sites Palii et al. (2015).

Breathing modes can be decoupled by canonical transformations of their creation and annihilation operators, leading to a system with an effective electron-electron interaction and independent dispersive phonons:

[TABLE]

where

[TABLE]

Here, the notation $j(i)$ and $\langle i,j\rangle$ denote the first neighbor $j$ for a given $i$ and pairwise neighbors, respectively and $z$ is the coordination number, which we assume equal to $6$ whenever is not explicitly specified. The actual derivation of Eqs. (5), (6), (7) and (8) is documented in SN 3.

III Mean field

In this section we present a mean field approximation of the electron-electron Hamiltonian produced in Eq. (5). The idea is to consider first order fluctuations in spins, denoted as $\delta n_{i}$, while neglecting second order fluctuations as in Ref. Cardy (1996). To this end, we write $n_{i}=n-\delta n_{i}$ where $n$ is the mean value of the spins. It is easy to show that the electron-electron Hamiltonian in Eq. (5) becomes

[TABLE]

Notice that under this approach, the spin-spin product no longer appears and, instead, spins interact with a mean field proportional to $\sim n$. Then, computing the partition function using the mean field Hamiltonian, using Eq. (9), is direct, and yields:

[TABLE]

Notice in the partition function (Eq. (10)) how the dependence on the index $i$ from the product falls with the long range interaction, i.e., $U_{ij}$. In the case for $N\gg 1$ we may approximate the sum of long range interactions as

[TABLE]

The derivation of Eq. (11) is shown in the Supplementary Note (SN) 1.

Hence, using the approximation presented in Eq. (11) we further simplify the partition function to,

[TABLE]

From Eq. (12) we may compute any thermodynamic quantity straightforward. Although we are using a mean field approach, the mean field free energy will depend upon the mean spin, which is not a proper free energy. Nevertheless, it may be proved that in the thermodynamical limit ($N\rightarrow\infty$) the minimum mean field free energy value as a function of $n$ equals the actual free energy Kardar (2007).

The mean field free energy per molecule, which we denote as $f_{mf}(n,T)$, is

[TABLE]

Similarly, in the case where we neglect long-range interaction, it is not difficult to realize that the mean field free energy, $f_{mf}^{SR}(n,T)$, yields:

[TABLE]

Notice from Eqs. (13) and (14) that for $\lambda\ll 2/z$, it follows $f_{mf}\approx f_{mf}^{SR}$. Conversely, under this approximation when $\lambda>2/z$ the long range interaction does not decay with distance. Hence, this imposes a constraint over the system which makes it impossible to draw conclusions from a mean field approach and, in fact, we shall see that the system is highly size-dependent (this is shown and discussed in SN 1).

Now, in the case of high temperatures, i.e., $T\gg\epsilon-2n\alpha^{2}/\left(1-\lambda z/2\right)$, the mean field free energy becomes

[TABLE]

with $\eta=\sum_{i=0}^{2}g_{i}n_{i}/\sum_{i=0}^{2}g_{i}$ being the high temperature mean spin. Notice that the previous Eq. (15) has a minimum at $n=\eta$ provided $\lambda<2/z$. This is shown in Fig. 1 where we present a contour plot of the mean field free Energy per particle (Eq. (13)) in a $T-n$ diagram. Moreover, this is in agreement with the fact that for small values of $\lambda$, $f_{mf}\approx f_{mf}^{SR}$ and in the case of short range interactions, at high temperatures the system must behave as a set of uncoupled spins. This implies that the mean spin is given by $\sum_{i=0}^{2}p_{i}n_{i}$, where $p_{i}=g_{i}/\sum_{m=0}^{2}g_{m}$ is the probability obtained simply from the degeneracy of each configuration neglecting any coupling in the model.

An interesting feature of the model is the fact that in the case where $\lambda>2/z$, the point $n=\eta$ becomes a maximum and the global minimum is located at $n=0$, i.e., at high temperatures the mean spin becomes zero. This is shown in Fig. 1 where we present a contour plot of the mean field free Energy per particle (Eq. (13)) in a $T-n$ diagram. Thus, the condition $\lambda>2/z$ preempts the maximization of entropy at high temperatures. In the SN 1 we show that, in fact, for $\lambda>2/z$ the long range interactions go as $N^{1/4}\exp(N)$, hence any mean field approach fails in describing the model in this parameter region. Hence, we should constrict the model to $\lambda<2/z$.

Critical Temperature

In this section we compute the critical temperature. In particular, we show that there are some parameter values in which at low temperatures the system is equally likely to have a transition where all spins are $2$ or [math] mimicking a two-level system. To understand this, let us consider the mean field free energy in the case of short range interactions (Eq. (14)), however, the procedure is the same in the case of long range interactions mean field free energy. Notice that the exponentials inside the logarithm may go to zero or infinity for low temperatures depending on the values of the parameters as well as the mean spin. Let us consider low temperatures and that $n<\epsilon/2\alpha^{2}\left(1+\lambda z/2\right)$. Then, for low temperatures, we may approximate Eq. (14) with solely the first term. Similarly, in the case where $n>\epsilon/2\alpha^{2}\left(1+\lambda z/2\right)$ we may approximate Eq. (14) thus to write:

[TABLE]

In Fig. 2 we show the comparison between Eq. (14) and the approximation (16) for some fixed parameter values. Now, given we are using mean field, the minima of the mean field free energy should correspond with the actual free energy. Thus, we should expect that the system is likely to choose low or high spin for the same given parameter-set when two things are fulfilled, namely,

• •

The mean field free energy at low spin is equal to the mean field free energy at high spin.

• •

The temperature is of the order of the maximum of the free energy.

The first condition translates into equating Eq. (16) for $n=0$ with itself for $n=2$. From this we obtain the temperature $T_{f}=2\left(\epsilon-2\alpha^{2}\left(1+\lambda z/2\right)\right)/\ln g_{2}$. The second condition implies that $T_{f}$ should be of the order of the maximum value in the free energy. In Fig. 2 we have plotted the mean field free energy (Eq. (16)) at temperature $T=T_{f}$ as well as the temperature $T_{f}$ for a given set of parameters (see caption), in particular, we fixed $\alpha\approx 1.5$. Notice that the maximum of the free energy is of the order of the temperature which is $\approx 1.5$. As we will see later, the numerical simulations yield a critical temperature of $\approx 1$ for the same parameter values.

We may further compute the critical temperature under this criteria. We do this as follows: first, we locate the intersection between Eq. (16) when $n<\epsilon/2\alpha^{2}\left(1+\lambda z/2\right)$ and itself when $n>\epsilon/2\alpha^{2}\left(1+\lambda z/2\right)$. Let us denote this intersection as $n_{int}(T)$. It is a feasible task to obtain

[TABLE]

Then, solving for $T$ the Eq.

[TABLE]

yields the critical temperature $T_{c}$. In Fig. 3 we have plotted the critical temperature vs $\alpha$ from which one may appreciate that $T_{c}$ decreases as $\alpha$ increases, which is qualitatively consistent with the numerical simulations to be discussed in the next section.

Additionally, we may predict if the transition is to LS or HS. This is done by noticing that if the intersection between Eq. (16) when $n<\epsilon/2\alpha^{2}\left(1+\lambda z/2\right)$ and itself when $n>\epsilon/2\alpha^{2}\left(1+\lambda z/2\right)$ happens at $n>1$ then the mean field free energy has a minimum at $n=0$. Conversely, if it happens at $n<1$ then the minimum occurs at $n=2$. In Fig. 3 we have plotted $n_{int}(T_{c})$ vs $\alpha$. Notice that for $\alpha<1.5$ the system goes to LS while for $\alpha>1.5$ the system goes to HS. However, for $\alpha\approx 1.5$ the system is likely to go HS or LS.

A similar analysis may be done in case of $\lambda<2/z$ with long range interactions and the outcome is essentially the same. However, the critical temperature is somewhat lower than when neglecting long range interactions.

IV Discussion

In this section we compare the mean field results with the numerical simulations. This were done in C++ and Julia using a Monte Carlo algorithm with a Metropolis test Newman and Barkema (1999) for system sizes $10^{3},\;15^{3},\;20^{3}$ and $100^{3}$ spins with periodic boundary conditions. We fixed the parameters $\epsilon=10$ and $z=6$ while considering, both, $\lambda=0.25(<2/z)$ and $\lambda=0.5(>z/6)$ and then for a given fixed set of the previous values we fixed $\alpha$ in the range from $0.1$ to $3.0$ in the simulations performed in Julia, while in the simulations performed in C++ we increased that range to $4.0$. We initiate the simulations at $T_{in}=15$ with all spins having a value equal to $2$. Then, we start lowering the temperature in steps of $\Delta T=0.25$ and $\Delta T=1.0$ in the numerical simulation done in C++ and Julia, respectively, until reaching $T\approx 0$.

Now, in Fig. 4 we show contour plots of the mean field Eqs. (13) and (14) for different values of $\alpha$ and $\lambda$ in a $T$-$n$ diagram. As was discussed in the previous section, mean field approach predicts that, in the thermodynamic limit, when long range interactions are neglected for $\alpha\approx\sqrt{5/2}$ the system is likely to be in a LS or HS state as shown in Fig. 4. In Fig. 5 we have plotted the mean spin obtained from the simulations, which we denote as, $\langle n\rangle$ vs temperature for the aforementioned case obtained from the simulations. Notice that the black data points in Fig. 5, corresponding to $\alpha=1.5$, show mean spin equal to [math] at low temperatures. Conversely, the black data points in Fig. 5 show mean spin equal to $2$ at low temperatures and also correspond to $\alpha-1.5$, i.e., for the same parameter values the system is equally likely in having mean spin equal to $2$ as well as equal to [math]. This is in agreement with the mean field Eq. (13) which has been plotted in Fig. 4. We have summarized all this in Table 1.

In Fig. 6 we show the mean spin $\langle n\rangle$ vs temperature obtained from our simulations for different system sizes and parameter values which are specified in the legend of the figure. Fig. 6 (a) to (c) correspond to $\epsilon=10,\;z=6$ and $\lambda=0.25$ neglecting long range interactions while considering different system sizes specified in the legends of each plot. Each curve corresponds to a fixed value of alpha that spans from $0.1$ to $4.0$ (Fig. 6 shows the color code we are using), as was described in detail at the beginning of this section. The first thing to notice is that for the system size we are considering, the results are robust. Secondly, notice that although initially all spins have values equal to $2$, after the system is equilibrated the mean spin value ranges from $\sim 1.25$ for low $\alpha$-value, to $2.0$ for high $\alpha$-value. Finally, the black data points correspond to $\alpha\approx 1.5$, that is the value predicted by mean field theory where the system can go to spin [math] or spin $1$. Fig. 6 (d) to (f) correspond to $\epsilon=10,\;z=6$ and $\lambda=0.5$ neglecting long range interactions while considering different system sizes specified in the legends of each plot. The black data points correspond to $\alpha=1.2$. The mean field approach predicts that for $1.2<\alpha<1.3$ the system may go to spin $2$ or [math] at low temperatures, which is consistent with these simulation results. Altogether, for the parameter values considered, the comparison between the mean field approach and the simulation results work seemingly well. Besides, the critical temperature predicted by the mean field approach is close to the values obtained in the Monte Carlo simulation, though it is advisable to not relay in mean field approaches to compute such quantities, in general.

Now, when we consider long range interactions in the Monte Carlo Simulations, the results are not as obvious as the predicted from mean field approach. Fig. 6 (g) to (i) correspond to $\epsilon=10,\;z=6$ and $\lambda=0.25$ considering long range interactions as well as different system sizes specified in the legends of each plot. The first thing to notice is that the results obtained in the case where $N=10^{3}$ does not match those obtained in the case where $N=20^{3}$, thus a bigger system size is required for results to converge to the thermodynamical limit. Furthermore, the black data points correspond to $\alpha\approx 1.5$, i.e., around this value the system is likely to go to a spin value equal to $2$ or spin [math]. However, the mean field approach predicts a lower $\alpha$-value and this has to do with the fact that we are relying on the critical temperature obtained by the mean field approach to obtain the $\alpha$-value where the system can go to spin $2$ or spin $1$.

Finally, Fig. 6 (j) to (i) correspond to $\epsilon=10,\;z=6$ and $\lambda=0.5$ considering long range interactions as well as different system sizes specified in the legends of each plot. Notice that the results are profoundly system size dependent, which is why mean field approaches are fruitless. Nonetheless, this is in agreement, in principle, with the theoretical prediction of the model for $\lambda<2/z$ where it is obtained long range oscillating interactions (see SN 1). Moreover, system size dependence may be related to cooperativity properties, a feature characteristic of SCO phenomena Murray and Kepert (2004).

V Conclusions

In this work we have studied the spin-crossover transition using a simple model derived from quite general assumptions. We have presented a thorough discussion on the parameter space using a mean field approach and Monte Carlo simulation. Using the mean field we were able to obtain the critical temperature. When long range interactions were neglected, the mean field critical temperature compared well with the Monte Carlo simulation results. We further showed, both, through the mean field approach and the Monte Carlo simulations that for a small region in the parameter space, at low temperature the model may be in LS or HS state. When long range interactions are considered, the mean field approach loses accuracy as expected since long range interactions facilitates cooperativity. We further showed that when this happens, the spin vs $T$ curve is highly system size dependent. Hysteric effects have not been addressed here, but we may speculate that it might be possible in the parameter region where the system is likely to be in a LS or HS state. However, we leave that for future work.

Acknowledgments

We gratefully thank DGAPA-PAPIIT project IN102717. J.Q.T.M. acknowledges a doctoral fellowship from CONACyT.

Supplementary Note 1

Here we derive a mean field approximation for the long range interaction $U_{ij}$ First, notice that for $\bm{q}\rightarrow 0$,

[TABLE]

Therefore, when $N\gg 1$, then the long range interactions become

[TABLE]

The way Eq. (11) is obtained seems rather naive. Thus, let us provide a different approach. Notice that in the case of a cubic lattice with lattice parameter $a$, in the thermodynamical limit we have

[TABLE]

Thus, the long range interaction may be expressed as,

[TABLE]

where $f(r,\kappa)$ is defined as

[TABLE]

and

[TABLE]

We may further approximate this. First, for consistency, let us define $V(r_{ij})$ as

[TABLE]

Now, notice that

[TABLE]

where $\rho(\bm{r})$ is the density, i.e.,

[TABLE]

Now, let us replace $\rho(\bm{r})\rightarrow\rho\equiv N/V$ and make a change in variables to the center of mass and the relative distance taken pairwise, namely,

[TABLE]

which has Jacobian equal to unity. Henceforth, it is straightforward to show that

[TABLE]

Applying these ideas to Eq. (S22) yields

[TABLE]

In the case $\kappa^{2}>0$ the result is straightforward and yields Eq. (11). However, for $\kappa^{2}<0$ one should consider a maximum length. Naturally, this maximum length should be of the order of $\sim N^{1/3}$, yet, this implies that $\sum U_{ij}\sim N^{4/3}e^{N^{1/3}}$ which is fine for small system size. But the whole approximation is based on the assumption of the thermodynamical limit ($N\gg 1$). Thus, mean field does not apply in this regime.

Now, by taking into account the long range approximation, notice that

[TABLE]

Supplementary Note 2

In this section we provide a decription on how to calculate the long range interacion $U_{ij}$ in square lattices following a weighted sum is reciprocal space, using the method proposed in Ref. Chadi and Cohen (1973). We need first to calculate,

[TABLE]

it follows,

[TABLE]

Observe that in $s(q)$, we need only to consider the lower triangle of the first quadrant where $\bm{q}$ has components $q_{x}>0$ and $q_{y}\leq q_{x}$, since $\cos(x)$ is an even function.

Now we need the sum,

[TABLE]

Again, due to symmetry reasons, we concentrate our attention into a lower triangle of the first quadrant where $\bm{q}$ has components $q_{x}>0$ and $q_{y}\leq q_{x}$. Suppose for example that our mesh has equally spaced 16 q points. Only three of them are nonequivalent, $\bm{q}_{1}=(1/8,1/8)$ appears 4 times, $\bm{q}_{2}=(3/8,3/8)$ appears 4 times and $\bm{q}_{3}=(3/8,1/8)$ 8 times. Thus each point in a sum has a weight,

[TABLE]

We now use our mesh in $\bm{q}$-space,

[TABLE]

where,

[TABLE]

or,

[TABLE]

where $\theta_{l}=\tan^{-1}(q_{l}^{x}/q_{l}^{y})$ and $q_{l}=\sqrt{(q_{l}^{x})^{2}+(q_{l}^{y})^{2}}$

The mesh can be further refined using this technique.

Let us calculate $U_{ij}$ in a cubic lattice. First,

[TABLE]

In a 3D cubic lattice, we use a mesh grid of 64 $q$ points. This can be reduced to compute only these $q$-points Chadi and Cohen (1973),

$\bm{q}_{1}=(1/8,1/8,1/8)$ with $W_{1}=1/8$

$\bm{q}_{2}=(3/8,3/8,1/8)$ and $W_{2}=3/8$

$\bm{q}_{3}=(3/8,1/8,1/8)$ and $W_{3}=3/8$

$\bm{q}_{4}=(3/8,3/8,3/8)$ , and $W_{4}=1/8$

[TABLE]

If required, consider points

[TABLE]

with $n$ an integer and then iterate. Then,

[TABLE]

where $n_{l}$ is the number of points obtained from $\bm{q}_{l}$ under all symmetry operations of the lattice point group.

Supplementary Note 3

In this section we show how to obtain the effective Hamiltonian (Eq. (6)). To this end, we consider a $d$-dimensional periodic array of $N$ mono-nuclear metal complexes which we tag with latin letters. Each site $i$ contains an ion $\text{Fe}^{+2}$ in an octahedral site surrounded by non-magnetic ligands. We denote the number of electrons occupying states eg at ion $i$th with $n_{i}$. These electrons are coupled to a local “breathing” vibration mode described by the creation and annihilation operators over site $i$ denoted as $\hat{a}_{i}^{{\dagger}}$ and $\hat{a}_{i}$, respectively. Breathing modes at neighboring sites are also coupled through their interaction with acoustic phonons. Thus, the effective Hamiltonian for this electron-local vibrations system may be expressed as:

[TABLE]

with

[TABLE]

In these expressions $\epsilon$ is the excitation energy per eg electron, while $\alpha$ and $\lambda$ are coupling parameters. Now, let us consider a transformation on the local phonon operators, such that,

[TABLE]

Therefore, Eq. (S44) becomes,

[TABLE]

Now, let us express the operators $\hat{b}_{i}^{{\dagger}}$ and $\hat{b}_{i}$ in Fourier space, i.e.,

[TABLE]

Notice that

[TABLE]

where we have used the identity,

[TABLE]

Here $\delta_{\bm{q}\bm{q}^{\prime}}$ is the Kronecker delta defined as

[TABLE]

Now, let us define $s(\bm{q})$ as

[TABLE]

where we denote $j(i)$ as the index $j$ running over the first-neighbor lattice sites of site $i$. Now, for a lattice with inversion symmetry, i.e., $s(\bm{q})=s(-\bm{q})$, the phonon part of the Hamiltonian, $\hat{H}_{ph}$, becomes

[TABLE]

In order to express Eq. (S53) in a canonical way, we further introduce the new operators $\hat{f}_{\bm{q}}^{{\dagger}}$ and $\hat{f}_{\bm{q}}$, such that,

[TABLE]

Substituting Eq. (S54) in Eq. (S53) and after some algebra, one obtains

[TABLE]

with

[TABLE]

Similarly, the electron-phonon coupled term may be decoupled by introducing the same transformations. It is not difficult to show that,

[TABLE]

Finally, we define $\hat{c}_{\bm{q}}^{{\dagger}}$ and $\hat{c}_{\bm{q}}$ as,

[TABLE]

and arrange terms to finally obtain,

[TABLE]