Spin-lattice dynamics of surface vs core magnetization in Fe nanoparticles

Magnetic nanoparticles (NPs) are fundamental components for applications such as catalysts or biomedical materials.^1–3 Developing numerical atomistic tools allowing for accurate predictions of the temperature dependence of the NPs' magnetization dynamics, as well as other relevant quantities of interest such as heat capacity, is of practical design interest.⁴ Such tools could bring insight into the fundamental processes at stake and help tailor NPs by locating optimum properties, including size and shape, for technological⁵ and biomedical applications.⁶ 

Our numerical effort is based on the classical spin-lattice dynamics (SLD) approach.⁷ SLD is an atomistic method assigning classical spins to the magnetic atoms in the simulated system.^8–10 The potential energy is computed by combining a purely mechanical potential with a Heisenberg exchange interaction. Through the interatomic dependence of the exchange interaction, the magnetic spins are coupled to the atomic potential energy. Thus, the corresponding molecular dynamics simulation of such a system allows us to calculate simultaneously the influence of (spin-derived) forces on the atoms' positions and the changes in the precession of the spins by the magnetic field set up by the surrounding atoms. Such spin-lattice simulations have been used until now to describe the temperature dependence of magneto-mechanical properties, phase transitions, phonon dispersion, demagnetization experiments, and other phenomena.^11–15 

Dos Santos et al. showed that SLD could be used to compute the temperature dependence of the magnetization for iron NPs.¹⁶ In this former study, the influence of thermal fluctuations on the total magnetization was accounted for both the lattice and the spin systems, since their coupling is evaluated in a time-dependent way by the SLD simulations. However, this work used fixed size for the atomic moments,¹⁶ thus ignoring their norm fluctuations due to surface effects.¹⁷ Those results yielded good qualitative agreement for the NPs' magnetization vs temperature trends, but quantitative agreement with experiments was obtained only if the magnetization was re-normalized. Since the total magnetization of NPs is a strong function of the NP diameter for small sizes, this situation is unsatisfactory.

In this work, we present an approach allowing us to account for surface effects on the magnetic moments. We explore its effects on the magnetization of iron NPs and display good quantitative agreement with experimental measurements. We also show that our approach can be straightforwardly used to compute the heat-capacity of magnetic NPs, accounting for both lattice and magnetic contributions.

We investigate spherical NPs with a diameter d between 1 and 8 nm. Considered NP diameters are small compared to typical domain-wall thickness measured in iron,¹⁸ allowing us to treat them as single-domain magnets. The spheres are cut out from bcc Fe crystals and relaxed for 10 ps. In our SLD implementation, this means that the atoms' positions are allowed to thermally expand and keep pressure to zero, while the spins respond to the changed environment. The atoms interact through forces computed from the spin Hamiltonian and the Chamati et al. potential¹⁹ (see the supplementary material for more details). This potential was designed to accurately compute surface energies, and former studies leveraged it to investigate finite size systems.^20,21

Initially, each classical spin vector $si$ (of unit length) points in the [001] direction. The total spin S of a NP is given by the magnitude of the vectorial sum over all atomic spins $si$ as

Thus, at temperature T = 0, S = 1.

Each atom i carries a local magnetic moment (whose magnitude is determined by the atom's environment) of size μ_i and direction $si$⁠. The average magnetic moment per atom in the NP is then given by

which gives us the magnetization of the considered NP.

The interaction between the atomic spins may change the spin direction and also influence the atomic motion. It is dominated by the exchange interaction, $−J(rij)sisj$⁠. The space dependence of the exchange interaction is described by a Bethe–Slater function, fit to the results by Pajda et al.,²⁵ since it was proven to describe well the magnetism of Fe NPs.¹⁶ J(r) is the only position dependent term in this interaction and is responsible for energy exchange between the spin and the atomic systems. Cubic magneto-crystalline anisotropy is also included as in Ref. 26, and details are given in the supplementary material [see Eqs. (S2) and (S3)].

Our time step is set to 1 fs. We use Langevin thermostats to equilibrate our NPs at the desired temperatures for a time of 0.5 ns. After this time, the magnetic properties are measured as an average over 0.3 ns. Further details on our SLD approach are provided in Ref. 16.

We assign a magnetic moment to each iron atom depending on its local environment. Figure 1(a) displays the magnetic moment of Fe atoms in compressed or expanded bcc lattices, and how magnetism vanishes at high atom densities and increases toward the free atom value (⁠$4μB$⁠) for expanded lattices.²⁷ In a zero-pressure bcc lattice, $μequ≈2.154μB$⁠. We thus use the local atomic volume in order to assign local magnetic moments to Fe atoms.²⁸ Our work only uses the local volume, but more detailed information about the local atomic coordination might be needed for a more precise assignment of local moments.²⁹ We also note that thermal expansion is more pronounced in NPs than in the bulk, and the average atomic volume increases strongly with decreasing diameter (see Fig. S1 in the supplementary material).

The local electronic density—as obtained from the EAM-type potential—carries information about the local atomic volume, see Fig. 1(b). We thus have an algorithmic mapping between the local electronic density (readily available at each step of the SLD simulation) and the size of the local moments. We fit our data to a function of the form^24,28

For the Chamati et al. potential,¹⁹ we obtain $C=(3.457±0.259)μB$ and $γ=0.414±0.073$⁠; ρ_c denotes the critical electron density at which ferromagnetism vanishes.

Figure 2 exemplifies our results on the radial dependence of magnetic moments $μ(r)$ for a NP of diameter 2.3 nm. The data can be divided in two categories,

• a core region extending up to $r≈9$ Å, in which the magnetic moment is close to the bulk value and

• a shell of width $Δr≈2.5$ Å, in which the moments increase linearly with r up to values of 2.6–2.65 μ_B.

This feature holds also at elevated temperature: at 600 K, the atomic moments calculated through Eq. (3) are only slightly higher (but within the uncertainty limits), since thermal expansion increases the atomic volume, as seen in Figs. S1 and S3 of the supplementary material. We note that this separation into core and shell regions was found for all NP sizes investigated here; in particular, the shell width is always equal to 2.5 Å. This is plausible, since it corresponds to a roughly monatomic shell. The second outer atomic shell feels an environment that resembles bulk Fe, and hence the magnetic moments take values close to bulk atoms.^30–32 All SLD calculations are performed leveraging the SPIN package of LAMMPS.^26,33 This SLD implementation considers normalized spin vectors when computing the exchange interaction. Therefore, in order to take into account the influence of the atomic moments' distribution, which varies with respect to the atomic volume and the neighboring spin configuration,^17,34 atoms were divided into two separate groups: (i) the core with moment $μcore$ and (ii) the shell with moment $μshell$⁠, where the thickness of the shell has been fixed to 2.5 Å. As can be seen in Fig. 2, this only allows us to approximate the continuous magnetic moment fluctuations. We choose $μcore=2.15μB$ and $μshell=2.45μB$ as the average over the core and shell atoms, respectively, in agreement with Fig. 2. Both the shell width and average shell magnetic moments proved valid for all cluster sizes (⁠$≤8$ nm) and temperatures (⁠$T≤900$ K) (see Figs. S2 and S3 in the supplementary material). NPs larger than 8 nm are expected to behave similarly, see Fig. S3.

This grouping influences the spin-lattice dynamics only via the exchange interaction. While core spins interact via $J(rij)sisj$⁠, the factor J is scaled by $μshell/μcore$ for core–shell spin interactions, and by $(μshell/μcore)2$ for shell–shell interactions. Former studies have been reporting bulk exchange interaction calculations for iron,^35,36 but we are not aware of similar computations for cluster shells.

Experimental data are available in Ref. 37 for Fe NPs at 120 K with a number of atoms below N = 800, corresponding to diameters below 2.6 nm. Figure 3 shows that our SLD simulation results—denoted by $μshell=2.45$—underestimate the experimental values by up to approximately 25%. While the magnetic moment increases with decreasing the NP size, it only leads to magnetic moments of around $2.25μB$⁠, even though the number of core atoms is smaller than the number of shell atoms for these small clusters, N < 100. An analysis of our data shows that this is caused by the misorientation of the atom spin for these small clusters: while the core spin assumes values of S = 0.94, the shell spin is below 0.92, and hence the total magnetization of NPs does not increase much at T = 120 K. Note that a simulation with all magnetic moments fixed to their bulk value $μequ$ shows a size independent magnetic moment of around $2.0μB≃0.92μequ$¹⁶ (open circles in Fig. 3).

Experimental data reach moments up to $3μB$ for small NPs (⁠$N≲100$⁠), whereas our larger magnetic moment values are approximately $2.8μB$⁠, obtained for atoms in edge positions on the NP surface. This demonstrates that our model underestimates the magnetic moment of shell atoms. Even deviations from spherical shapes (which may occur in experiments) would not produce larger magnetic moments in our approach. Those high moments in small clusters have been investigated using tight-binding calculations³⁹ and have been assigned to structural changes in small NPs favoring fcc-like and icosahedral coordinations connected to a re-ordering of atomic orbitals and associated spin changes. Such structural changes are beyond the scope of this work (where a bcc bulk-like structure is always assumed).

We therefore tested a second model, in which $μshell=3.0μB$⁠. Figure 3 shows that this calculation gives better agreement with the experimental data (reducing the error to approximately 10% for the smaller NPs). Figure 3 also shows reasonable agreement between our simulations and other experimental results for 1.6 nm clusters.³⁸ However, even for the smallest clusters, the magnetization remains below $3.0μB$ for the reasons discussed above: the core moments are only $2.1μB$ and the thermal fluctuations disorder the low-coordinated shell spins at 120 K.

A discussion of the temperature dependence of the magnetization brings us further insight. Figure 4(a) reproduces experimental data for 2.3 nm diameter clusters (N = 533). Our core–shell simulations with $μshell=2.45μB$ are in fair agreement with the experiments and only slightly overestimate the magnetic moment for intermediate temperatures (⁠$T=400−600 K$⁠). However, the calculations with $μshell=3.0μB$ strongly overestimate the magnetic moment for all temperatures, indicating that the spin-spin coupling is too strong in this model.

Motivated by the results of Billas et al.³⁷ (a large magnetic moment and a relatively low magnetization), we set a third model [denoted by $J*(r)$], in which we retain the high value of the shell moments, $μshell=3.0μB$⁠, but reduce the exchange coupling of shell atoms with other shell or core atoms to half their value. Thus, the function J(r) by Pajda et al.²⁵ is scaled by $12μshell/μcore$ for core–shell spin interactions, and by $12(μshell/μcore)2$ for shell–shell interactions.

When compared to our simulation results, this indicates that exchange between shell–shell and shell–core spins is weaker than we assumed before. Here, we find that the 1/2 re-scaling factor of the shell–shell and shell–core J(r) couplings gives excellent agreement with experiments. An alternative approach could tune the exchange coefficients, and advanced fitting methods could be employed.⁴⁰ Previous simulations use the same exchange for the core and the shell^41,42 but, given the current lack of theoretical or experimental guidance we take surface exchange as a free parameter, as in previous studies.^43,44

Figure 4(a) demonstrates that this choice nicely reproduces the experimental data for the temperature dependence of the 2.3 nm diameter clusters. Compared to the calculations with unchanged J, the magnetic moments of NPs at 120 K have been reduced for all NP sizes, thus deteriorating the comparison with experiments, see Fig. 3. As discussed above, this disagreement is likely caused by structural changes in the NPs, and therefore out of scope of the present calculations.

Locating the Curie transition on magnetic nanoclusters is challenging. Heat capacity (C_p) measurements give an indirect determination of this phase-transition temperature.^45,46 C_p is computed by taking the internal energy's derivative with respect to temperature. (Details are provided in the supplementary material.) Our simulation results for the 2.3-nm NPs, Fig. 4(b), show C_p peaks near 500 K (associated with the ferromagnetic-paramagnetic phase-transition). This is a considerable shift compared to the bulk value, consistent with our magnetization vs temperature results. The C_p curve for the NP with homogeneous values of μ and J(r) shows a similar peak, but shifted to higher temperatures since those shell spins require more energy to disorder than the shell–core case with $J*$⁠. We note that frozen-lattice spin dynamics simulations give different positions and shapes of the specific maximum, indicating that coupled spin-lattice simulations are required to correlate C_p and the Curie temperature. These results could be especially relevant for biomedical applications, since a precise estimation of C_p is essential to evaluate heating of NPs for applications like magnetic induction hyperthermia.⁴⁷ This is particularly important as lattice-only simulations do not account for magnetic degrees of freedom and thus cannot recover this intermediate temperature C_p peak.⁴⁸ Most “atomistic spin dynamics” (ASD) simulations use a fixed perfect lattice, ignoring the thermal motion of the atoms. Dos Santos et al.¹⁶ studied the importance of moving atoms for 3 nm NPs and found that frozen-lattice simulations lead to unphysical compressive stresses of 0.2 GPa. Figure S4 in the supplementary material compares our magnetic-moment results to those obtained by a frozen-lattice calculation and shows better agreement with the experiment for the moving lattice case.

Figure 5 displays how the NP spins become disordered with temperature. The total spin S averaged throughout the NP, and over the shell and the core only are plotted for 2.3 nm diameter clusters. The distinction between the core and shell atoms and their different magnetic moments and couplings has a strong influence on the spin disordering. The model with $μshell=3.0μB$ shows a strong coupling of the shell and core atoms, which is not unexpected, since for this cluster size, more than 51% of all atoms are in the shell. In this model, the low coordination of surface atoms is compensated with higher magnetic exchange coupling among them. On the other hand, spins show a clearer distinction between the shell and core values for the $μshell=2.45μB$ model and even more for the $J*(r)$ model, demonstrating the decreased coupling between the moments in these models.

Because of their higher coordination, the core atoms remain magnetically stiffer to higher temperatures than shell atoms. In previous work,¹⁶ it was shown that the main influence of the atomic thermal motion on the magnetization is felt at temperatures beyond 400 K and leads to a substantial decline in magnetization; this influence is included in our present calculations.

We used SLD to perform a core–shell investigation of magnetic NPs. Calculations of the influence of the atomic volume on the magnetic moment based on a local bcc-like coordination show that the surface shell, in which atomic moments deviate from their bulk value, is only 2.5 Å thick, and atomic moments there assume an average value of $μshell=2.45μB$⁠. We greatly improved agreement with available experimental data³⁷ by increasing the shell moment to $μshell=3.0μB$ but decreasing the exchange interaction among shell atoms and between the shell and core atoms by 50%. We showed that a core–shell model has sufficient flexibility to represent the magnetization of Fe NPs with SLD and to accurately compute quantities of practical design interest such as the temperature and the NP size dependence of the magnetization or heat-capacity.

We believe a major approximation of our classical model consisted in treating the exchange coupling J(r) as temperature independent; including this dependence might lead to more precise predictions,^36,49,50 improving agreement in Fig. 3, while retaining the excellent agreement in Fig. 4(a). In addition, we have considered the magnetic cubic anisotropy identical for all spins,¹⁶ although it should vary near the surface due to the coordination reduction.^51,52 Different values for NPs having different core/shell fractions should be considered (with an expected reduction of about 10% when halving the cluster size⁵³). Former studies also leveraged the Néel pair interaction to reproduce magneto-elastic and surface anisotropic effects^15,54 but using it for SLD simulations of NPs remains to be investigated.

In future work, the possibility of assigning arbitrary magnetic moment values for individual atoms during the calculation of exchange interaction will be investigated. This would allow us to reproduce the moment fluctuations displayed in Fig. 2 and should improve our agreement with experimental predictions. Complementary ab initio computations of finite size effects on both the size of moments and the effective exchange interaction between atoms for small clusters could also be designed.

See the supplementary material for additional information on the dependence of the average atomic volume of NPs on their size, the cubic magneto-crystalline anisotropy, a comparison of our calculations to frozen lattice simulations, the Hamiltonian used, and heat capacity calculation details.

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project No. 268565370–TRR 173 (Project A06). Simulations were performed at the High Performance Cluster Elwetritsch (Regionales Hochschulrechenzentrum, TU Kaiserslautern, Germany). E.M.B. thanks support from SIIP-UNCuyo 06/M104 and ANPCyT PICTO-UUMM-2019-00048 grants. J.T. acknowledges that Sandia National Laboratories is a multimission laboratory managed and operated by the National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.