Giant thermal hysteresis in Verwey transition of single domain Fe3O4 nanoparticles

Most interesting phenomena of condensed matter physics originate from interactions among different degrees of freedom, making it a very intriguing yet challenging question how certain ground states emerge from only a limited number of atoms in assembly. This is especially the case for strongly correlated electron systems with overwhelming complexity. The Verwey transition of Fe3O4 is a classic example of this category, of which the origin is still elusive 80 years after the first report. Here we report, for the first time, that the Verwey transition of Fe3O4 nanoparticles exhibits size-dependent thermal hysteresis in magnetization, 57Fe NMR, and XRD measurements. The hysteresis width passes a maximum of 11 K when the size is 120 nm while dropping to only 1 K for the bulk sample. This behavior is very similar to that of magnetic coercivity and the critical sizes of the hysteresis and the magnetic single domain are identical. We interpret it as a manifestation of charge ordering and spin ordering correlation in a single domain. This work paves a new way of undertaking researches in the vibrant field of strongly correlated electron physics combined with nanoscience.

Strongly correlated electron systems are a group of materials, where electron correlations determine a delicate balance among the competing phases of the systems [1][2][3] . Naturally, the challenge becomes ever so mounting when one tries to handle a case where there are correlations involving not just one but all four degrees of freedom of solid: charge, orbital, lattice, and spin. This is exactly the case for the Verwey transition of Fe 3 O 4 , where the electronic degree of freedom (charge) is tightly coupled to all three others [4][5][6][7][8] . As the temperature becomes lower than the Verwey transition temperature, T V , a number of events take place simultaneously: charge ordering of Fe 3+ and Fe 2+ cations, the lattice symmetry change from cubic (Fd 3 m) to monoclinic (C c ), metal-insulator transition, sudden drop in magnetic susceptibility, and the magnetic easy axis change from (111) to (100). Resonant X-ray scattering technique also showed a hint of orbital ordering 9 . Due to such complexity of these correlated phenomena, the origin of the Verwey transition still remains unresolved even after about 80 year since its first report 4 .
The central question of strongly correlated electron systems is how the strong electron correlations emerge out of the assembly of individual uncorrelated atoms. To be specific to the Verwey transition of Fe 3 O 4 , it will be a very interesting study to investigate the size dependence of the transition. Given the nature of the complex interactions, it may well offer a new window of opportunities looking into the intimate details of the Verwey transition. The recent advance in nanoscience has brought new methodology to solid state physics for studying the emergent properties of a limited number of atoms in assembly. Through solution-based synthesis, colloidal nanoparticles (NPs) can be prepared with a precise size control in the range from sub-nanometer to tens of nanometers, which amounts to the number of atoms ranging from 10 1 to 10 20 . In terms of electronic energy structure (molecular orbital to band) and symmetry of atomic arrangement (point group to space group), transition from molecule to solid takes place at around 10 2 -10 5 atoms 10 . Armed with the latest huge developments in the nanochemistry field, one can open a new, exciting opportunity of monitoring the emergence of electron correlations at the borderline between molecules and solids. For example, the critical pressure of the solid-solid phase transition for CdSe NPs significantly increases as their sizes are reduced to less than few nanometers 11 , whereas the melting temperatures of CdS and metal NPs show an opposite trend [12][13][14] . Nevertheless, studies are relatively rare on more complex transitions emerging from electron correlations for colloidal NPs.
Previously, we developed a fine synthetic method to prepare stoichiometric Fe 3 O 4 NPs with various sizes to observe that T V becomes size dependent for the NPs smaller than 40 nm 15 . Here, we show another, yet hitherto unreported feature of the Verwey transition. By probing all the four degrees of freedom using independent techniques with superconducting quantum interference device (SQUID) magnetometer, cryogenic X-ray diffraction (XRD), 57 Fe nuclear magnetic resonance (NMR), and resonant inelastic X-ray scattering (RIXS), we reveal that the Verwey transition of stoichiometric Fe 3 O 4 NPs exhibits large size-dependent thermal hysteresis and that the hysteresis width is correlated with the critical size of the magnetic single domain, a feature is not discovered yet.

Results
The Verwey transition is extremely sensitive to oxidation such that it is substantially suppressed when off-stoichiometry parameter δ, defined as Fe 3(1-δ) O 4 , is larger than only 1% 16 . Therefore, we conducted preparation and storage of our Fe 3 O 4 NPs with extreme care to exclude any possible cause of off-stoichiometry 15 Fig. S7 and Table S1 in SI. Obtained Fe 2+ /Fe 3+ ratio for 7 μm is 1/1.9994 and that for 42 nm is 1/1.9985, which corresponds to the off-stoichiometry parameter δ = 0.0002 for 7 μm and δ = 0.0005 for 42 nm. Based on these results, we confirmed that both samples are high-stoichiometric within the error of δ = ±0.0003.
The magnetization data of our NPs show a clear correlation between thermal hysteresis and the particle size in nanometer scale (Fig. 1). Bulk standard sample (7 μm) exhibits a characteristic drop in magnetization at the known value of T V with a thermal hysteresis of 1 K, consistent with those reported in literature 17 . For the NPs, on the other hand, as the size decreases from 390 to 120 nm, the width of the thermal hysteresis increases from 2 to 11 K before falling down until the Verwey transition itself disappears below 10 nm 15 . When the heating/cooling rate varies from 0.1 to 10 K min −1 , the value of T V as well as the hysteresis width do not show any significant change within error bars (see Fig. S3 in SI). In controlled oxidation experiment done on 7 μm and 16 nm samples, the thermal hysteresis width of 16 nm sample after oxidation significantly change from 4.6 K to 0.7 K and that of 7 μm sample decrease from 1 K to 0.3 K (see Fig. S4 in SI). Therefore, we can rule out any form of sample degradation from possible explanation, especially in the form of oxidation. It implies that the size-dependent thermal hysteresis observed from the NPs is most likely to be an intrinsic effect.
Our subsequent 57 Fe NMR measurements tell a similar story too. The NMR signal measured above T V (130 K) shows a very narrow single peak for the 7 μm sample (Fig. 2a). This peak progressively gets broadened and shifted to the lower frequencies with reducing size, indicating inhomogeneous spin dynamics in the NPs due to the high surface-to-volume ratio. When measured below T V (110 K), the single peak of the 7 μm sample gets split into several ones, which is a clear sign of a new magnetic phase with quite different internal magnetic field distribution 18 . This general tendency is consistently seen for the NP samples. It is in stark contrast with the almost flat signals measured on γ-Fe 2 O 3 (maghemite) under the identical condition, which is another evidence of our experimental observations being intrinsic effect of NPs on the Verwey transition, not due to oxidation. Yet a more drastic observation is made on the temperature dependent NMR data (Fig. 2b). For the 7 μm sample, the single peak splits into weakened multiple peaks and reappears at 120-121 K upon cooling and heating with almost the same peak widths. On the other hand, the NMR data collected on the 42 nm sample display a significant thermal hysteresis; the single peak is split at 108.5 K in cooling and recovered at 119.5 K in heating with a thermal hysteresis as large as 11 K. The size-dependent hysteresis of the NMR data becomes much clearer when the broadness of the NMR spectra are plotted against the temperature (Fig. 2c). Note that the spectra broadness present in Fig. 2c is defined at quarter maximum of the strongest NMR peak at each temperature in the same frequency region. The hysteresis decreases with the size of the NPs smaller than 77 nm, as we observed from the magnetization data in Fig. 1b.
The thermal hysteresis is also registered in the temperature-dependent XRD data (Fig. 3a,b). Upon the Verwey transition, Fe 3 O 4 (440) peak becomes significantly broadened due to the change in the lattice structure. The transition temperatures of the 42 nm sample measured in cooling and heating procedures differ by 7 K, in agreement with both the magnetization and the NMR data discussed above. We also undertook RIXS experiments at near Fe L 3 -edge with σ polarization on the samples at the temperatures above (300 K) and below (40 K) T V . As shown in Fig. 3c, there is a clear peak centered at 0.2 eV for the 7 μm sample at 300 K. This peak is assigned to a polaron excitation in a recent Fe 3 O 4 single crystal study 19 . With reducing the size, this peak gets progressively suppressed and almost disappears for the 7 nm sample. For comparison, our data collected on γ-Fe 2 O 3 do not show the peak at all. When we compare the RIXS data collected at 40 and 300 K, there is no visible temperature dependence as one can see in Fig. 3d (see Fig. S5 in SI for multiple peak fittings to obtain the peak area): a recent RIXS measurement made on single crystal Fe 3 O 4 showed that the polaron peak persists well into the high temperature phase with very small temperature dependence up to 550 K 19 . Similarly, the peak area does not show much of temperature dependence either. It is to be noted that the measurements were taken with the energy resolution of 100 meV and the spectra with π polarization incident beam also show the same excitation. With all the data taken together, the size dependence of the thermal hysteresis becomes much clearer as shown in Fig. 4a and b. The hysteresis width, ΔT V = T V (heating) − T V (cooling), shows clear size dependence with the critical size of 120 nm. For the sizes larger than or equal to 120 nm, ΔT V is proportional to D −1 where D is the mean size of the NPs. On the other hand, it shows −D −3/2 behaviour for D < 120 nm. We note that this size dependence of ΔT V is very similar to the size-dependent magnetic property of the NPs 20 . As shown in Fig. 4c, the ratio of the remnant and saturated magnetization, M r /M s , from isothermal magnetization at 20 K gradually increases with reducing size, before flattening off at a value of 0.4 far below 120 nm (see Fig. S6 in SI for the raw data). At 300 K, this parameter shows a peak at the same size. This observation is consistent with the theoretical prediction for Fe 3 O 4 that the value of M r /M s for a single domain is 0.5 21 , indicating that the NPs become magnetic single domain below the critical size of 120 nm. Transition from multi-to single domain is also corroborated by the size-dependent coercivity data, H c (D), in Fig. 4d that show the maximum at around 120 nm. Furthermore, the fitting curves from the ΔT V data is remarkably similar with the coercivity of the NPs. Based on this relationship between the size dependence of ΔT V and magnetism, we conclude that the Verwey transition and magnetic domain share the same critical size.

Discussion
Let us now turn our attention to the origin of the observations we made above: (i) the size-dependent thermal hysteresis of the Verwey transition, ΔT V (D), and ii) the coincidence of the critical size of ΔT V (D) and magnetic domain. In general, the width of the hysteresis reflects the kinetic energy barrier for the phase transition. The hysteresis of 11 K (0.95 meV) at its maximum gives a rough estimation of energy correlated with charge ordering which is considered to drive the Verwey transition 8 . As we discussed above, the Verwey transition has the composite nature of intercoupling between the four degrees of freedom in solid: lattice, orbital, spin, and charge. Among them, the energy scale of spin degree of freedom (magnetism) is the smallest and about a few meV, being closest to that of ΔT V , while the others usually have larger characteristic energy scales, often of a few eV. According to the previous literature, iron oxide NPs with the size of 160 nm or smaller have a single lattice domain 22,23 . Therefore, we can safely assume that our NPs form in a single domain of all four degrees of freedom below 120 nm, and the size of magnetic single domain has dominant effect on the kinetics of the Verwey transition.
Study on the relationship between the kinetics of phase transition and the size of colloidal NPs dates back to two decades ago 11,24 . In 1997, it was reported that the high pressure-induced solid-solid phase transition of CdS NPs show size-dependent hysteresis with the kinetic energy barrier rapidly increasing with the size of the NPs from 1.5 up to 4.3 nm 24 . Apparently, this observation has strong similarity with our result that the NPs have a higher energy barrier for the phase transition compared to bulk solid and that the hysteresis width passes a maximum as the size gets reduced. For CdS NPs, it was suggested that the size-dependent hysteresis is related with the single nucleation event within a nanoparticle upon the phase transition. As long as the size of a nanoparticle is smaller than the critical size of the nucleus of another phase, a phase transition is completed by a single nucleation and the energy barrier is linearly proportional to its volume. When the particle size is large enough (10 2 -10 3 nm), multiple nucleation and growth take place within the volume, which lowers the energy barrier with respect to that of the single nucleation (see Scheme S1 in SI for the detail of the model). However, this model does not apply to our data: ΔT V (D) is not proportional to D 3 but −D −3/2 when D < 120 nm, and the crossover from −D −3/2 to D −1 behaviour is abrupt and discontinuous, rather than gradual change from single nucleation to multiple nucleation.
The significant resemblance of ΔT V (D) and H c (D) in Fig. 4b and d suggests that the coupling between charge and spin degrees of freedom is manifested in the form of metastability in the Verwey transition. For a uniaxial magnetic single domain, H c is proportional to −D −3/2 due to a balance between the energy barrier for coherent spin rotation and thermal energy 20 . In a magnetic multi-domain regime, domain wall nucleation is at work to lower the energy barrier for spin rotation. In other words, H c reflects the kinetic stability of magnetization reversal. And, in the Verwey transition, ΔT V has the same meaning for transition between charge ordering and disordering. At this point, the mechanism underlying correlation of charge and spin ordering is not understood yet. Given the nature of strong correlations between the spin and charge degrees of freedom, one can understand this similarity in the way that when domains of the charge channel get modified upon reducing size naturally the corresponding domain of the spin channel will follow the suit through the anticipated strong coupling between the spin and charge degrees of freedom. Therefore, it seems to us natural to observe such a strong similarity between ΔT V (D) and H c (D). Interestingly, an increased thermal hysteresis was also observed in thin film when Fe 3 O 4 thin films were prepared under optimized oxygen partial pressure to match the stoichiometry 25 . In this case, they cited that this thermal hysteresis is not directly related to the average domain size of film, but instead they suggest that anti-phase boundaries may affect the thermal hysteresis.
To summarize, we observed that the Verwey transition of stoichiometric Fe 3 O 4 NPs shows a very large size-dependent thermal hysteresis. For our bulk standard sample, the hysteresis width is only 1 K. In nanometer scale, however, the hysteresis gets significantly enhanced reaching the maximum of 11 K for 120 nm sample. For the NPs smaller than 120 nm, the hysteresis width is well fitted with −D −3/2 curve, similar to the size dependence of the coercivity. Also, the size of the maximum hysteresis coincides with the critical size of magnetic single domain. Overall, ΔT V (D) and H c (D) show significant resemblance. We think that this close relationship suggests a much tighter coupling than have been thought of the charge and spin degrees of freedom in Fe 3 O 4 .

Synthesis of stoichiometric Fe 3 O 4 nanoparticles.
For the precise stoichiometry control, the whole procedures were performed under slightly reducing atmosphere by using a standard Schlenk technique under a gas flow of 4 wt% CO/CO 2 15 . All of the solvents were deaerated and stored in a glove box before use. In a typical synthesis, a mixture of Fe(acac) 3 and oleic acid in benzyl ether was degassed under vacuum for 1 hr. The mixture was heated to 290 °C at a rate of 20 K min −1 and then kept at the same temperature with vigorous stirring for 30 min to complete For the NPs smaller than 100 nm, a mixture of 33.9 g of oleic acid in 312 g of benzyl ether was used with varying amount of Fe(acac) 3 from 24.5 g for 77 nm NPs to 5.33 g for 7 nm ones. For the larger NPs, 3.39 g of oleic acid in 31.2 g of benzyl ether was used with 2.34 g and 1.92 g of Fe(acac) 3 for 390 and 120 nm NPs, respectively. To prepare bulk standard Fe 3 O 4 (7 μm particles), magnetite powder from a commercial source was annealed at 1300 °C for 24 hours under CO gas flow to obtain the exact stoichiometry. To see the effect of off-stoichiometry, 7 μm sample was deliberately oxidized under air at 200 °C for six days and the NP samples in ambient air at room temperature for four days.
Characterizations. Size distribution and morphology of Fe 3 O 4 NPs were analyzed by using a JEOL JEM-2010 transmission electron microscope operating at 200 kV. Magnetic susceptibility was measured with a Quantum Design SQUID magnetometer MPMS 5XL. To scan the thermal hysteresis, a sample was cooled down to 20 K at zero field first. Then measurement was carried out at 100 Oe from 20 to 200 K and from 200 to 20 K in a row. Note that a conventional zero field cooling and field cooling curve is not same with our experiment, because our purpose of measurements is to see the thermal hysteresis of Verwey transition. The value of T V was obtained from Gaussian fit of dM/dT curves. To define the values of the remanent and saturated magnetization, we also carried out isothermal magnetization measurements from −6000 to 6000 Oe at 20 and 300 K, respectively. Cryogenic XRD measurement was conducted using a Bruker D8 Discover System with Oxford cryosystems. NMR spectra were taken using a home-made solid state NMR spectroscopy instrument equipped with a cryostat. We swept the frequency from 67.5 to 71.5 MHz under zero field at different temperatures in the range from 90 to 130 K. For both XRD and NMR data, T V was calculated by fitting derivative of the peak width with respect to the temperature with a Gaussian curve. Measurement time for each temperature, including both temperature stabilization and scan time, was 2-3 min for magnetometer, 50 min for XRD, and 1 hour for NMR, respectively. RIXS experiments were performed at 05A1 beamline of the Taiwan Light Source at the National Synchrotron Radiation Research Center (NSRRC) in Taiwan 19,26 . We used an AGM-AGS spectrometer and the scattering angle was set at 90 degree. The measurements were conducted at zero field in a high vacuum chamber (~10 −9 Torr). The incident X-ray energy was defined by X-ray absorption spectrum and fixed at Fe L 3 edge minus 4 eV (~706 eV). All samples were prepared in pelletized form with the diameter of 5 mm. We carried out the measurements using both σ and π polarized X-ray at 300 and 40 K. The experimental resolution of energy loss is about 100 meV. The Mössbauer spectra were obtained by using a transmission mode at 295 K. The Mössbauer spectrometer of the electromechanical type consists of a fixed absorber and a moving source in constant-acceleration mode with a 57 Co source of 50 mCi in a rhodium matrix. In addition, Mössbauer spectrometer was calibrated by using an α-Fe foil. The Mössbauer spectra were analyzed by a least-squares fitting procedure and provided the magnetic hyperfine field (H hf ), Isomer shift (∆), electric quadrupole splitting (E Q ), and relative area the ratio of Fe sites.