Effect of uniaxial stress on the electronic band structure of NbP

The Weyl semimetal NbP exhibits a very small Fermi surface consisting of two electron and two hole pockets, whose fourfold degeneracy in $k$ space is tied to the rotational symmetry of the underlying tetragonal crystal lattice. By applying uniaxial stress, the crystal symmetry can be reduced, which successively leads to a degeneracy lifting of the Fermi-surface pockets. This is reflected by a splitting of the Shubnikov-de Haas frequencies when the magnetic field is aligned along the $c$ axis of the tetragonal lattice. In this study, we present the measurement of Shubnikov-de Haas oscillations of single-crystalline NbP samples under uniaxial tension, combined with state-of-the-art calculations of the electronic band structure. Our results show qualitative agreement between calculated and experimentally determined Shubnikov-de Haas frequencies, demonstrating the robustness of the band-structure calculations upon introducing strain. Furthermore, we predict a significant shift of the Weyl points with increasing uniaxial tension, allowing for an effective tuning to the Fermi level at only 0.8% of strain along the $a$ axis.


I. INTRODUCTION
The symmetry and topology of electronic band structures in crystalline materials is inherently connected to the symmetry of the underlying crystal lattice. A powerful experimental tool for the exploration of the symmetry dependence of electronic properties is the application of uniaxial stress, as it allows to reduce the symmetry of the crystal lattice and thus provides a rather different kind of information than application of hydrostatic pressure. In recent years, the interest in topological properties of electronic band structures surged, which also aroused the interest in studying the uniaxial strain response of materials with remarkable topological band-structure features, using both theoretical 1-12 and experimental [13][14][15][16][17] methods.
One of those materials is the semimetallic transitionmetal monopnictide NbP, which has recently stimulated numerous studies [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36] . It exhibits a very small Fermi sea consisting of two electron and two hole pockets 18,19 , whose fourfold degeneracy in k space is tied the rotational symmetry of the underlying tetragonal lattice. Cyclotron masses in the order of 5-15% of the free electron mass m 0 and large charge-carrier mobilities in the order of 10 6 − 10 7 cm 2 V −1 s −120 give rise to an extremely large magnetoresistance (MR) of up to 10 6 % at 1.8 K and 9 T in the cleanest NbP samples 20 , superimposed by pronounced Shubnikov-de Haas (SdH) oscillations. The small and highly degenerate pockets and the accompanied experimental access via SdH oscillations render NbP a promising platform for studying the symmetry-breaking of the electronic bands induced by uniaxial stress. As a member of the TaAs family, NbP is predicted 37,38 to exhibit 12 pairs of Weyl points, i.e., crossings of valence and conduction bands with nearly linear energy-momentum dispersion which are topologi-cally protected by the broken inversion symmetry of the crystal lattice. They can be classified into two fourand eightfold energy-degenerate groups, whose degeneracy is also tied to the rotational symmetry of the crystal lattice. The existence of Weyl points near the Fermi level is thought to account for the ultrahigh mobilities observed in NbP 20,21 . The Fermi surface of NbP has been studied 18 theoretically by means of ab-initio density functional theory (DFT) calculations as well as experimentally by analyzing the de Haas-van Alphen oscillations found in the magnetic-torque signal. Measuring magnetic quantum oscillations by means of SdH experiments, also the evolution of the Fermi surface of NbP with applying hydrostatic pressure has been probed 22 . These experiments showed a relatively stable behavior of the SdH frequencies but significant changes of their amplitudes, indicating slight deformation of the electron and hole pockets. The effect of uniaxial stress on the electronic band structure of NbP has not yet been investigated, however, it is particularly interesting as the stress-induced breaking of the fourfold rotational symmetry is expected to lift the degeneracies of the Fermi pockets and the Weyl points.
In this paper, we report on the effect of uniaxial stress on the Fermi surface of NbP, determined via SdH oscillations in the MR of two single-crystalline NbP samples (Sample 1 and 2) subjected to direct application of uniaxial tension in a three-piezo stack apparatus (CS100, Razorbill Instruments Ltd). For Sample 1, the magnetic field H was applied along a high-symmetry direction, resulting in a splitting of the SdH frequencies upon applying uniaxial stress due to the breaking of the fourfold rotational symmetry and thus lifting of the Fermi-pocket degeneracy. In contrast, the SdH frequencies in Sample 2, for which H has been applied along a low-symmetry direction, did not exhibit splitting and just shifted upon in-creasing strain, as in that direction no crystalline symmetry is broken. Both trends are consistent and reversible with varying ε 1 , and are in good agreement with our DFT calculations. NbP has a non-centrosymmetric tetragonal crystal lattice (space group I4 1 md, No. 109) with the lattice parameters a = b = 3.3324(2)Å and c = 11.13705(7)Å 39 . The crystal structure is shown in Fig. 1(a).
Single crystals of NbP were grown using chemical vapor transport reactions and oriented using X-ray diffraction. We cut two ∼3 mm-long bar-shaped samples along the a axis from the same single-crystalline piece. Sample 1 is 280(10) µm wide and 75(5) µm thick, with the shortest side along the c axis. Sample 2 is 310(10) µm wide and 140(5) µm thick, with the shortest side along the b axis.
In order to determine the elastic response of NbP to applied uniaxial stress, the elastic constants have been measured on a third sample of the same batch using an ultrasound pulse-echo technique at liquid-nitrogen temperature. For that, the main crystal planes as well as the [110]-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 elas-tic waves. Ultrasound propagation was measured at 27 − 100 MHz with sound-pulse durations ranging from 50 − 200 ns.
For studying the SdH oscillations, electrical-transport measurements have been carried out at temperatures between 2 and 10 K in a Quantum Design Physical Property Measurement system (PPMS) equipped with a 9 T magnet (results shown only for T = 2 K). In all measurements, current has been applied along the a axis, i.e., parallel to the uniaxial stress σ 1 . The resistivity ρ xx and Hall resistivity ρ xy have been measured in a 4-probe configuration using Stanford SR-830 lock-in amplifiers with a reference frequency of 117.23 Hz and an excitation current of 1 mA. Electrical contacts were fabricated using 25 µm Pt wire and silver paint.
After characterizing the samples in the unstrained state, they were subsequently mounted to the strain cell using titanium plates and spacers with the same height as the samples, glued from top and bottom with insulating expoxy (Stycast 2850FT, Loctite) and tightened with brass screws [see Fig. 1(c)]. The length of the strained part of the sample is approximately 1-1.5 mm. The displacement is measured via a pre-calibrated capacitor, using an LCR-Meter (E4980 AL, Keysight) at 27.7 kHz. The peculiarities of this strain cell are extensively covered in Ref. 40. Notably, the size and homogeneity of the effective strain sensitively depend on a number of parameters, most importantly the sample's spring constant k s as well as the thickness and stiffness of the epoxy. The setup with fixation from both bottom and top of the sample has been reported 40 to exhibit a reasonable degree of homogeneity in the center part of the strained sample, which is where we probed the voltage. At 2 K, the maximum capacitance range in our strain cell was 880-945 fF for Sample 1 and 925-950 fF for Sample 2; in that range the displacement (few µm) as a function of capacitance can be approximated as linear. Overreading of the displacement due to device deformation is taken into account with a correction factor (1 + G · k s /k τ ) −1 , where k τ is the torsional stiffness of the strain cell and G accounts for geometric details. Strain losses due to the finite shear stiffness of the epoxy are accounted for with another correction factor, which is taken to be 80% as estimated in previous works on this type of strain cell and same epoxy 41,42 . During the mounting procedure, heat curing of the epoxy and thus thermal expansion of the piezoelectric stacks and other parts of the chassis lead to a deviation of the sample from the zero-strain state with respect to the calibration when cooled back to room temperature. We account for this by subtracting a constant offset, such that a linear extrapolation of the SdH frequencies under strain intersects with the SdH frequency in the unstrained state. This yielded a moderate compression of −0.17 % for Sample 1 and −0.1 % for Sample 2 at room temperature without having any voltage applied to the piezoelectric stacks.
For our ab-initio investigations of the electronic band structure, we employ DFT as implemented in the VASP package 43 . Here, plane waves with pseudopotentials are used as a basis set and the generalized-gradient approximation (GGA) 44 is utilized for the description of the exchange-correlation potential. For calculating fine meshes in the Brillouin zone we create a tight-binding Hamiltonian from Wannier functions using the Wan-nier90 package 45 with initial projections to the d orbitals of Nb and to the p orbitals of P. For the evaluation of the extremal orbits a k mesh with a resolution better than 0.004Å −1 was used.

III. ELASTIC PROPERTIES
To determine the elastic response for a given applied uniaxial stress and to be able to estimate the effectively applied strain in the samples, the elastic stiffness tensor C ij of NbP needs to be considered. Using Voigt notation (xx → 1, yy → 2, zz → 3, xz → 4, yz → 5, xy → 6), the sample's strain components ε j are related to the applied stress components via We determined all elastic constants except for C 13 via measurement of the speed of sound v = C eff /ρ for longitudinal and transverse modes along the main axes as well as along [110], ρ being the mass density (ρ = 6.52 g cm −3 for NbP 39 ). The results are shown in Tab. I and show good agreement with the values calculated from DFT 46 .
Applying uniaxial tension along a (σ 1 > 0, σ i = 0 for i = 1) results in strain ε 1 = σ 1 /E along the same axis, where the Young's modulus E in that configuration holds S 11 being the first component of the elastic compliance tensor S ji . Given the good agreement between experimentally determined and theoretical elastic constants, we used the calculated value 46 C 13 = 130.8 GPa to estimate the Young's modulus, yielding E = 210(60) GPa (E = 251 GPa using only the calculated values). With k s = EA/l, A being the sample cross section area and l its length, the sample spring constant yields (4.4 ± 1.4) × 10 6 Nm −1 for Sample 1 and (6.3 ± 1.8) × 10 6 Nm −1 for Sample 2. This results in an about 15% and 25% too large displacement value due to cell deformation in Sample 1 and Sample 2, respectively.
Naturally, tension along a results in compression along b with the Poisson ratio and along c, respectively, with This yields ν 21 = 0.52 (16) and ν 31 = 0.23(11) (ν 21 = 0.36 and ν 31 = 0.29 from calculated values). In our DFT calculations, the latter were used. The effect of uniaxial stress σ 1 on the crystal symmetry is illustrated in Fig. 1(b): the fourfold rotational symmetry, more precisely the screw operation 4 1 parallel to the c axis of the tetragonal lattice, is broken and the crystal symmetry reduces to an orthorhombic lattice. In our setup, only the displacement ∆l along the a axis is measured, allowing for the determination of ε 1 = ∆l/l.

IV. RESULTS
The MR ρ xx (H) and Hall resistivity ρ xy (H) of both samples are shown at different (tensile) strains ε 1 > 0 in Fig. 2, whereby H has been applied along c for Sample 1 and along b for Sample 2. Both samples exhibit a large MR characteristic for NbP 20 , superimposed by SdH oscillations. In agreement with previous reports, the amplitude of the SdH oscillations strongly varies depending on the H direction and is most pronounced for H along c, which is also apparent comparing the MR in Sample The strain response of ρ xx and ρ xy was completely elastic, i.e., no hysteresis was observed upon increasing or decreasing the stress and all curves could be reproduced by applying the same stress again. This provides an ideal playground for in-situ manipulation of the band structure, without introducing additional disorder or facing the irreversibility of high-pressure experiments. Notably, both MR and Hall resistivity do not exhibit significant changes apart from the superimposed oscillations, which shows that the deformation of the Fermi surface does not have a strong effect on the general charge-carrier properties, e.g., density and mobility, but rather leads to a redistribution of the carriers in k space. This is in contrast to the effect of hydrostatic pressure on NbP 22 , where both MR and Hall resistivity change considerably.
By Fourier transformation ofρ xx /ρ xx (H −1 ), a characteristic spectrum is revealed. Via the Onsager relation 47 F = (Φ 0 /2π 2 )A ext , each SdH frequency F is related to an extremal orbit of the Fermi surface plane perpendicular to the applied H, where A ext is the area enclosed by this orbit and Φ 0 is the magnetic flux quantum. The SdH spectra of the two samples in the unstrained state are shown in Fig. 3(a), for H rotated from the c to the b axis. Noticeable differences in the intermediate range between 20 and 70 • occur, presumably due to a bigger sensitivity to misalignments in that range. However, for H along b and c the spectrum of the two samples is in good agreement and is roughly stable upon misalignments of ±5 • . This is expected, as both samples are cut from the same piece and thus the Fermi level should be nearly identical. The H direction in our strain setup is fixed, thus, the matching of the SdH spectra is a necessary check to justify the comparison of the two samples.
The Fermi surface of NbP is shown in Fig. 8(a). There are a number of extremal orbits [schematically shown in Fig. 3(b)], which we labeled as in Ref. 18, resulting in a complex spectrum. The two hole pockets, H1 and H2, exhibit only maximum orbits. For H along c, the orbit is termed as δ 1 (δ 2 ). For H along b or a, the fourfold rotational symmetry yields two maximum orbits per hole pocket, termed as δ 1 and δ 1 (δ 2 and δ 2 ). The two electron pockets, E1 and E2, exhibit a number of both maximum and minimum orbits. For H along c, E1 (E2) exhibits two maxima, α 1 (α 2 ) and γ 1 (γ 2 ), and a minimum orbit β 1 [β 2 yields a SdH frequency smaller than 1 T and is not shown in Fig. 3(b)]. For H along b or a, there are three possible orbits for E1 (E2), termed α 1 (α 2 ), α 1 (α 2 ), and γ 1 (γ 2 ). Variation of the Fermi level to 11 meV above the charge-neutral point resulted in the best matching between calculated and experimental SdH frequencies. However, we could not attribute each orbit to a peak in the experimental spectra. Also, we observe experimental peaks in a range where DFT yields no fundamental frequencies, some of which can be attributed to higher harmonics. Some orbits yield very similar SdH frequencies, as for example α 1 (32.8 T) and γ 1 (31.1 T), which falls below our experimental resolution. In case of two calculated frequencies lying close to a peak in the experimental spectrum, we chose the orbit with the highest congruency of the calculated and experimental dF/dε 1 . This procedure led to the same assignments of experimental and calculated values as in Refs. 18  F r e q u e n c y ( T ) ther challenges of disentangling the rich spectrum of NbP are thoroughly covered in Ref. 18. In the following, we focus on Fourier peaks which are distinct and could be attributed to a theoretical orbit with good agreement, marked in Fig. 3(a) and summarized in Table II. The Fourier spectra for H c (Sample 1) for different strains are shown in Fig. 4(a). A clear splitting of the α 1 and β 1 peaks into an increasing and a decreasing branch upon increasing ε 1 is observed. The evolution of α 2 is less pronounced, however, signatures of the α 2 peak found for some of the strain values are indicative of a splitting as well. As the increasing branch of α 2 is crossing the decreasing branch of β 1 , and the amplitude of β 1 surpasses that of α 2 , the increasing branch of the α 2 peak cannot be fully traced. The splitting of the peaks is consistent with our DFT calculations and is a strong evidence for the lifting of the degeneracy of the Fermi pockets via breaking the fourfold rotational symmetry of the tetragonal crystal lattice. As illustrated in Figs. 4(b), each fourfold degenerate orbit splits into two shrinking and two expanding orbits, resulting in two branches of SdH frequencies [see also Fig. 8(d)].
In contrast, the Fourier peaks for H b (Sample 2) do not exhibit a splitting or broadening and just shift upon increasing strain [ Fig. 5(a)]. This qualitative behavior is congruent with our DFT calculations and serves as a control experiment, confirming that the splitting of the SdH frequencies is tied to the orientation of H along the 4 1 screw axis of the crystal lattice.
The experimental and calculated SdH frequencies for the attributed orbits are plotted versus ε 1 in Fig. 6. All calculated trends are qualitatively reproduced by the experimental values, demonstrating that the strain evolution of the Fermi surface of NbP is well described by DFT. We note, however, that the absolute values of the majority of experimentally determined frequencies differ from the calculated values by ∼10%. This deviation is of similar magnitude as in the previous reports comparing quantum oscillations and DFT results for NbP 18,22 . The derivatives dF/dε 1 are extracted by linear fits of F (ε 1 ) and are given next to the DFT results in Table II. The sign and order of magnitude of the dF/dε 1 values are in good agreement between theory and experiment for most orbits, whereas matching in terms of absolute values for some orbits, e.g., α 1 and α 1 differs by more than 100%. Erroneous extraction of the experimental values might be caused by the smearing of the peaks due to the field-dependent damping of the SdH amplitude. The cyclotron masses m c = 2 /2π · ∂A/∂E of the corresponding orbits of area A = A ext at the Fermi energy E F are extracted from the T -dependent damping of the Fourier peaks via 47 where R T is the amplitude factor, the reduced Planck constant, e the electron charge, and k B the Boltzmann constant. 1/H is taken as the average (1/H) avg of the Fourier window ∆(1/H). We extracted m c for several orbits (Fig. 7) from the Fourier spectra at T = 2 − 10 K. Close to zero strain, the experimental values for α 1 and and similarly deviate from the higher calculated values. Cyclotron masses extracted from Fourier spectra, especially with frequencies below 100 T, are usually underestimated by up to 50-60 % 48,49 due to the strong field dependence of the SdH amplitude. The extracted mass heavily depends on ∆(1/H) and (1/H) avg , whereas decreasing ∆(1/H) generally yields m c values closer to the real one. However, as the sampling interval of the discrete Fourier spectrum is (∆(1/H)) −1 , decreasing ∆(1/H) drastically reduces the resolution and makes identifying distinct peaks for H c (Sample 1) impossible. Additional complication arises from the overlapping of the Fourier peaks. Therefore, the large error bars in Fig. 7 are estimated by extracting m c for different ∆(1/H) and (1/H) avg as well as fitting Eq. (5) to Fourier spectra of simulated data with similar frequencies and cyclotron masses. Notably, envelope extraction via inverse Fourier transformation of a frequency window (with rounded corners 49 ) enclosing a Fourier peak of the simulated data also did not result in reliable determination of m c . Due to the large uncertainty, we were not able to resolve the splitting of m c predicted by DFT; the extracted values can only be considered as a rough estimate. The cyclotron masses for H b (Sample 2) are presumably slightly more reli- able, as the peaks are more distinguishable and a smaller ∆(1/H) window could be used.
The overall qualitative (and in some parts quantitative) agreement of the DFT calculations with the SdH data might as well allow us to draw conclusions about the Weyl points in NbP under uniaxial stress. The 12 pairs of Weyl points in NbP fall in two groups: four pairs in the k z = 0 plane (WP1) and eight pairs in the k z ≈ ±π/c plane (WP2). Upon increasing ε 1 , the four (eight) WP1 (WP2) pairs split into two types of two (four) pairs each, labelled "high" and "low" [see Figs. 8(b),(e)], shifting up and down in energy, respectively. Notably, the energy shift is quite significant, with tens of meV (or hundreds of K in terms of E/k B ) at a few percent of strain [ Fig. 8(c),(f)]. Predicted by our calculations, the WP1 high pairs should be situated at the Fermi level at ε 1 ≈ 1.6 %, and WP2 low pairs at ε 1 ≈ 0.8 %, respectively. In contrast to the application of hydrostatic pressure 22 , where several GPa are required to tune WP2 a few meV closer to the Fermi level and WP1 does not exhibit any notable shift, uniaxial stress might allow for effective tuning of the Weyl points without having to deal with additional disorder as is the case with chemical doping. As epitaxial thin films of NbP have already been grown 50 , realization of strained thin films might render a promising opportunity towards "Weyl-tronic" devices.

V. SUMMARY
In summary, we studied the response of SdH oscillations in NbP to direct application of uniaxial stress and combined the experimental results with DFT calculations. By applying stress along one of the axes of the square ab plane of the tetragonal crystal lattice, the fourfold rotational symmetry is broken. This is reflected by a splitting of the peaks in the Fourier spectra of the SdH oscillations for H along the c axis. When H is applied along one of the equal axes, the symmetry of the orbits on the Fermi surface remains unchanged, thus only shifting of the Fourier peaks is observed. Our observations are qualitatively congruent with the DFT calculations, demonstrating their robustness for NbP in the strained state. Furthermore, we predict a significant energy shifting of the Weyl points, whereas increasing tensile strain of one percent leads to a shift of tens of meV. This might be an effective way to tune the Weyl points to the Fermi level in order to realize Weyl-based devices.