# Fileset

[EurJIC_2024_e202400656.pdf](https://mdr.nims.go.jp/filesets/6fc3f4a5-7073-47d1-8c9c-e9b2fdcce78a/download)

## Creator

[Seiji Miyashita](https://orcid.org/0000-0003-0681-3910), [Masamichi Nishino](https://orcid.org/0000-0002-2060-2303), [Cristian Enachescu](https://orcid.org/0000-0001-5173-5451)

## Rights

This is the peer reviewed version of the following article: S. Miyashita, M. Nishino, C. Enachescu, Eur. J. Inorg. Chem. 2024, e202400656., which has been published in final form at https://doi.org/10.1002/ejic.202400656. This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions for Use of Self-Archived Versions. This article may not be enhanced, enriched or otherwise transformed into a derivative work, without express permission from Wiley or by statutory rights under applicable legislation. Copyright notices must not be removed, obscured or modified. The article must be linked to Wiley’s version of record on Wiley Online Library and any embedding, framing or otherwise making available the article or pages thereof by third parties from platforms, services and websites other than Wiley Online Library must be prohibited.[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[Mechanisms of Separation of the Elastic Step and the Thermal Step](https://mdr.nims.go.jp/datasets/76c1e007-8606-4b31-a634-9f3b85dddaee)

## Fulltext

Mechanisms of separation of the elastic step and the thermal stepSeiji Miyashita∗JSR-UTokyo Collaboration Hub, CURIE, Department of Physics,Graduate School of Science, University of Tokyo, Bunkyo-Ku, Tokyo, 113-0033, JapanMasamichi Nishino†Research Center for Materials Nanoarchitectonics,National Institute for Materials Science, 1-1 Namiki,Tsukuba, Ibaraki 305-0044, JapanCristian Enachescu‡Faculty of Physics, Alexandru Ioan Cuza University of Iasi, Iasi 700506, Romania(Dated: December 9, 2024)The elastic interaction causes a two-step conversion from the low-spin (LS) to the high-spin state(HS) following a femtosecond laser pulse excitation. These steps occur on different timescales: aspin conversion (elastic step) in a short time due to pressure propagation and a late spin conversion(thermal step) resulting from heat diffusion. In this paper, we provide a brief review of modelsfor spin crossover phenomena, emphasizing the role of the difference of degeneracy of the HS stateand LS state and elastic interaction due to changes in local structures. We also review theoreticalapproaches to these effects. Additionally, to analyze the dynamical process of spin conversion afterphoto-irradiation, it is necessary to consider mechanisms of heat transfer among spins. Recently,we proposed a two-temperature model that incorporates both lattice and spin temperatures, whichenables us to reproduce the two steps under an early uniformization of the lattice temperature,highlighting the importance of the timescales of these processes. Furthermore, we explore possiblemechanisms to separate the two steps by examining the different timescales of processes of spinconversions by mechanical and entropic forces, and also by considering changes of local temperaturedue to spin conversion.I. INTRODUCTIONSpin-crossover (SC) phenomenon is a transition between bistable states, typically the high spin (HS) state and thelow spin (LS) states [1]. Due to the difference in degeneracies of the LS and HS states, various interesting thermalprocesses emerge, including different forms of hysteresis. This phenomena have been analyzed by the so called Isinglike model [2], where the effects of degeneracy are incorporated by introducing a temperature-dependent external fieldin the usual non-degenerate Ising model. The metastability of the phases results in a range of intriguing orderingprocesses [3–5].Apart from the degeneracy difference, it has been noted out that the volume difference of the two states, or moregenerally difference of local lattice structure, causes deformation of the lattice. The elastic energy resulting from thisdeformation induces an effective interaction among the spin states [6–9]. This elastic interaction is responsible for thelong range ferromagnetic interaction [10] which causes various interesting phenomena [11–18].The elastic interaction plays a significant role not only in static thermodynamic properties but also in dynamicsin ordering process. Recently, an intriguing dynamical ordering process following ultrafast photo-irradiation drivenby the elastic interaction was observed [19]. Specifically, the high pressure due to the elastic force at the excited HSpropagates in a short period corresponding to sound propagation, and expands the lattice, which causes an increase ofHS fraction. This phenomenon is known as the ”elastic step.” This step occurs much before further increase of the HSfraction due to heating of the whole system due to photo-irradiation, referred to as the ”thermal step.” Consequently,two peaks appear in the HS fraction’s evolution. This dynamical phenomenon is driven by propagation of elasticdeformation and thermal diffusion [19, 20]. The appearance of these two peaks depends on the sample size. In a smallsystem, the injected energy escapes to the thermal bath (sample environment) before the thermal step occurs, so onlythe elastic step is observed. while In a large sample, however both steps are observed.∗miyashita-seiji@g.ecc.u-tokyo.ac.jp@keio.jp†NISHINO.Masamichi@nims.go.jp‡cristian.enachescu@uaic.ro2One may attribute such separation to the properties of fast propagation of the elastic deformation and slow prop-agation of heat. To study this process of spin conversion with thermal diffusion, models of dynamic processes havebeen explored [20, 21]. Assuming that thermal diffusion is much slower compared to the pressure propagation, thetemperature profile in the system becomes highly inhomogeneous in the process, which can cause the two steps indevelopment of HS fraction. This model is called “single-temperature model”. However, it has reported that thethermalization of the lattice is very fast [22] and uniformization of the temperature happens in a short time before thethermal step. In this context, a two-temperature model has been introduced, considering two types of temperatures:the lattice temperature and the spin temperature [21]. The former represents the energy of the total lattice vibration,while the latter corresponds to an energy responsible to the spin conversion.To study these dynamical properties, detailed timescales of processes are important. In the present paper, we beginwith a brief review on modeling of the spin crossover phenomena and explore modeling of the dynamics of processes.First, we consider the timescales of relaxation in both of the original degenerate model and the Ising-like model,finding that both yield the same thermodynamic properties. Next, we consider dynamical features of processes, suchas the difference between the spin conversion process driven by a mechanical force (pressure) and that due to anentropic effect (difference of degeneracy). We also assess the effect of heat emission during the spin conversion. Withthese features in mind, we attempt to realize the two steps in a single-temperature model with rapid temperatureuniformization.II. MODELSA. Ising-like modelAs aforementioned, the degeneracy of the LS and HS states plays key role for the SC phenomena. The HS state hasthe energy D and degeneracy nHS while the LS state has the energy −D and degeneracy nLS. Electrical situations ofthe LS and HS are depicted in Fig. 1, where the difference of the spins is shown (SLS = 1/2 and SHS = 5/2 for the LSand HS states, respectively). It should be noted that the degeneracies nHS and nLS come from not only the spin statebut also from the lattice degree of freedom. Therefore, the ratio g = nHS/nLS is much larger than that of the spin partgspin = (2SHS +1)/(2SLS +1). To treat the effect of degeneracy, the so-called Ising-like model has been introduced [2–FIG. 1: Schematic picture of the LS and HS states of Fe3+ ion. (Courtesy of Professor Hiroko Tokoro)4]. Let HS spin states of the ith ion be noted by si = +1(i = 1, · · ·nHS) and LS by si = −1(i = 1, · · ·nLS). Then, thepartition function is given byZ =n∑si=−1, · · · ,−1︸ ︷︷ ︸nLS,1 · · · 1︸ ︷︷ ︸nHSe−βEi = nHSe−βD + nLSeβD =√nHSnLS(√nLSnHSeβD +√nHSnHLe−βD)=√nHSnLS(eβD−kBT ln(nHS/nLS))/2 + e−β(D−kBT ln(nHS/nLS))/2)=∑σ=±1e−βHeffσ, (1)whereHeff = −D +12kBT ln g, g =nHSnLS. (2)3Thus, the original model with degenerate states is expressed by an Ising model (σ = ±1 without degeneracy) withthe effective temperature-dependent field Heff . This transformation is also good even if the Hamiltonian containsinteractions among spin states as explained in Computational Details Part A. Thus, models with interactions amongspin states:HIsing = −J∑<ij>sisj +D∑isi (3)is expressed byHIsing = −J∑<ij>σiσj −Heff∑iσi. (4)For the phase transitions of SC systems, the cooperative interaction is inevitable. The importance of the cooperativityamong spins for the transition was pointed out from the experiments of the measurement of the heat capacity bySorai and Seki, and they proposed the so-called domain model [23]. The relation between the domain model andthe Ising-like model was clarified by a microscopic theoretical study [24]. This type of modeling is good for bistablesystems with different degeneracies, i.e., not only for SC system but also charge transfer model and charge transferSC model [25].Thermal properties of the spin crossover model are obtained from the phase diagram (T,H) of corresponding Isingmodel [3, 4]. There appear several different types of ordering processes as depicted in Fig. 2, i.e., (I) smooth crossover,(II) smooth transition with low-temperature metastable HS phase, (III) first-order phase transition, (IV) first-orderphase transition with low-temperature metastable HS phase, and (V) hidden LS phase. Temperature dependences of1 2−2−1012 H/TCT/TC(I)(II)(III)(IV)(V)coexistence linedisordered statespinodal fieldFIG. 2: Phase diagram of an Ising model and various cases of spin crossover system which are expressed by lines of the slopegiven by 12kB ln(nHSnLS).HS fraction fHSfHS =〈σi〉+ 12(5)to the Types I ∼ V are depicted in Fig. 3. One of characteristics of SC model is the so-called LIESST (light inducedexcited spin state transition) [26, 27]. That is, when the LS state is excited to HS state by an photo-irradiation, theexcited HS state has a long lifetime. There may be two cases of LIESST. One of them is caused by long relaxation time40 5 10 1500.510fHSTemperature/JType I (D=12 g=20)kBT/J=2D/ln g0 5 1000.510fHSTemperature/JType II (D=5 g=20)kB T/J=2D/ln g0 2 4 600.510fHSTemperature/JType III (D=6 g=20)kB T/J=2D/ln g0 2 4 600.510fHSTemperature/JType IV (D=5.5 g=20)kB T/J=2D/ln g0 2 4 600.510fHSTemperature/JType V (D=5 g=20)kBT/J=2D/ln gFIG. 3: Corresponding temperature dependence of HS fraction to Fig. 2.at low temperatures where all the processes change very slowly. The other is the case where the HS state is trappedat the metastable branch in case IV. In this case fHS first relaxes to the local stable value, although it finally relaxesto LS state after a long time. [27] Systematic changes from the type III to IV have been studied in Rb2Mn[Fe(CN)6families, where the case (IV) was found in the experiments [28] and also the case (V) in the experiment [29]. Whenthe interaction between the nearest neighbor sites is antiferromagnetic, the corresponding phase diagram of the Isingmodel becomes complicated and the SC ordering becomes two-step [30]. When the interactions become frustrated,e.g., the case in which signs of interactions of nearest neighbor and the next nearest neighbor pairs are different,two-step and more complicated ordering process appear which have been found in experiments [17, 31]B. Models for elastic interactionElastic interaction has been recognized as an ingredient of the SC transition and lead to various interesting proper-ties [32]. The volume, more generally local lattice structure, changes between the LS and HS states, and the changecauses deformation of the lattice [6–9]. This deformation causes a long-range interaction among spins used the elasticenergy of the lattice [10].Here we briefly review models so for studied for the elastic interaction in spin-crossover systems.1. Model AThe elastic interaction was firstly studied by the model A explained below. In the model, motions of quantities aretreated by a molecular-dynamics simulation with Nose-Hoover thermostat [6, 33]. This model is suitable to reproducevarious dynamical properties.HA =N∑iP 2i2M+N∑ip2i2m+N∑iV intrai (ri) +∑<ij>V interij (Xi,Y j , ri, rj). (6)5Here, the position and the momentum of the i-th molecule are given by (Xi,P i). The coordinate (molecular radius)and the corresponding momentum representing the intra-molecular motion (breathing mode) are denoted by (ri, pi).The elastic interaction among molecules forming a lattice structure (Fig. 4) is given by the excess energy due todeformation of nearest neighbor pairsV interij (Xi,Y j , ri, rj) = f(dij), (7)where f(x) is a function of the distance between molecules dij = |Xi −Xj | − (ri + rj). The concrete form whichwas used in the paper [6] has the form f(x) = D(ea′(x−x0) + e−b′(x−x0))with a′ = 0.5 and b′ = 1.0. The value x0 ischosen to make f(x) have the minimum at x = 0. An alternate form of f(x) could be the form of Hlattice (13) givenbelow.To keep the lattice to be a not titled square shape, we include interaction between the next nearest neighbor (nnn)pairs in addition, which is small compared to that of the nearest neighbor. If we use the triangular lattice, we do notneed the nnn interaction.FIG. 4: Elastic interaction V interij (Xi,Y j , ri, rj) among molecules (Courtesy of Professor Hiroko Tokoro).The energy for the internal degree of freedom is given by V intrai (ri) as a function of the radius of the molecule ri.The effective potential energy for the radius ri, Vintrai (ri), is given by an bistable form depicted in Fig. 5. Representingthe differences of energies and of degeneracies, the potential energies for the LS and HS states are given by harmonicpotentials byVLS =12aLS(ri − rLS)2, VHS =12aHS(ri − rHS)2 + 2D, (8)respectively. The coefficients aLS and aHS are taken to be proportional to the stiffness constants of breathing modesof the LS and HS states, respectively. We assume that the potential energies are hybridized by a mixing term J , andthen the hybridized potential energy E is given by∣∣∣∣ VLS − E JJ VHS − E∣∣∣∣ = 0, (9)which givesE± =VLS + VHS ±√4J2 + (VHS − VLS)22. (10)These energies are depicted in Fig. 5. In actual spin-crossover materials, there exists a large entropy difference betweenthe LS and HS states. As mentioned, it should be noted that the degeneracies nLS and nHS come from not only thespin part but also contributions from surrounding ligands. To treat such large entropy difference in the moleculardynamics model (6), an oscillator coupling to the spin state in each molecule was introduced, in which the frequencyof the oscillator varies depending on the spin state [6]. To control temperature and/or pressure, Nose-Hoover thermalbath and/or Andersen pressure bath were used [6] . With this modeling, various peculiar phenomena due to theelastic interactions, e.g., macroscopic nucleation phenomena [11], have been studied. This model is suitable to studyeffects of mechanical dynamics at zero and finite temperatures because the motion is given by molecular dynamics.6rLS rHSrV−(r)V+(r)VHS(r)VLS(r)FIG. 5: Schematic picture of the potential energies VLS (black), VHS (green), V+ (blue, dash), and V− (blue, solid). V− is usedfor V intrai (ri).2. Model BIn the model B, the motions of both spin and lattice are given by transition probabilities with the detailed balance [7].In this model, the energy of the positions ({Xi}) and the spin states {σi} is given byHB = Hspin +Hlattice. (11)The spin part isHspin =∑i(D − 12kBT ln g)σi, g =nHSnLS, (12)and the lattice part isHlattice =12∑<ij>k1 (|Xi −Xj | −R(σi)−R(σi))2+12∑<il>k2(|Xi −X l| −√2(R(σi) +R(σj)))2, (13)where < ij > denotes all the nearest neighbor pairs and < il > the next nearest neighbor pairs in a square lattice asdepicted in Fig. 4. Here, the second term is introduced to keep the lattice to be square shape. Thus, k2 is set to besmall compared to k1. R(σ) denotes the radius of the molecule of the spin state of σ. In the model B, both the motionof position (lattice) and the change of spin state are treated by transition probabilities (i.e., as usual Monte Carlomethods). By this model, various interesting properties due to elastic interaction have been studied. For example, themean-field type long-range ferromagnetic interaction among spins [10], distribution of local stress [13], a complicatedphase diagram due to competition between a short-range antiferromagnetic interaction and the long-range interactiondue the elastic model [14–16], and multistep transitions [17].Furthermore, this model can be applied to study dynamical properties. Using different timescales for the spin-state conversion by MC and lattice motion by MD, it was shown that the nature of interface growth (domain wallpropagation) varies depending on the relative timescales [12].3. Model CThermo-mechano-elastic method proposed by Enachescu consists of over-damped lattice dynamics and thermalconversion of spins [8]. In this model, the positions of the molecules evolve according to the equation:md2Xidt2= F i − µdXidt, (14)where µ is a damping constant. This equation brings the system to the low energy state. Because of the friction,in a long time, the state reaches the ground state (or at least local minimum state) and the motion does not take7the temperature effect onto account. But, when the mass of molecules is large, it gives a good description of latticemotion. Also, for practical calculation, the lattice dynamics is performed in a finite time and the lattice state doesnot go to the minimum energy position. To improve this point, one may use the Nose-Hoover formalism [33, 34]or introducing a noise to realize the equilibrium state by the fluctuation and dissipation relation [35]. For the spinconversion, Monte Carlo update with the following transition probability is used:P iHS→LS = 1τ exp(−EA − κpikBT),P iLS→HS = 1τ exp(−EA + κpikBT)exp(−D − kBT ln gkBT), (15)where pi is the sum of the pressures from the neighboring sites to the ith molecule. In the triangular lattice, weconsider only nearest neighbor pairs, and the elastic interaction and the force are given byHlatice =12kδx2ij , F = −kδxij , (16)where δxij is a change of the distance from the equilibrium length. Then the pressure is given bypi =∑j:neighboring siteskδxij . (17)The relation between models B and C has been studied in Ref. [9]. This model has been applied for several SCphenomena, e.g., photo-thermal swotching [18, 20, 21]. Using over-damped lattice motion in Model B is physicallyequivalent to Model C. In our previous papers, we utilized Model C; however, in the present paper, we employ ModelB with over-damped lattice motion, which produces qualitatively (and almost quantitatively) the same results asthose obtained from Model C.C. Models for dynamics after ultrafast pulse photo-irradiationTo study mechanism of the elastic and thermal steps after ultrafast pulse photo-irradiation, models with thermalconduction is necessary. Here, we review models so far studied.It should be noted that in the simulation, the Hamiltonian (11) is only for the potential energy which consists ofthe spin part (12) and elastic part (13). In classical systems, the kinetic part is independent of the potential part.The temperature somehow represents kinetic energy, and the transfer of kinetic energy is treated by the heat transferwhich is not included in the Hamiltonian. To study dynamics of systems with non-uniform temperatures, we needsome heat transfer equation in addition to the Hamiltonian, which will be explained in the following subsections.1. Single-temperature modelIn this model, we consider a local temperature, TL(x, t), for a molecule at a position x at a time t. Providing heattransfer between neighboring sites and also heat escape to external thermal bath with the temperature Tbath, the timeevolution of the local temperature, {TL(x, t)} is given byTL(x, t+ ∆t) = TL(x, t) + ∆t[λ∑nn:x′(T (x′, t)− T (x, t)) + κ∑x(Tbath − T (x, t))], (18)where λ gives the rate of energy diffusion to the neighboring sites, κ is the rate of heat escape to the external thermalbath. The unit of λ and κ is 1/MCS. In the present study, we demonstrate qualitative features of the processes, andthus here we choose the parameters to reproduce phenomena with typical features. In this respect, we study models intwo dimensions. Real systems are in three dimensions where we must consider the inhomogeneous irradiation from thesurface to the bulk, which is not considered in the present paper although it was studied the previous papers [20, 21].In the present paper, we use the following parameters: The lattice is a 60 × 60 square lattice, and D = 1064K,g = 1069, k1 = 80000K/nm2, k2 = 8000K/nm2 (the equilibrium distance rLS = R(−1) = 1 and the length is measuredby the unit nm), and the bath temperature Tbath = 280K. In Fig. 6, we show the temperature dependence of fHS inthe present model with the parameters. In the present study, 2D/kB ln g ' 305, the bath temperature is Tbath = 280K8which is set just below the hysteresis loop. The average temperature Tall just after the irradiation with TEX = 480Kis given byTall = Tbath × 0.77 + TEX × 0.23 ' 326K. (19)250 300 350−101TfHSTbathTall(480)T0=2D/kBln gFIG. 6: Temperature dependence of fHS of the present model.Here, we consider the following processes after photo-irradiation. Before the photo-irradiation, the system is inequilibrium at a temperature Tbath which is just below the coexistence (hysteresis) region and the spins are in LSstate. At the irradiation, some of molecules are excited and the local temperature increases to TEX as shown in Fig. 7.FIG. 7: (left) Before the photo-irradiation where almost all the local sites are in LS state (dark blue), and (right) a molecule(yellow) is excited by the photo irradiation.We consider a case in which some of molecules are excited. Because the excited molecules tend to expand to the sizeof HS state, very high pressure is placed on the neighboring sites. Microscopically, the Fe-ligand distance increasesin an excited molecule, but the volume (lattice parameter) of the whole system is little changed immediately afterirradiation. Thus, the high pressure is applied to the neighboring molecules. In the present paper, we do not refer tothe details of this fast process although it is important as a microscopic process.The energy diffuses to the neighboring sites, and the temperature of the irradiated site decreases heating up theneighboring sites (Fig. 8(a)). In this process, some of the irradiated sites are pushed back to LS state (Fig. 8(b)) dueto the high pressure, and the number of HS molecules decreases just after the irradiation. When the excited spin keepsHS state, high pressure remains, and the high pressure propagates over the system and the size of system expands.This expansion of the system causes reduction of the pressure, and makes the molecules to convert to HS, whichmakes an increase of the HS fraction. When the energy transfer is slow, after the increase the number of HS moleculesdecreases because the decrease of the excited HS molecules, which makes a peak of HS fraction. These processes givethe elastic step. This step occurs at the timescale of pressure propagation (an order of the sound velocity which is anorder of few nano second). After a long time (in timescale of heat diffusion of an order of µs), temperatures of all thesites become higher than the conversion temperature, and all the sites tend to become HS state, which is the thermal9(a) (b)FIG. 8: (a) Conversion neighboring sites to HS state (orange), and (b) a case where the excited site is pushed back to LS statewith high energy (right blue).step. After much longer time, energy escapes to the external thermal bath, and the temperature comes back to theoriginal one of the LS state where almost all the spins convert to LS again. Like Fig. 2 of Ref. [21] obtained by themodel C, we show an example of such process in Fig. 9. Time dependence of nHS (HS fraction) is depicted in Fig. 9.Lower panel shows nHS in logarithmic time. The change at the irradiation occurs at t = 100 which is very close tothe left axis in the upper panel. The resulting state depends on parameters of excitation, heat transfers, and the heatleaking to the bath. By choosing an appropriate set of parameters, we can reproduce the two steps similarly to thatobserved experiment [19].In the figure, we also plot time dependence of the quantities “vol(t)/vol(0)”, a(t)/a(0), 10×∆a/a(0) and 5×E(t)−E(0)/E(0). To see the change clearly, multiplied values are plotted. “vol” is the volume (area) of the lattice andvol(0) is the value in the LS state, and a is the average distance of the neighboring sites i and j,√(〈|Xi −Xj |2〉,∆a/a(0) = (a(t)− a(0))/a(0)) and E is the sum of local temperature of the current configuration. In ComputationalDetails Part B, a vs nHS is plotted in Fig. 15.nHSa(t)/a(0) x 10vol(t)/vol(0)nHSt(E(t)−E(0))/E(0) x5log(t)600001001201000200034000∆a(t)/a(0)1.000.5FIG. 9: Example of two steps in a single-temperature model with rHS/rLS = 1.102, λ = 0.0002, and κ = 0.000017 in a two-dimensional system with 60 × 60 sites. “vol” is the volume of the lattice, E the total energy, “a” the average lattice constantand nHS is HS fraction. The lower panel shows nHS in the logarithmic time. Temperature of the irradiated site is TEX = 480K,and the probability of irradiated sites is p = 0.23.However, it is found that the spin configuration is very inhomogeneous as depicted in Fig, 10 similarly to Fig. 2(right) of Ref. [21]. This inhomogeneity of local temperature contradicts experimental observation that the thermal-ization of the lattice is very fast around 10ns which is about the time of elastic step [22]. They found that the dynamicsspan from subpicosecond molecular photo switching followed by volume expansion (on a nanosecond timescale) andadditional thermal switching (on a microsecond timescale) from results of time-resolved X-ray diffraction. Therefore,uniformization of the temperature should be achieved much before the thermal step.If λ is large, local temperature becomes uniform fast, but the dominant spin conversions are given by thermalconversion (i.e., thermal step), and the elastic step does not appear.10time=120 time=500 time=1000time=2000 time=10000 time=34000FIG. 10: Spin configurations at times. The colors denote different temperatures: yellow : T ≥ 390, orange: 390 > T ≥ 370,dark brown: 370 > T ≥ 350, light brawn: 350 > T ≥ 330 green : 330 > T ≥ 310, light blue: 310 > T ≥ 290, dark blue :290 > T ≥ 280, black : 280 > T .2. Two-temperature modelTo overcome this difficulty, we introduced a two-temperature model [21]. In this model, we consider two kinds oftemperatures, i.e., the lattice temperature TL(x, t) which expresses the total energy at a molecule at a position x ata time t and the spin temperature TS(x, t) which is responsible for the spin conversion. The temperatures follw thefollowing equations:TL(x, t+ ∆t) = TL(x, t) + ∆t[λ∑nn:x′(TL(x′, t)− TL(x, t)) + κ∑x(Tbath − TL(x, t)) ,](20)andTS(x, t+ ∆t) = TS(x, t) + ∆tλ′∑x(TL(x, t)− TS(x, t)) , (21)where λ gives the rate of energy diffusion to the neighboring sites, κ is the rate of heat escape to the external thermalbath, and λ′ is the rate of conversion from the lattice temperature to the spin temperature. In our recent paper [21],we used slightly different model. But here we use this model which is based on essentially the same picture.We assume the spin state at each site is given by the equilibrium state with its spin temperature. Here we assumethat λ � λ′, so that the energy diffuses rapidly, and the spin temperature comes close to the lattice temperaturerather slowly. The escape rate to the external bath κ is assumed to be much smaller than λ′. In this case, thelattice temperatures become uniform in a short time, while the spin temperature remains inhomogeneous. Andthus, tuning these parameters, we can produce the elastic step and thermal step, separately in almost uniform latticetemperatures [21]. A two-step of HS fraction is depicted in Fig. 11 which corresponds to Fig. 3 in our recent paper [21].In the following section, we consider properties of the timescales of spin conversion in various cases.11nHSa(t)/a(0) x 10vol(t)/vol(0)nHSt(E(t)−E(0))/E(0) x5log(t)60000100120 300 75012000∆a(t)/a(0)1.000.5time=300(a) (b)FIG. 11: (a)Example of two steps in a two-temperature model with rHS/rLS = 1.1, λ = 0.02, λS = 0.001and κ = 0.00005 ina two-dimensional system with 60 × 60 sites. The notations are the same in Fig. 9. (b) Profile of the lattice temperature atthe time t = 300 denoted by the arrow in (a). The temperature of the irradiated site is TEX = 480K, and the probability ofirradiated sites is p = 0.23.III. TIMESCALES OF PROCESSESNow we consider the timescales of processes. For thermodynamics, the timescales of relaxation process are notimportant. But for dynamical properties, the timescales have crucial role. Here we consider the following points(1) timescales of the degeneracy model and the Ising-like model,(2) timescales of thermal process and mechanical process,(3) timescales of heat emission to the local bath in mechanical spin-conversion process.A. Timescales of the degeneracy model and the Ising-like modelLet us begin with (1) timescales of the degeneracy model and the Ising-like model. For spin-crossover phenomena,the degeneracy, i.e., the difference of degeneracy of the LS and HS states is important. For equilibrium thermodynamicproperties, the effect of the different degeneracies can be expressed by an effective temperature dependent field (2).At high temperatures, the system tends to be HS state because of the entropy effect, i.e., nHS > nLS as depicted inFig. 12(left). In the Ising-like model, the effective field (2) becomes positive and the energy of HS state in the Ising-like model becomes lower as depicted in Fig. 12(right). These two pictures have the same thermodynamic propertiesbecause their partition functions are the same.Now we study dynamical properties of the models. In particular, we obtain the relaxation rates of both models inMonte Carlo simulations which give qualitative nature of relaxations.Dynamics of Monte Carlo simulations are given by the so-called master equation which is equation of motion of theprobability distribution updated by the transition probabilities among the states. Giving the transition probabilityper a unit time from the state i to j by {wi→j}, the master equation is given byP (i, t+ ∆t) = P (i, t)−∆tn∑j 6=iwi→jP (i, t) + ∆tn∑j 6=iwj→iP (j, t), i = 1, · · · ,M. (22)Expressing the probabilities as a vector P (t), the master equation is given by in a matrix form with the time evolution12FIG. 12: (left)Comparison of spin conversion in the original degenerate model, and (right) that in the Ising-like model.operator L:P (1, t+ ∆t)P (2, t+ ∆t)...P (M, t+ ∆t) =w1→1 w2→1∆t · · · wM→1∆tw1→2∆t w2→2 · · · wM→2∆t....... . ....w1→M∆t w2→M∆t · · · wM→MP (1, t)P (2, t)...P (M, t)→ P (t+ ∆t) = LP (t), (23)Here, we assume the detailed balance among the transition probabilities to realize the thermal equilibrium state:wi→jwj→i=Peq(j)Peq(i). (24)In the present study, we adopt the choice of Glauber type transition probability for each spin:wi→j ∝ e−βEj , β =1kBT, (25)andwi→i = 1−M∑i 6=jwi→j∆t. (26)In the model with explicit degeneracies, the energy difference of the LS and HS states is 2D, and energetically LS isfavorable. But because nHS is much larger than nLS, HS state is realized at high temperatures. In the Ising-like model,the variable σi just takes ±1, but the effective energy includes the degeneracy effect as the temperature dependentfield (2). And thus, HS state is realized at high temperatures. These processes are different, and we need to check thedifference in dynamics. . The time-evolution operator for the degenerate model is expressed by a master equation isgiven by an n× n matrix, where n = nHS + nLS, while that of the Ising-like model is a 2× 2 matrix. The relaxationrate is obtained from the eigenvalues of the matrices.1. Timescales of the degeneracy modelThe time evolution operator for a step of ∆t in the original degenerate model is given byP (t+ ∆t) = LP (t) =α β · · · β γ′ · · · γ′β α................ . ..... . . β......β β β α γ′ · · · γ′γ · · · γ α′ · · · β′.........γ · · · · · · γ β′ · · · α′p1......pnHSq1......qnLS, (27)13whereβ = wHS→HS =1τHH∆t, γ = wHS→LS =1τHLeD/kBT∆t, (28)β′ = wLS→LS =1τLL∆t, γ′ = wLS→HS1τLHe−D/kBT∆t, (29)α = 1−∑j 6=iw1→j = 1− (nHS − 1)β − nLSγ, (i ≤ nHS), (30)andα′ = 1−∑j 6=iw1→j = 1− (nLS − 1)β′ − nHSγ′ (nHS < i ≤ n). (31)Here we assumeτLH = τHL = τX. (32)Now we obtain eigenstates and eigenvalues. In general, one of the eigenstates of the time evolution operator Lrepresents the steady (equilibrium) state with the eigenvalue λ0 = 1, and other eigenstates correspond to relaxationmodes with the eigenvalues λi < 1 (i = 1, 2, · · ·M − 1) where M = nHS + nLS is the dimension of the time evolutionoperator L.From the symmetry of the matrix L, first, we consider a uniform state in which all the HS states are uniformlypopulated with a probability p and all the LS states are uniformly populated with a probability q:p1......pnHSq1......qLS=p......pq......q. (33)In this case, the eigenvalue equation with the eigenvalue λ is given byp(α+ (nHS − 1)β) + qnLSγ′ = λp, q(α′ + (nLS − 1)β′) + pnHSγ = λq. (34)Substituting (30) and (31), we have the relations:p(1− nLSγ) + qnLSγ′ = λp, q(1− nHSγ′) + pnHSγ = λq. (35)Thus, the λ is given by∣∣∣∣ 1− nLSγ − λ nLSγ′nHSγ 1− nHSγ′ − λ∣∣∣∣ = 0→ (λ− 1)(λ− (1− nLSγ − nHSγ′)) = 0. (36)The eigenvalues areλ0 = 1, λ1 = 1− nLSγ − nHSγ′ = 1− nLS1τXeD/kBT∆t− nHS1τXe−D/kBT∆t = 1− ∆tτX(nLSeD/kBT + nHSe−D/kBT).(37)The eigenvalue λ0 is for equilibrium state and the eigenvalue λ1 gives the relaxation time.(λ1)n = e−t/τ , t = n∆t,→ 1τ= − lnλ1∆t(38)14= − 1∆tln(1− ∆tτX(nLSeD/kBT + nHSe−D/kBT))' 1τX(nLSeD/kBT + nHSe−D/kBT). (39)The ration of p and q for the case of λ0 = 1 is given by−pnLSγ + qnLSγ′ = 0→ p =γ′γq → p : q = nHSe−D/kBT : nLSeD/kBT . (40)This gives the equilibrium state which corresponds to the canonical distribution.The ration of p and q for the case of λ = 1− nLSγ − nHSγ′ is given bypnHS + qnLS = 0→ p = −nLSnHSq. (41)The sum of all the components of P 1 is zero (pnHS + qnLS = 0), which indicates that its eigenstate gives a relaxationmode giving a deviation from the equilibrium state.Next, we obtain non-uniform eigenstates which are orthogonal to the uniform eigenstates obtained above. Then,we consider states satisfying the following relations:nHS∑i=1pi = 0,n=nHS+nLS∑j=nHS+1qi = 0, (42)i.e., eigenvectors of the forms: p1p2...pnHS0...0, and00...0q1...qnLS. (43)The operation for the upper part (of HS states) is given byα β · · · ββ. . . · · · β....... . ....β · · · · · · αp1p2...pnHS = (α− β)p1p2...pnHS , (44)where we used the condition∑nHSi=1 pi = 0. Thus, the eigenvalues are (nHS − 1)-fold degenerate to beλ = α− β = 1− nHSβ − nLSγ. (45)Similarly, from the lower part, we find (nLS − 1)-fold degenerate eigenvalues:λ = α′ − β′ = 1− nHSβ′ − nLSγ′. (46)Because uniformizations within HS states or within LS states are expected to be fast, i.e., we assume these eigenvaluesfor non-uniform cases are small. Therefore, the relaxation process the uniform modes with slow relaxation time playsimportant role.2. Timescales of the Ising-like modelIn this case that the system has only two states (σ = ±1), the time evolution matrix is given byP (t+ ∆t) = LP (t) =(1− wHS→LS∆t wLS→HS∆twHS→LS∆t 1− wLS→HS∆t)(pHSpLS), (47)15wherewHS→LS =1τ0eβ(D−kBT ln g/2) =1τ0√1geβD =1τ0√nHSnLSnLSeβD,wLS→HS =1τ0e−β(D−kBT ln g/2) =1τ0√ge−βD =1τ0√nHSnLSnHSe−βD. (48)The eigenvalues are obtained by∣∣∣∣ 1− wHS→LS∆t− λ wLS→HS∆twHS→LS∆t 1− wLS→HS∆t− λ∣∣∣∣ = 0→ (λ− 1)(λ− (1− wHS→LS − wLS→HS)) = 0. (49)asλIsing−like0 = 1, λIsing−like1 = 1− wHS→LS∆t− wLS→HS∆t. (50)The relaxation time is given by(λIsing−like1 )n = e−t/τIsing−like, t = n∆t,→ 1τ Ising−like= − lnλIsing−like1∆t' wHS→LS + wLS→HS. (51)Comparing τ for the degenerate model and τ Ising−like, i.e.,1τ=1τX(nLSeβD + nHSe−βD) , 1τ Ising−like=1τ0√nHSnLS(nLSeβD + nHSe−βD) , (52)we find that the relaxation times of the Ising-like model and the degeneracy model have the same dependencies onthe temperature. Here τX is a typical time of contacts for each microscopic state si with the thermal bath, whileτ Ising−like is a typical time of contacts for one of HS state and LS state σi with the thermal bath. For the overalltimescale, they are different by the factor τ0√nHSnLS/τX. But this factor is constant throughout the process and thedifference does not affect temperature dependence of the relaxation process. Therefore, we conclude that both modelshave the same dynamics, and the use of the Ising-like model is verified to study the original degenerate model.B. Considerations on dynamical effects on the spin conversionsIn Sec. II C 2, we considered the two-temperature model to explain the two-step evolution of spin conversions in therapid uniformization of lattice temperature. We may have a question what the essential role of the two temperaturesis, and what physical meaning of the spin temperature is. The answer should be in the timescales of microscopicprocesses, e.g., energy transfer and spin conversion, and so on.The two-temperature model suggests the importance of the timescales of processes. Here, we explore some mech-anisms which make the elastic step in a rapid uniformization of the temperature in a single-temperature model bymodifying the transition probability considering dynamical effects on the spin conversion.As we discuss above, the dynamics of the spin conversion is not simply given by thermal process which is given byMC simulation with the detailed balance, and we need to consider properies of direct mechanical process. In principle,the latter cannot be taken into account in MC. But below, we try to argue modelings with modifying the transitionprobability.Here, we consider the timescales of thermal process and mechanical process. Because of the degeneracy, spin stateschange not only due to mechanical force (the pressure), but also due to the entropic force due to the difference ofdegeneracies g. When the mechanical deformation is small, spin conversion takes place as thermal process (Fig. 13(a)),which may be called entropy-driven conversion. Figure 13(a) schematically shows that the effective energy of HS partis lowered by the entropy effect D − kBT ln g/2 although the energy of HS part itself is higher by D. On the otherhand, when the mechanical deformation is large, spin conversion takes place by a mechanical force, e.g., the pressure(Fig. 13(b)), which may be called mechanical-driven conversion. Figure 13(b) schematically shows that HS state isforced to convert to LS state due to large pressure. In this process, heat emission takes place in a short time and thelocal heat bath may not absorb it immediately.16(a) (b)FIG. 13: (a) Spin conversion caused by the degeneracy difference (entropy-driven conversion), (b) Spin conversion caused bymechanical force (mechanical-driven conversion).1. Timescales of thermal process and mechanical processWe expect the relaxation time for the entropy driven process to be longer than that of mechanical driven processesbecause the entropic force becomes relevant after many contacts of the heat bath to recognize the degeneracy difference.Although such an effect is difficult to be expressed by a transition probability satisfying the detailed balance, we dareto introduce dependence of the relaxation time as a function of mechanical energy to take into account such aspectapproximately and try to find two-step conversion within a single temperature model with fast uniformization oflattice temperature. For example, we may introduce dependence of the relaxation time as a function of mechanicalenergy.1τ0=a tanh(EM0E0)+ ba+ b, (53)where E0 is the total energy and EM0 is the mechanical part. a and b are control parameters. By tuning parameters,we may have two steps (not shown), but this choice of transition probability violates the detailed balance. Thus, inprinciple, the difference of dynamics must be taken into account by a mechanical way but not by the form of transferprobability, which should be an important subject in the future. Some trials of such treatments will be discussedlater.2. timescales of heat absorption to the local bath in mechanical spin-conversion processNext, we consider the timescale of heat flow to the thermal bath in spin conversion. In the master equationapproach, the transitions occur in contact with the thermal bath. There, the change of energy in the transition isassumed to be absorbed to the thermal bath immediately. This assumption is good when the bath is very large whereexchanges of the energy are done immediately and the temperature of the bath does not change.However, when we consider the local temperature, the local bath which consists of various modes of a moleculeis not so large. Then, it is expected to take some time before the emitted energy is absorbed in the bath. Dueto the emitted energy, the local temperature would change in a short time. As a typical case, let us consider theprocess where the photo-irradiated excited spin is pushed back to LS state due to the high pressure (Fig. 13(b)).In this process, the energy (i.e., heat) is emitted to the local environment (Fig. 14(a)). Therefore, the local latticetemperature increases. We try to reproduce the two-step spin conversion to take the heating due to the emissionwithin the single temperature model with fast uniformization of the temperature.It is found that, if we take into account temperature increase due to the emission, we can reproduce two stepsin a single-temperature model with tuned parameters (at t = 700 and 13000 in Fig. 14(b)). When λ is large, thetemperature profile becomes homogeneous in an early stage even in the single-temperature model. In Fig. 14(c), thetemperature profile at t = 200 is depicted.In this section, we explored possible mechanisms to realize two steps in a single-temperature model taking intoaccount some dynamical mechanisms, and we reproduce two-steps in a single temperature by considering additionaleffects at photo-excited sites. These results could lead to the hypothesis that the spin-temperature in the two-temperature model may represent the difference types of spin conversion (i.e., mechanical or thermal), or heat emissionat the spin-conversion of irradiated molecule. However, more systematic investigations on dynamical effects arenecessary.17nHSavol.nHStElog(t)5000010020070013000time=200(a) (b) (c)FIG. 14: (a) Heat emission in the process of compression by the pressure. (b) Two steps in a single-temperature model withtuned parameters taking into account the temperature increases due to compression. Example of two steps in a two-temperaturemodel with rH/rL = 1.1, λ = 0.02 and κ = 0.00005 in a two-dimensional system with 60 × 60 sites. the notations are thesame in Fig. 9. Temperature increase by heat emission is 100K. (c) temperature profile at t = 200. The temperature of theirradiated site is TEX = 400K, and the probability of irradiated sites is p = 0.23.IV. SUMMARY AND DISCUSSIONIn this paper, theoretical models for spin crossover phenomena were reviewed. Characteristics of spin crossoversystems are the difference of the degeneracy of the LS and HS states, i.e., nHS and nLS, and the elastic interactiondue to the difference of local lattice structure, e.g., the volume of molecules. The former is treated by the Ising-likemodel [2–4], and the latter is treated by several models taking the elastic energy due to lattice deformation [6–9].Moreover, to explain the dynamical aspects of ordering processes, such as the two-stage ordering after photo-irradiation(elastic step and thermal step) [19, 20], mechanism of heat (energy) conduction is important. For this point, single-temperature model and the two-temperature model have been introduced. The conventional single-temperaturemodel studied in Sec III C-1 cannot reproduce the two-step transition in a rapid uniformization of temperaturefound experimentally [22], and the concept of the spin-temperature which is a part of the energy responsible to thespin conversion has been introduced. This two-temperature model succeeded to reproduce the two steps in rapiduniformization of the lattice temperature [21]. The idea of the spin temperature suggests the importance of thetimescale of the processes. In the present paper, we studied timescales of the original degenerate model and the Ising-like model. It was found that the relaxation time of both models have the same temperature dependence althoughthe overall factor is different. Thus, it is verified to use the Ising-like model to study dynamics. Moreover, we try torealize the two steps, i.e., the elastic step besides the thermal step in rapid uniformization of the temperature in asingle temperature model. We considered some modifications of transition probability, taking into account effects dueto the difference of thermal and entropic forces and of local heat emission in the spin-conversion by pressure. It wasfound that the elastic step appears in such modifications.For real time phenomena, more investigations of time dependent processes are necessary, which will be importantfor fast manipulation of spin states by ultrafast pulses. As for the discussion on the timescales related to the thermalexchange between the molecule and its local heat bath, concrete information of intramolecular and intermolecularvibrations and ligand distortions is inevitable. Such microscopic details depend on materials. In the present work, weonly discussed the importance of timescales, and the details will be postponed to future studies. Moreover, the presentpaper does not refer to the microscopic processes of irradiation where the interaction between the electromagneticfield and the spins is important. Such processes are relevant also for general manipulation of quantum systems. Mostof works focus on control of spin state [36], but the dynamics with abrupt strong pumping such as the present case ofirradiation should be an important issue for the future studies.AcknowledgmentsWe would like to thank all the collaborators. The present work was partially supported by JSPS KAKENHI GrantNumber JP18K03444. During the preparation of the present paper, one of the authors, our great collaborator and18friend, Cristian Enachescu passed away. We would like to express our gratitude for his essential contribution to thispaper and his great achievements in the field of spin transition materials.Computational DetailsPart A: Ising-like model with interactionThe partition function for the system with interactions is given byZ =∑all the states of sieβ(J∑ij sisj−D∑Ni=1 si). (54)We express the Boltzmann factor for each state of {si} by {σi}:eβ(J∑ij sisj−D∑Ni=1 si) = eβ(J∑ij σiσj−D∑Ni=1 σi) (55)Taking the degeneracy of si, the partition function is given byZ =∑all the states of σiN∏ini(σi)eβ(J∑ij σiσj−D∑Ni=1 σi), (56)where ni(σi) is nHS when σi = 1 and nLS when σi = −1, which is expressed asni(σi) =√nHSnLSe12σi lnnHSnLS ={nHS for σi = 1nLS for σi = −1(57)Then, the partition function is given byZ =∑{σi=±1}eβ(J∑ij σiσj−D∑Ni=1 σi)√nHSnLS)Ne12 lnnHSnLS∑Ni=1 σi= (√nHSnLS)N∑{σi=±1}eβ(J∑ij σiσj+(−D+ 12kBT lnnHSnLS)∑Ni=1 σi). (58)Thus, the model is regarded as an Ising model with the effective temperature dependent field:Heff = −D +12kBT lnnHSnLS. (59)For a ferromagnetic model which has a critical temperature TC, the system exhibits a first-order phase transitionas function of the magnetic field H below the critical point. Thus, the present model exhibits a first order phasetransition at Heff = 0 ifkBT0 =2Dln nHSnLS< kBTC. (60)Thus, the phase transition occurs whenD <12kBTC lnnHSnLS. (61)Here we do not use any approximation, e.g., the mean-field approximation, and the condition of the first-order phasetransition holds generally.190 111.11.2nHSavol/vol(0)FIG. 15: Plot of a (blue solid line) and vol (black dotted line) vs nHS.Part B: Relations ”a” and “vol” vs. ”nHS”In Fig. 15, we plot a (blue solid line) and vol (black dotted line) vs. nHS. The arrows shows the movements in time.At an early time, nHS shows the elastic step and changes nonmonotonically, while ”a” and “vol” change monotonically.[1] O. Kahn, Molecular Magnetism (VCII, 1993), M. A. Halcrow(ed.), Spin-Crossover Materials(Wiley, 2013). P. Gütlich andA. Goodwin, Spin Crossover in Transition Metal Compounds, Vols. I-III. Topics in Current Chemistry, Vols. 233, 234,235. Springer-Verlag, Berlin, Heidelberg, New York (2004). Kamel Boukheddaden; Seiji Miyashita; Smail Triki, J. Appl.Phys. 132, 220402 (2022).[2] J. Wajnflasz and J. Pick: Phys (France) 32, C1-91-C91-92 (1971). A. Bousseksou, J. Nasser, J. Linares, K. Boukheddaden,and F. Varret: J. Phys. I France 2, 1381 (1992).[3] S. Miyashita, Y. Konishi, H. Tokoro, M. Nishino,K. Boukheddaden, F. Varret: Prog. Theor. Phys. 114, 719 (2005).[4] S. Miyashita, Proceedings of the Japan Academy Series B 86, (643-666) (2010).[5] S. Miyashita, Collapse of Metastability, Springer-Nature (2022).[6] M. Nishino, K. Boukheddaden, Y. Konishi, and S. Miyashita, Phys. Rev. Lett. 98, 247203 (2007). K. Boukheddaden, M.Nishino and S. Miyashita, Phys. Rev. Lett. 100, 177206 (2008). M. Nishino, K. Boukheddaden and S. Miyashita, Phys.Rev. B 79, 012409 (1-4) (2009).[7] Y Konishi, H. Tokoro, M. Nishino and S. Miyashita, Phys. Rev. Lett. 100, 067206 (1-4) (2008).[8] C. Enachescu, L. Stoleriu, A. Stancu, and A. Hauser, Phys. Rev. Lett. 102, 257204 (2009). C. Enachescu, L. Stoleriu, A.Stancu, and A. Hauser, Phys. Rev. B 82, 104114 (2010).[9] C. Enachescu, M. Nishino, S. Miyashita, K. Boukheddaden, F. Varret, and P. A. Rikvold, Phys. Rev. B 91, 104102 (2015).[10] S. Miyashita, Y. Konishi, M. Nishino, H. Tokoro and P.A. Rikvold, Phys. Rev. B 77, 0144105 (2008).[11] M. Nishino, C. Enachescu, S. Miyashita, P. A. Rikvold, K. Boukheddaden, and F. Varret, Sci. Rep. 1, 162 (2011).[12] M. Nishino, T. Nakada, C. Enachescu, K. Boukheddaden, and S. Miyashita, Phys. Rev. B 88, 094303 (2013).[13] A. Slimani, K. Boukheddaden, F. Varret, H. Oubouchou, M. Nishino and S. Miyashita, Phys. Rev. B 87, 014111 (2013),A. Slimani, K. Boukheddaden, F. Varret, M. Nishino and S. Miyashita, J. Chem. Phys. 139, 194706 (2013).[14] T. Nakada, T. Mori, S. Miyashita, M. Nishino, S. Todo, W. Nicolazzai and P. A. Rikvold, Phys. Rev. B 85, 054408(1-8)(2012).20[15] P. A. Rikvold, G. Brown, S. Miyashita, C. Omand, and M. Nishino, Phys. Rev. B 93, 064109 (2016). C. H. Chan, G.Brown, and P. A. Rikvold, Phys. Rev. E 95, 053302 (2017), Phys. Rev. B 96, 174428 (2017). M. Nishino, S. Miyashita,and P. A. Rikvold, Phys. Rev. B 96, 144425 (2017). M. Nishino, P. A. Rikvold, C. Omand, and S. Miyashita, Phys. Rev.B 98, 144402 (2018).[16] Yogendra Singh, Hassane Oubouchou, Masamichi Nishino, Seiji Miyashita, and Kamel Boukheddaden, Phys. Rev. B 101,054105 (1-15) (2020).[17] Masamichi Nishino, Cristian Enachescu, and Seiji Miyashita, Phys. Rev. B 100, 134414 (2019). M. Nishino, Y. Singh, K.Boukheddaden, and S. Miyashita, Tutorial on elastic interaction models for multistep spin-crossover transitions, J. Appl.Phys. 130, 141102 (2021).[18] Y. Hu, M. Picher, N. M. Tran, M. Palluel, L. Stoleriu, N. Daro, S. Mornet, C. Enachescu, E. Freysz, F. Banhart, G.Chastanet, Adv. Materials 33, 2105586 (2021).[19] R. Bertoni, M. Lorenc, H. Cailleau, A. Tissot, J. Laisney, M. Boillot, L. Stoleriu, A. Stancu, C. Enachescu, and E. Collet,Nat. Mater. 15, 606 (2016).[20] C. Enachescu, L. Stoleriu, M. Nishino, S. Miyashita,3 A. Stancu, M. Lorenc, R. Bertoni, H. Cailleau, and E. Collet, Phys.Rev. B 95, 224107 (2017).[21] Laurentiu Stoleriu, Masamichi Nishino, Seiji Miyashita, Alexandru Stancu, Roman Bertoni, Eric Collet, Maciej Lorenc,and Cristian Enachescu, Phys. Rev. B 108, 014306(1-11) (2023).[22] H. Cailleau, M. Lorenc, L. Guerin, M. Servol, E. Collet, and M. Buron-Le Cointe, Acta Crystallogr A Found Crystallogr.66, 189 (2010).[23] M. Sorai and S. Seki, J. Phys. Chem. Solids, 35, 555 (1974), M. Sorai, Bull. Chem. Soc. Jpn. 74, 2223 (2001).[24] M. Nishino, S. Miyashita and K. Boukheddaden, J. Chem. Phys. 118, 4594-4597 (2003).[25] S. Miyashita and N. Kojima, Prog. Theor. Phys. 109, 729-739 (2003).[26] S. Zerdane, L.Wilbraham, M. Cammarata, O. Iasco, E. Riviere, M. L. Boillot, I. Ciofini, and E. Collet, Chem. Sci. 8, 4978(2017). A. Marino, P. Chakraborty, M. Servol, M. Lorenc, E. Collet, and A. Hauser, Angew. Chem. Int. Ed. 53, 3863(2014). S. Decurtins, P. Gütlich, C. P. Kohler, H. Spiering, and A. Hauser, Chem. Phys. Lett. 105, 1 (1984). A. Hauser,Top. Curr. Chem. 234, 155 (2004).[27] Y. Konishi, H. Tokoro, M. Nishino and S. Miyashita, J. Phys. Soc. Jan. 75, 114603 1-9 (2006).[28] H. Tokoro, S. Miyashita, K. Hashimoto, and S. Ohkoshi, Phys. Rev. B73, 172415 (2006).[29] H. Tokoro and S. Ohkoshi, Appl. Phys. Lett., 93, 021906 (2008).[30] M. Nishino and S. Miyashita, Phys. Rev. B 88, 014108 (1-14) (2013).[31] N. Bréfuel, H. Watanabe, L. Toupet, J. Come, M. Kojima, N. Matsumoto, E. Collet, K. Tanaka, and J. P. Tuchagues,Ang. Chem. Int. Ed 48, 9304 (2009). H. Watanabe, K. Tanaka, N. Brefuel, H. Cailleau, J. F. Letard, S. Ravy, P. Fertey,M. Nishino, S. Miyashita and E. Collet, Phys. Rev. B 93, 014419 (2016).[32] N. Willenbacher and H. Spiering, J. Phys. C 21, 1423 (1988). H. Spiering and N. Willenbacher, J. Phys.: Condens. Matter1, 10089 (1989). H. Spiering, Top. Curr. Chem. 235, 171 (2004).[33] S. Nosé, J. Chem. Phys. 81, 511 (1984). W. G. Hoover, Phys. Rev. A 31, 1695 (1985).[34] C. Enachescu and W. Nicolazzi, Comptes Rendus Chimie 21, 1179 (2018).[35] M. Nishino and S. Miyashita, Phys. Rev. B 91, 134411(1-13) (2015).[36] S. Miyashita, T. Shirai, T. Mori, H. De Raedt, S. Bertaina and I. Chiorescu,, J. Phys. B Atomic Mol and Opt Phys, 45,124010 (2012). T. Shirai, S. Todo, H. De Raedt, and S. Miyashita, Phys. Rev. A98, 043802 (2018). N. Lambert, et al.,Phys. Rev. Research 5, 013181 (2023), and X. Wang et al. AVS Quantum Sci. 5, 043801 (2023).