Modeling the Transient Behavior of Ionic Diodes with the Nernst–Planck–Poisson Equations

Ionic circuitry formed from soft charged polymers promises to enable a new kind of analog computing and usher in a new age of human–computer interaction, but devices remain in their infancy. Ionic diodes, created by the interface between two oppositely charged polymers, form one of the most fundamental elements of ionic circuitry. However, it is currently not well understood how the boundary conditions and geometry of a diode influence its transient response and performance. A consistent set of metrics with which to analyze these diodes is first developed, then a continuum model based on the Nernst–Planck–Poisson equations is used to study how device construction and three types of boundary conditions affect performance. This model is solved using the finite volume method and gives insight into transient diode behavior not described by previous analyses. It is demonstrated how different parts of a diode's construction, thickness, and operating voltage act as the limiting factors for different regimes. To conclude, it is demonstrated how blocking electrodes limit both diode lifetime and rectification behavior and how neutral bath boundary conditions play an important role in the performance of some diodes, with recommendations for building higher performance devices in the future.


Introduction
Polyelectrolytes-polymers with charged backbones-are used to construct a wide range of ionic sensors and circuits.Using the rich library of organic polyelectrolytes available to the research community, researchers have been able to use additional functionalities of these materials to imbue ionic devices with multi-physics interactions.10][11] At the core of many of these devices are ionic diodes-interfaces between two differently charged materials that rectify currents. [1,3,11]Because DOI: 10.1002/adsr.202300131these junctions arise at the interface between any such polymers, they are also present in most ionic circuits formed from multiple elements, [12][13][14] as well as a variety of other ionic devices.23] As fundamental as these junctions are to many devices, the mechanisms that occur at these junctions are not completely well modeled.Researchers have modeled junctions in a variety of ways, including continuum [24] and moleculardynamics [25] based formulations.Most continuum models start with the Nernst-Planck-Poisson equations (NPP, briefly explained in Section 3).These equations are notoriously computationally expensive to solve in their full form, and most computational models to date have made simplifying assumptions. [26]The most common simplifying assumption is to solve only for the equilibrium state of the system. [24,25,27]However, most existing ionic devices are relatively slow to respond to ionic stimuli, with response times on the order of tens to hundreds of seconds. [12,13,21,28]This slow transient response is relevant because it limits the applications of ionic circuits to domains where a slow response is acceptable, and impedes the goal of interfacing these devices with biological signals, which often occur much faster.This slow timescale is ignored by existing models that focus on the equilibrium state.In simulating the NPP equations more generally, a number of papers have solved the transient time problem, [29][30][31][32] but as far as we are aware, these methods have not been applied to understanding the transient response of polyelectrolyte based ionic diodes specifically.Sabbagh et al. performed a transient simulation of an ionic diode, but did not share transient results for the diode itself, limiting their discussion to the equilibrium state. [23]Carneiro-Net et al. simulated the transient response of a micropore-based diode, with interesting conclusions, but did not discuss the internal mechanisms in detail; micropore diodes operate in both a different regime and with a different principle than bulk polyelectrolyte diodes. [33]A fundamental understanding of these dynamics that occur inside a polyelectrolyte diode is critical to identifying and shortening the response time of ionic circuits and sensors.
In this paper, we perform a time-dependent numerical analysis of the NPP equations as applied to a variety of ionic diode constructions.We use a 1D discretization based on the finitevolume method with a fully-implicit time-stepping scheme.We begin by introducing a set of four performance metrics in Section 2 that are closely tied to how ionic diodes are used in real applications.Next, in Section 3 we briefly describe the NPP governing equations and solution method, with more details presented in Experimental Section.In Section 4 we choose a baseline diode construction that demonstrates the main features possible in a diode transient response and then detail the key material and geometry features that influence each stage of this response.We also compare this model to experimental data from literature.Finally in Section 5 we examine the influence that microporous electrodes and electrode electrochemical reactions have on the transient diode response.

Performance Metrics
In order to discuss the performance of diodes, we must first define a set of performance metrics.Each of these metrics captures an independent aspect of how well a particular diode will perform when used in the context of a larger device.The metrics are intended to capture the deviation from ideal behavior.The primary goal of a diode is to rectify flows of species.An ideal diode allows arbitrarily large ionic currents to flow when a driving force is applied in one direction and zero current when a driving force is applied in the other direction.Deviation from this ideal behavior leads to the first two performance metrics: forward bias current capacity and current rectification ratio, denoted I f and  I respectively.An ideal diode also switches instantly between its conducting and non-conducting states.In contrast, physical diodes require finite time to switch between these states.In that time, some additional charge leaks back across the junction.This leads to the third and fourth performance metrics that we will use: the switching time and reverse excess charge, denoted t s and  q respectively.These metrics are explained in more detail in the list below and demonstrated graphically in Figure 1.
• I f is defined as the ionic current that flows through the diode when it is forward biased, normalized by the area of the junction interface.Like all the metrics, I f depends on the inherent physical properties of the diode as well as the measurement conditions.The former includes the material dielectric constant and charge concentration, while the latter includes driving voltage and measurement time.At long measurements times, I f approaches I ∞ f , the steady state forward bias current.•  i is the ratio of final forward and reverse bias currents when the diode is subjected to a square wave driving voltage.This driving voltage is centered around the voltage for which no current flows in the diode."Final" means the currents at the end of each period of forward and reverse bias.This metric was chosen to match the way that ionic diodes are often characterized experimentally.Like I f , it is a function of both inherent diode properties and measurement details.If the period of the square wave is very large compared to the relevant diode timescales, then this metric approximates the steady state rectification ratio.Top: Applied square wave voltage with a peak-peak amplitude of V peak-peak and a half cycle switching period of T switch .Bottom: Diode current over time.I f is steady state forward bias current.I r is the steady state reverse bias current. I is the current rectification ratio.t s,r is the reverse bias settling time. q is the excess charge metric that quantifies how much charge returns in reverse bias, which is defined in terms of the ratio of the lightblue and dark-blue-striped areas.The red and blue areas correspond to charge.
• t s,r is the time required for the current to settle within 5% of the steady state reverse bias current after a step change in applied voltage from forward to reverse bias.•  q is the excess charge passed before t s,r , divided by the charge that would have flowed in steady state for the same period of time.Calculating  q does require either the voltage switching period to be long or an independent measurement of the steady state reverse bias current.

Modeling Approach
In order to explore the effect of carrier properties, boundary conditions and bulk material properties on the performance metrics, we implemented a continuum model based on the standard NPP governing equations for ion transport.As the development of these equations is discussed extensively in previous work, we will only briefly discuss them. [34]We model all of the diodes in this paper using a 1D framework for which the diode is assumed to be uniform in the xy plane and vary with z coordinate.The NPP equations model the combined diffusion and electromigration of charged species.In order to efficiently model the electric fields at the timescales we are interested in (μs to s), we assume that the speed of light is instantaneous.This simplifies Maxwell's equations to their electrostatic form, which gives us the Poisson part of the NPP equations [Equation (1)].In order to derive the Nernst-Planck part [Equation (2)], we assume a form for the free energy of our ions that combines ideal solution mixing with an electrical potential energy.Finally, we assume that the flux of species is linearly proportional to the gradient of the electrochemical potential, which ensures a minimization of the free energy of the system. [34] wherec (i) and z (i) are the concentration and valence of the ith mobile species, ϕ is the electric potential, ϵ is the material dielectric constant, F is Faraday's constant, R is the universal gas constant, and T is the material temperature.We model three kinds of boundary conditions: capacitive blocking electrodes, infinite concentration baths, and Butler-Volmer type surface chemical reactions.The details of boundary condition implementation are presented in Experimental Section.
In order to solve these coupled non-linear PDEs, we discretize space using the finite volume method and discretize time using a second order backward differencing formulation (BDF2). [35]This scheme is implemented using FiPy, an open source finite volume python library. [36]We use the fully implicit forms of the discretization in order to take larger timesteps than an explicit scheme would allow.We additionally use a modified Powell's method provided by MINPACK through SciPy instead of the fixed point iterations that FiPy uses by default. [37,38]The specifics of this scheme are discussed in Experimental Section.A table with representative mesh point counts, timesteps counts, and simulation time is given in Table S1, Supporting Information.

Baseline Ionic Diode and Variations
In order to provide a baseline against which to compare the ionic diodes variants, we will first outline the construction and operating principles of a diode connected to salt baths, then detail the sensitivity of this diode to its construction and operating conditions using the continuum model.We choose this setup as our baseline because it avoids limits imposed by other diode constructions.The mechanisms and explanations introduced will largely apply to the subsequent sections as well, where we will explore additional boundary conditions and diode constructions.
The baseline diode (Figure 2) is formed from four regions.From one side to the other these are: a neutral salt bath, a polycation with a non-mobile positive charge background, a polyanion with a non-mobile negative charge background, and finally a second neutral salt bath.There are two mobile ions, an anion and a cation, both monovalent and with opposite charge to each other.We set z = 0 at the junction between the two polyelectrolytes.The left and right boundaries have constant concentration boundary conditions enforced for the two mobile ions.The function of this baseline diode is best understood as three diodes/interfaces in series.In addition to the rectifying interface formed between the two polyelectrolytes, two more rectifying interfaces are formed between the polyelectrolytes and the neutral baths.We will refer to the bath-polyelectrolyte junctions as the outer junctions and the polyelectrolyte-polyelectrolyte junction as the central junction.The equilibrium voltage and ion concentrations consist of constant-value regions in the bulk of each material separated by steep changes in value at each of these three junctions (Figure 2).These equilibrium features match those found in prior work. [24,25,39]The sharp transitions in concentration and voltage occur over the Debye length.The flat regions are nearly exactly electro-neutral, with concentrations that obey Donnan exclusion. [40]o explain how diodes are able to rectify currents and detail their transient response, we analyze the difference between what occurs in the diode in forward and reverse bias.In this discussion we will use the terms "majority carrier" and "minority carrier."A majority carrier is an ion in a region with background charge (the polyelectrolyte regions) where the ion has opposite charge to that background charge.For example, chlorine ions in a quaternary-ammonium based polyelectrolyte would be majority carriers.Minority carriers are ions that have the same charge as the background charge in a given region.At equilibrium, with no applied voltage, the concentrations of mobile ions within each region are constant and no current flows anywhere in the diode.In order to satisfy the Donnan equilibrium between the neutral bath regions and each polyelectrolyte region, the ratio of minority carrier concentration within each electrolyte to the bath concentration must be the same as the ratio of bath concentration to majority carrier concentration. [40]Since the concentration of majority carriers is larger than the concentration of minority carriers in a given region, the response of the majority carriers to an electric field determines/dominates the concentration profile inside a polyelectrolyte at lengths larger than the Debye length.Majority carriers experience a larger electromotive flux than minority carriers for a given and shared electric field, so they will tend to accumulate at interfaces where they become minority carriers (and thus less driven).Likewise, minority carriers deplete at places where they transition to majority carriers (and are thus driven more).The currents that flow through each section are controlled by the minority carrier currents.The concentrations of minority carriers are much smaller than that of the majority carriers.This implies that the conductivity of these carriers is smaller.As every ion passes through the diode it must be a minority carrier in one region.In dynamic equilibrium, the current for each species is uniform across the device (implied by the lack of continued ion accumulation).Together, these two facts imply that the currents are limited by the regions of lowest conductivity for each species, which typically is where the ions are minority carriers.
The set of constants for our baseline simulation is listed in Table 1.The values for these constants are chosen based on two criteria.First, values are chosen to be within the range of constants for experimental systems.Second, the values have been chosen such that each of the transient response and spatial features are visible on the same figure.To satisfy this second constraint, L ionomer , L bath , C ionomer , and C bath are on the smaller side for experimental systems, while ϵ is on the larger side. [13,41,42]R and F are universal constants.

Response to Applied Voltage Square Wave
When a voltage is applied across the diode, the currents that pass through the diode depend on the polarity of the voltage applied.By looking at how this current changes with time and the internal concentration profiles inside the diode, we can begin to understand the transient response of the diode (Figure 3).The transients shown in this paper are taken from the second voltage square wave cycle to represent a typical cycle.
When a positive voltage is applied to the side of the diode with a polyanion, the central junction is placed into forward bias while the bath-polyelectrolyte junctions are placed in reverse bias.Since the central junction is a much better rectifier than the bath junctions if the polyelectrolyte concentrations are large (because of the much larger difference in background charge), we denote this as forward bias for the overall device.This for-ward bias proceeds in two stages.In the first stage of forward bias, the majority carriers in each polyelectrolyte are drawn toward the central junction, slightly depleting the outer junctions and slightly enriching the central junction as the imposed electric field is minimized by charge rearrangement (Figure 3b, red).In the second stage of forward bias, the minority carriers that crossed the central junction begin to diffuse across the polyelectrolyte regions until they reach the bath-polyelectrolyte interfaces.As they diffuse, these minority carriers increase the conductivity of the polyelectrolytes.Once the bath-polyelectrolyte junction is reached, these carriers re-enrich the outer junctions with carriers (Figure 3b, gray) and raise the overall current flowing through the diode.
In reverse bias, that is, when a negative voltage is applied to the polyanion side, the opposite transport process occurs.Majority carriers flee the central junction and accumulate at the outer junctions, placing the central junction into reverse bias and the outer junctions into forward bias (Figure 3c).Minority carriers slowly cross the central junction and become majority carriers, while the opposite occurs at the outer junctions.If the reverse bias is applied after a period of forward bias, then this process is again broken into two separate stages.First, the concentration of minority carriers is quickly pulled to near zero close to the central junction and increases at the outer junctions (Figure 3c, gray).Then minority carriers that enriched each polyelectrolyte in forward bias slowly diffuse back toward the central junction (Figure 3b, blue).Note that this diffusion is not limited to minority carriers near the central junction.Minority carriers throughout the whole polyelectrolyte diffuse back and contribute to this second stage.As the enrichment in each polyelectrolyte decreases, the current carried by these minority carriers eventually falls, eventually reaching a reverse bias steady state current smaller than the forward bias steady state current.If there is no preceding forward bias, the first stage proceeds as before, but, since there are no excess minority carriers to extract, there is no second stage.The lack of this second extraction stage means that the diode reaches equilibrium much faster.

Comparison with Experiments
In order to validate this model, we compared our simulation with an experimentally realized diode with a similar structure, that of Sabbagh et al. [13] This diode is constructed from poly(diallyldimethylammonium chloride) (pDAD-MAC) and poly(2-acrylamido-2-methyl-1-propanesulfonic acid) (pAMPSA) hydrogels on a planar microfluidic chip, with counter ions K + and Cl − respectively.These hydrogels have a very high ionic content since they are polymerized from ≈5M monomer solutions, with every monomer contributing one counter-ion.The diode is formed from two trapezoidal-shaped pieces of polyelectrolyte, with the narrow edges of the trapezoid in contact with each other and the wide edges of the trapezoid in contact with a neutral salt bath (Figure 4a, upper).The junction area is ≈230 μm by 25 μm, while the bath junctions are ≈400 μm by 25 μm.The width of each polyelectrolyte from inner junction to outer junction is ≈150 μm.
The parameters we used to simulate this diode were either derived from information provided in the paper, or hand fit to the shown response while staying within reasonable limits.A couple of simplifying assumptions had to be made to represent the diode effectively using our model.First, since our model is 1D, we replace the trapezoidal polyelectrolytes with a uniformly wide polyelectrolyte at the average trapezoid width of 315 μm (Figure 4a, lower).The baths in the physical device are also much wider than the polyelectrolytes, meaning that their conductivity is higher than a 1d model would predict.To compensate for this decreased width, we use a short bath thickness of only 8 μm, which was chosen to match the experimental data.The only other free parameter was the diffusivity, D ionomer which was set to the same value for both polyelectrolytes, and matched the data best with a value of 1× 10 −10 m 2 s −1 .The full set of diode parameters used in the simulation is presented in Table S2, Supporting Information.In order to improve numerical convergence, we scaled all concentrations down 1000 times lower than in the experimental diode, and then scaled the currents up by 1000 times.This scaling has minimal effect on the results since the currents scale linearly with concentration (Figure S1, Supporting Information).Our model is able to capture many of the key features of the transient response of this diode (Figure 4b).Specifically, the simulation matches the magnitude of response and the reverse bias timescales.It also predicts the gradual increase in forward bias current over time, the origin of which is explained in the previous section.
The main feature that our simulation does not capture well is the fast forward bias timescale.We predict a much faster fall in current than is present in the experimental data.It is likely that at least some of this difference is explained by our model being 1D.The drop in current is controlled by the outer interface, which, in the experimental device is much larger than the junction interface width.This size difference means that there are more ions to deplete at the outer junction, which might increase the time required to finish this depletion.We also predict a much smaller reverse bias current, and thus a much larger  I than present in experiments (5900 vs. 52).Part of this current disparity might be explained by the current leakage mentioned in the paper. [13]oth the disparity in current and the disparity in timescale might additionally be explained by the assumption of our model that the activity of each ion is equal to its concentration.If the activity coefficient is instead much lower than one, this will reduce the Donnan exclusion between the baths and the polyelectrolyte regions, increasing the minority carrier concentrations.This increase in concentration will increase the time required to deplete the junctions and thus the fast forward bias timescale.This concentration increase will additionally reduce  I because it will increase the reverse bias current.

Parameter Dependence
The lengths of the different regions that comprise a diode strongly influence diode performance.In Section 4 we discussed  [13] (upper) and our simulated geometry (lower).We simplify the 3D experimental geometry (effectively 2D) to a uniform 1D system in order to simulate it using our numerical code.b) Diode current versus time for the experimental diode (dotted) and our simulation (blue).The simulated reverse and forward bias steps shown are taken in sequence after one initial forward bias step.The data from Sabbahh et al. are taken in an unknown order and with unknown prior voltage history, and are presented in the order shown for comparison with our simulation.how in the second stage of both reverse and forward bias, ions must fully diffuse across the polyelectrolyte layers before equilibrium is reached.This diffusion process is controlled by the width of the polyelectrolyte layers.By comparing this diffusion time to the desired switching period of the diode we can gain insight into when this diffusion process controls performance.The natural non-dimensional number relating the switching time and polyelectrolyte thickness is given in Equation ( 3), where T switch is the voltage switching period.
When D is much larger than 1, the minority carriers waves diffuse all the way to the outer junctions during one forward bias step.When D is much smaller than 1, the wave fails to reach the other junction and the diode current never reaches equilibrium within the switching period.
The chosen switching time and polyelectrolyte width have a strong effect on the current and charge rectification metrics we defined in Section 2 (Figure 5).The performance plots are divided into two regions: One where D is small and the other when D is large, separated by D ≈ 0.25.We discuss each of these regions in turn.When the switching time is large or the polyelectrolyte width is small then diffusion is fast compared to the device geometry ( D > 0.25).In this regime, increasing the polyelectrolyte width greatly increases  I .This increase mainly comes from a large drop in the steady state reverse bias current, as evidenced by that fact that I f also drops.At the same time, increas-ing L ionomer also increases  q and t s,r .Both these metrics increase because more minority carriers are stored in each polyelectrolyte during forward bias.Most of these charges then return in reverse bias, increasing the time it takes for the diode to settle.When the switching time is fast or the polyelectrolyte width is large then diffusion is slow compared to the device geometry ( D < 0.25).In this case we observe a different effect: minority carriers from the central junction do not reach the outer junctions during the switching period.As L ionomer is increased, both  I and  q fall.The equilibrium diode state is never reached, and the forward bias current never has the opportunity to reach its steady state value.This can also be observed on the t s,r plot, where the settling times plateau at a maximum of ≈9 s.This plateau occurs because the true settling time is now longer than the switching period of the diode.Notably, most polyelectrolytes have a range of diffusivity from 1×10 −10 to 1×10 −12 m 2 s −1 , which combined with a typical switching frequency of 10s, gives an acceptable polyelectrolyte thickness of 10 to 100 μm before D falls under 0.1.To switch faster requires a smaller thickness.
This switching time and polyelectrolyte width versus rectification ratio trade-off influence diode design decisions.If the goal for a diode is to have the highest rectification ratio, then thick polyelectrolyte layers are desirable, but come at the cost of a slow switching time.If, instead, the goal is to create a diode that switches quickly or has a high current density, then thin polyelectrolyte regions are required.
The driving force, or applied electrochemical potential, also influences diode performance.In our model, because the baths are at equal concentrations, the driving force across the diode comes entirely from the applied voltage.As the applied voltage increases, the forward bias current increases dramatically and the reverse bias current stays roughly the same.The forward bias current increases largely due to the increase in equilibrium minority carrier concentration everywhere, and especially near the outer junction (Figure 6a).This concentration increase leads to a nearly linear increase in rectification ratio with applied voltage (Figure 6b, top).This concentration increase also sub-linearly in-creases the reverse bias settling time of the diode (t s,r ) as more ions must be extracted (Figure 6b, bottom).
As discussed previously, the concentration of minority carriers in the polyelectrolyte regions depends on the Donnan exclusion between those regions and the neutral baths.The higher the bath concentration (c bath ), the more minority carriers occupy the diode at equilibrium and the larger current that passes in both forward and reverse bias.For our baseline diode, we use  the ratio between equilibrium minority and majority carrier concentrations to quantify the diode's minority carrier content.When this ratio is small, and for monovalent mobile ions, it can be expressed in terms of the bath and polyelectrolyte background charge concentrations as in Equation ( 4).
As the bath concentration increases, this ratio increases, and both forward and reverse bias currents increase by similar amounts as each other, which hurts both  I and  q (Figure 7, left).As c minority c majority approaches one, the diode ceases to function as a diode, with the current rectification ratio approaching 1.One might expect that the reverse bias settling time, t s,r , would greatly increase with increasing bath concentration, because many more ions must be extracted from the junction in order for it to turn off.This does happen, but is partially compensated by the higher currents that are able to flow in a diode with more ions inside (Figure 7, lower right).Together, these two effects somewhat cancel out and the settling time is not strongly unaffected by the bath concentration (Figure 7, upper right).
While decreasing the bath concentration is effective at increasing diode rectification performance, it has another consequence: it decreases the conductivity of the baths that connect the diode to the rest of the circuit.For our simulations, the width of neutral bath is the same order of magnitude as the polyelectrolyte width.As the bath concentration decreases, more voltage must drop over the bath to flow the same current.This means that in order to achieve the same voltage drop across the polyelectrolytes, larger voltages must be applied.This problem is most relevant when the conductance of the baths is equal or smaller to that of the polyelectrolytes.If the baths are much more conductive than the polyelectrolytes, then the majority of the voltage applied to the system will drop over the polyelectrolytes.

External Interface Performance Limits
The performance of an ionic diode can depend strongly on what its boundaries are connected to, and how those boundaries are controlled.Up to this point, we have simulated only concentration boundary conditions, which place the least constraints on diode performance.In this section, we explore the effects of two other types boundary conditions: blocking electrodes and electrochemical reactions.Blocking electrodes are electrodes that do not allow charge to pass through them.Instead, they provide a capacitive coupling between the electrode and the diode material.Redox electrodes use a chemical redox reaction to generate and consume mobile species at the electrodes.

Blocking Electrodes
Blocking electrodes are a simple way to connect a diode to a larger electrical circuit.They have no chemical reactions, and rely on a capacitive interaction between an electrical conductor and the solution ions to pass current.Despite their simplicity, their use imposes two main constraints on the function of a diode.
First, for a given voltage, the electrodes can store a finite quantity of charge (which is determined by their areal capacitance, denoted as C here in units of C m −2 ).This finite capacitance limits the total charge that the diode can pass in a given bias state, which, when combined with the current that flows across the diode during that bias state, yields two additional system time constants,  c,f and  c,r .These time constants are inversely proportional to the current that flows and proportional to the electrode capacitance.Since the current in reverse bias is smaller than in forward bias by a factor of  I , to first order, the reverse bias capacitance timescale,  c,r , will be  I times larger than the forward bias capacitance timescale,  c,f .Thus,  c,r ≈  c, f  I .This relationship implies that in order to continue functioning in an ionic circuit long term, a diode with capacitive boundaries must be placed in reverse bias for  I times longer than it is placed in forward bias.The better a diode is at rectifying current, the less fractional time it can be used in forward bias.Assuming that a diode is run in an overall balanced fashion, the finite capacitance can still affect performance in a single diode cycle if the capacitive timescale is short enough compared to the switching time.As the capacitive boundaries charge up, the voltage that drops across those boundaries increases.The larger this voltage drop, the less voltage that drops across the diode itself.In order for this voltage to remain less than one tenth of the applied voltage over one forward bias interval,  c,r must be greater than ten times the switching time.With existing hydrogel-based diodes, forward bias current densities of 200 A m −2 [13,28] are reasonable.In order for  c,f to be at least ten times a switching time of 10 s, C must be greater than approximately This capacitance is far larger than any existing material.For comparison, microporous carbon electrodes might more reasonably have a capacitance of 8 F m −2 . [43]Alternatively, using a microporous material, the electrode area could be 2500 times larger than the junction area.This might be reasonable for extremely small junctions, but poses an almost insurmountable limitation on hydrogel diodes with blocking electrodes.For other polyelectrolyte materials with lower current densities, the situation is more favorable.For example, with a current density of 2 mA m −2 and a V peak-peak of 700 mV, the water-free ionomer diode sensor by Kim et al. requires an areal capacitance of only ≈0.3 F m −2 , which is eminently achievable. [41]econd, blocking electrodes can be used either in combination with a neutral salt bath or directly in contact with the polyelectrolytes.Existing devices that use blocking electrodes tend to use them directly in contact with the polyelectrolytes, which eliminates the function of the outer junctions in keeping ions that have already crossed the entire diode from coming back and keeping the concentration of minority carriers low.The removal of the salt baths simplifies device construction by removing two layers, but significantly increases the reverse bias current and decreases the forward bias current, decreasing the performance of the diode (Figure 8).The decrease in forward bias current compared to the diode with baths occurs because the microporous region has a much lower ion concentration and thus conductivity than the baths do.

Redox Electrodes
Finally, we implemented electrochemical boundary conditions that represent the Ag-AgCl electrodes that are commonly used in studying diodes.These electrodes are generally used in combination with a neutral salt bath, [13,14] and as such, their performance is comparable to that case.The case where electrodes are directly connected to the polyelectrolytes was previously studied analytically by Yamamato and Doi [24] and experimentally by Nyamayaro et al. [44] When the salt baths are large compared to the polyelectrolyte length and the electrochemical reaction is fast, the performance of each diode is indistinguishable from that of a concentration boundary (Figure S4, Supporting Information).When the chemical reaction is slower, more of the applied voltage must drop over the electrode interfaces in order to supply the same current (Figure S5, Supporting Information).If the reaction is slow enough, this becomes the limiting process and hurts both apparent rectification and current density performance.Thus a reaction with fast kinetics like Ag-AgCl is ideal.Previously, we showed that for concentration Dirichlet boundary conditions, there is a linear increase in steady state current with applied voltage.With chemical reaction boundary conditions, we also observe a linear increase in current with voltage (Figure S6, Supporting Information).Our formulation uses the Frumkin-Butler-Volmer (FBV) formulation, which explicitly incorporates the diffuse ion layer and relates the reaction rate to the electric field at the electrode surface. [35,45]Because our concentrations are relatively small, our Debye length is large compared to ion size, which places us toward the Gouy-Chapman limit of the FBV equations.Since the voltages we are applying are also small, our current will scale approximately linearly with voltage.This stands in contrast with the exponential current scaling reported by Yamamato and Doi for chemical reaction boundary conditions.Their model operates with large voltages and in the Helmholtz limit, which leads to an exponential current scaling. [45]sing chemical reaction boundary conditions also places constraints on diode lifetime, mainly in terms of the quantity of reactant available.For example, with Ag-AgCl electrodes, there is only a finite amount of oxidized silver present on each electrode.This layer ranges from nm to μm, depending on the electrode.A typical electrode preparation procedure [46] of 10 A m −2 for a few minutes forms an AgCl layer with a thickness of ≈1 μm, which stores ≈3.7× 10 3 C m −2 of charge. [47]Running at this same current density in a diode results the same lifetime of a few min-Like with the capacitive electrodes, this run time can be improved by either increasing the electrode area compared to the junction area or decreasing the current that flows through the diode.The electrode area can either be increased geometrically, or by using a porous structure similar to that of a microporous capacitive electrode. [48,49]

Conclusions
In this paper we have introduced a set of metrics that allowed us to discuss the transient response of a polyelectrolyte diode under a variety of parametric regimes and boundary conditions.Existing literature mainly discusses the rectification ratio, which ignores the transient response that is fundamental to the real world application of these diodes.We propose t s,r and  q to serve as a set of standard metrics by which the transient response of future experimental and simulated diodes can be compared.Additionally, we discussed in detail for the first time how the observed current response of an ionic diode is created from the internal movement of ions.Previous literature has focused on the role of boundary conditions and the central junction in determining the diode behavior, whereas we showed that the outer junctions play an important role as well in controlling both the short and intermediate timescale transient response for many physically relevant parts of the parameter space.We then showed how device construction could be effectively summarized through a non-dimensional diffusivity, D, the minority carrier content, and the driving force, V peak-peak , that determine the nature of the diode transient.
Our work has a direct impact on the construction and optimization of ionic diodes in many applications.For example, if trying to quickly and accurately control the release of a neurotransmitter or implement digital logic, the most important metric might be switching time.Prior work discussing switching time has recommended sharpening the inner junction. [42]Our work provides a basis for further recommendations.To optimize for that case, a diode should have as thin polyelectrolyte regions as possible, use a small driving voltage, and use small bath ion concentrations (listed in order of importance).If the diode is instead used for low frequency energy harvesting or osmotic energy harvesting, where the rectification ratio may be more important than the switching time, then a small bath concentration, higher voltage, and larger polyelectrolyte region should be used (also in order of importance).
One known limitation of the model is that in using the Nernst-Planck equations, we implicitly assume that the activities of the polymers are equal to their concentrations.This is a standard assumption in related literature, [24,50] but is known to be inaccurate for the concentrated polyelectrolytes present in most ionic devices. [51,52]The difference between concentration and activity will influence the potentials that drop across each of the junctions and the strength of the Donnan exclusion of minority carriers.This will not significantly change the way our performance metrics vary with device construction, but will change the absolute values of those metrics and what behaviors are possible in a given parameter range.
Future work can also likely improve on the experimental agreement by extending the model to three dimensions.This extension does not require any additional theory to be developed, but does pose a significant challenge in terms of computational complexity.We solved the full form of the NPP equations in 1D, and doing so for our systems under study was at the limit of easily available computational resources.Rewriting the performance critical parts of our code in a lower level language than Python would likely bring some more complex systems within computational reach, but a change in approach is needed to simulate large scale ionic circuits.Of specific concern to improving ionic sensors is extending this model to ionic transistors formed from multiple junctions, which will be critical to amplifying small ionic sensor signals.As we have shown, the transient responses of ionic circuits will continue to play a large role in their performance, and thus computational and design tools must be developed that reduce the cost of analyzing transients for larger systems.

Experimental Section
Voltage Boundary: A pseudo square wave Dirichlet voltage boundary V square (t) was applied according to ϕ(L, t) = V square (t) as shown in Equation (5).In this equation, A is the amplitude of the wave,  is the period of the wave, and k is a sharpness parameter that determines how fast the square wave switches sign.For the simulations, k was chosen to be 0.01, which resulted in a voltage transition time of 0.1 s.A pseudo-square wave was used instead of a sharp square wave to improve numerical stability around the transition.
Capacitive Boundary: In physical diodes, high electrode capacitances are achieved by using microporous electrodes.This microporous effect is fundamentally a higher dimensional phenomenon.In order to model it to first order, the effect was coarse-grained by creating a region at the diode boundaries with an artificially large dielectric constant.This allowed to model a variety of electrode capacitances using the 1D simulation.From Guy-Chapman theory, the capacitance of a double layer per area for small voltage-drops is   d where ϵ is the dielectric constant and  d is the Debye length.The increased dielectric constant captured the increased polarizability of the region due to the proximity of the volume to the electrode surface.For the simulations shown, the relative dielectric constant was increased to 1×10 6 , which is 1.28×10 4 times that of water.This corresponded to a capacitance of ≈8.2 F m −2 , which was on par with some commercially available microporous electrodes. [43]lectrochemical Reaction Boundary: Butler-Volmer kinetics were included by including additional source terms to Equation (2) as detailed in Equation ( 6), where  is a characteristic length for the charge transfer reaction ( ≈ 1× 10 −10 m),  is the charge transfer number for the reaction (with positive numbers transferring electrons into the electrode), n is the surface normal vector, and j 0 and c (i) 0 are the equilibrium exchange current and concentration respectively for involved species i. Equation ( 6) is based on the Frumkin-Butler-Volmer formulation. [45] This differs from the standard Butler-Volmer form because the space charge boundary layers of the system were explicitly simulated.The electrochemical reaction proceeded according to the local change in free energy associated with an electron hopping between the electrode surface and an ion right at the inner Helmholtz plane.Simulation Stages: Every simulation was calculated in two time stages.First, the system was equilibrated for 500 s with ϕ(0, t) = 0 and d dz (L, t) = 0.In the second stage of the simulation, the pseudo-square wave voltage boundary [Equation ( 5)] was applied with first forward and then reverse bias.Each of the rectification metrics was calculated using the second voltage cycle.
Time Discretization: A second-order variable-step backwards-finitedifference time-stepping (VSBDF2) was used in the simulations. [53]This time stepping scheme is described by Equation (7), where u is a degree of freedom, and  = Δt now /Δt old is a timestep factor, Δt now , Δt old are the current and previous timesteps respectively, and g(u n + 1 ) is an implicit formulation of the right hand side of the Nernst-Planck equation [Equation ( 2 Mesh: The solution was obtained over a highly non-uniform mesh which was chosen for computational efficiency.This non-uniformity was efficient because there are large spatial changes in the degrees of freedom (concentrations c (i) and electric potential ϕ) near the material interfaces.As such, the meshwas highly refined near these interfaces.The mesh spacing was computed by solving the differential equation in Equation (8), where (x) is a mesh density function.(x) is chosen to be the maximum of a baseline mesh density and a set of Gaussian peaks centered at each interface and with standard deviation a few times the Debye length (Figure S7, Supporting Information).Then, a mesh face was chosen to be placed at every z(n) where n is a positive integer.dz() d = 1 (x) (8)   Current Computation: The displacement and drift currents were summed and averaged over the simulation domain to produce the currents shown in each every figure.For the bath concentration parameter sweep, the median of the current was taken instead of the average to eliminate minor current outliers that arise from discretization errors at the junctions.The currents were computed after running the implicit solver for each time-step using an explicit formulation of the NPP equations.
Preprocessing: Data was extracted from Figure S8A, Supporting Information in the paper by Sabbagh et al. [13] using WebPlotDigitizer (https:// automeris.io/WebPlotDigitizer).The experimental current data were normalized by the average area of the device.
Statistical Methods: No statistical analyses were performed for this paper.

Figure 1 .
Figure 1.Diagram of the a typical diode response with annotated metrics.Top: Applied square wave voltage with a peak-peak amplitude of V peak-peak and a half cycle switching period of T switch .Bottom: Diode current over time.I f is steady state forward bias current.I r is the steady state reverse bias current. I is the current rectification ratio.t s,r is the reverse bias settling time. q is the excess charge metric that quantifies how much charge returns in reverse bias, which is defined in terms of the ratio of the lightblue and dark-blue-striped areas.The red and blue areas correspond to charge.

Figure 2 .
Figure 2. Equilibrium voltage and mobile ion concentration versus space with diode constituent regions annotated.From left to right these regions are a neutral salt bath, a polycation with mobile anions, a polyanion with mobile cations, and a second neutral salt bath.

Figure 3 .
Figure 3. Transient response for the baseline diode in forward and reverse bias when subject to a 200 mV square wave at 0.1 Hz. t = 0 represents the start of the second voltage cycle.a) Time series plot of the electrical current flowing through the diode over time.The cycle shown is the second cycle of the applied square wave.b) Voltage and anion concentration versus z coordinate during forward bias.The colors correspond to the colored points in (A) from 0 to 10 s.A wave of minority carrier anions travels from the central junction to the outer junctions.The outer junction is initially depleted of ions and then re-enriched when the wave arrives.c) Voltage and anion concentration versus z coordinate during reverse bias.The colors correspond to the points in (A) from 10 to 20 s.Minority anions that crossed the junction in forward bias are re-extracted from the poly-anionic region across the central junction.The outer junction is initially enriched, then returns to normal.

Figure 4 .
Figure 4. a) Diagram of the construction of a diode from Sabbagh et al.[13] (upper) and our simulated geometry (lower).We simplify the 3D experimental geometry (effectively 2D) to a uniform 1D system in order to simulate it using our numerical code.b) Diode current versus time for the experimental diode (dotted) and our simulation (blue).The simulated reverse and forward bias steps shown are taken in sequence after one initial forward bias step.The data from Sabbahh et al. are taken in an unknown order and with unknown prior voltage history, and are presented in the order shown for comparison with our simulation.

Figure 5 .
Figure 5.The four performance metrics versus the width of the polyelectrolyte (ionomer) regions in the diode.The blue, orange, and green lines respectively represent increasing diffusivity of the ions in the entire device.The black stars indicate the place along each respective curve where the nondimensional parameter D is 1.Larger L ionomer corresponds with smaller D. Upper left: Current rectification ratio,  I .Upper right: Reverse bias settling time t s,r .Note that here we use the current at the end of the square wave period to approximate the steady state current, which is why the settling time reaches a plateau of the square wave period.Lower left: Reverse bias excess charge,  I .Lower right: Forward bias current at the end of a forward bias cycle, I f .

Figure 6 .
Figure 6.a) Forward bias equilibrium anion concentration in the polyelectrolyte regions versus z coordinate.The colors correspond to voltages and the color of points on the right plots.The higher voltages greatly increase the concentration of minority and even majority carriers in their respective regions.b) Current rectification (top) and reverse bias settling time (bottom) versus peak-peak applied boundary voltage. I grows nearly linearly with voltage, while t s,r grows sub-linearly with voltage.The metrics not shown here are shown in Figure S2, Supporting Information.

Figure 7 .
Figure 7.The four performance metrics versus the ion concentration of the bath, c bath .Upper left: Current rectification ratio,  I . I upper right: reverse bias settling time t s,r .Lower left: reverse bias excess charge,  I .Lower right: forward bias current at the end of a forward bias cycle, I f .

Figure 8 .
Figure 8. a) Anion concentration in the polyanionic region of the a diode for both high capacitance and concentration (baseline) boundary conditions.The lack of outer junction for the capacitive case increases the minority concentrations near the outer boundary.b) Current rectification ratio (top) and forward bias current (bottom) versus the thickness of the polyelectrolyte layer.Removing the outer junction greatly decreases both metrics.