Atomistic study on the pressure dependence of the melting point of NdFe₁₂

The transition to green energy relies on the improvement of two areas of technology, the harnessing of energy from natural sources, and the application of that energy. In both cases electrical engines have a significant role to play, for example in the electrical engines of wind turbines, and power trains of electric cars. The most effective electrical engines contain permanent magnets on their internal rotor, which require a high magnetisation, coercivity, and curie temperature in order to operate efficiently over the range of conditions the motor may be placed in. The curie temperature is of particular importance for electrical engines in cars, which are prone to heating up to a temperature significant enough to decrease the magnetisation of the magnets and hence lower the efficiency of the engine.

Currently, in high performance electrical engines Nd₂Fe₁₄B is the magnet of choice with some of the neodymium replaced with dysprosium to increase curie temperature giving (Nd, Dy)₂Fe₁₄B.¹ Whilst the highest performing magnet currently on the market, its requirement of dysprosium to increase its curie temperature to an acceptable level significantly increases its costs due to the scarcity of dysprosium. Therefore, there has been continued research into new magnetic materials which can provide a cheaper alternative with a higher curie temperature. Recently, Miyake et al.² theoretically investigated a new magnetic structure, NdFe₁₂N, showing it to have improved magnetic properties over Nd₂Fe₁₄B, confirmed later on by Hirayama et al.³ when grown epitaxially before having its properties investigated.

However, the base structure, NdFe₁₂, is generally unstable in binary form, unable to form without a supporting structure of similar lattice properties providing it a base off of which to grow. This is obviously a significant issue for manufacturing as although Hirayama et al.³ noted that once a thin layer of the phase is deposited it continues grow without further stabilisation (at least up to 320nm), epitaxial deposition is currently a very inefficient way to bulk manufacture materials.

Therefore, in order to gain insight into what causes the instability of this material, and similar promising 1-12 phases (SmFe₁₂, SmCo₁₂), a better understanding of the thermal evolution of these structure is necessary. One of the key structural changes that any material goes through, which is particularly relevant for manufacturing, is melting, and this is the subject investigated in this paper.

To simulate the melting point of NdFe₁₂ we used the solid-liquid coexistence method⁴ in LAMMPS,⁵ using Morse potentials⁶ for our interatomic interactions. The solid-liquid coexistence method is a molecular dynamics technique which allows you to simulate accurately the melting temperature of a material by bringing a melted and solid phase of the material together at a constant temperature and pressure in an NPT ensemble, before allowing the system to evolve. The systems evolution gives an indication of its proximity to the melting temperature (MT), liquid phase growth indicating > MT, a stable equilibrium in which the atoms of the liquid phase have become an amorphous solid indicating <MT, and a stable equilibrium in which the liquid phase remains liquid indicating =MT (close to the correct temperature). Figure 1 shows the base system and the three evolved states.

Specifically, in our simulations, a complete supercell of the simulated structure containing ∼10,000 atoms is instantiated with periodic boundary conditions, and allowed to equilibrate at 1000K in an NPT ensemble for ∼80 picoseconds or 80,000 time steps.

One half of the total supercell has the velocity of, and force on its atom set to zero, effectively freezing the atoms in place, and making them non-interactive with the other half of the cell, whilst the other is heated well above its melting point in an NVT ensemble to give one totally melted and one frozen half of the structure.

From this state, the two halves are brought back together. However, in order to avoid the characteristic shock waves which, occur when the two halves are instantaneously allowed to interact again, the atoms present in the regions, in which the two halves are in contact, have extra spring forces attached to them. These spring forces tether the atoms near to the positions they are currently in, by exerting a force on them toward their original position when they stray from it. This force is proportional to the distance they are from their original position.

In this way, the atomic movement gained as a consequence of the two halves being brought together can be damped out of the simulation, resulting in an increase of accuracy. Using this technique, the system is equilibrated in the NPT ensemble at the desired pressure and temperature, as the springs are gradually removed over 135 picoseconds or 135,000 time steps. Once the springs are removed the whole system is allowed to evolve over a further 250 picoseconds, it’s evolution can then be analysed to find out whether the melting point is at above or below the chosen temperature.

Using the method above we began by using a set of Morse potentials derived from Chens inversion lattice technique.⁷ However, on investigation, we found the potentials gave highly inflated results of the well-known melting temperatures for the constituent elements of NdFe₁₂, neodymium and iron, along with incorrect values for various well known constants such as the cohesive energy, and Bulk modulus. In order to correct this, we used genetic algorithm coding and the General Utility Lattice Program (GULP)⁸ to refit the potentials to a series of experimentally derived constants found in the literature, see Table I.

An iterative process was used to fit the potentials, proceeding in generations in which a series of genes containing values for the base parameters of the Morse potential are varied via random single percentage point scaling, and interchanging of similar parameters – known as mutation, and cross-breeding respectively in the genetic algorithm literature. Each gene has its fitness determined by calculating the percentage difference between the constants gained from its use in simulation, and those expected from the literature.

These differences are weighted and summed to give an overall value of fitness, the lower the value, the fitter the gene, the more accurate the potential. Weighting of the differences allows the preferential selection of potentials that closely match constants considered of greater importance, in this case cohesive energy. In order to converge to a solution, the algorithm discards genes of poor fitness at the end of each iteration, and in the next iteration replaces them with modified versions of the fitter genes before repeating the process. Over a number of generations this provides potentials that more accurately fit the desired constants, the constants gained from the final potentials chosen as a result of this method are shown in Table I, with comparison to their expected values.

The Morse potentials for the Nd-Nd and Fe-Fe interactions were fitted first, using the base elemental unit cells of neodymium and iron respectively. This has two benefits, firstly the structures have both been thoroughly studied, and therefore accurate constants with which to fit the potentials can be found for them, and secondly the only interatomic interaction that occurs between atoms in either structure is the one we are interested in fitting. As mentioned in order to fit the potentials we used constants from the literature shown in Table I. The genes which contain parameter sets for the potential were instantiated so that each parameter varies over a wide range of possible values. As each parameter in the set varies independently a very large amount of the potentials parameter space can immediately be explored by instantiating the genes in this way. Parameters for each potential were varied by cross breeding and mutation, then used for the interatomic potentials in GULP simulations. Constants resulting from these simulations are used to calculate potential fitness. Both cross breeding and mutation occurred pseudo-randomly across the generations. The probability of a cross breeding event followed a sinusoidal curve within an exponential envelope function which slowly reduced to zero. Mutation event probability remained constant, and instead the percentage by which a mutation varied a genes potential parameters decreased exponentially towards a constant none zero value. Across hundreds of generations, the algorithm slowly moved towards regions of good fitness within the parameter space.

At the end of each generation a king, which is the potential of highest fitness for that generation, was chosen, and if no potential from the proceeding generations improved on its fitness it was kept as the best fit potential. After hundreds of generations the final king was used as the potential for future simulations.

In order to fit the final Nd-Fe potential two reference structures Nd₂Fe₁₇,¹⁴ and NdFe₁₂² were used. However, as these structures have not been studied as thoroughly, only their lattice constants, comparative cohesive energies per atom, and a reasonable estimate of a similar Bulk modulus to Nd₂Fe₁₄B were used for parameter fitting. Comparative cohesive energy is the per atom energy difference between Nd₂Fe₁₇ and NdFe₁₂, and is set up to be slightly in favour of Nd₂Fe₁₇, the more stable of the two structures. Identical atomic interaction is modelled by the previously fit potentials, so that only the potential of interest was varied to produce a good fit.

The results of the subsequent melting simulations are shown in Figure 2. As it shows, our calculated melting temperatures of the pure elemental forms of neodymium and iron are ∼1700K, and ∼3000K, overestimated by a factor 1.32, and 1.66 respectively, indicating firstly that the use of simplistic Morse potentials is not effective for quantitative melting analysis, and secondly that we should expect a similar scaling of the NdFe₁₂ melting temperatures. Accounting for this and bearing in mind that interactions within the NdFe₁₂ structure are dominated by Fe-Fe interactions, we estimate the ambient pressure melting temperature of NdFe₁₂ to be ∼1200K, reduced from the calculated value of 2040K. As can be seen in Figure 2, increasing the pressure of the simulation increases the structures melting temperature to 2110K at 500MPa, which would reduce to ∼1271K.

Whilst there are no experimental melting temperatures to compare our results to, it is interesting to compare our calculated melting point to the reported annealing temperatures of similar ternary structures NdFe₁₁Ti, and NdFe₁₁Mo, which are 1373K,¹⁵ and >1173K¹⁶ respectively. Indicating that the melting temperatures of these materials is likely on the order of hundreds of Kelvin greater than that calculated for NdFe₁₂. This is not surprising given that the addition of either ternary element to the binary NdFe₁₂ structure is known to stabilise the phase.

Although these comparisons of melting temperatures are interesting, in reality it is likely that NdFe₁₂ would decompose into alpha-iron and a stable rare earth intermetallic (Nd₂Fe₁₇), at a lower temperature than the melting points indicated here, due to NdFe₁₂ being a metastable structure. Consequently, another drawback of our model, is its inability to predict such decomposition.

Bearing in the mind the aforementioned problems, it is important to clarify this investigation was performed both as a first step towards a fuller understanding of the phase behaviour of NdFe₁₂, and the development of a methodology that could be used to computationally investigate a whole range of magnetic structures. Development of such a methodology is of high importance, as it would allow inexpensive computational investigation to aid in the formation and refinement of manufacturing protocols for new materials.

Currently however, it is clear that our method does not provide wholly accurate results, likely derived from our choice of the Morse interatomic potential, which simplistically assumes that the force between two atoms is independent of the local environment in-between them. Therefore, to improve this methodology it is necessary to develop, based on the harmonic potentials derived here, a set of MEAM potentials, which will include a correction field for each position and interaction in the 1:12 phase. In particular, further research into the NdFe₁₂ structure will look at how the melting temperature, and structural properties of the material change with substitution of some of the iron atoms for titanium.

This work was funded by the ESPRC grant EP/P02047X/1 “Picosecond Dynamics of Magnetic Exchange Springs.” Thank you to the Toyota Motor Corporation for sponsoring this work.