Hysteresis dynamics of rare earth nickelates: unusual scaling exponent and asymmetric spinodal decomposition

Understanding the dynamics of phase-transitions, interpretations of their experimental observations and their agreement with theoretical predictions continue to be a long-standing research interest. Here, we present detailed phase-transition dynamics of rare earth nickelates associated with its first-order metal–insulator transition. The thermal hysteresis shows absence of training effect and defies the Preisach model. A large phase-coexistence in insulating state during cooling suggests kinetically arrested glassy dynamics of the phase-transition. Experimentally derived hysteresis scaling exponent is much larger than the mean-field predicted universal value of 2/3. In the phase-coexistence region, the quench and hold measurement depicts higher stability of the metallic state compare to that of the insulating one; highlighting the manifestation of phase-coexistence via asymmetric spinodal decomposition. All these observations for nickelates are in stark contrast to the phase-transition dynamics of canonically similar vanadates but are closer to those of glasses, alloys. A substantial disagreement between the experiment and theory emphasizes the necessity to incorporate system-dependent details for the accurate interpretation of the experimental results.


Introduction
Under constraints of many-body interactions, the response of a system to a time varying external perturbation is not instantaneous due to relaxational delay. This makes the response hysteretic in nature; a non-equilibrium and nonlinear phenomenon. A variety of physical processes such as ferroic orderings, spin crossover, first-order phase-transitions, etc, occur with prominent hysteresis [1][2][3][4][5]. Knowledge of the dynamics associated with the hysteresis can help understand the underlying mechanisms of such phases as well as of the states which are far from equilibrium. Further, controlling the shape and width of the hysteresis has various technological aspects too; for example, in data storage, numerical display, memristive switching devices, sensors, modulators, etc, [6,7]. Given the fundamental interest in hysteretic phenomena and its abundance in nature, there has been immense research interest in past, to model the hysteretic processes.
Preisach presented a mathematical model for the hysteresis and highlighted its wiping out and congruency properties [8]. Ising model was developed originally to explain magnetic ordering and transitions which considers magnetic spins to interact only with its nearest neighbor spins [9]. This oversimplified model could explain qualitatively magnetic transition in two-dimensional systems, however, it is numerically tough to solve Ising model for three dimension which is the case for real systems. To make the calculation easy, mean-field approximation came into picture which replaces all the nearest neighbor interactions with an effective field and neglects effect of fluctuations about its mean value [10]. One important feature of the hysteretic process was revealed that the change in the area of the hysteresis loop follows scaling relation and the derived scaling exponent is universal, independent of the material systems [11,12]. After rigorous works of Zhong and co-workers, applicability and universal nature of hysteresis scaling were also suggested to be equally valid in the evolution of linearly driven first-order hysteretic phase-transitions [13][14][15][16]. Similarly, a number of experimental methods and tools were developed to extract the microscopic details of hysteretic processes and validate it with the theoretical predictions, for example, first-order reversal curve (FORC) and quench and hold measurements which tell about coexistence region at microscopic level which may not be observable from experimental hysteresis loop and stability of metastable state, respectively [5,[17][18][19].
The success of these existing theoretical models and experimental methods in the study of hysteresis dynamics of many material systems suggests that these models and experimental methods can have great potential to unveil the microscopic mechanisms of a variety of features exhibited by complex systems, particularly 3d-transition metal oxides (TMOs). Due to subtle competition among lattice, charge, spin and orbital degrees of freedom, in general, the electronic states of TMOs are fragile. As a result, multiple phases co-exist in degenerate conformations [20,21]. A slight change in external perturbation drastically modifies the strength of different phases and their spatial distribution in the phase-coexistence region. Hence, it is important to understand the dynamics of phase-coexistence around the phase-transition. The studies carried out on thermal hysteresis associated with first-order metal-insulator transition (MIT) in vanadates indeed show a good agreement between the theoretical predictions and the experimental results. For example, experimental hysteresis was reproduced by incorporating mean-field approximation and effective medium theory [5]. The hysteresis scaling exponent deduced from experimental data also lies on the close limit of theoretically predicted universal value of 2/3 [18]. However, the universal nature of the phase-transition and the hysteresis dynamics of the entire range of TMOs are not explored yet.
Rare earth nickelates (RNiO 3 ) are one of the contemporary systems which exhibit first-order MIT with prominent hysteresis. The recent studies on nickelates highlight its several peculiar behaviors, to name a few: (i) nickelates (a charge-transfer insulator), similar to vanadates (a Mott-Hubbard insulator), have been found to exhibit carrier mass divergence type of MIT which was initially believed to be carrier number vanishing type [22]. (ii) Nickelates show the evidence of continuous domain wall formation (a second-order phase-transition) along with abrupt transition of insulating clusters to metallic clusters (a first-order transition) at the onset of MIT [23]. This is a rare example of coexistence of first and second-order phase-transitions and so far, only has been found in nickelates. (iii) A deeper investigation on the metallic nature of nickelates reveals that the temperature exponent of resistivity is extremely sensitive to epitaxial strain and can be tuned from n = 1 to 3 [24]. Thus, nickelates are actually not a bad metals (which is most of the time loosely spoken) but their metallic states can be tuned from Planckian type to Fermi liquid to bad metal type [25][26][27][28] and thus multiple mechanisms simultaneously play roles to determine its metallic nature. Nevertheless, nickelates, within the strongly correlated 3d-TMOs, are closest to vanadates as they share several similarities: both are Mott insulators and exhibit similar phase diagrams incorporating structural, magnetic and electronic transitions [29]. So, it is natural to expect, like vanadates, the experimental hysteresis dynamics of nickelates to be in accordance with theoretical predictions. However, we do not find this true as our detailed experimental investigation shows substantial deviation of experimental results from the theoretical predictions thus unravelling one more peculiar behavior of nickelates.
We prepared two nickelate films: NdNiO 3 (NNO) film fabricated on single crystal LaAlO 3 (LAO) (100) substrate and La 0.5 Eu 0.5 NiO 3 (LENO) film on NdGaO 3 (NGO) (001) substrate which has nearly same lattice constant to that of NNO (∼3.807Å) but has quenched disorder due to different R-site cation sizes of La 3+ (1.216Å) and Eu 3+ (1.120Å) [30]. We observed qualitatively similar hysteresis dynamics in both the films; a generic nature of the hysteresis dynamics of nickelate system. Both the films show no training effect, violation of Preisach's congruency theorem and the phase-transition is dominated by kinetically arrested like glassy dynamics during cooling cycle. The obtained hysteresis scaling exponent is much larger than the mean-field predicted universal value of 2/3. The phase-coexistence manifests via asymmetric spinodal decomposition. These properties of nickelates are closer to structurally complex glasses, alloys rather than vanadates. The deviations of experimental results from theoretical predictions emphasize the importance of incorporating system-specific details in the development of a theoretical model.

Experimental details
The epitaxial nickelate thin films were fabricated using pulsed laser deposition technique. The deposition parameters used for film growth were: 710 • C deposition temperature, 20 Pa deposition oxygen pressure and 2.5 mJ cm −2 (in case of NNO film)/3.0 mJ cm −2 (in case of LENO film) laser pulse fluence with 8 Hz pulse repetition rate. After deposition, the films were annealed at the deposition temperature in 200 Pa of oxygen pressure for five minutes. Finally, the films were cooled down to room temperature at the rate of 5 • C min −1 . Structural characterizations were performed using PANalytical x-ray diffractometer. Electronic transport properties were investigated in the temperature range of 2-300 K with the help of Physical Property Measurement System in four probe arrangement. See supplementary material, sections S1 and S2 for the details of structural analysis and transport measurements.

Training effect and Preisach model
The temperature (T ) dependent DC resistivity (ρ) of the NNO film, as plotted in figure 1(a), shows first-order MIT with a pronounced thermal hysteresis. To explore the possibility of training effect, we repeated the measurement of ρ-τ loop with a fixed temperature ramp rate of 3 K min −1 ( figure 1(a)). Overlapping of all the loops within the experimental resolution rules out the possible training effect. About five orders of change in resistivity across MIT and absence of training effect confirm the film of high quality. Figure 1(b) shows construction of minor hysteresis loops inside major hysteresis loop in the temperature range of 120-130 K (chosen arbitrarily) as follows: the film was cooled from 300 K to 115 K then warmed back to 130 K. From 130 K, a closed loop was formed by traversing the temperature 130 K → 120 K → 130 K. Now, the film was cooled from 130 K to 110 K and then the process was repeated to construct another closed loop in the range of 120-130 K. Subsequently cooling the film up to 105 K and 100 K, two more minor loops were constructed. We chose subsequent cooling temperatures with 5 K difference so that there was a sufficient vertical shift between any two consecutive minor loops (see figure 1(b)). The wiping out and congruency properties are necessary and sufficient conditions for a hysteresis loop to follow Preisach model [8]. The former makes sure that the major hysteresis loop is well-defined irrespective of its past history while the latter ensures that the shape and area of all the minor hysteresis loops are same [4]. Here, the NNO film exhibits well-defined thermal hysteresis and lacks training effect, thus, satisfies the wiping out property. While, a significant difference in the areas of minor hysteresis loops clearly shows the failure of congruency property (see figure 1(b)). This suggests that the thermal hysteresis of the NNO film does not follow the Preisach model. The reason of the non-applicability of Preisach model lies in its inherent assumption that the major hysteresis is an ensemble of a number of small hysteresis, called hysterons, which do not interact with each other [8]. However, in presence of long range interaction in case of nickelates (strongly correlated system), such hysterons certainly feel the presence of each other, thus, resulting in the failure of Preisach model in nickelate films.

FORC analysis
Unlike the Preisach model which is entirely a mathematical construction and is based on some hypothetical assumptions, the FORC distribution is of completely analytical in nature and is applicable to any type of hysteretic process [5]. A FORC data is collected as follow: starting from a high temperature metallic state (here 300 K) the sample is cooled down to a lower temperature, called reversal temperature (T R ) and then warmed up again to 300 K. The ρ-τ curve, starting from T R to 300 K is defined as FORC (figure 2(a)). By cooling down to discrete T R values in the range of 2-300 K and subsequently warming up to 300 K, a large number of FORCs were collected. These data are called 'heating FORCs' as the temperature is always raised from each T R to higher temperatures. In similar pattern, cooling FORCs were also collected by varying T R and cooling the sample to 2 K from each T R ( figure 2(b)). Finally, heating and cooling FORC distributions were calculated by taking the mixed second order derivative of the reversal curves: where r(T R , T) = log 10 [ρ(T R , T)] and plotted in figures 2(c) and (d), respectively for NNO film. The second order derivative eliminates those resistivity data-points which remain constant on varying T or T R . Thus, the non-zero value of Ω(T R , T) represents irreversible parts which are introduced by coexistence of insulating/metallic clusters during the phase-transition. Larger magnitude of Ω(T R , T) corresponds to larger irreversibility or stronger phase-coexistence. In both the heating and cooling FORCs, the strength of phase-coexistence, as expected, is maximum around the T MI and along the T = T R line (figures 2(c) and (d)). Further, in both the cases, Ω(T R , T) has some non-zero values in low temperature insulating region. This suggests that the system is not in a single electronic state even at much lower temperatures far away from the hysteretic region. A small fraction of metallic clusters exists in the insulating region. Most importantly, the nickelate film shows a clear difference between heating and cooling FORCs. The amount of low temperature phase-coexistence in insulating region is very large in cooling FORC compared to that in heating FORC (see figures 2(c) and (d)). This signifies that the phase-transition during cooling protocol remains largely incomplete which is a well-known signature of kinetically arrested glassy dynamics [31].

Hysteresis scaling
To investigate the universal nature of the MIT dynamics of the nickelate film, we performed the linear temperature ramp rate dependent resistivity measurements ( figure 3). In general, the delay in phase-transition (T MI ) (or equivalently hysteresis loop area) increases with increasing linear temperature ramp rate (R) and is predicted to follow a universal scaling form, given by: where T MI (R 0 ) is the free parameter corresponding to the transition temperature in quasistatic condition, 'a' is a constant and 'α' is the scaling exponent [13]. '+' sign corresponds to heating protocol and '−' sign corresponds to cooling protocol. Under the mean-field model, the value of 'α' is predicted to be 2/3 in both the heating and cooling protocols, a universal value independent of material system [18,32].  shows the evolution of hysteresis loop of the NNO film as a function of temperature ramp rate while insets show the magnified view of the shift in hysteresis curve in the direction of arrows on increasing the ramp rate for both the protocols. We extracted the hysteresis scaling exponent by scaling the T MI with respect to ramp rate using equation (2) (inset in figure 3). The extracted values of scaling exponent (α) are 0.94 ± 0.07 and 0.98 ± 0.04 in heating and cooling protocols, respectively. This is clear deviation from the mean-field predicted universal scaling exponent of 2/3. Note that such large values have been reported so far only in case of glasses and alloys [14,33,34], nickelate is the first known oxide system which exhibits such a large scaling exponent.

Quench and hold measurement
The stability of metastable insulating/metallic clusters in the phase-coexistence state during phase-transition can give more concrete evidences for the observed kinetically arrested glassy dynamics and large hysteresis scaling exponent. The stability of a metastable state can be estimated by quench and hold experiment. When a system is quenched and then is held, it relaxes from initial metastable state to a stable state. Hence, the time-evolution of metastable state into stable state provides an estimation of the stability of such metastable state. We rapidly cooled (heated) the NNO film slightly below (above) the T MI up to different temperatures and observed the time-dependent change in resistivity (Δρ) of the film at those temperatures (figures 4(a) and (b)). In both the heating and cooling quenches, at any temperature, Δρ has two parts: sharp fall/rise for a small time-interval, say, for 3 min followed by a slow exponential fall/rise at later times. The sharp fall/rise in Δρ is known to occur due to nonlinear effect of phase-transition, lattice distortion and sudden growth of metallic/insulating clusters [19]. This also suggests that a major fraction of metallic (insulating) clusters converts into insulating (metallic) clusters as the sample temperature crosses the T MI upon cooling (heating). The slow exponential fall/rise indicates the slowing down of the phase-transition as metastable state converts into stable state. It is also seen that Δρ does not saturate even after much longer time. The non-saturation is more prominent in cooling quench than in heating quench. This clearly shows that the time-scale of metal to insulator transition in cooling is much larger than that of insulator to metal transition in heating sequence, indicating occurrence of kinetic arrest of phase-transition upon cooling. For quantitative estimation, we fitted the slow part of Δρ to an exponential decay function: where ρ 0 and A are constants and τ is relaxation time for the respective cluster growths (figure 4) [19]. The upper signs in right-hand side of equation (3) correspond to heating protocol while bottom signs are for cooling protocol. In both the cases, τ peaks around the T MI (insets in figures 4(a) and (b)). This behavior of τ is in accordance with the fact that the stability of metastable state is highest around T MI and decreases away from the T MI . Further, the value of τ is 4-5 times higher in cooling protocol than in heating protocol. This suggests that much of the metallic clusters do not transition into the insulating clusters but remain in metastable metallic state during cooling. Thus, the phase-transition exhibits glassy dynamics. It also reveals that the stability of metastable metallic phase in phase-coexistence region is much higher than that of metastable insulating phase. In general, stability of different metastable phases in phase-coexistence region is nearly the same and this coexistence is defined to occur through 'symmetric spinodal decomposition'. On the other hand, when phase-coexistence region consists of metastable phases of different stability limit then such coexistence is defined to occur through 'asymmetric spinodal decomposition' [35]. Clearly, the phase-coexistence associated with nickelates occurs through asymmetric spinodal decomposition. As the hysteresis dynamics of a system is significantly affected by the presence of disorder in the system, we carried out preceding experimental investigations on the quenched disorder LENO film [30], to unveil the role of disorder in our case. In general, presence of a disorder hinders the onset of any phase ordering and increases electronic inhomogeneity. In our case, it is expected that the quenched disorder would reduce the difference between the observed dynamics of phase-transitions in heating (a clean insulator to metal transition) and cooling (a kinetically arrested metal to insulator transition). Hence, one can expect that the phase-transition dynamics data might match with the theoretical predictions. On the contrary, despite of having much reduced sharpness of the transition and strength of irreversible part in FORC, the hysteresis dynamics of LENO film is qualitatively similar to that of NNO film (see supplementary material sections S3-S5). It also shows clear difference between heating and cooling FORCs and that the hysteresis scaling exponent is greater than 2/3 (supplementary material section S5). The quench and hold measurement (figures 4(c) and (d)) shows no time-dependence of Δρ in heating protocol while a weak time-dependence (one order lesser compare to that of NNO film) exists in cooling protocol. Overall, it is surmised that the phase-coexistence manifests via asymmetric spinodal decomposition for disordered LENO film as well. Due to disorder induced electronic inhomogeneity, Δρ in cooling protocol for the LENO film does not have any specific temperature dependence as reflected in its relaxation time pattern (inset of figure 4(d)). We did not try to fit Δρ and extract relaxation time in heating protocol due to its negligible time-dependence in this case. These observations of LENO film confirm the kinetically arrested glassy dynamics of phase-transition during cooling as well as a generic nature of hysteresis dynamics of entire nickelate system.

Theoretical consideration: compressible Ising model
Finally, we present theoretically simulated resistivity curve obtained after using compressible Ising model and compare it with experimental resistivity curve to see how much the experimental hysteresis of nickelate Figure 5. Variation of order parameter (insulator fraction, φ) with temperature: a comparison between transport measurement and theoretical result generated from compressible Ising model. The quasi-static order parameter (derived from quench and hold data) is also shown with navy blue circles (cooling) and orange circles (heating). matches with or deviates from the theoretically reproduced one. The compressible Ising model, introduced by Domb [36] and Henriques [9], is an extension of Ising model where mean field Hamiltonian is modified by lattice distortion. The effect of lattice distortion on the Hamiltonian of the system is included by an additional harmonic elastic energy ∝ N 2 (v − v 0 ) 2 ; N is number of sites, v 0 is equilibrium volume and v is volume after lattice distortion. It is also assumed that the lattice distortion affects the spin-spin exchange interaction in the Ising model given by; ; J 0 is exchange interaction constant in equilibrium and J d is the change in exchange interaction strength due to lattice distortion. With these inclusions, this model provides description of structural phase-transition and first order temperature driven phase-transition. The final expression for the free energy per site (f ), in this case, is given by [18]: where ξ = j 2 d q 8j 0 K and T c = qj 0 K B ; q is the coordination number of lattice, K is the positive parameter related to the inverse of compressibility and Φ is the order parameter (here insulator volume fraction). See supplementary section S6 for detailed derivation of equation (4). The evolution of order parameter Φ in mean field picture can be determined by dissipative time-dependent Landau equation or model A [37,38]: One of the two free parameters of the model, T C is already specified by the cooling spinodal temperature, i.e., T C ≈ 128 K and the other free parameter ξ is numerically calculated to be 0.1696 from the value of the heating spinodal temperature 139 K. The sample is considered as an inhomogeneous ensemble with a Gaussian distribution of mean 128 K with standard deviation of 9 K. The value of the parameter λ is 3.5 s −1 .
In figure 5, we show the temperature evolution of order parameter Φ calculated from the resistivity measured with temperature ramp rate of 0.5 K min −1 for the NNO film (as described in supplementary section S7). We also calculated quasi-static order parameter from quench and hold measurement. The quasi-static order parameter corresponds to value of resistivity after a long wait (∼45 min) once the sample is quenched and then held at the given temperature (as shown in figure 4). The orange and navy blue circles in figure 5 show temperature evolution of thus obtained quasi-static order parameter at a few temperatures. It is evident that the quasi-static order parameter differs marginally from the dynamic order parameter (i.e., order parameter evolution at temperature ramp rate of 0.5 K min −1 ). That indicates the athermal character of the transition i.e., metastable phase of the system is arrested if the temperature is kept fixed. In spinodal decomposition formalism, system can sustain in its metastable phase up to spinodal (when the activation barrier vanishes). Further, in figure 5, the dotted purple line shows the simulated curve obtained after using equations (4) and (5) with choosing the free parameter properly. We clearly see that the simulated curve significantly deviates from the experimental curve. Unlike the case for vanadate system, the theoretical hysteresis does not match properly with the experimental one for nickelate system. This highlights the crucial role played by disorder and thermal fluctuations in hysteresis dynamics which equations (4) and (5) ignore. Depending on the system, the source of disorder and thermal fluctuations can be different. In the current case, electronic complexity, kinetic arrest and quenched disorder are the main source for disorder and thermal fluctuations. Also, equation (5) assumes that each metallic/insulating domain has its own transition temperature and the individual transition temperature follows Gaussian distribution centered around the mean transition temperature. However, in real scenario, the distribution can deviate from Gaussian behavior.

Discussion
Now, we point out stark contrast in heuristically observed hysteresis dynamics of nickelate with that of vanadate and its similarity with glass or alloy systems which is the main highlight of this work. The experimentally observed hysteresis dynamics and phase-transition in vanadates closely follow the theoretical predictions and universal scaling; no difference between the phase-transitions in heating and cooling is observed [5,18]. The compressible Ising model is able to reproduce experimental hysteresis [18]. Even in its pressure-temperature phase diagram in the vicinity of critical endpoint, the critical exponents were found to be consistent with those were extracted using mean-field theory [39]. However, as described in preceding sections, dynamics of phase-transition in nickelates in heating and cooling protocols clearly differ from each other, hysteresis scaling exponent deviates from the mean-field predicted value, exhibiting large value of scaling exponent. Further, the time evolution of metastable state into stable state in quench and hold measurement is much slower in case of nickelate (order of 10s of minutes) compared to in case of vanadate (order of 10s of seconds). Finally, the compressible Ising model is not able to reproduce completely the experimental hysteresis of nickelate. Interestingly, all these features of hysteresis dynamics of nickelate closely resembles with that in case of glass or alloy systems. Thus, our study highlights unusual phase-transition dynamics of rare earth nickelates.
Kinetically arrested like glassy dynamics of phase-transition is a normal feature for alloys and glasses, e.g. FM-AFM transition in doped CeF 2 alloys, magneto-structural transition in Gd 5 (Ge x Si 1−x ) 4 , FeRh alloy, magnetic shape memory alloys, geometrically frustrated spin chain compounds (e.g. Ca 3 Co 2 O 6 ), etc, [40][41][42][43][44]. Kinetic arrest in such systems is a result of structural complexity. While the glassy dynamics in relatively structurally simple compounds such as manganites, cobaltites, cuprates, inter-metallics and multiferroics is argued to occur due to electronic complexities [31]. At low temperatures, orderings of more than one electronic phase and phase-coexistence make homogeneous phase unstable, leading to competing ground states and kinetically arrested glassy dynamics of phase-transition during cooling.
The occurrence of kinetic arrest of phase-transition in nickelates while cooling increases the stability of metastable metallic clusters and causes the phase-transition to occur through asymmetric spinodal decomposition. For the same reason, the hysteresis scaling exponent also deviates significantly from the theoretically predicted one. Such deviations in scaling exponent have been reported in some other systems too (see supplementary table S1) which put validity of universal nature of hysteresis scaling exponent in question. Kagawa et al [45] observed unusual critical behavior in quasi two-dimensional organic superconductor (a system quite similar to cuprate superconductors) and suggested the need of new universality class. In their case, though the experimentally deduced critical exponents follow the scaling relation, their values significantly deviate from the theoretical ones. Note that the hysteresis scaling exponent indirectly tells that how much dissipative a system is [14]. Higher the dissipative a system is, larger will be the energy dissipation in driving the system and so the larger the hysteresis scaling exponent. The large energy dissipation in driving a hysteretic process may occur due to various reasons depending on the system which mean-field approximation does not consider. For example, structural complexity in glass and alloy systems, magnetic anisotropy, domain wall movements and thermal fluctuation (random noise) etc, in a magnetic hysteresis [46] and electronic inhomogeneity in complex systems where fundamental degrees of freedom are strongly coupled and give rise to coexistence of multiple phases and kinetically arrested phase-transition. The hysteresis scaling value in such system would be larger than the mean-field predicted universal value of 2/3. Our study invokes the necessity to include such details in the existing models to appropriately predict the experimental results. Such rigorous experimental studies of hysteresis dynamics on a variety of systems have not been performed sufficiently yet. In this scenario, the present study was much desired to guide the theoretical approaches and particularly, to validate the universal nature of the scaling exponents.

Conclusions
To sum up, our study reveals that the hysteresis mechanism of nickelates is dominated by kinetically arrested glassy dynamics and the phase-transition occurs through asymmetric spinodal decomposition. These findings put nickelate system closer to structurally complex glasses and alloys rather than vanadate system. The deviation of experimental hysteresis scaling exponent from theoretical prediction creates a huge scope for refinement of existing theoretical framework by incorporating the system-dependent details. Such elaborated experimental investigation of hysteretic behavior of a variety of systems can serve as a testing ground for further theoretical improvements.