# Dynamical Properties of Spin-Crossover Solids Within the Kinetic Spin-1 BEG Model in the Presence of a Time-Dependent Magnetic Field

Saliou Bolarinwa Ogou, Djidjoho Toussaint Oke, Félix Hontinfinde, Kamel Boukheddaden

► To cite this version:

Saliou Bolarinwa Ogou, Djidjoho Toussaint Oke, Félix Hontinfinde, Kamel Boukheddaden. Dynamical Properties of Spin-Crossover Solids Within the Kinetic Spin-1 BEG Model in the Presence of a Time-Dependent Magnetic Field. Advanced Theory and Simulations, 2019, 2 (5), pp.1800192. 10.1002/adts.201800192. hal-03043796

HAL Id: hal-03043796

https://hal.science/hal-03043796

Submitted on 25 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamical Properties of Spin-Crossover Solids Within the Kinetic Spin-1 BEG Model in the Presence of a Time-Dependent Magnetic Field

Saliou Bolarinwa Ogou, Djidjoho Toussaint Oke, Félix Hontinfinde, and Kamel Boukheddaden

Spin-crossover (SCO) and Prussian blue analogs (PBAs) materials are investigated in 2D with a three-state Blume–Emery–Griffiths (BEG) model where each spin interacts with its nearest neighbors (nn) and may be either in high-spin (HS) or low-spin (LS) state. The interactions through the system lattice are temperature-dependent to account for spin-phonon interactions. The system is also in contact with an oscillating magnetic field energy. The generated numerical results by the dynamic mean field theory (DMFT) study approach are consistent with those derived by kinetic Monte Carlo (KMC) simulations with Glauber dynamics and Arrhenius transition rates. First-order transitions with thermally induced hysteresis phenomena have been observed. Near the hysteresis loops, the model exhibits throughout relaxation curves, some fluctuations in the LS phase, strengthened by increasing temperature where this phenomenon becomes temperature- and magnetic field-dependent.

## 1. Introduction

Because of the needs for our companies in the field of the treatment and storage of informations, transportation of art-works or other valuable objects, evaluation of the degree of alteration of materials after collisions,[1] the design of miniaturized devices with fast response is a major stake. Switchable molecular materials with thermodynamic bistability are one of the solutions to this requirement. Among possible systems, spin-crossover (SCO) materials with a central metal ion with the property to switch between two different spin states, constitute

serious and original candidates. Their lay-out leads interactions between the compo-nent molecules of materials, which cause the appearance of unique or multiple phys-ical properties (conduction, magnetism, photochromism, nonlinear optics ...).[2–6] Furthermore, their molecular aspect makes them easily flexible chemically to optimize their properties or to combine them with other molecular systems to generate multi-functional materials. Therefore, under the influence of various stimuli such as light, pressure, temperature, magnetic and electric fields, etc.[7–10] phase transitions be-tween the low-spin (LS) diamagnetic state and the high-spin (HS) paramagnetic state are triggered.[7,9,10] The change of induced

unit of memory can be reduced to a molecule, thus allowing to reach storage capacity more important than those of conventional materials. The thermally induced spin transition leads to both electric and structural changes, often observed as a color and magnetic moment changes.[2,11,12] When the interactions between molecules are weak, the HS fraction changes smoothly with the temperature, whereas when it becomes strong, the system exhibits cooperative phenomena,[13–15] which manifest through the existence of first-order phase transitions accompanied with thermal hysteresis. An interesting example is that of the [Fe(NH<sub>2</sub>tr z)<sub>3</sub>](NO<sub>3</sub>)<sub>2</sub>, (NH<sub>2</sub>tr z = 4 – amino – 1, 2, 4 – tr iazole) which exhibits a spin transition with a hysteresis loop near the room temperature region.[16] Then, the change in HS fraction becomes sharper and sharper with increasing interaction strength between molecules. Of course, the interaction in SCO solids is dominated by the variations of unit-cell volume and bond length, that are considerably larger in the HS state. These result, in addition to the larger electronic degeneracy, also in a larger phonon density for this state.[13,14,17] At the atomic scale and in the case of Fe(II), the SCO phenomenon is the result of the redistribution of the electrons between the bonding t<sub>2g</sub> and the antibonding e<sub>g</sub> orbitals. In the diamagnetic (σ = 0) LS state, only bonding orbitals are populated (t<sub>2g</sub> e<sub>g</sub><sup>0</sup>), while in the paramagnetic (σ = 2) HS state, the electronic configuration becomes (t<sub>2g</sub><sup>4</sup> e<sub>g</sub><sup>2</sup>). As demonstrated in several works,[13,14,18] the elastic in-teractions are at the hearth of the existence of cooperative effects in SCO materials and play a crucial role in the existence of the

Dr. S. B. Ogou, Dr. D. T. Oke, Prof. F. Hontinfinde
Département de Physique (FAST) et Institut des Mathématiques et de Sciences Physiques (IMSP)
Université d'Abomey-Calavi
01 BP 613 Porto-Novo, Benin
E-mail: felix.hontinfinde@fast.uac.bj

Dr. D. T. Oke, Prof. K. Boukheddaden
Groupe d'Etudes de la Matière Condensée
Université Paris-Saclay, Université de Versailles/St. Quentin en Yvelines-CNRS
45 Avenue des Etats Unis, F78035 Versailles Cedex, France

<figure><img src="image_2.png" /></figure>

The ORCID identification number(s) for the author(s) of this article can be found under https://doi.org/10.1002/adts.201800192
first-order transitions and the thermally induced hysteresis loops observed experimentally.[13] Experimental results allow to establish that SCO transitions involve both electronic transformation (spin and orbital) and structural modifications.[18–23]

From the theoretical point of view, many interesting works have addressed static and dynamic properties of SCO materials. Most of them are based on Ising-like models,[24–26] atom-phonon coupling descriptions,[12,27–32] or more recently on complex elastic descriptions[13,23,33–36] or taking into account the volume change at the transition. These models are solved analytically using mean-field approaches or by Monte Carlo simulations[13,28,33] or using molecular dynamics.[14,37]

In the present investigations, the SCO is described using the Blume–Emery–Griffiths (BEG) model that allows to include not only the interactions between the SCO sites but also to take into account, possible magnetic interactions or magnetic field effect in the high-spin state. In the past, several methods were used to solve the Ising spin-1 model under an oscillating magnetic field[38–44] or mixed spins kinetic models,[45–47] which exhibited interesting physical properties. In this work, we consider the BEG model for SCO or PBAs[24–36] and crossover solids under oscillating magnetic field to describe their thermal and time-dependent properties. As in our previous works,[48–50] the quadrupolar coupling parameter K between SCO units is assumed to depend linearly on the absolute temperature T in the form, $K = \alpha k_B T$, in order to take into account acoustic phonons origin of the interactions between SCO units. The effective “ligand-field” energy D depends on the absolute temperature and contains the combined effects of the ligand-field strength Δ and that of the degeneracy ratio g between LS and HS states which result in an entropic term, stabilizing the HS state at high temperature. The model is solved in the framework of the DMFT approach. One gets the motion equations of the order parameters $m = \langle \sigma \rangle$ (magnetization) and $n_{HS} = \langle \sigma^2 \rangle$ (HS fraction), which describe nonequilibrium properties of the system. On the other hand, dynamical properties of the model are computed by using KMC simulations[48] with Glauber dynamics, where the time scale is given with suitable Arrhenius transition rates.[49–51] The relaxation of HS fraction shows oscillating nonlinear shapes[33,52] accompanied by some fluctuations in the stationary state, in which a residual HS fraction remains due to the effect of the magnetic field. The obtained results of DMFT simulations are discussed in relation with those by KMC calculations. Analysis of spatio-temporal configurations have been performed to explain the spin transition during the relaxation process.

The paper is organized as follows. Section 2 is devoted to the presentation of the model and the description of the used calculation methods: DMFT and KMC simulations. In Section 3, we present and discuss the obtained results. Section 4 contains the conclusion and outlines some possible developments of this work.

## 2. Hamiltonian Model

The Hamiltonian of spin-1 BEG model adapted for SCO materials with magnetic interactions and applied magnetic field is written as follows:

$$H = -J \sum_{\langle i,j \rangle} \sigma_i \sigma_j - K \sum_{\langle i,j \rangle} \sigma_i^2 \sigma_j^2 + D \sum_i \sigma_i^2 - h(t) \sum_i \sigma_i \quad (1)$$

where $\sigma_i = \pm 1, 0$ are fictitious spin values located at site i of the square lattice. The spins $\sigma_i = \pm 1$ describe the magnetic HS spin state and $\sigma_i = 0$ is associated with the diamagnetic LS state. According to the 2D character of the system, and to the square symmetry of the lattice in which only nearest-neighbor (nn) interactions are considered, the coordination number is $z = 4$. The magnetic interactions between the magnetic states, $\sigma_i = \pm 1$, are taken into account through the exchange term J and the quadrupolar interactions (between the SCO sites) are introduced through K. Due to the elastic nature of the SCO transition,[33–37] this K term that takes into account the phonon contribution which also originates from the intra-molecular vibrations, lattice distortion, and/or acoustic phonons, provides the elastic long-range interactions between the SCO units, and is written here as $K = \alpha k_B T^{[48–50]}$ with the ratio $\gamma = J/K$ taken as a tunable parameter. $D = \Delta - k_B T \ln(g)$ is the effective ligand-field strength and h is the energy associated with the application of an external magnetic field. In the following, a radio frequency magnetic field $h(t) = h_0 \cos(\omega t)$ is considered where $\omega$ is the oscillation frequency.

### 2.1. Dynamic Mean-Field Theory (DMFT) Approach

Throughout the spin lattice, each spin site i feels the following one-site Hamiltonian, $H_i$, in the mean field approximation

$$H_i = -zJ m \sigma_i - zK n_{HS} \sigma_i^2 + D \sigma_i^2 - h(t) \sigma_i \quad (2)$$

where $m = \langle \sigma_i \rangle$ and $n_{HS} = \langle \sigma_i^2 \rangle$ are considered as invariant by translation over the lattice.

The associated mean-field free energy per site is given by

$$F(m, n_{HS}, T) = U - TS \quad (3)$$

where U and S are the internal energy and entropy per site of the system, respectively, given by

$$U = -\frac{z}{2} J m^2 - \frac{z}{2} K n_{HS}^2 + D n_{HS} - h(t) m \quad (4)$$

$$S = k_B \beta \langle H_i \rangle + k_B \ln \sum_{\sigma} e^{-\beta H_i} \quad (5)$$

Here, $\beta = \frac{1}{k_B T}$ ($T$ is the temperature) and $\langle H_i \rangle = -zJ m^2 - zK n_{HS}^2 + D n_{HS} - hm$ is the average value of $H_i$. After some simple calculations, the variational free energy is given by

$$F(m, n_{HS}, T) = \frac{z}{2} J m^2 + \frac{z}{2} K n_{HS}^2 - k_B T \ln[1 + 2e^{\beta(zK n_{HS} - D)} \cosh \beta(zJm + h(t))] \quad (6)$$
### 2.1.1. The Dynamic Choice

The dynamical properties of Hamiltonian 1 are investigated in the frame of a microscopic master equation which governs the time evolution of the probability $P(\{\sigma_i\}, t)$ to occupy the spin configuration $\{\sigma\}^{[49-51,53]}$ at time $t$. The flux of probability accounts for transitions from the configuration $\{\sigma\}_{i \rightleftharpoons \sigma'}_{i}$ with transition rates $W$. Following the Glauber-type stochastic dynamics, the master equation reads

$$ \frac{\partial P(\{\sigma\}, t)}{\partial t} = -\sum_{i=1}^{N} W_i(\sigma_i \rightarrow \sigma_i') P(\{\sigma\}_{i}, \sigma_i, t) + \sum_{i=1}^{N} W_i(\sigma_i' \rightarrow \sigma_i) P(\{\sigma\}_{j}, \sigma_i', t) \quad (7) $$

where $W_i(\sigma_i \rightarrow \sigma_i')$ is the transition rate of the $i$th spin from the value $\sigma_i$ to $\sigma_i'$ defined by

$$ W_i(\sigma_i \rightarrow \sigma_i') = \frac{1}{3\tau} \frac{e^{-\beta H_i}}{\sum_{\sigma_i'} e^{-\beta H_i}} \quad (8) $$

$\tau$ stands for the Arrhenius time scale and $\frac{1}{\tau}$ denotes the effective intramolecular frequency associated with the “spontaneous spin reversals”$^{[49-51]}$

$$ \frac{1}{\tau} = \frac{1}{\tau_0} e^{-\beta E_0^a} \quad (9) $$

where $1/\tau_0$ is the “intrinsic” frequency spin-flip process between HS and LS states, taken as constant and $E_0^a$ denotes the intramolecular vibronic energy barrier.

In the equilibrium state, the probabilities satisfy the detailed balance condition

$$ \frac{W(\sigma_i \rightarrow \sigma_i')}{W(\sigma_i' \rightarrow \sigma_i)} = \frac{P(\{\sigma_i\}, \sigma_i', t)}{P(\{\sigma_i\}, \sigma_i, t)} \quad (10) $$

Thus, the probability per unit time is given by

$$ W_i(\sigma_i \rightarrow \sigma_i') = W_i(\sigma_i') = \frac{e^{-\beta H_i}}{3\tau \left(1 + 2e^{\beta(zK_{HS} - D)} \cosh \beta(zJm + h)\right)} \quad (11) $$

Thus, these transition rates $W$, must fulfill the detailed balance condition.

### 2.1.2. The Motion Equations

The mean value of spin $\sigma_k$ and $\sigma_k^2$, respectively, associated to the magnetization and HS fraction calculations are

$$ m = \sum_{\{\sigma\}} \sigma_k P(\{\sigma_i\}, t) \quad \text{and} \quad n_{HS} = \sum_{\{\sigma\}} \sigma_k^2 P(\{\sigma_i\}, t) \quad (12) $$

and $\sum_{\sigma_k' \neq \sigma_k} W_k(\sigma_k') = \frac{1}{3\tau}$. After some calculations, one gets the time evolution of the magnetization $m$ and the fraction $n_{HS}$ in

the forms

$$ \left\{ \begin{aligned} \frac{dm}{dt} &= \frac{e^{-\beta E_0^a}}{3\tau_0} \left( -m + \frac{2 \sinh \beta(zJm + h(t))}{e^{-\beta(zK_{HS} - D)} + 2 \cosh \beta(zJm + h(t))} \right) \\ \frac{dn_{HS}}{dt} &= \frac{e^{-\beta E_0^a}}{3\tau_0} \left( -n_{HS} + \frac{2 \cosh \beta(zJm + h(t))}{e^{-\beta(zK_{HS} - D)} + 2 \cosh \beta(zJm + h(t))} \right) \end{aligned} \right. \quad (13) $$

These coupled nonlinear differential equations (13) are numerically solved by using fourth order Adams–Moulton predictor-corrector methods.$^{[54-57]}$ Thermal properties are got by integrating on $[0, 2\pi]$ in time. Then, the dynamic order parameters are defined as the time-average magnetization $m$ and $n_{HS}$ fraction over a period of the oscillating magnetic field energy

$$ \left\{ \begin{aligned} m &= \frac{1}{2\pi} \int_{0}^{2\pi} m(\xi) d\xi \\ n_{HS} &= \frac{1}{2\pi} \int_{0}^{2\pi} n_{HS}(\xi) d\xi \end{aligned} \right. \quad (14) $$

where $\xi = \omega t$. Here, $\omega$ becomes a simple parameter and we considered $\xi$ as the control parameter in time. The method used to solve the system equations (14) is that of Romberg$^{[54-57]}$ with Richardson extrapolation to obtain thermodynamic quantities such as: magnetization $m$, fraction $n_{HS}$, and the magnetic susceptibility given by

$$ \chi = \frac{n_{HS} - m^2}{k_B T} \quad (15) $$

Within the system dynamics, Van Hove-like equations$^{[58]}$ give the same system equations through the free energy per site. One gets

$$ \left\{ \begin{aligned} \frac{dm}{dt} &= -\frac{1}{3zJ\tau} \frac{\partial F}{\partial m} \\ \frac{dn_{HS}}{dt} &= -\frac{1}{3zJ\tau} \frac{\partial F}{\partial n_{HS}} \end{aligned} \right. \quad (16) $$

## 2.2. Kinetic Monte Carlo Simulations

Here, the $n$-fold algorithm developed by Boltz, Kalos, and Lebowitz (BKL)$^{[59]}$ is used to investigate the model properties. In the Monte Carlo simulations$^{[60,61]}$ procedure, an initial configuration $\{\sigma\}$ (of linear size $L$) of the square system is chosen with periodic boundary conditions. Then, the total number of possible spin-flip processes that lead to another configuration $\{\sigma'\}$ is calculated. Let us denote by $W(\{\sigma\} \rightarrow \{\sigma'\})$ the transition rate from $\{\sigma\}$ to $\{\sigma'\}$. For the Glauber spin-flip dynamics, $W(\{\sigma\} \rightarrow \{\sigma'\})$ becomes

$$ W(\{\sigma\} \rightarrow \{\sigma'\}) = \frac{1}{3\tau(1 + e^{\beta\Delta E})} \quad (17) $$

where $\Delta E = E(\{\sigma'\}) - E(\{\sigma\})$ denotes the change in the system energy associated to the spin-flip move and $\frac{1}{\tau}$ is defined as in
Equation (9). Then, one calculates the total evolution rate $Q(\{\sigma\})$ of $\{\sigma\}$ by considering all possible processes

$$Q(\{\sigma\}) = \sum_{a=1}^{N} W_a \qquad (18)$$

where $a$ stands for a given possible spin-flip process. Two random numbers $0 \le r_1, r_2 \le 1$ are chosen to calculate the lifetime $t(\{\sigma\}) = -\ln(r_1)/Q(\{\sigma\})$ of configuration $\{\sigma\}$ and a random evolution rate $\eta(\{\sigma\}) = r_2 Q(\{\sigma\})$. The first process $a$ where the following condition

$$\sum_{a=1}^{b} W_a > \eta(\{\sigma\}) \qquad (19)$$

is fulfilled is realized with probability 1. The total evolution time $\bar{t}$ is given by $\bar{t} = \sum_{\{\sigma\}} t(\{\sigma\})$ where the sum runs over all spin configurations generated up to the steady state. After a suitable number $n_{MC}$ of Monte Carlo steps, stationary values of physical quantities of interest are calculated. An adapted trapezoidal integration method is used after relaxation of physical quantities in real time. Furthermore, 20 – 25 independent runs are considered in the averaging procedure.

# 3. Results and Discussions

## 3.1. Thermal Equilibrium Properties of the System in Mean-Field Approximation

First, in the stationary case, the system's Hamiltonian is analyzed without the external magnetic field ($h_0 = 0$ K) and without magnetic exchange interaction ($\gamma = 0$) as well. In this case, the model is solved in the framework of mean-field theory (MFT). Then, the system is reduced to two-dependent equations, solved numerically. By raising the temperature sequentially from 20 to 80 K with step $\Delta T = 0.1$ K, some thermal behaviors are got and displayed in Figure 1 for selected values of model parameters. Figure 1a displays the thermal behaviors of $n_{HS}$ fraction for specific values of $\alpha$. When $\alpha \ge 1$, first order transitions take place (see Figure 1b). Then, the reduced equilibrium temperature is a decreasing function of $\alpha$. Up to $\alpha = 1.05$, the analytical expression of $T_{eq}/\Delta^{[48]}$ coincides with the one obtained by the DMFT approach and discrepancy only appears beyond. Figure 1c shows the behavior of the corresponding magnetic susceptibility with temperature and Figure 1d depicts the $\alpha$-dependence of the free energy per site for selected values of $\Delta$. The crossing point, $F^{eq} = 0$, corresponds to the onset of the first-order transition as already obtained in ref. [49]. Figure 2 displays two-temperature dependencies of the free energy. With increasing temperature, at each temperature ($T^*$) and equilibrium temperature $T_{eq}$, the free energy is a decreasing function. Up to $T_{eq} = T^*$, the transition is of first-order and at high temperature, gradual spin-conversion occurs. For weak free energy, the phase is of diamagnetic (LS phase), and between $T_{eq}$ and $T_c$ is of paramagnetic. The ferromagnetic phase prevails when the free energy becomes higher than in other phases.

<figure><img src="image_3.png" /><figcaption>Figure 1. Some features of the model as functions of temperature $T$ and parameter $\alpha$ in mean-field approximation. a,c) Gradual and first-order transitions are depicted for selected values of $\alpha$ and fixed other parameters. The curves shift to the left when $\alpha$ increases and jumps appear in their behaviors. b) The reduced equilibrium temperature $T_{eq}/\Delta$ as function of $\alpha$ obtained by two methods (see text). d) Equilibrium free energy as function of $\alpha$. The point $F_{\alpha=1.05}^{eq} = 0^{[49]}$ is the onset of first-order transitions for selected values of $\Delta$ (400, 450, and 500 K). $n_{HS}$ and $m$ are solutions at thermodynamic equilibrium, that is, corresponding to the minimum of the free energy. Other parameters values are written in the panels.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. Thermal behavior of the free energy per site through different phase transitions. At low temperatures, free energy is obtained of $T_{eq} = T^*$. Beyond that, values of the free energy are quite different and paramagnetic phase is favored between them. Over $T^*$, ferromagnetic phase prevails. Under the equilibrium temperature, the phase is diamagnetic (LS phase). Other values of used parameters are $\gamma = 0$, $\Delta = 400$ K, $h_0 = 0$, K, and $g = 100$.</figcaption></figure>

## 3.2. Nonequilibrium Properties of the System in the DMFT

### 3.2.1. Temperature-Dependent System and Phase Diagram

The effects of external magnetic field constraint on spin-crossover compounds have been thoroughly studied in the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Time-dependent behaviors of magnetization *m* and *n*<sub>HS</sub> fraction for selected values of *α* and *γ*. a,b) For *γ* = 0, the effect of *α* shows two-dependent temperature regimes: for weak interactions in contribution of phonon, “aperiodic” regime appears and when these interactions are strong enough, the regime is of “pseudo-periodic”. The same behaviors are obtained at *α* = 0.5 and varying values of the parameter *γ*. Other values of model parameters considered are: *T* = 40 K, *Δ* = 400 K, *g* = 100, and *h*<sub>0</sub> = 150 K.</figcaption></figure>

literature.[2,62,63] Static and pulsed fields have been considered. In most cases, small effects are got and a high magnetic field (30 – 100 T) is required to obtain meaningful results.[64–66] At so high field, it is not easy to measure experimentally the magnetic properties and experiments with optical reflectivity detection technique are often used.[62,63] To our knowledge, besides pulsed time-dependent magnetic fields, experiments with time-dependent sinusoidal fields are lacking in the literature. Motivated by the aim to generate qualitative and quantitative predictions on this interesting case, a sinusoidal field with a high amplitude *h*₀ = 150 K is considered in the present work. The individual spin-flip frequency is set to the value 1/3τ₀ = 1s⁻¹ and the activation energy to E⁰ = 0 K. The relaxation of the magnetization *m* and the HS fraction *n*<sub>HS</sub> is obtained in time for selected values of *α* and *γ* (see **Figure 3**). These results are numerically obtained at fixed temperature (T = 40 K) and for initial conditions (±1, 0). Two regimes appear when increasing *α* values (Figures 3a,b): “pseudo-periodic” for high values of *α* and “aperiodic” for low *α* values. The same phenomena have been observed for increasing *γ* (see Figures 3c,d). One can suspect that in the model, *γ* and *α* play a leading role on the amplitude of the oscillations obtained during the relaxation. Now, we could notice that according to *γ* and *α* values, the system may undergo a transition of the first-order, second-order, or gradual spin-conversion.[48,49] A deep analysis of the relaxation in the vicinity of the thermal hysteresis loop could help to know about such a behavior of the model system. Thus, we display the results in the next section where the relaxation is made near the hysteresis loop at low temperatures. After that, the thermal behaviors are obtained on *m* and *n*<sub>HS</sub> for some selected parameter values (**Figure 4**). First-order spin transitions characterized

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Some thermal behaviors of the model as function of parameter *α*. Gradual spin-conversion and first-order transition are depicted for selected values of *α* and fixed other parameters. a) magnetization *m*, b) fraction *n*<sub>HS</sub>, and c) magnetic susceptibility *χ*. The curves shift to the left for low temperatures when *α* increases. d) The reduced equilibrium temperature as function of *α* at fixed *z*. Other values of used parameters are: *γ* = 0, *Δ* = 400 K, *h*₀ = 150 K, and *g* = 100.</figcaption></figure>

by jumps in *m* and *n*<sub>HS</sub>, occur at around 40 K when raising consequently *α* values (see Figure 4a–c). Such a low transition temperature range has been observed during the investigation of the Mn<sup>111</sup>[(pyrol)<sub>3</sub>tren] in the presence of a low magnetic field of the order of 23 T.[67] The (pyrol)<sub>3</sub>tren is the trianionic Schiff base obtained from the condensation of pyrolle – 2 – carboxaldehyde with 2, 2', 2"tris(ethylamino)amine. The reduced equilibrium temperature *T*<sub>eq</sub>/Δ is a decreasing function of *α* as already obtained in refs. [48,49] In **Figure 5**, the same tendencies are obtained for varying values of *γ* at fixed values of *α* with, however, a clear saturation of the net magnetization *m* at high temperature. Here, the HS units are considered to be created and interact magnetically in the system. Phase boundaries are determined, and the corresponding phase diagrams constructed (**Figure 6**). The second-order phase boundaries are obtained when the total lattice magnetization vanishes. At the first-order transition temperatures, both *m* and *n*<sub>HS</sub> show discontinuities in their behaviors. Corresponding discontinuities are observed in the behavior of the magnetic susceptibility *χ*. The positions of the TCP (tricritical point) differ slightly from those found by the corrected effective field theory (CEFT)[48] and by exact recursion relations on the Bethe lattice (BL).[49] Positions of the TCP for varying values of *Δ* indicate the existence of tricritical lines. BL approach, CEFT, and DMFT methods yield similar temperature phase diagrams. But by the DMFT approximation, at relatively low *γ* values (*γ* < *γ*<sup>T</sup>), two temperatures may exist: *T*<sub>eq</sub> and *T*<sub>c</sub> with *T*<sub>eq</sub> < *T*<sub>c</sub>. For this range of temperatures, dia-para- and para-ferro-magnetic phase transitions can be possible and the HS fraction reaches its saturated value, *n*<sub>HS</sub> = 1. At the vicinity of the first-order transitions, we displayed in **Figure 7** the associated thermal hysteresis loops for some selected values of model parameters. When raising the temperature from 15 to 50 K with 0.1 temperature step, for fixed *γ* = 0 and *α* = 3, the equilibrium temperature is about *T*<sub>1/2↑</sub> = 35.80 K on the upward branch and *T*<sub>1/2↓</sub> = 24.90 K on the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Some thermal behaviors of the model as function of parameter $\gamma$. For $\gamma \neq 0$, second- and first-order transitions occurred when increasing $\gamma$ values. a) Magnetization $m$, b) fraction $n_{HS}$, and c) magnetic susceptibility $\chi$. The curves shift to the left again for low temperatures when $\gamma$ increases as already obtained.[48,49] Other values of model parameters are considered: $\alpha = 0.5$, $\Delta = 400$ K, $h_0 = 150$ K, and $g = 100$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Thermal phase diagram of the model in $(\gamma, T_c)$ plane for selected values of $\alpha$. Three phases are found as in refs. [48,49]: dia-para- and ferromagnetic phases. The triangle point is the tricritical point (TCP) at $\gamma_t$. These values tend to be zero when the phonon contribution is important. Tricritical line is obtained for varying $\Delta$ values (dots line in d). Dashed lines correspond to critical temperatures $T_c$ and continuous lines are associated to equilibrium temperature $T_{eq}$. Other values of model parameters are written in the panels.</figcaption></figure>

backward branch (Figures 7a,b). When $\gamma = 2$ and $\alpha = 1$, the width of the hysteresis loop is slightly reduced (5.6 K in order) but the whole cycle is shifted to relatively high temperatures, with 41.80 K and 36.2 K for the upward and backward transition temperatures, respectively (see Figure 7c).

### 3.2.2. Isothermal Relaxation Near the Hysteresis Loops under Oscillating Field

It is interesting to investigate the relaxation properties of the system according to the model parameters. The method used is based on kinetic Monte Carlo (KMC) simulations, using the BKL algorithm which was described in the previous section. A sample of $L \times L = 100 \times 100$ lattice sites is considered and some of the obtained results are compared with those by DMFT approach. The relaxation processes take place from the HS metastable state to the LS stable state at fixed temperature. The final data are averaged over 25 independent simulation runs. As reported in **Figure 8** for both methods, one gets from the initial stage, concave curves of the order parameters in real and arbitrary time which somewhat coincide in the LS phase. For $\gamma = 0$ and temperature set to $T = 40$ K, the system relaxes slowly by the DMFT method compared to that of KMC simulations for varying values of $\alpha$. This can be related to an invariant space introduced to calculate the mean value of spin $s_i$ at site $i$. The relaxation tail presented some oscillations around 0 (see Figure 8a,c).
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Thermal hysteresis loops of fraction $n_{HS}$ and magnetization $m$ at the vicinity of first-order transitions. a) and b) are in same captions: $\alpha = 3$, $\gamma = 0$, $\Delta = 400$ K, $g = 100$, and $h_0 = 150$ K. In panel c): $\alpha = 1$ and $\gamma = 2$, the loop area is reduced and obtained at relatively high temperatures. The red down arrow indicates the starting point of isothermal relaxation.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Compared results of isothermal relaxation curves from HS metastable state to LS state and for the magnetization $m$ by both calculation methods and for selected values of $\alpha$ and $\gamma$. Oscillations are visible in the LS phases and become important in magnitude with the KMC method (red curves). Faster is the relaxation with KMC than by the DMFT (black curves) approach where the system takes a long time accompanied by some fluctuations from initial stage. a,b) $\gamma = 0$ and $\alpha = 0$; 1. c,d) $\alpha = 0.5$ and $\gamma = 0$; 1. Other values considered for model parameters are: $T = 40$ K, $\Delta = 400$ K, $h_0 = 150$ K, and $g = 100$.</figcaption></figure>

for both order parameters. In other panels (Figure 8b,d), $n_{HS}$ decreases in time with some undulations from initial stage when the system reaches the LS phase. It is interesting to notice that while the magnetization oscillation frequency is exactly that of the applied field $h(t) = h_0 \cos(\omega t)$, the double frequency is obtained for $n_{HS}$ as it clearly appears in Figure 8b,d. The main reason of this frequency doubling has to be related to the quadratic nature of $n_{HS} = (\sigma^2)$. Now for the validity limits of the model, the relaxation curves are integrated in time by an adapted trapezoidal method at each temperature and some KMC results are compared to those obtained by DMFT (Figure 9). For gradual transitions, all results are similar. Discrepancy appears for high interaction strength, that is, in the first-order transition domain (Figure 9a,b). The obtained results are averaged over 20 independent numerical experiences. The hysteresis loops are consequently obtained for several values of $\alpha$ and $\gamma$ (Figure 9c,d). At $\gamma = 0$ and from $\alpha = 3$ to $\alpha = 4$, the equilibrium temperatures $T_{eq}$ shift to low temperatures as already obtained in our previous works.[48,50]

Let us examine the relaxation phenomena near the thermal hysteresis loops of Figure 9c. The “intrinsic” frequency flip-flip process between HS/LS states is $\frac{1}{3\tau_0} = 10^2 s^{-1}$ and the intramolecular vibronic energy barrier $E_0^2 = 70$ K. In Figure 10a,b, we notice that the relaxation curves are sensitive to the temperature in both calculation methods (e.g., $T = 23.9$ and $T = 24$ K). The oscillations appear near the loops’ area, except for low-temperature from the loop (Figure 10a,b for KMC results). Near the loop, fluctuations become important and the system
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Some thermal behaviors of the system as functions of model parameters by both methods. Agreement results are observed in gradual spin transition domains (low values of α or γ). a,b) The results are totally different in the first-order transition domain (by increasing α or γ values). Corresponding hysteresis loops are displayed in c) and d). The loop area increases when interactions become important and shifts to low temperatures. The red down arrows indicate the starting point of isothermal relaxation. Other parameters are written in panels.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Isothermal relaxation curves near hysteresis loops. a,b) The temperature is sensitive in the relaxation processes for both methods ($T = 23.9$ K and $T = 24$ K). c,d) Near the loops, the system oscillates in the HS and LS phases when the relaxation occurs at $T = 27.4$ K and $T = 24$ K, respectively. Some snapshots are displayed in next panels (Figures 11 and 12). Values considered for other parameters are: $\Delta = 400$ K, $g = 100$, $\frac{1}{3\tau_0} = 10^2 s^{-1}$, and $E_0^a = 70$ K.</figcaption></figure>

oscillates in the HS and LS phases where some spins ($\sigma = -1$) are favored in the process at $T = 27.4$ K (Figure 10c) and then, spin-transition occurs for trending spin values $\sigma = 0$. Then, the induced magnetization $m$ has fluctuations through negative values in the LS phase. Hence, the temperature is an important factor in this range. In Figure 10d, at $T = 22$ K near the loop (for $\alpha = 4$. and $\gamma = 0$) similar situation is presented. Further investigations of the dynamics are needed to substantiate the role of the magnetic field in the course of the relaxation. For a better understanding of these phenomena, some associated snapshots of spin configurations are displayed in Figures 11 and 12 obtained at 24 and 27.4 K, and show that throughout the relaxation processes, LS units are formed and get bigger to coalesce with some few HS units. One gets, at time $t_1 = 1.55$ s, 10% of compact LS nucleated units. At time $t_2 = 1.65$ s (Figure 10b), similar behaviors are

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Snapshots of the system configuration near the hysteresis loops at $T = 24$ K at different relaxation times: $t_1 = 1.55$; $t_2 = 1.62$; $t_3 = 1.82$; and $t_4 = 1.92$ s. White clusters correspond to those of LS units whereas red (spin +1) and the blue (spin -1) dots ones consist of HS unit clusters. $n_{LS}$ increases in the course of the time. In a–d), respectively, 10%, 15%, 40%, and 55% of the lattice are occupied by LS units. Values considered for other parameters are: $\alpha = 3$, $\gamma = 0$, $\Delta = 400$ K, $g = 100$, $\frac{1}{3\tau_0} = 10^2 s^{-1}$, and $E_0^a = 70$ K.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Snapshots of the system configuration near the hysteresis loop at $T = 27.4$ K at different relaxation times: $t_5 = 1.90$; $t_6 = 1.96$; $t_7 = 2.22$; and $t_8 = 5.01$ s. Corresponding percentages of LS units are also: a) 10%, b) 15%, c) 40%, and d) 55%. Values considered for other parameters are: $\alpha = 3$, $\gamma = 0$, $\Delta = 400$ K, $g = 100$, $\frac{1}{3\tau_0} = 10^2 s^{-1}$, and $E_0^a = 70$ K. Colors have the same meaning as in Figure 11.</figcaption></figure>
observed for the fixed parameters values used with the presence of some few LS units at their edge, as well as sparse small LS domains inside the lattice. When we approach the loop enough, the time increases and the nucleation phenomenon is emphasized (see figures for more details). LS aggregates are rapidly formed everywhere and become sparse, surrounded with some HS units of spin values $\sigma = -1$. Through these panels, we have concluded that in the relaxation, we have denoted: magnetic relaxation process in the HS phase and spin-transition between para- and diamagnetic phases with competition of temperature and magnetic field.

These studies allow us to substantiate the features of LS phase growth involved during the relaxation process of an SCO system under oscillating magnetic field. Two interacting remarks emerge. The first one is the key role played by the temperature T as a driving force which induces the thermal spin-transition. The second one is that of the magnetic-transition played by the oscillating magnetic-field when magnetic energy propagates within the system.

## 4. Conclusion

The investigation of time-dependent systems has attracted much interest in recent years not only for their fundamental aspects but also for their applicability, such as quantum transport, quantum optics, quantum information, spintronics. In this work, the theoretical study of a 2D spin-1 BEG model under the constraint of a time-dependent oscillating magnetic-field energy is performed to investigate SCO and PBAs materials by means of dynamic mean field theory (DMFT) calculations and kinetic Monte Carlo (KMC) simulations with Glauber dynamics. The fascinating results obtained by both methods bear some resemblance and are obtained according to Arrhenius' time scale. Usually, gradual and first-order transitions are obtained with varying values of model parameters in the phase diagram. Around first-order transitions, the system exhibits hysteresis loops whose area depends on the interactions strength as already obtained in our previous works. Near the loop, the interactions increase and thermal fluctuations are important too, causing oscillations in HS and LS phases where the magnetic-field strength and the temperature play a crucial role. The relaxation process from metastable HS state occurs with two processes: magnetic-transition in HS phase and spin-transition between para- and diamagnetic phases. The analysis of the relaxation process strengthens the transitions that occurred through the stochastic distribution of the LS island species. One could think that the relaxation processes follow a 2D-nucleation as expressed in ref. [50].

## Conflict of Interest

The authors declare no conflict of interest.

## Keywords

2D-nucleation, dynamic mean field approach, kinetic Monte Carlo simulations, phase diagram, relaxation dynamics, spin-crossover compounds

[1] C. M. Jureschi, J. Linares, A. Rotaru, M. Ritti, M. Parlier, M. M. Dîrtu, M. Wolff, Y. Garcia, *Sensors* **2015**, *15*, 2388.

[2] P. Gütlich, H. A. Goodwin, *Spin-Crossover in Transition Metal Compounds I, II, III*, Springer, Berlin **2004**, pp.233–235.

[3] J. H. Ammeter, *Nov. J. Chem.* **1980**, *4*, 631.

[4] S. Ohkoshi, K. Hashimoto, *J. Am. Chem. Soc.* **1999**, *121*, 10591.

[5] S. Ohkoshi, S. Ikeda, T. Hozumi, T. Kashiwagi, K. Hashimoto, *J. Am. Chem. Soc.* **2006**, *128*, 5320.

[6] S. Ohkoshi, K. Imoto, Y. Tsunobuchi, S. Takano, H. Tokoro, *Nat. Chem.* **2011**, *3*, 564.

[7] H. Tokoro, S.-I. Ohkoshi, K. Hashimoto, *Appl. Phys. Lett.* **2003**, *82*, 1245.

[8] F. Varret, K. Boukheddaden, C. Chong, A. Goujon, B. Gillon, J. Jeftic, A. Hauser, *Eur. Phys. Lett.* **2007**, *77*, 30007.

[9] D. A. Pejaković, J. L. Maison, C. Kitamura, J. S. Miller, A. J. Epstein, *Polyhedron* **2001**, *20*, 1435.

[10] K. Kato, Y. Moritomo, M. Takata, M. Sakata, M. Umekawa, N. Hamada, S. Ohkoshi, H. Tokoro, K. Hashimoto, *Phys. Rev. Lett.* **2003**, *91*, 255502.

[11] H. Banerjee, S. Chakraborty, T. Saha-Dasgupta, *Inorganics* **2017**, *5*, 47.

[12] A. Gîndulescu, A. Rotaru, J. Linares, M. Dimian, J. Nasser, *J. Phys: Conf. Ser.* **2011**, *268*, 012007.

[13] M. Nishino, S. Miyashita, P. A. Rikvold, *Phys. Rev. B*, **2017**, *96*, 144425.

[14] C. Enachescu, L. Stoleriu, A. Stancu, A. Hauser, *Phys. Rev. B*, **2010**, *82*, 104114.

[15] M. Sorai, S. Seki, *J. Phys. Chem. Solids* **1974**, *35*, 555.

[16] M. M. Dîrtu, C. Neuhausen, A. D. Naik, A. Rotaru, L. Spinu, Y. Garcia, *Inorg. Chem.* **2010**, *49*, 5723.

[17] W. Nicolazzi, J. Pavlik, S. Bedoui, G. Molnár, A. Bousseksou, *Eur. Phys. J. Special Topics* **2013**, *222*, 1137.

[18] M. Paez-Espejo, M. Sy, K. Boukheddaden, *J. Am. Chem. Soc.* **2016**, *138*, 3202.

[19] M. A. Halcrow, *Spin-Crossover Materials: Properties and Applications*, Wiley, New York **2013**.

[20] P. Gütlich, A. B. Gasper, Y. Garcia, *Beilstein J. Org. Chem.* **2013**, *9*, 342.

[21] C. M. Quintero, G. Félix, I. Suleimanov, J. S. Costa, G. Molnár, L. Salmon, W. Nicolazzi, A. Bousseksou, *Beilstein J. Nanotechnol.* **2014**, *5*, 2230.

[22] E. König, *Struct. Bond.* **1991**, *76*, 51.

[23] H. Spiering, N. Willenbacher, *J. Phys. Condens. Matter* **1989**, *1*, 10089.

[24] I. Gudyma, V. Ivashko, J. Linares, *J. Appl. Phys.* **2014**, *116*, 173509.

[25] M. Nishino, S. Miyashita, K. Boukheddaden, *J. Chem. Phys.* **2003**, *118*, 4594.

[26] J. Wajnflasz, *Phys. Status Solidi B* **1970**, *40*, 537.

[27] J. A. Nasser, *Eur. Phys. J. B* **2001**, *21*, 3.

[28] J. A. Nasser, S. Topçu, L. Chassagne, M. Wakim, B. Bennali, J. Linares, Y. Alayli, *Eur. Phys. J. B* **2011**, *83*, 115.

[29] A. Rotaru, A. Carmona, F. Combaud, J. Linares, A. Stancu, J. Nasser, *Polyhedron* **2009**, *28*, 1684.

[30] A. Rotaru, J. Linares, S. Mordelet, A. Stancu, J. Nasser, *J. Appl. Phys.* **2009**, *106*, 043507.

[31] J. A. Nasser, *Eur. Phys. J. B* **2009**, *48*, 19.

[32] A. Rotaru, J. Linares, E. Codjovi, J. Nasser, A. Stancu, *J. Appl. Phys.* **2008**, *103*, 07B908.

[33] A. Slimani, K. Boukheddaden, K. Yamashita, *Phys. Rev. B* **2015**, *92*, 014111.

[34] M. Nishino, S. Miyashita, *Phys. Rev. B* **2013**, *88*, 014108.

[35] K. Boukheddaden, M. Nishino, S. Miyashita, *Phys. Rev. B* **2007**, *75*, 094112.
[36] J. A. Nasser, K. Boukheddaden, J. Linares, Eur. Phys. J B **2004**, *39*, 219.

[37] K. Boukheddaden, S. Miyashita, M. Nishino, Phys. Rev. Lett. **2008**, *100*, 177206.

[38] S. A. Deviren, E. Albayrak, Physica A **2011**, *390*, 3283.

[39] M. Acharyya, Phys. Rev. E **1999**, *59*, 128.

[40] T. Tomé, M. J. de Oliveira, Phys. Rev. A **1990**, *41*, 5252.

[41] R. Erdem, M. Keskin, Phys. Rev. E **2001**, *64*, 026102.

[42] O. Canko, B. Deviren, M. Keskin, J. Phys. Condens. Matter **2006**, *18*, 6635.

[43] M. Keskin, O. Kanko, Ü. Temizer, Phys. Rev. E **2005**, *72*, 036125.

[44] S. W. Sides, P. A. Rikvold, M. A. Novotny, Phys. Rev. Lett. **1998**, *81*, 834.

[45] M. Keskin, O. Canko, J. Korean, Phys. Soc. **2009**, *55*, 1344.

[46] M. Keskin, M. Ertaş, O. Canko, Phys. Scr. **2009**, *79*, 025501.

[47] M. Bati, M. Ertaş, J. Supercond. Nov. Magn. **2016**, *29*, 2835.

[48] T. D. Oke, F. Hontinfinde, K. Boukheddaden, Eur. Phys. J. B **2013**, *86*, 271.

[49] T. D. Oke, F. Hontinfinde, K. Boukheddaden, Appl. Phys. A **2015**, *120*, 309.

[50] T. D. Oke, F. Hontinfinde, K. Boukheddaden, Comput. Conden. Matter **2016**, *9*, 27.

[51] a) M. Nishino, K. Boukheddaden, S. Miyashita, F. Varret, Polyhedron **2005**, *24*, 2852.; b) M. Nishino, K. Boukheddaden, S. Miyashita, F. Varret, Phys. Rev. B **2003**, *68*, 224402, and references therein.

[52] C. Enachescu, L. Stoleriu, A. Stancu, A. Hauser, Phys. Rev. Lett. **2009**, *102*, 257204.

[53] S. Mouri, K. Tanaka, S. Bonhommeau, N. O. Moussa, G. Molnár, A. Bousseksou, Phys. Rev. B **2008**, *78*, 174308.

[54] T. Sauer, Numerical Analysis, 2nd ed., Pearson Education, Inc., Hoboken, NJ **2012**.

[55] R. L. Burden, J. D. Faires, Numerical Analysis, 9th ed., Boston, Cengage, USA **2010**.

[56] A. Quarteroni, R. Sacco, F. Saleri, Méthodes Numériques: Algorithmes, Analyse et Applications, Springer-Verlag, Milano **2007**.

[57] S. D. Conte, C. de Boor, Elementary Numerical Analysis: An Algorithmic Approach, 3rd ed., McGraw-Hill, New York, USA **1980**.

[58] M. Nishino, K. Boukheddaden, S. Miyashita, F. Varret, Phys. Rev. **2005**, *72*, 064452.

[59] A. B. Bortz, M. H. Kalos, J. L. Lebowitz, J. Comput. Phys. **1975**, *17* L733.

[60] S. Jain, Monte Carlo Simulations of Disordered Systems, World Scientific Publishing Co. Pte. Ltd, Singapore **1992**.

[61] K. Binder, D. W. Heermann, Monte Carlo Simulation in Statistical Physics, 5th ed., Springer-Verlag, Berlin Heidelberg **2010**. https://doi.org/10.1007/978-3-642-03163-2

[62] A. Bousseksou, F. Varret, M. Goiran, K. Boukheddaden, J.-P. Tuchagues, Top. Curr. Chem. **2004**, *235*, 65.

[63] P. Gütlich, H. A. Goodwin, Spin-Crossover Transition Metal Compounds III, Springer-Verlag, Heidelberg, Berlin **2004**.

[64] A. Bousseksou, N. Nègre, M. Goiran, L. Salmon, J. P. Tuchagues, M.-L. Boillot, K. Boukheddaden, F. Varret, Eur. Phys. J. B **2000**, *13*, 451.

[65] S. Bonhommeau, G. Molnár, M. Goiran, K. Boukheddaden, A. Bousseksou, Phys. Rev. B **2006**, *74*, 064424.

[66] a) S. Kimura, T. Otani, Y. Narumi, K. Kindo, M. Nakano, G. Matsubayashi, J. Phys. Soc. Jpn. **2003**, 75SB, 122; b) P. G. Sim, E. Sinn, J. Am. Chem. Soc. **1981**, *103*, 241.

[67] Y. Garcia, O. Kahn, J.-P. Ader, A. Buzdin, Y. Meurdesoif, M. Guillot, Phys. Rev. Lett. A **2000**, *271*, 145.