Surface-gate-defined single-electron transistor in a MoS2 bilayer M. Javaid,1, ∗ Daniel W. Drumm,1, 2 Salvy P. Russo,1 and Andrew D. Greentree1, 2 1 arXiv:1608.07894v1 [cond-mat.mes-hall] 29 Aug 2016 Chemical and Quantum Physics, School of Sciences, RMIT University, Melbourne VIC 3001, Australia 2 The Australian Research Council Centre of Excellence for Nanoscale BioPhotonics, School of Sciences, RMIT University, Melbourne, VIC 3001, Australia (Dated: August 30, 2016) We report the multi-scale modeling and design of a gate-defined single-electron transistor in a MoS2 bilayer. By combining density-functional theory and finite-element analysis, we design a surface gate structure to electrostatically define and tune a quantum dot and its associated tunnel barriers in the MoS2 bilayer. Our approach suggests new pathways for the creation of novel quantum electronic devices in two-dimensional materials. I. INTRODUCTION Novel two-dimensional (2D) materials such as graphene and transition-metal dichalcogenides (TMDCs) have attracted significant interest due to their unique electronic and optical properties [1–3]. The ultrathin geometry and dangling-bond-free interfaces of 2D materials make them good candidates to integrate on various substrates [4]. Among these 2D materials, MoS2 (a TMDC) is promising for quantum electronics; it possesses interesting layer-dependent properties. For example, by decreasing the thickness of MoS2 from bulk to a single layer, its band gap switches from indirect to direct and increases by more than 0.6 eV, leading to strong photoluminescence from a single layer [2, 5]. The single-electron transistor (SET) is a device where electrons tunnel one by one to and from a small island through tunnel barriers [6]. The tunneling of an electron and the quantization of charge are controlled by gate electrodes. SETs have extensive applications in nano-electronic devices [7]. They have been proposed as a future alternative to conventional CMOS transistors [8], and have been used as nano-scale electrometers [9], capable√of measuring sub-electron charge variations to 10−6 qe / Hz. They also have many applications in single-electron logic circuits [10] and single-electron turnstile devices [11]. Many different techniques are used to fabricate SETs based in materials such as Si [12–14], GaAs [15] and carbon nanotubes [16]. There is a significant body of literature on fabricating SETs. For example: Thelander et. al. have designed gold-nanoparticle SETs using carbon nanotubes as leads at 200 K [17], Klein et. al. reported a colloidal chemistry technique to create cadmium-selenide nanocrystals of varying sizes [18], and Kim et. al. used a focused-ion-beam technique to fabricate SETs operating at room temperature [19]. Recently, there has been extensive work to define quantumdot devices in 2D materials including graphene [20] and ∗ maria.javaid@rmit.edu.au TMDCs [21]. Surface gates have been used to create tunable tunnel barriers in Si interfaces [22], but to date there have been no reported designs of electrostatically tunable tunnel barriers and tunable quantum dots in 2D MoS2 . It is known that the MoS2 monolayer band structure shows no response to applied electric fields [23]. In contrast, MoS2 bilayers exhibit significant band-structure modification under perpendicular electric field [23, 24]. Here we present a lithographically appealing design of a surface-gate-defined SET in a MoS2 bilayer, and model its physical characteristics. The structure of our proposed device is a MoS2 bilayer sandwiched between a HfO2 substrate and cap layer, with metallic surface gates, as shown in Fig. 1(a). We use the surface-gate electric fields to modify the local band structure in the MoS2 bilayer. We present a fully self-consistent design of the device obtained through multi-scale modeling. Recently, a more conventional gate design has been used to demonstrate quantum confined structures on a few layers of WSe2 with tunnel barriers defined by electric fields at a temperature of 240 mK [21]. We modeled a MoS2 bilayer using density-functional theory (DFT) to study the effects of electric field on the band structure of MoS2 bilayers. We performed numerical simulations to calculate the potential at the MoS2 bilayer due to the surface gates and used this potential to study the modification in the lowest unoccupied molecular orbital (LUMO) of the MoS2 bilayer. We modeled our device geometry (obtained from LUMO bending in the MoS2 bilayer) with the commercial finite-element analysis simulation tool comsol [25] to obtain the selfcapacitance of the SET island and the capacitances between the island and each of the electrodes. We used these capacitances in numerical simulations of transport through the SET at a temperature of 1 K and studied its physical characteristics. We iteratively cycled through comsol and numerical simulations to obtain a consistent picture of the device. This paper is organized as follows: we begin our discussion with the DFT modeling of the MoS2 bilayer followed by the comsol modeling. Then we present the influence of the potential due to surface gates in the MoS2 bilayer plane and numerical transport simulations of the resulting SET device. 2 (a) G1 G3 G4 G2 Cap layer Substrate (b) (c) CisR(Vg ) CidR(Vg ) Island Cig Vs Vd Vg −0.16 −0.22 0 −200 −200 (e) −0.28 Potential (V) Potential (V) X (nm) (d) 200 −0.28 0 Y (nm) 200 −0.34 −0.3 Island Source Drain −0.32 −0.34 −200 0 Y (nm) 200 (f ) X (nm) Z (nm) Y (nm) FIG. 1. (Color online) (a) Schematic showing a MoS2 bilayer sandwiched between HfO2 cap layer and substrate, with metallic surface gates above the cap layer. Gates G3 and G4 are shown separated from G1 for visibility, but are actually collinear. (b) Perspective view of the black layer in (a) showing A-A′ stacked, 2H MoS2 bilayer, Mo (large, green) atoms, S (small, yellow) atoms. (c) Equivalent circuit diagram of SET showing the capacitance/resistance between the island and the source/drain and capacitance between the island and the top gates. (d) Top view of the potential created on top of the the MoS2 bilayer sheet, showing the creation of well-defined source, drain, and island regions. (e) Potential slice along bilayer directly underneath gates G1, G3, and G4, showing source, island, and drain for a surface-gate voltage of 0.27 V. (f) 3D view of the geometry used in the comsol modeling. The (blue) highlighted regions are the source, island and drain in the MoS2 bilayer plane. Lengths are in nm. II. DFT MODELING OF MOLYBDENUM-DISULPHIDE BILAYER STRUCTURES Although the effects of electric field on the band structure of MoS2 bilayers are available in the literature, the focus has either been on behavior over extreme field strength ranges (with sampling too coarse for our device) [23] or on the band-gap response without specific discussion of the LUMO physics [24]. Thus, to achieve the insight necessary to design our device, we explore the LUMO physics directly in the target operating regime. We investigated the band structure of an optimized, A-A′ stacked (Mo atoms in the top monolayer above the S atoms of the bottom monolayer), 2H-phase MoS2 bilayer for varying electric fields. A-A′ is the most-reported stacking order for bulk MoS2 and it is very close in energy to the most stable stacking order A-B [24]. A perspective view of the MoS2 bilayer is shown in Fig. 1 (b). To model the effects of electric field on the MoS2 bilayers, we used DFT in crystal09 [26, 27]. MoS2 has hexagonal symmetry and belongs to the P 63 /mmc space group, with lattice parameters a = 3.17 Å and c = 12.324 Å [28]. We created the bilayer unit slab by cutting a (001) plane from a bulk MoS2 model, and including vacuum to a total cell height of c = 500 Å. The exchange and correlation terms were described via the PBEsol functional [29], which generally predicts lattice constants more accurately than PBE and LSDA (thus improving the equilibrium properties of solids), and also handles the electronic response to potentials better than most GGA functionals [29]. We used Gaussian basis sets; Mo SC HAYWSC311(d31)G cora 1997 [30] for Mo atoms (a Hay-Wadt effective-core pseudopotential [31] combined with a valence-electron basis set) and S pob TZVP 2012 [32] for S atoms (an all-electron basis set). We set a 8×16×1 Monkhorst-Pack [33] k-point mesh. We optimized the geometry of a MoS2 bilayer unit slab under zero electric field and then used this optimized geometry to study the effects of electric field on the band structure of MoS2 bilayers. We calculated the band structure along the high-symmetry path Γ-M-K-Γ. Band structures (shown in the Appendix) were calculated for varying electric fields applied perpendicular to the MoS2 bilayers. The responses of the LUMO and the highest occupied molecular orbital (HOMO) of the MoS2 bilayer to applied electric fields are shown in Fig. 2(a). The HOMO remains static with electric fields of up to ±0.2 V/Å (in agreement with [23]) while there is a significant response in the LUMO. Fig. 2(a) also shows some inconsistent points in the LUMO and HOMO under electric field. We excluded these points from further calculations as the band structure at these points is qualitatively different and we believe that these inconsistent points do not represent the underlying physics. We applied a linear fit to the filtered LUMO and HOMO data, and set zero energy at the LUMO intercept (at zero field). As an electric field is applied up to ±0.2 V/Å, the LUMO bends through 138 meV. Thus electron confinement is achievable in a MoS2 bilayer via electric-field modification of the LUMO, which occurs at a rate of 690 meV / (V/Å). III. CIRCUIT SIMULATIONS IN COMSOL We modeled our device geometry in comsol [25] to get the capacitances between all the device components. The 3 to the total capacitance of the island. The ground plane is 11 nm (thickness of the HfO2 substrate) below the island. The wire-frame rendering of the geometry used in comsol modeling is shown in Fig. 1(f), where all the lengths are in nm. For an island of diameter 12 nm, the self-capacitance is CiΣ = 8.45 aF which corresponds to a charging energy of e2 /CiΣ = 18.9 meV, where CiΣ is defined as (a) LUMO Energy (eV) 0 −0.7 CiΣ = Cis + Cid + Cig1 + Cig2 + Cigp (1) HOMO −1.4 0 0.1 Abs. electric field (V/Å) 0.2 (b) 0.6 Island 0.45 0.45 0.3 0.3 Source Drain 0.15 0.15 0 0 Barrier Barrier −0.15 −0.3 −200 −0.15 −100 0 Position (nm) 100 Relative LUMO bending (meV) Relative electric Field (mV/Å) 0.6 −0.3 200 FIG. 2. (Color online) (a) Responses of the MoS2 bilayer LUMO (grey circles) and HOMO (blue squares) to applied electric field, with first-order fits to LUMO (red line) and HOMO (black line) data. Negative electric-field values are folded over to the positive axis. (b) Electric field (solid blue line) relative to the dashed line in (a), created across the MoS2 bilayer for 0.27 V surface-gate voltage. This electric field is used to calculate the LUMO bending (dashed red line) via the line-fit parameters from (a). The electric field and the LUMO bending corresponding to the intersection of dashed line with LUMO fit in (a) is set at zero. The Fermi level (black dashed line), barrier and island regions (blue dashed arrows), and region boundaries (black vertical bars) are marked. geometry (obtained through potential simulations discussed later and converged by iteratively cycling through numerical potential and transport simulations and comsol) used for comsol modeling is an island of diameter 12 nm and the same thickness as the MoS2 bilayer. There are 16 nm wide MoS2 barriers between the island and inplane source and drain. The surface gates are modeled as two intersecting wires each of 5 nm radius sitting atop the 7 nm HfO2 cap layer. We excluded gates G3 and G4 of Fig. 1 (a) from our comsol modeling as they are collinear with gate G1 and do not contribute significantly Here Cis = −0.252 aF is the capacitance between the island and the buried source, Cid = −0.253 aF is the capacitance between the island and the buried drain, Cig1 = −2.2 aF is the capacitance between the island and gate G1, Cig2 = −2.38 aF is the capacitance between the island and gate G2, and Cigp = −3.35 aF is the capacitance between the island and the ground plane. IV. SINGLE-ELECTRON TRANSPORT SIMULATIONS In the design of our proposed MoS2 bilayer SET, there are two longer surface gates G1 and G2, perpendicular to each other. Two shorter gates, G3 and G4, are collinear with the G1 gate as shown in Fig. 1(a). All of these surface gates are electrostatically insulated from each other; this can be done by thermally growing a thin oxide layer between them. The equivalent electrical circuit diagram is shown in Fig. 1(c), where R(Vg ) is the tunnel resistance between the island and the source/drain. This resistance is a function of the gate voltage Vg applied to the surface gates. Vs and Vd are the source and drain potentials applied directly to the buried source and drain in the MoS2 bilayer plane (created by the surface gates). We modeled our device design through in-house numerical simulations. Gates G3 and G4 were placed on top of gate G1, with their finite edges 32 nm apart from each other. The thicknesses of the cap layer and the substrate were 7 nm and 11 nm respectively. A varying DC potential of 0.2 V to 0.34 V was applied to all surface gates. All of these parameters and the DC potential range were finalized through iterations between the comsol modeling and the numerical simulations. The top view of the potential created in the MoS2 bilayer plane due to the surface gates is shown in Fig. 1(d). The intersection of gates G1 and G2 generates a well-defined island in the MoS2 bilayer plane. Gates G1, G3, and G4 define the source/drain wiring configuration and tunable tunnel barriers in the MoS2 bilayer plane. An analytical expression for the potential created due to one of the surface gates is given in the Appendix. Fig. 1(e) shows a slice of the potential on top of the MoS2 bilayer plane and underneath gates G1, G3, and G4. Note that there is a small local minimum in the potential (and the electric field) before the barriers. The (a) RT (arb. units) 4 7,500 5,500 3,500 0.2 (b) 0.24 0.28 Vg (V) 0.32 Current (arb. units) 0.02 electron tunnels to the island and current flows through the SET which is shown by the peaks in Fig. 3(b). On the other hand, when none of the island states are available to the incoming electron, it cannot tunnel to the island and the current through the SET drops to zero, defining the Coulomb-blockade region between the peaks. The height variation between the current peaks arises from the changing response of the MoS2 -bilayer density of states in the source/drain and island, in contrast to an ideal metallic SET. V. 0.01 0 0.2 0.24 0.28 Vg (V) 0.32 FIG. 3. (a) Variation of tunnel resistance with gate voltage, Vg . (b) Characteristic Coulomb blockade peaks for a bias voltage of 100 µV and at a temperature of 1 K, showing a variable current through the MoS2 bilayer SET. potential difference across the MoS2 bilayer is proportional to the electric field (Fig. 2(b)). We used the linear fit to Fig. 2(a) to obtain the LUMO bending at fields shown in Fig. 2(b). This LUMO bending as a function of the position across the MoS2 bilayer plane is shown in Fig. 2(b). With the above dimensions of the device, a 12 nm diameter island and 16 nm wide tunnel barriers are created in the MoS2 bilayer plane. The dimensions of the SET are comparable with experimental structures in GaAs/ AlGaAs [34]. The size of the island and the width of the tunnel barriers can be controlled through the substrate and the cap layer thicknesses. The LUMO bending is used as an input to numerical transport simulation. Fig. 3(a) shows the variation of tunnel resistance with Vg . Increases to the gate voltage cause the barrier height to rise, reducing the transmission coefficient. Due to scaling of input electron energy, there is a minimal change in transmission coefficient with gate voltage but density of states in both source/drain and island increases predominantly and that leads to the observed decrease in tunnel resistance. (For details of the tunnel resistance and transport calculations, see the Appendix) The behavior of the SET can be described by varying Vg . Fig. 3(b) shows the characteristic current passing through the SET at a temperature of 1 K and a bias voltage of 100 µV (In this paper, we are only considering the low-bias regime.) By varying the surface gate voltage, the overall potential landscape through the MoS2 bilayer changes, thus modifying the barrier height and alignment between incoming electron energy and island states. When the energy of an incoming electron resonates with an unoccupied energy level in the island, that CONCLUSIONS The surface gate approach to quantum devices is quite flexible and amenable to scalable integrated devices. With identical gate voltages, we already have significant control over barrier height, the energy states in the island and hence, electron tunneling. One could achieve further independent control of the incoming electron energy or island states by applying different voltages to each gate. Moreover, due to the small size of the island, our preliminary modeling holds promise for relatively hightemperature single-electron devices. The most important feature of the gate design shown here, as opposed to more conventional gate-defined quantum dot designs (such as, for example, [21]) is that the creation of potential wells directly underneath and in the near field of the metallic surface gates should allow moreprecise lithographic definition of the quantum features in the active layer. In the conventional approach to creating sub-surface quantum dots, the relatively large distance from the gate layer to the active layer gives rise to a complicated potential landscape that must be solved for self-consistently with Schrödinger-Poisson solvers. Conversely, as our active layer is in the near field of the surface gates, the potential should more precisely follow the metallisation on the surface layer. This in turn means that more complicated integrated circuits should be realisable, for example multi-island dot structures or integrated qubit/readout type circuits. This design ethos is exemplified by the surface gate defined SETs fabricated at a silicon/silica interface by Angus et al. [22]. Here we have studied an embedded MoS2 bilayer system. However, we would expect our results to be easily applicable to other bilayer two-dimensional materials that exhibit a LUMO response to electric field (e.g. MoSe2 , MoTe2 and WS2 [23]), with only minimal modification of the applied surface potentials. VI. ACKNOWLEDGEMENTS MJ and ADG acknowledge funding by ARC Discovery Grant No. D130104381. DWD acknowledges the support of the ARC Centre of Excellence for Nanoscale BioPhotonics (CE140100003). This work was supported by computational resources provided by the Australian 5 (a) Government through the National Computational Infrastructure under the National Computational Merit Allocation Scheme. E-Ef (eV) 2.72 0 −2.72 Here we present a brief overview of the design of the MoS2 bilayer single-electron transistor, followed by details of our modeling. We show our analytic analysis of the electric field across the bilayer plane, the response of the MoS2 bilayer band structure to electric fields, the spatial behavior of the LUMO under different gate voltages, and details of our calculations of the resistance of (and current through) the SET. Gates G1 and G2 of Fig. 1(a) are approximated as infinitely long lines of charge, the intersection of which defines our island via the doubly-strong steep potential. To establish source and drain leads, we use gates G3 and G4 to create similarly strong potential landscapes directly underneath gate G1. For convenience, gates G3 and G4 are approximated as one infinite line of charge with a missing, finite segment centred on the intersection of gates G1 and G2. For our treatment, this dictates V3 = V4 , but that condition is not necessary if the semi-infinite gates G3 and G4 are treated separately and identical source and drain behavior is not required. It should also be noted that V2 ≈ V3 = V4 achieves comparable island and lead depths below the barrier peaks; V2 = V3 = V4 is not a constraint. Gate G1 is therefore effectively a plunger gate controlling the depth of all components, which could in principle be operated in isolation from the rest (although V1 +V3 > V2 is required to avoid an unintended perpendicular lead under gate G2). For simplicity here, however, we set V g = V 1 = V2 = V3 = V 4 . Defining x as the distance in the bilayer plane transverse to gate G1 (with zero value directly underneath G1), (y similarly with reference to G2), and z as the bilayer-ground plane separation, the electric potential difference created across the MoS2 bilayers is the sum effect from all gates, which is analytic. For example, the potential due to gate G1 is:  2  λ x + (z − d)2 V = log 2 , (A1) 2πǫo ǫr x + (z + d)2 where λ is the linear charge density of gate G1, ǫr is the dielectric constant of HfO2 [35] and d is the distance of gate G1 from the ground plane. The other potentials are similarly derived, and the full potential is the sum effect from all of the gates. The finite difference between the top and bottom edges of the bilayer is evaluated to provide the field strength across the bilayer. Having obtained a two-dimensional map of the potential difference across the bilayer plane – e.g., Fig. 1(d) in the main paper – we then consider its effect on the MoS2 band structure at fields from 0 to 0.2 V/Å in steps no larger than 0.01 V/Å, and smaller as required to verify −5.44 EF=0.1 V/Å 5.44 2.72 E-Ef (eV) Appendix A (b) EF=0 V/Å 5.44 0 −2.72 Γ M K Γ −5.44 Γ M K Γ FIG. 4. (a) Band structure of MoS2 bilayer under zero electric field. There is an indirect band gap of 1.30 eV. (b) Band structure of MoS2 bilayer under finite electric field. There is an indirect band gap of 1.23 eV. linear behavior. Fig. 4 shows two sample band structures. We are particularly interested in the minimal substrate footprint, and the constraint of reasonable applied field strength. Therefore in the main paper, without loss of generality in the method, we explicitly consider the high charging energy case, when the classical occupancy of the SET island can change by at most one electron at a time. Thus, the extent of the LUMO bending created must match well to the charging energy of the system; ∆L ≈ E = 1 T −1 Q C Q, 2 (A2) where Q is the charge on the SET, defined by Q=[Qi Qs Qg1 Qg2 Qd Qgp ]T . C is the capacitance matrix between all the device components and is defined as:   CiΣ Cis Cig1 Cig2 Cid Cigp  Csi CsΣ Csg1 Csg2 Csd Csgp    C Cg1Σ Cg1g2 Cg1d Cg1gp  C C =  g1i g1s  (A3)  Cg2i Cg2s Cg2g1 Cg2Σ Cg2d Cg2gp   C  di Cds Cdg1 Cdg2 CdΣ Cdgp Cgpi Cgps Cgpg1 Cgpg2 Cgpd CgpΣ the terms of which are obtained via comsol modeling as described in the main paper. The tunnel resistance cannot be evaluated exactly. The usual approach to an arbitrary barrier, the WKB approximation, is not directly applicable if the incoming energy of the electron matches the barrier energy at any point. In fact, not only does the approximation not hold at those points, but also in their vicinities. The connection formulas, using an Airy function ansatz to patch difficult regions, are the standard way to circumvent this issue; however, the conditions for their use are not met by this barrier. Accordingly, we make the approximation that the barrier is finite and rectangular, with the identical width (16 nm) and height equal to the peak value of LUMO bending, called Lpeak while setting L(−200) = 0. Since such a barrier is definitively larger, and no thinner, than the original barrier, its transmission coefficient is guaranteed to be lower than any experimental device. 6 Island LUMO bending (meV) 0.45 Vg =0.34V Vg =0.305V Vg =0.27V 0.3 Vg =0.235V Vg = 0.27 V, where µ is set such that ε is equal to the black-dashed line. There, the island has 12 nm diameter and 16 nm wide tunnel barriers. L is the barrier height created by LUMO bending due to surface gate potentials. The tunnel resistance is given by [37] Vg =0.20V 0.15 RT−1 = 0 −0.15 −100 0 Position(nm) 100 200 FIG. 5. (Color online) LUMO bending as function of position along the source-island-drain axis for varying gate voltages within the device target operating regime. The barrier function behaves linearly with Vg . Vertical marks show the boundaries of the labelled regions (black-dashed lines with arrow heads). ε is not shown since it varies with Vg , but its values for each Vg can be identified as the intersections of the corresponding LUMO data traces with the region boundaries. The transmission coefficient of a finite square barrier is: T = e−2γ , (A4) where 1 γ= h̄ (A6) where | T (ε) | is the transmission co-efficient, ρ1 is the density of states in source/drain electrodes and ρ2 is the density of states in the island. We approximate the density of states by assuming constant band bending over the leads (and island), giving: Barrier Barrier −200 2πe2 ρ1 ρ2 | T (ε) |2 , h̄ Za q q a 2m(Lpeak − ε)dx = 2m (Lpeak − ε), h̄ 0 (A5) h̄ = h/2π is the reduced Planck’s constant, m = 0.38mo is the effective mass of electrons in the MoS2 bilayer [36], a is the width of the tunnel barrier, and ε is the energy of an incoming electron in the MoS2 bilayer. Here ε can also be regarded as ε = Vsd + Eth + µ, where Vsd is the source-drain bias, Eth is the thermal energy described by the Fermi smearing incorporated below, and µ is the chemical potential of the electron which is set by the absolute source and drain voltages Vs and Vd . (Recall Vsd = Vs − Vd ; modifying the absolute value of each in tandem has no effect on Vsd .) We assume that the incoming electron has µ such that ε ≈ 0.38 [Lpeak − L (−200)] = 0.38Lpeak ; this defines the points where L = ε (the region boundaries) at consistent locations in space – enforcing the desired geometry with large charging energy. (The left edge of L has been set as zero energy, the absolute value of µ is somewhat different to the marked ε.) Fig. 5 shows L (x) for several values of Vg , together with the constant region boundaries. In principle, this assumption could be relaxed as needed to describe a particular physical device. As an example, Fig. 2(b) in the main text shows the particular case of the optimised structure under ρ= Zεf D(ε) dε = 2 D(ε − εo ) × (ε − εo ); 3 (A7) εo Here εf is the Fermi level and εo is the energy of the lowest level in the source/drain leads of MoS2 bilayer plane. (This corresponds to making the assumption that electron-hole recombination is negligibly slow compared to the tunneling rate). D(ε) is the density of states per unit energy in a solid, given by √ 8π 2 3 √ D(ε) = (A8) m 2 ε × η for ε > εo , h3 where η is the volume. Hence, √ 16π 2 3 √ m 2 ε − εo (ε − εo ) × ηs for ε > εo ρ1 = 3h3 (A9) and √ 16π 2 3 √ ρ2 = m 2 ε − εo (ε − εo ) × ηi for ε > εo (A10) 3h3 Here ηs and ηi are the volumes for the source and island respectively. Volumes were estimated assuming the bilayer height of 0.94 Å, and island and lead in-plane areas 2 of π (6 nm) and 400 × 20 nm2 , respectively which leads to tunnel resistance of 3500 MΩ to 8500 MΩ for this idealised design in perfect MoS2 in HfO2 . Practically, many factors involved in experiment (e.g., fabrication imprecision, defects, etc.) often accumulate to affect resistances by several orders of magnitude. Although we use a semiconductor material (MoS2 bilayer) in the design of our SET, after connecting it with the external batteries the Fermi level rises into the conduction band in the source/island/drain regions while sitting below the conduction band in the tunnel barrier regions. Thus, the SET behaves effectively like a metalinsulator-metal junction and all the standard approximations for metallic SET are still applicable to our SET. To determine the current through the SET, we need the tunneling rates on and off the island for different 7 charge configurations of the SET. These tunneling rates are given by [7] Γnχi = n ∆Eχi 1 , n qe2 RT exp( ∆Eχi ) − 1 kB T (A11) where χ denotes the source or drain and i the island. n ∆Eχi is the energy difference between different charge configurations of the SET and is given by Here Vχ is the chemical potential applied to the source/drain lead. The current through the SET is sum of the probabilities for all possible tunneling rates on and off the island. This current is given by I = qe ∞ X X n=−∞ χ=s,d pn (Γnχi − Γniχ ) (A13) (A12) where pn is the probability that island is in a state with n excess electrons. These probabilities are calculated by the master equation as described in [38]. [1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005). [2] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Phy. Rev. Lett. 105, 136805 (2010). [3] G. Fiori, F. Bonaccorso, G. Iannaccone, T. Palacios, D. Neumaier, A. Seabaugh, S. K. Banerjee, and L. Colombo, Nat. Nanotechnol. 9, 768 (2014). [4] K. Bolotin, K. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, and H. Stormer, Solid State Commun. 146, 351 (2008). [5] A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C. Y. Chim, G. Galli, and F. Wang, Nano Lett. 10, 1271 (2010). [6] V. Singh, A. Agrawal, and S. Singh, IJSCE 2, 502 (2012). [7] H. Grabert and M. H. Devoret, (Plenum Press, New York, 1992). [8] S. Mahapatra, V. Vaish, C. Wasshuber, K. Banerjee, and A. M. Ionescu, IEEE Trans. Electron Devices 51, 1772 (2004). [9] M. Devoret and R. Schoelkopf, Nature 406, 1039 (2000). [10] S. J. Kim, J. J. Lee, H. J. Kang, J. B. Choi, Y.-S. Yu, Y. Takahashi, and D. G. Hasko, Appl. Phys. Lett. 101, 183101 (2012). [11] Z. L. Yuan, Science 295, 102 (2002). [12] M. Fuechsle, S. Mahapatra, F. A. Zwanenburg, M. Friesen, M. A. Eriksson, and M. Y. Simmons, Nat. Nanotechnol. 5, 502 (2010). [13] M. Hofheinz, X. Jehl, M. Sanquer, G. Molas, M. Vinet, and S. Deleonibus, Appl. Phys. Lett. 89, 143504 (2006). [14] S. J. Shin, C. S. Jung, B. J. Park, T. K. Yoon, J. J. Lee, S. J. Kim, J. B. Choi, Y. Takahashi, and D. G. Hasko, App. Phys. Lett. 97, 103101 (2010). [15] R. G. Knobel and A. N. Cleland, Nature 424, 291 (2003). [16] P. Stokes and S. I. Khondaker, Appl. Phys. Lett. 92, 262107 (2008). [17] C. Thelander, M. H. Magnusson, K. Deppert, L. Samuelson, P. R. Poulsen, J. Nygard, and J. Borggreen, Appl. Phys. Lett. 79, 2106 (2001). [18] D. L. Klein, R. Roth, A. K. L. Lim, A. P. Alivisatos, and P. L. McEuen, Nature 389, 699 (1997). [19] T. W. Kim, D. C. Choo, J. H. Shim, and S. O. Kang, Appl. Phys. Lett. 80, 2168 (2002). [20] C. Stampfer, E. Schurtenberger, F. Molitor, J. Güttinger, T. Ihn, and K. Ensslin, Nano Lett. 8, 2378 (2008). [21] X.-X. Song, D. Liu, V. Mosallanejad, J. You, T.-Y. Han, D.-T. Chen, H.-O. Li, G. Cao, M. Xiao, G.-C. Guo, and G.-P. Guo, Nanoscale 7, 16867 (2015). S. J. Angus, A. J. Ferguson, A. S. Dzurak, and R. G. Clark, Nano Lett. 7, 2051 (2007). A. Ramasubramaniam, D. Naveh, and E. Towe, Phys. Rev. B 84, 205325 (2011). Q. Liu, L. Li, Y. Li, Z. Gao, Z. Chen, and J. Lu, J. Phys. Chem. C 116, 21556 (2012). comsol multiphysics R version 5.1, COMSOL, Inc., Burlington, MA, USA, 2008. R. Dovesi, R. Orlando, B. Civalleri, C. Roetti, V. R. Saunders, and C. M. Zicovich-Wilson, Zeitschrift für Kristallographie 220, 571 (2005). R. Dovesi, V. R. Saunders, C. Roetti, R. Orlando, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. J. Bush, M. Llunell, C. Science, and A. Technologies, Crystal09 User’s Manual (2013). R. G. Dickinson and L. Pauling, J. Am. Chem. Soc. 45, 1466 (1923). J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100, 136406 (2008). F. Corà, A. Patel, N. M. Harrison, C. Roetti, and C. R. A. Catlow, J. Mater. Chem. 7, 959 (1997). A. Ehlers, M. Böhme, S. Dapprich, A. Gobbi, A. Höllwarth, V. Jonas, K. Köhler, R. Stegmann, A. Veldkamp, and G. Frenking, Chem. Phys. Lett. 208, 111 (1993). M. F. Peintinger, D. V. Oliveira, and T. Bredow, J. Comput. Chem. 34, 451 (2013). H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976). W. G. van der Wiel, S. De Franceschi, J. M. Elzerman, T. Fujisawa, S. Tarucha, and L. P. Kouwenhoven, Rev. Mod. Phys. 75, 1 (2002). M. Jones, Y. Kwon, and D. Norton, Appl. Phys. A 81, 285 (2005). T. Cheiwchanchamnangij and W. R. L. Lambrecht, Phys. Rev. B 85, 205302 (2012). P. Koppinen, Bias and temperature dependence analysis of the tunneling current of normal metal-insulatornormal metal tunnel junctions, Master’s thesis, University of Jyvaskyla, department of physics (2003). V. I. Conrad, A. D. Greentree, D. N. Jamieson, and L. C. L. Hollenberg, J. Comput. Theor. Nanos. 2, 214 (2005). n ∆Eχi = E(n − 1) − E(n) + Vχ qe [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38]