Electronic Transport in 2D-Based Printed FETs from a Multiscale Perspective

www.advelectronicmat.de 2100972 (1 of 10) © 2022 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH ReseaRch aRticle Electronic Transport in 2D-Based Printed FETs from a Multiscale Perspective Marta Perucchini, Damiano Marian,* Enrique G. Marin, Teresa Cusati, Giuseppe Iannaccone, and Gianluca Fiori* DOI: 10.1002/aelm.202100972 dispersions can be deposited layer-by-layer on a wide range of substrates to define the different components in electronic circuits.[2,3] The recent experimental realization of inkjet-printed nanosheet-network transistors has already demonstrated the potential of this approach in the work of Kelly et al.,[4] where all device elements (electrodes, channel and insulator) are made of printed 2D materials flakes. Similarly, Carey et al.[5] realized a fully inkjet-printed field-effect transistor (FET) on a textile substrate, moving a step forward in the applications of such structures in wearable and flexible electronics, whereas Jewel et al.[6] and Monne et al.[7] improved the on/off current ratio of FETs by mixing graphene and MoS2 inks. By means of most common magnitudes such as resistance and mobility, electrical characterization of solution-processed devices is paramount for assessing and benchmarking their performance with respect to more traditional architectures based on single-flake channels. For instance, printed transistors based on transition metal dichalcogenides (TMDs) in ref. [4] showed a network-channel mobility around 0.1 cm Vs 21 1 2 MoS µ = −− , while Lin et al.[8] reported room temperature mobility of about 10 cm2 V−1 s−1 with MoS2 spin-coated thin-film transistors. Graphene conductive lines have instead been characterized in the latest work realized by some of the authors,[9] demonstrating μ in the order of tens of cm2 V−1 s−1. However, the main transport mechanisms at play in nonhomogeneous structures are still far from being understood, so that a rigorous theoretical framework able to provide physical insights as well as fabrication guidelines is needed and currently lacking. The printed structures are, indeed, made of 3D arrangements of 2D flakes, which form an interconnected network (i.e., a heterogeneous material), substantially different from thin films generally used in the fabrication of transistors and electronic circuits with conventional micro and nanoelectronic techniques. For this reason, the study of electronic transport in these devices requires a physical description able to capture the phenomena involved both at the nanometer scale of the flakes and at the micrometer scale of the networks. These resemble to carbon nanotube (CNT) or metal nanowire networks, for which some models, mainly based on percolation theory, have already been proposed in the last decade.[10–12] For the sake of clarity, with percolation one refers to the onset of conduction across a region previously showing an insulating As 2D materials (2DMs) gain the research limelight as a technological option for obtaining on-demand printed low-cost integrated circuits with reduced environmental impact, theoretical methods able to provide the necessary fabrication guidelines acquire fundamental importance. Here, a multiscale modeling technique is exploited to study electronic transport in devices consisting of a printed 2DM network of flakes. The approach implements a Monte Carlo scheme to generate the flake distribution. By means of ab initio density functional theory calculations together with non equilibrium Green’s functions formalism, detailed physical insights on flake-to-flake transport mechanisms are provided. This later feeds a 3D drift-diffusion and Poisson solution to compute self-consistently transport and electrostatics in the device. The method is applied to MoS2 and graphene-based dielectrically gated FETs, highlighting the impact of the structure density and variability on the mobility and sheet resistance. The prediction capability of the proposed approach is validated against electrical measurements of in-house printed graphene conductive lines as a function of film thickness, demonstrating its strong potential as a guide for future experimental activity in thefield. M. Perucchini, D. Marian, T. Cusati, G. Iannaccone, G. Fiori Dipartimento di Ingegneria dell’Informazione Università di Pisa Pisa 56122, Italy E-mail: [email protected]; [email protected] E. G. Marin Dpto. Electrónica, Fac. Ciencias Universidad de Granada Granada 18071, Spain The ORCID identification number(s) for the author(s) of this article can be found under https://doi.org/10.1002/aelm.202100972. 1. Introduction In recent years there has been a widespread effort to integrate 2D materials (2DMs) in the most diverse and innovative electronic applications, so as to exploit at best their unique properties like atomic thickness, transparency, and flexibility. One of the biggest opportunities is represented by printable electronics, as recently demonstrated in ref. [1]. Through this technology, inks prepared with various 2D-materials flakes © 2022 The Authors. Advanced Electronic Materials published by WileyVCH GmbH. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Adv. Electron. Mater. 2022, 2100972 www.advancedsciencenews.com www.advelectronicmat.de 2100972 (2 of 10) © 2022 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH behavior, when the number of conducting connections reaches a density that exceeds a critical value, the so-called percolation threshold. Above this value, the more links are added, the more the conductivity increases thanks to the multiplication of conducting paths, up to the point that the whole structure is interconnected in a continuous network, showing a bulk-like behavior.[13] The modeling problem was early faced for CNT-based 2D networks by Kumar et al., and Pimparkar et al.[10–12] with approximated tube-to-tube coupling described in a driftdiffusive regime, and followed by the works in refs. [14,15]. More recently, Colasanti etal.[16] proposed a more realistic 3D arrangement of the CNTs where the modeling of transport across single-point junctions was developed further, whereas the variability in the DOS of the carbon nanotubes was also investigated by Schießl etal.[17]. These models, however, cannot be directly applied to the case of arrangements of quasi 2D objects, like graphene or transition metal dichalcogenide flakes, since additional complexities are introduced, such as the need of understanding vertical transport across overlapping areas between different flakes.[18] Very recently, Rizzi etal.[19,20] have proposed a computationally efficient resistor network model for graphene-based conductor materials to calculate the overall electrical conductivity. Nonetheless, this and other methods require strong assumptions in the consideration of vertical transport: for example, the overlapping regions are collapsed into single points or external/fitting parameters are needed to model conductance between overlappingnanotubes. In this work, we leverage the most powerful available simulation approach, based on multi-scale modeling and able to deal with the different levels of abstraction needed to understand electronic transport in micrometer-sized devices made of 2DM flakes printed on top of dielectric substrates. We first define a randomly generated network, with arbitrarily oriented flakes and variable overlapping areas and connections through a Monte Carlo approach. Second, since the current flow in the structure is strongly dependent on the charge carrier transport along the in-plane and out-of-plane directions, we derive an anisotropic conductance, exploiting accurate ab initio density functional theory (DFT) calculations combined with nonequilibrium Green’s functions formalism.[21] Finally, we self-consistently solve the 3D Poisson and continuity equations in a device structure. Moreover, to take into account the intrinsic variability of the results due to the heterogeneous channel material, we provide a statistical analysis of the relevant quantities such as mobility, sheet resistance, and transfer characteristics. The model requires as input only the in-plane mobility, which can depend on the quality of the ink and on the experimental procedure after the printing process, and the density of the flakes, what we call the filling factor (FF), that is, the volume occupied by the flakes with respect to the overall channel volume. In this scheme we are able to compare our simulations with graphene-based printed conductive lines realized in-house, and evaluated as a function of the film thickness.[1,9] This is a first, but important step toward assessing the capabilities of the proposed approach to predict the electronic behavior of inkjet-printed devices, contributing to the efforts of narrowing the gap between theory and experiments in printed electronics, and that can possibly be extended to other materials and structures. By means of such approach, we can suggest possible routes to improve the performance of graphene and MoS2 dielectrically gated printedFETs. 2. Results and Discussion To exemplify the capabilities of the proposed multiscale approach as well as to provide guidelines regarding the main variables that determine electrical transport in 3D networks of 2D printed flakes, we have opted for a structure adapted from a 2D inkjet-printed FET.[1] In particular, we study the behavior of dielectrically gated FETs with the channel made of networks of flakes, in the three-terminal device architecture schematically represented in Figure1a,b. A gate voltage VG is applied to the top metal to electrostatically control the channel, while a drain-to-source voltage difference VDS is set between channel ends. The channel is composed by a randomly generated 2DMbased network (see Experimental Section) of flakes with a thickness of 0.5 nm and an average lateral dimension of 150 nm, the same order of magnitude of the flakes dispersed in the experimentally realized water-based inks,[1] randomly oriented on the x–y plane. A typical section of the flake distribution is shown in Figure1c for two overlapping planes parallel to the gate. In particular, we have selected two different materials, namely graphene and MoS2, although the proposed approach can be extended to any arbitrary 2D material ink or even combinations of them. The channel region is embedded in 10 nm of top and bottom insulator with dielectric constant εr= 3.0, the typical value for hBN.[22,23] The channel extends for 3 μm in length and 1 μm in width with, otherwise stated, a total thickness of 100 nm (see Section S8, Supporting Information for a summary of the quantity used in the simulations). Being the channel material heterogeneous, the number of flakes as well as their connectivity play a major role in determining the device behavior. Quantitative estimation of these two magnitudes is given by the FF. To gain a visual insight of the phenomena we have plotted the current isosurfaces for two different filling factors FF = 0.5 and FF = 0.7 in Figure1d under two different gate biases VGS= 1.0V and VGS= 2.0V. As can be seen, the current flow along the device channel is not only dependent on the applied gate voltage but also hindered (facilitated) by the absence (presence) of 2DM flakes. In the following, the impact of the FF on the current will be discussed in moredetail. In a 2DM network of flakes the in-plane and out-of-plane conductance largely differ, with fundamental consequences in the transport properties of a printed device. For a correct modeling of the current flow in these structures, it is fundamental to determine through accurate calculations the relation between in-plane and out-of-plane transport in the network, and its dependence on the overlapping area between flakes or on the type of stacking. While in-plane transport is characterized by a semi-classical diffusive mobility, the out-of-plane transport is mediated by tunneling through a narrow (few angstroms) spatial gap between flakes, whose coupling is governed by van der Waals forces. In order to connect, the diffusive/semi-classical/ in-plane and ballistic/ab initio/out-of-plane transport, we first calculate both in-plane and out-of-plane ballistic transmission Adv. Electron. Mater. 2022, 2100972 www.advancedsciencenews.com www.advelectronicmat.de 2100972 (3 of 10) © 2022 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH and obtain the conductance ratio between them. This conductance ratio can be interpreted as the probability for a charge carrier to keep on moving within the flake or to tunnel into another flake. In this way the mobility controls the overall transport, according to the solution of the drift–diffusion equations, while the out-of-plane conductance is treatedaccurately. In order to determine the out-of-plane versus in-plane conductance in overlapping graphene and MoS2 flakes, we have performed ab initio DFT calculations by means of the Quantum Espresso suite, that we have combined with Maximally localized Wannier functions, nonequilibrium Green’s functions and Landauer–Buttiker theory [24,25] (details can be found in Experimental Section). In particular, we have first defined, for both materials, two crystal structures to compute in-plane (Figure2a,b top) and out-of-plane transmissions (Figure 2a,b bottom), mimicking the situation of electrical transport along a single flake or, vertically, across two flakes. The influence of the overlapping area between the flakes as well as their stacking arrangement (AA and AB for graphene, and AA′ and AB for MoS2) have been considered in the study. From the charge density difference in Figure2c,d obtained as described in the Supporting Information, it can be seen that a redistribution of charge happens in the overlapped region. For graphene, as expected, the absolute value of the calculated in-plane transmission (visible in Figure2e) increases linearly with energy. The out-of-plane transmission also increases with energy, although non-linearly, saturating with the overlapping area, while relevant differences (close to a factor × 2) are observed between AA and AB arrangements. MoS2 instead, with a 1.78eV gap, shows a region of null transmission in all curves as depicted in Figure2f: noticeably the relationship between transmission and energy is nonlinear and it is asymmetric for positive and negative energy values. The differences between AA′ and AB stacking arrangements are also asymmetric for electrons and holes, being pronounced in the valence band and small in the conduction band. The overlapping area has in both cases a small impact. These differences in flake coupling, at the nanoscopic scale, are expected to be blurred up by the randomness of the printed structure, when moving to the device microscopic scale. As a consequence, we have taken an average transmission coefficient and used it to characterize the tunnel transport between flakes. Next, applying the Landauer–Buttiker linear response theory,[26,27] we can obtain the conductance G versus Fermi level position, and, to assess transport anisotropy, we calculate the ratio of the out-of-plane to the in-plane conductance (G1/G0), together with their fitting curves shown Figure2g,h (more details can be found in the Supporting Information). For MoS2 we have only considered energies above the conduction band level, which correspond to n-type conduction, the one typically obtained by unintentional doping, and therefore analyzed in the presentwork. Figure 1. a) Sketch of a 2DM network based device and b) schematic of a top-gated FET structure considered for the calculations; c) a typical simulated flake distribution showing two overlapped sections along the transistor vertical direction in which the flakes with average lateral dimension L= 150nm occupy 50% of the volume; d) Isosurfaces of the current density absolute value as a function of the applied gate voltage (VG= 1.00V on the left, VG= 2.00V on the right) for filling factor equal to FF = 0.5 (top) and FF = 0.7 (bottom). Adv. Electron. Mater. 2022, 2100972 www.advancedsciencenews.com www.advelectronicmat.de 2100972 (4 of 10) © 2022 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH In order to validate the prediction capability of the proposed model and to provide insights into the mechanisms at play in fabricated devices, we refer to the graphene-based structures realized in refs. [1] and [9]. These samples are inkjet-printed onto SiO2/Si substrates controlling the material deposition by changing the number of print passes and therefore resulting in different channel thickness. Figure3a shows the experimental measurements of the sheet resistance versus thickness before (diamonds) and after (squares) annealing, which are compared with the performed simulations (circles and triangles for before and after annealing cases, respectively). For these structures it has been observed that environmental n-type doping is present, to the amount of around 1012 carriers cm−2 for graphene layer as measured in ref. [9]. The origin of this unintentional doping is attributed to residuals from ink solvents and chemical components, whereas humidity and water molecules would otherwise Figure 2. Schematics crystal structures of a) graphene and b) MoS2 for the calculation of in-plane (top figure) and out-of-plane (bottom figure) transmission, respectively. Charge density difference calculated for a section perpendicular on the y-axis of the overlapping structures: c) for graphene and d) for MoS2; transmission coefficient for e) graphene and f) MoS2 calculated for the monolayer and the partially overlapped structures and for different stacking and overlapping areas. g,h) Average conductance ratio (out-of-plane/in-plane) and fitting curves for the two materials. Adv. Electron. Mater. 2022, 2100972 www.advancedsciencenews.com www.advelectronicmat.de 2100972 (5 of 10) © 2022 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH lead to p-type doping. Coherently, we have assumed an unintentional doping density in graphene of Nd= 1012 cm−2 for the calculations reproducing the results in ref. [9] and Nd= 1011 cm−2 for those obtained in ref. [1]. The in-plane mobility of graphene in both cases is set to 200cm2 V−1 s−1.[28,29] Two different filling factors FF = 0.35/0.3 (ref. [1] or [9]) and FF = 0.7 are considered for the samples before and after annealing respectively, since the network of flakes is more compacted after the annealing and thus its density is increased. For each thickness, a set of 30 devices with different actual flake distribution is considered. As can be seen, the simulations are in very good agreement with the experimental results. We have sought to perform the same extensive validation for the extracted mobility in printed devices, but mobility measurements in printed lines and in the range of thicknesses considered in this study are hardly available in the literature. In ref. [9] some of the authors measured mobility in the range of 10–15cm2 V−1 s−1 (for thicknesses 125 ± 45nm and 235 ± 70nm) similar to the ones we obtain from our calculations (see Supporting Information). Following the same procedure, we calculated the sheet resistance also for MoS2-based network, for which experimental measurements as a function of thickness are however not available (Figure3b). Given the early stage of the printed technology and the complexity of the nonhomogeneous materials studied here, it is important to understand how the channel structure (described through the FF) impacts physical variables such as mobility and sheet resistance, that are key figures of merit for the device. It has to be mentioned that the FF in the simulations is linked, from a fabrication point of view, to process parameters like the number of print passes, the ink concentration or viscosity, and the drop spacing although a quantitative relationship cannot be easilyestablished. So that, we keep the device structure and the physical variables previously validated and we study the printed FET mobility and sheet resistance for five sets of 30 devices with a fixed thickness of 100nm as a function of the network density (Figure4). In particular, for the extracted mobility (Figure 4a) we observe values from ≈1 to ≈16cm2 V−1 s−1 increasing almost linearly with FF. The best value of the estimated device mobility is, therefore, around one-tenth of the in-plane mobility (μGr= 200cm2 V−1 s−1) that has been considered, pointing out that the connectivity of the network and the out-of-plane/in-plane conductance impact highly on the actual effective mobility of the device. Improving both the in-plane mobility and also assuring a good interaction between flakes within the inks are the objectives to pursue in the search for better experimentalperformance. By fitting the points with linear regression, it has been possible to extrapolate the x-axis intercept at FF ≈ 0.3, which can be considered as a threshold below which the effective mobility goes to zero since almost no conductive paths exist. A similar condition was also found in ref. [9] for the printed samples with a number of print passes around 10. Following the same procedure, we calculated the mobility also for MoS2-based network structures, for which experimental data in literature report a value of the order fo 0.1 cm2 V−1 s−1.[4] The simulated (undoped) devices show a linear relationship between FF and the mobility (see Figure4b) with values ranging from almost 0 to ≈1.0~cm2 V−1 s−1 (considering 2 2 µ = MoS cm2 V−1 s−1 as in refs. [4,30]), however highlighting a similar percolating threshold as for graphene. Indeed, the extremely low mobility of MoS2 inks is the main limiting factor jeopardizing the experimental successful realization of MoS2 inkjet-printed devices, with some remarkable exceptions as in the work by Lin etal.,[8] where MoS2 mobility values of the order of 10 cm2 V−1 s−1 are demonstrated. From our calculations, we observe that in the best scenario (FF=0.7) the obtained device mobility does not suffer a higher degradation with respect to the in-plane mobility (μ≈ 1.0 cm2 V−1 s−1 vs 2 2 µ = MoScm2 V−1 s−1) than in graphene (μ≈ 16 cm2 V−1 s−1 vs μGr= 200cm2 V−1 s−1). These results emphasize FF as a critical parameter for MoS2 printing, inferring that particular attention has to be paid not only to the good quality of the inks, but also to the printing parameters as the drop spacing, the number of print passes, and the ink concentration. In the inset of Figure4b, the cumulative occupation of the network is shown for the various FF. The experimental activity might seek higher FF in order to improve performance. Techniques as spin-coating 2DMs[8] Figure 3. a) Extracted sheet resistance for graphene as a function of the channel thickness compared with experimental results obtained in refs. [1] and [9] before (diamonds) and after (squares) annealing. Two different filling factors FF = 0.35/0.3 for ref. [1] or [9] and FF = 0.7 are considered in the simulations to reproduce the preand after annealing scenarios. b) Extracted sheet resistance MoS2 as a function of the channel thickness for two different filling factors FF = 0.4 and FF = 0.7. Adv. Electron. Mater. 2022, 2100972 www.advancedsciencenews.com www.advelectronicmat.de 2100972 (6 of 10) © 2022 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH might help to this purpose as well as to improve the interaction between flakes in the out-of-planetransport. As for the resistivity, we have calculated Rsheet for each device in a region for which the slope of the IDS–VDS curve is constant. The results are shown in Figure4c,d together with a fitting curve as a guide for the eye (details can be found in the Supporting Information). Rsheet is in the kΩ per square range for graphene and up to 108~Ω per square for MoS2-based FETs, increasing for lower FF. The different order of magnitude for the two materials can be expected due to the semiconducting electronic nature opposite to a semi-metallic one and their different anisotropic conductance. In both cases, however, the trend follows a percolating-type performance: if we consider approximately FF = 0.3 as the threshold, it is possible to see an abrupt decrease in resistance from FF = 0.3–0.4, which continues (even though to a lesser extent) for even larger FF values. From a qualitative perspective, we remark that the percolative behavior is a common feature of both experimental and simulated devices. As discussed in the mobility curves, the percolation threshold in the sheet resistance is found around a FF = 0.3 in our simulated structures. The predicted decaying trends for both materials highlight the impact of the FF on the device sheet resistance and draw attention to the necessity of having high flakes density in the inks so to avoid experimental realizations close to the percolating threshold. Above this point, the sheet resistance decreases approaching the bulk values, and the improvements in the sheet resistance figures of merit are moreincremental. In order to complete the view about the impact of the network structure and channel thickness on the electrical response, we have computed the transfer characteristics as a function of FF ranging from 0.3 to 0.7 and thicknesses from 50 up to 150nm for 30 different network distributions. The results for graphene as a function of FF and thickness are shown in Figure5a, where the solid lines indicate the mean current of the distribution Figure 4. Effective mobility and sheet resistance for a,c) graphene and b,d) MoS2 printed FETs, respectively as a function of the filling factor. The insets in (b) show the cumulative density extracted from a sample structure for each filling factor. Adv. Electron. Mater. 2022, 2100972 www.advancedsciencenews.com www.advelectronicmat.de 2100972 (7 of 10) © 2022 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH and the shaded area takes into account its standard deviation. The semi-metallic nature of graphene results in a conductive behavior along the whole channel thickness, regardless of the applied VGS. For a fairer comparison we have normalized IDS in Figure5a not only by the device width, but also by the channel thickness. While the IDS per unit width scales with the film thickness (see Figure S3, Supporting Information), the differences in IDS per unit area between thicknesses (Figure5a) become relevant only for larger FFs and for higher VGS biases, as the gate-induced charge fills with more carriers for transport the layers closer to the gate. On the other hand, although the modulation is very poor, even lower than what is expected in a graphene monolayer FET,[31] a well-defined trend correlates the increase in the current with FF. This can be explained by the presence of a large number of paths along the plane with preferential transport due to the higher conductance values. It is interesting to see that the IDS–VGS characteristics seem more affected by the connectivity among flakes, thus by their density (FF), rather than by the applied VDS considered. In fact, by comparing Figure5a where the FET transfer characteristics are shown as a function of FF for the same VDS= 0.3V and Figure 5b where, on the contrary, the curves are a function of the drain-to-source voltage for FF = 0.5 (and a thickness of 100nm), it can be seen that it is not possible to reach currents larger than 160 μA μm−2 modulating VDS up to =0.5V, whereas the value can be overtaken increasing FF to 0.7. With the same device structure, we have performed transport calculations for n-type FETs based on MoS2 networks. The Figure 5. Transfer characteristics for 30-samples simulations as a function of a) the flakes filling factor FF (different colors) and of channel thickness (empty, colored or black symbols) and as function of b) the drain-to-source voltage for graphene FETs. Transfer characteristics as a function of the flakes filling factor FF (different colors) and of channel thickness (empty, colored or black symbols) in c) linear and d) semilogaritmic scale for MoS2 FETs. The solid line indicates the mean value and the shaded areas include the standard deviation. Adv. Electron. Mater. 2022, 2100972 www.advancedsciencenews.com www.advelectronicmat.de 2100972 (8 of 10) © 2022 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH transfer characteristic curves in Figure 5c,d are obtained for VDS= 0.2 V and again FF ranging from 0.3 to 0.7 and thicknesses from 50 up to 150nm, each case consisting of a set of 30 different random distribution of flakes. For MoS2, IDS negligibly depends on the film thickness (from 50 to 150 nm). Indeed, the carriers contributing to IDS in MoS2 printed devices are mainly induced close to the gate, and more importantly the degradation of the transport through multiple tunnel transmissions constraints the possibility of reaching higher currents. A good modulation of the current with the positive applied gate bias can be appreciated, and in particular, the sub-threshold swing values reach, as expected, 60 mV decade−1 for all network density, except FF = 0.3 for which is equal to 73 mV decade−1. As can be seen, IDS is very sensitive to FF, with an increase of more than two orders of magnitude in the range of densities considered, a factor which may be of paramount importance when such devices are experimentally fabricated. The On/Off ratio has also been extracted. Considering the small applied voltage, still, it has been possible to achieve On-to-Off current ratios for FF=0.7 of more than two orders of magnitude, with an average of ≈558. This On-to-Off current ratio can be expected to reach larger values, in line with IRDS specifications,[32] when higher drain-to-source voltage isconsidered. The performance of transistors in terms of On/Off current ratio however degrades to 51, when reducing the flake density to 0.3. It has to be noticed however that the ratio reported for experimentally realized MoS2 on SiO2 printed transistors is in the range of 3 to 6,[33] a low value which can be explained mainly by the presence of defects and traps which deteriorate the electronic properties of the material,[33] and that here have beenneglected. Examining the fabricated devices and comparing their parameters with those extracted from the ones we simulated by means of the modeling approach proposed, it is possible to highlight some fundamental aspects of transport in network-based channels. On one hand we observed how the out-of-plane tunnel transmission is detrimental for the overall network-mobility. In this sense, attention should be drawn in optimizing the interface between overlapping flakes, for instance by means of solvents or chemical additives in the ink formulation or by modulating the flakes thickness and phase as in ref. [8]. On the other hand, by increasing the amount of flakes in the channel it is possible to increase the network mobility while simultaneously reducing the sheet resistance. From our understanding, this trend can be attributed to the multiplication of connections between flakes and consequently of the number of percolative paths. From a fabrication perspective, the effect can be augmented by means of strategies designed to improve the network compactness, such as—but not limited to— annealing processing steps. The impact of the channel thickness is very dependent on the channel material: while the semi-metallic nature of graphene results in a proportional scaling of the current with the channel thickness, in MoS2 the current barely changes with it, as most of the conductive paths are determined by the region of the channel close to the gateterminal. 3. Conclusion In this work, we have presented a simulation platform to study the electrical properties and behavior of FETs based on 2DM flakes printed network, intended as a support for the inkjet printing process, a development that constitutes an advancement to be used in the realization of flexible electronic components and for the interpretation of the device operation. To this purpose, we have exploited the multi-scale modeling approach combining different levels of abstraction in order to comprehend the electronic transport in micrometer-sized devices made of nanometric 2DM flakes. We believe it is an important step toward the prediction of the electronic behavior of inkjetprinted devices, helping in narrowing the gap between theory and experiments in printed electronics. The model involves the random definition of the network structure, the evaluation of the out-of-plane versus in-plane conductance from ab initio DFT and nonequilibrium Green’s functions simulations and the electronic transport calculations through a self-consistent solution of Poisson and continuity equations under a drift–diffusion framework. The results for the transfer and output characteristics for different sets of structures show a scenario typical of percolation, where the structure density plays a fundamental role in the resulting device properties together with the electronic structure of the material. In particular, graphene flake networks show very poor modulation of the current with the applied gate voltage, whereas MoS2 flake networks show ideal SS values when defects and traps are neglected. By extracting the effective mobility for both graphene and MoS2 n-FET devices from the computed IDS–VG curves it has been possible to define a percolating threshold (FF ≈ 0.3). Given as input the in-plane mobility of the flakes and the filling factor, parameters that can depend on the particular experiment taken into account, the model is able to capture the behavior of available experimental results for graphene inkjet-printed structures for sheet resistance as a function of the film thickness. Moreover, the present model indicates the out-of-plane conductance—which embodies the flake-to-flake resistivity—and the filling factor, which is linked to the network connectivity, as key parameters for the design of new devices, being the film thickness relevant in graphene channels and negligible in MoS2 ones. Finally, the good predictability of the methodology used in this work can be applied to other 2DMs networks and even to mixtures of such materials, providing insights into the behavior of nonhomogeneous channels. Moreover, the presented platform can be exploited to study the impact of different geometrical variables, and their variability, into the device electrical performance, including, for example, 2D flakes orientation with respect to the transport direction, relative angle orientation between flakes, flake size, flakes tilt, interface traps, grain boundaries, etc. These parameters can eventually have a relevant impact as it happens, for example, for transport across grain boundaries in graphene.[34–37] Finally, we aim at extending our study exploring different device architectures andgeometries. 4. Experimental Section Transmission Calculations: The model proposed for investigating the transport between flakes started from a DFT-level simulation, (using the Quantum Espresso suite[24]) of an bilayer system comprising the two flakes of interest. The structure was made of two graphene layers, overlapped in the AB or AA stacking fashion with an interplanar distance Adv. Electron. Mater. 2022, 2100972 www.advancedsciencenews.com www.advelectronicmat.de 2100972 (9 of 10) © 2022 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH of 3.3 Å. This value was obtained from the equilibrium configuration spacing in ref. [38]. Periodic boundary conditions were assumed in the x and y direction and to avoid interactions due to the repetition of the cell in x direction, 20Å of vacuum was considered perpendicularly to the 2D material layers. Similarly, for MoS2 two layers overlapped in a AB or AA′ fashion, with the experimental interplanar distance of 2.98Å.[39] Monkhorst–Pack 1 × 12 × 12 and 1 × 6 × 6 k-meshes has been used for the Brillouin-zone integration in graphene and MoS2,respectively. From the DFT band-structure we extract a Hamiltonian in a tightbinding-like form in a basis of maximally localized wannier functions (MLWF) through the Wannier90 software,[25] that later allowed a fine manipulation to build systems of arbitrary size. Inspired by the procedure employed in refs. [40,41], a Wannier tight-binding-like Hamiltonian is achieved (through carefully projecting the band-structure on a precise selection of atomic orbitals), where one can clearly identify two diagonal blocks in correspondence to the top and bottom flakes, and two off-diagonal blocks that explain the coupling between the layers (see Figure S1a, Supporting Information and ref. [18]). With this procedure one can recover from each of the diagonal blocks the bands of the single flake, while the full Hamiltonian still contained all the information about the tunneling coupling coming directly from first principles. One can easily: i) study the heterostructure of interest by extending the Hamiltonian with single-layer blocks at both sides, so to build a single-flake/overlapped region/single flake geometry; or ii) define the size of the overlapping area by replicating the full Hamiltonian (see Figure S1b, Supporting Information and ref. [18]). Next, the Wannier Hamiltonian and nonequilibrium Green’s functions were employed to obtain the out-of-plane tunneling transmission for practical printed systems, that is, scaling the flakes to hundreds of nanometers using NanoTCAD ViDES.[21,42] Transmission was calculated considering 360 k−points in the ydirection. Monte Carlo Network Generation: In order to provide an accurate description of the materials structure for printed devices, an algorithm to create a random arrangement of flakes in a 3D grid was defined, extending the framework of the NanoTCAD ViDES software.[21] With this method it was possible to individually assign the properties (such as orientation and dimension) to the flake objects to have full customization. Once the shape and thickness of the flakes were set, their size and orientation with respect to the horizontal axis can be randomly generated following a normal distribution with the desired mean and variance. In this manuscript the flakes were oriented parallel to the main transport axis connecting source and drain. Similarly, a set of (x, y) coordinate was randomly extracted for each flake (on a pool covering the device area on the x–y plane) to represent its origin and central point, from which the desired shape was constructed. Afterward, each plane in the third direction (z) was simultaneously filled with the proper number of flakes to reach the desired density, that is, the network FF. The flake planes were thus assumed to be parallel to each other. The flakes in each single layer were instead sequentially positioned: in order to avoid nonphysical situations those flakes falling on already occupied grid points were, in fact, to beexcluded. Electrostatics and Transport Self-Consistent Solution: The algorithm coupled the 3D electrostatics solution and the 3D continuity equation under a drift-diffusion approach. As initial guess, the Poisson’s equation was self-consistently solved by removing the drain-to-source voltage and fixing the proper Dirichlet conditions for the desired VG with an input quasi-Fermi potential equal to 0V. To evaluate the charge in the MoS2 flakes the analytical carrier density expression based on the parabolic effective mass approximation was used although any arbitrary analytical density of states distribution can be considered instead (as it had been done for graphene). The grid resolution was 0.5nm along the z direction, that is, treating the flakes as purely bidimensional, and 5nm along the x and y directions. The channel region was embedded in 10 nm of top and bottom insulator with dielectric constant εr = 3.0, the typical value for hBN.[22,23] The same homogeneous permittivity was assumed within the channel region, since the impact of a heterogeneous ε was observed to be negligible. The charge and the quasi-Fermi level were obtained via the conjugate gradient iterative method. The relationship between carrier concentration and in-flake and inter-flake transport, was solved self-consistently considering properly both of the potential and the quasi-Fermi level spatial distributions in the device, and the transport anisotropy was considered in the different in-plane and outof-plane mobility. The error in potential between two iterations was checked against the convergence criteria (5 × 10−3), until convergence was reached. In such case, the current density components can be calculated for each direction as well as the total output current collected at the drain contact. For what concerns the mobility values, they were defined independently for the six faces of each grid elemental volume so to have μ= 0 cm2 V−1 s−1 for those channel points facing an “empty” regions, otherwise equal to μx, μy for neighboring points in the horizontal planes and μz(EF) as discussed in the main text for vertically adjacent channelpoints. Supporting Information Supporting Information is available from the Wiley Online Library or from the author. Acknowledgements The authors gratefully acknowledge the ERC PEP2D (contract no. 770047) and H2020 WASP (contract no. 825213) and the Crosslab Department of Excellence project for financial support. E.G.M. also acknowledges the support by the Spanish MCI through the project PID2020-116518GB-I00 and by the Junta de Andalucía-Consejería de Economía y Conocimiento/ FEDER-EU through the Project A-TIC-646-UGR20. Open access funding provided by Universita degli Studi di Pisa within the CRUI-CARE Agreement. Conflict of Interest The authors declare no conflict ofinterest. Data Availability Statement The data that support the findings of this study are available from the corresponding author upon reasonablerequest. Keywords 2D materials, drift–diffusion, field-effect transistor, multi-scale simulations, network, printable electronics Received: September 9, 2021 Revised: November 19, 2021 Published online: [1] D. McManus, S. Vranic, F. Withers, V. Sanchez-Romaguera, M. Macucci, H. Yang, R. Sorrentino, K. Parvez, S. K. Son, G. Iannaccone, K. Kostarelos, G. Fiori, C. Casiraghi, Nat. Nanotechnol. 2017, 12, 343. [2] J.Wang, Y.Liu, Z.Fan, W.Wang, B.Wang, Z.Guo, Adv. Compos. Hybrid Mater. 2019, 2, 1. [3] G. L.Goh, H.Zhang, T. H.Chong, W. Y.Yeong, Adv. Electron. Mater. 2021, 7, 2100445. Adv. Electron. Mater. 2022, 2100972