Precise Fermi level engineering in a topological Weyl semimetal via fast ion implantation

,

The precise controllability of the Fermi level and carrier doping in a solid-state material holds paramount significance in condensed matter physics and materials science.In microelectronics such as high-electron-mobility transistors, the ability to fine-tune the Fermi level is crucial to control the operation thresholds and high-frequency responses [1][2][3].In strongly correlated systems such as high-temperature cuprate superconductors, achieving optimal hole doping is essential for enhancing the superconducting critical temperature [4][5][6].As for energy harvesting applications such as thermoelectrics, a well-adjusted Fermi level strikes a delicate balance between electrical conductivity and thermopower, which leads to an optimized power factor [7][8][9][10].A particularly noteworthy utilization for Fermi level fine-tuning lies in topological materials, where only a precisely positioned Fermi level can unveil the topological properties.For instance, in topological insulators, dissipationless electronic states such as quantum anomalous Hall states can only manifest when the Fermi level is located precisely within the small surface bandgap.This can only be achieved through a strategic blend of chemical doping [11][12][13] and electrostatic gating in thin heterostructures [14][15][16][17].For topological Weyl semimetals (WSM), exotic phenomena like large charge-to-spin interconversion [18,19], strong higher-order photoresponse [20,21], anomalous Nernst effect [22,23], and quantized thermoelectric Hall effect [24,25] emerge due to a divergent Berry curvature, which can only happen when the Fermi level is tuned to be exactly at the Weyl nodes.However, unlike in lower dimensional materials, in which the Fermi level can be continuously fine-tuned through a variety of accessible techniques like electrostatic and ionic gating, achieving this level of control in three-dimensional bulk materials has presented a longstanding challenge.Chemical doping can only be performed during the synthesis and offers limited precision, which does not meet the stringent requirement for topological semimetals.The electronic states in topological semimetals, such as from the Weyl and Dirac bands, have a linear dispersion with zero density-of-states at the topological singularities, making the Fermi level location extremely sensitive to the carrier concentration [26,27].Therefore, there is an urgent need to attain enhanced fine-tunability of topological semimetals for more practical energy and information applications.
In this work, we use high-energy, accelerator-based hydrogen implantation to achieve ultra-fine carrier doping in WSM of tantalum phosphide (TaP) [24].Although hydrogenation and low-energy ion implantation have been used to tune materials properties [28][29][30][31][32][33], there are a few key distinct features presented in the current setup.Foremost, the fine controllability of the advanced accelerator technique enables an ultrafine tuning of the Fermi level down to the meV regime, far beyond the previous approaches.In addition, the high energy beam also enables the doping of a bulk crystal beyond a thin film with tens of nanometers [30,31,34].Moreover, thanks to the rapid development of the accelerator technology with high precision control of the accelerator beam energy and flux of ions, we were able to achieve a priori doping planning through a combination of DFT and Monte Carlo calculations, which eventually shows quantitative agreement with actual experiments.This sets up an example for ultrafine doping planning experiments.The Weyl nodes are preserved The doping experiment uses the Tandem accelerator located at the Center for Science and Technology with Accelerators and Radiation (CSTAR).The schematic used in our experiment is shown in Fig. 1a.This type of accelerator consists of two successive linear accelerators with the same power sup-ply.A cesium sputtering source is used as the ion source to generate negative ions from most of the elements with low vapor pressure (excluding noble gases).This source works by cesium sputtering a solid target, which produces negative ions that are then subjected to electron exchange with neutral cesium.Negative ions are initially introduced into the first stage (blue arrow), where they are accelerated toward the positive high-voltage (HV) terminal.Subsequently, depending on the experimental preferences, nitrogen gas can be injected into the HV terminal to strip some of the electrons from the ions, which transforms them into positive ions.As a result, the Tandem accelerator enables a remarkable capability to dope tens of different elements in the period table of elements.We have chosen negative hydrogen as the dopant, which offers several advantages in materials science due to its small atomic size, i.e., it does not significantly alter the structure of pristine materials [35].Electrochemical hydrogenation can also be used for electron doping of semiconductors with good doping regime controllability, but not high-precision doping level control [36].In addition, hydrogen can exhibit a range of exotic properties as a dopant that makes it particularly interesting for research and applications, including metal-to-insulator transitions [37], magnetostriction [38], enhanced superconductivity [39], and paraelectric-to-ferroelectric phase transitions [40].
The earlier study [24] highlights TaP as a remarkable topo- logical Weyl semimetal with an inversion symmetry-breaking crystal structure (Fig. 1b).It exhibits a non-saturating thermopower, a giant magnetoresistance, and quantized thermoelectric Hall effect.Shubnikov-de Haas oscillations (SdH) are observed in the background-subtracted MR (∆MR) data at temperatures below 25 K.The Fourier transform of ∆MR showed four pockets, with two small carrier pockets, one with a low-frequency F α = 4 T and the other with F β = 18 T indicating the presence of the W2 Weyl node, in addition to the electron pocket from the W1 Weyl node contributing to F β (Fig. 1c) [24].Introducing H − ions through irradiation is expected to shift the Fermi level closer to the Weyl node by carrier doping, as illustrated in the schematic diagram in Fig. 1d.By employing the Vienna ab initio simulation package (VASP), we conducted density functional theory (DFT) calculations.These calculations provide a qualitative insight into the energy redistribution of Weyl node energy levels upon H doping.The energy of the Weyl nodes with respect to the Fermi level (E f ) is shown in Fig. 2a for TaP and Fig. 2b for the doped TaP system (one H in Ta 16 P 16 ).Here, the horizontal axis is the index of the Weyl nodes.The Weyl nodes are ranked according to their energy, so the index itself does not carry additional information.Notably, the doped TaP system exhibits Weyl nodes that are distributed somewhat con-tinuously.Comparing the doped TaP to the pristine TaP, our calculations indicate that certain Weyl nodes in the doped system lie merely 0.5 meV away from the Fermi level, whereas, in the pristine system, they are much further at 10 meV from the Fermi level.A theoretical model and ab initio calculations estimating the correlation between hydrogen doping concentration in TaP and the shift in the Fermi level are discussed in Supplementary Information (SI).
To calculate the implanted ion energy, penetration depth within the crystal, and damage to the target material, the stopping and range of ions in matter software (SRIM/TRIM) [41] was employed.The detailed calculations are presented in the SI.Single crystals of TaP were irradiated with 20 keV H − ions with varying durations to achieve the desired doping concentration.The summary of different samples with their respective dose of irradiation and time duration is presented in Table I.
The irradiated samples named S 3m and S 20m exhibit a significantly high magnetoresistance (MR) on the order of 10 4 % and accompanied by SdH oscillations at temperatures below 15 K (see Fig. 3a, and Fig. 3f, respectively).In contrast, the S 2h sample does not display such oscillations even up to 9 T (details in SI).It is important to note that the analysis of carrier pockets based on quantum oscillations can be influenced by the choice of the background for MR.Therefore, we calculated the background-subtracted MR (∆MR) using different background subtraction methods (details in SI).Fig. 3g illustrate the ∆MR obtained with a polynomial-fitted background subtraction method, while Fig. 3c, and Fig. 3h show the corresponding fast Fourier transform (FFT) results for S 3m and S 20m , respectively.These FFT plots reveal a noticeably lower oscillation frequency compared to the pristine TaP [24].Furthermore, Fig. 3d, and Fig. 3i present the results of the Landau level (LL) fan diagram analysis which indicates the presence of more oscillations with lower frequencies.The oscillation frequencies have been repeatedly extracted using various background subtraction methods, and scatter plotted together to show agreement between methods (see Fig. 3e, and Fig. 3j).This finding supports and validates our DFT calculations concerning the continuous energy distribution of Weyl nodes in proximity to the Fermi surface compared to pristine TaP.
To find the effect of irradiation on the carrier concentration and the mobility Fig. 4, the longitudinal and transverse conductivity data were fitted with a two-band model.The details of the data extraction from these measurements are discussed in the SI.The compensation temperature increased to 160 K and 155 K for S 3m and S 20m , respectively, whereas it is 110 K for pristine TaP [24].Fig. 5 shows the normalized carrier concentration n T extracted from the simple extended Kohler's rule.The consistent temperature-dependent trend observed in our analysis, in accordance with the fitted two-band model, serves as compelling evidence for the accuracy and validity of our findings (see SI for detailed calculation).
In summary, we have addressed a novel method for ultrafine Fermi level tuning in the prototype bulk WSM TaP using high-energy hydrogen implantation facilitated by acceleratorbased techniques.Our transport measurements demonstrate a successful increment in the charge neutral point temperature through precise tuning of the Fermi level in proximity to the Weyl nodes.These experimental results are achievable by the guidance from DFT calculations.The approach presented in this study represents a highly controllable and universally applicable technique for fine-tuning the Fermi level and quantum orderings in bulk quantum materials.By enabling precise control over the Fermi level in bulk crystals, the work could open up new possibilities for exploring and manipulating the unique properties of topological WSM and other quantum materials.This advancement holds significant promise for the development of future quantum technologies and applications.(ARPES) experimental results, and the order of doping magnitude is what we are concerned about (We ignored the red lines in these figures since they are from the ab initio calculation).The nodes were approximated with upright but distorted cones, and their dispersion slopes were extracted from the available ARPES figure of each node type.With this approximated model, the constant energy (E) surfaces are ellipsoids with semi-principal axes Here, E 0 is the node's energy level, k 0 is the node's location, and d is the dispersion slope.This implies that, for each node, the energy increment from the Fermi level (E f ) to E would lead to an expansion in the volume of momentum space located below the Fermi level by The absolute value is omitted from energy differences with E 0 , resulting in negative volumes for energies below E 0 .This adjustment ensures the accurate representation of momentum space volume increments.Given the standard practice of expressing energy relative to the Fermi level, the formulas were subsequently simplified using this notation as follows: Consider energy levels E 1 and E 2 corresponding to the nodes of W1 and W2, respectively, where E 1 = 24 meV and E 2 = −40 meV.Assuming that E is close to 0 (E f ) enough such that only the nodes contribute to the total momentum space volume increment.The total increment in momentum space volume (sixteen from W1 and eight from W2) would be where d 1 and d 2 are average dispersion slopes of W1 and W2 nodes respectively.This signifies that the number of electronic states increases by In an alternative perspective, the requirement arises for TaP to undergo doping with an additional electron concentration of It can be observed from the extracted figures that Therefore, the doping electron concentration in cm −3 of TaP for the desired Fermi level shift, E, in meV is The plot depicting the equation ( 7) for E ∈ [−50, 50] meV is shown in Fig. S1c.By substituting E with E 1 , and E 2 , the model predicts the doping electron concentration of 5.3×10 19 cm −3 for W1, and doping hole concentration of 8.6×10 19 cm −3 for W2, respectively.

Ab initio Calculation of H-doped TaP
To estimate the correlation between hydrogen doping concentration in TaP and the shift in the Fermi level, ab initio calculations were done.It involves determining the count of unoccupied states situated between the initial and adjusted Fermi levels.In this Subsequently, structure relaxation was performed on the structure from the MaterialsProject to get a stable structure according to the pseudopotentials used.This stable structure was used in Quantum Espresso to calculate the self-consistent function, and non-self-consistent calculation with high k-point density grids [52].Then, we perform density of states and band structure calculation pathway.The band structure and density of the states plot are shown in Fig. S3a and Fig. S3b, respectively.The doping relation was achieved by integrating the density of states starting from the Fermi level.This allowed the determination of the number of states per unit cell that needed to be filled with the doped electron to effectuate the shift of the Fermi level (shown in the bottom of Fig. S3b).Finally, the doping electron concentration was calculated by dividing the unit cell volume.Our calculation predicts the doping electron concentration of 9.7×10 19 cm −3 for W1, and doping hole concentration of 1.4×10 20 cm −3 for W2.

SRIM and TRIM Calculation
We employed the SRIM (stopping and range of ions in matter) and TRIM (transport of ions in matter) software [41] to simulate realistic irradiation conditions and investigate the defect profile.TRIM simulations provided us with the flexibility to adjust parameters such as target type, thickness, and incident ion beam energy.Through the TRIM simulations, we obtained valuable results that revealed the depth profile generated by both the incident ions and recoils.For our study, we deliberately chose H − ions with an energy of 20 keV to irradiate the TaP crystal.This selection was driven by our specific interest in doping the ions near the surface of the sample (∼ 5000 Å).By considering the beam current and aperture size, we calculated the incident ion flux.Furthermore, we accounted for the subtraction of secondary electrons and the total duration of irradiation, enabling us to determine the number of implanted ions in each sample.In our experimental setup, we utilized a 4 mm diameter beam to irradiate samples of size 2×2 mm 2 .The beam current for S 3m was 140 nA, and for S 20m and S 2h was 170 nA.We varied the irradiation time accordingly to achieve the targeted doping concentration of ≃ 10 19 cm −3 , 10 20 cm −3 and 10 21 cm −3 for S 3m , S 20m and S 2h samples as per our mathematical and computational calculations.The total doping concentration in the sample is determined by multiplying the irradiation dose for each sample by the penetration depth.The irradiation dose is summarized in the main text (Table I).The electrical measurements of the samples were carried out using the electrical transport option (ETO) of the physical property measurement system (PPMS).We utilized the ETO with a six-probe configuration to obtain the required data, enabling simultaneous measurements of both ρ xx and ρ xy .We employed a symmetric probe configuration to minimize the influence of manually made contact misalignment.In this setup, the longitudinal resistivity ρ xx exhibited symmetry with respect to the applied magnetic field, while the transverse resistivity ρ xy displayed antisymmetry.We employed the following equations to achieve the symmetric probe: The magnetic field dependence of ρ xx and ρ xy at different temperatures is illustrated in Fig. S4.Shubnikov-de Haas (SdH) oscillations can be observed in both ρ xx and ρ xy at low temperatures, up to 15 K for S 3m , and 10 K for S 20m .However, no such oscillations are observed for S 2h sample, indicating that the presence of high disorder suppresses the quantum oscillations.Thus, our further experimental analysis primarily focuses on the S 3m , and S 20m samples as our main interest lies in precise Fermi level engineering in a topological Weyl semimetal.and transverse conductivity, σ xy , using the following equations: The field dependence of σ xx and σ xy at different temperatures are illustrated in Fig. S5a and Fig. S5b for S 3m , respectively, and Fig. S6a and Fig. S6b for S 20m , respectively.To extract the carrier concentration and mobility, we simultaneously fit the σ xx and σ xy data using two-band model defined by: Fig. S5c-d, and Fig. S6c-d show example fitting results of the curves for 2 K temperature for S 3m , and S 20m , respectively.The root mean squared relative error (RMSRE) linked to the fitting procedure at different temperatures using the two-band model is displayed in Fig. S5e-f, and Fig. S6e-f for S 3m , and S 20m , respectively.The main text includes the temperature variation of electron and hole concentration, along with their corresponding mobility, for each doping level and compares them with pristine TaP (Fig. 4).Our findings reveal that H − -ion implantation significantly impacts the Fermi surface by considerably elevating the electron-hole compensation temperature despite the mobility being reduced by an order of magnitude.

Analysis of Quantum Oscillation
Magnetoresistance (MR) data was extracted from ρ xx data, using the following relation: The quantum oscillation data, ∆MR, was obtained by subtracting background signals from the MR data.To ensure robustness, we employed three independent methods to extract and analyze ∆MR: (a) background-free curvature (Fig. S7 and Fig. S8), (b) polynomial background (partially in Fig. 3, and fully in Fig. S9 and Fig. S10), and (c) piece-wise polynomial background (Fig. S11 and Fig. S12).In (a), a second-order derivative with respect to the magnetic field B is applied to the MR data, automatically eliminating all constant, linear, and quadratic terms, which mostly contribute to background instead of oscillation, without the need for actual background selection.Notably, the fast Fourier transform (FFT) analysis unveiled distinct oscillation frequencies.After applying a standard signal filtering process through FFT against 1/B, we isolated the oscillation components from the data and determined the corresponding Landau levels (LLs).The LL index fan corresponding to the different frequencies for each fitting model has been plotted.All three methods consistently support the existence of a low-frequency carrier pocket, F α 1.87 T -3.94 T for S 3m and 0.35 T -1.14 T for S 20m .In contrast, the pristine TaP exhibits the same behavior in the range F α 2.3 T -4 T and F β 18 T -19 T [24].Our comprehensive analysis confirms the presence of new oscillation frequencies, suggesting that the carrier pockets of the Weyl node can indeed reach the desired 0 Landau level, and additional Weyl nodes appear with H − ion irradiation.

Extended Kohler's rule of magnetoresistance
The MR curves do not coincide onto a single curve when plotted against B/ρ xx (0), indicating that the MR = f (B/ρ xx (0)) rule proposed by Kohler is not applicable to S 3m and S 20m , as observed in our study.However, an intriguing observation is that all the curves in Fig. S13a-b are nearly parallel with each other.This suggests the possibility of a single temperature-dependent multiplier affecting MR (y-axis) or B/ρ xx (0) (x-axis), causing them to coincide or merge into one curve.To understand this behavior further, MR curves are normalized by a scaling procedure using the 300 K curve, and adjust the parameter n T individually for each curve [53].The outcome is clearly visible in Fig. S13c-d, where all the curves seamlessly coincide onto the 300 K curve when the data is scaled using B/ρ xx (0)n T .The corresponding values of n T for the different curves at various temperatures can be found in Fig. S14.This scaling procedure sheds light on the temperature-dependent behavior of the MR and provides valuable insights into the underlying physical mechanisms governing the observed phenomena.We have further normalized n T and compared it with the data extracted by using a two-band model from the relation [53]:

FIG. 1 .FIG. 2 .
FIG. 1. a) Tandem accelerator schematic.It includes an ion source to produce negative ions.Selected ion species are injected toward the terminal using a low-energy magnet.The ion beam is focused by a magnetic quadrupole lens and the desired ion species are directed into the beamline for the experiment.b) Crystal structure of the Weyl semimetal TaP.c) Fast Fourier transform of the longitudinal magnetoresistance of TaP reveals quantum oscillations with specific frequencies (reproduced from [24]).d) Illustration of the irradiation effect on the Fermi level of TaP.The red arrow indicates the shift of the Fermi level towards a Weyl node.

FIG. 3 .
FIG. 3. a) Variation of magnetoresistance (MR) with magnetic field (B) at different temperatures for sample S3m.b) Background-subtracted MR (∆MR) of a versus 1/B for a few temperatures that exhibit Shubnikov-de Haas (SdH) oscillations.This ∆MR is calculated using q polynomial-fitted background subtraction.c) Fast Fourier transform (FFT) of b. d) Example of a Landau level fan diagram analysis which indicates the presence of lower frequencies in b at 5 K than TaP [24].e) Scatter plot of the oscillation frequencies extracted from b using various background subtraction methods at 5 K. f)-j) are the same as a)-e), but for the S20m sample.

FIG. 4 .
FIG. 4. a) Temperature-dependent profiles of electron (ne) and hole concentrations (n h ).b) Temperature-dependent trends of electron (µe) and hole mobilities (µ h ).The top profiles serve as references for TaP.Our analysis involves the extraction of equivalent profiles from the two-band model fitting, applied to the S3m and S20m samples as annotated in the plots.

FIG. 5 .
FIG. 5. Direct comparison of the temperature dependence of the normalized adjusted parameter nT , as extracted from both the extended Kohler's rule and the two-band model for a) S3m, and b) S20m.The plots show good agreement between simple extended Kohler's rule and fitted two-band model.
FIG. S2.Convergence test on simulation's hyperparameters: a kinetic energy cutoff for wavefunction (ecutwfc), b kinetic energy cutoff for charge density (ecutrho), c Gaussian spreading for Brillouin-zone integration in metals (degauss), and d number of k-point grids on each axis.

FIG
FIG. S4.Magnetic field dependent of longitudinal resistivity, ρxx, at different temperatures up to 300 K for a S3m, b S20m, and c S 2h .d-f The same as a-c, respectively, but for ρxy.

FIG
FIG. S5.Temperature variation of a longitudinal conductivity, σxx, and b transverse conductivity, σxy, for S3m.The fitting results of the curves in a, and b at 2 K with two-band model are shown in c, and d, respectively.The fitting RMSRE (Root Mean Squared Relative Error) values at different temperatures for e σxx, and f σxy.
FIG. S7.MR analysis of S3m with background-free curvature method.MR plots and their corresponding second derivatives which work as ∆MR at a 2 K, b 5 K, and c 15 K. d-f Fast Fourier transform (FFT) of a-c, respectively.g-i display the Shubnikov-de Haas (SdH) oscillation and Landau level fan plot of d-f, respectively.