Temperature-dependent phonon spectra of magnetic random solid solutions

A first-principles-based computational tool for simulating phonons of magnetic random solid solutions including thermal magnetic fluctuations is developed. The method takes fluctuations of force constants due to magnetic excitations as well as due to chemical disorder into account. The developed approach correctly predicts the experimentally observed unusual phonon hardening of a transverse acoustic mode in Fe–Pd an Fe–Pt Invar alloys with increasing temperature. This peculiar behavior, which cannot be explained within a conventional harmonic picture, turns out to be a consequence of thermal magnetic fluctuations. The proposed methodology can be straightforwardly applied to a wide range of materials to reveal new insights into physical behaviors and to design materials through computation, which were not accessible so far. Researchers have developed a numerical method to predict the atomic vibrations of disordered magnetic solids above zero temperature. Accurate methods for computing the vibrational modes of materials at different temperatures are needed to determine thermodynamic properties from first principles. However, for magnetic alloys the effects of the non-uniform distribution of atoms and quantum magnetic behaviour each require special techniques for calculations to be tractable. Yuji Ikeda and colleagues from Japan, Germany and the Netherlands have designed a scheme that combines an approach for ordered magnetic systems with established approximations for disordered materials. Their method successfully captures the unusual thermal expansion properties of Invar alloys, which are strongly affected by the interplay of atomic and magnetic fluctuations, and should provide insights into the physics of other complex materials such as high entropy alloys.


INTRODUCTION
Magnetic random solid solutions represent a large and important class of crystalline materials ranging from structural materials such as steels, [1][2][3][4] including Invar alloys, [5][6][7][8][9] up to multicomponent magnetic high-entropy alloys. [10][11][12] The simultaneous presence of chemical disorder and thermal magnetic fluctuations as well as their couplings to lattice vibrations play pivotal roles in many of these alloys. Lattice vibrations largely dominate thermodynamic properties of materials 13 and contribute to phase stability, 14 which is a key parameter for the computational design of new and innovative materials. A computational scheme that simulates the lattice vibrations in magnetic random solid solutions by properly taking into account both, magnetic fluctuations as well as chemical disorder (as sketched in Fig. 1), is therefore of genuine importance.
The delicate interactions between lattice vibrations, chemical disorder, and thermal magnetic fluctuations can cause extreme and unusual physical properties. A prominent example is the hardening of a transverse acoustic phonon mode and elastic constants with increasing temperature in Invar alloys. 15,16 Since thermal expansion-usually dominating the temperature dependence of phonon modes-is in such alloys negligible, [5][6][7][8][9] the inclusion of explicit temperature-dependent excitations, such as magnetic fluctuations, is critical to resolve such peculiarities.
We therefore propose a first-principles-based method to calculate lattice vibrations of magnetic random solid solutions, which addresses both, thermal magnetic fluctuations as well as chemical disorder simultaneously. This is achieved by extending and combining the above-mentioned disjunct approaches into a methodological framework, which allows to predict the temperature-dependent phonon spectra of magnetic random solid solutions. To achieve this goal, we apply a two-step procedure and adiabatically decouple the (fast) magnetic and the (slow) chemical degrees of freedom. To include thermal magnetic fluctuations, we utilize force constants (FCs) which implicitly depend on the temperature-dependent magnetic state. 22 These FCs are obtained using the spin-space averaging method 21 in combination with quantum Monte Carlo simulations for an effective Heisenberg spin Hamiltonian. To account for chemical disorder, which induces variations of atomic masses and FCs among atomic sites, we employ the ICPA 31,32 and the band unfolding. [33][34][35][36]

RESULTS
We first calculate the FC s for the two extreme magnetic states: (i) the low-temperature ideal ferromagnetic (FM) state, where all magnetic moments point to the same direction, and (ii) the hightemperature ideal paramagnetic (PM) state, where the magnetic moments are fully disordered. The FCs in the ideal PM state are obtained using the spin-space averaging method, 21,28 which consists in a statistical average over a large set of randomly distributed collinear magnetic moments (see Sec. B in the Supplemental Material for further details). The FCs in the magnetic state at temperature T, which has magnetic short-range order, are then calculated considering the impact of the magnons as 22 where Φ M−M′ (T) denotes the FCs between the pair of chemical elements M and M′ at T, and Φ FM MÀM 0 and Φ PM MÀM 0 denote the FCs in the ideal FM and in the ideal PM states, respectively. Here, the original formalism for chemically ordered systems 22 is extended to random solid solutions by considering the element-resolved FCs. The mixing parameter α(T) is directly related to the magnetic energy 22 (see Sec. C in the Supplemental Material). Spin quantization effects for the magnetic energy and hence for α(T), being critical below T C , 47 are incorporated by performing numerically precise quantum Monte Carlo simulations for an effective nearest-neighbor Heisenberg spin Hamiltonian. 22,[47][48][49] Here the solid solution is modeled by randomly distributing the chemical elements having different spin values onto the magnetic sites. Note that the reverse impact of vibrations on the spin Hamiltonian, and hence on the magnetic energetics and α(T), is considered to be small and, therefore, not included in this approach. It is also worth emphasizing that classical Monte Carlo simulations, which do not include spin quantization effects, would result in inaccurate magnetic energetics below T C 47 and hence inaccurate α(T).
The obtained Φ M−M′ (T) are next employed to derive the phonon spectrum of the magnetic random solid solution at any given T. In principle both the ICPA 31,32 and the band unfolding [33][34][35][36] can address the variations of atomic masses and FCs among atomic sites due to chemical disorder. The ICPA is a Green's-functionbased method and analytically incorporates the variations of atomic masses and FCs; a fully random solid solution is rigorously modeled in the ICPA. This method, however, requires the numerical solution of a relatively complex set of equations, 31,32 making its extension to multicomponent systems challenging. In contrast, the band unfolding is a supercell-based method and can therefore be straightforwardly extended to the multicomponent systems. This method suffers, however, from the undesired periodicity originating from the limited size of the supercell, which may cause spurious features in the computed phonon spectra. To eliminate them, a much larger effective supercell is constructed, and the FCs Φ M−M′ (T) computed from the original (smaller) supercell are assigned to the atomic sites of the effective supercell. Such an effective supercell has longer periodicity than the original supercell and therefore includes a larger number of distinct local configurations of chemical elements.
We apply the developed approach to two experimentally wellstudied magnetic alloys, namely to disordered face-centered cubic (fcc) Fe 0.72 Pd 0.28 and Fe 0.72 Pt 0.28 alloys. Both Invar alloys reveal the aforementioned characteristic phonon hardening of a 110 h i (in fractional coordinates for the conventional fcc unit cell) transverse acoustic mode when heated up above the Curie temperature, T C . 15,50 In the following we focus on (i) the impact of mass and FC variations among atomic sites as well as the performance of the ICPA and the band unfolding to incorporate it and, in particular, (ii) the impact of thermal magnetic fluctuations on the phonon spectra.
In order to discuss the impact of mass and FC variations among atomic sites separately from the impact of thermal magnetic fluctuations, we first focus on the ideal FM state [α(T = 0) = 1], i.e., in the absence of thermal magnetic fluctuations. Figure 2 shows the phonon spectra of Fe 0.72 Pd 0.28 and Fe 0.72 Pt 0.28 in the ideal FM state calculated using the band unfolding (top panels) and the ICPA (middle panels). The results obtained from the two different methods are very similar to each other for both alloys, although the formalisms of the ICPA and of the band unfolding are different as detailed above. Both the methods probably incorporate the impact of chemical disorder into the calculations of phonon spectra in a similar quality. The peak positions of the computed phonon spectra are also found to be in good agreement with experimental phonon frequencies.
The spectra also show phonon broadening originating from the variations of atomic masses and FCs among atomic sites due to the chemical disorder. For Fe 0.72 Pd 0.28 this effect is relatively small, while for Fe 0.72 Pt 0.28 large phonon broadening particularly in the frequency region around 3-5 THz is observed. The large phonon broadening in Fe 0.72 Pt 0.28 occurs probably due to the large differences of atomic masses and FCs among the chemical elements. The atomic mass of Pt relative to Fe (≈3.5) is much larger than that of Pd (≈1.9). The FCs in Fe 0.72 Pt 0.28 are also largely different among the distinct combinations of the chemical elements compared with those in Fe 0.72 Pd 0.28 (see Sec. E in the Supplemental Material).
To elucidate the impacts of the variations of atomic masses and FCs among atomic sites more clearly, the phonon dispersion relations are also calculated in the absence of FC and mass fluctuations using the concentration-weighted average atomic mass m and the crystallographically-symmetrized FCs Φ. The FCs Φ are obtained by first applying each symmetry operation of the fcc structure to the original FCs and then taking the average of the transformed FCs irrespective of the chemical elements. The results are shown in the bottom panels in Fig. 2. Note that by construction, no phonon broadening is obtained in this case. For the Fe-Pd alloy, the phonon frequencies derived from m and Φ still agree reasonably well with those in experiments. For the Fe-Pt alloy, however, the deviations are significant, in particular around the X point. This indicates the importance of taking the variations of atomic masses and FCs among atomic sites into account for accurate phonon computations of random solid solutions. Having verified the importance of an appropriate treatment of chemical disorder in random solid solutions, we next analyze the impact of thermal magnetic fluctuations on the phonon spectra from the viewpoint of their temperature dependence. The temperature-dependent spectra are obtained using Φ M−M′ (T) from Eq. (1). Calculations are carried out at three representative characteristic temperatures: below, near, and above T C (575 K for Fe 0.72 Pd 0. 28 15 and 367 K for Fe 0.72 Pt 0. 28 50 ). We focus in particular on the lowest-frequency phonon branch along the 110 h i direction, which shows the unusual hardening with increasing temperature. Note that within the harmonic approximation and in absence of explicit temperature-dependent excitations, phonon spectra only implicitly depend on temperature via the volume expansion. However, as mentioned above, the considered Invar alloys reveal negligible expansion below T C . A conventional harmonic approximation would therefore predict a temperature-independent phonon spectra and cannot explain the experimental data. Figure 3 shows the results employing the band unfolding (Note that similarly to the discussed ideal FM state (Fig. 2), the ICPA results are again in excellent agreement with the band-unfolding results for the intermediate temperatures and therefore shown only in the Supplemental Material.) For comparison, the results for the ideal FM and the ideal PM states are also shown. The peak positions of the computed temperature-dependent phonon spectra are in agreement with experimental phonon frequencies at the corresponding temperatures, while the results for the ideal FM and the ideal PM states provide the lower and the upper bounds, respectively. This indicates that thermal magnetic fluctuations are crucial to reproduce the experimental data.
Our simulations further indicate the presence of considerable magnetic short-range order (i.e., a finite value of α(T) in Fig. 3) even above T C , similar to the one observed for pure Fe. 22 Without taking magnetic short-range order into account, i.e., by limiting the simulations solely to the ideal PM state, the agreement with experiments is lost. Magnetic short-range order hence plays a crucial role, not only for pure Fe, but also for Fe-based alloys. We also observe slight deviations between our simulations and the experiment above T C . Since the Invar effect is lost at the temperature above T C , these deviations may be attributed to the increase of the lattice constants, which is not yet included in our calculations.
The temperature dependence of the phonon modes may be also related to the martensitic transformation in these alloys. As observed in Fig. 3, the slope of the lowest-frequency branch along the 110 h i direction around the Γ point largely depends on temperature for both Fe 0.72 Pd 0.28 and Fe 0.72 Pt 0.28 alloys and decreases with decreasing temperature. The slope is more than four times smaller in the ideal FM state than in the ideal PM state for both the alloys. Since the slope of this branch is proportional to ffiffiffiffi ffi C 0 p , 51 where C′ is the elastic stiffness constant defined as C′ ≡ (C 11 − C 12 )/2, it is actually found that C′ changes dramatically due to the magnetic ordering. C′ is related to the dynamical stability along the Bain path, 52 and therefore the softening of the 110 h i phonon branch is commonly considered as a precursor of the martensitic transformation. In experiments, temperature-induced martensitic transformations are indeed observed in disordered fcc Fe 1−x Pd x and Fe 1−x Pt x (x = 0.25-0.33) alloys, e.g., at 314.5 K 53 for Fe 0.72 Pd 0.28 and approximately 125 K 54 for Fe 0.73 Pt 0.27 . Having identified magnetic fluctuations as a driving force for the softening of the 110 h i phonon branch, our results suggest that the martensitic transformation in these Fe-based alloys is driven by magnon-phonon contributions.
Finally, we show the supercell-size dependence of the phonon spectra calculated using the band unfolding approach. We considered three supercell models with different sizes, namely the 32-atom 2 × 2 × 2, the 256-atom 4 × 4 × 4, and the 846-atom 6 × 6 × 6 supercells of the conventional fcc unit cell. The atomic configurations of the 2 × 2 × 2 supercell models were the same as those of the special-quasirandom-structure (SQS) 55 models used to extract the FCs, while the atomic configurations of the larger supercell models were determined using a pseudorandom number generator. Figure 4 shows the unfolded phonon spectra of the Fe 0.72 Pd 0.28 and Fe 0.72 Pt 0.28 in the ideal FM state obtained from these supercell models. Even for the smallest 2 × 2 × 2 supercell models (top panels), the phonon spectra are qualitatively similar to those obtained from the larger supercell models. The spectra obtained from the 2 × 2 × 2 supercell models are, however, less smooth than those obtained from the larger supercell models. The ruggedness in the spectra obtained from the 2 × 2 × 2 supercell models is due to the small number of local environments sampled around the atoms. In contrast, the spectra obtained from the 4 × 4 × 4 supercell models (middle panels) and from the 6 × 6 × 6 supercell models (bottom panels) are smooth and very similar to each other for both alloys. This indicates that the spectra obtained from the 6 × 6 × 6 supercell models are almost converged with respect to the supercell size.  16 At these temperatures both these disordered alloys are in the FM phase Temperature-dependent phonons of magnetic random alloys Y Ikeda et al.

DISCUSSION
We propose a first-principles-based method to incorporate both thermal magnetic fluctuations and chemical disorder into a unified computational framework of phonon calculations for magnetic random solid solutions. Chemical disorder, which leads to variations of atomic masses and FCs among atomic sites, is taken into account using the ICPA and the band unfolding in combination with SQSs. Thermal magnetic fluctuations are incorporated using the spin-space averaging and quantum Monte Carlo simulations using an effective Heisenberg spin Hamiltonian.
The proposed methodology is applied to Fe-Pd and Fe-Pt Invar alloys. Both the ICPA and the band unfolding are found to be equally capable for computing phonon spectra of the chemically disordered alloys. Both methods also reveal phonon broadening due to the chemical disorder. Taking thermal magnetic fluctuations into account, the developed approach shows excellent agreement between the peak positions of the computed temperature-dependent phonon spectra and experimental phonon frequencies. In particular, the approach correctly reproduces the experimentally observed unusual hardening in iron-based Invar alloys at high temperatures and shows that magnetic fluctuations are responsible for it. Finite-temperature magnon-phonon contributions also cause the experimentally observed softening of the elastic stiffness constant C′ and hence trigger the martensitic transformation in these alloys.
The proposed approach can be straightforwardly applied to other complex magnetic random solid solutions such as magnetic high-entropy alloys, which are expected to reveal complex physical mechanisms caused by the interplay between thermal magnetic fluctuations, chemical disorder, and lattice vibrations. It may be also possible to apply the proposed approach to magnetic random solid solutions with chemical short-range order; e.g., a supercell model of a disordered alloy with chemical short-range order can be used for the band unfolding.

METHODS
Chemical and magnetic (to obtain the FCs in the ideal PM state) disorder in these alloys were simulated by SQSs 55 based on the 32-atom 2 × 2 × 2 supercell of the conventional fcc unit cell (see Sec. A in the Supplemental Material). Chemical compositions for the SQS supercell models were chosen to be Fe 0.75 Pd 0.25 and Fe 0.75 Pt 0.25 , being close to the experimental ones, namely, Fe 0.72 Pd 0.28 and Fe 0.72 Pt 0.28 , respectively.
The mixing parameter α(T) in Eq. (1) was determined using an effective Heisenberg spin Hamiltonian for a 7 × 7 × 7 supercell including 1372 magnetic sites (see Sec. C in the Supplemental Material). The FCs of the SQS models were calculated using finite displacements of 0.01 Å. The FCs in the ideal PM state, Φ PM MÀM 0 , were obtained using the spin-space averaging method 21 by applying all inequivalent symmetry operations of the fcc structure to the SQS models with a disordered magnetic configuration (See also Sec. B in the Supplemental Material). Applying the symmetry operations to the SQS disordered-magnetic configuration, multiple different disordered magnetic configurations are obtained. At the magnetic high-temperature limit, such disordered magnetic configurations are dominant and have equal weights. 21,22 The averaging of the FCs was performed for each chemically inequivalent pair in the disordered-alloy model (e.g., Fe-Fe, Fe-Pd, Pd-Pd). For the chosen SQS model there are 1536 (=48 rotations × 32 translations) inequivalent symmetry operations. To be specific, taking the first nearest neighbor FCs of the Fe-Fe pair in the ideal PM state as an example, the statistical average with equal weights was computed over the FCs of 864 (=1536 × 0.75 × 0.75) atomic pairs in the SQS models. The FCs up to the fourth nearest neighbors were then employed in the ICPA and the band unfolding. For the band-unfolding procedure, the obtained Φ M−M′ (T) were assigned to the atomic sites of effective 864-atom 6 × 6 × 6 supercell models, where the composition ratios were fixed to Fe 0.75 Pd 0.25 and Fe 0.75 Pt 0.25 to be close to the experimental ones.
Electronic structures were calculated in the framework of densityfunctional theory within the generalized gradient approximation of the Perdew-Burke-Ernzerhof from 56 using the plane-wave basis projector augmented-wave method 57 as implemented in the vasp code. [58][59][60] A plane-wave energy cutoff of 350 eV was used. Since thermal expansion is small for the chosen Invar alloys, fixed lattice constants were applied for the computer simulations, which were taken from experimental data at room temperature; 3.755 Å 15 (disordered fcc Fe 0.72 Pd 0.28 ) and 3.749 Å 61 (disordered fcc Fe 0.72 Pt 0.28 ). Internal atomic positions were fully relaxed while keeping the lattices of the SQS models fixed.

Data availability
The authors declare that all the data supporting the findings of this study are available within the paper and its Supplemental Material.