Phonons in magnetically disordered materials Magnetic versus phononic time scales

The lattice dynamics in magnetic materials, such as Fe depends on the degree of disorder of the atomic magnetic moments and the time scale of spin ﬂuctuations. Using ﬁrst-principles methods, we have studied this effect by determining the force constant matrix in two limits: (i) When spin ﬂuctuations are much faster than the atom vibrations, their combined impact is captured by a spin-space averaged force constant matrix, (ii) when individual spin ﬂuctuations are sufﬁciently slow to scatter the phonon modes, the itinerant coherent potential approximation with spin-pair resolved force constants (i.e., (cid:2) ↑↑ , (cid:2) ↓↓ , and (cid:2) ↑↓ ) is employed in this paper. The physical consequences for the vibrational spectral functions are analyzed by systematically modifying the input parameters (magnetization and ratio of force constants betweens atoms with equal and opposite spin directions) and by deriving them for the prototype material system bcc Fe from ﬁrst-principles calculations. In the paramagnetic regime, the two limits yield identical phonon spectra. Below the Curie temperature, however, there are regions in the parametric study that show qualitative differences, including a broadening of the phonon peaks. For bcc Fe, however, the quantitative modiﬁcations of phonon frequencies turn out to be small.


I. INTRODUCTION
Lattice vibrations are typically the dominant temperatureinduced excitations in solids and, therefore, play a pivotal role in state-of-the-art materials science research [1]. Important insights about phase stability [2,3], ordering [4], and elastic properties [5] can be obtained from lattice dynamical studies. For chemically ordered systems, the calculation of quasiharmonic phonon frequencies requires the diagonalization of the dynamical matrix, which can be straightforwardly performed with standard techniques. The presence of disorder, however, introduces additional complexities. The disorderinduced scattering is determined by fluctuations in both the masses and the interatomic force constants [6][7][8][9][10], leading to off-diagonal disorder in the dynamical matrix. Moreover, the sum rule obeyed by the diagonal and off-diagonal parts of the force constants ensures that the disorder at a certain site depends upon its neighborhood-the so-called environmental disorder [11]. Due to this complex nature of the phonon problem, the calculation of the dynamical matrix in disordered systems is challenging.
The effect of large mass disorder between the components has meanwhile been investigated for several alloys [12][13][14].
Neutron-scattering measurements revealed disorder-induced asymmetric broad phonon peaks in these systems. The special features related to disorder further include the splitting of phonon modes into two branches [6,7]. If one studies the thermodynamics of pure Fe, however, one is not confronted with mass disorder due to the presence of only one chemical species. Nevertheless, there can still be disorder in the force constants between atoms with differently aligned magnetic moments. This situation is comparable to a Ni-Co alloy, which has negligible mass disorder but force constants that are significantly different from pure Ni. As a consequence of the disorder in the force constants, appreciable modifications in the phonon dispersions of a Ni-Co alloy as compared to bulk Ni has been observed in neutron-scattering experiments [15].
The theoretical investigation of vibrational effects in disordered systems was for a long time confined to chemically random alloys. Only recently, significant improvements in methodology and computational facility have provided access to systems with magnetic disorder as well [16][17][18][19][20]. The delicate interplay between magnetic and lattice degrees of freedom forms the basis for several investigations and interesting insights. For example, in CrN, the contribution of magnetic disorder to the vibrational free energy is essential for achieving transition temperatures from the cubic paramagnetic phase to the orthorhombic antiferromagnetic phase [21] in good agreement with experiment [22,23]. In Ni-Mn-Ga Heusler alloys, the approximate treatment of magnetic excitations by the fixed-spin moment approach removes the dynamical instability (soft phonon modes) in the cubic austenite thereby stabilizing this phase [24,25]. Furthermore, spin-phonon coupling is of tremendous importance for the interpretation of various phenomena in pure Fe and Fe-based materials [26,27]. Recently, it has been shown that the incorporation of a temperature-dependent magnetic energy yields reasonable agreement between theoretical phonons and experimental data for the temperature dependence of phonon frequencies in bcc Fe [28].
To meet these challenges, several theoretical approaches have been suggested in recent years to calculate phonons in paramagnetic systems with fully disordered spins: disordered local moments and spin molecular dynamics [29,30], a spin-spiral approach [31], dynamical mean-field theory [32], and a spin-space averaging (SSA) procedure [17]. All of them have specific approximations, particularly, connected to the handling of the magnetic vs phononic time scales. The SSA allows us to describe paramagnetic phonon modes and even their temperature dependence across the magnetic transition temperature reasonably well [28,33]. The approach has the underlying assumption that spin directions are changing so fast that the various local spin configurations become indistinguishable for the forces on the atoms. In other approaches, such as the above-mentioned molecular dynamics [29], the atomic forces depend on quasistatic magnetic configurations, allowing a modification of both on comparable time scales. It has recently been shown that a regular update of these magnetic configurations on the same time scale and dependent on the present displacement in the MD has a strong impact on phonon broadening [34].
The present study of magnetic disorder distinguishes force constants between pairs of atoms with equal and with opposite spins, whereas averaging force constants over magnetic configurations in the neighborhood of these pairs. The idea behind this approach is that a disordered arrangement of atoms with spin-up and spin-down changes sufficiently slow to resolve the interaction of lattice vibrations with these magnetic fluctuations. For this purpose, we make use of the coherent potential approximation (CPA) [35] and its generalizations. Similar to the application of CPA for chemical disorder, the traveling phonon experiences a nonhomogeneous background due to the fluctuations, which results into scattering effects and finite phonon lifetimes.
A straightforward application of CPA normally fails to provide an accurate description for phonons mainly due to two reasons: (i) It is a single-site theory, i.e., it cannot capture nonlocal fluctuations, and (ii) the resulting medium is structureless which prohibits incorporation of structural relaxations. Various generalizations of CPA have been suggested [36][37][38] of which two approaches based on the augmented space theorem of Mookerjee [39], have been particularly successful in the case of substitutionally disordered alloys. These are the itinerant coherent potential approximation (ICPA) of Ghosh et al. [40] and the augmented space recursion of Saha et al. [41] and Alam and Mookerjee [42]. Both of them have proven to provide almost identical results for several systems [43].
These methods have to date solely been applied to chemically disordered systems, and the generalization to magnetic disorder performed in this paper is so far missing. Beyond the relevance for magnetism, the study is also valuable for determining the impact of force constant fluctuations since no contributions of mass disorder are overshadowing the effect in this case. We demonstrate and discuss the salient features for the resulting phonon spectral density as a function of magnetization, including modified dispersions and the broadening of phonon linewidths. After having identified these trends, we apply our methodology to bcc Fe and clarify in this way the temperature-dependent impact of magnetic fluctuations on an intensively studied material system. The methodology, in the future, can be extended to study vibrational properties of magnetically disordered random solid solutions. Such an extension would require consideration of both chemical and magnetic disorders within the disorder model along with an expansion of the configurational space within the ICPA to specify both degrees of freedom.

II. METHODOLOGY
Our formalism consists of two major steps: First, it is necessary to derive the force constants of magnetically disordered Fe using density functional theory (DFT). Only if they are available, ICPA can be employed in a second step to perform the configuration averaging and to consider the magnetic fluctuations.

A. Density functional theory
The DFT calculations for Fe have been performed with a plane-wave basis set as implemented in the Vienna ab initio simulation package (VASP) [44,45] in order to obtain accurate forces based on the Hellmann-Feynman theorem. The ion-electron interactions are treated with the projectoraugmented wave method [46] together with the generalized gradient approximation for the exchange-correlation potential parametrized by Perdew-Burke-Ernzerhof [47]. A large basis with a plane-wave cutoff of 340 eV and a k-point grid of 7 × 7 × 7 are used. For the force calculations, the Methfessel-Paxton scheme [48] with a smearing of 0.1 eV is used. The ab initio force constants are calculated using a direct approach in which each atom in a 4 × 4 × 4 supercell (of the primitive unit cell) is moved by 0.015 Å along three Cartesian axes, and the forces on the atoms are calculated. The lattice parameter used for bcc Fe is 2.871 Å. The ICPA calculations described below are performed with a 25 × 25 × 25 k mesh and 2000 energy points.

B. Magnetic disorder
For the calculation of the force constant matrix in the paramagnetic state, we use a randomly disordered collinear spin configuration constructed using the concept of special quasirandom structures (SQSs) [49] as obtained from the ATAT package [50]. Within these structures, the N atoms in a given periodically repeated cell are distributed such that their distinct correlation functions (k,m) best match the ensembleaveraged correlation functions (k,m) for the infinite perfectly random configuration. Here, (k, m) corresponds to the vertex (structural figure) defined by the number of involved atoms k (pairs, triples, ...) and the order of separation of the atoms m [first, second, ..., nearest neighbors (nns)].
Because of the low symmetry of the SQSs, the force constant tensors between a given pair of atoms do not have the symmetry of the underlying lattice (here: bcc). Therefore, the displacement of each atom gives rise to a different force constant matrix. Inelastic neutron scattering, on the other hand, uses the crystal symmetry of the ideal lattice to 094201-2 determine phonon dispersions. In order to achieve force constant tensors with such a symmetry also for the disordered alloy, we perform a two-step averaging procedure [51].
(1) We average for each atomic pair individually (here: up to second-nearest neighbor) over the force constants, to achieve transformed 3 × 3 force constant matrices with the rotational symmetry of the underlying crystal lattice. This step is identical to SSA where it is motivated by a much shorter magnetic as compared to the phononic time scale.
(2) We average over the force constant matrices of all atomic pairs in the supercell belonging to one of the sets ↑↑, ↓↓, or ↑↓ to eventually achieve three 3 × 3 force constant matrices ↑↑ , ↓↓ , and ↑↓ that also obey the translational symmetry of the crystal. The distinction between spin pairs allows for a partial consideration of magnetic disorder and their fluctuations on a time scale comparable to phonons.
Apart from the fully paramagnetic state, also unequal concentrations of ↑ and ↓ spins as expected below the Curie temperature are considered. Furthermore, to mimic magnetic short-range order (SRO), the overall concentrations are kept equal, but the second averaging is performed only over force constants that have been determined in a configuration with an unequal number of ↑ and ↓ spins for the nearest-neighbor shell.

C. Itinerant coherent potential approximation (ICPA)
The central quantity for the determination of phonon energies and linewidths in magnetically disordered Fe is the spectral density belonging to the Green's function of lattice excitations, which itself is the Fourier transform of a displacementdisplacement Green's function. Here, is again the 3 × 3 force constant matrix as determined above, ω is the frequency, containing a small imaginary part, and m is the mass operator. In the case of paramagnetic Fe, the two kinds of atoms ↑ and ↓ have identical mass, i.e., m = m1. This simplifies the formalism as compared to a situation with chemical disorder [40]. However, since the force constants ↑↑ , ↓↓ , and ↑↓ still represent a disorder, the averaging · · · is performed over spin-configurational degrees of freedom.
The disorder is characterized by the concentrations c ↑ and c ↓ of atoms with collinear spins in both directions. Within the space of all possible spin configurations of the system, the state of site i is defined by |↑ i or |↓ i . The Hamiltonian of the system, however, operates in the Hilbert space where mass m, force constant , and Green's function G are defined. Using the augmented space formalism (ASF) [39], we expand the Hilbert space with the configuration space to take into account the statistics of random site occupancy. To this end, one of the bases representing the site-average state in the configuration space, which is augmented to the real space, is given as If we would only consider the projection¯ of the force constant operatorˆ on this subspace, then the corresponding Green's function was that of the SSA approach, This is because the corresponding matrix element in Eq. (A5) completes the averaging over force constants such that no distinction between spin pairs is present anymore.
Within the ICPA method [40], corrections are employed before the averaging over spin directions happens. Based on the profound experience with chemical disorder, the singlefluctuation approximation is used as Thus, | j f k is a state in which fluctuation in the average state |0 is present only on site j. The states in the ASF can, therefore, be expressed as |i0 and |i f k where i specifies the site index of the dynamic variable, position or momentum, and k indicates the presence of fluctuation on site i. Within the ASF, the force constant operator becomes in a block representation,ˆ with block matrix elements, where i and j specify site indices whereas k and l indicate the presence of fluctuation on these sites. The ASF is a smart bookkeeping technique to keep track of the fluctuations in the environment due to disorder. For details of the performed projection, the reader is referred to Ref. [40]. A similar ASF representation could be also introduced for the Green's function. The relevant propagator is only the configurationaverage part, and due to the inversion in (5), the self-energy (ω 2 ) contains fluctuations of the form with where1 is the identity matrix associated with two sites with fluctuation present on both sites. In particular, and † in (8) can be considered as annihilation and creation operators of single spin fluctuations in the material.
We follow the philosophy of CPA and replace the Dyson equation for the itinerator, by its self-consistent extension, where G (i) is identical to the full propagatorḠ in (7) except that all irreducible scatterings beginning or ending on site i are omitted. These equations need to be resolved iteratively. The concept corresponds to the traveling cluster approximation (TCA) [52] due to the "itineration" of fluctuations caused by the off-diagonal character of the force constant matrix involving, e.g., nearest-neighbor interactions. Note that other cluster generalizations of CPA are typically difficult to implement beyond pair scattering and diagonal disorder since they do not provide any precise guidance regarding the diagrams that are to be summed. Due to the mapping onto matrix elements in the augmented space [53] as performed here, the enlarged Hilbert space to accommodate configuration fluctuations due to spin disorder, and the overlapping sets with scattering are straightforwardly mapped onto the set of sites with "fluctuation states" [54]. Thus, the ideas in the TCA are implemented with the advantage of a tractable approach to treat off-diagonal disorder.
Since the Green's function in the augmented space is, by construction, site translation invariant, it can also be Fourier transformed to q space as where N defines the number of neighbors perturbed due to the single-site fluctuation in a lattice. For the investigation of phonon spectra throughout the Brillouin zone, we are interested in the peak position and width of the spectral function of this Green's function, with G σ 1 ,σ 2 ( q, ω 2 ) being the conditional or partial Green's function in q space. The advantage of the method is that these results are obtained spin resolved with σ 1 , σ 2 representing the spin directions. The required matrix elements in the special case of a bcc crystal structure are provided in the Appendix. Particularly important is , which directly determines the self-energy in (8) and is of the principal shape, Using these equations, one can quickly convince oneself that some special cases are fulfilled. First of all, the impact of fluctuations disappears by definition, if ↑↑ , ↓↓ , and ↑↓ are identical, i.e., if the forces do not depend on the spin. The other way around, fluctuation effects become more dominant the larger the difference ↑↑ − ↑↓ of the system is since Less obvious is the fact that the fluctuation effects also disappear in the ideal paramagnetic limit, where c ↑ = c ↓ and ↑↑ = ↓↓ . This can be considered as a first major result of the ICPA approach. It represents the fact that single-site fluctuations of spins do not change the ideal paramagnetic state for the force constants. This is schematically demonstrated in Fig. 1. In this respect, the impact of the force-constant disorder is qualitatively different from mass disorder. The situation can be also different for two-site spin fluctuations. In our evaluation of this limit not only vanishes, but also Eq. (9) becomes identical with the SSA propagator G SSA , and the spectral density will have sharp peaks without any broadening.
The self-energy (8) obviously also disappears in the case of magnetic saturation c ↑ = 1 and c ↓ = 0. The fluctuations are, therefore, expected to be strongest close to the magnetic transition temperature where a long-range or a short-range order is present.

A. Force constant matrix
We take paramagnetic bcc Fe described by a SQS model as a starting point for our discussion to get realistic input parameters for the force constant matrix. For this purpose, we have performed collinear DFT calculations (Sec. II A) and have averaged the force constants (Sec. II B). Tests with 094201-4   Table I contains the results for ↑↑ and ↑↓ for the paramagnet (c ↑ = c ↓ ). In this case, the symmetry requires ↑↑ = ↓↓ . We observe that the ratio in the nn force constants between ↑↑ and ↑↓ is for paramagnetic bcc Fe below 1.3. This is even true if the averaging is limited to those local configurations where the number of nearest-neighbor ↑ and ↓ spins is not equal, labeled as SRO in Table I (second  block).
It is also evident from Table I that, in bcc Fe, the nnn force constants are-in contrast to fcc Fe [40]-still significant. This is a consequence of the more open structure of bcc with a comparable distance of the first and second atomic shells. We further note the strong impact of the magnetic configuration on these force constants. Therefore, the nnn shell has been taken into account. In contrast to nn and nnn force constants, third-and fourth-neighbor force constants are orders of magnitude smaller and, hence, have negligible impact on phonon frequencies.
In the case that the magnetic moments are partially ordered-such as c ↑ :c ↓ = 3:1 (see the third block of Table I)-already the ratio between ↑↑ and ↓↓ in the nn shell can be as high as 2.0, whereas even the sign of the force constants in the nnn shell can be different. Furthermore, ↑↑ and ↓↓ do not need to be identical anymore. These observations raise the expectation that significant deviations from a SSA approach can be expected below the Curie temperature.

B. Model parameters
To qualitatively investigate the impact of the resulting fluctuations on the spectral density, we have simplified the DFT values for the force constants of bcc Fe to model parameters and have systematically modified them. Using Table I as a guideline for realistic sizes and ratios of the force constants, the actually chosen model parameters are given in Table II Fig. 2 is close to the homogeneous case of ↑↑ = ↓↓ = ↑↓ for which no fluctuations are present, and, therefore, a δ peak in the spectral function is expected. The small broadening of the peak is due to the ratio nn ↑↑ / nn ↓↓ = 2.0, which has, however, little impact since only 1/4 2 = 6% of the nn force constants are described by nn ↓↓ . In contrast to this, the ratio nn ↑↑ / nn ↑↓ has a stronger impact since nn ↑↓ describes 6/16 = 37.5% of the nn force constants. The modifications have two consequences for the spectral function: On one hand, the position of the phonon peak is shifted. This seems to be mainly caused by the increase in nn ↑↑ as indicated by the partial spectral function for the ↑↑ contribution to the Green's function (the red line in Fig. 2).

The increase in nn
↑↑ is chosen such that the SSA result (the vertical dashed line) is not changed by the modification of the force constants. Nevertheless, the difference of ICPA and SSA phonon energies seems to scale with the ratio nn ↑↑ / nn ↑↓ . On the other hand, with the increase in the ratio nn ↑↑ / nn ↑↓ , the spectral lines become broad and asymmetric. In particular, a strong shoulder towards lower frequencies develops. In cases of extreme disorder, the spectral function can even split into two peaks. The main reason seems to be that the peak of the partial spectral function for the ↑↓ contribution (the green line in Fig. 2) does not follow the trend of the ↑↑ partial spectral function and remains close to the SSA value. It is additionally broadened towards lower frequencies and even shows negative values. The emergence of negative partial spectral densities is mathematically correct since only the total spectral function has a physical meaning. In fact, the negative contributions are mainly compensated by the enhanced values of the ↓↓ partial spectral function [8,9].
The behavior demonstrates the averaging of spin degrees of freedom in the ICPA: The initial assumption is that the force constants of the majority and minority as well as the equal and opposite spin can be distinguished. The SSA part (3) of the TABLE II. Real-space nearest-neighbor force constants αβ (notation as in Table I) chosen for the systematic study in Fig. 2. The units are 10 3 dyn cm −1 . Green's function (7) Table I have a very limited impact on the ICPA spectral function. This is mainly due to the fact that the symmetry constraints of the bcc lattice yield only diagonal contributions in nnn [see Eq. (A3)] and, consequently, limit the number of terms entering¯ ,˜ , and .
After having studied the impact of the force constants on the phonons, we next investigate the impact of the magnetic temperature by changing the concentrations c ↑ :c ↓ in the ICPA (which was fixed to 3:1 in Fig. 2). In Fig. 3, the dependence along the [ζ , ζ , 0]L phonon branch for the fully disordered (c ↑ = 0.50), a partially ordered (c ↑ = 0.75) and a fully ordered (c ↑ = 1.0) magnetic state is provided, whereas keeping the interatomic force constants unchanged. Taking the ferromagnetic state as a reference, the phonons are significantly softened when decreasing the magnetization (i.e., the concentration c ↑ ). This tendency is clearly visible throughout the whole phonon branch but is particularly pronounced in the region 0.2 ζ 0.4.
In order to understand the wave-vector dependence along [ζ , ζ , 0]L, we have plotted the partial spectral functions of c ↑ = 0.5 for three ζ values in the insets of Fig. 3. It can be seen that, at ζ = 0.15, all three contributions to the total spectral function are peaked at almost the same frequency, which is identical to the SSA result (the dashed line). At ζ = 0.3, a separation of the ↑↑ contribution from the ↓↓ and ↑↓ contributions can be noted because the stronger force constant gives rise to higher frequencies. Since the latter contributions (the blue dashed dispersion line) still dominate Frequency (meV) Total spectral function (arbitrary unit) the total spectral density, the phonon frequency is clearly shifted below the SSA result. At the same time, we observe a broadening of the phonon peak. As the wave vector is increased, the respective peaks further separate, and the distribution of spectral weight of the ↑↑ contribution (the red dashed dispersion line) is changed such that the twopeak structure changes to a pronounced peak well above the SSA position. This high-frequency peak dominates the total spectral function close to the edge of the Brillouin zone, therewith reducing the softening effect as compared to SSA.
The dependence of the full spectral function on the concentration ratio c ↑ :c ↓ is shown in the lower panel of The concentration ratios c ↑ = c ↓ (ideal paramagnet) and c ↑ :c ↓ = 3:1 (ferromagnet close to saturation) are considered. The SSA approach with magnon-phonon coupling is an interpolation to the experimental temperature (1043 K) [28]. The "6up2dn" result corresponds to force constants in a paramagnetic supercell with SRO such that locally c ↑ :c ↓ = 3:1. The black filled circles correspond to experimental data [55].
(c ↑ = 1.0). The opposite limit of c ↑ = c ↓ also yields to a sharpening of peaks due to reduced fluctuation effects. However, it does not correspond to the perfect paramagnet because the force constants are fixed to nn ↑↑ / nn ↓↓ = 2. At ζ = 0.45, remainders of the two-peak structure for c ↑ = 0.5 are visible, but the peaks are rather flat. The behavior for intermediate concentrations is for all ζ values a smooth transition between these two extreme cases. The distance between the extrema depends again on the wave vector but seems to be largest when the spectral function becomes broadest, i.e., when the strongest fluctuation effects occur.

C. Application to bcc Fe
To quantitatively evaluate the performance of the ICPA for bcc Fe, we now return to the force constants presented in Table I. All three scenarios considered in the table are plotted in Fig. 4. As indicated earlier, the ICPA result for c ↑ = c ↓ (the brown lines) is formally identical with the SSA result published in Ref. [17], except the restriction of the force constants to first and second neighbors. Therefore, no fluctuation effects are expected in this high-temperature limit.
Taking into account that the experiments have been performed at the Curie temperature of Fe (1043 K), a SRO of localized moments instead of a fully random distribution can be expected [56]. To evaluate its impact, we have selected only those nearest-neighbor configurations with six ↑ and two ↓ spins. The resulting force constants (the second block of Table I) are approximately the same as for all the other configurations (the first block of Table I). Nevertheless, due to the different concentrations c ↑ and c ↓ used in this case (considering the SRO as a total magnetization), the corresponding phonon energies are significantly shifted upwards (the red line in Fig. 4). Although we can now also formally expect deviations of the ICPA phonon energies from a similar SSA procedure, the relatively small difference between φ ↑↑ and φ ↑↓ will make the quantitative phonon corrections also small. We, therefore, compare the approach with the interpolation scheme of SSA and ferromagnetic force constants to describe the spin behavior at 1043 K as published in Ref. [28] (the orange line in Fig. 4). We note that the soft transversal branches are better captured by the SRO assumption, whereas the interpolation schemes works better for the stiffer modes.
A direct comparison of ICPA and SSA phonon energies can be performed if c ↑ :c ↓ = 3:1 globally throughout the crystal. As compared to the scenario with only SRO configurations of c ↑ :c ↓ = 3:1, the phonon energies are again shifted upwards. However, the force constant disorder caused by the difference between φ ↑↑ and φ ↑↓ (the third block of Table I) is in this case even smaller than in the previous discussion. As a consequence, the fluctuation effects on the phonon energies (compare the blue and green lines in Fig. 4) are below the experimental resolution. Hence, bcc Fe seems not to be a material for which a slower time scale for magnetic fluctuations as considered in ICPA has a noticeable impact on phonon energies. The resulting phonon scattering becomes, however, apparent if one studies phonon linewidths. As discussed in the previous subsection, they are strongly influenced by force constant disorder. Figure 5 demonstrates that the relatively small spin-induced force constant disorder of bcc Fe has a noticeable impact on phonon lifetimes. The broadening of the phonon spectrum seems to be strongest for transversal branches close to high-symmetry points in the Brillouin zone.

IV. CONCLUSIONS
In conclusion, we have investigated the impact of time scales on the lattice dynamics in systems with magnetic fluctuations. To this end, the ICPA has been employed to describe the magnetic disorder. In contrast to chemically disordered systems for which ICPA has been originally developed, no mass disorder is present in this case. The method averages over different local spin configurations contained in a fully disordered spin distribution for a supercell in such a way 094201-7 that a temporal distinction between force constants belonging to the three different spin pairs ↑↑, ↓↓, and ↑↓ of the neighboring atoms remains possible. This is different from the SSA approach where a time scale of spin fluctuations is considered that is too fast to allow for such a distinction.
In the paramagnetic case, characterized by equal concentrations of up and down spins, both approaches yield identical results. This indicates that, instead of a consideration of individual fluctuations, an average over all spin configuration as performed in SSA is sufficient in this temperature regime. Below the Curie temperature, however, the relevance of magnetic fluctuations for phonon line shapes and frequencies has been demonstrated.
Nevertheless, when calculating the phonon dispersion in the prototype material system bcc Fe, the results of ICPA and SSA are very similar even below the Curie temperature. The main difference is the observation of a phonon broadening, which is predicted by ICPA but does not exist on the time scale considered by SSA. In addition, the impact of SRO in the paramagnetic case has been studied. The comparison of the obtained results with the experimental data indicates the validity of our approach. These findings, thus, suggest that magnetic disorder as present near the Curie temperature can influence the lattice dynamics in magnetic systems.

ACKNOWLEDGMENT
Funding by the Deutsche Forschungsgemeinschaft (DFG) within the priority Programme No. SPP 1599 is gratefully acknowledged.

APPENDIX: MATRIX ELEMENTS
In the following, the expressions for the operators G −1 SSA , , and (˜ −¯ ·1), which are required as input for the propagator (7), are provided for a bcc crystal structure. Similarly, the expressions for fcc crystals, given in the Appendix of Ref. [40], can be reduced to the case of vanishing mass disorder.
The nn and nnn sites are (± 1 2 , ± 1 2 , ± 1 2 ), (±1, 0, 0), (0, ±1, 0), and (0, 0, ±1), respectively, and corresponding force constant matrices are 00 = (000,000 The force constants between the reference atom and the other neighbors can be calculated by applying cubic symmetry operations to the above two matrices. We use the notation for the Cartesian com-ponents α, β of relative lattice vectors R α/β i and R α/β j for the two neighboring atoms at sites i and j. We then obtain the following matrix elements for¯ i j (with α = β): where R nn and R nnn specify distance between sites at nn and nnn positions, respectively. The matrix elements for i j are only nonzero if either i or j is a fluctuation site. In this case, we obtain, similar to (A4), the following matrix elements (α = β): where the last three matrix elements are given for the special case of φ ↑↑ = φ ↓↓ . In the case of˜ i j ( q), a distinction between the positions of the fluctuation sites is required. If i and j are fluctuation sites (corresponding to creation and annihilation sites) at nn or nnn positions, then the following equations hold: ↑↓ ]e i q· R nnn for R α nnn = ±1, (˜ ( q)) αα nnn = [2c ↑ c ↓ φ 2yy ↑↑ + (c 2 ↑ + c 2 ↓ )φ 2yy ↑↓ ] − 2c ↑ c ↓ [φ 2yy ↑↑ − φ 2yy ↑↓ ]e i q· R nnn for R α nnn = 0, with φ ↑↑ = φ ↓↓ for the last three terms. If only i or j is a fluctuation site, then it is required that another site k is involved, at which the second fluctuation takes place. If k is located at a 094201-8 nn position of i and j such that R nn = R ki = − R k j , then, If k is located at a nnn position of i and j such that R nnn = R ki = − R k j , then, ↑↓ ]e i q· R nnn for R α nnn = ±1, ↑↓ ]e i q· R nnn for R α nnn = 0, (A8) with φ ↑↑ = φ ↓↓ for the last three terms. If none of the above conditions are fulfilled, then˜ i j ( q) becomes a null matrix. If φ ↑↑ = φ ↓↓ = φ ↑↓ , i.e., the dynamical matrix does not depend on the spin configuration, then, the system is free from any kind of fluctuation. This is consistent with the mathematical observation that and (˜ −¯ )( q) reduce to null matrices. Subsequently, the self-energy contribution due to scattering vanishes. This even happens if φ ↑↑ = φ ↑↓ , but φ ↑↑ = φ ↓↓ and c ↑ = c ↓ in the present case. Under such circumstances, the ICPA and the SSA yield identical results.