Simulation of Z-Shaped Graphene Geometric Diodes Using Particle-in-Cell Monte Carlo Method in the Quasi-Ballistic Regime

Geometric diodes are planar conductors patterned asymmetrically to provide electrical asymmetry, and they have exhibited high-frequency rectification in infrared rectennas. These devices function by ballistic or quasi-ballistic transport in which the transport characteristics are sensitive to the device geometry. Common methods for predicting device performance rely on the assumption of totally ballistic transport and neglect the effects of electron momentum relaxation. We present a particle-in-cell Monte Carlo simulation method that allows the prediction of the current–voltage characteristics of geometric diodes operating quasi-ballistically, with the mean-free-path length shorter than the critical device dimensions. With this simulation method, we analyze a new diode geometry made from graphene that shows an improvement in rectification capability over previous geometries. We find that the current rectification capability of a given geometry is optimized for a specific mean-free-path length, such that arbitrarily large mean-free-path lengths are not desirable. These results present a new avenue for understanding geometric effects in the quasi-ballistic regime and show that the relationship between device dimensions and the carrier mean-free-path length can be adjusted to optimize device performance.


Introduction
We defined geometric diodes as planar conductors patterned with a geometric asymmetry that gives rise to a preferred current direction in 2008 [1], and we have demonstrated that they rectify at DC and infrared frequencies [2,3]. Since these diodes are planar and therefore have extremely low capacitance, they have the potential to provide ultra-fast rectification. They can be used in rectennas, diodes coupled to antennas, to rectify infrared signals and produce DC power [3], and they can also be used in other applications that demand high-speed electronics. Other high-frequency diodes, such as metal-insulator-insulator-metal diodes, have substantially higher capacitance due to their parallel-plate configuration and hence lower operating frequencies [4] even for ideal materials [5]. In addition, geometric diodes may be able to achieve higher current asymmetries than other high-frequency diodes [6].
The mean-free-path length (λ MFP ) of charge carriers in the geometric diode material is critical. For a geometric diode to exhibit rectification, charge carriers must travel ballistically (λ MFP greater than device dimensions) or quasi-ballistically (λ MFP near device dimensions) through the patterned material, which requires that device dimensions be on the order of λ MFP or lower. This imposes difficult fabrication requirements for devices made of conventional metals, as typical mean-free-path lengths are below 100 nm. This difficulty can be relaxed somewhat by using more exotic materials, such as graphene, which has shown mean-free-path lengths on the micron scale [7]. Such graphene devices have been demonstrated experimentally [2,3,8] and are capable of infrared rectification. Methods of modeling geometric diodes often incorporate the Landauer-Buttiker formalism to obtain device current-voltage (I(V)) characteristics, using a transmission function obtained theoretically [9] or through simulation [10,11]. This approach is elegant but so far has been applied only to purely ballistic transport and does not include the momentum relaxation that is present in quasi-ballistic devices, which can still exhibit rectifying behavior. Another simulation method incorporates a Monte Carlo algorithm with Drude-like carrier transport through the devices, and it includes momentum relaxation by way of scattering events [3,12]. In this approach, the device geometry is established, and the input voltage is applied across the device at the two electrodes in order to compute the geometry-dependent electric field, which is assumed to be constant throughout the simulation time. Then, a certain number of charge carriers, each of which represents multiple electrons, are generated at random locations in the device, each with a randomly directed velocity with magnitude equal to the Fermi velocity (v F ). Then, the carriers are allowed to move under the influence of the electric field over a small timestep, with the electric field perturbing their initial randomly directed momentum. If a charge carrier collides with a device boundary, it undergoes specular reflection and is allowed to continue propagating. If it crosses either of the electrodes, then a counter is incremented before another carrier is injected from the opposite electrode. Once a certain carrier has traveled a distance equal to the mean-free-path length, it is given a new randomly directed velocity, with magnitude v F , and the process repeats for a predetermined number of time steps. After a chosen time has elapsed, the current is derived from the counter and the total simulation time.
We present a new simulation method that incorporates a particle-in-cell algorithm (PIC) [13] to determine the I(V) characteristic for a given device. This method, which builds on the previous Monte Carlo method [3], decreases the computation time and includes a self-consistent electric field computation at every time step. Our application of this method also includes a more accurate implementation of the collision time statistics from the Drude model. Since this algorithm includes momentum relaxation, it is useful for the simulation of both ballistic and quasi-ballistic devices. Our previous simulations and demonstrations made use of an inverse-arrowhead design [3]. Here, we use this simulator to demonstrate a new diode geometry, the Z-diode, which shows improved current asymmetry over previous designs.

Materials and Methods
The simulation begins by establishing the diode boundaries in MATLAB [14] and creating a triangular mesh within them. The nodes of the mesh are the discrete points at which the charge concentration, electric potential, and electric field will be computed later in the algorithm. Then, a certain number of macroparticles, each of which can represent multiple electrons, are scattered within the device at random locations. The use of macroparticles, where each one typically represents 100 or 200 individual electrons depending on the carrier concentration of the material, facilitates quicker computation than simulating individual carriers. Each macroparticle is initially distributed to the three mesh nodes surrounding it using a weighting determined from the relative distances between the particle and each of the nodes. This establishes the initial spatial charge density. There is assumed to be a homogenous distribution of background positive charge that maintains charge neutrality in the device as a whole. Each of the macroparticles is initiated with a randomly directed velocity vector with magnitude equal to v F . The probabilistic nature of electron collisions is incorporated by assigning each macroparticle a time interval in the future at which it will undergo a momentum relaxation event. These times are sampled from an exponential decay distribution where e −t/τ represents the probability that an electron picked at random will have no collision during the next time interval t, and τ = λ MFP /v F represents the average momentum relaxation time of the choice material [15].
The input voltage is set as a boundary condition between the two electrodes. MATLAB's PDE package is used to solve Poisson's equation for the potential at every mesh node. Then, the electric field is computed from the gradient of the solution. Then, the field at each particle is determined from the field at the surrounding three nodes using the distance weighting that was used for the charge distribution. The field exerts a force on the particle, changing its velocity over a small time-step dt, and the particles hop forward in time along the new velocity vector. The counter that tracks the time at which a given particle will undergo a momentum-destroying collision is decremented by dt for each particle. The new locations of the particles are used to compute the charge density again, and the process repeats itself many times. This methodology, which is usually incorporated in PIC algorithms, computes inter-particle interactions through the spatial charge density rather than through Coulombic forces between particle pairs, which would be more computationally expensive. Figure 1 shows the setup and results graphically for an inverse-arrowhead diode shape [3]. Note that for clarity, the figure contains only 20 macroparticles, but in full simulations, hundreds are used. For the simulations described herein, we always ensured that a sufficient number of macroparticles were used so that the results did not change with a further increase in the number of macroparticles. For a lower number of macroparticles, each macroparticle represents a larger number of individual electrons. Since inter-electron forces are not considered among electrons within the same macroparticle, more electron interactions are ignored and care must be given to ensure error is not introduced. The input voltage is set as a boundary condition between the two electrodes. MATLAB's PDE package is used to solve Poisson's equation for the potential at every mesh node. Then, the electric field is computed from the gradient of the solution. Then, the field at each particle is determined from the field at the surrounding three nodes using the distance weighting that was used for the charge distribution. The field exerts a force on the particle, changing its velocity over a small time-step dt, and the particles hop forward in time along the new velocity vector. The counter that tracks the time at which a given particle will undergo a momentum-destroying collision is decremented by dt for each particle. The new locations of the particles are used to compute the charge density again, and the process repeats itself many times. This methodology, which is usually incorporated in PIC algorithms, computes inter-particle interactions through the spatial charge density rather than through Coulombic forces between particle pairs, which would be more computationally expensive. Figure 1 shows the setup and results graphically for an inverse-arrowhead diode shape [3]. Note that for clarity, the figure contains only 20 macroparticles, but in full simulations, hundreds are used. For the simulations described herein, we always ensured that a sufficient number of macroparticles were used so that the results did not change with a further increase in the number of macroparticles. For a lower number of macroparticles, each macroparticle represents a larger number of individual electrons. Since inter-electron forces are not considered among electrons within the same macroparticle, more electron interactions are ignored and care must be given to ensure error is not introduced. At any step, if the time-to-next-collision counter goes below zero for any particle, that particle's momentum is reset by giving it a new randomly-directed velocity with magnitude vF, and its counter is reset from a sampling of the exponential decay distribution. The new counter value is reduced by the amount of time the particle traveled after the collision event, i.e., the previous negative value of the counter. This ensures that the momentum relaxation statistics are obeyed throughout the simulation and, on average, the particles will travel approximately a distance of λMFP between collisions. If a particle collides with a device boundary, it is reflected specularly and continues propagating. When a particle crosses one of the electrode boundaries, a counter, Q, which counts the total number of charges to pass through the device throughout the entire simulation, is incremented or decremented accordingly. Then, the total current, I, is computed as, At any step, if the time-to-next-collision counter goes below zero for any particle, that particle's momentum is reset by giving it a new randomly-directed velocity with magnitude v F , and its counter is reset from a sampling of the exponential decay distribution. The new counter value is reduced by the amount of time the particle traveled after the collision event, i.e., the previous negative value of the counter. This ensures that the momentum relaxation statistics are obeyed throughout the simulation and, on average, the particles will travel approximately a distance of λ MFP between collisions. If a particle collides with a device boundary, it is reflected specularly and continues propagating. When a particle crosses one of the electrode boundaries, a counter, Q, which counts the total number of charges to pass through the device throughout the entire simulation, is incremented or decremented accordingly. Then, the total current, I, is computed as, where n is the material's free carrier concentration per unit area, A is the total device area, N macro is the number of macroparticles used in the simulation, q is the elementary charge, and N t is the total number of simulated time steps.
It is important to note that dt must be chosen carefully with regard to τ. Since a particle's collision counter is decremented by dt at every timestep, for larger dt, there is a larger margin by which a particle's collision counter can drop below zero. This corresponds to the particle overshooting its actual collision location, which results in an unintentional increase in the distance traveled between collisions. Then, the simulated current will be overstated due to an artifactual increase in mobility. To avoid this problem, mesh elements must be small enough to capture the large potential gradients that occur at geometric constrictions in the device (see the large electric field acting on the particle near the constriction at x = 0 in Figure 1b). If the mesh is too coarse, large potential gradients (electric fields) will be understated and yield a lower simulated current. On the other hand, arbitrarily small timesteps and mesh sizes require long simulation times. For each geometry, we repeated simulations for decreasing dt and mesh size to find values for each where the resulting I(V) remained approximately constant with further reductions. Once we determined a sufficiently small timestep and mesh size to give a grid-independent solution for a given geometry, we used them for all subsequent simulations of that geometry. A time step of dt = τ/10 was used for all simulations herein. For devices with geometric constrictions, the mesh element size varied throughout with a larger mesh size near the electrodes and a smaller mesh size near the constrictions.

Rectangular Strip
To verify that the simulator accurately predicts current-voltage characteristics, it was run on a rectangular sample with no constrictions and then compared to the simple analytical solution from the Drude model [15]. Figure 2 shows the simulated and modeled current, I, for applied voltage V 0 . The mesh size was set to 250 nm and constant throughout the device, since the electric field was uniform due to the absence of constrictions. The simulated result agrees well with the analytical solution, which supports the accuracy of the simulation. where n is the material's free carrier concentration per unit area, A is the total device area, Nmacro is the number of macroparticles used in the simulation, q is the elementary charge, and Nt is the total number of simulated time steps.
It is important to note that dt must be chosen carefully with regard to τ. Since a particle's collision counter is decremented by dt at every timestep, for larger dt, there is a larger margin by which a particle's collision counter can drop below zero. This corresponds to the particle overshooting its actual collision location, which results in an unintentional increase in the distance traveled between collisions. Then, the simulated current will be overstated due to an artifactual increase in mobility. To avoid this problem, mesh elements must be small enough to capture the large potential gradients that occur at geometric constrictions in the device (see the large electric field acting on the particle near the constriction at x = 0 in Figure 1b). If the mesh is too coarse, large potential gradients (electric fields) will be understated and yield a lower simulated current. On the other hand, arbitrarily small timesteps and mesh sizes require long simulation times. For each geometry, we repeated simulations for decreasing dt and mesh size to find values for each where the resulting I(V) remained approximately constant with further reductions. Once we determined a sufficiently small timestep and mesh size to give a grid-independent solution for a given geometry, we used them for all subsequent simulations of that geometry. A time step of dt = τ/10 was used for all simulations herein. For devices with geometric constrictions, the mesh element size varied throughout with a larger mesh size near the electrodes and a smaller mesh size near the constrictions.

Rectangular Strip
To verify that the simulator accurately predicts current-voltage characteristics, it was run on a rectangular sample with no constrictions and then compared to the simple analytical solution from the Drude model [15]. Figure 2 shows the simulated and modeled current, I, for applied voltage V0. The mesh size was set to 250 nm and constant throughout the device, since the electric field was uniform due to the absence of constrictions. The simulated result agrees well with the analytical solution, which supports the accuracy of the simulation.

Inverse-Arrowhead and Z-Diode
After demonstration of the inverse-arrowhead diode, Zhu proposed that alternative geometries could present increased current asymmetry, with one in particular being the Z-diode [12], which is compared to the inverse-arrowhead in Figure 3. In both the inverse-Nanomaterials 2021, 11, 2361 5 of 9 arrowhead (Figure 3a) and the Z-diode (Figure 3b), under forward bias, charges drift from left to right, and the geometry is such that reflections off the device boundaries funnel the charges through the constriction, or neck. However, in reverse bias, there is a notable difference in the geometries of the two devices. As charges drift right to left, the inversearrowhead geometry allows for easier transmission through the neck because the opening is perpendicular to the general drift direction. This reverse leakage is reduced in the Z-diode because the neck opening is oriented perpendicular to the general drift direction.

Inverse-Arrowhead and Z-Diode
After demonstration of the inverse-arrowhead diode, Zhu proposed that alternative geometries could present increased current asymmetry, with one in particular being the Z-diode [12], which is compared to the inverse-arrowhead in Figure 3. In both the inversearrowhead (Figure 3a) and the Z-diode (Figure 3b), under forward bias, charges drift from left to right, and the geometry is such that reflections off the device boundaries funnel the charges through the constriction, or neck. However, in reverse bias, there is a notable difference in the geometries of the two devices. As charges drift right to left, the inversearrowhead geometry allows for easier transmission through the neck because the opening is perpendicular to the general drift direction. This reverse leakage is reduced in the Zdiode because the neck opening is oriented perpendicular to the general drift direction. To compare the rectification capabilities of the two geometries, their I(V) characteristics were simulated. The relevant geometric variables are illustrated in Figure 4. The shoulders, i.e., the electrodes from which charges are injected, have width ds, while the neck, or the minimum constriction in the device, has width dn. The voltage, V0, is applied between the left and right electrodes to drive the current, I, through the device. The strength of the geometric rectification effect depends on the device dimensions relative to λMFP [12]. Both simulations used λMFP = 500 nm, dn = 100 nm, and ds = 1 μm. Voltages were swept from −0.5 to 0.5 V.  The material was assumed to be graphene for the simulation, such that the material properties can be tuned by applying a gate bias. Much graphene processing is carried out on thermally oxidized silicon wafers. In this case, the graphene device sits atop the oxide, To compare the rectification capabilities of the two geometries, their I(V) characteristics were simulated. The relevant geometric variables are illustrated in Figure 4. The shoulders, i.e., the electrodes from which charges are injected, have width d s , while the neck, or the minimum constriction in the device, has width d n . The voltage, V 0 , is applied between the left and right electrodes to drive the current, I, through the device. The strength of the geometric rectification effect depends on the device dimensions relative to λ MFP [12]. Both simulations used λ MFP = 500 nm, d n = 100 nm, and d s = 1 µm. Voltages were swept from −0.5 to 0.5 V.

Inverse-Arrowhead and Z-Diode
After demonstration of the inverse-arrowhead diode, Zhu proposed that alternative geometries could present increased current asymmetry, with one in particular being the Z-diode [12], which is compared to the inverse-arrowhead in Figure 3. In both the inversearrowhead (Figure 3a) and the Z-diode (Figure 3b), under forward bias, charges drift from left to right, and the geometry is such that reflections off the device boundaries funnel the charges through the constriction, or neck. However, in reverse bias, there is a notable difference in the geometries of the two devices. As charges drift right to left, the inversearrowhead geometry allows for easier transmission through the neck because the opening is perpendicular to the general drift direction. This reverse leakage is reduced in the Zdiode because the neck opening is oriented perpendicular to the general drift direction. To compare the rectification capabilities of the two geometries, their I(V) characteristics were simulated. The relevant geometric variables are illustrated in Figure 4. The shoulders, i.e., the electrodes from which charges are injected, have width ds, while the neck, or the minimum constriction in the device, has width dn. The voltage, V0, is applied between the left and right electrodes to drive the current, I, through the device. The strength of the geometric rectification effect depends on the device dimensions relative to λMFP [12]. Both simulations used λMFP = 500 nm, dn = 100 nm, and ds = 1 μm. Voltages were swept from −0.5 to 0.5 V.  The material was assumed to be graphene for the simulation, such that the material properties can be tuned by applying a gate bias. Much graphene processing is carried out on thermally oxidized silicon wafers. In this case, the graphene device sits atop the oxide, The material was assumed to be graphene for the simulation, such that the material properties can be tuned by applying a gate bias. Much graphene processing is carried out on thermally oxidized silicon wafers. In this case, the graphene device sits atop the oxide, and a gate voltage is applied to the silicon layer below the oxide. The dependence of carrier concentration on gate voltage can be described by n = 0 r V g /tq [16], where 0 is the permittivity of free space, r = 3.9 is the relative permittivity of SiO 2 , V g is the gate voltage relative to the Dirac point where the carrier concentration is minimum, t is the thickness of the SiO 2 substrate, and q is the elementary charge. V g can be directly applied by biasing the substrate, but is often nonzero even without bias due to the unintentional impurity doping that occurs in fabrication [3,16,17]. In such cases, it is not uncommon for V g to be in the order of 10 V, so we set V g = 10 V for our simulations. The oxide thickness is set to t = 300 nm, since this is a common SiO 2 thickness used in graphene processing. For these parameters, electrons are the majority carriers, and their concentration is n = 7.2 × 10 11 cm −2 . The effective mass, m * , follows the relation, where h is Planck's constant, and m 0 is the electron mass [18]. For these simulations, this gives m * = 0.02. For all the simulations given here, we used a time step of dt = τ/10. The maximum mesh size was set to 100 nm at the shoulders and 40 nm at the neck. Reducing either the time step or the maximum mesh size further resulted in longer computation times for the same I(V) curves. The resulting I(V) curves for both diodes, as well as their current asymmetries, are shown in Figure 5. Both diodes show rectification behavior in the quasi-ballistic regime, where the charge carriers undergo momentum relaxation events. As expected, the Z-diode exhibits a more nonlinear I(V) curve and a higher current asymmetry above 0.1 V. This is predominantly due to the stronger suppression of the reverse current due the geometry being more effective at reflecting charge carriers away from the neck in reverse bias. Although it is not yet clear why, this Z-diode is less asymmetric at 0.1 V and below, indicating the existence of a turn-on voltage that must be supplied for substantial rectification. For the given dimensions, the Z-diode shows a higher resistance at low biases, which must be considered when impedance matching devices to other circuit components, such as antennas in rectenna circuits. In a recent study by Wang et. al, an inverse-arrowhead diode was fabricated and measured with results showing correspondence to the simulator [18].
rier concentration on gate voltage can be described by n = ϵ 0 ϵ r V g /tq [16], where ϵ 0 is the permittivity of free space, ϵ r 3.9 is the relative permittivity of SiO2, V g is the gate voltage relative to the Dirac point where the carrier concentration is minimum, t is the thickness of the SiO2 substrate, and q is the elementary charge. Vg can be directly applied by biasing the substrate, but is often nonzero even without bias due to the unintentional impurity doping that occurs in fabrication [3,16,17]. In such cases, it is not uncommon for Vg to be in the order of 10 V, so we set Vg = 10 V for our simulations. The oxide thickness is set to t = 300 nm, since this is a common SiO2 thickness used in graphene processing. For these parameters, electrons are the majority carriers, and their concentration is n = 7.2 × 10 11 cm -2 . The effective mass, m * , follows the relation, where h is Planck's constant, and m 0 is the electron mass [18]. For these simulations, this gives m * = 0.02. For all the simulations given here, we used a time step of dt = τ/10. The maximum mesh size was set to 100 nm at the shoulders and 40 nm at the neck. Reducing either the time step or the maximum mesh size further resulted in longer computation times for the same I(V) curves.
The resulting I(V) curves for both diodes, as well as their current asymmetries, are shown in Figure 5. Both diodes show rectification behavior in the quasi-ballistic regime, where the charge carriers undergo momentum relaxation events. As expected, the Z-diode exhibits a more nonlinear curve and a higher current asymmetry above 0.1 V. This is predominantly due to the stronger suppression of the reverse current due the geometry being more effective at reflecting charge carriers away from the neck in reverse bias. Although it is not yet clear why, this Z-diode is less asymmetric at 0.1 V and below, indicating the existence of a turn-on voltage that must be supplied for substantial rectification. For the given dimensions, the Z-diode shows a higher resistance at low biases, which must be considered when impedance matching devices to other circuit components, such as antennas in rectenna circuits. In a recent study by Wang et. al, an inversearrowhead diode was fabricated and measured with results showing correspondence to the simulator [18].

Asymmetry Variation with λ MFP
In the fully ballistic regime, λ MFP is assumed to be larger than all device dimensions. In the quasi-ballistic regime, this is not necessarily the case, and so changing the device size or λ MFP changes the prevalence of momentum-destroying collisions during electron transport. To understand this effect, the previous Z-diode was simulated over a range of λ MFP . For each value of λ MFP , the currents were simulated for applied voltages of V 0 = ±0.5 V, and the resulting asymmetry was calculated, with the results shown in Figure 6.
As λMFP approaches zero, the asymmetry approaches 1, which is what is expected for macroscopic conductors with large critical dimensions relative to λMFP that do not exhibit quasi-ballistic transport. In such a case, the carrier momentum is destroyed too often for transport to be dominated by the device geometry. However, the most interesting revelation is that the current asymmetry peaks for a finite mean-free-path length, which has previously not been discussed in literature. In the case of the Z-diode with a 100 nm neck, the peak asymmetry appears to occur around λ peak ≈ 600 nm, with asymmetry dropping for λ MFP > λ peak . This effect may be due to carriers with increasingly large λMFP undergoing more reflections at device boundaries under reverse bias before having their momenta reset. A charge carrier that is initially deflected away from the neck can reflect off the bottom boundary and be directed toward the neck again by the strong electric field. If the carrier undergoes no momentum relaxation in the time frame, then its velocity will not be randomized before it travels through the neck and contributes to an increased reverse current. The value of λpeak is likely to change based on the applied voltage and the diode shape, and it may not exist at all for certain geometries.

Discussion
The simulation method described here presents several advantages compared to existing simulation methods. It is useful in studying the rectification behavior of diodes that operate in the quasi-ballistic regime and does not assume fully ballistic transport as do the other methods discussed in the introduction. In addition, it includes the effects of electric As λ MFP approaches zero, the asymmetry approaches 1, which is what is expected for macroscopic conductors with large critical dimensions relative to λ MFP that do not exhibit quasi-ballistic transport. In such a case, the carrier momentum is destroyed too often for transport to be dominated by the device geometry. However, the most interesting revelation is that the current asymmetry peaks for a finite mean-free-path length, which has previously not been discussed in literature. In the case of the Z-diode with a 100 nm neck, the peak asymmetry appears to occur around λ peak ≈ 600 nm, with asymmetry dropping for λ MFP > λ peak . This effect may be due to carriers with increasingly large λ MFP undergoing more reflections at device boundaries under reverse bias before having their momenta reset. A charge carrier that is initially deflected away from the neck can reflect off the bottom boundary and be directed toward the neck again by the strong electric field. If the carrier undergoes no momentum relaxation in the time frame, then its velocity will not be randomized before it travels through the neck and contributes to an increased reverse current. The value of λ peak is likely to change based on the applied voltage and the diode shape, and it may not exist at all for certain geometries.

Discussion
The simulation method described here presents several advantages compared to existing simulation methods. It is useful in studying the rectification behavior of diodes that operate in the quasi-ballistic regime and does not assume fully ballistic transport as do the other methods discussed in the introduction. In addition, it includes the effects of electric fields that result from inhomogeneous charge distributions and is computationally more efficient than existing Monte Carlo simulation methods.
It is important to note several factors and assumptions involved in the simulation that can compromise its validity in certain situations. This method assumes that the Fermi level of the charge carriers relative to the Dirac point, E F , is much larger than k B T. Here, k B is the Boltzmann constant and T is the ambient temperature. In this assumption, all of the carriers are assumed to have energy given by E F with negligible spread. In the case of graphene, the Fermi level is given by E F = m * m 0 v F 2 [19], which in our case gives E F = 100 meV, which is larger than k B T by a factor of 4 at room temperature.
This simulation method also neglects quantum effects due to the wave nature of electrons. The de Broglie wavelength can be used to determine the length scales below which wave-like properties of the charge carriers become important. From the relation E F =hkv F [19], whereh is the reduced Planck constant, and k is the carrier wave vector, we have k = 150 µm −1 for the above simulations, which translates to a de Broglie wavelength of 40 nm. This, along with the large coherence lengths reported for graphene [20], suggests that interference effects due to the wave nature of the electrons become important for device dimensions below the 100 nm scale. The de Broglie wavelength varies with the gate voltage, decreasing for greater V g , and can therefore be controlled to some extent.
At high electric fields, saturation effects due to high-energy charge carriers exciting optical phonons need to be accounted for [21]. For the parameters assumed in these simulations, i.e., V g = 10 V for graphene on a 300 nm thick SiO 2 substrate at room temperature, saturation effects become strong for applied fields around 0.5 V/µm [21]. For micron-sale devices, this simulation method should not be used with applied voltages exceeding the 1 V scale. The resulting current will be overstated as velocity saturation is ignored.

Conclusions
In our simulation approach, macroparticles are used to determine the potential and electric field at every mesh node in the device, and the particles move under the influence of the electrostatic forces. Each macroparticle's momentum is periodically destroyed, and its direction is randomized before moving under the influence of the electric fields. The simulator is run over many collisions, and the current is derived from the rate at which the macroparticles traverse the device. By including momentum relaxation, it is applicable to quasi-ballistic transport as well as ballistic transport, and therefore, it can be used to predict the behavior of a wider range of devices than if the momentum relaxation was not included. It is also computationally more efficient than existing Monte Carlo simulation methods. The simulator has certain limits to its applicability. For very small device features, below 100 nm for the material properties used here, quantum effects can begin to take hold and cause erroneous simulation results. For graphene, the simulator will also overstate the current for electric fields over 0.5 V/µm due to its neglect of saturation effects.
Through simulation, we have studied a new geometric diode design, the Z-diode, which exhibits improved rectification over 0.1 V compared to the previous inverse-arrowhead design. It was found that the current asymmetry of such diodes is sensitive to the device dimensions relative to the carrier mean-free-path length. In particular, for a given Z-diode size, there exists a finite mean-free-path length that maximizes the current asymmetry at a particular voltage, which is a result that has previously not been discussed in the literature. Equivalently, this suggests that given a material with a specified mean-free-path length, the critical dimensions of the diode can be tuned to maximize asymmetry.