Towards a DFT-based layered model for TCAD simulations of MoS2

Solid-State Electronics 197 (2022) 108437 Available online 23 August 2022 0038-1101/© 2022 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Towards a DFT-based layered model for TCAD simulations of MoS 2 L. Donetti, C. Marquez, C. Navarro, C. Medina-Bailon, J.L. Padilla, C. Sampedro, F. Gamiz Departamento de Electronica and CITIC, Universidad de Granada, Granada, Spain ARTICLE INFO The review of this paper was arranged by “S. Cristoloveanu” Keywords: Density functional theory MoS 2 TCAD Dielectric constant Density of states ABSTRACT In this work, we employ the results of atomistic DFT calculation to extract useful parameters for the simulation of few-layers MoS 2 structures with traditional TCAD tools. In particular, we focus on the charge distribution, which allows us to obtain a layered model for the dielectric constant, and on the effective densities of states in the conduction and valence bands taking into account the full 2D density of states. Using this model, we compute the capacitance of a metal–oxide–semiconductor structure and compare it to the one obtained employing a uniform model with averaged effective parameters. 1. Introduction The layered structure of 2D semiconductors such as MoS 2 presents several challenges for their accurate simulation employing standard TCAD tools. Apart from the need of thickness-specific material parameters [1–3], descriptions based on bulk parameters cannot accurately reproduce the layered spatial charge distribution and the 2D density of states. Therefore, even if the properties of MoS 2 are being successfully studied through ab initio atomistic methods [4], it is not straightforward to include this knowledge in classical TCAD tools [5]. In this paper, we model the layered structure of MoS 2 through a stack of alternating semiconductor and insulating layers representing the Van der Waals (VdW) gaps [6], with material parameters extracted directly from Density Functional Theory (DFT) calculations (Section 2). In particular, in Section 3, we obtain a layered dielectric model through the analysis of charge distribution in monoand few-layer structures with an applied out-of-plane electric field, and extract layer-dependent Nc and Nv values from the calculation of the full Density Of States (DOS). Then, in Section 4, we implement our model and compare the results of the layered model with those of a uniform model with averaged values of the parameters. Finally, we draw our conclusions in Section 5. 2. DFT calculations DFT calculations are performed with QuantumATK (version S2021.06) [7], employing its LCAO calculator with GGA exchange/correlation, PBE functional and Grimme DFT-D2 VdW correction. Spin–orbit interaction is included to obtain a better description of the band structure of the valence band around its maximum. We consider few-layer MoS 2 structures with a number of layers, N, between 1 and 10, keeping the experimental value of the in-plane lattice constant. For each structure in equilibrium, we compute the band structure that will be employed to obtain the DOS. Then, we apply an electric field perpendicular to the semiconductor layers by forcing a potential difference between the boundaries of the simulation cell. For the biased structures, we compute the induced electron density difference Δn, and the electrostatic potential difference ΔV, taking in-plane averages as a function of the out-of-plane position z [8]. 3. Layered model Some examples of DFT calculations with an applied bias are shown in Fig. 1. The behavior of the total charge density and the dipole moment density in each layer as a function of the external electric field suggests a layered dielectric model, where MoS 2 layers are separated by VdW gaps with a different dielectric constant, as the one shown in Fig. 2. By fitting the DFT results, the two dielectric constants and the thickness of the different layers can be extracted [8]. The resulting parameter values (see Table 1) are slightly different from those of Ref. [8] because here the spin–orbit interaction is included in the DFT calculation and a slightly different lattice constant is employed. The results suggest that the thickness of the external layers is different from that of the internal ones and also that an empty gap is needed to “fill” the thickness of the multilayer structure, as shown by the empty rectangles at the sides of the E-mail address: [email protected] (L. Donetti). Contents lists available at ScienceDirect Solid State Electronics journal homepage: www.elsevier.com/locate/sse https://doi.org/10.1016/j.sse.2022.108437 Received 4 July 2022; Received in revised form 12 August 2022; Accepted 15 August 2022 Solid State Electronics 197 (2022) 108437 2 layered structure in Fig. 2: the thickness of such gaps is given by d/2−ts/2, where d=6.15 Å is the interlayer distance, equal to half the c lattice constant of bulk MoS 2 . It is not straightforward to model the DOS of 2D materials through the usual 3D expressions involving the effective mass because of several reasons. First of all, it must be taken into account that the position of the Conduction Band Minimum (CBM) and Valence Band Maximum (VBM) in reciprocal space vary with the number of layers, N: the valley multiplicities are modified and also the corresponding band curvatures which define the effective masses change. For example, as we can see in Fig. 3, the gap is direct in the mono-layer case, while it becomes indirect for N>1. Then, even assuming that the material parameters can vary with N, the DOS of 2D bands presents a functional dependence on energy which is different from the one for 3D materials. As an example, in Fig. 4, we show the DOS, g(E), of monoand some fewlayers structure. In each case, the DOS is approximately constant near the extrema of each band and discrete steps can be observed when the minimum (or maximum) of different bands are reached. The narrow peaks near the jumps are artifacts of the numerical procedure (to improve the smoothness of the resulting DOS, we employ the tetrahedron method [9] for Brillouin zone integration, but some small fluctuations are still present). For non-degenerate bulk materials, the electron density, n, can be obtained by: n=Ncexp(−Ec−EF kT )(1) where Nc is an effective density of states in the conduction band, Ec is the CBM, EF is the Fermi level, k is Boltzmann constant, T is the absolute temperature. In this case, Nc is proportional to m3/2 DOS (where mDOS is the DOS effective mass), because of the energy dependence of Fig. 1. Excess electron density Δn (a, c) and electric field E (b, d) for monolayer (a and b) and 4-layer (c and d) MoS 2 , with 2 V bias. In (a) and (c), closed and open circles represent the z positions of Mo and S atoms, respectively. Fig. 2. Layered model for MoS 2 , showing the meaning of the different thickness parameters of Table 1 and the interlayer distance, d. Table 1 Parameters of the layered dielectric model of Fig. 2. ∊s ∊b ts (Å) ti (Å) 22.8 2.2 5.15 4.53 Fig. 3. Band structure of MoS 2 from mono-layer to 8 layers. The light blue areas denote the band gap and it is shown to help visualize the band extrema (both CBM and VBM are at K point for mono-layer, CBM at Q and VBM at Γ points for multi-layer structures). Energies are referenced to the Fermi level. Fig. 4. DOS of MoS 2 structures with different number of layers. L. Donetti et al. Solid State Electronics 197 (2022) 108437 3 g(E)∝(E−Ec)1/2. However, as we can observe in Fig. 4, g(E)shows neither a bulk nor a pure 2D behavior (whose DOS is a single step function). Therefore, we employ a numerical approach to compute the effective 2D DOS in the conduction band, Nc,2D, taking into account the whole DOS, with the following expression: Nc,2D=∫∞ Ec g(E)exp((E−Ec)/kT)dE (2) which is related to the 2D electron density by n2D= Nc,2Dexp( − (Ec−EF)/kT). The results of Eq. (2) at room temperature (T=300k) and an analogous one for Nv,2D are shown in Fig. 5(a). The jump between mono-layer and bi-layer is due to the fact that the position of the conduction band minima and valence band maxima change: the DOS reflects the change in the effective mass and in the degeneracies of the extrema position in the reciprocal space. For example, regarding Nc,2D, in Fig. 4 we can see that the height of the first step of the DOS of the conduction band is much smaller in the mono-layer case than in the multi-layer case. The density of states is larger for N>1 because the CBM is at Q point which has a larger multiplicity than the K point (6 vs. 2) and also because of a different band curvature. Then, the increase for larger N is caused by the fact that the subsequent bands get closer in energy to the first one. To obtain 3D effective densities of states Nc and Nv, we normalize Nc,2D and Nv,2D by the corresponding semiconductor thickness: Nc/v=Nc/v,2D ts+ (N−1)ti (layered model)(3) These bulk parameters show a pronounced decrease as N grows, as depicted in Fig. 5(b), because the increase of the 2D quantities (Fig. 5(a)) is much smaller than the increase of the thickness. Finally, we compute the temperature dependence of the effective densities of states. In Fig. 6(a), we represent normalized Nc and Nv as a function of T: for N=2 we can observe a linear dependence on T, as expected for 2D materials, and when N increases the behavior becomes more similar to the bulk one, with a T3/2 dependence (especially at larger temperatures). Indeed, in thick 2D systems, the sum of many step functions corresponding to the different sub-bands can give rise to an overall shape similar to the 3D case. This fact should not be surprising, because also other properties, such as the size of the band gap Eg, approach the bulk value already for N=8 [2]. A behavior corresponding to a relatively larger exponent is also observed for the monolayer case: this can be explained taking into account the two close conduction band minima which give rise to quite close steps in the DOS (Fig. 4). 4. TCAD simulations We implement the layered model for MoS 2 structures of different thickness in Sentaurus TCAD [5], with parameters taken from Table 1. For comparison, we also introduce a uniform model with an effective averaged dielectric constant [8]. The electron and hole densities are computed employing the values of Nc and Nv extracted in the previous section, which depend on the sample thickness. In the case of the uniform model, these values have been corrected to take into account the different effective semiconductor volume, that is: Nc/v=Nc/v,2D Nd (uniform model)(4) As shown in the previous section, the DOS behavior does not fully correspond to either 2D or 3D model, so that the values of Nc and Nv are directly set through the parameters Nc300 and Nv300, in eDOSMass and hDOSMass, respectively (with Formula =2). We first simulate isolated MoS 2 structures with an applied bias, in order to compare the electrostatic potential profile and electric field with the output of the DFT calculations of Section 2. The results for the case of a 4-layer structure and an applied bias of 2 V are shown in Fig. 7, exhibiting good agreement, at least for the layered model. Next, we simulate the capacitance of metal/oxide/semiconductor structures composed by a 1 nm SiO 2 layer and MoS 2 with different number of layers (inset in Fig. 8). The semiconductor is considered slightly doped (with ND equivalent to a 2D density of 1 ×1010 cm2), the metal workfunction corresponds to MoS 2 midgap, the bottom electrode is kept at Vbot =0V, while the top electrode bias, Vtop, is varied. In Fig. 8, we show the simulated capacitance compared to the one obtained for similar structures where MoS 2 layers are substituted by a uniform Fig. 5. 2D (a) and 3D (b) effective DOS of MoS 2 structures as a function of the number of layers N. Fig. 6. Temperature dependence of N c (a) and Nv (b), normalized to 1 at T= 150 K in log–log scale. The expected behavior for bulk materials (∝T3/2) and for strictly 2D materials (∝T) are also shown as dotted lines. L. Donetti et al. Solid State Electronics 197 (2022) 108437 4 material with averaged parameters [8]. At zero bias, the results obtained with the layered and uniform models are the same, while the results differ for large positive and negative values of Vtop . The main difference for all values of N are the smaller capacitance of the layered models for large absolute values of Vtop and the larger values of Vtop needed to reach the inversion/accumulation regions. These results can be explained by noting that the layered model of Fig. 2 includes small dielectric regions before the first semiconductor layer and after the last one. 5. Conclusions Employing DFT simulations, we have extracted dielectric constants and effective DOS in the conduction and valence bands for the implementation of a layered MoS 2 model to be employed in TCAD simulations. We have shown the relevance of taking into account the layered structure of few-layer MoS 2 by comparing the TCAD results with those obtained with a uniform model. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments The authors would like to thank the financial support of projects H2020-MSCA-IF-2019 Ref. 895322 (EU Horizon 2020 programme), TEC2017-89800-R (Spanish State Research Agency, AEI), Juan de la Cierva Incorporaci´ on Fellowship scheme 307 under grant agreement No. IJC2019-040003-I (MICINN/AEI). P18-RT-4826 (Regional Government of Andalusia) and B-TIC-515-UGR18 (University of Granada). Funding for open access charge: CBUA/Universidad de Granada. References [1] Mak KF, Lee C, Hone J, Shan J, Heinz TF. Atomically thin MoS2: a new direct-gap semiconductor. Phys Rev Lett 2010;105(13). https://doi.org/10.1103/ physrevlett.105.136805. [2] Kuc A, Zibouche N, Heine T. Influence of quantum confinement on the electronic structure of the transition metal sulfide T)§2. Phys Rev B 2011;83(24). https://doi. org/10.1103/physrevb.83.245213. [3] Kang Y, Jeon D, Kim T. Local mapping of the thickness-dependent dielectric constant of MoS2. J Phys Chem C 2021;125(6):3611–5. https://doi.org/10.1021/ acs.jpcc.0c11198. [4] Roldan R, Silva-Guillen JA, Lopez-Sancho MP, Guinea F, Cappelluti E, Ordej´ on P. Electronic properties of single-layer and multilayer transition metal dichalcogenidesMX2(M=Mo, W and X=S, Se). Annalen der Physik 2014;526(9–10): 347–57. https://doi.org/10.1002/andp.201400128. [5] Synopsys Sentaurus Device User Guide (T-2022.03); 2022. [6] Mirabelli G, Hurley PK, Duffy R. Physics-based modelling of MoS2: the layered structure concept. Semicond Sci Technol 2019;34(5):055015. https://doi.org/ 10.1088/1361-6641/ab121b. [7] Smidstrup S, Markussen T, Vancraeyveld P, Wellendorff J, Schneider J, Gunst T, Verstichel B, Stradi D, Khomyakov PA, Vej-Hansen UG, Lee M-E, Chill ST, Rasmussen F, Penazzi G, Corsetti F, Ojanper¨ a A, Jensen K, Palsgaard MLN, Martinez U, Blom A, Brandbyge M, Stokbro K, et al. QuantumATK: an integrated platform of electronic and atomic-scale modelling tools. J Phys: Condens Matter 2019;32(1):015901. https://doi.org/10.1088/1361-648x/ab4007. [8] Donetti L, Navarro C, Marquez C, Medina-Bailon C, Padilla J, Gamiz F. DFT-based layered dielectric model of few-layer MoS2. Solid-State Electron 2022;194:108346. https://doi.org/10.1016/j.sse.2022.108346. [9] Bl¨ ochl PE, Jepsen O, Andersen OK. Improved tetrahedron method for brillouin-zone integrations. Phys Rev B 1994;49:16223–33. https://doi.org/10.1103/ PhysRevB.49.16223. Fig. 7. Electrostatic potential (a) and electric field (b) compared between DFT calculation and TCAD simulations employing the layered model (solid curve) and the uniform model (dashed curve), for the 4-layer MoS 2 with an applied bias of 2V. Fig. 8. Capacitance of the structure depicted in the inset obtained with TCAD simulations, for the layered (solid lines) with different number of layers N and uniform (dashed lines) corresponding models. The oxide capacitance C ox is shown as a reference. L. Donetti et al.