Slow oxidation of magnetite nanoparticles elucidates the limits of the Verwey transition

Magnetite (Fe3O4) is of fundamental importance for the Verwey transition near TV = 125 K, below which a complex lattice distortion and electron orders occur. The Verwey transition is suppressed by chemical doping effects giving rise to well-documented first and second-order regimes, but the origin of the order change is unclear. Here, we show that slow oxidation of monodisperse Fe3O4 nanoparticles leads to an intriguing variation of the Verwey transition: an initial drop of TV to a minimum at 70 K after 75 days and a followed recovery to 95 K after 160 days. A physical model based on both doping and doping-gradient effects accounts quantitatively for this evolution between inhomogeneous to homogeneous doping regimes. This work demonstrates that slow oxidation of nanoparticles can give exquisite control and separation of homogeneous and inhomogeneous doping effects on the Verwey transition and offers opportunities for similar insights into complex electronic and magnetic phase transitions in other materials.

O xidation is a fundamental chemical process used to tune many material properties such as changes of correlated electron behavior, as typified by the Verwey transition in magnetite. The signatures of this transition-drastic changes in the crystal structure, magnetization, electric resistivity, thermal conductivity, and other properties near T V = 125 K 1 -have been studied for over eight decades, but many aspects remain mysterious. Among these is the great sensitivity of the transition to oxidation by incorporating tiny amounts of excess oxygen as Fe 3(1−δ) O 4 or equivalent cation doping, e.g., with zinc. Oxygen doping of bulk samples suppresses T V to minimum reported values near 80 K for δ = 0.012 [2][3][4][5][6] . A crossover from a sharp firstorder transition to a much broader second-order change was observed near T V ≈ 100 K for a critical doping δ c = 0.0039 2,6 . The firstto second-order change was thought to correspond to disruption of the long-range electronic order. However, recent studies have shown that the complex arrangement of charge, orbital, and three Fe-site trimeron 7 states observed by X-ray crystallography below T V is preserved in the second-order regime and selective doping of one site was proposed to account for the changeover 8 .
Control of the tiny nonstoichiometry in bulk Fe 3(1−δ) O 4 samples is very challenging, as high temperatures have to be used to give appreciable oxygen diffusion with samples quenched to ambient conditions for characterization [2][3][4][5][6] . A further issue is that any inhomogeneity in the doping may lead to further suppression of T V as the Verwey transition is also known to be sensitive to non-hydrostatic stresses at constant doping 9 . In general, it is challenging to disentangle the intrinsic doping effect from any contribution due to inhomogeneity when correlated electron systems are chemically doped. However, nanoparticles of magnetite have been developed extensively in recent years for biomedical applications ranging from magnetic resonance imaging contrast agents to thermal therapy, to destroy cancer cells 10 . Chemical control of nanoparticle growth has also enabled fundamental aspects of the Verwey transition, such as the intrinsic domain size dependence, to be explored at the nanoscale [11][12][13] . As oxygen diffusion into nanoparticles occurs far more rapidly than into bulk material, this has enabled us to explore the effects of room-temperature oxidation of magnetite.
Here we show that slow oxidation of highly stoichiometric and monodisperse Fe 3 O 4 nanoparticles reveals an intriguing evolution of the Verwey transition. By exposure to several oxygen partial pressures at ambient but controlled temperature of 30°C, we found an initial drop of T V down to 70 K after 45-75 days, followed by a recovery up to 95 K after 160 days. We confirmed that the T V persists up to 1071 days in all experiments. A simple physical model based on Fick's law and basicassumptions about doping and doping-gradient effects can reproduce the intriguing time evolution. This demonstrates that the persistent 95 K value corresponds to the lower limit for homogenously doped magnetite and hence for the first-order regime. In comparison, further suppression down to 70 K results from inhomogeneous strains that characterize the second-order region. We expect that this work will establish an important way to control chemical doping in nanoparticles and can be further applicable to understand related phenomena for other materials having complex phase transitions.

Results and discussions
High-quality stoichiometric Fe 3 O 4 nanoparticles were prepared by a method used in previous studies 11 . The average size of nanoparticles used here is 44 ± 3 nm. We note that this sample consists of single magnetic domain particles based on comparison of the coercivity to that reported in a previous study 13 . Samples were exposed to controlled oxygen partial pressures P(O 2 ) from 0.2 to 2 atm, as described in "Methods", with magnetization, X-ray diffraction (XRD), nuclear magnetic resonance (NMR) spectra, and heat capacity measurements used to monitor changes. The magnetization for freshly prepared samples shows a sharp drop at 120 K with visible hysteresis demonstrating a first-order Verwey transition. Initial oxidation in air (P(O 2 ) = 0.2 atm) leads to broadening and suppression of the Verwey transition. The thermal hysteresis decreases and is suppressed after 13 days when T V = 103.7 K (see Supplementary Fig. 1), indicating a change from first-order to second-order behavior. Further oxidation leads to increased suppression of T V and broadening of the Verwey transition, quantified as ΔT V from the full-width at halfmaximum (FWHM) of the peak in the first derivative of the magnetization (dM/dT) up to 78 days ( Fig. 1a and see also Supplementary  Fig. 2a, b). However, further oxidation beyond 78 days led to a recovery of the Verwey transition, with T V increasing and ΔT V decreasing, until a constant value of T V = 95 K was obtained at around 160 days. This persists thereafter for measured times up to 1071 days. The variations of T V and ΔT V are replicated in measurements of other quantities such as the coercive field (H c ) (Supplementary Fig. 2h) and the heat capacity, as plotted in Fig. 1b. The heat capacity shows that the sharp transition seen in the fresh sample becomes broadened upon oxidation up to 84 days but then sharpens as T V increases until 140 days and the Verwey transition is still visible after 1071 days of oxidation with T V = 95 K. The entropy changes through the Verwey transition were estimated and showed a very similar trend to T V as shown in Supplementary Fig. 3. Similar variations are observed in XRD, through variations in the (440) peak broadening ( Supplementary Fig. 4a-c) and in Fe NMR spectra ( Supplementary Fig. 5a, b).
Temporal evolutions of these measurements are shown in Fig. 1, and magnetization, XRD, and heat capacity are shown at four different stages of oxidation in Fig. 2, to illustrate the principal changes clearly. The time dependences of T V and ΔT V are summarized in Fig. 3a, showing good agreement between T V values estimated from magnetization, XRD, NMR, and heat capacity measurements. Transmission electron microscopy images show that nanoparticle size and morphology have not changed after 92 days of oxidation ( Supplementary Fig. 6), whereas Mössbauer spectra show that the characteristic magnetite pattern of Fe 2+ and Fe 3+ signals is preserved after 102 days (Supplementary Fig. 7). Although magnetization and heat capacity measurements on the 1071 day sample confirm that the Verwey transition is still present with T V = 95 K, the changes in magnetization and entropy at T V are much smaller than for the 150-day sample, indicating that further oxidation of the nanoparticles leads to a shell of γ-Fe 2 O 3 (which has no Verwey transition) surrounding a core of The discovered variation of the Verwey transition temperature with oxidation time is surprising as a gradual oxidation process would be expected to lead to a monotonic decrease in T V to a minimum value at the maximum doping accommodated by the Fe 3 O 4 lattice. The observation that T V is initially suppressed by as much as 25 K below the final persistent value of 95 K demonstrates that the transition is suppressed not only by the doping effect from the concentration of added oxygen C but also by inhomogeneous strains created by the oxygen concentration gradient dC/dr between the oxygen-rich exterior and oxygen-poor interior of the nanoparticles. This accounts for the maximum ΔT V transition width being observed at the minimum T V . We have developed a simple model to quantify these effects using Fick's law of diffusion to simulate dC/dr, which can readily be transformed to strain (see "Methods"). The calculated total strain σ tot variation is well-matched with that of ΔT V (Fig. 4a). Furthermore, the overall variation of T V can be modeled by assuming that T V is decreased by both a doping term, proportional to C, and to a strain term that scales as dC/dr. Using parameters fitted to the ΔT V and T V values, an excellent fit to the time dependence of T V is obtained (see Fig. 4b and "Methods"). The fitted diffusion coefficient D = 2.4 × 10 −19 cm 2 /s for these Fe 3 O 4 nanoparticles is consistent with D~10 −20 cm 2 /s values reported in the literature 14,15 .
The above results demonstrate that ambient temperature exposure of magnetite nanoparticles to oxygen gas provides a simple way to slow down the kinetics of oxidation, enabling the effects of oxygen doping and oxygen-doping gradients on the Verwey transition to be separated. A striking observation is that the persistent value of ARTICLE T V = 95 K observed at the upper homogenous doping limit for magnetite is not at the lower end of the~80-120 K range found from previous studies of bulk magnetite and the present nanoparticle results. Instead, it lies close to the T V ≈ 100 K crossover between the firstand second-order Verwey transitions. This demonstrates that the crossover corresponds to the intrinsic lower temperature limit of the Verwey transition in homogeneously doped magnetite. Hence, the reported critical doping δ c = 0.0039 corresponds to the true upper limit for homogenous oxygen doping. Lower T V values down to 70 K in the second-order regime result from the additional effects of strain gradients on the transition. In the present experiments, these result from the oxygen concentration gradients between the surface and center of the nanoparticles. However, in previous experiments on bulk samples where second-order transitions were reported for higher oxygen-doping levels (0.0039 < δ < 0.012), the incorporation of additional oxygen likely led to the formation of small oxygen-rich regions within the magnetite that gave rise to the strain gradients and nucleate the formation of a secondary γ-Fe 2 O 3 phase on further oxidation beyond δ = 0.012.
This experimental approach provides a simple way for elucidating doping effects in solid oxides and related materials. Oxygen diffusion into bulk samples at room temperature is prohibitively slow, but the use of nanoparticles where diffusion lengths are reduced to tens of nanometers have enabled the inhomogeneous and homogenous doping of magnetite, and subsequent doping to maghemite, to be separated on timescales of 10's to 100's of days. The rate of oxidation here to homogenous doping level δ c = 0.0039 in 2t min ≈ 140 days corresponds to the incorporation of only~70 oxygen atoms per magnetite nanoparticle per day, leading to slow changes in both electronic doping and strain effects. Elastic strain couples to most phase transitions in solids and can lead to large differences in structural phase transition temperatures depending on whether strain is constrained by external stresses ("clamped") or relaxed to thermodynamic equilibrium ("unclamped"). In Landau theory, such changes in strain coupling can drive phase transitions between the second and first order, as observed here for magnetite. Our method of following the dynamics of slow chemical reactions may offer further new insights into how doping and strain effects tune correlated-electron properties in other materials, e.g., transformations of antiferromagnetic LaMnO 3 to ferromagnetic  T Aging Time (Days) magnetoresistive LaMnO 3+δ and of insulating YBa 2 Cu 3 O 6 to superconducting YBa 2 Cu 3 O 7 .
In summary, investigation of the Verwey electron-ordering transition during slow oxidation of Fe 3 O 4 nanoparticles up to 1071 days reveals an unusual behavior, with initial suppression and broadening of T V followed by recovery to a persistent value near 95 K. This variation is explained quantitatively by a simple diffusion model in which both doping concentration and concentration gradient effects are essential. The results account for the previously reported firstand second-order regimes of the Verwey transition as being homogenously and inhomogeneously doped, respectively. Such slow experiments, corresponding to the incorporation of only tens of atoms per nanoparticle per day, are likely to give further insights into other chemical tuning processes of correlated-electron materials.

Methods
Sample preparation and oxidation condition. Stoichiometric and homogenous Fe 3 O 4 nanoparticles were prepared using a standard Schlenk technique 11 . We used thermal decomposition of iron acetylacetonate precursor in the oleic acid surfactant. We varied the precursor-to-surfactant ratios to control the nanoparticles size and we chose 44 nm-diameter nanoparticles with a 3 nm distribution for this work. Samples were stored in an incubator for oxidation by ambient air. The temperature inside the incubator was maintained at around 30 ± 1°C and the humidity was fixed at~20%. For the experiments performed at other oxygen pressures, samples were sealed in quartz glass tubes under pressures of P(O 2 ) = 1.0, 1.1, and 2.0 atm.
Characterizations of the physical properties. Magnetic susceptibility measurements were conducted using a vibrating sample magnetometer (Lakeshore 7300 series). The samples were cooled to 30 K under zero field and then measured at 100 Oe during heating up to 150 K and cooling down to 30 K in sequence. The isothermal magnetization measurements were also carried out from −6000 to 6000 Oe at 30 and 150 K.
To further support the observations, we examined the structural change induced by the Verwey transition. The XRD experiments were performed using a commercial diffractometer (Bruker D8 Discover System) equipped with an Oxford cryosystem. We mainly focused on the (440) Bragg peak, as it is most drastically modified upon oxidation. The single (440) peak of the cubic phase splits into two peaks when transformed into the monoclinic phase below the Verwey transition. However, it can only be seen as a broad peak for nanoparticle cases, because the XRD peaks are significantly broadened due to size effects ( Supplementary Fig. 4a,  b). As shown in Fig. 1c, the drop of peak height is no longer observed after 50 days of oxidation but appears again in the data taken after 115 days of oxidation. This is reflected by the FWHM of the (440) peak ( Supplementary Fig. 4c).
To examine the oxidation effect on the spin dynamics, we measured the 57 Fe NMR spectra, which are very sensitive to the local field distribution from the Fe ions. NMR spectra of all samples were measured using a home-made solid-state NMR spectroscopy instrument with a cryostat as in the ref. 12 . The single peak located at 69.5 MHz splits into multiple peak structures through the Verwey transition for the fresh sample 12 , but the split NMR peaks merge and become a broad peak in the spectra of aged samples as oxidation progresses (Supplementary Fig. 5a, b). The entire oxidation dependence of the NMR spectra taken at 80 K is summarized in Fig. 1d, where we found that the peak splitting is not observed from 54 days to 78 days. After~2 years of oxidation, we measured the NMR spectrum and found that it was slightly more structured than that taken on a sample with 78 days of oxidation.
To measure the heat capacity of the Fe 3 O 4 nanoparticles, we used both the Physical Property Measurement System 9 from Quantum Design and a home-built cryostat system. The powdered sample was pelletized before the measurements and mounted onto the sample stage after measuring the addenda. The heat capacity was measured from 30 to 150 K under zero field as shown in Fig. 1b. To define the T V , we first subtracted the background signals using the sum of the polynomial functions. The T V was defined as the maximum of the remaining heat capacity.
The size distribution and morphology of the samples were characterized by using a JEOL JEM-2010 TEM operating at 200 kV. We could not find any significant difference between the fresh and oxidized samples from the transmission electron microscopy images (Supplementary Fig. 6).
Mössbauer spectra of selected samples were measured below and above T V in transmission geometry. We used a 57 Co radiation source and the spectra were fitted with Lorentzian functions as described elsewhere [16][17][18] . Spectra taken after 102 days of oxidation shows a similar structure to that of the fresh sample ( Supplementary Fig. 7) and fitted hyperfine fields (H hf ), quadrupole splitting (E Q ), and isomer shift are summarized in Supplementary Tables 1 and 2. These confirm that the oxidized sample is still magnetite and not maghemite or other iron oxides.
Determination of T V . We defined T V from the center of the transition shown in several measurements. In the case of magnetization, the first derivative of the magnetization as a function of temperature (dM/dT) was fitted with a Gaussian function giving T V and ΔT V as the center and the FWHM of the peak, respectively. For the XRD, we obtained the (440) peak amplitude as a function of temperature by fitting the data with a Gaussian function and estimated T V from the center of the peak amplitude drop. T V was also determined as the mid-point in the curve showing the NMR peak height drop without any fitting. T V from the heat capacity measurements was taken as the center position of the peak appeared in C/T using Gaussian function. Before the peak fitting, the background of C/T was subtracted based on the background fitting using a polynomial with the formula of a 0 T 2 þ a 1 T 3 þ a 2 T 4 þ a 3 T 5 .
Strain σ calculation from a diffusion model. In our theoretical model, the concentration gradient dC/dr is related to oxidation-induced strain. To calculate the total strain σ tot , we first simulate dC/dr from a simple diffusion model using Fick's law. The diffusion equation in a sphere is 19 : To solve the above partial differential equation, we applied the initial and boundary conditions as follows.
Here, C is the amount of oxygen diffused into the sphere, with C = 1 being fully oxidized. The boundary condition is time dependent and the C at the outer surface is assumed to increase linearly until t min and remains constant after that. By solving the diffusion equation, we obtained C and dC/dr as a function of radius r ( Supplementary Fig. 9a). For the calculation results shown in Fig. 4 and Supplementary Fig. 9a, we assumed a = 20 nm, D = 10 −19 cm 2 /s, and t min = 70 days.
As described in the main text, using the simulated dC/dr, we calculated the total strain σ tot inside nanoparticles. We assume that the dC/dr applies pressure on adjacent oxidized shells. Based on the theory of elasticity 20 , the strain σ i applied in the interface between i-th and j-th shells with hydrostatic pressure p i is written as below with the assumption that p i = p j .
We assume that the hydrostatic pressure p i is linearly proportional to dC/dr in the nanoparticle so the total strain σ tot can be calculated σ tot ¼ ∑ i σ i , as shown in Fig. 4a and Supplementary Fig. 9b. Based on the D dependence of σ tot , the diffusion coefficient was taken to be D = 10 −19 cm 2 /s. Fitting the oxidation time dependence of T V . It is known that T V is linearly dependent on off-stoichiometry parameter δ(t) in the bulk case 2,17 . Here we assume that in nanoparticles, the additional term dδ(t)/dr also affects T V (t), as the strain builds up during the oxidation as described in the main text. As δ(t) is equivalent to C(t) from the above diffusion model, we write T V (t) as: Here, a 1~a3 are arbitrary constants. As dC(t)/dr scales as strain σ tot and transition width ΔT V (Fig. 4a), we replace dC(t)/dr with ΔT V (t) and C(t) can be assumed to be an exponential decay function, so the equation becomes: Using the ΔT V (t) values from magnetization measurements, we fit the oxidation time dependence of T V (t) from Eq. (6) as shown in Fig. 4b. We found that a 0 1 ¼ 95:47 K, a 0 2 ¼ 23:01 K, a 0 3 ¼ 0:5, and τ ¼ 19:44 days. The time exponent should be related to the diffusion coefficient D as 1=τ = −Dπ 2 /a 2 from the general solution of the diffusion equation in a sphere 19 , giving D = 2.41 10 −19 cm 2 /s in excellent agreement with the assumed value of D = 10 −19 cm 2 /s above.

Data availability
The data that support the findings of this study are available from the corresponding author J.-G.P. upon reasonable request.