# 1D Spin-Crossover Molecular Chain with Degenerate States

Andrii Gudyma  Max Planck Institute of Microstructure Physics, Weinberg 2,
06120 Halle/Saale, Germany  Iurii Gudyma<EMAIL_ADDRESS>Physical,
Technical and Computer Sciences Institute of Yuriy Fedkovych Chernivtsi
National University, 58012 Chernivtsi, Ukraine

###### Abstract

A study of the one-dimensional molecular chain (MC) with two single-particle
degenerate states is presented. We establish connection of the MC with the
Ising model with phononic interactions and investigate properties of the model
using a transfer matrix method. The transfer matrix method offers a promising
pathway for simulating such materials properties. The role of degeneracy of
states and phononic interaction being made explicit. We analyze regimes of the
system and parameters of the occurring crossover. Here, we present exact
results for the magnetization per spin, the correlation function and the
effective volume of the system. We demonstrate possibility of existence of two
peaks in the specific heat capacity thermal behavior.

††preprint: AIP/123-QED

##  I Introduction

Quasi-one-dimensional systems play an essential role in nanotechnology.
Usually phase transitions do not take place in these systems. It is well known
that, in general, one dimensional systems with short range interaction do not
undergo a phase transition, with possible exception when the studied model is
considered at zero or infinite temperatures. The absence of phase transition
in such systems has been determined by wide variety of conditions. It means
that the manifestation of each case of phase transition in one-dimensional
(1D) systems demands a detailed study. For example, presence of hysteresis and
many-step crossover Spin-Crossover1d-1  were found in SCO materials. These
phenomena appear as a result of the spin-elastic frustration.  Depending on
the particular configuration of the magnetic molecules, a variety of
magnetization behaviors can be described even with a relatively simple model.
The importance of such studies is related to the occurrence of first-order
phase transition in molecular chains, that could find a wide variety of
practical applications, including further miniaturization of nanodevices and
nanomaterials Nanotechnology . Molecular building blocks are among the most
promising candidates for the future electronics SENTHILKUMAR2017176 .

When the energy gap between the ground state and the first excited state
becomes sufficiently perceptible and distances to the rest of the excited
states are much greater, then the two-state model may be used as the
theoretical framework. The two-state interaction model is formulated in terms
of Ising spin variables, $s_{n}=\pm 1$. Mathematically, this corresponds to
the application of the quasi-spin-$\frac{1}{2}$ Ising model for molecular
chains. If the lowest states of the system are degenerate, then we have a
degenerate Ising model. Note that, as a rule, the degeneracy of the high-
energy spin state is more significant. Ising model with the degenerated states
was introduced by Huang Huang1975 , and developed and studied numerically
HARRIS1985299 ; bousseksou1992ising ; boukheddaden2000one ; Schebarchov2014 ;
sienkiewicz2014finite ; cajahuaringa2019nonequilibrium . The equivalence of
the states’ degeneracy and temperature dependent effective magnetic field
acting on spin-crossover molecules was established bousseksou1992ising ;
boukheddaden2000one . Various mechano-elastic nasser2005diluted ;
konishi2008monte ; Lecomte2008Two-variable ; enachescu2009model ; ye2015monte
and Ising-like models linares1999analytical ; Boukheddaden_2007 ;
boukheddaden2008molecular  were developed to describe effects which appear in
spin crossover materials. Model which deserve special focus of interest is the
compressible Ising model zagrebnov1972spin ; Salinas_1973 ; Henriques_1987 ;
lehmann2012anomalous ; Spin-Crossover1d-1 . For this model in 1D the elastic
interactions were taken into account using transfer matrix technique
Linares2004 ; Boukheddaden2007  and this model was applied to spin crossover
materials Spin-Crossover1d-1 ; Spin-Crossover1d-2 ; Spin-Crossover1d-3 .

For example, in coordination Fe (II) linear chain compounds exhibiting thermal
spin crossover (SCO) transitions are associated with diamagnetic-to-
paramagnetic switching between low-spin (LS) and high-spin (HS) state
nebbali2018one ; lan2019thermal ; vinogradova2020cooperative . In other words,
SCO chains undergo drastic changes in the spin configuration, leading them to
switch in a reversible way from a low-spin state (LS) to a high-spin state
(HS). The SCO phenomenon represents the paradigm of bistability at molecular
level which emerges at the macroscopic scale, offering potential applications
in the development of new generations of electronic devices such as molecular
spintronics devices, nonvolatile memories, molecular sensors, displays, and
reversible switches raman2013interface ; matsumoto2014programmable ;
bairagi2016molecular ; mullaney2017spin ; kipgen2018evolution ;
ridier2020unprecedented ; kobke2020reversible . The energy difference between
HS and LS states is due to the competition between the crystal field
splitting, which prefers doubly occupied $d$ orbitals and, hence, LS, and
Hund’s first rule, which favors the HS state. At lower temperature the LS
state is dominant, while the HS state is preferred at higher temperature
halcrow2013spin ; gudyma2015kinetics . From this point of view, the study of
chain Fe (II) complex by x-ray-absorption fine-structure spectroscopy (XAFS)
afford opportunity to obtain detailed information on the phase transition from
the local-structure aspect yokoyama1998spin . The extended XAFS analysis
showed the abrupt change of the interatomic distance at the transition
temperature. Large distance differences between the LS and HS states are
ascribed to the electron occupation of the $3d$ levels in the HS state. The LS
state shows the ${}^{1}A_{1g}$ electronic configuration of
$(3dt_{2g})^{6}(3de_{g})^{0}$ where no electrons occupy the
$\sigma$-antibonding $3de_{g}$ orbitals, while the HS state exhibits the
${}^{5}T_{2g}$ configuration of $(3dt_{2g})^{6}(3de_{g})^{2}$, this leading to
a significant weakening of the metal-ligand bond. Videlicet, the coupling of
the molecular electronic state with lattice transformation together with
intermolecular interactions result in very cooperative transitions.

In this work we consider molecular chain of particles with two inner states.
In our model these two inner states may be degenerated. The particles are a
subject to a pair potential of general type which differentiates the inner
states of particles. We establish the connection between molecular chain and
the Ising-like model with phonon interactions. We treat the model analytically
using transfer matrix formalism. Here our goal is first to show the role of
degeneracy of states for scenario of spin transition and analyze parameters of
the spin crossover. We focus our attention on the thermal behavior of the
magnetization and specific heat capacity.

The outline of this work is as follows. Sec. II defines the model’s formalism.
In Sec.III we introduce the transfer matrix formalism and make finite $N$
analytic calculations of the effective volume, average magnetization and the
correlation function. In Sec. IV we present our analytical and numerical
results for spin-crossover molecular chain. Finally, results and discussions
are given in Sec. V.

##  II The system

During the past two decades molecular magnetism has experienced significant
advancements, and molecular materials have been able not only to display the
different type of magnetic behaviors, initially identified in simpler solids
including ferromagnets, but they have also been able to provide examples of
materials exhibiting novel magnetic phenomena molnar2018spin ;
coronado2020molecular . In order to characterize the behavior of 1D molecular
chain considering both, structural and magnetic properties, we are basing on
the simple conception in which the molecules compound a coordination polymer.

We denote higher energy state as pseudo-spin $+1$ state and lower energy state
as pseudo-spin $-1$ state. Spin $+1$ state has degeneracy $g_{+}$ and spin
$-1$ has degeneracy $g_{-}$. We assume pair interactions depend on the
particles pseudo-states, and the potentials of the corresponding pairs of
pseudo-spin states interactions are $V_{--}(r)$, $V_{+-}(r)$ and $V_{++}(r)$.
Schematic view of the three microscopic configurations of two neighboring
atoms in the chain is shown in the Fig. 1(a). The Hamiltonian of the system is
a sum of the pair potentials and single particle field Linares2004 ;
Boukheddaden2007 ; Lecomte2008Two-variable

$\hat{H}=\sum_{n=1}^{N-1}V_{s_{n}s_{n+1}}(x_{n}-x_{n+1})+\sum_{n=1}^{N}W_{s_{n}},$
(1)

where $N$ is the total number of molecules and $W_{s_{n}}$ is the energy of
the single-molecule pseudo-state. The difference of the pseudo-state energies
$\Delta=W_{+}-W_{-}$ is the external ligand field acting on a single molecule.
We apply an harmonic approximation for the pair nearest-neighbor potential
$V_{s_{n}s_{n+1}}(r)$ at the potential minimum

$V_{s_{n}s_{n+1}}(r)=V_{s_{n}s_{n+1}}^{(0)}+\frac{1}{2}K_{s_{n}s_{n+1}}\left(r-a_{s_{n}s_{n+1}}\right)^{2},$
(2)

where $a_{s_{n}s_{n+1}}$ is a distance where potential has minimum,
$V_{s_{n}s_{n+1}}^{(0)}$ is the potential depth and $K_{s_{n}s_{n+1}}$ is an
elastic constant coupling $n$-th and $(n+1)$-st molecules in the pseudo-states
$s_{n}$ and $s_{n+1}$. Schematic treatment of the potentials is shown in Fig.
1(b). In Eq. (2) the first term cannot be nullified as
$V_{s_{n}s_{n+1}}^{(0)}=V_{s_{n}s_{n+1}}(a_{s_{n}s_{n+1}})$ and it is
different for the $V_{--}(r)$, $V_{+-}(r)$ and $V_{++}(r)$ potentials.

Figure 1: Schematic interactions of the pseudo-spin states and treatment of
the inter-particle potential. (a) All possible configurations of the nearest
pseudo-spin states. Interaction potentials and average distances between
particles depend on the pseudo-spin states. (b) Interaction potential and
harmonic approximation. We consider possible displacement of the particles
from the equilibrium position for the given pseudo-state configuration to be
small.

Let’s introduce variables $q_{n}=x_{n}-x_{n+1}$. We split initial Hamiltonian
(1) to two parts

$\hat{H}=\hat{H}_{1}+\hat{H}_{2},$ (3)

where

$\hat{H}_{1}=\sum_{n=1}^{N-1}V_{s_{n}s_{n+1}}^{(0)}+\sum_{n=1}^{N}W_{s_{n}},$
(4)

and

$\hat{H}_{2}=\frac{1}{2}\sum_{n=1}^{N-1}K_{s_{n}s_{n+1}}\left(q_{n}-a_{s_{n}s_{n+1}}\right)^{2}.$
(5)

We rewrite first part of the Hamiltonian (4) in terms of pseudo-spins

$\hat{H}_{1}=E_{0}\\!-\\!\sum_{n=1}^{N-1}\\!J\\!s_{n}s_{n+1}\\!-\\!\sum_{n=1}^{N-1}\\!B\\!\frac{s_{n}+s_{n+1}}{2}\\!-\\!W(s_{1})\\!-\\!W(s_{N})\\!,$
(6)

where

$E_{0}=\frac{N-1}{4}\left(V_{--}^{(0)}+V_{++}^{(0)}\right)+\frac{N-1}{2}V_{+-}^{(0)}+N\frac{W_{+}+W_{-}}{2},$
(7a)
$J=-\frac{1}{4}\left(V_{--}^{(0)}+V_{++}^{(0)}\right)+\frac{1}{2}V_{+-}^{(0)},$
(7b) $B=\frac{1}{4}\left(V_{++}^{(0)}-V_{--}^{(0)}\right)-\frac{\Delta}{2},$
(7c) $W(s_{n})=-\frac{\Delta}{4}s_{n}.$ (7d)

The spins on the edge of the system experience additional external field
$W(s)$ due to the geometry of the system.

###  II.1 Partition function

The statistical properties of the model are completely determined by the
partition function

$Z\\!=\\!\sum_{\langle s_{1},\ldots,s_{N}\rangle}\\!\iiint dq_{1}\cdots
dq_{N-1}g_{s_{1}}\cdots g_{s_{N}}e^{-\beta
E\left(q_{1},\ldots,q_{N-1},s_{1},\ldots,s_{N}\right)}\\\ =\sum_{\langle
s_{1},\ldots,s_{N}\rangle}g_{s_{1}}\cdots
g_{s_{N}}\left(\prod_{n=1}^{N-1}\sqrt{\frac{2\pi}{\beta
K_{s_{n}s_{n+1}}}}\right)e^{-\beta E_{1}},$ (8)

where the sum is over all states $\langle s_{1},\ldots,s_{N}\rangle$ with
energy $E\left(s_{1},\ldots,s_{N}\right)$, $g_{s_{n}}$ represents the
degeneracy of the state $s_{n}$ and $\beta$ denotes the inverse of the
Boltzmann constant times temperature. It is convenient to consider an ensemble
in which $Z$ depends on the temperature and the field. Certainly, the choice
of boundary conditions becomes irrelevant in the thermodynamic limit.

The coefficients we get during the integration over the phonon degrees of
freedom may be expressed in terms of the spin operators as follows

$\sqrt{\frac{2\pi}{\beta K_{s_{n}s_{n+1}}}}=e^{\delta\epsilon+\delta
b(s_{n}+s_{n+1})/2+\delta js_{n}s_{n+1}},$ (9)

where

$\delta\epsilon=-\frac{1}{8}\ln\left(\frac{\beta^{4}}{(2\pi)^{4}}K_{+-}^{2}K_{--}K_{++}\right),$
(10a) $\delta j=\frac{1}{8}\ln\left(\frac{K_{+-}^{2}}{K_{--}K_{++}}\right),$
(10b) $\delta b=\frac{1}{4}\ln\left(\frac{K_{--}}{K_{++}}\right).$ (10c)

We present spin state degeneracies as follows

$g_{s_{n}}=e^{\frac{1}{2}\left(\ln g_{+}+\ln g_{-}\right)+\frac{1}{2}\left(\ln
g_{+}-\ln g_{-}\right)s_{n}}.$ (11)

We express expression under the exponent of the (8) as a sum of pair terms and
single-particles terms acting only on the boundary. Therefore we have

$Z=\sum_{\langle
s_{1},\ldots,s_{N}\rangle}e^{\epsilon+\sum_{n=1}^{N-1}v(s_{n},s_{n+1})+w(s_{1})+w(s_{N})}$
(12)

where

$w(s_{n})=w\frac{s_{n}}{2},$ (13)

with $w=\frac{1}{2}\ln g-\frac{\beta\Delta}{2}$, and effective two-particle
energy terms

$v(s_{n},s_{n+1})=js_{n}s_{n+1}+b(s_{n}+s_{n+1})/2,$ (14)

and

$\epsilon=-\beta E_{0}+\frac{N}{2}\left(\ln g_{+}+\ln
g_{-}\right)+(N-1)\delta\epsilon,$ (15a) $j=\beta J+\delta j,$ (15b) $b=\beta
B+\frac{1}{2}\ln g+\delta b.$ (15c)

We make notation $g=\frac{g_{+}}{g_{-}}$. As we will see later, it is not the
degeneracy of each energy state that is essential, but the ratio of HS and LS
states degeneracies.

The partition function (12) can be expressed as the partition function of the
Ising model with the effective Hamiltonian

$\hat{H}_{eff}=E_{0,eff}-\sum_{n=1}^{N-1}J_{eff}\hat{s}_{n}\hat{s}_{n+1}-\sum_{n=2}^{N-1}B_{eff}\hat{s}_{n}\\\
+\frac{B_{boundary}+B_{eff}}{2}(\hat{s}_{1}+\hat{s}_{N}),$ (16)

where $E_{0,eff}=E_{0}-\frac{Nk_{B}T}{2}\ln g_{+}g_{-}-(N-1)\delta\epsilon
k_{B}T$, $J_{eff}=J\\!+\\!\delta jk_{B}T$,
$B_{\text{eff}}=B+\frac{k_{B}T}{2}\ln g+\delta bk_{B}T$ and
$B_{boundary}=-\frac{\Delta}{2}+\frac{k_{B}T}{2}\ln g$. The effective
Hamiltonian coincides with the Ising model Hamiltonian with the reference
energy, effective magnetic field and ferromagnetic interaction constant being
functions of temperature. This dependence roots from the taking into account
pseudo-states degeneracy and phononic interactions. Eqs. (7d) and (10c)
establish correspondence between the parameters of the actual inter-molecular
potentials and parameters of the effective Ising model.

The effective temperature dependent magnetic field bousseksou1992ising ;
boukheddaden2000one  $B_{\text{eff}}=B+\frac{k_{B}T}{2}\ln g+\delta bk_{B}T$
acts on the bulk, while molecules on the edges are subject of the external
field $\frac{1}{2}(B_{boundary}+B_{eff})$, reflecting their boundary
distinction. Effective Ising-like systems with the phononic interactions were
introduced in works Linares2004 ; Boukheddaden2007 . In our paper we continue
investigation of this model with lifting down two significant restrictions: a)
all phononic interactions arise from the intermolecular interactions have same
potential depth; b) average distances of the intermolecular potential should
form arithmetic progression. We note, this results in having different form of
the effective ferromagnetic interaction constant comparing to conclusions
Linares2004 ; Boukheddaden2007 . We consider phononic interactions in the
scope of corrections to the parameters of an effective Ising model while
coefficients $\delta\epsilon$, $\delta j$ and $\delta b$ define corrections to
the effective energy reference, ferromagnetic interaction strength and
external field. Accordingly the fictitious spins interact with their nearest
neighbor in the Ising way (parametrized by the phononic interactions). In
addition, we have an external field (B) parametrized by the phononic
interactions and the degeneracy between the two lower states. The spin-
independent part cannot be discarded since it is a function of temperature.
Altogether, our model contains the following parameters: the degeneracy
between the two lower states $g$, the inter-molecular coupling $J$, the
external ligand field $\Delta$ and the external field $B$, the elastic
constants $K_{s_{n}s_{n+1}}$, the temperature $T$.

##  III Transfer-matrix formalism

Thermodynamic properties of the system are completely described with the
partition function. Here we use the transfer matrix formalism Linares2004 ;
Boukheddaden2007  to calculate the partition function.  We rewrite the
partition function (12) as

$Z=e^{\epsilon}\operatorname{Tr}\hat{T}^{N-1}\hat{R},$ (17)

where transfer matrix is

$\hat{T}=e^{v(s_{n},s_{n+1})}=\begin{pmatrix}e^{j+b}&e^{-j}\\\
e^{-j}&e^{j-b}\end{pmatrix},$ (18)

and matrix $\hat{R}$ is accounting effects of the field acting on the surface
spins

$\hat{R}=e^{w(s_{N})+w(s_{1})}=\begin{pmatrix}e^{w}&1\\\
1&e^{-w}\end{pmatrix}.$ (19)

For calculating $\operatorname{Tr}\hat{T}^{N-1}\hat{R}$ we should make a
rotation of the basis to one where the transfer matrix is diagonal

$Z=e^{\epsilon}\operatorname{Tr}\hat{U}\hat{U}^{-1}\hat{T}^{N-1}\hat{U}\hat{U}^{-1}\hat{R},$
(20)

where the rotation matrix

$\hat{U}=\begin{pmatrix}\cos\phi&\sin\phi\\\ -\sin\phi&\cos\phi\end{pmatrix},$
(21)

and the rotation angle $\phi$ is given by the equation

$\cot 2\phi=e^{2j}\sinh(b).$ (22)

The eigenvalues of the transfer matrix are

$\lambda_{\pm}=\left(e^{j}\cosh b\pm\sqrt{e^{2j}\sinh^{2}b+e^{-2j}}\right).$
(23)

Therefore, we obtain the partition function

$Z=e^{\epsilon}\left(c_{+}\lambda_{+}^{N-1}+c_{-}\lambda_{-}^{N-1}\right),$
(24)

where the coefficients are

$c_{+}=\cosh w+\frac{e^{-2j}+\sinh b\sinh w}{\sqrt{\sinh^{2}b+e^{-4j}}},$
(25a) $c_{-}=\cosh w-\frac{e^{-2j}+\sinh b\sinh
w}{\sqrt{\sinh^{2}b+e^{-4j}}}.$ (25b)

The free energy density is given by

$f=-\frac{1}{\beta N}\ln Z.$ (26)

In the thermodynamic limit we have the desired result:

$f=-\lim_{N\to\infty}\frac{1}{\beta N}\ln
Z=-\frac{\epsilon}{N\beta}-\frac{1}{\beta}\ln\lambda_{+}.$ (27)

Average fictitious magnetization per spin can be calculated directly by the
formula

$\left\langle s\right\rangle=\frac{1}{N}\sum_{n=1}^{N}\left\langle
s_{n}\right\rangle=\frac{1}{N}\sum_{n=1}^{N}\frac{1}{Z}\mathrm{Tr}\hat{s}_{n}e^{-\beta\hat{H}}=\frac{1}{N}\frac{\partial\ln
Z}{\partial b}.$ (28)

The magnetization per spin $m=\left\langle s\right\rangle$ at nonzero $T$ and
$B$:

$m=\frac{\sinh(b)}{\sqrt{\sinh^{2}b+e^{-4j}}}.$ (29)

Results for symmetric degeneracies case $g_{+}=g_{-}$ and without phononic
part repeat well-known behavior of conventional Ising model.

###  III.1 Volume of the system and the correlation function

Let’s calculate average effective volume of the molecular chain

$L=\sum_{n=1}^{N-1}\langle x_{n+1}-x_{n}\rangle=\sum_{n=1}^{N-1}\langle
q_{n}\rangle.$ (30)

By the definition

$\langle q_{n}\rangle\\!=\\!\frac{1}{Z}\sum_{\langle
s_{1},\ldots,s_{N}\rangle}\\!\iiint dq_{1}\cdots dq_{N-1}q_{n}g_{s_{1}}\cdots
g_{s_{N}}e^{-\beta E}$ (31)

Integrating over the phonon degrees of freedom we get

$\langle q_{n}\rangle=\frac{\sum_{\langle
s_{1},\ldots,s_{N}\rangle}a_{s_{n}s_{n+1}}e^{-\beta H_{eff}}}{\sum_{\langle
s_{1},\ldots,s_{N}\rangle}e^{-\beta H_{eff}}}.$ (32)

Thus, the volume of the system $L=\sum_{n=1}^{N-1}\langle
a_{s_{n}s_{n+1}}\rangle$. We rewrite later expression as follows

$L=\sum_{n=1}^{N-1}\left(a_{\epsilon}+a_{J}\langle
s_{n}s_{n+1}\rangle+\frac{a_{B}}{2}\langle s_{n}+s_{n+1}\rangle\right),$ (33)

where $a_{\epsilon}=\frac{1}{4}\left(a_{--}+a_{++}\right)+\frac{1}{2}a_{+-}$,
$a_{J}=\frac{1}{4}\left(a_{--}+a_{++}\right)-\frac{1}{2}a_{+-}$ and
$a_{B}=\frac{1}{2}\left(a_{++}-a_{--}\right)$. Thus, the effective volume of
system is connected with the average magnetization and the correlation
function. Average fictitious spin at the position $n$ corresponds to the local
magnetization

$\langle\hat{s}_{n}\rangle=\frac{\operatorname{Tr}\hat{T}^{n-1}\hat{\sigma}_{z}\hat{T}^{N-n}\hat{R}}{\operatorname{Tr}\hat{T}^{N-1}\hat{R}}.$
(34)

For calculating this expression we go to the eigenbasis of the transfer
matrix. All matrices except of $\hat{U}^{-1}\hat{\sigma}_{z}\hat{U}$ already
have been calculated above. Therefore

$\hat{U}^{-1}\hat{\sigma}_{z}\hat{U}=\begin{pmatrix}-m&-1-m\\\
-1+m&m\end{pmatrix},$ (35)

where $m$ is the expression on the right side of the Eq. (29). We get local
magnetization

$\langle\hat{s}_{n}\rangle=m+\frac{C_{+-}e^{-\frac{n-1}{\xi}}+C_{-+}e^{-\frac{N-n}{\xi}}}{{c_{+}+c_{-}e^{-\frac{N-1}{\xi}}}},$
(36)

where coefficients are

$C_{-+}=C_{+-}=(m^{2}-1)(-\sinh w+e^{2j}\sinh b),$ (37)

and the correlation length is $\xi=-\ln\frac{\lambda_{-}}{\lambda_{+}}$. It is
easy to see that since $\lambda_{-}<\lambda_{+}$, $\xi>0$. The average
magnetization is

$\langle\hat{s}\rangle=m+\frac{C_{+-}+C_{-+}}{N\left(1-e^{-\frac{N}{\xi}}\right)\left(c_{+}+c_{-}e^{-\frac{N-1}{\xi}}\right)}.$
(38)

In the thermodynamic limit average magnetization (38) goes to magnetization
given by the Eq. (29). We note that only average over all spins magnetization
coincides with the classic Ising model result, while average of the individual
spin distinct from the classic result due to the system boundary. We see
boundary effects do not vanish even in the thermodynamic limit.

The local correlation function $G_{n}(r)$ by definition is

$G_{n}(r)=\langle\hat{s}_{n}\hat{s}_{n+r}\rangle=\frac{\operatorname{Tr}\hat{T}^{n-1}\hat{\sigma}_{z}\hat{T}^{r}\hat{\sigma}_{z}\hat{T}^{N-n-r}\hat{R}}{\operatorname{Tr}\hat{T}^{N-1}\hat{R}}.$
(39)

Therefore

$\langle\hat{s}_{n}\hat{s}_{n+r}\rangle=m^{2}+(1-m^{2})\frac{c_{+}e^{-\frac{r}{\xi}}+c_{-}e^{-\frac{N-r-1}{\xi}}}{c_{+}+c_{-}e^{-\frac{N-1}{\xi}}}\\\
+mC_{+-}\frac{e^{-\frac{n-1}{\xi}}-e^{-\frac{n-1+r}{\xi}}+e^{-\frac{N-n-r}{\xi}}-e^{-\frac{N-n}{\xi}}}{c_{+}+c_{-}e^{-\frac{N-1}{\xi}}}.$
(40)

We get the correlation function in the thermodynamic limit

$G(r)=\frac{1}{N}\sum_{n=1}^{N-r-1}\langle\hat{s}_{n}\hat{s}_{n+r}\rangle=m^{2}+(1-m^{2})e^{-\frac{r}{\xi}}.$
(41)

The average magnetization given by the Eq. (36) and the correlation function
given by (40) are exact. We see the average correlation function matches with
the classic Ising model result bellucci2013correlation  in the thermodynamic
limit. Local correlation function (see Eq. (40)) has information about the
edges of the system even in the thermodynamic limit.

Finally, we get expression

$L=(N-1)(a_{\epsilon}+a_{B}m+a_{J}(m^{2}+(1-m^{2})e^{-\frac{1}{\xi}}))\\\
+a_{B}\frac{C_{+-}+C_{-+}}{{c_{+}+c_{-}e^{-\frac{N-1}{\xi}}}}\left[\frac{1}{1-e^{-\frac{N}{\xi}}}-\frac{1}{2}(1+e^{-\frac{N-1}{\xi}})\right]\\\
+a_{J}m\frac{C_{+-}(1-e^{-\frac{1}{\xi}})}{c_{+}+c_{-}e^{-\frac{N-1}{\xi}}}\left[\frac{1}{1-e^{-\frac{N-1}{\xi}}}+\frac{e^{-\frac{N-1}{\xi}}}{1-e^{\frac{N-1}{\xi}}}\right].$
(42)

Expression (42) is exact and defines density $\rho^{-1}=L/N\to
a_{\epsilon}+a_{B}m+a_{J}(m^{2}+(1-m^{2})e^{-\frac{1}{\xi}})$ in the
thermodynamic limit.

##  IV Spin-crossover chains

Here, we consider a general case where degeneracies of the pseudo-spin states
are not equal $g_{+}\neq g_{-}$ and phonon corrections (10c) are non-zero.
This situation is typical for a number of molecular crystals, including spin-
crossover materials. Note that for spin-crossover molecular states, the
relative degeneracy of the high-energy state in comparison with the low-energy
one is large.

We introduce the equilibrium temperature $T_{eq}$ as a temperature when
pseudo-spin states have equal occupations. This occurs when effective field
vanishes $b=0$. Therefore we write

$T_{eq}=-\frac{B}{k_{B}\left(\frac{1}{2}\ln g+\delta b\right)}.$ (43)

We note that for certain values of the external field $B$ and pseudo-spin
degeneracies $g$, temperature $T_{eq}$ can be negative what means that for the
given field and degeneracies there is no such temperature that pseudo-spin
states would have equal occupations.

Figure 2: Average magnetization $m(T,T_{eq})$ for $\frac{1}{2}\ln g+\delta
b=-0.5$ and $\delta j=0$. Dashed line indicates condition $T=T_{eq}$.

Figure 3: (a)-(c) Average magnetization as a function of temperature for
$T_{eq}<T_{crossover}$, $T_{eq}=T_{crossover}$ and $T_{eq}>T_{crossover}$ and
various values of the $\frac{1}{2}\ln g+\delta b=0.1,0.5,1,2$ and $\delta
j=0.35$. In the case $T_{eq}<T_{crossover}$ average magnetization has maximum
at $T=T_{0}$. For $T_{eq}=T_{crossover}$, maximum is reached at
$T_{0}=\infty$. When $T_{eq}>T_{crossover}$ average magnetization is a
monotonous function of temperature and has no extrema.

A spin-crossover material may abruptly change its macroscopic behavior, for
example magnetization, when the external conditions, such as temperature, are
varied. The region where this happens is called abrupt crossover, and it marks
a transition from one state to another. Dependence of average magnetization
per spin as a function of temperature and equilibrium temperature
$m\left(k_{B}T/J,k_{B}T_{eq}/J\right)$ for the fixed value $\frac{1}{2}\ln
g+\delta b=0.5$ is illustrated graphically in Fig. 2. The dashed line
corresponds to the condition $T=T_{eq}$, or equivalently $B_{eff}=0$. In this
case, the magnetization per spin $m=0$. The average magnetization is positive
at any temperature when $T_{eq}$ is negative. For a positive $T_{eq}$, $m>0$
when $T>T_{eq}$, and $m<0$ when $T<T_{eq}$, respectively. We see that for a
large $T_{eq}$, transition from $m=-1$ to $m=1$ is wide and smooth while at a
small $T_{eq}$ this transition happens in the small region close to the point
$T=T_{eq}$. Such behaviour of the average magnetization in the limit
$T_{eq}\to 0$ gives discontinuous behavior, and we have the phase transition
at $T=0$. Therefore we talk about the crossover from smooth magnetization
behaviour to jump-like.

We identify the crossover temperature in the following way. At zero
temperature the average magnetization is defined only by the sign of the
effective magnetic field $m(T=0)=\operatorname{sign}(b)$, while in the high
$T$ regime the average magnetization goes to value

$m(T\to\infty)=\frac{\sinh\left(\frac{1}{2}\ln g+\delta
b\right)}{\sqrt{\sinh^{2}\left(\frac{1}{2}\ln g+\delta b\right)+e^{-4\delta
j}}}.$ (44)

In the vicinity of $T=T_{eq}$, symmetry $\hat{s}\to-\hat{s}$ with $B_{eff}\to-
B_{eff}$ exists and thus the average magnetization is an anti-symmetric
function of $t=(T-T_{eq})/T_{eq}$. Consequently for a small $T_{eq}$ the
average magnetization changes abruptly from $m=-1$ to nearly $1$ and then goes
to the high temperature asymptotic, while for large values of the $T_{eq}$
transition from $m=-1$ to the high $T$ asymptotic is smooth and has no maximum
in-between. Therefore we identify the abrupt regime by existence of the local
maximum of the average magnetization, and the temperature $T_{0}$ as the
temperature at which the average magnetization has maximum. The gradual regime
we identify as a regime with no extrema in the average magnetization. Thus
condition for the existence of abrupt change in the average magnetization is

$\left.\frac{\partial m}{\partial T}\right|_{T=T_{0}}=0.$ (45)

Maximal equilibrium temperature $T_{eq}$ at which Eq. (45) has finite
solutions for the $T_{0}$ is the crossover temperature. The derivative is
always positive and the magnetization $m$ is a monotonous function of
temperature when $T_{eq}<T_{crossover}$. For identifying the crossover
temperature exactly we write

$\operatorname{sign}\left(\frac{\partial m}{\partial
T}\right)=\operatorname{sign}\left(\frac{k_{B}T_{eq}}{2J}\left(\frac{1}{2}\ln
g+\delta b\right)\right.\\\ \left.-\tanh\left(\left(\frac{1}{2}\ln g+\delta
b\right)\left(1-\frac{T_{eq}}{T}\right)\right)\right).$ (46)

Abrupt crossover doesn’t occur when $\frac{\partial m}{\partial T}$ has same
sign at any temperature. This might happen only when function
$\tanh\left(\left(\frac{1}{2}\ln g+\delta
b\right)\left(1-\frac{T_{eq}}{T}\right)\right)$ is smaller than
$\frac{T_{eq}\left(\frac{1}{2}\ln g+\delta b\right)}{T_{c}(1-\delta j)}$ at
any finite temperature and equal to it at the infinite temperature. Thus we
get the crossover temperature

$T_{crossover}=\frac{2J}{k_{B}}\frac{\tanh\left(\frac{1}{2}\ln g+\delta
b\right)}{\frac{1}{2}\ln g+\delta b}.$ (47)

We note that the function $\tanh(x)/x$ for all $x$ is smaller than one, and
thus $T_{crossover}\leq\frac{2J}{k_{B}}$. The maximal crossover temperature
$T_{crossover}=\frac{2J}{k_{B}}$ is reached when $\frac{1}{2}\ln g+\delta b\to
0$. In Fig. 3(a-c) dependence of the average magnetization on temperature for
fixed $T_{eq}$ is plotted. Three different cases are considered: (a)
$T_{eq}<T_{crossover}$, (b) $T_{eq}=T_{crossover}$ and (c)
$T_{eq}>T_{crossover}$.  At high temperatures the magnetization goes to the
asymptotic value given by the Eq. (44). At all panels (a-c) evolution of the
magnetization is shown for four values of the $\frac{1}{2}\ln g+\delta
b=0.1,0.5,1,2$ and $\delta j=0.3$ being constant. Corresponding curves of the
same color from different panels share zero- and high temperature asymptotes.
The high temperature asymptotic value of magnetization $m(T\to\infty)$ goes to
$0$ when $\frac{1}{2}\ln g+\delta b\to 0$ . Asymptotic value $m=1$ is reached
in the case of the large difference in the pseudo-state degeneracies ($g\gg
1$).

Condition defining maximal magnetization is given by the Eq. (45). Therefore
we get

$T_{0}=\frac{T_{eq}}{1-\frac{1}{\frac{1}{2}\ln g+\delta
b}\operatorname{arctanh}\left(\frac{k_{B}T_{eq}}{2J}\left(\frac{1}{2}\ln
g+\delta b\right)\right)}.$ (48)

Figure 4: $T_{0}$ as a function of ratio $T_{eq}$ for various values of the
expression $\frac{1}{2}\ln g+\delta b=0.1,0.5,1$ and $2$ (right to left: red,
green, orange, blue). Vertical dashed lines are asymptotes of the
corresponding $T_{crossover}$ given by Eq. (47).

Dependence of $T_{0}$ from the $T_{eq}$ for various values of $\frac{1}{2}\ln
g+\delta b$ are plotted in Fig. 4. In the limit $T_{eq}\to 0$, $T_{0}$ goes to
zero as well. With increase of the $T_{eq}$ goes to infinity as $T_{eq}\to
T_{crossover}$. In the region over the line defined by the Eq. (48) the
derivative $\frac{\partial m}{\partial T}>0$ and in the region below
$\frac{\partial m}{\partial T}<0$. Eq. (48) can be rewritten explicitly in
terms of the effective Hamiltonian parameters as following

$T_{0}=-\frac{B}{k_{B}}\frac{1}{\frac{1}{2}\ln g+\delta
b+\operatorname{arctanh}\frac{B}{J}}.$ (49)

It is interesting to note that the average magnetization does depend on the
phonon correction $\delta j$, while $T_{0}$ and $T_{crossover}$ do not depend
on it.

Figure 5: Average magnetization as a function of temperature for
$T_{eq}=0.3T_{crossvoer}$ and various values of $\delta
j=0.5,0.1,-0.1,-0.5,-0.9$. Maximum of magnetization is at $T=T_{0}$.

In Fig. 5 dependence of the average magnetization on temperature for fixed
$T_{eq}=0.3T_{crossvoer}$ is depicted. The phonon interaction constant $\delta
j$ was chosen from the set $0.5,0.1,-0.1,-0.5,-0.9$. With the increase of
$\delta j$ we observe expected change of the high temperature asymptotes. This
results in the difference of the magnetization change magnitude near the
$T=T_{eq}$, but location of all characteristic temperatures remains the same.

The internal energy is $U=-\frac{\partial}{\partial\beta}\ln Z$

$U=E_{0}+N\left(\frac{1}{2}k_{B}T-B\langle s\rangle-JG(1)\right),$ (50)

where the average local magnetization and the correlation function are given
by the Eqs. (36) and (40). In the thermodynamic limit these expressions reduce
to the average magnetization $\langle s\rangle=m$ and the correlation function
given by Eq. (41).

Figure 6: Specific heat capacity per particle $c_{P}$ (solid line) as a
function of temperature and the $-B\frac{\partial m}{\partial T}$ derivative
(dashed line). (a) $T_{eq}=3T_{crossover}$. Parameters are chosen to be same
as Fig. 3(c). Maximum of the heat capacity is shifted from the maximum of the
derivative $\frac{\partial m}{\partial T}$. (b) $T_{eq}=0.3T_{crossover}$.
Parameters are chosen to be same as in Fig. 3(a).

Therefore the specific heat capacity per particle at constant pressure in the
thermodynamic limit is

$c_{P}=\frac{1}{2}k_{B}-B\frac{\partial m}{\partial
T}-J\frac{\partial}{\partial
T}\left(m^{2}+(1-m^{2})e^{-\frac{1}{\xi}}\right).$ (51)

In Fig. 6 the specific heat capacity is plotted for two cases: (a)
$T_{eq}=3T_{crossover}$ and (b) $T_{eq}=0.3T_{crossover}$. These cases
correspond to Figs. 3(c) and (a). Other parameters are fixed and chosen to be
the same as for the orange curves in Fig. 3. At the left panel we observe the
Schottky anomaly which transforms into large narrow peak at the right panel
with increase of intermolecular interaction stiffness. In the models with
$J=0$ maximum of the heat capacity corresponds exactly to the maximum of
function $\frac{\partial m}{\partial T}$ Boukheddaden2007 . We note that in
our model the location of the maximum of heat capacity is different from the
maximum of the function $\frac{\partial m}{\partial T}$. This difference
originates from the non-zero ferromagnetic interaction $J$ and is defined by
the correlation function. Another important effect we observe is the existence
of the second maximum of the heat capacity when the crossover is abrupt. In
this case second peak has much smaller high and broader comparing to the
narrow peak around $T_{0}$. This effect is novel for 1D systems and was not
present in similar 1D Ising-like models.

##  V Summary and conclusions

The recent progress of molecular magnetism towards low dimensions, with an
increasing focus on applications, led to our research has been accented on
peculiarity of 1d spin crossover compound. We have systematically studied the
one-dimensional molecular chain with degenerate states and phononic
interaction. Exact solutions by the method of transfer matrix modified for
free boundary were obtained. We have shown that there is an exact mapping
between the elastic molecular chain and Ising models of SCO materials. We have
shown this explicitly for a one-dimensional model, but our approach can be
extended to other elastic models and geometries. We carefully consider the
possible problems involving finite size effects and boundary conditions on the
properties of the model which might be investigated in molecular chain.

An effective Ising-like Hamiltonian corresponding to the exact partition
function of the system is obtained. Moreover, this Hamiltonian and its
structure have physical consequences, which implies that the effective
Hamiltonian is much more than just a convenient mathematical construct. The
pivotal role in the transition in low-dimensional molecular structure can be
physically attributed to the relative degeneracy of states. The significant
volume change in the cores molecules (SCO-molecules in the present case) leads
to essential renormalization of the correlation function. We have analyzed the
regimes of the LS–HS crossover and identified if the crossover is abrupt or
gradual for the specific parameters. In the case of abrupt crossover we have
showed possibility of existing two peak thermal dependence of the specific
heat capacity.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

## References

  * [1] R Traiche, M Sy, and K Boukheddaden.  Elastic frustration in 1d spin-crossover chains: Evidence of multi-step transitions and self-organizations of the spin states.  The Journal of Physical Chemistry C, 122(7):4083–4096, 2018.
  * [2] G. Mohan Kumar.  Nanotechnology: Nanomaterials and Nanodevices.  Narosa Publishing House, New Delhi, 2016.
  * [3] Kuppusamy Senthil Kumar and Mario Ruben.  Emerging trends in spin crossover (sco) based functional materials and devices.  Coordination Chemistry Reviews, 346:176–205, 2017.
  * [4] H. W. Huang.  Asymmetrical Ising model.  Phys. Rev. B, 12:216–218, 1975.
  * [5] R. Harris.  A study of first- and second-order phase transitions using Monte Carlo simulations in the micro-canonical ensemble.  Physics Letters A, 111(6):299 – 303, 1985.
  * [6] A. Bousseksou, J. Nasser, J. Linares, K. Boukheddaden, and F. Varret.  Ising-like model for the two-step spin-crossover.  J. Phys. I France, 2(7):1381–1403, 1992.
  * [7] K Boukheddaden, J Linares, H Spiering, and F Varret.  One-dimensional ising-like systems: an analytical investigation of the static and dynamic properties, applied to spin-crossover relaxation.  The European Physical Journal B-Condensed Matter and Complex Systems, 15(2):317–326, 2000.
  * [8] D. Schebarchov, T. P. Schulze, and S. C. Hendy.  Degenerate Ising model for atomistic simulation of crystal-melt interfaces.  The Journal of Chemical Physics, 140(7):074704, 2014.
  * [9] J. Sienkiewicz, K. Suchecki, and J. A Hołyst.  Finite size induces crossover temperature in growing spin chains.  Physical Review E, 89(1):012105, 2014.
  * [10] S. Cajahuaringa and A. Antonelli.  Nonequilibrium free energy methods applied to magnetic systems: the degenerate ising model.  Journal of Statistical Physics, 175(5):1006–1021, 2019.
  * [11] JA Nasser.  Diluted spin conversion compounds behaviours in the atom-phonon coupling model: case of not too large dilution.  The European Physical Journal B-Condensed Matter and Complex Systems, 48(1):19–27, 2005.
  * [12] Y Konishi, H Tokoro, M Nishino, and S Miyashita.  Monte carlo simulation of pressure-induced phase transitions in spin-crossover materials.  Physical Review Letters, 100(6):067206, 2008.
  * [13] W. Nicolazzi, S. Pillet, and C. Lecomte.  Two-variable anharmonic model for spin-crossover solids: A like-spin domains interpretation.  Phys. Rev. B, 78:174401, Nov 2008.
  * [14] Cr Enachescu, L Stoleriu, Al Stancu, and A Hauser.  Model for elastic relaxation phenomena in finite 2d hexagonal molecular lattices.  Physical review letters, 102(25):257204, 2009.
  * [15] Hong-Zhou Ye, Chong Sun, and Hong Jiang.  Monte-carlo simulations of spin-crossover phenomena based on a vibronic ising-like model with realistic parameters.  Physical Chemistry Chemical Physics, 17(10):6801–6808, 2015.
  * [16] J Linares, H Spiering, and F Varret.  Analytical solution of 1d ising-like systems modified by weak long range interaction.  The European Physical Journal B-Condensed Matter and Complex Systems, 10(2):271–275, 1999.
  * [17] K Boukheddaden, J Linares, R Tanasa, and C Chong.  Theoretical investigations on an axial next nearest neighbour ising-like model for spin crossover solids: one- and two-step spin transitions.  Journal of Physics: Condensed Matter, 19(10):106201, feb 2007.
  * [18] K Boukheddaden, M Nishino, and S Miyashita.  Molecular dynamics and transfer integral investigations of an elastic anharmonic model for phonon-induced spin crossover.  Physical review letters, 100(17):177206, 2008.
  * [19] VA Zagrebnov and BK Fedyanin.  Spin-phonon interaction in the ising model.  Theoretical and Mathematical Physics, 10(1):84–93, 1972.
  * [20] S R Salinas.  On the one-dimensional compressible ising model.  Journal of Physics A: Mathematical, Nuclear and General, 6(10):1527–1533, oct 1973.
  * [21] V B Henriques and S R Salinas.  Effective spin hamiltonians for compressible ising models.  Journal of Physics C: Solid State Physics, 20(16):2415–2429, jun 1987.
  * [22] A Lehmann-Szweykowska, M Kurzyński, R Wojciechowski, M Wiesner, and B Mróz.  Anomalous phase transitions in licsso 4 in the compressible hcp ising model.  Acta Physica Polonica, A., 121, 2012.
  * [23] JA Nasser, K Boukheddaden, and J Linares.  Two-step spin conversion and other effects in the atom-phonon coupling model.  The European Physical Journal B-Condensed Matter and Complex Systems, 39(2):219–227, 2004.
  * [24] K. Boukheddaden, S. Miyashita, and M. Nishino.  Elastic interaction among transition metals in one-dimensional spin-crossover solids.  Phys. Rev. B, 75:094112, 2007.
  * [25] W Nicolazzi, J Pavlik, S Bedoui, G Molnár, and A Bousseksou.  Elastic ising-like model for the nucleation and domain formation in spin crossover molecular solids.  The European Physical Journal Special Topics, 222(5):1137–1159, 2013.
  * [26] T D Oke, F Hontinfinde, and K Boukheddaden.  Bethe lattice approach and relaxation dynamics study of spin-crossover materials.  Applied Physics A, 120(1):309–320, 2015.
  * [27] K. Nebbali, C. D. Mekuimemba, C. Charles, S. Yefsah, G. Chastanet, A. J. Mota, E. Colacio, and S. Triki.  One-dimensional thiocyanato-bridged fe (ii) spin crossover cooperative polymer with unusual fen5s coordination sphere.  Inorganic chemistry, 57(19):12338–12346, 2018.
  * [28] W. Lan, F. J. Valverde-Muñoz, Y. Dou, X. Hao, M. C. Muñoz, Z. Zhou, H. Liu, Q. Liu, J. A. Real, and D. Zhang.  A thermal-and light-induced switchable one-dimensional rare loop-like spin crossover coordination polymer.  Dalton Transactions, 48(45):17014–17021, 2019.
  * [29] K. A. Vinogradova, D. P. Pishchur, V. Y. Komarov, L. G. Lavrenova, and M. B. Bushuev.  Cooperative spin transition in a 1d-polymeric complex [fe (4-ethyl-1, 2, 4-triazole) 3] sif6· nh2o.  Inorganica Chimica Acta, page 119560, 2020.
  * [30] K. V. Raman, A. M. Kamerbeek, A. Mukherjee, N. Atodiresei, T. K. Sen, P. Lazić, V. Caciuc, R. Michel, D. Stalke, S. K. Mandal, et al.  Interface-engineered templates for molecular spin memory devices.  Nature, 493(7433):509–513, 2013.
  * [31] T. Matsumoto, G. N. Newton, T. Shiga, S. Hayami, Y. Matsui, H. Okamoto, R. Kumai, Y. Murakami, and H. Oshio.  Programmable spin-state switching in a mixed-valence spin-crossover iron grid.  Nature Communications, 5(1):1–8, 2014.
  * [32] K. Bairagi, O. Iasco, A. Bellec, A. Kartsev, D. Li, J. Lagoute, C. Chacon, Y. Girard, S. Rousset, F. Miserque, et al.  Molecular-scale dynamics of light-induced spin cross-over in a two-dimensional layer.  Nature Communications, 7(1):1–7, 2016.
  * [33] B. R. Mullaney, L. Goux-Capes, D. J. Price, G. Chastanet, J.-F. Létard, and C. J. Kepert.  Spin crossover-induced colossal positive and negative thermal expansion in a nanoporous coordination framework material.  Nature Communications, 8(1):1–6, 2017.
  * [34] L. Kipgen, M. Bernien, S. Ossinger, F. Nickel, A. J. Britton, L. M. Arruda, H. Naggert, C. Luo, C. Lotze, H. Ryll, et al.  Evolution of cooperativity in the spin transition of an iron (ii) complex on a graphite surface.  Nature Communications, 9(1):1–8, 2018.
  * [35] K. Ridier, A.-C. Bas, Y. Zhang, L. Routaboul, L. Salmon, G. Molnár, C. Bergaud, and A. Bousseksou.  Unprecedented switching endurance affords for high-resolution surface temperature mapping using a spin-crossover film.  Nature Communications, 11(1):1–9, 2020.
  * [36] A. Köbke, F. Gutzeit, F. Röhricht, A. Schlimm, J. Grunwald, F. Tuczek, M. Studniarek, D. Longo, F. Choueikani, E. Otero, et al.  Reversible coordination-induced spin-state switching in complexes on metal surfaces.  Nature Nanotechnology, 15(1):18–21, 2020.
  * [37] M. A. Halcrow, editor.  Spin-crossover materials: properties and applications.  John Wiley & Sons, 2013.
  * [38] Iu. Gudyma, C. Enachescu, and A. Maksymov.  Kinetics of nonequilibrium transition in spin-crossover compounds.  In Nanocomposites, Nanophotonics, Nanobiotechnology, and Applications, pages 375–401. Springer, 2015.
  * [39] T. Yokoyama, Y. Murakami, M. Kiguchi, T. Komatsu, and N. Kojima.  Spin-crossover phase transition of a chain Fe(II) complex studied by x-ray-absorption fine-structure spectroscopy.  Phys. Rev. B, 58:14238–14244, 1998.
  * [40] G. Molnár, S. Rat, L. Salmon, W. Nicolazzi, and A. Bousseksou.  Spin crossover nanomaterials: from fundamental concepts to devices.  Advanced Materials, 30(5):1703862, 2018.
  * [41] E. Coronado.  Molecular magnetism: from chemical design to spin control in molecules, materials and devices.  Nature Reviews Materials, 5(2):87–104, 2020.
  * [42] S. Bellucci and V. Ohanyan.  Correlation functions in one-dimensional spin lattices with Ising and heisenberg bonds.  The European Physical Journal B, 86(11):446, 2013.