Strong anisotropy of electron-phonon interaction in NbP probed by magnetoacoustic quantum oscillations

In this study, we report on the observation of de Haas-van Alphen-type quantum oscillations (QO) in the ultrasound velocity of NbP as well as `giant QO' in the ultrasound attenuation in pulsed magnetic fields. The difference of the QO amplitude for different acoustic modes reveals a strong anisotropy of the effective deformation potential, which we estimate to be as high as $9\,\mathrm{eV}$ for certain parts of the Fermi surface. Furthermore, the natural filtering of QO frequencies and the tracing of the individual Landau levels to the quantum limit allows for a more detailed investigation of the Fermi surface of NbP as was previously achieved by means of analyzing QO observed in magnetization or electrical resistivity.


I. INTRODUCTION
Probing the propagation of ultrasound in the quantum regime of electrons yields detailed information on the nature of electron-phonon interactions. The ultrasound velocity in such regime exhibits quantum oscillations (QO), which can be understood both from a thermodynamic argument [1, 2] as well as from a self-consistent treatment of ultrasound propagation as a stream of acoustic phonons interacting with an electron gas that is quantized into Landau levels (LL) [3][4][5][6]. Both approaches yield the same result, namely, the amplitude of the QO being dependent on the (effective) deformation potential Ξ k i = dE k /dε i , which is a measure for the change of energy E k of an electronic band k at a given strain ε i . The connection to the microscopic picture can be understood intuitively by recalling that the probability for an electron in the k-th band of being scattered by a phonon-mode corresponding to ε i is proportional to (Ξ k i ) 2 [3][4][5][6][7][8][9]. Employing measurements of magnetoacoustic QO, the deformation potential and its anisotropy have been experimentally determined for many metals and semimetals (see for example Refs. 3,[9][10][11][12][13][14]. Recently, the semimetallic transition-metal monopnictide NbP is of great interest, mainly due to its symmetryprotected crossings of conduction and valence bands which potentially host Weyl fermions [15][16][17]. It exhibits a very small and highly anisotropic Fermi surface, consisting of intercalated spin-split pairs of electron and hole pockets due to spin-orbit coupling [18]. The small Fermi surface gives rise to pronounced QO of relatively low frequencies, which have so far been observed in magnetization [18][19][20], electrical resistivity [21][22][23][24], Hall resistivity [21,23], thermal conductivity [19], thermopower [19], and * clemens.schindler@cpfs.mpg.de † johannes.gooth@cpfs.mpg.de heat capacity [19]. The superposition of QO originating from different extremal Fermi-surface orbits yield a rich Fourier spectrum, especially when H is aligned along the c axis of the tetragonal lattice and the extremal orbits are the smallest. The peaks in the Fourier spectra could be assigned to orbits via comparison of experimental data to ab initio density functional theory (DFT) calculations [18], however, ambiguities due to the limited resolution and the broadness of the Fourier peaks remained. In a recent study by some of the authors [23], the evolution of the Fermi surface upon direct application of uniaxial stress along the a axis has been probed by means of Shubnikov-de Haas (SdH) oscillations in the electrical resistivity. These experiments revealed a strong strain dependence of the SdH oscillations, which, besides the additional information regarding the orbit assignments, also render NbP a promising platform for studying magnetoacoustic QO. Furthermore, the strong anisotropy of the Fermi surface is suggestive of a highly anisotropic electron-phonon interaction as well, which can be most conveniently investigated via ultrasonic measurements.
In this paper, we report on the measurements of QO in the ultrasound velocity and attenuation in a NbP ) corresponding to the elastic moduli C 11 , C 33 , C 44 , C 66 , and (C 11 − C 12 )/2 (using Voigt notation). Here, u is the displacement vector and q is the direction of propagation of the acoustic wave. Significant differences of the individual QO amplitudes between the modes were revealed. A large signal-tonoise ratio, the usage of pulsed magnetic fields beyond the quantum limit, the high quality of our sample resulting in peak-shaped QO (presence of higher harmonics of the Fourier series), and the natural filtering of certain QO frequencies due to the anisotropic electron-phonon interaction allowed for a detailed analysis of the QO fre-quencies and amplitude ratios. Thereby, the anisotropy of Ξ k i and partially also the cyclotron masses, cyclotron mobilities, and phase factors for several extremal Fermisurface orbits were determined. The QO frequency spectrum could be analyzed via direct assignments of the LL peaks rather than Fourier analysis as in previous studies, which allowed for the assignment of formerly elusive orbits. In addition, the extremal nature (maximum or minimum) of the individual orbits could be deduced from the asymmetric shape of the LL peaks.

II. METHODS
NbP has a tetragonal crystal lattice (space group I4 1 md, no. 109) with the lattice parameters a = b = 3.3324(2) A and c = 11.13705(7) A [25]. A singlecrystalline sample of NbP was grown using chemical vapor transport reactions; the sample has also been used in our previous work [23] for the determination of the elastic moduli. For acoustic modes propagating along one of the main axes, the sample was cut accordingly to a cuboid-shape of dimension 1.92×1.80×0.88 mm 3 . For the (C 11 − C 12 )/2 mode, two cuts parallel to the (110) plane were subsequently added. The crystal planes were carefully polished and two lithium-niobate (LiNbO 3 ) transducers (Z-cut for longitudinal waves and X-cut for transverse waves) were glued to opposite parallel surfaces for excitation and detection of acoustic waves. The relative ultrasound-velocity changes ∆v/v and attenuation changes ∆α were measured using an ultrasound pulseecho phase-sensitive detection technique [9,26] in pulsed magnetic fields up to 38 T (test pulses up to 56 T) at temperatures ranging from 1.35 to 30 K. Excitation frequencies were varied from 27 to 100 MHz with pulse durations ranging from 50 − 200 ns. Strain-induction coupling, i.e., the Alpher-Rubin effect [2], may be safely neglected at the used frequencies as the large magnetoresistance in NbP even at moderate magnetic fields (µ 0 H > 1 T) prevents from a strong skin effect.

III. RESULTS
The change of sound velocity ∆v/v and the change of sound attenuation ∆α vs magnetic field at T = 1.35 K are shown for different acoustic modes in Fig. 1. Here, ∆v/v refers to the change compared to the sound velocity at zero magnetic field v = C eff /ρ, where C eff is the effective elastic constant governing the respective mode [27] and ρ is the mass density (ρ = 6.52 g cm −3 for NbP [25]). ∆v/v shows pronounced QO with high harmonic content, whereas dominant frequencies and size of the oscillation amplitudes strongly vary between the modes. Strikingly, the QO amplitude in the C 66 mode is smaller by a factor of ≈ 20 compared to the other modes, where for the last few LL changes in v by more than one part in a thousand are observed. ∆α exhibits QO with a characteristic spike-like shape, also varying in terms of amplitude and dominant frequencies depending on the mode. We recall that the physical mechanism responsible for the QO in ultrasound attenuation, which are commonly termed as 'giant QO' [1,8], is not related to the Landau tubes passing through the extremal parts of the Fermi surface as in the de Haas-van Alphen (dHvA)-type oscillations. Instead, spikes in ∆α occur when the Landau tubes pass through the Fermi-surface section, where the component of the Fermi velocity parallel to q is equal to the phase velocity of sound [1, 3,8,28]. This resonance condition is the reason for the spike-like shape, as it is only fulfilled for particular values of the wavevector in contrast to the contribution of many wavevectors in the dHvA-type oscillations. Notably, the resonant Fermi-surface orbits can differ substantially from the extremal orbits, especially when q ⊥ H. Hence, the position of the observed spikes in ∆α do not necessarily coincide with the LL peaks in ∆v/v. Above 30 T, all electrons and holes are confined to their lowest LL; and v(H) and α(H) exhibit a steady slope in the investigated field (measured up to 56 T for C 44 ) and temperature range, showing no signatures for correlationdriven charge instabilities. Such correlation-driven phase transitions, e.g., a charge density wave, would manifest in a slope change of ∆v/v and a peak in ∆α [29], and have been predicted to occur in the extreme quantum limit of Weyl semimetals [30,31]. Notably, there have been observations of indicative features in the extreme quantum limit in the electrical resistivity and in the sound velocity and attenuation in the related compound TaAs [32,33]. However, in case of pristine NbP the interaction strength presumably is too feeble as to allow for experimental access to these energy scales within our achievable field and temperature range.
A. Quantum oscillations in the velocity of sound

Frequency analysis and orbit assignment
To analyze the QO in the ultrasound velocity, −∆v/v is plotted against 1/H (Fig. 2). The ultrasound velocity, just as any thermodynamic property of a material, exhibits singularities upon increasing magnetic field whenever a cyclotron orbit corresponding to a LL is exactly equal to an extremal orbit of the Fermi-surface sheet perpendicular to the applied H. According to the Onsager relation [1], these singularities are periodic in 1/H with the frequency F = ( /2πe)A ext , where A ext is the area enclosed by the corresponding extremal orbit, is the reduced Planck constant and e the electron charge. Plotting LL number vs 1/H, F can then be extracted using a linear fit [see Fig. 2 For a maximum orbit, −∆v/v will increase with (1/H) −1/2 approaching a LL singularity from a lower field and then decrease steeply, once the area of the corresponding cyclotron orbit exceeds that of the maximum orbit [28]. Accordingly, for a minimum orbit these slopes are reversed and the steep rise appears on the low-field side of the LL peak. If smearing due to finite temperature and electron scattering is sufficiently suppressed, the QO retain a high harmonic content and approach a sawtooth-like shape. The asymmetry of the individual LL peaks then allows for identifying whether the corresponding peak is arising from a maximum or minimum orbit of the Fermi surface.
Clearly, the dominant frequency of 30.89 T in C 11 and C 33 (also very well distinguishable in the (C 11 − C 12 )/2 mode) stems from a maximum orbit [most apparent for the last three LL, see Fig. 2(b)]. It is also the most pronounced frequency in the SdH oscillations in magnetoresistance [ Fig. 2(b) top], whose shape resembles that of the C 11 mode. As assigned in Ref. 18 based on DFT calculations and further indicated by comparing experimental and calculated strain dependences [23], this frequency is likely stemming from the α 1 rather than the γ 1 orbit [hereafter, we use the same labeling for the extremal orbits of NbP as in these Refs., see Fig. 2(a)]. The α 1 oscillation is much less pronounced in C 44 [see Fig. 2(c)], allowing for a clear identification of the 14.74 T oscillation as a minimum orbit, assigned to β 1 . After having identified the LL peaks for α 1 and β 1 , the remaining peaks in the high-field range might be assigned to the γ 1 orbit and possibly also the δ 1 orbit [see Fig. 2(e)]. The assignment to δ 1 is thereby rather speculative; the second peak at approx. 0.06 T −1 might also stem from the last LL of δ 2 . At low fields, a 0.9 T oscillation with minimum-orbit characteristics is visible in C 44 , assigned to β 2 [ Fig. 2(d)]. Furthermore, by applying a low-pass Fourier filter to C 11 an oscillation of 6.81 T is singled out, which was also identified in the Fourier spectra from previous QO studies [18,22,23] and assigned to the α 2 orbit [ Fig. 2(f)]. The extracted frequencies are summarized in Table I. We note that we did not observe additional QO patterns predicted to occur in Weyl semimetals when the Fermi level is near the Weyl points [5].

Lifshitz-Kosevich fit
The actual shape of the QO in ∆v/v can be described by a Fourier series taking finite-temperature smearing of the Fermi-Dirac distribution and LL broadening due to electron scattering into account. After Lifshitz and Kosevich [1], the oscillatory part of ∆v/v for a single QO frequency without spin degeneracy holds where m c denotes the effective cyclotron mass, V the real space volume, A ext " is the curvature of the Fermi surface at the extremal orbit and ϕ is the phase factor. The ±π/4 phase shift accounts for whether the orbit is maximum (−) or minimum (+). Damping of the QO due to thermal smearing of the Fermi distribution is accounted for by the factor [1] Damping due to electron scattering is taken into account by the Dingle damping factor [1] where T D is the Dingle temperature and µ c is the mobility of an electron exerting a cyclotron motion in an applied magnetic field (not to be confused with the zerofield transport mobility, which, depending on the current direction, can significantly differ from µ c in case of a large band anisotropy [35]). The β 1 , β 2 , and α 1 oscillations were clearly distinguishable in C 44 and C 33 , respectively, ( d ) and could be approximated using the first 20 harmonics of Eq. (1). From fits to the QO for different temperatures (Fig. 3), the damping factors R D and R T could be extracted, allowing for the determination of m c , ϕ, µ c , and T D (summarized in Table I). The fitting procedure was performed globally for all temperatures with the shared parameters F (fixed), m c , ϕ, and µ c , and an independent amplitude prefactor. We note that the direct fitting of the naturally filtered QO yields a greater reliability for the m c values compared to the analysis of Fourier spectra, as there the field-dependent amplitude damping usually leads to a systematic underestimation  (7) - Min The calculated frequencies were obtained from density functional theory in our previous study [23]. b Ξ 1 has been estimated with Eq. (5)

Discussion of the phase factor
The phase factors extracted from fitting Eq. (1) to the ∆v/v data are around 0.5 for the extremal orbits α 2 and β 2 on the electron pocket E2, and vary from 0.27 to 0.20 for the orbits α 1 , β 1 and γ 1 on E1. According to recent theory works by Alexandradinata et al. [37,38], the phase factor generally consists of three contributions where ϕ M is the Maslov correction (ϕ M = 1/2 for orbits that are compressible to a circle, which is the case for all orbits in NbP), ϕ B is the geometric phase, i.e., Berry phase [39], that an electron acquires upon encircling the orbit in reciprocal space, and ϕ d is the dynamic phase factor which accounts for the generalized Zeeman interaction of the intrinsic and orbital magnetic moment. The main interest in analyzing the phase contributions lies in the extraction of ϕ B , as it potentially allows to identify topologically non-trivial bands, such as Weyl or Dirac bands [39]. Indeed, under certain symmetry constraints (for details, see Refs. 37 and 38) ϕ d vanishes or can only take quantized values ±1/2, which then allows to draw conclusions about ϕ B . As all orbits in NbP for H c can be mapped onto themselves in k space upon applying a mirror operation (mirror planes k x = 0 or k y = 0, see Ref. 34), they belong to the classification (II-A, u = 1, s = 0) of Tab. I in Ref. 37, and ϕ d can be either 0 or 1/2 depending on details of the band structure. Hence, at first glance it may seem that a deviation of ϕ from 0 or 1/2 can be regarded as a signature of a non-zero ϕ B . However, it was shown by Klotz et al. [18] that the Fermisurface pockets in NbP intersecting with the Weyl bands, E1 and H1, always encompass a pair of Weyl points and should thus exhibit a trivial phase shift of ϕ B = 1 or 0. Hence, the extracted phase factors of α 1 , β 1 and γ 1 are at odds with the possible values predicted by theory. It is rather speculative why this is the case, the reason might be slight misalignment of the magnetic field, wrong orbit assignment or, more generally, inaccuracy of the DFT calculations, although the latter two are highly improbable given the otherwise good agreement. The extracted ϕ of E2 are not contradicting theory, but are also not particularly informative regarding the topological nature of the bands.

Extraction of the deformation potentials
Comparing the amplitudes of the same orbit for different modes, the ratio of the C −1 ii (dF/dε i ) 2 values can be extracted. With the known elastic constants from our previous study [23], the ratio of the effective deformation potentials can then be calculated via [2] The amplitude ratios for the individual orbits have been extracted by selecting well distinguishable LL peaks (near the quantum limit) and divide their top-to-bottom heights. In case there was no separate LL peak, as for example for the β 1 orbit in C 11 and C 33 , the height was estimated via fitting of two Lorentzian functions with fixed centers (Fig. 4), whereas the center positions were extracted from comparison with other modes (see Fig. 2). The resulting deformation potentials w.r.t. Ξ 1 are summarized in Table I. They are strongly anisotropic -measurable Ξ values vary by up to a factor of ≈ 8 depending on the direction of strain -which reflects the anisotropy of the electronic bands in NbP (see DFT calculations in Refs. 34 and 40). In contrast to the isotropic behavior in conventional metals, the electron-phonon scattering in NbP [and transferably other (Weyl) semimetals with anisotropic bands] is highly selective.
With the ∂F/∂ε 1 values gathered from Ref. 23, Ξ 1 can be estimated via Eq. (5) to be 2.1 eV (2.5 eV) for α 1 and 1.4 eV (2.2 eV) for β 1 taking experimentally (calculated) values, respectively. For β 1 , this results in an effective deformation potential of 9 eV (14 eV) for shear strain along c. This potential is among the highest reported values [10,41,42] and illustrates how electrons in the narrow part of the electron pocket are extremely susceptible to interaction with phonon modes corresponding to such shear strain. We note that upon applying strain along an axis perpendicular to the c axis, the breaking of the rotational symmetry leads to a degeneracy lifting of the Fermi pockets and ∂F/∂ε 1 actually splits into a positive and a negative branch [23]. As in Eq. (1) the sign of ∂F/∂ε 1 is canceled due to the square, we took the average of the absolute values in order to estimate Ξ 1 . The 'giant QO' in ∆α are less straightforward to analyze, as the position of the resonant orbits in reciprocal space is rather complicated to determine for each corresponding phonon mode. If plotted against 1/H [ Fig. 5(a)], two periodic series of spikes are very well distinguishable, labeled as F 1 and F 2 . The Onsager relation is valid for the 'giant QO' as well; linear fits to the spike positions vs LL number yield F 1 = 29.8 T and F 2 = 14.5 T [ Fig. 5(b)]. The areas enclosed by the resonant orbits are thus close to those of α 1 and β 1 . A puzzling feature is the observation of the same frequencies in two modes with perpendicular q, e.g., F 1 in both C 11 and C 33 . This observation might be explained by the peculiar shape of the Fermi surface in NbP, where fourfold degenerate sickle-like pockets are located near the edges of the first Brillouin zone. In this particular case, the resonant condition might be fulfilled for the same orbit for elastic waves propagating both along a and c.
In contrast to the QO in ∆v/v, the exact shape of the spikes in ∆α is rather difficult to fit. Each δ function corresponding to a spike must be convoluted with various distribution functions accounting for the effects of finite temperature and electron scattering [1]. In our case, this did not seem viable as multiple frequencies superimpose each other and similar information on the electronic properties has already been extracted from the QO in ∆v/v, where also the signal-to-noise ratio was more favorable. Nevertheless, the slight asymmetry of the spikes can be attributed to an indirect effect of electron scattering, where the smearing of the LL relaxes the resonance condition [1]. The spikes of F 1 and F 2 are broader towards the low-field side [see Fig. 1(b)], which is indicative of a convex curvature of the Fermi surface at the resonant orbit (A" < 0).

IV. SUMMARY
In summary, we studied the QO in ultrasound velocity and attenuation in NbP in pulsed magnetic fields. Thereby, fields with H c beyond the quantum limit were applied. We compared the QO for several acoustic modes, revealing significant differences as to which orbits are dominant. By extracting the amplitudes of the QO in the ultrasound velocity, the anisotropy of the deformation potentials has been determined for several extremal orbits. Thereby, a large deformation potential of approximately 9 eV for the minimum orbit β 1 under shear strain along the c axis has been revealed, suggesting that electrons in this part of the Fermi surface are very susceptible to interactions with the phonon modes corresponding to C 44 . Furthermore, the high harmonic content of the QO and the large field range allowed for a more reliable determination of the frequencies, effective cyclotron masses, and mobilities as was previously achieved by means of Fourier analysis. On a side note, we did not find any signatures for correlated electron states in the quantum limit of (pristine) NbP.