Evidence of the Coulomb gap in the density of states of MoS$_2$

$\mathrm{MoS_2}$ is an emergent van der Waals material that shows promising prospects in semiconductor industry and optoelectronic applications. However, its electronic properties are not yet fully understood. In particular, the nature of the insulating state at low carrier density deserves further investigation, as it is important for fundamental research and applications. In this study, we investigate the insulating state of a dual-gated exfoliated bilayer $\mathrm{MoS_2}$ field-effect transistor by performing magnetotransport experiments. We observe positive and non-saturating magnetoresistance, in a regime where only one band contributes to electron transport. At low electron density ($\sim 1.4\times 10^{12}~\mathrm{cm^{-2}}$) and a perpendicular magnetic field of 7 Tesla, the resistance exceeds by more than one order of magnitude the zero field resistance and exponentially drops with increasing temperature. We attribute this observation to strong electron localization. Both temperature and magnetic field dependence can, at least qualitatively, be described by the Efros-Shklovskii law, predicting the formation of a Coulomb gap in the density of states due to Coulomb interactions. However, the localization length obtained from fitting the temperature dependence exceeds by more than one order of magnitude the one obtained from the magnetic field dependence. We attribute this discrepancy to the presence of a nearby metallic gate, which provides electrostatic screening and thus reduces long-range Coulomb interactions. The result of our study suggests that the insulating state of $\mathrm{MoS_2}$ originates from a combination of disorder-driven electron localization and Coulomb interactions.


I. INTRODUCTION
The resistivity ρ of some semiconductors shows a metal-insulator transition as a function of the electron density n [1]. For densities larger than a critical value n c the resistivity shows a metallic temperature dependence (dρ/dT > 0), while below n c it shows an insulating temperature dependence (dρ/dT < 0). This metal-insulator transition attracted great interest in the late 1990s [2][3][4][5]. In two-dimensional (2D) semiconductors the origin of the metallic phase is controversial [6][7][8], as it was predicted that any amount of defects would inexorably lead to electron localization at zero temperature in a 2D system [9,10]. The insulating phase at low densities can be due to either intriguing correlated states, like Wigner crystals [11] or disorder-induced electron localization [5], as well as a combination of the two effects.
In highly disordered systems, charge transport at low temperatures occurs via electron hopping between localized states [12], known as variable-range hopping (VRH). The conductivity in hopping transport at zero magnetic field is usually described by an exponential dependence on the temperature of the form where T 0 and p ≤ 1 are constants that depend on the hopping mechanism. In a non-interacting system, the density of states close to the Fermi energy is constant (but finite) and the conductivity is described by Mott's law [13], for which p = 1/3 (for two-dimensional systems). When electrons are strongly localized, the long-range Coulomb potential is not efficiently screened. Electron correlations result in a Coulomb gap in the density of states close to the Fermi energy [14,15]. The modified density of states changes the temperature dependence of the hopping conductivity, which is now characterized by the parameter p = 1/2, as described by the Efros-Shklovskii (ES) theory [16]. The insulating phase of MoS 2 has been experimentally studied in monolayers [17] and multilayers [18,19], where both thermally activated transport at intermediate temperatures (T ∼ 50 K to 100 K) and Mott VRH transport at lower temperatures have been observed. In addition, it is expected that electron-electron interactions play a major role in determining the electronic properties due to the large electron effective mass [m * ≈ (0.4 − 0.6)m 0 ] of MoS 2 , especially at low densities. Indeed, signatures of interaction effects have already been reported in the literature [20,21], among which there was also the observation of a Wigner crystal in MoSe 2 [22]. Therefore, MoS 2 , and in general, semiconducting transition metal dichalcogenides (TMDs), are good candidates for the observation of the Coulomb gap in the density of states. However, the observation of interaction effects is restricted to low densities, where the Coulomb energy dominates over the kinetic energy of electrons. Transport experiments in this density range are challenging in most materials and require low defect densities [23]. The observation of the Coulomb gap in MoS 2 remains to date elusive [19] due to the large density of intrinsic defects.
Here, we investigate magnetotransport in bilayer MoS 2 encapsulated in hexagonal boron nitride (hBN). We first demonstrate the high quality of our device and a complete understanding of the bands that contribute to electron transport. Then we tune the density, such that electron transport occurs in a single (twofold degenerate) band. In this regime, we observe a metal-insulator transition at the electron density n c ≈ 1.7×10 12 cm −2 . Below this transition, the zero-field resistance and the magnetoresistance show an exponential decay with increasing temperature, which is qualitatively consistent with the Efros-Shklovskii law [16], suggesting that electron correlations open a gap in the density of states. However, the localization length obtained by fitting the temperature dependence exceeds by more than one order of magnitude the one obtained from the magnetic field dependence. We attribute this discrepancy to the presence of a nearby metallic gate, which provides electrostatic screening and thus reduces long-range Coulomb interactions.

II. RESULTS AND DISCUSSION
In this study, a bilayer MoS 2 with dual-gated architecture is employed to study electron transport. Figure 1(a) shows optical images (upper panel) and a schematic side view (lower panel) of the device. The fabrication starts by assembling a thin hBN and graphite layers onto a silicon/silicon dioxide chip (285 nm of oxide layer). The graphite serves as the bottom gate, while the hBN is the gate dielectric material. The layers are stacked together using a dry-transfer technique (see Refs. [21,24,25] for details). We pattern metallic contacts (Cr/Au: 5 nm/10 nm) with standard electron beam lithography and electron beam evaporation. The region of the contacts is then cleaned thoroughly with the tip of an atomic force microscope in contact mode. This step is crucial to remove the residues of the lithography process. The bilayer MoS 2 is obtained by mechanically cleaving a bulk MoS 2 crystal from natural sources (SPI supply) and identified by optical contrast, which has proven to be a very reliable method [26][27][28][29]. A second stack with a layer of hBN and the bilayer MoS 2 is assembled, then aligned and deposited on top of the contacts. Both exfoliation and assembling take place in a glove box with argon atmosphere (H 2 O, O 2 < 0.1 ppm). The metallic top gate is patterned with standard lithography processes and covers the entire MoS 2 flake. In the final step, we vacuum anneal the sample at 250 • C for 4 h to improve the contact interface. All measurements presented in this work (if not explicitly stated otherwise) are performed at a temperature of 1.3 K at low frequency (∼ 30 Hz) and low excitation voltage (V rms SD = 100 µV) with standard lock-in techniques.
First, we characterize the contacts to ensure Ohmic behavior and low contact resistance. The dc output characteristics are shown in Fig. 1(b) for different top gate voltages (V TG ). The two-terminal resistance changes from < 1 kΩ to ∼ 1 MΩ as a function of V TG . The current shows a linear dependence on the applied source-drain voltage V SD down to the lowest V TG (inset), indicating a vanishing Schottky barrier. The average contact resistance of source and drain contacts (1 and 2) is < 300 Ω for V TG ≥ 6 V and does not depend on the applied bottom gate voltage (V BG ). To the best of our knowledge, this is one of the lowest values reported in the literature and is comparable to bismuth contacts [30]. The vanishing Schottky barrier and the low contact resistance allow us to study electron transport at low densities (down to 1.4 × 10 12 cm −2 ) at low source-drain voltage (100 µV).
The four-terminal conductivity (σ) as a function of V TG is shown in Fig. 1(c). We identify two voltage ranges, where the conductivity features a linear dependence on V TG , which differs by the slope. The two slopes of σ(V TG ) yield the mobilities ∼ 1000 cm 2 V −1 s −1 (dotted line) and ∼ 2400 cm 2 V −1 s −1 (dashed line). This specific shape of the four-terminal conductivity is a general property of high-quality single-gated MoS 2 devices [see Fig. 1(c) in Ref. [20] and Fig. 1(d) in the supplemental material of Ref. [21] for a direct comparison]. The kink marked with the number 2 is related to the population of the upper spin-orbit (SO) split band, as will be demonstrated below. Figure 1(d) shows the derivative of the conductivity dσ/dV BG as a function of V TG and V BG . For V BG < 0.7 V the V BG dependence of the conductivity is monotonic and resembles its V TG dependence. There are two distinct rates dσ/dV BG that are separated by a black dashed line (number 2). This line follows constant density conditions, separating the single band from the two-band regime. From this result, we conclude that the separation between the bands is not displacement field dependent (i.e., the SO gap is not tunable with the electric field). In contrast, for V BG > 0.7 V the dependence of the conductivity on V BG is nonmonotonic. The conductivity features a local maximum (number 3) that is almost independent of V TG . The origin of the negative dσ/dV BG is related to the population of the bottom layer and will be discussed in Sec. II A.
So far we have characterized the field effect transistor at zero magnetic field. Understanding the contribution of the bands in multilayer MoS 2 devices is an essential step for interpreting their electronic properties. In the following, we will demonstrate that all features observed in the conductivity can be attributed to the population of different bands in the MoS 2 bilayer. To verify this hypothesis, we now turn our attention to magnetotransport measurements, from which we determine the electron density and the band occupation.

A. Band and layer occupation
The aim of this part is to determine how the different bands are filled with electrons as we change the gate voltages. For this purpose we measure the magnetoconductivity σ xx , applying a perpendicular magnetic field B. The conductivity σ xx is obtained by tensor inversion of the two-dimensional resistivity with components ρ xx = W V 3,4 /(IL 2 ) and ρ xy = V 3,5 /I. From the Shubnikov-de Haas oscillations (SdHO) we determine the density n i of band i as we did in previous works [21,24,25]. The total density is given by where C T = 146 nF/cm 2 and C B = 225 nF/cm 2 . The total density is related to the band densities via n tot = i n i . We do not determine the total density from the Hall effect, because the contacts extend across most of the conducting channel, leading to an overestimation of the density. Figure 2(a) shows the derivative of the magnetoconductivity dσ xx /dV BG as a function of B and V BG for V TG = 11.5 V. The conductivity features oscillations periodic in B −1 , as we expect for SdHO. From the Fourier spectrum of the dσ xx /dV BG (1/B) (see Ref. [24,25] for details) we determine the densities n i shown in Fig. 2(b). As in Ref. [24], we attribute a twofold valley degeneracy to each band i, which accounts for the two K valleys (see schematics in the lower right panel of Fig. 2). In the regime V BG < 0.7 V, where we find only two frequencies in the Fourier spectrum, there is remarkable agreement between the calculated n tot and the experimentally defined density n 1 + n 2 . In our interpretation, the densities n 1 and n 2 belong to the top MoS 2 layer. In the regime V BG > 0.7 V, the density n 1 + n 2 saturates, because the bottom MoS 2 layer becomes conducting, screening the field effect of the bottom gate. This behavior is in complete agreement with the behavior found and explained in previous studies [24,25].
The layer occupation can be controlled in dual-gated devices. For V TG > 0 and V BG < 0 only the top layer is occupied by electrons and we can tune between singleand double-band transport. The onset of n 2 (marked by the number 2) corresponds to the population of the upper SO split band of the top layer (see schematics of the band structure in Fig 2). To support this interpretation we determine the density n SO ≈ 3.5 × 10 12 cm −2 that is required to start filling the upper SO split band. In an effective mass approximation, this density corresponds to an energy ∆ SO ∼ 14 meV, in agreement with the value reported in Ref. [24]. Figure 2(c) shows dσ xx /dV BG at B = 7 T as a function of the voltage applied to the top and bottom gates. This measurement defines a phase diagram for the bilayer MoS 2 that is divided into five different regimes by the dashed (dotted) lines. Each line indicates the onset of a specific band (from 1 to 4 increasing V BG ). The bilayer is tuned from an insulating into a conducting phase, where up to four bands contribute to electron transport. In the remaining part of this paper, we focus on the narrow single-band transport regime enclosed by the dashed lines in Fig. 2(c).

B. Strong electron localization
We now turn our attention to the magnetoresistance in the single-band transport regime. According to the Drude model, the longitudinal resistance does not de-pend on the magnetic field. However, we observe a nonsaturating positive magnetoresistance in the low-density range (< 3.5 × 10 12 cm −2 ). The magnetoresistance increases by more than one order of magnitude at B = 7 T and n tot = 1.4 × 10 12 cm −2 , as shown in Figs. 3(a)-3(c). A qualitatively similar positive magnetoresistance was observed in monolayer MoS 2 [31], where a negative magnetoresistance (due to weak localization) gradually turned into a positive magnetoresistance as the temperature and density were lowered. This effect was attributed to the transition from weak localization to weak antilocalization, despite the authors realizing that the shape of the magnetoresistance was not well described by the theory of weak antilocalization. Also, in another bilayer [24,32] and three-layer MoS 2 devices [25] a similar behavior has been observed but was not further investigated. The observation of a positive magnetoresistance in samples prepared and studied in different research groups suggests that there is a common origin for this effect. Owing to the relatively low electron mobility in MoS 2 and the large intrinsic defect density [33], we consider the role of disorder [34] and electron localization to describe the nonsaturating positive magnetoresistance.

Figures 3(a)-3(c)
show the four-terminal longitudinal resistivity ρ as a function T for three different densities (from 1.4 × 10 12 cm −2 to 3.5 × 10 12 cm −2 ) and various magnetic fields. At zero magnetic field (violet curve) we observe a metal-insulator transition as we lower the density. It is remarkable that this transition occurs at relatively high densities ∼ 1.7 × 10 12 cm −2 , where the ratio between Coulomb and kinetic energy is still moderate r s = E C /E kin ∼ 7. Therefore, the transition is more likely to be a manifestation of strong electron localization when the Fermi energy approaches the bottom of the conduction band, rather than an insulating state due to correlations. At the lowest temperature, we estimate the product (k F × ℓ e ) between Fermi wave vector (k F ) and electron mean free path (ℓ e ). Strong electron localization occurs for k F × ℓ e ≤ 1 [as shown in Fig. 3(a)], while band conduction occurs for k F × ℓ e ≫ 1. Applying a perpendicular magnetic field breaks this condition, inducing a metal-insulator transition even for k F × ℓ e > 1 [see Fig. 3 The resistivity is plotted on a logarithmic scale as a function of T −1/2 in Fig. 3(d) for the lowest density and various magnetic fields. At low temperatures (< 10 K) all the curves follow a linear dependence, as predicted by the ES theory. The slope of the curves increases with increasing magnetic field. Therefore, we define a parameter T * (B) that depends on the magnetic field, such that the resistivity is described by This function seems to capture the temperature dependence of the resistivity for T < 10 K, while above this temperature the resistance shows a weaker temperature and magnetic field dependence. The parameter T 1/2 * obtained from the fit is shown in Fig. 3(f). This parameter depends roughly linearly (at least above B ≈ 1 T) on the applied magnetic field and the slope decreases with increasing electron density. Similarly, we plot the resistivity for different densities [see Fig. 3(e)] and extrapolate the parameter T 1/2 * as a function of densities. The density dependence of T 1/2 * is shown in Fig. 3(g), where we see T 1/2 * increasing rapidly for decreasing density. We considered other models to describe our data, such as Mott's law (p = 1/3) and thermally activated nearestneighbor hopping (p = 1). While thermally activated hopping clearly fails in describing the temperature de- pendence of the resistivity, we obtain reasonable fits also with Mott's law, which yields a slightly larger least-meansquare error. Based on the temperature dependence of the resistivity, we cannot reliably distinguish between the two models because the resistivity does not change by several orders of magnitude. For this reason, we also consider the role played by the magnetic field.
At finite magnetic fields, the tails of the electron wave functions decay faster and the overlap between the localized wave functions decreases. This leads to a reduction of the tunneling probability and thus to an increase in the resistance. This effect results in an exponential increase of the resistance of the form [12] where m depends on the range of magnetic field and the assumptions made in the theoretical model [12,35,36].
First, we focus on the low magnetic field range (ℓ B ≫ a, ℓ B being the magnetic length and a the localization length), where the effect of B can be treated as a small correction to T 0 . The low magnetic field correction of the percolation parameter ξ is given by [37] ∆ξ(B) = ξ(B) − ξ(0) = C 2 a 4 ℓ 4 where C 2 = 0.002 is a numerical parameter, and ℓ B = h/eB. This expression yields the magnetoresistivity which can be expanded in a quadratic expression for A(T )B 2 << 1. Figures 4(a) and 4(b) show the low magnetic field range of the resistivity at the density 1.54 × 10 12 cm −2 (i.e., below n c ) and for two exemplary temperatures. At T = 10 K the resistivity shows a quadratic dependence on B, as predicted by Eq. (2). At this temperature, the model fits the data well in the entire magnetic field range probed in our experiments, as shown by the inset of Fig. 4(a). The range of magnetic field for which this model is able to describe the data shrinks with lowering the temperature because the parameter ∆ξ(B) grows with lowering the temperature (∝ T −3/2 ), thus limiting the magnetic field range for which our approximation is valid. Therefore, it is not surprising that at T = 1.3 K [ Fig. 4(b)] the model deviates from the data at B c ≈ 1 T where ∆ξ/ξ 0 ≈ 0.5 becomes a large correction.
Fitting the magnetoresistance at different temperatures provides the temperature dependence of the curvature, which is the only fitting parameter. The result of the fit is shown in Fig. 4(c) on a log-log scale for the density 1.54 × 10 12 cm −2 . The data follow the temperature dependence T −(1.6±0.1) . For comparison, we show the temperature dependence (T −1 ) predicted by the Mott theory, which clearly deviates from the trend of our data, while the data are in good agreement with ES theory (T −3/2 ). Based on this observation, we conclude that, in "diluted" MoS 2 , long-range Coulomb interactions lead to the formation of a gap in the density of states at the Fermi energy. This conclusion is further supported by the interaction parameter r s ∼ 7, which is modest but indicates that the Coulomb energy is significantly larger than the kinetic energy of free electrons.
At this point, we would like to discuss our quantitative results. By fitting the data with ES law we can estimate the localization length a according to where ϵ r ≈ 7 is the relative dielectric constant of MoS 2 , and C 1 = 6.2 is a numerical constant [37]. This equation, however, provides localization lengths of a few micrometers, which is surprisingly large considering the interparticle spacing (n −1/2 ∼ 10 nm). To verify the validity of this result, we estimate the localization length from the curvature of the magnetoresistance and compare the two results. By inserting T 0 (obtained from the temperature dependence) in the definition of A(T ) we obtain a ≈ 100 nm at a density of 1.54 × 10 12 cm −2 . This result differs by more than one order of magnitude from the one obtained by Eq. (3).
Here, we offer an argument that might explain the discrepancy between these two results. We note that in our device the metallic gate is only separated from the 2D electron system by a thin hBN layer (d = 13 nm). This distance is comparable to the interparticle distance n −1/2 ∼ 10 nm. As a consequence, the metallic gate screens Coulomb interactions at large distances (r ≫ d), where the Coulomb potential becomes dipolelike (∝ r −3 ). The presence of the metal plate partially suppresses the Coulomb gap due to its screening effect [38][39][40][41][42][43]. For small hopping distances the density of available states still has a gap, while for long hopping distances the density of states is constant. At lower temperatures, long-distance hopping becomes favorable and VRH is determined by the constant density of states (like in the Mott VRH). Therefore, it is expected that in the presence of a nearby metallic gate, there is a critical temperature below which the system shows a transition from ES to the "screened" Mott VRH [41]. This VRH is different from the usual Mott VRH, which is characterized by a different density of states.
We estimate this temperature to be around 20 K by assuming a ≈ 100 nm, placing our system just below the transition temperature. In contrast, the hopping distance is of the order of d, thus still in the Coulomb gap regime. Therefore, it is not clear from these characteristic quantities in which regime our sample is. As we argued above, we cannot, based on the temperature dependence, clearly distinguish between ES and Mott VRH. Under this condition, Eq. (3) may not be valid, as we are close to the transition between the two hopping transport regimes.
On the other hand, the magnetic field promotes electron localization and reduces the probability of long-range hopping (less overlap of wave function tails). As shown in Fig. 3(f), the magnetic field dependence of the fitting parameter T * (B) ∝ B 2 is compatible with the wavefunction shrinking effect. In addition, the curvature A(T ) clearly establishes that at finite field VRH follows ES theory. For this reason, we consider our estimation of the localization length from the analysis shown in Fig. 4 to be valid, from which we estimate a ∼ 100 nm. Finally, we consider the asymptotic limit at high magnetic fields (a ≫ ℓ B ). The typical hopping distance is strongly reduced by the effect of the magnetic field and the typical hopping energy increases. As a consequence, the Coulomb gap might be negligible at high B-fields, as proposed by Nguyen [37]. The resistivity is expected to follow the asymptotic behavior described by Eq. (4), with C ∝ T −1/2 being a temperature-dependent coefficient. This equation captures the temperature dependence of the resistivity at high magnetic fields. However, the resistivity seems to follow exp(CB) instead of Eq. (4), as shown in Fig. 4(d). Given the divergent susceptibility of MoS 2 at low density [20], we speculate that the spin polarization might contribute in increasing the magnetoresistance. As it was reported in Ref. [44] the mechanism that provides an enhanced magnetoresistance is related to the blocking of hopping due to spin polarization. Disentangling spin from orbital effects requires in-plane magnetic fields. However, the in-plane magnetic field does not couple to the electron spins in the K valleys of MoS 2 , because the spin is locked out-of-plane due to SO coupling [20,45].

III. CONCLUSION
The low-temperature resistivity of bilayer MoS 2 undergoes a transition from metallic to insulating temperature dependence at a critical density n c ≈ 1.7 × 10 12 cm −2 . This density is one order of magnitude larger than in silicon metal-oxide-semiconductor field-effect transistors [46] and three orders of magnitude larger than in GaAs quantum wells [47,48]. We attribute this transition to a disorder-induced transition in agreement with other metal-insulator-transition observed in disordered materials [46]. In fact, our observation is in line with the proposal of Klapwijk and Das Sarma [5], which states that the transition should be observed when the electron density is close to a few electrons per ionized impurity. The ionized impurities may be sulfur vacancies, which are known to have an inhomogeneous distribution with an average density of 1 × 10 12 − 1 × 10 13 cm −2 [49], thus the same order of magnitude as n c in our sample.
In the insulating phase, the resistance drops exponentially with increasing temperature, compatible with variable-range hopping. The limited temperature range considered in our experiments does not provide a conclusive distinction between Mott and Efros-Shklovskii laws. On the other end, the magnetic field dependence at low fields closely follows an Efros-Shklovskii law. The Coulomb gap is likely to appear in the density of states of MoS 2 , as the interaction parameter (r s ∼ 7) suggests that Coulomb energy is the dominant energy scale. However, the presence of a nearby metallic gate could contribute to the suppression of the gap in the density of state at large distances (r h > 2d) [39,[41][42][43]. Despite the screening effect, we observe the presence of the Coulomb gap by applying a perpendicular magnetic field. We interpret this result based on the wave function shrinking, which reduces long distant hopping.