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

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.


Introduction
Because of the needs for our companies in the field of the treatment and storage of informations, transportation of artworks 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 layout leads interactions between the component molecules of materials, which cause the appearance of unique or multiple physical properties (conduction, magnetism, photochromism, nonlinear optics . . . ). [2][3][4][5][6] Furthermore, their molecular aspect makes them easily flexible chemically to optimize their properties or to combine them with other molecular systems to generate multifunctional materials. Therefore, under the influence of various stimuli such as light, pressure, temperature, magnetic and electric fields, etc., [7][8][9][10] phase transitions between the low-spin (LS) diamagnetic state and the high-spin (HS) paramagnetic state are triggered. [7,9,10] The change of induced spin state allows to obtain devices for high density information storage in which the 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][14][15] which manifest through the existence of first-order phase transitions accompanied with thermal hysteresis. An interesting example is that of the [F e(N H 2 tr z) 3 ](NO 3 ) 2 , (N H 2 tr z = 4 − amino − 1, 2, 4 − tr iazol e) 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 F e(I I), the SCO phenomenon is the result of the redistribution of the electrons between the bonding t 2g and the antibonding e g orbitals. In the diamagnetic (σ = 0) LS state, only bonding orbitals are populated (t 6 2g e 0 g ), while in the paramagnetic (σ = 2) HS state, the electronic configuration becomes (t 4 2g e 2 g ). As demonstrated in several works, [13,14,18] the elastic interactions are at the hearth of the existence of cooperative effects in SCO materials and play a crucial role in the existence of the www.advancedsciencenews.com www.advtheorysimul.com 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][19][20][21][22][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][25][26] atom-phonon coupling descriptions, [12,[27][28][29][30][31][32] or more recently on complex elastic descriptions [13,23,[33][34][35][36] or taking into account the volume change at the transition. These models are solved analytically using meanfield 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][39][40][41][42][43][44] or mixed spins kinetic models, [45][46][47] which exhibited interesting physical properties. In this work, we consider the BEG model for SCO or PBAs [24][25][26][27][28][29][30][31][32][33][34][35][36] and crossover solids under oscillating magnetic field to describe their thermal and time-dependent properties. As in our previous works, [48][49][50] the quadrupolar coupling parameter K between SCO units is assumed to depend linearly on the absolute temperature T in the form, K = α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 = σ (magnetization) and n HS = σ 2 (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][50][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 spatiotemporal 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.

Hamiltonian Model
The Hamiltonian of spin-1 BEG model adapted for SCO materials with magnetic interactions and applied magnetic field is written as follows: where σ i = ±1, 0 are fictitious spin values located at site i of the square lattice. The spins σ i = ±1 describe the magnetic HS spin state and σ 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, σ i = ±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][34][35][36][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 longrange interactions between the SCO units, and is written here as K = αk B T [48][49][50] with the ratio γ = J /K taken as a tunable parameter. D = − 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(ωt) is considered where ω is the oscillation frequency.

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 where m = σ i and n HS = σ 2 i are considered as invariant by translation over the lattice.
The associated mean-field free energy per site is given by where U and S are the internal energy and entropy per site of the system, respectively, given by Here, β = 1 k B T (T is the temperature) and H i = −zJ m 2 − zK n 2 HS + Dn HS − hm is the average value of H i . After some simple calculations, the variational free energy is given by

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({σ i }, t) to occupy the spin configuration {σ } [49][50][51]53] at time t. The flux of probability accounts for transitions from the configuration {σ } i {σ } i with transition rates W. Following the Glauber-type stochastic dynamics, the master equation reads where τ stands for the Arrhenius time scale and 1 τ denotes the effective intramolecular frequency associated with the "spontaneous spin reversals" [49][50][51] where 1/τ 0 is the "intrinsic" frequency spin-flip process between HS and LS states, taken as constant and E a 0 denotes the intramolecular vibronic energy barrier.
In the equilibrium state, the probabilities satisfy the detailed balance condition Thus, the probability per unit time is given by Thus, these transition rates W, must fulfill the detailed balance condition.

The Motion Equations
The mean value of spin σ k and σ 2 k , respectively, associated to the magnetization and HS fraction calculations are and σ k =σ k W k (σ k ) = 1 3τ . After some calculations, one gets the time evolution of the magnetization m and the fraction n HS in the forms These coupled nonlinear differential equations (13) are numerically solved by using fourth order Adams-Moulton predictorcorrector methods. [54][55][56][57] Thermal properties are got by integrating on [0, 2π ] 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 where ξ = ωt. Here, ω becomes a simple parameter and we considered ξ as the control parameter in time. The method used to solve the system equations (14) is that of Romberg [54][55][56][57] with Richardson extrapolation to obtain thermodynamic quantities such as: magnetization m, fraction n HS , and the magnetic susceptibility given by Within the system dynamics, Van Hove-like equations [58] give the same system equations through the free energy per site. One gets

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 {σ } (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 {σ } is calculated. Let where [48] denotes the change in the system energy associated to the spin-flip move and 1 τ is defined as in where a stands for a given possible spin-flip process. Two random numbers 0 ≤ r 1 , r 2 ≤ 1 are chosen to calculate the lifetime is fulfilled is realized with probability 1. The total evolution timē t is given byt

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 (γ = 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 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 α. When α ≥ 1, first order transitions take place (see Figure 1b). Then, the reduced equilibrium temperature is a decreasing function of α. Up to α = 1.05, the analytical expression of T eq / [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 α-dependence of the free energy per site for selected values of . 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. .05 = 0 [49] is the onset of first-order transitions for selected values of (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.

Temperature-Dependent System and Phase Diagram
The effects of external magnetic field constraint on spincrossover compounds have been thoroughly studied in the  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][65][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 0 = 150 K is considered in the present work. The individual spin-flip frequency is set to the value 1 3τ 0 = 1s −1 and the activation energy to E a 0 = 0 K. The relaxation of the magnetization m and the HS fraction n HS 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 HS for some selected parameter values (Figure 4). First-order spin transitions characterized by jumps in m and n HS , 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 I I I [( pyr ol) 3 tr en] in the presence of a low magnetic field of the order of 23 T. [67] The ( pyr ol) 3 tr en is the trianionic Schiff base obtained from the condensation of pyr olle − 2 − car boxaldehyde with 2, 2 , 2 tr is (ethylamino)amine. The reduced equilibrium temperature T eq / 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 HS 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 (γ < γ T ), two temperatures may exist: T eq and T c with T eq < T c . 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 HS = 1. At the vicinity of the first-order transitions, we displayed in  backward branch (Figures 7a,b). When γ = 2 and α = 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).

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 × L = 100 × 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 γ = 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 α. 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) www.advancedsciencenews.com www.advtheorysimul.com  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(ω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 = σ 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 α and γ (Figure 9c,d).
At γ = 0 and from α = 3 to α = 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 spinflip process between HS/LS states is 1 3τ 0 = 10 2 s −1 and the intramolecular vibronic energy barrier E a 0 = 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   standing 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 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 σ = −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.

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 spintransition 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].