Vibrational and magnetic signatures of extended defects in Fe

Defects change the phonon spectrum and also the magnetic properties of bcc-Fe. Using molecular dynamics simulation, the influence of defects – vacancies, dislocations, and grain boundaries – on the phonon spectra and magnetic properties of bcc-Fe is determined. It is found that the main influence of defects consists in a decrease of the amplitude of the longitudinal peak, PL, at around 37 meV. While the change in phonon spectra shows only little dependence on the defect type, the quantitative decrease of PL is proportional to the defect concentration. Local magnetic moments can be determined from the local atomic volumes. Again, the changes in the magnetic moments of a defective crystal are linear in the defect concentrations. In addition, the change of the phonon density of states and the magnetic moments under homogeneous uniaxial strain are investigated.


Introduction
Since point defects are the most common defect type in metals, their influence on the phonon density of states (pDOS) was studied intensely. This influence is of basic importance because it enters the calculation of defect migration energies and defect formation entropies as well as the temperature dependence of the free energy of Fe. Early calculations of the vibrational spectra of vacancycontaining crystals were limited to pair potentials [1][2][3][4][5]. In view of the known limitations of pair potentials to describe the elastic properties of metals properly [6], later approaches turned to many-body potentials. Kislov and Mazurenko [7] calculated the vibrational spectrum of bcc-Fe using the Finnis-Sinclair potential [8]. They found a shift of the transverse phonon peaks towards lower frequencies which they ascribed to a weakening of the interatomic interaction in the near-neighbor environment. Gairola et al. [9] calculated the local density of states of the neighbors of a vacancy in bcc Fe using a Green's function method and a self-determined interaction potential. They found that changes in the vibrational spectrum are restricted to first neighbors; the pDOS features a decrease towards higher frequencies and a small shift towards lower frequencies. Self-interstitials show localized modes (related to the stretching of the dumbbell bond), which are absent for vacancies; their vibrational properties a e-mail: urbassek@rhrk.uni-kl.de have been calculated more rarely [10,11]. Also localized low-frequency modes exist that correspond to easy glide directions of the interstitial [11].
In contrast to point defects, the effect of extended defects on the vibrational spectrum has been determined only rarely. An exception is provided by interstitial clusters, since these are important in radiation damage studies [11].
Derlet et al. [12] studied the low-frequency behavior of nano-crystalline Cu. They found that in the lowfrequency limit, the grain-boundary component of the pDOS exhibits an increase ∝ ω rather than ∝ ω 2 . They also found a lowering of the longitudinal peak in nanocrystalline matter.
Later Meyer et al. [13] studied how the pDOS of nanoclusters evolves towards that of the bulk material. They found that the contribution of atoms that are not in a well-characterized crystal environment shows a shift towards lower frequencies and a lowered longitudinal peak amplitude.
In the present paper, we use atomistic simulation to study the effect of defects on the pDOS and the magnetic moments in bcc Fe. Isolated vacancies are studied for comparison purposes, while the study of extended dislocation networks is novel. Poly-crystalline specimens are also investigated due to their principal importance for experiment.

Methods
Our Fe crystals are cubic with a side length of 14.4 nm and contain N = 250 000 atoms. We equilibrate them at a temperature of 300 K for 50 ps, before evaluating the pDOS.
Defects are introduced into the sample as follows.
(i) Vacancies are introduced by randomly deleting atoms in the crystal. Vacancy concentrations up to n v = 0.03 are considered. (ii) Dislocations are produced according to the recipe of [28,29] by introducing straight edge dislocations, running in 100 direction, into the crystal. The dislocations were annealed according to the recipe of [30] by heating to 1230 K, equilibrating for 100 ps and then cooling down to 300 K. We produce 9 samples with dislocation densities ρ between 0 and 5.959 × 10 16 m −2 . The resulting structures contain a dislocation network consisting of dislocations with Burgers vectors b = 1 2 111 and b = 100 , in which the latter are remnants from the production process. Here, the dislocation density is calculated as ρ = L disl /V , where L disl is the total length of dislocations in the sample, and V is the sample volume. For purposes of comparison, we can relate the dislocation density to a defect concentration n c as follows.
If Ω 0 = 11.78Å 3 denotes the atomic volume of bcc Fe, the number of atoms in the core of a dislocation line of length L disl is given by (1) (iii) As an alternative method to introduce dislocations into the sample, we strained the sample uniaxially along a 100 direction with a strain rate of 10 10 s −1 .
In the orthogonal directions, the boundary conditions guaranteed zero stress. Before evaluating the data, we equilibrated the system for 10 ps. (iv) Grain boundaries are generated using the free software 'atomsk' [31]. The nanocrystalline specimen created contains 6 grains with an average size of 9.75 nm. They are equilibrated by annealing at 1230 K, as described above for dislocated crystals, item (ii) above.
The pDOS, γ(ω), is calculated from the velocity autocorrelation function γ(τ ) of the iron atoms Fig. 1. pDOS as calculated using the Ackland et al. [33] potential, compared to experimental data at 295 K [55]. The quadratic dependence at small phonon energies is highlighted.
where the sum runs over all iron atoms i in the sample and . . . t0 denotes the average taken over the reference times t 0 . All calculations are performed at 300 K, unless noted otherwise. The velocity autocorrelation function is calculated for times up to τ max = 50 ps; this appears sufficient since it corresponds to a frequency resolution of 1/τ max = 0.02 THz or 0.67 cm −1 [32]. The pDOS is obtained as the Fourier transform of the velocity autocorrelation function [32]. We normalize it to area 1 such that γ(ω) = 1. We use the Fe potential developed by Ackland et al. [33] in this study, since it provides a fair representation of the Fe pDOS when comparing to experimental data. A variety of other potentials were investigated by us, but they did not reproduce this quantity as satisfactorily, see Appendix A.
The software LAMMPS [34] is used for performing the simulations. OVITO [35] is employed to render snapshots and to calculate the dislocation length.

pDOS
Experimental data of the Fe pDOS have been available since the late 1960s where the phonon dispersion was measured by neutron scattering [36][37][38]. We reproduce these data in Figure 1 for 300 K and compare them with our calculated data. Deviations between experiment and simulation show up; as discussed in the Appendix A, these are typical for all interatomic potentials.
The shape of the pDOS, g(ω), is characterized by a quadratic increase at low frequencies according to where v D is the Debye velocity of sound, and a highfrequency cut-off. In between, the shape is dominated by 3 peaks; the two lower-frequency peaks correspond to mainly transverse-polarized phonons, while the last peak is due to longitudinal phonons.
In the present study, we base our analysis on the peak height of the longitudinal peak, since this is a quantity that is easily obtained also from an experimentally measured pDOS. Certainly, from a theoretical point of view, it might be preferable to base the analysis on the fraction of modes present under the entire longitudinal peak by integrating the pDOS around a suitable interval around the peak maximum. However, this interval might be difficult to define as it changes with the defect types and concentrations considered. We therefore restrict the discussion to the peak height itself.
At low temperatures, the spectrum is shifted to higher frequencies and the relative position of the peaks changes. For this reason, we restrict our analysis to one temperature, 300 K.

Defective crystals: vacancies
The influence of vacancies on the vibrational spectrum of crystals has been investigated since long. As outlined in the Introduction, Section 1, the effect of vacancies is localized to the first two neighbor shells around it and can be characterized as a shift of the spectrum to smaller frequencies.
In Figure 2a, we investigate the effect on the pDOS of a random assembly of vacancies in a bcc-Fe crystal. Clearly, the vacancy densities applied in our study are considerably larger than in experiment; they are used in order to bring out the influences more strikingly. The main effect is a decrease of the height of the longitudinal peak situated at around 37 meV. To a lesser extent, also the second transverse peak, at around 27 meV, is reduced. In compensation, the minima between the peaks are slightly filled up, and also the pDOS at low frequencies, between 8 and 18 meV, is slightly increased. A change of the low-frequency dependence, g(ω) ∝ v 3 D ω 2 , equation (3), can hardly be detected below 5 meV. This independence implies that the elastic moduli of Fe are only negligibly changed by the introduction of vacancies at the densities applied here.
We conclude that the change in the longitudinal phonon peak is the most pronounced feature of the influence of vacancies on the phonon spectrum. Figure 2b shows that up to a vacancy concentration of 2%, the longitudinal peak amplitude decreases linearly with the vacancy density as where P L denotes the longitudinal peak height relative to that of an ideal crystal, and n v is the vacancy concentration. We find α v = 5.083. For higher densities, the peak amplitude appears to stabilize.

Defective crystals: dislocations
The change of the pDOS upon the insertion of dislocations is shown in Figure 3a. The effects parallel closely those in a vacancy-filled crystal; the main effect is again a decrease in the longitudinal peak amplitude. As shown in Figure 3b, the effect is again approximately linear in the defect density, in particular for dislocation densities ρ 3 × 10 16 m −2 . Here, we find in analogy to equation (4) where α d = 0.03 (10 16 m −2 ) −1 at low ρ 3 × 10 16 m −2 and α d = 0.013 (10 16 m −2 ) −1 at higher ρ.
In order to compare the changes induced by dislocations with those induced by vacancies, we convert the dislocation densities to the concentration of core atoms, n c , via equation (1). We note that this comparison only allows to assess the strength of the respective influence of vacancies and dislocations; the physical origin of the changes in the pDOS are caused by the structure of the defects -in particular the local environment around the defects -and will not be assessed here. In terms of the concentration of core atoms, n c , the change in peak height reads with . Hence α c = 60 at low ρ 3 × 10 16 m −2 and α c = 25 at higher ρ. These values are a factor of 5-12 higher than the vacancy values α v . This discrepancy can be interpreted such that the concentration of core atoms is larger than the estimate equation (1) by a factor of around 10. In other words, the cross-sectional area around the dislocation line where the density is changed with respect to the ideal-lattice density is larger than an atomic cross section. This is in line with atomistic calculations of the core structure of dislocations in Fe [17].
The effect of a single dislocation on the frequency spectrum of a crystal does not appear to have been studied previously in such detail as that of point defects. It may however be assumed that the effect is based on two aspects: (i) the underdense dislocation core, (ii) the farreaching strain field around a dislocation. Atoms in and around the dislocation core are expected to experience similar effects as around vacancies, since also here the local atom density is reduced. The strain field around dislocations, on the other hand, is complex and has both tensile and compressive components [39,40]; it is therefore hard to predict its influence on the height of the longitudinal peak, P L . If one assumes that the close environment of the dislocation dominates the effect on P L , the similarity of the vibrational spectra of vacancy-and dislocation-filled crystals in decreasing the longitudinal-peak height becomes understandable. Note, however, that dislocations also populate high-frequency modes -in contrast to vacanciesand hence lead to a broadening of the longitudinal peak.

Defective crystals: strain
An alternative method to introduce dislocations into a crystal consists in straining the crystal until it responds by dislocation formation. Figure 4a shows that dislocations start to develop at a strain of around 18%; this late start of defect formation is caused by our use of an ideal crystal. Then the dislocation density increases more or less monotonically until a strain of around 40%, and then saturates at values of around 10 17 m −2 . Clearly, such high strains can only be applied to ideal crystals without fracturing them. Figure 4b shows that under strain the spectra change more pronouncedly than in the cases considered above. In particular, a clear increase in the curvature at low frequencies is seen. In agreement with equation (3), this increase is related to the increase of elastic constants in the strained material. Concomitantly, a slight increase of the pDOS at the highest frequencies is observed; in particular, also the frequency of the highest phonon mode is lifted. Again, this is caused by the high strain in the crystal.
Similarly as for the defects described previously -vacancies and unstrained dislocations -the amplitude of the longitudinal peak decreases. However, the changes in the minima and also in the second transverse peak are now more pronounced than in the previous cases. Figures 4c and 4d quantify the decrease of the longitudinal peak amplitude with increasing strain and dislocation density, respectively. The dependence now differs strongly from the approximately linear dependences found above.
(i) Up to a tensile strain of 5%, the pDOS and P L are unchanged. (ii) At 10% strain, P L is significantly reduced, by around 25%, even though no dislocations were generated. Figure 4b reveals that at this strain, the spectrum already considerably changes. Phonon modes are re-distributed towards the low-and high-frequency regions, signaling the strain-induced increase of the elastic moduli and the speed of sound on the one hand and the increase of the maximum phonon frequency on the other hand. As a result, the entire phonon peak region between 20 and 40 meV is depopulated of modes. Hence the change occurring between 5 and 10% strain is not caused by defect formation, but by a re-distribution of modes from the peak region to the edges of the spectrum. (iii) At even higher strains, the relative changes are minor and mostly produced by noise. (iv) However, for strains between roughly 20 and 30%, where the dislocation density increases strongly, an approximately linear decrease of P L with ρ can be found, similar to equation (5): where α s = 0.013 (10 16 m −2 ) −1 , in excellent agreement with the value of α d given in Section 3.2 above. (v) In this same regime, the dependence of P L on strain is with α p = 0.244.

Defective crystals: grain boundaries
In Figure 5, we compare the pDOS of a single-crystalline sample with that of the nanocrystalline sample. Again we observe a relocation of modes from the central peak region of the spectrum towards the high-and low-frequency tails. This feature is in agreement with previous experimental [41] and simulational [12] studies of nanocrystalline Ni which showed that with decreasing grain size, this mode relocation becomes more pronounced. The increase of both low-and high-frequency modes is attributed to modes localized at grain boundaries; this redistribution leads to a 'broadening' of the pDOS towards low and high frequencies [13]. When straining the poly-crystalline sample, the strain induces complex changes both in the grain-boundary regions and in the bulk of the grains; such changes may occur already at considerably lower strains than in a single-crystalline specimen [42]. However, an analysis of defect generation in poly-crystalline materials is nontrivial, since grain boundaries are easily mis-identified as dislocations by defect-identification algorithms such as OVITO; in addition, it is hard to monitor changes in the grain boundaries during straining. We therefore refrain from correlating strain with defect densities in the polycrystal.  Figure 6 shows that -in contrast to the single-crystal, Figure 4 -the change in P L depends rather monotonously on strain. A roughly linear relationship can be read off Figure 6b, which we write as with α p = 0.43. Here we excluded the initial value at zero strain from our fit, since the initial decrease of P L with strain is stronger. This may be due to the fact that initially the grain boundaries were not sufficiently relaxed such that the initial straining leads to more complex grain-boundary changes. The value of α p for straining a poly-crystalline sample is almost double as large as the value for single-crystal straining. This demonstrates that grain boundaries and their changes during strain have a major effect on the evolution of the pDOS.

Magnetic moments
Bcc-iron is ferromagnetic; in the ideal defect-free lattice, at zero temperature, each atom carries a magnetic moment of µ equ = 2.154 µ B [43], where µ B denotes the Bohr magneton. In the vicinity of defects, the magnetic moment of an atom i, µ i , will change. The physical origin of the change is complex. A major influence is exerted by the so-called magneto-volume effect which predicts that the magnetic moment of an Fe atom depends on its local atomic volume. For a free atom (infinite atomic volume), it is µ = 4.0 µ B ; with shrinking volume the moment monotonically decreases until it vanishes at a critical atomic volume [14]. Since the atomic volumes are changed in the vicinity of defects, the magneto-volume effect influences the magnetic moments locally. A second influence is provided by the local coordination of iron atoms. The coordination effect has been investigated recently by Zhang et al. [44] in Fe. They found that this effect becomes dominant in particular for compressed iron with atomic volumes smaller than the equilibrium volume in bcc-Fe, Ω 0 = 11.78Å 3 .
Our atomistic results on the structure of defective Fe may be used to assess the changes in magnetism via the magneto-volume effect; we ignore the coordination effect, since the local volumes in defect structures are as a rule enhanced with respect to the equilibrium volume, such that the coordination effect will lose significance. In the embedded-atom-model type potential that we use here, the local volume of atom i may be most easily assessed by the value of the electron density n e i at this site. The electron density n e i is obtained as the sum of the contributions of all surrounding atoms, where j denotes the neighbor atoms and r ij their distance to atom i. The function f (r) denotes the electron density provided by a neighbor atom at distance r and is part of the description of the interatomic potential function. Figure 7a displays the dependence of the electron density on atomic volume for a homogeneously expanded/contracted lattice. Derlet and Dudarev [26,43,45,46] used this idea to correlate the magnetic moment of atom i, µ i , directly with its electron density, n e i . They proposed a function It has been argued that the calculation of the local electron density to some extent also incorporates the coordination effect, since it is based on the individual bond lengths r ij [25]. Figure 7b shows ab-initio data of the dependence of µ on the atomic volume and demonstrates that equation (11) is able to fit this dependence well via the local electron density. In the present study, the range of atomic volumes tested is 11-16Å 3 , such that our fit is well justified. Our fit parameters are C = (3.423 ± 0.234) µ B and γ = 0.278 ± 0.046. Here, the error denotes the accuracy of the fit procedure. n e c denotes the electron density Fig. 7. Dependence of (a) the electron density and (b) the local magnetic moment on the atomic volume. Data were calculated for a homogeneously expanded / contracted bcc structure. Ab-initio data obtained by Moruzzi et al. [61] and Herper et al. [62] are taken from reference [43]. The cross marks the equilibrium magnetic moment, µequ = 2.154 µB, at the equilibrium atomic volume of bcc Fe, 11.78Å. at the critical volume, Ω c = 7.750Å 3 , and is only used for normalization purposes. equation (11) results in an equilibrium magnetic moment of µ equ = 2.149 (2.166) µ B for the single-crystalline (poly-crystalline) structure. We note that the dependence of the magnetic moment µ on the atomic volume Ω can be well approximated by a linear dependence in the vicinity of the equilibrium volume, with k v = 0.770. We show in Figure 8a that the magnetic moment depends on the vacancy concentration in a defective crystal in a linear way, with k v = 0.477. This value can be compared to the value of k given above, equation (12), if we assume Ω − Ω 0 ∼ = n v Ω 0 . k v is of similar size, but somewhat smaller than k. This can be justified by the fact that the volume change caused by insertion of a vacancy is somewhat smaller than an atomic volume, due to relaxation effects [47]. For dislocations, the linear relationship is more noisy, see Figure 8b: with k d = 0.388 × 10 −3 (10 16 m −2 ) −1 . Rewriting this expression in terms of the concentration of atoms in the core of a dislocation line, n c , equation (1), we obtain with k d = 0.75. This value is somewhat larger than k v , equation (13), implying -analogous to the effect seen for the pDOS -that the number of core atoms in a dislocation underestimates the effects on the changes in magnetic moments. The influence of grain boundaries on the magnetism can be assessed by comparing the magnetic moments in our single-crystalline (2.149 µ B ) and in a poly-crystalline (2.166 µ B ) specimen. The increase is caused by the open volume generated in the vicinity of the grain boundaries.
The magnetic moments change in a strained specimen. We first discuss the changes under small strains < 1%, where the specimen responds in a linear elastic way. The relative volume change then is given by where ν = C 12 /(C 12 + C 12 ) = 0.37 is the Poisson number when the uniaxial strain is exerted in [001] direction, and the C ij are the elastic moduli, which for our potential are tabulated in [33]. Inserting this dependence in equation (12) leads to with k = 0.16. This dependence is well fulfilled in our simulation (not shown). The data for larger strains, however, deviate from this behavior, see Figure 9. Here, a value of k = 0.045 (0.039) is obtained for the single-crystalline (poly-crystalline) case, considerably smaller than in the linear-elastic case. This behavior is mainly caused by a considerably slower increase of the atomic volume than that predicted by equation (16), since at large strains, the strain-induced stresses relax to decrease the atomic volumes.
In addition, the effect of strain is weaker in the poly-crystalline specimen, presumably since already the grain-boundary-filled polycrystalline material has a higher initial magnetic moment than the single-crystalline specimen. Note that the dislocations created in the singlecrystalline material hardly allow to reach the unstrained value of the polycrystalline material, 2.166 µ B . This comparison thus demonstrates the strong influence that grain boundaries have on the magnetic moment.
For a uniaxial strain, atomic volumes increase as = (Ω − Ω 0 )/Ω 0 and we can therefore immediately compare the values of k under strain with the value of k, equation (12), by more than an order of magnitude.

Conclusions
Using atomistic simulation, we studied the influence of defects on the pDOS and the magnetic moments in bcc-Fe. We found the following.
1. Defects lead to a decrease of the amplitude of the longitudinal phonon peak, P L . In compensation, the pDOS increases in the minima between the peaks, and also in the low-frequency region. 2. The changes induced by different types of defects -vacancies, dislocations, grain boundaries -to the longitudinal peak height are similar. In other words, it appears difficult to use the changes in the longitudinal peak height as a tool to identify the defect types. 3. In contrast to vacancies, both dislocations and grain boundaries lead to an increased population of high-frequency modes, i.e., a broadening of the longitudinal peak. 4. At low defect densities, the peak of the longitudinal phonons decreases in proportion to the defect density. 5. Changes in the magnetic moment of defective crystals are also linear in the defect density. These changes can be explained since defects lead to an increase in the atomic volume of neighboring atoms, and hence in their magnetic moment.  Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix A: Interatomic potentials
A considerable number of potentials have been derived for bcc Fe. We used here the potential by Ackland et al. [33]. However, we also tested the potentials of Chamati et al. [48], Etesami [17,54] for the quality of their pDOS. The data are shown in Figure A.1. Interestingly, none of the potentials reproduces the experimental data [55] accurately.
Since in this study, we focus on the changes in the longitudinal phonon peak, are interested in a fair representation of the position and relative height of the longitudinal and the two transverse phonon peaks. These appear to be best represented by the Ackland et al. potential [33] (and also in the less well kown Olsson potential), and we selected this potential for our study.
The pDOS of the Mendelev et al. [50] potential has also been provided by Talati et al. [56], in fair agreement with our results.
Dragoni et al. [57] provide a comprehensive review of the thermoelastic properties of bcc-Fe in four potentials: The Mendelev potential [50], the Meyer-Entel potential [51], the Ouyang et al. potential [58] and the Marchese et al. potential [59]. Among others they also calculate the phonon dispersions, but not the density of states. They conclude that all considered potentials have their weaknesses; for example, the Mendelev potential [50] shows anomalies in the softening of the C 44 elastic constant.
Previously, Marinica and Willaime [10] calculated the bulk phonon dispersions and phonon densities of state for the potentials by Ackland et al. [33], Mendelev et al. [50] and its improved version [60], and Dudarev and Derlet [45] and found good performance for all potentials at zero pressure. At higher pressure (10 GPa) -which may be relevant in locally strained defective regions of a metalthe performance is not so good.
It should be noted that potentials are usually developed by fitting to basic metallurgical quantities, such as the cohesive energy, lattice constant, elastic moduli, and point defect energies. Depending on the motivation of the developers, other quantities may enter, such as surface energies and relaxations, melting or phase transformation temperatures, etc. The pDOS usually does not belong to the quantities used for potential development. This may explain why the potentials considered in this Appendix show such a variation in the description of the pDOS.
Malerba et al. [17] report on the ability of potentials to describe dislocation cores in Fe correctly. They note that among the potentials studied here Ackland et al. [33], Mendelev et al. [50], and Marinica et al. [17,54] perform well in this regard. This justifies our use of the Ackland potential for modeling plasticity in Fe.