Modeling Hydrodynamic Charge Transport in Graphene

Graphene has exceptional electronic properties, such as zero band gap, massless carriers, and high mobility. These exotic carrier properties enable the design and development of unique graphene devices. However, traditional semiconductor solvers based on drift-diffusion equations are not capable of modeling and simulating the charge distribution and transport in graphene, accurately, to its full extent. The effects of charge inertia, viscosity, collective charge movement, contact doping, etc., cannot be accounted for by the conventional Poisson-drift-diffusion models, due to the underlying assumptions and simplifications. Therefore, this article proposes two mathematical models to analyze and simulate graphene-based devices. The first model is based on a modified nonlinear Poisson’s equation, which solves for the Fermi level and charge distribution electrostatically on graphene, by considering gating and contact doping. The second proposed solver focuses on the transport of the carriers by solving a hydrodynamic model. Furthermore, this model is applied to a Tesla-valve structure, where the viscosity and collective motion of the carriers play an important role, giving rise to rectification. These two models allow us to model unique electronic properties of graphene that could be paramount for the design of future graphene devices.


Introduction
Graphene has been the focus of various research topics for the last two decades and is pushing toward the market [1,2], and even new physical effects such as Moiré materials with superconductivity [3,4] and chiral plasmons [5] are still being discovered. For the future of graphene and for the development of graphene-based devices, it is crucial to be able to understand its many effects and to model the electronic transport in graphene. Thanks to its 2-D structure, graphene has exceptional electronic transport characteristics such as zero-mass carriers at the Dirac point and very high charge mobility. These properties make it challenging to use traditional semiconductor-modeling techniques with graphene. Additionally, due to the collective behavior of charge carriers and the dominant electronelectron (e-e) collisions occurring in graphene, considerably high viscosity levels in carriers of graphene have been observed. The viscosity of electrons in graphene has been shown to lead to phenomena such as whirlpools and negative local resistance [6], superballistic transport [7], and Poiseuille flow [8,9]. When wisely utilized, these hydrodynamic effects could also be utilized to implement electronic devices that show rectifying properties. A 2-D Tesla valve has been implemented [10] that acts like a viscometer and demonstrates rectification. Graphene is considered as a semi-metal or a zero-gap semiconductor; however, its special attributes and zero band gap make it unique. Thus, the use of traditional semiconductor models, based on classical drift-diffusion equations, is challenging. There have been multiple attempts to extend the drift-diffusion solvers to model graphene-based devices [11][12][13][14][15], and they are valid for a range of applications; however, they lack the ability to fully capture the aforementioned phenomena, due to viscosity and the many-body interactions occurring with the carriers in graphene. Moreover, the existing hydrodynamic models for charge transport in graphene are limited to the relativistic case at 0 K [16]. Therefore, a generalized hydrodynamic study is required to model the motion of electrons in graphene and graphene-based electronic devices accurately. By considering mass, momentum, and energy conservation of the charge carriers in graphene, the charge-transport dynamics can be predicted [17][18][19]. In this work, we aim to demonstrate novel computational solvers and algorithms to model graphene's electronic-transport properties, and we report two separate solvers based on the Finite Element Method (FEM). The first solver is to model the gating and contact doping in graphene, which solves the arising nonlinear Poisson's equation iteratively to determine the Fermi level and charge densities in graphene layer. The second reported solver models the hydrodynamic charge transport in the time domain, based on the Discontinuous Galerkin Time Domain (DGTD-FEM) algorithm, and delivers dynamic charge transport in graphene in the time domain.
This article is structured as follows. Section 2 focuses on the theory and physical model to cover the nonlinear electrostatic equation to model the electrostatic gating and contact doping in graphene as well as the hydrodynamic transport. Section 3 shows the results arising from the nonlinear electrostatic solver, as well as the transient hydrodynamic transport solver that demonstrates Poiseuille flow and rectification in a graphene-based Tesla valve. Section 4 discusses the results, possible future work, and concludes the paper.

Nonlinear Electrostatic Approach: Modeling Gating and Contact Doping
The transport properties of graphene depend on many internal and external parameters, such as its deposition technique, impurities, type of materials in contact with graphene, operating temperature, etc. These aspects need to be accurately considered for the design of next-generation graphene-based electronic devices in order to account for the effects, such as nonlinear contact doping, due to metals. To exploit the full potential of graphene's exotic properties, one should be able to access and manipulate the Fermi level of graphene. Making use of graphene in electrical devices requires shifting the Fermi level, by gating and making contacts with metals, which inherently dopes the graphene locally and results in Fermi level shifts. Both of these effects, gating and contact doping, can already be analyzed under static conditions where there is no net current flowing through the graphene.
The key relations that determine the electrostatic charge and Fermi level distribution are the Thomas-Fermi equation and Poisson's equation. These equations can be written for the analysis of charge distribution of graphene under electrostatic conditions. Using the equilibrium Fermi-Dirac function f FD (E − E F ), and 2-D density of states in graphene g 2D (E) = 2|E| π 2 v 2 F , the electron density per area n 2D e at a given Fermi Level E F for nonzero temperature T, is written as where k B is the Boltzmann constant and Li n (x) is the poly-logarithm function of the n-th order [12]. Electron-hole symmetry in graphene's band diagram yields a similar relationship for the hole density n 2D h (E F ) = n 2D e (−E F ). For devices where the graphene layer is sandwiched between dielectric layers, as in Figure 1, the following nonlinear Poisson's equation can be written [11,14,15]: where ϕ is the electrostatic potential, ε is permittivity, q denotes elementary charge, and t gr is the graphene's thickness, taken as 0.35 nm. Dividing by graphene's thickness on the right-hand side of (3) makes sure that the volumetric charge density is used for the charge density term in Poisson's equation. Here, Equation (2) represents the nonlinear Poisson's equation (NLP), as the charge distribution in graphene is also a function of the potential. The obtained potential on the graphene layer shifts the band diagram of graphene based on E = −qϕ, and the corresponding Fermi level with respect to the Dirac point can be updated. Therefore, solving this equation would require an iterative approach between (1) and (2), as in [20]. Gate potentials are applied by enforcing Dirichlet boundary conditions onto the model. For the other boundaries, a homogeneous Neumann boundary condition (∇ϕ·n = 0, where n is the outward-pointing unit normal vector), is used. The obtained solution directly gives the charge distribution in graphene under the gate potential (without contacts present directly on the graphene), since the scenario describes the electrostatic case without any current flowing. where is the electrostatic potential, is permittivity, denotes elementary charge, and is the graphene's thickness, taken as 0.35 nm. Dividing by graphene's thickness on the right-hand side of (3) makes sure that the volumetric charge density is used for the charge density term in Poisson's equation. Here, Equation (2) represents the nonlinear Poisson's equation (NLP), as the charge distribution in graphene is also a function of the potential. The obtained potential on the graphene layer shifts the band diagram of graphene based on = − , and the corresponding Fermi level with respect to the Dirac point can be updated. Therefore, solving this equation would require an iterative approach between (1) and (2), as in [20]. Gate potentials are applied by enforcing Dirichlet boundary conditions onto the model. For the other boundaries, a homogeneous Neumann boundary condition ( • = 0, where is the outward-pointing unit normal vector), is used. The obtained solution directly gives the charge distribution in graphene under the gate potential (without contacts present directly on the graphene), since the scenario describes the electrostatic case without any current flowing. In this work, Equation (2) has been discretized using the Finite Element Method (FEM) with a triangular mesh and with scalar (nodal) shape functions and is solved iteratively together with Equation (1), by using successive under-relaxation method as in [20,23]. Equation (2) is strongly coupled with Equation (1), where the potential is one coupling channel, and the concentration of carriers is the other one. In each iteration, for a given potential , carrier densities are computed by (1), and then they are transferred into  [21,22] (discussed further in the supplementary notes and Figure S1 In this work, Equation (2) has been discretized using the Finite Element Method (FEM) with a triangular mesh and with scalar (nodal) shape functions and is solved iteratively together with Equation (1), by using successive under-relaxation method as in [20,23]. Equation (2) is strongly coupled with Equation (1), where the potential is one coupling channel, and the concentration of carriers is the other one. In each iteration, for a given potential ϕ, carrier densities are computed by (1), and then they are transferred into the Poisson's Equation (2) to update the potential distribution again. This iterative scheme continues until the potential distribution stabilizes. The model in Figure 1a also assumes that graphene layer is in contact with a source of carrier, which does not affect the Fermi level position and does not cause doping of the graphene.
In addition to the analysis of gating on graphene, the contact doping and Fermi level pinning due to the contact of metals on graphene [24] can also be incorporated into this model, by applying additional Dirichlet boundary conditions as in Figure 1c. Since the potential ϕ is directly related to the Fermi level energy of graphene (E F ), it is possible to enforce the Fermi level pinning directly, by applying additional Dirichlet boundary conditions for the potential under equilibrium, and this leads to carrier doping due to the metal contacts. With these results, it is possible to visualize built-in fields, calculate the conductance of a device, and, ultimately, design contacts and dielectrics to find the ideal operation ranges for the gate voltages.

Hydrodynamic Charge Transport in Graphene 2.2.1. Hydrodynamic Transport Model
The accurate analysis of semiconductor devices requires advanced models for the charge carrier transport. The traditional drift-diffusion models (DDM) start losing their validity as the operation frequencies increase [10], and the oscillation periods of the carriers become comparable with the momentum and energy-relaxation-time constants. Therefore, hydrodynamic models (HDM) have been employed to solve mass, momentum, and energy conservation equations for the carriers in the semiconductors, to provide better accuracy. Especially for reaching THz operation speeds and developing sub-micron devices, modeling inertia effects and ballistic transport with HDM is vital [17,18,[25][26][27][28][29][30]. Applying a number of moment-expansion operations on the Boltzmann Transport Equation (BTE) would lead to transport models with increasing accuracy and complexity. One of the widely accepted set of equations for HDM are given for unipolar charge transport (assuming only electrons are present), as in (4)-(6).
where the unknowns of the system are electron concentration n e , momentum density of electrons p, and the energy density of electrons w. Additionally, κ stands for the heatconduction coefficient (κ = κ 0 µ n0 nk 2 B T 0 as in the Wiedemann-Franz law), T n is the carrier temperature, v stands for the electron velocity, and E is the electric field (E = −∇ϕ) that exerts electric force on the electrons. A set of additional equations are given to relate the parameters as p = m e nv, and w = 3 2 nk B T + 1 2 m e n|v| 2 , where m e is the effective carrier mass. The complementary collision terms can be written as ∂n e ∂t c where τ p and τ w represent the momentum and energy relaxation times, respectively, and T 0 stands for ambient temperature. By assuming the geometry is 2-D, and separating the components of the vectorial quantities, this set of HDM equations can be written in a more compact way: F(u) = a x v x u + [0, n e k B T n , 0, v x n e k B T n ] T + a y v y u + 0, 0, n e k B T n , v y n e k B T n T , for the unknown vector u = n e , p x , p y , w T . Solutions of this set for conventional semiconductors have already been demonstrated by the use of DGTD-FEM in 1-D [25] and 2-D [18]. However, due to the exotic properties of graphene, direct application of this model on graphene is challenging.

Graphene Properties and Parameters
The hydrodynamic properties of charge transport in graphene have been demonstrated multiple times, such as its viscosity effects [8,10], Poiseuille flow [9], etc. Moreover, macroscopic models based on the Boltzmann equation are developed for the analysis of dynamic effects in graphene, such as thermal transport [31]. However, a generalized HDM solver to model charge-carrier transport has been missing. The HDM-transport solver would also complete the charge-transport analysis in graphene under non-equilibrium conditions, when it is coupled with the nonlinear Poisson's solver described in the previous section. In other words, the full HDM solver, as described in (7)-(9), can be coupled with (2) and (3) to provide full analysis. Another approach could also be solving (2) and (3) initially to obtain the Fermi level in the graphene sheet, due to contact doping and gating, and (7)-(9) can be solved together with Poisson's equation for the transport analysis, assuming contact doping stays steady. However, the introduced set of HDM Equations (7)-(9) are derived for traditional semiconductors with well-defined effective carrier mass. Additionally, for traditional semiconductors, the momentum and energy-relaxation mechanisms are well understood, and the corresponding time constants can be found in the literature. Therefore, in order for HDM transport to be applied on graphene, material specific parameters, carrier effective mass, saturation velocity, momentum, and energy-relaxation-time constants need to be known. As is widely known, one of graphene's extreme properties is its exceptional carrier mass, due to its conical Dirac band structure. This massless characteristic of electrons and holes also results in very high mobility and electrical conductivity, making graphene a unique material. In this work, we follow the classical approach given in [32], without invoking the Dirac equation to obtain an average effective electron (and also hole) mass in graphene: where v F denotes Fermi velocity (= 8.96 × 10 4 m/s). For E F = 0, average effective carrier mass is also given with the following formula: Another very important parameter for hydrodynamic transport is the momentumrelaxation-time constant that reveals itself in the right-hand-side collision term of the momentum-conservation Equation (5). This parameter is primarily dependent on the scattering mechanisms taking place, and for high-temperature applications, phonon scattering would be the primary scattering mechanism. The momentum-relaxation-time constant for phonon scattering in graphene as a function of Fermi level, phonon velocities, mass density, and acoustic deformation potential can be obtained by following the formalism in [33,34]. For T = 300 K and E F = 100 meV, τ p corresponds to 0.27 ps, while it can potentially reach up to 2 ps for the undoped graphene with a Fermi level near the Dirac point.
Moreover, carrier-saturation velocity also plays a significant role for the transport and can be written as [35,36]: where ω OP ∼ = 200 meV is the inelastic optical phonon (OP) scattering rate, observed when SiO 2 is used as the gate-insulator material. The last missing parameter for the HDM of graphene is energy-relaxation time that appears in the collision term of the energy conservation Equation (6). For this work, it is taken as 1 ps for high-temperature operation [37,38]. Having obtained all the material-related transport parameters for graphene, the HDM set can be solved, and the hydrodynamic transport effects can be observed with the appropriate set of boundary conditions.

Results
In order to test and analyze the developed 2-D hydrodynamic solver for graphene, a 2-D device geometry, as in Figure 2, is chosen. The orange part in Figure 2 denotes the graphene-based device as well as the full computation domain for the HDM simulation. The geometry is inspired from a mechanical device designed by Nikola Tesla, also referred to as a Tesla valve [39]. The main purpose and advantage of this device is that when a fluid flows through this device, it encounters different resistances depending on the flow direction. In other words, in one direction, (from a to b in Figure 2) the fluid is divided into two branches, and they recombine again with a large angle (opposing the main flow direction) creating whirlpools, loss of energy, and higher resistance for the flow. However, for the other direction (fluid flow from b to a), the flow again is divided into two branches, but it recombines later with a smoother angle without causing too much loss of energy. The authors of [10] have fabricated a similar Tesla-valve-like graphene device operating at low temperature, based on this phenomenon, where electron flow is more resistive in one direction with diode-like electrical rectification. Since a commercial hydrodynamic carrier transport solver for graphene is not available, the authors of [10] have performed simulations using fluid-dynamics solvers to show the effects qualitatively and also to verify their experimental findings. In this work, we also choose a similar structure with smaller dimensions to model the hydrodynamic charge transport in graphene at room temperature.
where ℏ ≅ 200 meV is the inelastic optical phonon (OP) scattering rate, observed when SiO2 is used as the gate-insulator material. The last missing parameter for the HDM of graphene is energy-relaxation time that appears in the collision term of the energy conservation Equation (6). For this work, it is taken as 1 ps for high-temperature operation [37,38]. Having obtained all the materialrelated transport parameters for graphene, the HDM set can be solved, and the hydrodynamic transport effects can be observed with the appropriate set of boundary conditions.

Results
In order to test and analyze the developed 2-D hydrodynamic solver for graphene, a 2-D device geometry, as in Figure 2, is chosen. The orange part in Figure 2 denotes the graphene-based device as well as the full computation domain for the HDM simulation. The geometry is inspired from a mechanical device designed by Nikola Tesla, also referred to as a Tesla valve [39]. The main purpose and advantage of this device is that when a fluid flows through this device, it encounters different resistances depending on the flow direction. In other words, in one direction, (from a to b in Figure 2) the fluid is divided into two branches, and they recombine again with a large angle (opposing the main flow direction) creating whirlpools, loss of energy, and higher resistance for the flow. However, for the other direction (fluid flow from b to a), the flow again is divided into two branches, but it recombines later with a smoother angle without causing too much loss of energy. The authors of [10] have fabricated a similar Tesla-valve-like graphene device operating at low temperature, based on this phenomenon, where electron flow is more resistive in one direction with diode-like electrical rectification. Since a commercial hydrodynamic carrier transport solver for graphene is not available, the authors of [10] have performed simulations using fluid-dynamics solvers to show the effects qualitatively and also to verify their experimental findings. In this work, we also choose a similar structure with smaller dimensions to model the hydrodynamic charge transport in graphene at room temperature.  Ohmic contacts are depicted with "a" and "b", the charge carriers flow depending on the applied potential, either from "a" to "b" or from "b" to "a". The boundaries are insulating, except for the ohmic-metallic contacts. Hard and easy flow directions are given by the arrows.

Rectification: Tesla Valve
The developed HDM solver for graphene is applied for the device geometry given in Figure 2, which represents the fluidic device, called a Tesla valve, in nanoscale for electron transport in graphene. In order to observe the pure viscosity effects and the working principles of a Tesla valve, a homogeneous Fermi level distribution is assumed on the graphene, which can be provided by gating from underneath. Moreover, ideal ohmic contacts without contact doping effects are employed for the boundary conditions at the "a" and "b" boundaries in Figure 2. Assuming homogeneous E F = 100 meV on the graphene layer, a unipolar HDM solver (assuming only electrons exist) is employed, together with the Poisson's equation to simulate the structure. The electric field (driving force for the carriers) between the ideal ohmic contacts are obtained by solving the Poisson's equation, with respective Dirichlet boundary conditions. In other words, the model (7)-(9) is coupled with a Poisson's equation solver in lateral dimensions, for the computation of the hydrodynamic charge carrier transport in graphene. This analysis is done independently of the transverse potential computation, thanks to the homogeneous Fermi level assumption. For a more accurate analysis, Poisson's equation can be considered in 3-D, to account for both lateral and transverse variations at the same time.
The narrow width, in the order of 10 nm, of the graphene device in Figure 2, might also lead to edge effects and confinement effects, which could potentially change the electron velocities [31,40,41]. In this work, we focus on macroscopic carrier transport and collective carrier motion, and the mentioned effects are omitted by using average values. However, thanks to the flexible nature of HDM equations, the varying carrier velocity and different phonon-scattering times can also be incorporated into the model for more accurate graphene-nanoribbon simulations.
There are two types of boundary conditions used in this structure, the first one is the ohmic contact condition on the boundaries 'a' and 'b', to make sure the potential is applied on electrons, so they can enter and leave the structure without resistance. The other boundary condition is applied on the edges of the Tesla valve structure. For that purpose, an isolation boundary is used to eliminate the possible flow outside of the computational domain. This isolation condition is implemented weakly, by assuming zero velocity just outside of the graphene device. In other words, usage of the DGTD method allowed us to assign zero velocity for the carriers on the "ghost elements", just outside the device domain, and their presence is felt through the numerical-flux term in the DGTD-FEM formalism [18]. For the numerical flux, local Lax Friedrichs flux has been utilized [42], and time steps as small as 0.01 fs are needed to provide the stability, due to high carrier velocities and concentrations. The developed solver takes around 290 ms CPU time for the computation of each time step on an Intel i7-7700, requiring less than 100 MB memory for the device geometry in Figure 2, with triangular mesh consisting of 21,000 triangular elements. To sum up, for simulating 1 ps of dynamic transport, usually there is a required duration for the steady state to be formed: the solver computes 100,000 time steps. Thanks to the discontinuous nature of the employed DGTD method, a parallel approach can be developed for the next-generation implementation, and the resulting computation cost and required CPU time can be relieved. This would also lead to acceptable computation durations for the simulations of bigger structures and devices made out of graphene.
Electron flow directions and velocities can be seen in Figure 3, for the cases where the contact on the right is grounded, and 50 mV is applied to the left contact, and vice versa. As it can be seen from Figure 3a, electrons meet again without too much resistance, when they flow in the easy direction, whereas electrons with opposing flow directions cause loss of energy in the hard-flow direction, as can be seen from Figure 3b. This different resistance experienced by the electrons, depending on the flow direction, also shows itself in the total current flowing through the device. Figure 4 depicts the transient current vs. time and the transient potential on the left and right contacts. When the left contact is kept at 50 mV, electrons flow toward the left, and they generate a higher absolute current, compared to the case where the right contact is kept at 50 mV and the left one is grounded. The diodicity parameter can be defined for this device similar to [10] as The diodicity parameter can be defined for this device similar to [10] as For the same level of bias applied from the contacts in both the easy-and hard-flow cases, the diodicity yields to = 1.0655. Higher diodicities can be easily achieved by cascading several of these Tesla valves elements.  . Transient voltage and current profiles of the Tesla valve during "easy" and "hard" flow modes. Between 1-2 ps, a bias in the 'hard' flow direction is applied, i.e., the right contact is biased The diodicity parameter can be defined for this device similar to [10] as For the same level of bias applied from the contacts in both the easy-and hard-flow cases, the diodicity yields to = 1.0655. Higher diodicities can be easily achieved by cascading several of these Tesla valves elements.  . Transient voltage and current profiles of the Tesla valve during "easy" and "hard" flow modes. Between 1-2 ps, a bias in the 'hard' flow direction is applied, i.e., the right contact is biased Figure 4. Transient voltage and current profiles of the Tesla valve during "easy" and "hard" flow modes. Between 1-2 ps, a bias in the 'hard' flow direction is applied, i.e., the right contact is biased to 50 mV while the left contact is kept grounded. Between 3-4 ps, the biased is reversed to facilitate current in 'easy' direction. The computed current level between 3-4 ps is higher than the current levels observed in the 'hard' flow direction between 1-2 ps.
For the same level of bias applied from the contacts in both the easy-and hard-flow cases, the diodicity yields to D = 1.0655. Higher diodicities can be easily achieved by cascading several of these Tesla valves elements.

Poiseuille Flow
The classical Poiseuille flow is a well-known phenomenon, due to hydrodynamic effects and the collective motion that occurs when fluids flow in a long duct, such as a pipe. It has also been shown that under certain circumstances, the electrons in graphene also demonstrate the Poiseuille flow, while carrying a current [8,9]. In order to account for this effect, modified isolation boundary conditions have been used in the HDM solver, by enforcing the velocity of carriers to be zero, just outside the computational domain (in the 'ghost' elements). Therefore, the velocity profiles in the narrow channels, as in Figure 5, can be obtained; this ensures that the current is the maximum in the middle of the narrow graphene channels, and the electron flow tries to avoid the edges. Moreover, the current levels are higher when the electrons are flowing in the easier direction, which is labeled as 'forward', and maximum of the current is less in the 'backward' or harder-flow direction. current in 'easy' direction. The computed current level between 3-4 ps is higher than the current levels observed in the 'hard' flow direction between 1-2 ps.

Poiseuille Flow
The classical Poiseuille flow is a well-known phenomenon, due to hydrodynamic effects and the collective motion that occurs when fluids flow in a long duct, such as a pipe. It has also been shown that under certain circumstances, the electrons in graphene also demonstrate the Poiseuille flow, while carrying a current [8,9]. In order to account for this effect, modified isolation boundary conditions have been used in the HDM solver, by enforcing the velocity of carriers to be zero, just outside the computational domain (in the 'ghost' elements). Therefore, the velocity profiles in the narrow channels, as in Figure 5, can be obtained; this ensures that the current is the maximum in the middle of the narrow graphene channels, and the electron flow tries to avoid the edges. Moreover, the current levels are higher when the electrons are flowing in the easier direction, which is labeled as 'forward', and maximum of the current is less in the 'backward' or harder-flow direction.

Discussion
Two separate solvers to model and simulate the charge behavior in graphene have been developed. The first simulator combines the nonlinear Poisson's equation with the Thomas-Fermi equation in graphene, for the computation of the Fermi level and carrier distribution under the gate potential and contact doping, due to the metal layers. An example scenario is given in Figure 1 that analyzes the static fields in the transverse direction. The second solver considers the dynamic behavior of electrons and holes in graphene, by solving mass, momentum, and energy-conservation equations together. The proposed transport solver reveals the hydrodynamic charge transport in the lateral dimension and is qualitatively validated by simulating a graphene Tesla valve device. The modified HDM in graphene can account for the viscosity effects and collective carrier motion, and this enables modeling the hydrodynamic effects occurring in graphene-based devices, together with the Poisson's solver. Thanks to the modified isolation boundary condition, the Poiseuille flow can also be modeled, which is effective especially for the

Discussion
Two separate solvers to model and simulate the charge behavior in graphene have been developed. The first simulator combines the nonlinear Poisson's equation with the Thomas-Fermi equation in graphene, for the computation of the Fermi level and carrier distribution under the gate potential and contact doping, due to the metal layers. An example scenario is given in Figure 1 that analyzes the static fields in the transverse direction. The second solver considers the dynamic behavior of electrons and holes in graphene, by solving mass, momentum, and energy-conservation equations together. The proposed transport solver reveals the hydrodynamic charge transport in the lateral dimension and is qualitatively validated by simulating a graphene Tesla valve device. The modified HDM in graphene can account for the viscosity effects and collective carrier motion, and this enables modeling the hydrodynamic effects occurring in graphene-based devices, together with the Poisson's solver. Thanks to the modified isolation boundary condition, the Poiseuille flow can also be modeled, which is effective especially for the narrow carrier flow channels. Finally, the results obtained by the proposed models match, qualitatively, with the experimental findings of the literature for both the electrostatic case with contact doping and the hydrodynamic transport seen in the Tesla valve.
The developed solvers can be adopted for the development of next-generation graphene devices, especially for THz applications, by modeling and making use of the carrier inertia effects, viscosity effects, and possible ballistic charge transport. Furthermore, graphene nanoribbons (GNR) can also be modeled with this approach, by using appropriate carrier velocity and phonon-scattering times, emerging due to edge effects. Moreover, a parallel solver for the HDM can be implemented in the next generation, thanks to the employed DGTD algorithm.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/ma15124141/s1, Figure S1: 2-D device geometry (cross-section) with two gates at the top and the bottom. Graphene is depicted as orange layer in the middle, the green layers represent dielectrics, the gate contacts are shown in golden, and Dirichlet conditions are applied to those boundaries for solving the nonlinear Poisson's equation; Figure S2. Fermi level and charge carrier profile on graphene layer with respect to changing top gate potentials while the bottom gate is kept grounded. (a,c,e,g,i) show the Fermi level profile for the varying gate potentials (50 mV, 100 mV, 500 mV, 1 V, 5 V respectively) computed by using Equations (S1) and (S4). The full Thomas-Fermi dynamics in (S1) considers temperature dependence (black curves) while (S4) assumes 0 K (purple curves). (b,d,f,h,j) depict the carrier distribution on the graphene layer for the varying gate potentials. Blue, orange and dashed black curves are obtained from the full model (1)-(3), the purple curve for the electron concentration assumes 0 K based on (S4). As the gate potential, and the Fermi level, increases the computed distributions from both approaches come closer to each other, moreover hole concentration becomes negligible implying the 0 K approach could be used and poly-logarithm function can be avoided. Funding: The work presented in this paper is funded partially from the European Union's Horizon 2020 rEuropean Union's Horizon 2020 (H2020-NMBP-07-2017) under grant agreement MMAMA n • 761036.