Ambipolar device simulation based on the drift-diffusion model in ion-gated transition metal dichalcogenide transistors

Ionic gating is known as a powerful tool for investigation of electronic functionalities stemming from low voltage transistor operation to gate-induced electronic phase control including superconductivity. Two-dimensional (2D) material is one of the archetypal channel materials which exhibit a variety of gate-induced phenomena. Nevertheless, the device simulations on such ion-gated transistor devices have never been reported, despite its importance for the future design of device structures. In this paper, we developed a drift-diffusion (DD) model on a 2D material, WSe2 monolayer, attached with an ionic liquid, and succeeded in simulating the transport properties, potential profile, carrier density distributions in the transistor configuration. In particular, the simulation explains the ambipolar behavior with the gate voltage comparable to the band gap energy, as well as the formation of p-n junctions in the channel reported in several experimental papers. Such peculiar behavior becomes possible by the dramatic change of the potential profiles at the Schottky barrier by the ionic gating. The present result indicates that the DD model coupled to the Poisson equation is a fascinating platform to explain and predict further functionalities of ion-gated transistors through including the spin, valley, and optical degrees of freedom.

In particular, TMDs have gained interest thanks to the intriguing phenomena coupled the valley and spin degrees of freedom. Thus, TMDs can be applied to electronic and spintronic applications 16,17 . The TMD transistors are good candidates for the nextgeneration transistors since TMD semiconductors have enough band gap to switch ON/OFF states and high electron-hole mobilities 16,18 . Moreover, the short channel effect is suppressed compared to the Si devices since the 2D structure allows easier gate control 19 . The characteristics of the conventional metaloxide-semiconductor field-effect transistors (MOSFETs) are discussed mainly on MoS 2 18-24 , WSe 2 25-30 , and MoTe 2 31-34 . Ionic gating assists in realizing phenomena which are difficult to achieve in the conventional MOSFETs. The ambipolar behavior is one of those observed in ion-gated TMDs. In particular, hole transport in nominally n-type MoS 2 is observed only in EDLTs 24 . This ambipolar nature allows us to fabricate lateral p-n junctions in the channel without chemical doping and asymmetric metal electrodes 11 . By applying the source-drain bias voltage higher than the gate bias, the ionic distribution changes, allowing the ambipolar operation. Such ambipolar operation can be achieved also in the conventional transistors 25,[35][36][37] , but the operation voltage is significantly lower than that of the conventional FETs, and is in the voltage range similar to the bandgap energy. This offers an opportunity to determine the band gap directly from the transistor characteristics 38 . It is also known that the p-n junctions formed in the channel emits circularly polarized light, and furthermore, the circular polarization can be flipped by changing the direction of the injected current 13 . This enables the generation of the circularly polarized light using the electric field. This circular dichroism is ascribed to the imbalance of K and K′ points in the band structure caused by the voltage 39,40 . Moreover, the strong electric field induced by the ionic gating yields superconductivity in TMDs due to the high accumulation of carrier 1 .
To achieve more benefit of the ionic gating, the device operation mechanism of EDLTs should be understood. For theoretical studies of transistor operations, the drift-diffusion (DD) method is known as a powerful tool for device simulation to calculate the transport properties, band profile, and spatial distribution of carriers inside devices. As the down-scaling is employed to realize high performance in conventional transistors, the DD simulation is used to design ideal device structures 41 . However, there has been no DD simulation reported on the iongated transistors of TMDs.
In this work, we developed a 2D layer transistor model including an ionic liquid (IL) as a gate dielectric, based on the DD method and succeeded in simulating the distribution of ions, and the dynamics of electron and holes passing through the metal contacts. We examined the characteristics of ion-gated TMD transistors whose device size is about submicron while the several micron is the common channel size reported in experiments. The simulation shows that the down-scaled ion-gated WSe 2 transistor holds ambipolar behavior and the formation of the p-n junctions reported in larger-scale devices 13,15,38 . We clarified the transport 1 mechanism and the fundamental physics underneath the transistor operation using the band profile and the spatial distribution of ions and carriers obtained by the calculation.
The mechanism of the ambipolar transistor has been discussed theoretically on organic and amorphous materials because it is easy to inject both electrons and holes due to the in-gap states 42,43 . In single crystalline TMDs, such mid-gap states are absent. On the other hand, previous simulations on ion-gated organic transistors used a fixed capacitance induced by the ionic liquids 44,45 , which cannot include the variation of the capacitance due to the modification of the ion distribution by the source-drain bias voltage and the back action from the carriers in semiconductors. The general perspective of the ambipolar transport in MOSFETs has been discussed in ref. 46 .
Some specific implementations are required to adopt the DD model to the ion-gated transistors. We explain the points which differ from the DD model used for the conventional transistors in the next section. After describing the model, we show the results of device simulations on the ion-gated WSe 2 transistors. The current, band profile and ion/carrier density distributions are examined for the unipolar and ambipolar operation.

RESULTS AND DISCUSSION
Drift-diffusion model for ion-gated transition metal dichalcogenide transistors We start with the explanation of the DD model. We consider the 2D layer in (x, y) plane where the transport direction of the channel is set to x. The IL and oxide layer are piled up in the z direction as depicted in Fig. 1. We solve the continuity equation of each currents, 1 q ∇J n ðxÞ À RðxÞ ¼ 0; (1) À 1 q ∇J p ðxÞ À RðxÞ ¼ 0: Here, the electron and hole current densities are determined by well-known 1D Drift-Diffusion equations; J n ðxÞ ¼ Àqμ n nðxÞ∇Φðx; 0Þ þ qD n ðxÞ∇nðxÞ; J p ðxÞ ¼ Àqμ p pðxÞ∇Φðx; 0Þ À qD p ðxÞ∇pðxÞ; where J n(p) denotes the electron (hole) current densities. q is elementary charge, n(x) [p(x)] is the electron [hole] density at location (x, 0) in WSe 2 . n(x) and p(x) are defined in units of cm −2 . μ n(p) is the mobility of the electron (hole) and D n(p) (x) is the diffusion coefficient. Φ(x, 0) is the potential at position x. R(x) denotes the recombination rate, R(x) = R dir (x) + R SRH (x) + R n(p)−con (x). We consider two types of the recombination rate; the direct recombination rate, R dir (x), and the Shockley-Read-Hall (SRH) recombination rate, R SRH (x), where and R SRH ðxÞ ¼ nðxÞpðxÞ τ½nðxÞ þ pðxÞ : Here, β is the direct recombination parameter and τ is the recombination time of the SRH recombination. n in is the intrinsic carrier density of electrons and holes. For WSe 2 , n in ¼ where m e(h) is the effective mass of electrons (holes), m 0 is the mass of free electrons, and k B is the Boltzmann constant. R n(p)−con (x) indicates the generation or annihilation rate of the carriers which enter or exit from the channel at x through the thermionic or tunnel current. Thus, the boundary conditions for n(x) and p(x) are determined by the continuity equation in Eqs. (1), (2). The detailed description of R n(p)−con (x) is written in the Methods section.
We also solve the Poisson equation in the (x, z) plane in the whole system (WSe 2 , ILs, SiO 2 , and the metal contact) ∇ϵ α ϵ 0 Á ∇Φðx; zÞ ¼ ÀqCðx; zÞ; where Here, i + (x, z) and i − (x, z) are cation and anion concentrations, respectively (see Methods section for detail). d is the thickness of the WSe 2 layer, ϵ α is the relative permittivity of the material α (α ∈ WSe 2 , SiO 2 , and IL), and ϵ 0 is the permittivity of vacuum. ϵ WSe2 is direction dependent as shown in Table 1. The boundary condition is set to satisfy the continuity of the normal vector of the electric flux density at the interface between materials. The boundary condition at the metal  contact is set as Φ(x, t IL ) = V gs + Φ(x s , 0), where t IL is the thickness of IL. At the other edges, we adopt the Neumann's boundary condition, ∂Φð0;zÞ where L x is the length of WSe 2 layer in x direction, and t SiO2 is the thickness of the SiO 2 . The carrier densities [n(x) and p(x)] in the DD equation and the potential [Φ(x, z)] in the Poisson equation are determined self-consistently in the numerical calculations.
To adopt the above DD model to the ion-gated transistors, there are three points to be noticed regarding the difference between the conventional and ion-gated transistors. Firstly, in Eqs.
(3) and (4) the ordinary Einstein relation between the mobility and the diffusion coefficient is not appropriate for ion-gated transistors since the ions attract the high concentration of carriers above 10 13 cm −2 while the carrier concentration is usually below 10 12 cm −2 in the conventional transistors. Instead, we introduce the generalized Einstein relation (as explained in the Methods section) for the DD model at high carrier concentrations 47 .
The second point is the treatment of the contact. It is known that metal-TMD junctions usually show the Schottky behavior in the conventional TMD transistors [48][49][50] with some exceptions such as ones using graphene or Scandium contacts 20,51 . The mechanism of the Schottky transport has been precisely studied theoretically 52 . The conventional TMD transistors also have the problem of the Fermi-level-pinning due to the defect at the TMDmetal junctions 50 . Thus, only n-type operations have been reported in many works for conventional MoS 2 FETs. However, the Fermi-level pinning seems to be weak in the ion-gated TMDs, according to the experimental results 13,30 . We assume the contact made of an alloy of Ti and Au attached to semiconductor TMDs such as MoS 2 and WSe 2 acts as the n-type Schottky contact which only injects electrons. However, under the IL, as reported in ref. 13 , the carrier of unipolar transport changes from electrons to holes by changing the sign of the gate voltage for small source-drain voltage. Furthermore, the electron current shows the linear increase by applying the source-drain voltage, which indicates the Ohmic contact. Thus, the boundary condition for the Schottky contact is not appropriate for the DD model. Moreover, the boundary condition for the Ohmic contact which can inject either electron and hole is also not appropriate. To realize both Ohmic and ambipolar behavior, we adopt a practical Schottky contact model, which realizes not only the Schottky contact but also the Ohmic contact for both electron and hole carriers 53 . We expand the model in ref. 53 from 3D structure to the 2D structure that can be applied to thin layer devices.
Thirdly, since the ion has the finite size, there is a limitation of ions per volume in the liquid. Note that the potential becomes too sharp around the interface compared to the real ionic potential if we consider the infinitesimal size of ions. The finite size of ions results in the screening between ions that affect the potential. We adopt the model including the screening effect of IL which can be applied beyond the classical Gouy-Chapman-Stern model 54 . Note that the model is not applicable to the crowding situation of ions which occurs in the high gate bias voltage 55 . We explain the concrete description of the model in the Methods section.

Device simulation of WSe 2 transistors
The schematic view of ion-gated monolayer WSe 2 transistors considered is depicted in Fig. 1. The employed parameters for the calculation are written in Table 1. The length of WSe 2 in the x (channel) direction is set at 100 nm and the source and drain contacts with 10 nm length are attached to the top of WSe 2 monolayer. We assume that electrons tunnel between metal contacts and the channel at x = 10 nm and x = 90 nm so that the real channel length is 80 nm. We set the energy difference between the work function of TMD and metal contact to 0.6 eV considering the metal contact made of the alloy of Ti and Au 56 . In this calculation, we set the Fermi energy of the metal source contact to be zero. The energy level of the bottom of the conduction band is 0.6 eV at the source contact surface while the top of the valence band is −1.0 eV. The thickness of the IL is 10 nm and the metallic gate contact is attached to the top of the IL. The thickness of the SiO 2 layer at the bottom of the WSe 2 is 10 nm with the Neumann condition. We assume ions such as DEME-TFSI (N,N-diethyl-N-methyl-N-(2-methoxyethyl) ammonium bis (trifluoromethylsulfonyl)-imide) whose size is approximately 1 nm so that we fix the maximum ion concentration on the ion layer at the electric double layer to n max ¼ 1 10 21 cm −3 . The average ionic concentration is set to n i = 2 × 10 20 cm −3 . Though the employed n i is lower than the real concentration due to the restriction in the numerical calculation (abrupt change of the carrier and ion concentrations), the approximate maximum of the IL capacitance is 40 μF/cm 2 , which is large enough to bend the potential and make the Schottky barrier sufficiently thin so that the source and drain contacts show the Ohmic behavior. We use the constant mobility extracted from the experiment in ref. 15 . We neglect the V ds dependence of mobility since the potential is relatively flat inside the channel region beside the metal-TMD surface, as described later in the next section. Furthermore, we assume that the concentration of defects inside WSe 2 is negligible so that mobility is not affected by carrier concentration tuned by the gate voltage.
We start from examining the transport characteristics of iongated WSe 2 transistors using the theoretical model. As shown in Fig. 2a, b, the drain-source current, I ds , is calculated as a function of the voltage determined from the potential difference between the source contact and the center of the ionic liquid, Fig. 2b]. We should mention that V ref is adopted in Fig. 2a, b instead of the gate voltage between the source contact and the metal gate, V gs , to subtract the potential drop occurs at the surface of the gate contact for the considered system. Note that IL can be considered as two capacitors in series; one is the electric double layer capacitor (EDLC) of the gate electrode and the ionic liquid, and the other is the EDLC of WSe 2 and the ionic liquid. The surface area of the metal gate attached to the IL is usually much more substantial in experiments. Since the number of net ions at both capacitors should be the same, the ion concentration at the substantial metal gate is much smaller than the one at the WSe 2 surface. Thus, potential drops at the metal gate is negligible. However, in our calculation we set the surface area of the metal gate 100(nm)/80 (nm) larger than the one of the WSe 2 due to the constraints of the numerical calculation. Thus, approximately half of the potential drops at the metal gate. So we define V ref to compare the results with experiments. The bias voltage between the drain and source electrodes is set at V ds = 0.2 V. Note that the current flows from the drain contact to the source contact when V ds is positive. The center of the off-gate voltage slides to the minus (~−0.2 V) from the zero point since the energy difference between the metal contact and WSe 2 is Φ SB = 0.6 V which is 0.2 V above the center of the band gap. As shown in Fig. 2b, using the linear-scale plot, we estimate the threshold voltage from the intersection point of the line I ds = 0 and the dotted linear line of the current. Here, the dotted line is determined from the least-square method using the data from 1.2 V < V ref <1.5 V for the electron current, and −1.5 V < V ref < −1.3 V for the hole current. The determined threshold voltage is V th−e = 0.94 V for the electron current, and V th−h = −1.17 V for the hole current, respectively. The threshold voltages for electrons and holes have been used to measure the bandgap of the semiconductor as reported in ref. 38 . In our numerical calculation, the energy difference between the threshold voltages of electron and hole currents is e(V th−e − V th−h ) = 2.11 eV, which is 0.51 eV larger than the real band gap, E g . Figure 2c shows  Fig. 3. At V gs = −2 V, the band bending is extremely sharp at the contact, and both bottom of the conduction band and top of the valence band are shifted above the quasi-Fermi energy, φ n(p) , in the channel. Around the contact-channel junctions, the potential bends downward which induces the Schottky barrier. The IL makes the thickness of the Schottky barrier of the valence band at Fermi level about 0.9 nm for the holes to tunnel through. Note that the quasi-Fermi energy is plotted only for the main carrier since the minor carrier is almost zero due to the rapid reduction of the carrier inside the band gap for WSe 2 . At V gs = 0.2 V, the quasi-Fermi level is located inside the band gap. Thus, the tunnel of the carrier between the contact and channel is prohibited. Only the Schottky (thermionic) transport is possible, but the simulation shows that it is negligible at this gate voltage. At V gs = 2 V, the potential energy of the channel area becomes downward convex. Thus the quasi-Fermi energy is located at the conduction band which results in electron current. The carrier-type switching is possible with much smaller V gs compared to the case in the conventional transistors. The ions play a role of dopant (surface doping) and induce the ntype contact by making the Schottky barrier of the conduction band at the Fermi level about 0.4 nm so that carriers can tunnel through. The Ohmic behavior is confirmed by the I d -V ds characteristics discussed below. Here, we show that both n-type and p-type operation are possible depending on the direction of the Schottky barrier controlled by the metal gate attached to the IL even if the work function of the metal contact is located at the upper part of the band gap as we set in this simulation.
Next, we calculate the I d -V ds characteristics at V gs = 1.5 V where the carrier type is n-type. We plot the electron (dotted line), hole current (dot-dashed line), and the total current (solid line) as shown in Fig. 4. The I ds -V ds characteristics can be classified to the three regions: the linear (0 ⪅ V ds ⪅ 1.2 V), the saturation (1.2 V ⪅ V ds ⪅ 2.8 V), and the ambipolar region (V ds ⪆ 2.8V). In the linear region, the current increases with V ds approximately linearly around V ds = 0 V, followed by the saturation region starting at V ds 1.2 V. Above V ds~2 .8 V, I ds takes off again, indicating the onset of hole current injected from the drain electrode (ambipolar transport). This feature is similar to the experimental output curve of WSe 2 device for larger devices 13 . The increase of electron current in the ambipolar region is due to the recombination of electrons and holes inside the channel. Figure 5 shows the 2D plot of the net ion density of the IL (top), and band profile of the WSe 2 channel and the metal contacts (bottom) for the linear region (a: V ds = 0.5 V), saturation region (b: V ds = 2.5 V) and ambipolar region (c: V ds = 4 V). The figures from d to f are the enlarged views of the potential profile around the contacts. The open arrows denote the tunnel transport while the solid arrows denote the thermionic transport where carriers go beyond the barriers.
In the linear region (Fig. 5a), the cations are dominant at the interface attached to WSe 2 . Therefore, the current is dominated by electrons since the band gap allows only electrons to tunnel through the metal-WSe 2 junctions as depicted in Fig. 5d. By increasing the bias voltage V ds , the electrons tunneling from the channel to the drain contact gains. The inclination of the current in the linear region is estimated to be 7.5 × 10 −5 Ω −1 . Thus, the estimated contact resistance is R con~1 0 kΩ while the channel resistance is R ch = 1∕qμ n n~4 kΩ (n~4 × 10 13 cm −2 as discussed below) so that the contact resistance is comparable to the channel resistance, which results in the linear behavior.
In the saturation region, cations are dominant at the interface attached to WSe 2 while anions assemble at the drain contact as shown in the top panel of Fig. 5b Fig. 5b, e). The Fermi-level of the drain contact is located at the energy inside the band gap of the channel so that hole tunneling from the drain contact is limited by the Schottky junction as shown in Fig. 5e. The current is saturated even with increasing the V ds . This pinch-off behavior is also seen in conventional transistors.
Finally, in the ambipolar region, the net ion density changes gradually from the cation dominant to the anion dominant from the source to the drain at the interface attached to WSe 2 as shown in Fig. 5c. We find that anions assembled at the drain contact pull down the potential energy around the contact. By increasing V ds , the anions screens hole inside WSe 2 which makes the Schottky barrier thinner enough for holes to tunnel through. Furthermore, the Fermi energy of the drain contact goes below the top of the valence band in the channel at the metal-WSe 2 junction. Then, there are enough holes to tunnel through the Schottky barrier from the drain contact as depicted in Fig. 5f, which results in the p-n junction. In the p-n junction induced by IL, the drift and diffusion currents always move in the same direction. It means that the p-n junction exists only for the forward bias voltage higher than the built-in potential.
The ambipolar behavior does not occur for the conventional metal-oxide-semiconductor (MOS) transistors with the same V ds . In the Section III of the Supplementary information, we show the current and band profile of the WSe 2 MOS transistor, whose oxide layer has the same permittivity with IL ( Supplementary Fig. 2a). To fix the carrier concentration to the same order to the ion-gated transistors, we set the gate voltage to V gs = 10 V so that n~5 × 10 13 cm −2 at the equilibrium condition. The current increases monotonously with V ds and the transistor is always n-type as shown in Supplementary Fig. 2b. Supplementary Fig. 3 is the band profile for the MOS transistor. The potential energy decreases linearly at the channel region by applying V ds and the Schottky barrier at the drain area becomes perfect contact for electrons to pass through from WSe 2 to the metal. The Fermi energy of the metal contact is located inside the band gap so that the holes cannot pass through the barrier. Thus, the drain contact keeps acting as n-type contact and the current holds unipolar up to V ds = 4 V. Figure 6 displays the spatial distributions of the net carriers, electrons, and holes at the channel area for the linear region (a: V ds = 0.5 V), saturation region (b: V ds = 2.5 V), and ambipolar region (c: V ds = 4 V). In the linear region, the electron density slightly varies around n~4 × 10 13 cm −2 . In the saturation region, on the other hand, the electron densities are not constant anymore, and a depletion point appears at the drain contact. There are a few numbers of holes as shown in the inset of Fig. 6b (p~10 8 cm −2 ) which is injected from the drain contact by the thermionic transport although the hole tunneling is prohibited in the saturation region. In the ambipolar region, electrons and holes are spread throughout the whole channel area so that small difference between the electron and hole concentrations make the p-n junction. Therefore, the electrons and holes can recombine in the whole area of the channel as shown in the Supplementary Fig. 1 of the Supplementary information (Section II). There is no depletion region when the applied voltage is higher than the built-in potential. Then, the potential drop at the p-n junction is determined by the recombination zone. The size for the recombination zone is estimated to be 1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi β=D n max½pðxÞ p $ 1 μm in our case. Thus, the channel size we consider is much smaller than the recombination region. It leads to no potential change at the cross-section between n and p regions.
The application of the ion-gated TMD transistors can be very broad, stemming from light emitting devices 13,57 ,  Table 1 for other parameters.
A. Ueda et al.
superconductivity 58 , to thermoelectric properties 59 . Here, we pick up transistors that emits circularly polarized light owing to the peculiar electronic structure of TMDs. Since there has been no simulation model constructed in this area, the essential physical factor in determining the performance of the light-emitting device is yet not known. This simulation method is advantageous to explore the crucial factors and design the ideal device structure. We can examine the characteristics of light emission by calculating the recombination rate, as demonstrated in Section II of Supplementary Information. In the simulation, the parameters presented in Table 1 can be tuned to explore the crucial factors.
These factors can also be realized in experiments by changing the materials and device structures.
In conclusion, we have developed the DD model for the iongated TMD transistors and explained the characteristic transport behavior of the ion-gated WSe 2 transistors. We adopt the DD model coupled to the Poisson equation for calculating the band profile, carrier densities, and current. In the Poisson equation, the model including the screening effect of charge is employed to derive the potential profile of the IL. In the DD model, the generalized Einstein relation is considered for high carrier densities and the model of the practical Schottky contact is  employed for considering the Schottky and Ohmic behavior of the transistors. One of the peculiar features of ion-gated transistors of 2D materials is the ambipolar behavior with the application of the gate voltage comparable to the band gap energy. The simulation explains the fundamental physics underneath for ambipolar transport and the formation of p-n junctions in the channel. The present result indicates that the DD model coupled with the Poisson equation is a fascinating tool to study ion-gated transistors. By including the spin, valley, and optical degree of freedom to the DD model, we can examine the functionality and to design the device structure of ion-gated transistors.

METHODS
We introduce the concrete explanation of the models. In the DD equation, Eqs. (3) and (4), D n(p) (x) is not described using the ordinary Einstein relation since the Boltzmann distribution is not satisfied in the considered system due to the high carrier concentration above 10 13 cm −2 . Then, we adopt the generalized Einstein relation for the DD model 47 . The generalized Einstein relation which obeys the Fermi distribution is expressed as where Here, φ n(p) is the quasi Fermi potential of electrons (holes). k B is the Boltzmann constant, and T is the temperature. Note that g n(p) = 1 for the ordinary Einstein relation. For realizing the practical Schottky contact, we consider the thermionic current and tunnel current passing through the metal-TMD junctions. For the thermionic current density of electrons at the source (drain) contact is written as while the thermionic current density of holes is Here, F 1=2 is the Fermi integral where Γ(x) is the gamma function. A 2D is the Richardson constant for the 2D materials 60 , E c (x) and E v (x) denote the energy at the bottom of the conduction band and the top of the valence band, respectively. φ metal−s(d) is the Fermi energy of the source (drain) metal contact. Two is the factor stemmed from the valley degree of freedom of the monolayer TMD. We should mention that the thermionic current is considered only at the junction between the source (or drain) contact and channel, x s (x d ). The tunnel current densities at x from source (drain) contact for electrons and holes are written as Ecðx sðdÞ Þ dϵ F À1=2 À ϵ À φ n ðxÞ k B T À F À1=2 À ϵ À φ metalÀsðdÞ k B T ! T eÀsðdÞ ðϵÞ; respectively. Here, E(x) = −dΦ(x, 0)∕dx is the electric field in the x direction. Note that T e−s(d) (ϵ) and T h−s(d) (ϵ) are finite only when the arguments inside the exponential functions are negative. In the numerical calculation, the Fermi integral F 1=2 and F À1=2 has been obtained by using the approximation method based on the Sommerfeld expansion 61 The total recombination term is written as R(x) = R dir (x) + R SRH (x) + R n(p)−con (x). The ion distributions in the Poisson equation, Eq. (7), are expressed as i þ ðx; zÞ ¼ n s exp½À qΦ 0 ðx; zÞ kBT 1 þ 2γ sinh 2 qΦ 0 ðx; zÞ 2kBT ; (19) and i À ðx; zÞ ¼ n s exp½ qΦ 0 ðx; zÞ kBT 1 þ 2γ sinh 2 qΦ 0 ðx; zÞ 2kBT : (20) γ is the screening parameter of the ions stemmed from the finite size of ions. γ ¼ n i =n max where n max is the maximum concentration of both cation and anion which can pack in the liquid. Here we define shifted potential, Φ 0 ðx; zÞ, since Φ(x, z) is set to be zero reference to the Fermi energy of the source contact. Φ 0 ðx; zÞ and the coefficient of carrier concentration, n s , are determined to satisfy the condition Z Aion i þ ðx; zÞdxdz ¼ where A ion is the area of the ion. Note that the average carrier concentration, n i , is constant since there is no electrodes to supply ions in the system. In the numerical calculation, the Scharfetter-Gummel discretization is used for solving the DD equation which is described in the Section I of the Supplementary information.