A Gaussian Thermionic Emission Model for Analysis of Au/MoS2 Schottky Barrier Devices

Schottky barrier inhomogeneities are expected at the metal/TMDC interface and this can impact device performance. However, it is difficult to account for the distribution of interface inhomogeneity as most techniques average over the spot-area of the analytical tool, or the entire device measured for electrical I-V measurements. Commonly used models to extract Schottky barrier heights (SBH) neglect or fail to account for such inhomogeneities, which can lead to the extraction of incorrect SBH and Richardson constants. Here, we show that a gaussian modified thermionic emission model gives the best fit to experimental I-V-T data of van der Waals Au/p-MoS2 interfaces and allow the deconvolution of the SBH of the defective regions from the pristine region. By the inclusion of a gaussian distributed SBH in the macroscopic I-V-T analysis, we demonstrate that interface inhomogeneities due to defects are deconvoluted and well correlated to the impact on the device behavior across a wide temperature range from room temperature of 300 K down to 120 K. We verified the gaussian thermionic model across two different types of p-MoS2 (geological and synthetic), and finally compared the macroscopic SBH with the results of a nanoscopic technique, ballistic hole emission microscopy (BHEM). The results obtained using BHEM were consistent with the pristine Au/p-MoS2 SBH extracted from the gaussian modified thermionic emission model over hundreds of nanometers. Our findings show that the inclusion of Schottky barrier inhomogeneities in the analysis of I-V-T data is important to elucidate the impact of defects (e.g. grain boundaries, metallic impurities, etc.) and hence their influence on device behavior. We also find that the Richardson constant, a material specific constant typically treated as merely a fitting constant, is an important parameter to check for the validity of the transport model.

of van der Waals Au/p-MoS2 interfaces and allow the deconvolution of the Schottky barrier heights of the defective regions from the pristine region. By the inclusion of a gaussian distributed Schottky barrier height in the macroscopic I-V-T analysis, we demonstrate that interface inhomogeneities due to defects are deconvoluted and well correlated to the impact on the device behavior across a wide temperature range from room temperature of 300 K down to 120 K. We verified the gaussian thermionic model

I. Introduction
Schottky barrier inhomogeneities are present at metal/transition metal dichalcogenide (TMDC) interface due to potential fluctuations and defects in the material and this can impact the device performance [1][2][3][4]. However, the most commonly used thermionic emission transport model modified with a simple ideality factor correction [5] does not account sufficiently for non-ideal effects such as inhomogeneities, often leading to the extraction of an apparent Schottky barrier height (SBH) convoluted with other factors such as defects (inhomogeneity) and temperature.
The extracted apparent SBH from the simple model does not represent the true band alignment at the metal/semiconductor interface and is counterproductive to the correct understanding of the energetics of the interface. Several other models have been proposed to account for these non-ideal effects in the current-voltage (I-V) behavior of Schottky barrier devices [1][2][3][4][5][6][7][8]. In this paper, we compare the effectiveness of four types of commonly used transport models to extract correct values of the SBH and the material specific Richardson constant of MoS2. We verified the SBH extracted from the transport models against the SBH obtained from ballistic hole emission microscopy (BHEM) [9][10][11], which a direct measurement method for the SBH at the nanoscale, and verified the Richardson constants extracted from the models with theoretical calculated values based on the electron (hole) effective mass of MoS2 [12]. The accurate extraction of the SBH depends greatly on the successful determination of Is from the forward bias slope of the experimental I-V curves (Eq. 1). However, the thermionic emission model does not fit well to experimental I-V curves at low temperatures, and the modified (non-ideal) thermionic emission models [6,7], the thermionic field emission (thermally assisted tunneling) model or the generationrecombination model are often used instead [13,14]. While these modified models seem to fit certain experimental datasets well, the numbers that have been extracted have not always been reliable. For instance, negative SBHs have been reported [15,16] [23][24][25], where the Schottky barrier at the metal/semiconductor contact is modulated by the gate bias, the correct analysis of the Schottky barrier is crucial to the understanding of the subthreshold behavior of these 2D transistors, especially in the presence of defects [26][27][28].

II. Methods to extract the Schottky Barrier Height (SBH)
The thermionic emission model, Eq. 1, predicts that for V > 3kT/q, a plot of In I against V will be linear with a slope of 1 and its intercept at V = 0 will give Is. From Is, the Schottky barrier height can be extracted. A direct reading of Is from the experimental I-V curve is typically not used as the experimental reverse biased saturation current, as it also contains the image force lowering and other minority carrier effects [5]. However, the thermionic emission model is inadequate for a realistic metal/semiconductor interface especially at low temperatures. To account for the nonideal transport mechanisms, and series resistance in real devices, the ideality factor and the series resistance Rs, are empirically added to the model [5,29] and Eq. 1 becomes: Using the modified thermionic emission model, diode parameters such as the ideality factor n and barrier height ф can be plotted and were found to be dependent on the temperature. At low temperatures, the thermionic emission model does not fit well to the experimental data and increases greatly beyond 1, signifying non-ideal transport.
While the non-ideal transport has been attributed to additional current contributions from thermally assisted tunneling (thermionic field emission) across the Schottky barrier, generation-recombination current in the depletion region and image force lowering, it is not clear how the empirically modified thermionic emission model can account for these effects. From these modifications, a few versions of the Richardson plot will be analyzed.

A. Method 1: Ideal Richardson plot (In Is/T 2 vs 1/T)
This is the simplest method and is derived directly from the saturation current term of the thermionic emission model, Eq. 2. When temperature dependent plots can be obtained, a plot of ln(IS/T 2 ) vs. I/T, called a Richardson plot, will be a straight line where the slope and intercept at 1/T = 0 will give ф, and A** respectively. The empirically added ideality factor is not included in Eq. 2, hence this method is the ideal Richardson plot analysis.

B. Method 2: Ideality factor modified Richardson plot I (In Is/T 2 vs 1/nT)
To account for effects that cause deviations from ideal (n = 1) behavior, such as the image force and surface charges, which they argued to be also present at zero bias, Hackam and Harrop proposed a modified Richardson plot from Eq. 3 to include the ideality factor in the Is term [8]. The forward current now looks: The addition of n to the Is term in Eq. 4 now gives a modified Richardson plot Eq. 5 from which the SBH can be extracted from the gradient of the straight line and A** from the y-intercept.

C. Method 3: ideality factor modified Richardson plot II (nIn Is/T 2 vs 1/T)
Bhuiyan, Martinez and Esteve found that the Hackam and Harrop model does not work for them due to the presence of a strongly temperature dependent SBH and ideality factor measured [6] and that the A** extracted from using the Hackam and Harrop method is too large. Hence, they empirically proposed Eq. 6: Following their modification, the modified Richardson plot now reads:

D. Method 4: Inhomogeneous Gaussian barrier modified Richardson plot
Two different inhomogeneous Schottky barrier models have been proposed independently by Werner and Güttler [3,4], and Tung and coworkers [2,30]. Werner and Güttler used a gaussian approximation of the SBH distribution to account for the potential fluctuations at the interface, while Tung used a generalized model. While Tung's model is more rigorous, Werner and Gutter's model is simpler and can be placed into the context of BHEM measured SBHs. Hence in this paper, we focused on the Werner and Güttler model of gaussian SBH [4], which is given by: where is the apparent Schottky barrier height obtained as a result of the convolution of the gaussian distributed SBH with temperature in the thermionic emission model, Φ is the mean Schottky barrier height and σ is the standard deviation of the gaussian distribution. To obtain the σ of the Gaussian, a plot of against 1/T can be used. The gaussian standard deviations (σ) extracted from Eq. 8 can then be used to correct for the gaussian distributed SBH to obtain the gaussian corrected Richardson plots Eq. 9.
Here, a plot of ln ( 2 ) − 2 2 2 2 2 against 1 will give the A** in the intercept and Φ in the gradient.
The temperature dependence of the ideality factor and SBH, initially viewed as empirical inconsistencies in many experiments, is now well explained by Werner and Güttler to arise from the inhomogeneous SBH and that capacitance-voltage (C-V) measurements give Φ. In our experiments, we did not perform C-V measurements as capacitance measurement is not typically used in the operation of devices, but the current as a function of applied voltage (I-V) is used and is more common for analysis.
Werner and Güttler also showed that for lightly doped (10 15 to 10 17 cm -3 ) semiconductors, thermionic emission dominates carrier transport, even at low temperatures down to 77 K.
We demonstrate in our Au/MoS2 sample that by using a gaussian distributed SBH to account for these inhomogeneities, a more reliable value of the SBH can be obtained. We verify the gaussian thermionic emission model systematically using two different types of MoS2 (geological and synthetic crystals, from different suppliers) [31,32] by performing temperature dependent current voltage measurements

IV. Experimental
The p-type geological MoS2 crystal was obtained from Ward's Science [31], the p-type synthetic MoS2 crystal (intrinsic) was obtained from 2D semiconductors [32].
The ohmic contacts to the MoS2 crystals were deposited in a high vacuum e-beam evaporator system (AJA International) after cleaving off the top surface using sticky tape to obtain a fresh surface for the evaporation of Ti (5 nm Figure 2 shows

Methods 1, 2 and 3.
To analyze the temperature dependence of the SBH and ideality factor, Figure 3 shows the Richardson plots analyzed using used Methods 1, 2 and 3.
By plotting the temperature dependent Richardson plots, one can extract the effective Richardson constant and the SBH of the device from the gradient and the y-intercept respectively. However, the ideal Richardson plots and the ideality factor modified Richardson plots [6,7] yield A** values which are orders of magnitude away from the theoretical values and expected values valid for the model, or give unrealistic SBH values. Table 1 [4], which is given by Eq. 8. Plotting against q/2kT in Figure   4 shows the presence of two different regimes where the two diode parameters dominate.
We collected a few datasets over a few different locations on the same sample and   Table 2. Schottky barrier heights of Au/geological p-MoS2 and Au/synthetic p-MoS2 measured using temperature dependent currentvoltage (I-V-T) measurements and ballistic hole emission microscopy (BHEM). Defects 1 and 2 denote the two Schottky barrier heights obtained from the double gaussian distribution model of Schottky barrier heights that could arise from different kind of defects.

VI. Discussion
The variation across the same sample due to impurities [26,27], which suggests that these impurities are important contributors to the electrical behavior of the devices [28]. We have also observed similar defects and impurities in our samples using STM imaging (Supplemental Information, Figure S3), but similar to their analysis, we are unable to to identify the chemical composition of the defects in STM/BEEM and I-V measurements as these techniques does not have chemical sensitivity. We have shown that the double gaussian model proposed in this paper can be used to deconvolute the defective SBH from the pristine SBH. This method is also useful for research on engineering ohmic contacts to these materials.
Our MoS2 crystals are p-type crystals. In recent literature, most MoS2 crystals reported are n-type semiconductors [37][38][39][40][41][42]. However, early reports have noted that some geological samples are intrinsically p-type [43][44][45] and for synthetic crystals, depending on the growth method, the TMDC crystals can be unintentionally doped ntype or p-type due to doping from the chemical transport agent or impurity inclusion [46,47]. The Gaussian distributed model should also work for n-type  Figure S4), suggesting that the strain did not play a significant role in our measurements. Figure 6 shows the STM image of epitaxial Au films on MoS2 grown using slow deposition method which is an indirect evidence of a clean abrupt interface Au/MoS2 [17,19]. Our results lend support to the presence of an unpinned Fermi level for a well prepared van der Waal's Au/MoS2 interface that is deposited slowly (~0.2 Å/s), consistent with a few recent reports [18,19], and an old photoemission study [49]. We note that the gaussian distributed SBH analysis method we propose in this paper has been applied also to monolayer MoS2 devices [50]. In their report, Moon et al. analyzed the top and edge contact of Au/MoS2 n-type FET devices using the gaussian distribution model at different gate bias. They observed the top contact has a larger SBH and larger σ, showing more inhomogeneity than in the edge contact, which has a lower SBH and lower σ. However, they did not use the gaussian corrected Richardson plot, but they indicated that the standard Richardson plot is not valid due to an observed temperature dependence of the SBH. Our method of analysis can be used to bridge the gap between real 2D devices and theoretically proposed models [23] by deconvoluting SBH inhomogeneity from intrinsic material transport behaviors.

VII. Conclusion
We We used BHEM, which is a more tedious but direct method to measure the zero bias SBH without the need to rely on the validity of temperature dependent models, and to experimentally validate the importance of including a gaussian distributed SBH at the nanoscale in conventional I-V-T analysis. Our results provide the basic framework for extracting the pristine SBH from temperature dependent I-V data and demonstrate that with careful use of the dual parameter (A** + SBH) analysis, we avoid obtaining unphysical numbers that are counter-productive for understanding such interfaces. This implies that the I-V-T analysis can yield important insights on the SB inhomogeneities even though it might be a macroscale measurement.

VIII. Acknowledgements
We acknowledge funding from the A*STAR Pharos Grant No. 1527000016.

IX. Supplemental Information
Description of the Werner method to correct for the high series resistance found in the

Werner method
In the main text, we mentioned that the Werner method can be used to correct for the high series resistance of the Au/MoS2 diodes. The Werner method is described in detail here in the context of our experimental data.
Under forward bias and with series resistance contribution, the voltage across the diode, Vd = V -IRs and Vd >> kT, the thermionic emission current (Eq. S1) is given by the simplified form: where Id is the thermionic diode current. Differentiating Eq. S1 gives the small signal conductance G = dId/dV and one obtains: Werner showed that by plotting G/Id against G, named hereafter as the Werner plot, will give a straight line with y-axis intercept of q/nkT where n can be extracted, and xaxis gives the intercept of 1/Rs.  Figure S2 shows that for our experimental data, the P-L model fits better than the B-K model in the fitting range of 0.3 V to 1.3 V, about 0.4 V above the SBH. Therefore, the P-L model is used in our analysis. of Au/geological MoS2 (674 points) to increase signal to noise ratio, spectrum fits better to the P-L model (inset).

Defect imaging using scanning tunnelling microscopy
In the main text, we mentioned that we have also seen defects that were previously reported by Addou et al. in our geological MoS2 samples. [4,5] Figure S3 shows the STM images that we obtained from our MoS2 samples. Figure S3a and Figure S3b shows the same region scanned at positive and negative biases respectively. We  To study the effect of strain at our Au/MoS2 interface, we used Raman spectroscopy to compare the relative shifts of the vibrational modes. Figure S4 shows the Raman spectrum obtained on our Au/synthetic MoS2 device using a 532 mn laser at 120 K using a liquid nitrogen flow stage (Linkam). No shift of the Raman peaks was detected when we collected the spectrum on the bare MoS2 surface and on the Au/MoS2 although a decrease in peak intensity can be detected due to attenuation by the 15 nm Au layer. The Raman peaks the all match within 0.2 cm -1 which is the resolution limit of our Raman spectrometer. Hence, we can conclude that we did not detect any strain in our Au/MoS2 devices.  Figure S4. Raman spectrum collected on the Au/MoS2 surface and the bare MoS2 surface at 120 K using liquid nitrogen cooling and 532 nm laser excitation.