Reconfiguring crystal and electronic structures of MoS2 by substitutional doping

Doping of traditional semiconductors has enabled technological applications in modern electronics by tailoring their chemical, optical and electronic properties. However, substitutional doping in two-dimensional semiconductors is at a comparatively early stage, and the resultant effects are less explored. In this work, we report unusual effects of degenerate doping with Nb on structural, electronic and optical characteristics of MoS2 crystals. The doping readily induces a structural transformation from naturally occurring 2H stacking to 3R stacking. Electronically, a strong interaction of the Nb impurity states with the host valence bands drastically and nonlinearly modifies the electronic band structure with the valence band maximum of multilayer MoS2 at the Γ point pushed upward by hybridization with the Nb states. When thinned down to monolayers, in stark contrast, such significant nonlinear effect vanishes, instead resulting in strong and broadband photoluminescence via the formation of exciton complexes tightly bound to neutral acceptors.

S ubstitutional doping of bulk semiconductors, the atomistic substitution with non-isoelectronic impurities, allows to define the type of majority charge carriers and modulate their concentrations to a wide extent, such that they can electrically functionalize as the key component in electronic and optoelectronic devices. For the emerging two-dimensional (2D) or layered semiconductors such as transition metal dichalcogenides, MX 2 (M = transition metal and X = chalcogen), substitutional doping is highly desired to overcome their natively unipolar conduction 1 and substantial contact resistance 2 , despite recent efforts based on surface molecular doping 3,4 and phase transition 5,6 . In this context, substitutional doping in MX 2 has been recently experimentally realized, mainly by replacing the host cation M with other transition metal elements [7][8][9] , often leading to degenerate doping levels of free carriers or even switch in conduction type that was hardly achieved by the other doping techniques 3 . Also, recent theoretical studies report dopinginduced modifications in magnetic 10,11 and catalytic 12 properties of MX 2 .
The atomic d-orbitals in transition metal dopants can exhibit different extents of localization 13 . Namely, they are originally spatially localized with a constant energy level relative to the vacuum level, but these discrete levels may interact with each other to gain dispersion as their wavefunctions overlap, particularly at high doping concentrations. It ultimately results in spatial delocalization of the wavefunctions. Consequently, impurityrelated sub-bands may emerge near the band edge of the host, and at high densities they may eventually hybridize with the host bands, causing bandgap modification and redistribution of the density of states (DOS) 14,15 . As expected from the nearly linear, virtual crystal approximation widely adopted for semiconductor alloys 16 , however, such band restructuring effects often require considerable concentrations of dopants exceeding a few atomic percent to reach the alloying level, with exceptions identified in the so-called highly mismatched alloys 15 . Such dopants-induced band restructuring effects are believed to become stronger at reduced dimensionalities. Heavy doping may also lead to exotic behavior in reduced dimensionalities. For instance, quantumconfined Urbach tail 17 and dynamic surface exciton quenching 18 were observed in heavily doped zero-dimensional (0D) nanocrystals and one-dimensional (1D) carbon nanotubes, respectively. Yet, delocalization of dopant wavefunction and the resultant band restructuring in 2D semiconductors are still experimentally unresolved thus far, despite the fact that quite unique electronic and optical characteristics [19][20][21] were observed in heavily doped 2D materials.
Among the suggested cation dopants substituting molybdenum (Mo) in molybdenum disulfide (MoS 2 ), niobium (Nb) has one less d-electron, and is of particular importance from the thermodynamics point of view: it is predicted to have a negative formation energy 22 . In addition, Nb doping in MoS 2 was suggested to induce considerable degree of charge delocalization as studied in electron paramagnetic resonance experiments 23 and first-principles calculations 22 , as opposed to the case of native point defects (e.g., sulfur vacancies) which only induce localized and nondispersive (flat-energy level) states deep inside the bandgap 24 .
In this work, we report modifications in both crystal structure and electronic structure caused by Nb dopants in monolayer to bulk MoS 2 . We show doping-induced structural conversion from 2H to 3R stacking in the layered structure, accompanied by a renormalization of the valence band structure. In the multilayer or bulk, the wavefunctions of Nb dopants hybridize with the host valence band at the center of the Brillouin zone (Γ V ), thereby dramatically reducing the indirect bandgap. In stark contrast, such electronic restructuring is greatly suppressed in the ultimate 2D limit where Γ V of MoS 2 natively moves down, separating out the Nb states as isolated impurity states, i.e., ionization energy of the impurity states becomes greater due to intrinsically larger bandgap and reduced screening, eventually causing the formation of impurity bound excitons.

Results
Crystal structure of Nb-doped MoS 2 . Both Nb-free and Nbdoped MoS 2 single crystals, Nb x Mo 1-x S 2 (x = 0-0.01), were synthesized using the chemical vapor transport methods (see Methods). For the Nb-doped MoS 2 (MoS 2 :Nb), three representative doping concentrations were prepared to reach the degenerate doping: 0.1, 0.5, and 1% (where the % is defined as the atomic percent of Nb with respect to the density of host Mo sites). In our previous study, X-ray absorption and structural analysis have been employed to verify the substitutionality of the Nb doping as well as uniform distribution of dopants with no phase separation 7 . In addition, Hall-effect technique confirms that free hole concentrations are consistent with the synthesis condition (Supplementary Figure 1). The polymorphs of the synthesized bulk crystals were determined by convergent beam electron diffraction (CBED), and typical images are presented in Fig. 1c. With an obvious difference in symmetries between the MoS 2 and MoS 2 :Nb, the CBED patterns confirm the 2H and 3R stacking for the Nb-free and Nb-doped MoS 2 , respectively. The 2H and 3R are two different types of stacking order of MX 2 in the trigonal prismatic coordination, described by the space groups of P6 3 / mmc and R3m (Fig. 1a), respectively. Unlike the 2H phase, owing to the non-centrosymmetric crystal symmetry and the resultant pseudo-spin polarization in the valence bands, the 3R phase recently attracts great interests for study of valleytronics 25 . While the 2H phase is found more frequently for semiconducting MX 2 , the 3R phase is thermodynamically also stable thanks to the nearly negligible difference between their ground-state total energies, and is observed in most of Nb-doped samples we tested (38 out of 41 crystals; see additional CBED images in Supplementary Figure 2 and Note 1). A triangular shape of screw dislocation spiral is also often observed from the surface of the MoS 2 :Nb crystals, confirming their non-centrosymmetric stacking 26 (Supplementary Figure 3). Density functional theory (DFT) calculations were performed to confirm the stacking reconfiguration of MoS 2 from 2H phase in the updoped to 3R phase in the Nb-doped case. In Fig. 1b, the total energy difference between bilayer 2H and 3R MoS 2 :Nb (ΔE 0 ) turns into negative values, −0.13 to −0.41 meV per atom, depending on the Nb doping concentrations. We note that the 3R phase is energetically more stable than the 2H phase once the Nb dopant is added to MoS 2 supercell (the 6 × 6 (4 × 4) supercell is employed for Ⓐ and Ⓑ (Ⓒ) doping configurations). And, it becomes increasingly more stable upon further Nb incorporation. In short, it is the thermodynamic driving force that reorders the stacking stability, i.e., the 3R phase will naturally form once the Nb dopants are included in the crystal growth at elevated temperatures of~1000°C for > 2 weeks. The relative energetic stability can be related to free-carrier screening by holes residing in the d z 2 bands with some delocalization 23 in Nb-doped MoS 2 , effectively lowering the total energy of non-centrosymmetric 3R-MoS 2 . We note that the 3R phase was also experimentally observed in other heavily doped MoS 2 27 . Furthermore, the dynamic stability of such 3R MoS 2 :Nb bilayer is confirmed by the phonon dispersion showing no imaginary part, a criteria employed for other poly types of MoS 2 28,29 (Supplementary Figure 4). Consequently, we conclude that Nb doping provides a reliable route of acquiring the uncommon, non-centrosymmetric structure of MX 2 bulk crystals for potential valleytronic applications, circumventing the need for large-area monolayers, nor delicate adjustment of the growth temperature gradient 25 (the routinely adopted means to stabilize 3R MoS 2 ).
The rearrangement of stacking order by the substitutional Nb doping also renormalizes the A and B exciton transitions at the K point of the Brillouin zone, as evidenced in optical reflectance spectra for MoS 2 and MoS 2 :Nb bulk crystals measured at 4.5 K (Fig. 1d). First, the A exciton feature redshifts slightly by~10 meV in the MoS 2 :Nb crystal, determined by the lineshape fits to a Lorentz oscillator model 30 (Supplementary Figure 5 and Note 2). Second and more notably, the A-B energy splitting decreases from 0.19 to 0.15 eV, i.e., a much greater redshift of the B exciton peak. Also, the flatter dispersion along the K-H path in the 3R system 25 eliminates a distinctive H-point excitonic transition that is normally observable in bulk 2H-MoS 2 at low temperatures 30 . Quantitatively consistent with earlier reports on undoped 3R-MoS 2 25 , it can be seen that these renormalization effects in the A and B exciton energies come from the 2H-3R structural transformation in the MoS 2 :Nb, rather than a direct electronic effect associated with the Nb doping. Aside from these structural effects, strong and broad absorption is observed at 20-100 meV below the A exciton peak in the MoS 2 :Nb crystal (Fig. 1d), and is attributed to Nb impurity states located above the valence band edge at the K point (K V ). Yet, as for what we confirm at the K point, the direct electronic influence of these Nb impurity states on the host MoS 2 bands is negligible, akin to the case of conventional bulk semiconductors in the heavy doping regime.
Electronic restructuring of valence bands in the bilayer case. Few-layer structures were obtained by mechanical isolation from the bulk crystals, which enables further visualization of effects of the different stacking sequences. Figure 2a displays highresolution transmission electron microscopy (HRTEM) images of bilayers (2L) of undoped MoS 2 and MoS 2 :Nb. Here the 2H stacking of undoped MoS 2 (top panel) is easily recognized with no atom located at the center of each hexagon. In the MoS 2 :Nb bilayer, on the contrary, additional atomic columns appear in the center of each hexagon, with a slight but noticeable intensity difference between nearest neighbor columns. Atomically precise positioning and stack ordering of layers are possible in such 3Rstacking bilayer with the aid of line intensity profile 31 (Supplementary Figure 6). With extensive HRTEM investigation, we confirm that structurally the Nb doping only induces the 3R polytypism in MoS 2 , preserving its single crystallinity without causing extended defects and local lattice distortion.
Bulk MoS 2 is an indirect-bandgap semiconductor with negligible photoluminescence (PL), thus prevents access to its electronic bands in light emission spectroscopy. When the MoS 2 crystal is thinned to a few layers, however, noticeable PL emerges with two major features. A main peak at~1.85 eV is associated with (hot-carrier) transitions across the direct bandgap at the Kvalley, whereas a relatively weak PL peak at lower energies ( Fig. 2b) arises from transitions across the indirect bandgap involving the Γ V , which constitutes a combination of Mo À d z 2 and S À p z orbitals. Since these orbitals extend along the zdirection with strong overlap between neighboring layers, the indirect PL peak acts as an indicator of interlayer coupling: the lower the indirect PL peak energy is, the stronger the interlayer coupling 32 . 3R stacking is known to have a slightly shorter interlayer spacing compared to that of 2H, so the indirect PL of 3R bilayers is expected to redshift. However, to separate the structural (2H vs. 3R) effect from electronic (doped vs. undoped) effect in Nb-doped 3R bilayers, it is necessary to also include Nbfree 3R bilayers in the comparison. Indeed, the indirect PL of Nbfree 3R bilayers redshifts slightly by~30 meV, as shown in Fig. 2b. In contrast, the observed redshift of indirect PL from the Nbdoped 3R bilayers is much greater, up to~140 meV, and it indeed monotonically redshifts further with Nb fraction (x). These effects suggest additional mechanism of Nb doping beyond the mere 2H-3R structural conversion.
The unusual evolution of indirect optical transition in 2L-MoS 2 :Nb is attributed to a valence band modification by the Nb impurity states. As seen in Fig. 3, the impurity level of Nb replacing Mo (denoted as E I ) is theoretically known to be located below the Γ V in MoS 2 for both bulk and bilayer cases 33 , also judged from the reflectance spectra in Fig. 1d  in resonance with the valence bands. In this case, an upshift of Γ V (hence a reduction of the indirect bandgap) is expected from a simple two-level repulsion model, as exemplified in similar cases of valence (conduction) bands of GaAs where As is partially substituted with Bi or Sb (N) 14,34 . More generally, this two-level band-anticrossing model has successfully described the electronic band restructuring in a wide range of substitutional semiconductor alloys 15,35,36 . It predicts the restructuring and hybridization of host bands (E V (k)) with impurity states (E I ), giving rise to the two newly formed sub-bands expressed as where C is a parameter describing the hybridization strength, and x is the fraction of the impurities. For the case of bilayer MoS 2 : Nb, E V (k) and E I are the original valence band dispersion of MoS 2 and the originally non-dispersive Nb energy level, respectively, and E ± (k) is the two newly formed sub-bands as a result of the anticrossing interaction between E V (k) and E I , as shown in Fig. 3a. Even a sub-1% doping causes a significant band restructuring in MoS 2 :Nb, suggesting a strong band-anticrossing interaction with a high interaction parameter, C.
We also performed ab initio calculations on the electronic structures using DFT within the local density approximation and the results confirmed this model. Figure 2c shows the partial density of states (PDOS) plot for 2L undoped and Nb-doped MoS 2 in the 3R-type stacking order. Upon Nb doping, the indirect bandgap becomes narrower and the Fermi level (E F ) downshifts consistent with our experimental observation. As seen from the undoped MoS 2 , the band edges of both conduction and valence bands are mainly composed of Mo 4d states. For 2L-MoS 2 :Nb, the Nb 4d state contributes most significantly to the valence band maximum, Γ V , while its influence on the conduction band minimum is nearly negligible. Notably, contribution from the host Mo 4d state to the valence band edge is found to be sensitive to its distance to the Nb dopant; the closer to the Nb atom, the greater contribution of the Mo to the valence band maximum. This is not the case for the conduction band, hence it supports our band-anticrossing-type hybridization between the Nb and Mo 4d states at Γ V , as outlined above. Atomic contribution-resolved band structures further elucidate the case by visualizing the hybridization at Γ V (Supplementary Figure 7). That is, as the doping concentration increases, the valence band maximum (Γ V ) gains more E I character from the Nb dopants simply because of the greater fraction of Nb, x; in the meantime, our calculation suggests that the interaction parameter, C, can be simultaneously enhanced due to the closer physical distance to the dopants as x increases, all ultimately driving the indirect PL peak to monotonically redshift with x.
Bound excitons in the monolayer limit. When further thinning down to the monolayer limit, 2H-and 3R-MoS 2 no longer have a structural distinction (Supplementary Figure 8). From the bulk to monolayer, MoS 2 layers transit from indirect to direct bandgap, as well as experience reduced dielectric screening. These effects lead to well-known strong PL and large excitonic binding energy in monolayer MoS 2 . For monolayer MoS 2 , a drastic selfdownshift of the Γ V point makes the K V pockets the valence band maximum, across which direct optical transitions take place. In such a system, the Nb impurity states are now located above the valence band maximums, no longer crossing the valence bands of MoS 2 (Fig. 3b) and hence becoming isolated in-gap impurity states. Such reconfiguration of relative location between the host valence bands and the impurity level results in a few notable changes.
First, we observed an enhanced, broad-band and redshifted PL emission at room temperature (RT) from the monolayer MoS 2 : Nb, as shown in Fig. 4a. Combined with the corresponding absorption spectra (Supplementary Figure 9 and Note 3), it can be attributed to optically excited excitons (X) binding to neutral acceptor (A 0 ) resulting in the formation of A 0 X complexes. While an ionized acceptor (A − ) does not usually bind an exciton, the binding energy of exciton is in general the highest for a neutral acceptor according to the Haynes rule 37 . As illustrated in Fig. 3, the former (A − X) and latter (A 0 X) cases correspond to the bilayer and monolayer MoS 2 :Nb, respectively, under the assumption that E I becomes a deep impurity level and a portion of acceptors are no longer thermally ionized for the monolayer MoS 2 :Nb according to Fermi-Dirac statistics. In this sense, a direct PL of 2L-MoS 2 :Nb is only weakly affected by the presence of E I (Fig. 2b) but that of monolayer MoS 2 :Nb substantially changes. It is also noted that such impurity-bound PL becomes even broader and more redshifted upon increase in the Nb doping concentration. Next, PL of the monolayer MoS 2 :Nb is seen to be more than an order of magnitude stronger than that of the undoped monolayer  Fig. 4b. That is, while the decay dynamics of the exciton states at~1.88 eV in undoped MoS 2 is rather fast (lifetime, τ, on the order of 100 ps 38 ) approaching the instrumental response function (IRF), the PL at~1.68 eV in the monolayer MoS 2 :Nb shows significantly slower bi-exponential decay dynamics (τ 1 = 0.7 ± 0.05 ns and τ 2 = 4.3 ± 0.3 ns). This is in great contrast with the earlier PL studies on native defects like chalcogen atomic vacancies where PL decay times are measured as a few hundred ps 39,40 . Moreover, although such defect-related PL feature is observed only at cryogenic temperatures 24 , here the Nb-induced sub-band PL is much brighter and stable even at RT. Doping could indeed, in some cases, enhance radiative efficiency and could, for example, even enable RT lasing in GaAs nanowires 41 , yet such remarkable luminescence enhancement has not been reported in 2D semiconductors. Last, PL from undoped and Nb-doped monolayer MoS 2 shows a different laser excitation power dependence (Fig. 4c). As laser power steadily increases, PL in the monolayer MoS 2 :Nb, stemming from the impurity-bound excitons as described above, monotonically shifts toward higher energy, as similarly observed for localized excitons 42 whilst the undoped MoS 2 monolayer only displays a slight decrease presumably due to laser heating effect. Photoluminescence excitation (PLE) spectroscopy was also performed over a wide range of excitation energies to further understand this unusual PL behavior. It can be seen from Fig. 5 that the PLE signal appears when the excitation energy is above the A exciton absorption, and reaches a maximum when the excitation energy is equal to the B exciton absorption, as similarly reported in undoped monolayer MoS 2 43 . However, this enhancement of emission is not observed for excitation at C resonance indicating that an exciton first needs to form to bind to the impurities. Again, this PLE behavior is different from what is known for defect-related emission 44 . The fact that the PLE emission starts only after the excitation energy is higher than the A exciton energy (rather than when the excitation energy is equal to the Nb impurity energy level), is also an indication that the Nb impurity states no longer form a dispersive band in this case (unlike at the Γ V band in multilayer or bulk). Therefore, the Nb states stay as non-dispersive energy levels not hybridizing with the K V host bands (Fig. 3b), instead contributing to the formation of bound excitons. Indeed, hybridization of impurity states with the host bands is expected to be much weaker at the Brillouin zone edge (K) than at the center (Γ). This is because the interaction strength (C) in the band anticrossing model becomes the maximum at k = 0 (i.e., at Γ V ), with a rapidly decreasing k-dependence that follows the Fourier transform of the spatially-decaying overlap integral between the impurity wavefunction and the Wannier function of the host bands 45 . Moreover, the hybridization at K V is further weakened due to the farther energy separation of E I from E V (k) at K V in the monolayer (than Γ V in multilayer/bulk).

Discussion
We have demonstrated restructuring of layer stacking and electronic bands of MoS 2 by introducing substitutional Nb dopants. It drives the stacking sequence from ABABAB (2H) to ABCABC (3R), resulting in broken inversion symmetry even for multilayer MoS 2 . Our study reveals that the Nb impurity states in the bulk are quite dispersive and strongly interact with the extended states of host valence band, which is also supported by theoretical analysis. Strong interaction mainly occurs at the Γ V of the Brillouin zone for multilayer or bulk MoS 2 , suggesting a valence band hybridization. It is highlighted that all the structural and electronic modifications in multilayer MoS 2 :Nb happen still within the relatively low substitutional fractions of doping (rather than

Methods
Materials preparation. Undoped and Nb-doped MoS 2 single crystals were grown by a chemical vapor transport (CVT) method with the use of I 2 as a transport agent. For comparison, natural 2H-MoS 2 crystals were also purchased from SPI Supplies. High purity elements including Mo, S, Nb (all 99.99% purity), and the iodine transport agent were used for the crystal growth. The molar ratio of Mo/S was kept to 1:2 and the substitutional doping concentrations of niobium were designed to become 0.1, 0.5% and 1%, respectively. Afterwards, all reactants were placed in quartz tubes that were evacuated below~2 × 10 −5 Torr and sealed by an oxyhydrogen flame. A horizontal three-zone furnace was utilized to maintain an optimal temperature gradient for the diffusion of iodine transport agent. The quartz tubes were placed into the furnace where the high temperature zone was set at 1050°C and the low temperature zone was set at 935°C for 720 h.
Structural characterization. Both convergent beam electron diffraction (CBED) patterns and high-resolution transmission electron microscopy (HRTEM) images were acquired on a FEI Titan environmental TEM 80-300 at 80 kV accelerating voltage. A negative spherical aberration imaging technique with monochromated electron beam was used for HRTEM images of bilayer samples.
Optical measurements. The reflectance measurements on the bulk crystals used light from a 75 W Xenon lamp dispersed by a 1/2 m grating monochromator and a photo-multiplier tube detector, with the sample cooled using a continuous flow liquid He cryostat. Time-resolved photoluminescence measurements were performed using the time-correlated single-photon counting technique at room temperature. MoS 2 samples were excited with an ultrafast optical pulse of photon energy of 2.06 eV and power of 4.5 μJ/cm 2 . Pulse duration and repetition rate were 50 ps and 20 MHz, respectively. The PL spectra were obtained with Raman spectrometers in the back scattering geometry and continuous wave 488 nm laser as the excitation source. For photoluminescence excitation (PLE) measurements, a Fianium supercontinuum white laser (SC450) coupled to a laser line tuneable filter (LLTF) was used to provide excitation sources. The excitation lasers, ranging from 460 nm to 760 nm, have a bandwidth of 1−2 nm and their power was kept below 0.5 μW for all the measurements in order to keep the excitation rate in the linear regime and to avoid possible damage to the samples. Multiple measurements were made for each sample to check the reproducibility of the results.
Electronic structure calculations. Density functional theory (DFT) calculations were performed on the bilayer MoS 2 structures to elucidate their bandstructures before and after Nb doping. Geometry optimizations were performed using DFT implemented in the Vienna Ab initio Simulation Package (VASP) 46,47 within a Projected Augmented Wave (PAW) 48 basis and with the Perdew, Burke and Ernzerhof (PBE) functional 49 , ensuring sufficient vacuum (at least 20 Å) between periodic images along the z-direction (perpendicular to the MoS 2 plane) to minimize spurious interactions. During structural optimization, all atomic coordinates and lattice vectors were fully relaxed until the absolute value of the forces acting on each atom was less than 0.01 eV/Å. We checked that sufficient vacuum remains after relaxation. Plane-wave cutoffs were set to 400 eV and van der Waals interactions were accounted for via the DFT-D2 scheme 50 . For both undoped and Nbdoped MoS 2 , a 6 × 6 supercell was used along with a 4 × 4 × 1 k-point grid, unless otherwise stated, with the Monkhorst-Pack sampling 51 .
Data availability. The data that support the findings of this study are available from the corresponding authors on request.  5 Photoluminescence excitation (PLE) data of the monolayer MoS 2 :Nb (1%). Left and right y axis correspond to PL emission energy at 1.43 eV and differential reflectance spectra (ΔR/R), respectively, as a function of excitation energy. Top inset displays PLE intensity 2D map where the color scale represents emission intensity, and was collected with an excitation power of~0.5 μW at 300 K