Lattice dynamics in the double-helix antiferromagnet FeP

We present a comprehensive investigation of lattice dynamics in the double-helix antiferromagnet FeP by means of high-resolution time-of-flight neutron spectroscopy and ab-initio calculations. Phonons can hybridize with the magnetic excitations in noncollinear magnets to significantly influence their properties. We observed a rich spectrum of phonon excitations, which extends up to $\sim$50 meV. We performed detailed analysis of the observed and calculated spectra for all high-symmetry points and high-symmetry directions of the Brillouin zone. We show that the DFT calculations quantitatively capture the essential features of the observed phonons, including both dispersions and scattering intensities. By making use of the detailed intensity comparison between the theory and the data, we were able to identify displacement vectors for the majority of the observed modes. The overall excellent agreement between the DFT predictions and the experimental results breaks down for the lowest mode at the $Y$-point, whose energy is lower than calculated by $\sim$13%. The present study provides vital information on the lattice dynamics in FeP and demonstrates applicability of the DFT to novel pressure-induced phenomena in related materials, such as MnP and CrAs.

We present a comprehensive investigation of lattice dynamics in the double-helix antiferromagnet FeP by means of high-resolution time-of-flight neutron spectroscopy and ab-initio calculations. Phonons can hybridize with the magnetic excitations in noncollinear magnets to significantly influence their properties. We observed a rich spectrum of phonon excitations, which extends up to ∼50 meV. We performed detailed analysis of the observed and calculated spectra for all high-symmetry points and high-symmetry directions of the Brillouin zone. We show that the DFT calculations quantitatively capture the essential features of the observed phonons, including both dispersions and scattering intensities. By making use of the detailed intensity comparison between the theory and the data, we were able to identify displacement vectors for the majority of the observed modes. The overall excellent agreement between the DFT predictions and the experimental results breaks down for the lowest mode at the Y -point, whose energy is lower than calculated by ∼13%. The present study provides vital information on the lattice dynamics in FeP and demonstrates applicability of the DFT to novel pressure-induced phenomena in related materials, such as MnP and CrAs.

I. INTRODUCTION
Transition-metal monophosphates AP (A = Fe, Mn) are metallic binary compounds that attracted significant attention due to their magnetic and electronic properties. Both materials are itinerant helimagnets featuring an unusual "double helix" type of magnetic ordering with the propagation vector aligned along the c-axis of the orthorhombic crystal structure [1][2][3]. In addition, an unusual orbital density-wave ordering was recently observed in MnP by means of resonant inelastic x-ray scattering. The orbital ordering has half the period of the magnetic helix and shows a similar temperature dependence which implies a coupling between the two orders [4]. This feature is reminiscent of the nematic orbital order in the iron-based superconductors LaFeAsO and BaFe 2−x Co x As 2 [5][6][7].
Application of hydrostatic pressure in MnP suppresses the magnetic transition temperature down to 0 K and induces a quantum phase transition toward a nonmagnetic state. Pressure-induced variation of the magnetic spiral structure of MnP has been studied by nonresonant magnetic x-ray diffraction [8], muon-spin rotation [9], and neutron diffraction [10]. It was shown that the ambient-pressure incommensurate helical structure shortens, becomes commensurate, and rotates its propagation vector to the b-axis at a pressure of ∼2 GPa. This suggests a delicate balance of competing magnetic interactions [8][9][10]. Surprisingly, an unconventional superconducting phase was observed below T c ≈ 1 K at P ≈ 8 GPa, very close to the quantum critical point, making MnP the first ever observed Mn-based superconductor [11].
Even though the electronic [12][13][14][15][16], transport [17][18][19][20], and magnetic [1-3, 21, 22] properties of MnP and FeP have been intensively studied by different experimental techniques, there are no reports on the structural dynamics in these materials. Structural vibrations-phonons-are known to serve as a gluing mechanism for Cooper pairing in BCS theory of superconductivity. Moreover, phonons can hybridize with magnetic excitations in noncollinear magnets [23,24], which may significantly influence their properties. In the isostructural material MnAs (in its orthorhombic phase, which is stable between 315 and 393 K), giant coupling between a phonon soft mode and magnetic moments plays a crucial role in the magnetostructural phase transition, as evidenced by first-principles calculations [25][26][27].
Thus, knowledge of phonon dynamics can provide essential information on unusual physics in these materials.
Unlike MnAs, which is paramagnetic in its orthorhombic phase, another 3d-metal pnictide isostructural to FeP and MnP-CrAs-orders in a double-helix magnetic structure. [28]. The magnetic and lattice degrees of freedom were also found to be coupled in this material. The magnetic transition in CrAs is accompanied by a first-order isostructural transition which manifests itself by a large abrupt expansion of the b-axis and a slight reduction of a and c-axes [29]. Similarly to MnP, pressure-induced superconductivity was discovered in CrAs in the vicinity of the helical phase [30][31][32][33][34], and the magnetic propagation vector in CrAs undergoes the reorientation from c to b upon application of moderate pressure [35].
In FeP, the magnetic, lattice, and electronic degrees of freedom seem to be coupled to a lesser extent. Thus FeP can be considered as a model system for the phonon dynamics in the 3d-metal monopnictides family. In this work, we present arXiv:2009.07267v2 [cond-mat.str-el] 18 Sep 2020 the results of a comprehensive investigation of structural dynamics in FeP. We combine thermal-neutron time-of-flight spectroscopy and ab-initio phonon calculations in order to fully characterize the phonons over the entire Brillouin zone. The remainder of the manuscript is organized as follows: Section II describes the experimental details and the firstprinciples calculations. In Section III we discuss results of the measurements and compare the obtained neutron spectra with the calculations. The analysis is first presented for all high-symmetry points in the Brillouin zone, where the majority of lattice vibrations are resolved and identified. Then, the data on the phonon dispersions along high-symmetry directions are shown. In section IV, we summarize the results.

II. METHODS
The crystal structure of FeP is depicted in Figs. 1(a)-1(c). It is described by the space group P nma (no. 62) with the lattice parameters a = 5.197 Å, b = 3.099 Å, and c = 5.794 Å at room temperature [1,36]. Both Fe and P occupy the 4c Wyckoff position with the parameters x = 0.002 and z = 0.200 for Fe and x = 0.191 and z = 0.569 for P. It is worth noting that the P nma crystal structure of FeP can be considered as a distorted hexagonal NiAs-type structure (P6 3 /mmc space group), where the distorted hexagonal layers of Fe ions are stacked along the orthorhombic a-axis. Figure 1(c) depicts the crystal structure of FeP as viewed from [100]. Red circles highlight the Fe ions that have approximately the same coordinate along the a-axis, namely, x 1 = 0.048 and x 2 = 0.502. The highlighted sites form the buckled triangular layers that are contracted along the c-axis such that the isosceles triangle of neighbouring ions has the base angle of ∼56 • . The next Fe layer is then elongated along c with a ∼66 • angle. Green circles in Fig. 1(c) highlight the P ions that form a triangular layer with a larger buckling along a (the relative difference of the x coordinates for P layers x 1 − x 2 = 0.118), but with a smaller distortion along c.

A. Experimental details
High quality single crystals of FeP up to 500 mg in mass and 80 mm 3 in volume were grown by chemical vapor transport with iodine as a transport agent. The quality of the crystals were confirmed by means of EDX, XRD, magnetic, transport, high resolution TEM and single-crystal neutron diffraction experiments [36]. The temperature dependence of the resistivity shows the typical metallic behavior with a residual resistivity ratio RRR = 566. Below the the Néel temperature T N = 119.3 K, analysis of the neutron diffraction data yields an incommensurate magnetic structure with the propagation vector Q = (0, 0, ±δ), where δ ≈ 0.2 [1,36]. For neutron spectroscopy measurements, we used a single crystal with a mass of ∼0.5 g.
The experiment was conducted at the thermal-neutron direct-geometry time-of-flight spectrometer Merlin [37] lo- The sample was oriented with its [100] axis vertical, thus the (0K L) reciprocal plane was horizontal throughout the measurements. In addition to a large 2θ angular range of 135 • in the horizontal plane, the detector banks of Merlin also allow for detection of scattered neutrons up to 30 • in 2θ in the vertical plane. This enabled us to collect the data on the elemental excitations for all principal directions in reciprocal space within a single experimental setup. The incident neutron energy E i was set to 60 meV to achieve a reasonable experimental resolution while sampling a sufficiently large volume the momentum-energy space. The trade-off between the resolution and intensity was optimized by setting the chopper frequency to 450 Hz. In this setup, the resulting energy resolution at the elastic line is estimated as ∼3 meV. An additional dataset with a reduced counting time was collected at E i = 40 meV. During the measurements, the crystal was gradually rotated over 70 • around the [100] axis with a step of 0.5 • . To minimize anharmonic effects on the phonon spectra, all the measurements were performed at T = 6 K. The collected data were reduced and analysed using the HORACE software [38]. A symmetrization procedure was applied during data reduction, which means that the data from equivalent Q-directions in momentum space [for example, (H00) and (−H00)] were averaged to improve statistics. Thus, the covered Q-space was folded down to a 90 • sector in the (0K L) plane, which is irreducible for the orthorombic system. To plot the I(Q, E) intensity colormaps for different momentum directions, the data were integrated over ±0.2 Å −1 along the perpendicular momenta. The scattering intensity profiles in the high-symmetry points of the BZ (I(E)) were integrated over ±0.1 r.l.u. along each principal momentum direction [(H00), (0K0), and (00L)].
Phonon energies were obtained with the Phonon Explorer software, which works with unsymmetrized raw data [39]. It effectively combines statistics from all Brillouin zones where the phonon of interest has appreciable intensity as described in Sec. IIID.

B. First-principles calculations
Lattice dynamics calculations were carried out by means of the projector-augmented wave (PAW) method [40] and density functional theory (DFT) as implemented in the VASP software [41,42]. The generalized gradient approximation (GGA) functional with Perdew-Burke-Ernzerhof (PBE) parametrization [43] was used. The plane-wave cutoff was set to 600 eV. An 8 × 13 × 3 (Monkhorst-Pack scheme) [44] k-point mesh was used for Brillouin zone integration. The following electronic configurations was chosen: 3d 7 4s 1 for Fe, and 3p 3 3s 2 for P. Magnetic moments of Fe ions were taken into account by performing spin polarized calculations. The phonons were calculated by constructing a supercell (2×2×2) and calculating the force constants as implemented in PHONOPY [45]. Classification of vibrational modes was performed by means of PHONOPY. The notation used for the irreducible representations corresponds to ISOTROPY [46].
In an INS experiment, the dynamic structure factor S(Q, E) is measured, which takes the following form for the lattice vibrations: where E and ε d,s are the energy and the polarization of the phonon with the mode index s and the wave-vector q, τ is a reciprocal lattice vector, r d is the position vector for the atom d in the unit cell, M d its mass, and b d its coherent cross-section; W d is the Debye-Waller factor, and n s is the Bose-Einstein distribution. The INS intensity was simulated using the OCLIMAX software [47] which uses vibrational frequencies and polarization vectors from first-principles calculations (PHONOPY in this work) as the input. The simulated spectra were convolved with a Gaussian function with the standard deviation σ = 2 meV to model the experimental resolution.

A. Dynamic structure factor
We first discuss the INS spectra of FeP for momentum transfer along the [00L] (c*) direction. Figure 2(a) shows such a spectrum for energies E up to 35 meV and momenta |Q| up to ∼10 Å −1 . The spectrum reveals a clear sine-shaped mode that disperses within the ∼(25-35) meV range and acquires a noticeable intensity at |Q| greater then ∼4 Å −1 , such that two periods of the dispersion are included within the covered momentum transfer. The maxima of the dispersion coincide with |Q| ∼4.32 and ∼6.48 Å −1 , which corresponds to (004) and (006) reciprocal-lattice points, respectively, whereas the minima are found at the (00L) points with L = 5 and 7 r.l.u. A more detailed shape of the same mode is demonstrated in Fig. 2(b), where the data collected with a higher resolution are presented. The spectrum in Fig. 2(a) also features another mode with a much steeper dispersion emanating from the elastic (E = 0) line and reaching an energy of ∼27 meV at (00 15 2 ). The first mode is identified as a longitudinal optical phonon branch. The second mode represents the longitudinal acoustic vibrations.
The elementary excitations of the helimagnet FeP at low temperatures consist of phonons and spin excitations, which can be separated by analysis of the Q-dependence of the intensity of a particular excitation. The intensity of the magnetic excitations is proportional to the magnetic form-factor, which drops rapidly with increasing Q, whereas the phonon scattering is proportional to Q 2 and leads to a larger crosssection for the increasing momentum transfer. the I ∝ Q 2 trend drawn by the dashed line. Small deviations from the quadratic function are due to additional modulation by the structure factor, which alters the spectral weight of the optical phonon branch according to L = 2n + 1 for even/odd values of n at the (00L) point. The observed spectral weight of each phonon branch at a specific reduced wave-vector q can significantly vary at different absolute momenta Q. This enables optimizating the scattering geometry for the excitations of interest and performing mode assignments based on comparison of the measured and calculated intensities. Figures 3(a)-3(d) demonstrate the change of the INS intensity in the E-Q maps due to the phonon polarization factor, I ∝ (Q · ε) 2 (ε -the displacement vector), and the structure factor, I ∝ d e iQ·r d 2 , of the multiatomic unit cell. In Fig. 3(a), the data were plotted along the closed path between the Γ -points of the BZs adjacent to (H K L) with H = 1 and K and L changing between 1-3, and 5-7 respectively, such that the momentum path reads: (125)-(136)-(127)-(116)-(125). The interval (3) is a replica of the interval (1) and features a steep acoustic mode dispersing up to ∼33 meV, whereas a mode with a bandwidth of ∼16 meV is observed at the intervals (2) and (4). Because the reduced wave-vector is almost parallel to the absolute momentum (q Q) along the (125)-(136) and (127)-(116) paths, the intensity of the longitudinal acoustic mode is maximized. The opposite holds for the momenta connecting (136) with (127) and (116)

B. Spectra at the high-symmetry points
It is essential to take into account the variation of the dynamic structure factor when the experimental data are compared with the results of the lattice-dynamics calculations. In order to resolve different modes in the experiment, it is also required to analyze the INS intensity in many different BZs. In Figs. 4(a)-4(d) and 5(a)-5(a)-5(d), and in Table I, we summarize the comparison between the experimental and calculated energies of the phonons at all eight highsymmetry points of the BZ. The constant-Q cuts through the data (the intensity as a function of energy) at the Γ −point are shown in Fig. 4(a). To extract the exact peak positions, the experimental data (symbols) were fitted by a sum of Gaussian functions (solid curve) with three free parameters per peak: the peak center, the peak width (FWHM), and the peak amplitude. The results of the simulations are plotted along with the data for a direct comparison (shaded area).
As can be seen, the theory predicts a strong peak at ∼35 meV and a weaker one at ∼40 meV at Γ (017). Accord- ing to their displacement vectors, these phonons correspond to the vibrations with the A g symmetry. Consistent with the predictions, the experimental profile shows two peaks with correct intensity ratio and positions. The first peak is at E exp = 33.6 meV (E calc = 35.3 meV) and the second is at E exp = 39.6 meV (E calc = 41.5 meV). The energies of both peaks are slightly lower than the calculated values, by 1.7 and 1.9 meV, respectively, which translates to the relative deviations, ∆E/E = E exp − E calc /E exp , of −5.1% and −4.8%. The spectrum at Q = Γ (016) exhibits a strong peak at E ∼25 meV that is in a good agreement with the simulated spectrum, which predicts only one peak at the zone center in the (016) BZ. According to the calculated eigenvectors, the observed excitation can be identified as a B 3u phonon. Overall, the profiles shown in Fig. 4(a) for seven different BZs demonstrate good agreement with the DFT results. Every INS spectrum at Γ is characterized by dominant intensity of only a few modes, which can be well resolved in energy in different zones, except for the Γ (134) spectrum where two modes have close energies and similar intensities merging into a single peak (E calc = 41.6 and 41.9 meV) at E exp = 42.5 meV. The complete list of the observed and calculated mode energies is given in Table I.  These peaks were identified as the vibrations with S 1 + S 2 symmetry observed at ∼17 and ∼30 meV. Three modes are observed at T (0 3   It is worth noting that the lowest-energy phonon at S and T is softer than the calculated value by (8-12)%. On the contrary, the relative disagreement between E calc and E exp for the vibrations that have higher energies at S and T amounts to only (1-7)%. At the U-point, the low-lying mode is observed at the energy that is by 6% higher than the predicted value. The softening (hardening) of the above mentioned modes at the corners of the BZ corresponds to the softening (hardening) of a transverse acoustic mode with momentum along Γ -Y (Γ -Z).
The INS spectrum at the diagonal R-point is presented in Fig. 5(a). The R and S-points are the only wavevectors at which the modes have four-fold degeneracy (there are six four-fold-degenerate modes at R and S). All modes have two-fold degeneracy at all the other zone-edge highsymmetry points. The spectrum at R ( 1   2   3  2   11 2 ) reveals two excitations with E exp = 16 and 29 meV. These are the lowest two four-fold-degenerate modes with R 1 + R 2 symmetry. The second mode closely matches the calculated energy (E calc = 30.7 meV) but the first mode demonstrates a slightly larger discrepancy (E calc = 17.2 meV). The mode softening is less pronounced than the softening of the lowest-energy mode at S and T . This is expected as the R-points is further away from Y , where the largest mode softening (with respect to our DFT results) is observed.
The low-energy lattice vibrations at X can be observed at a number of BZs (e.g. X ( 1 2 24), X ( 1 2 34), or X ( 1 2 25)) as drawn in Fig. 5(b) where the highest INS intensity of these excitations (with X 1 and X 2 symmetry) is found at the (025) BZ center. The 17.8 meV peak at X ( 1 2 25) represents a mixture of two modes with close energies. Similarly, two phonons with energies of 37.2 and 39.8 meV are predicted to have spectral weight at X ( 1 2 24), where they are not resolved and seen as a single 37.1-meV peak.
The intensity profile of Y (0 5 2 5) [ Fig. 5(c)] reveals a transverse acoustic mode along the Γ -Y path (the (0K0) direction) that exhibits large softening compared with the DFT results. The experimental energy E exp = 16.6 meV deviates by 2.6 meV from the simulated value (E calc = 19.2 meV). The other modes at Y are in a better agreement with the theory. The transverse acoustic mode at Z reaches 18.4 meV, which is 1.9 meV larger than E calc . (See the spectrum at Z(13 9 2 ) in Fig. 5(d)). A number of the other modes were observed at the Z-point gthat are in a good agreement with the E calc values.
The values of E exp and E calc for all the experimentally determined vibrations at the high-symmetry points of the BZ are listed in Table I. If the relative discrepancies ∆E/E of all the observed modes are considered together, it can be found that the average ∆E/E amounts to ∼ −3% and its standard deviation is ∼4.5%. This indicates that our DFT calculations tend to slightly overestimate the phonon energies. If the calculated dispersions are renormalized by a factor of 0.97, the agreement between E exp and E calc are found within ±5% for all the modes, except for the transverse mode at T , S, and Y (softening, ∆E/E = 5%-13%, from T to Y ), and at X , U, and Z (hardening, ∆E/E = 7%-13%, from X to Z).

C. Displacement patterns
In the previous section we have shown that our calculations quantitatively reproduce the dispersion and intensity distribution of the observed phonon spectrum at all high symmetry points of the BZ. This allowed us to identify the   character of the lattice vibrations by analyzing the calculated eigenvectors of the observed phonons. Figures 6 and 7 show selected sketches of the lattice vibrations resulting from the calculated eigenvectors at the Γ -and other highsymmetry points of the BZ, respectively. Note that at the Γ -point we present only those phonons which were identified from the experiment [see Tab. I]. Fig. 7 shows the lowest-energy phonons for each high-symmetry point at the BZ boundary. The vibrations that have energies E < 40 meV mainly involve Fe ions, whereas P ions contribute mostly to the higher-energy phonons. These results agree well with naive expectations, because Fe atoms are approximately twice as heavy as the P atoms. It is worth noting that the phonons depicted in Figs. 6 and 7 mainly involve atomic displacement orthogonal to the [100] crystal direction. That is because our experimental geometry we had limited access to the (H00) direction of the reciprocal space. Therefore, we could hardly detect the phonons with polarization (displacement vectors) parallel to the a axis due to the polarization factor (Q · ε) 2 in the INS cross-section Eq. (1).
It is convenient to consider the displacements of Fe ions of the first five vibrations of Fig. 6 with respect to the quasihexagonal crystal structure (see Fig. 1(c)). The 23.4 meV excitations (B 1g ) involve the displacements of the Fe ions along the b-axis, such that the three neighbouring sites within the triangular layers cause the triangles to rotate around the axis perpendicular to the layers (the a-axis). The optical phonons of 26.5 and 33.5 meV represent the parallel translation of the adjacent triangular layers (in the opposite direction) along b and c, respectively. The B 3u and B 2g vibrations (24.9 meV and 37.8 meV) correspond to a contraction/elongation of the triangular layers along b. The 24.9 meV phonon recovers the hexagonal symmetry for the one Fe layer, but leads to a larger distortion of the other layer, whereas the 37.8-meV excitation leads to the symmetry recovery for both layers T (1 0.5 6.5) (the elongated layer is contracted and the contracted layer is elongated along b). The A g vibrations (39.6 and 47.1 meV) involve a complex displacement of P ions in the ac plane, and the phonons of 42.5 meV are characterized by a shift of P ions along b, which corresponds to a parallel translations of the P triangular layers along the b-axis (either opposite to the adjacent layers or in the same direction).
The lowest-energy phonons at the BZ boundary can be described as follows. The excitations at the symmetry points X , Z, and U correspond to the displacements only along b (of both Fe and P ions), whereas the vibrations at Y , T , and S have complex displacement patterns , where different sites move along different crystallographic directions. The 16.6 meV phonon of the Y -point (Y 2 symmetry) is of particular interest as this is the excitations that shows softening as compared to the results of our DFT calculations. It involves a larger displacement of two Fe ions along b, and a smaller displacement of two P ions in the ac plane. The movement of the Fe sublattice corresponds to the rotation of triangles within the layers where the adjacent layers rotate in the same direction.

D. Phonon dispersions
Having discussed the phonon modes at the high symmetry points of the BZ, we turn to following the phonon dispersion along the high symmetry directions. Figure 8 summarizes the experimental data for the Γ -X , Γ -Y , and Γ -Z paths in the BZ, which coincide with the (H00), (0K0), and (00L) directions, respectively. Because the INS spectral weight for different modes varies drastically as reciprocal space is traversed, one has to consider many different BZs, as shown in Fig. 8.
Only the 24.9 meV mode bears a significant spectral weight at Γ (014). At small reduced momenta along (H14), the spectral weight is transferred to the two lowest energy optical modes, which disperse downwards and reach the energies of ∼17.8 meV at the X -point at the BZ boundary. At the momenta between ( 1  2 14) and (114), the entire observed spectral weight transfers to the TA modes, which disperse in a good agreement with our ab initio calculations. As can be seen in Fig. 8, the spectral weight is altered when the momentum is changed from L = 4 to L = 5 for the same H and K. In the (015) BZ, the A 2 g mode has the highest intensity. It weakly disperses downwards and fades out on approach to the zone boundary. The spectral weight then transfers to the weakly dispersing B 1 2g mode in the (115) zone.
The INS intensity follows the structure factor. Thus, the same modes that are seen along the (H14) for 0 < H < 1 can be observed in (H25) and (H34), which yields the rule K, L± 1 for a given H [ Fig. 8]. In addition, the polarization factor makes the 39.6 meV mode visible, as it has a longitudinal polarization character. Further, the dispersion of the B 2 1g and B 1 2u modes can be followed in the (134) BZ. The dispersion relation of the B 1 3g and B 1 1g modes can be followed in the (H24) and (H35) BZs for 0 < H < 1, where they carry the maximal spectral weight. No spectral weight is present for the acoustic mode at these zones due to the structure factor.
The acoustic phonon intensity obeys the following systematic absence rules for the dispersion along the Γ -Y path, as we show in Fig. 8. The acoustic modes have nonzero intensity at the (125), (136), and (107), which is K + L odd for H = 1, whereas only the optical modes are observed in the (0K L) zones. The same applies for the acoustic dispersions along Γ -Z.
As mentioned in Sec. IIIB, the Y 2 mode shows a notable softening compared to the prediction of our ab initio calculations. The dispersion of this mode in the vicinity of the Y -point is clearly observed in the (0 2 + K 5) and (0 2 + K 6) zones for 0 < K < 1, where it shows a much steeper upward dispersion than predicted by the calculations. On the contrary, the Y 1 mode has a shallower dispersion in perfect agreement with the theory, as can be seen in the vicinity of the points (1 5 2 5), (1 5 2 6), and (1 1 2 7) in Fig. 8. The Γ -Z dispersion of the LA mode is characterized by large intensity in the (116) zone due to the polarization factor. Weaker intensities of both the LA and TA modes are observed at (134). The calculated energies and dispersions of the optical phonons are in good agreement with the INS data along the Γ -Z path as evidenced in Fig. 8.
We also resolved the phonon dispersions at momenta connecting all the high-symmetry points of the BZ boundary. Figure 9(a) shows a cut through the INS data along the momentum path Z-T -R-U-Z. To account for the spectral weight alteration in different BZs, we combined the spectra at the (015) and (125) BZs for Z-T and the (107) and (007) BZs for U-Z. The INS spectra in Fig. 9(a) cover the first five modes dispersing between ∼15 and ∼40 meV. All observed modes exhibit strong dispersion except along the T -R path, where a flat 16 meV mode is observed. This is in agreement with the lowest energy mode predicted by the DFT calculations. The data for the X -U-R-S-X and Y -S-R-T -Y paths are compared with the calculated dispersions in Figs. 9(b) and 9(c), respectively, which also demonstrate a relatively good agreement for all momenta except for the lowest-lying mode along T -Y . As can be seen from Fig. 9(c), the Y 2 mode softens (as compared to the DFT predictions) not only along Γ -Y (Fig. 8), but also along T -Y , where it is flat. Meanwhile, the closely-spaced Y 1 modes are in good agreement with the theoretically-predicted dispersion along T -Y . All the modes retain the four-fold degeneracy along R-S as a result of a pair of the nonsymmorphic symmetries of the P nma space group [14].
To compare the observed optical phonon energies with the results of our first-principles calculations in detail, we performed multi-zone fitting uning the entire unsymmetrized INS dataset [48]. In multi-zone fitting, the energy and linewidth of each mode at a given reduced momentum is extracted from the experimental data in many BZs simultaneously (for the same reduced momentum). This allows one to significantly improve statistics and to resolve modes that are close in energy but whose spectral weights differ in many BZs. During the multi-zone fitting, the intensity of each peak is allowed to vary in every BZ as each mode has different spectral weight in different BZs. However, the positions and widths of the peaks are fixed to be the same in every zone and refined globally. As a result, the extracted peak positions from all the BZs for the same crystal momentum direction can be plotted together, as shown in Figs. 10(a)-10(c) for the Γ -X , Γ -Y , and Γ -Z paths, respectively. To perform the multi-zone fitting, we used the PHONON EX-PLORER software [39,49] to fit constant-momentum cuts assuming Gaussian peak shapes for every phonon. Linear background was subtracted from every constant-momentum cut before the multi-zone fit was performed.
The positions of the acoustic modes at low energies below 10 meV were extracted from constant-energy cuts. For this, the I(Q) profiles were fitted with a Gaussian function for every fixed energy.
As can be seen in Figs. 10(a), a good agreement between the DFT results and experimental data is found for all re-solved modes for Γ -X . The slope of the TA modes matches the predicted energies for all momenta from the center of the zone to the zone boundary. However, due to a finite coverage of the (H00) direction in our experimental setup, it is not possible to verify the DFT predictions of the dispersion of the LA mode for Γ -X . A good agreement is also seen for many weakly dispersing optical phonon branches in the 20-40 meV range.
Because the LA and TA modes along the Γ -Y path are close in energy, we fit the constant-momentum INS data for E < 10 meV with only one peak function. Figure 10(b) shows that the experimental energies of the acoustic modes for the reduced momentum below 0.1 r.l.u. agree well with the DFT predictions. As was discussed in relation to Fig. 8, the Y 2 mode exhibits a significant softening of its dispersion in the vicinity of the Y -point. The softening of this lowenergy mode as compared to the DFT results is clear in Fig. 10(b). The discrepancy starts at ∼0.3 r.l.u. and increases as the momentum approaches the Y -point at 0.5 r.l.u. The other optical modes demonstrate a good agreement at all momenta. Figure 10(c) summarizes the extracted peak positions for Γ -Z. In this reciprocal-space direction, the TA and LA modes have significantly different slopes and could be resolved seperately in our INS data. The experimental and theoretical dispersions of the LA mode follow the same trend, but the experimentally obtained mode has slightly higher energies for all the momenta. The dispersions of the TA branches are found to be in good agreement. Similarly to Γ -X and Γ -Y , the optical dispersions for Γ -Z show good agreement between the DFT results and the INS data.

IV. DISCUSSION AND CONCLUSIONS
Since the magnetic structure of FeP is noncollinear, a significant renormalization of its magnon and phonon spectra can be expected in the presence of magnetoelastic coupling. A linear magnon-phonon term, which is absent in the case of collinear spin order, becomes nonzero in the presence of a noncollinear magnetic structure. This is a consequence of the absence of a global spin-quantization axis which allows for a one-magnon term in the spin Hamiltonian [23,24]. This also means magnon-phonon hybridization causes stronger INS spectral renormalization than the magnon-magnon and the phonon-phonon interactions [50][51][52].
The present comprehensive investigation of the lattice dynamics shows that the spin-lattice coupling is seemingly not playing a major role in FeP. This is evidenced by a very good overall agreement between the observed phonon spectra and the DFT calculations that do not take into account the magnetoelastic coupling. The ordering temperature of FeP yields an approximate expected magnon bandwidth of an order of a few dozen meV, or more, if a significant frustration of the exchange interactions is present. An estimate of the magnon bandwidth in FeP can also be made by comparison with CrAs, which is isostructural and orders into the same double-helix magnetic structure. Recent INS measurement of a CrAs powder sample showed that the magnetic excitations reach at least ∼100 meV. Taking into account that T N of CrAs is twice higher than that of FeP, an estimate of a 50 meV magnon bandwidth can be made for the latter. Since the influence of the magnon-phonon coupling is typically stronger in the vicinity of the crossing point of the bare modes [52], it is unlikely that the observed softening of the Y 2 phonon is caused by the overlapping magnon branches as those are expected at a higher energy at the BZ boundary. However, detailed INS measurements of the spin excitations of FeP are required to exclude the magnetoelastic origin of the softening of the Y 2 mode, which will be addressed in future studies.
Measurements of the temperature dependence of the phonon energies in the vicinity of T N might also provide useful information on the magnon-phonon coupling in this compound. For example, the lattice vibrations of the isostructural CrAs at the zone center were recently studied by Raman scattering [53]. Sen et al. [53] resolved all four fullsymmetric A g modes with the energies of ∼14, 23, 28, and 32 meV. The energies of the A g modes of FeP determined in our study (∼22-47 meV) are higher due to a lighter mass of the P ions, in agreement with the theoretical expectations. It was shown that the lowest-energy A g mode in CrAs undergoes a large (∼12%) renormalization at T N . It is therefore very intriguing to address the energy shift and the broadening of the phonon modes in FeP across the magnetic ordering temperature in future studies, for which the full characterization of the phonons presented in our work is indispensable. The detailed phonon analysis of FeP should serve as a foundation for further studies of the lattice dynamics of other representatives of the 3d-metal monopnictide family, where a more complex interplay between the charge, orbital, spin, and lattice degrees of freedom might be present.
The overall success of the DFT calculations to capture the lattice dynamics in FeP, as presented in our study, demonstrates that the DFT method may be applicable for investigations of the pressure-induced phenomena in related materials, such as the strong anisotropic compression and isostructural transitions in MnP [15] and CrAs [54]. This may allow one to address the important problem of the pressureinduced superconductivity in these compounds.
To conclude, we have performed INS measurements of lattice dynamics in FeP. The phonon spectra were collected in a large part of reciprocal space for energies up to 50 meV, which allowed us to resolve the majority of the acoustic and optical modes. We discussed momentum-dependence of the phonon spectral weight across many BZs and demonstrated the role of the nuclear structure and polarization factors in the variation of the observed INS intensities. With the help of accurate ab initio calculations, we were able to assign specific symmetries to the observed excitations and determine the displacement patterns for all modes. A detailed comparison between the theoretical phonon energies and the experimental data was considered at all high-symmetry points in the BZ, as well as all momenta along the highsymmetry directions in the entire BZ. The experimentally resolved phonon dispersions and the observed intensities of the specific modes showed a good agreement with the DFT predictions, except for the lowest-energy mode at the Y -point of the BZ boundary, which exhibit a noticeable softening of ∼13% as compared to the calculations. The origin of this discrepancy will be addressed in future studies.