Impact Statement
Permselective nanoporous materials are ubiquitous in desalination, energy harvesting and bio-sensing systems. Of particular importance are bipolar membranes and nanochannels that comprise two oppositely charged permselective regions. While a plethora of experimental works have characterized the electrical response of these systems, a fundamental understanding of the underlying physics determining the response is still missing. To address this knowledge gap, we have systematically simulated different bipolar nanofluidics systems subject to varying potential drops and characterized their electrical response to reveal signatures that are unique to every system. Our findings contribute to a more profound understanding of the various control parameters and mechanisms that determine the time-transient dynamics and the steady-state current–voltage response in bipolar systems and provide a valuable tool for interpreting experimental and numerical data of such systems. The insights from this work can be used to improve the design of fabricated bipolar devices.
1. Introduction
The transport of ions across permselective materials is immensely important for numerous applications, including, but not limited to, desalination (Reference Nikonenko, Kovalenko, Urtenov, Pismenskaya, Han, Sistat and PourcellyNikonenko et al. 2014; Reference Marbach and BocquetMarbach & Bocquet 2019), energy harvesting (Reference Siria, Poncharal, Biance, Fulcrand, Blase, Purcell and BocquetSiria et al. 2013; Reference Wang, Zhang, Zhu, Zhang, Liu, Fu, Zhu, Pu and QuWang et al. 2023), biomolecule sensing (Reference Vlassiouk, Kozel and SiwyVlassiouk, Kozel & Siwy 2009; Reference Slouka, Senapati and ChangSlouka, Senapati & Chang 2014) and fluidic-based electrical circuits (Reference Vlassiouk, Smirnov and SiwyVlassiouk, Smirnov & Siwy 2008; Reference Yossifon, Chang and ChangYossifon, Chang & Chang 2009; Reference Lucas and SiwyLucas & Siwy 2020; Reference Noy and DarlingNoy & Darling 2023; Reference Sebastian and GreenSebastian & Green 2023). However, the electrokinetic transport of ions depends on the geometry of the system, the distribution of the surface charge density, the electrolyte concentration, the applied voltage drop and many other parameters, such that it is virtually impossible to completely characterize the electrical response of the system for the full parameter space. The difficulty of characterizing the system is further frustrated by the fact that the governing equations of the transport are a set of coupled, nonlinear partial differential equations that require numerical evaluation.
Naturally, this has led to the reliance on either experiments or simplified numerical simulations (that focus on very specific systems with a small parameter space). The goal of this two-part work has been to provide a robust and systematic scan of a large parameter space to demonstrate how the electrical response of a permselective system varies as one particular parameter is varied. To that end, we divided the results and analysis into two. In Part 1 (Reference Abu-Rjal and GreenAbu-Rjal & Green 2024), we focused on the much more studied and intuitive steady-state characteristics. Here, in Part 2, we will focus on the time-transient result leading up to the steady-state response. Similar to Part 1, here too, we will show that every system considered has a different signature.
As in Part 1, we distinguish between three different systems – the unipolar system and the bipolar system where the bipolar system is then subdivided into two: the ideal and non-ideal. The difference between the unipolar and bipolar systems is the distribution of the surface charge density. Figure 1 shows a schematic of a bipolar system comprising two diffusion layers and two oppositely charged selective regions. When considering a unipolar system, one needs to remove one of the charged regions (without loss of generality, we will remove the positively charged surface by setting 
${L_3} = 0$ or 
${N_3} = 0$, where 
${L_3}$ and 
$N_3$ represent the non-dimensional length and volumetric excess counterion charge density in region 3, respectively). At the two ends of the system are two electrodes under an applied voltage drop of V, which measure a current, I (or vice versa, an applied current and measured voltage). Each of these systems is subjected to a wide range of voltages without and with the effect of electroosmotically driven advection while we characterize 
$I(t)$. To remove the dependency of the current, I, on the width, W, and height, H, of the system, we will consider the current density, 
$i(t)$ (discussed further below).
Figure 1. Schematic of a three-dimensional four-layered system comprising two diffusion layers connected by two permselective mediums under an applied voltage drop, V. The length of each of the four regions (
$k = 1,2,3,4$) is given by 
${L_k}$, the height is H and the width is W. The two outer regions are uncharged such that the concentrations of the positive ions (purple spheres) and the negative ions (green spheres) are the same. The two middle regions are charged with either a negative or positive surface charge density, leading to a surplus of counterions over coions. In the negatively charged region, the positive ions are the counterions, while in the positively charged regions, the negative ions are the counterions.
In Part 1, we demonstrated that the steady-state current density–voltage response 
$\langle \bar{i}\rangle - V$ of unipolar and bipolar systems differed from each other such that each configuration had a unique steady-state signature. In particular, the differences are by far more drastic when the effects of electroosmotic flows (EOF) are accounted for. Specifically, we showed that in contrast to unipolar systems, where over-limiting currents (OLCs) are always observed, in ‘ideal’ bipolar systems (the definition for ‘ideal’ is given below), OLCs are not observed at all. Further, any ‘non-ideal’ bipolar system displays characteristics of both unipolar and ideal bipolar systems, where the exact response is dependent on the degree of non-ideality. The three curves are shown in figure 2. In Part 2, we will explain these steady-state results by considering the time-transient response of the current subjected to a constant voltage applied at 
$t = 0$. We will show that the appearance of the electroosmotic instability (EOI) and whether it is sustained for long times, will depend on whether at least one initial depletion layer is formed in the system, and whether this initial depletion layer can be sustained. If double enrichment layers are formed, which is the case for negative voltages, the instability will not appear.
Figure 2. Steady-state current density–voltage, 
$\langle \bar{i}\rangle - V$, results without EOF (dashed lines) and with EOF (solid lines) for three scenarios: unipolar ( 
${N_3} = 0$), non-ideal bipolar ( 
${N_3} ={-} 10$) and ideal bipolar ( 
${N_3} ={-} 25$) systems. The inset is a zoomed view of the negative voltage near 
$\langle \bar{i}\rangle = 0$, showing that the 
$\langle \bar{i}\rangle - V$ curves of bipolar systems do not exhibit OLCs there. The error bars denote one standard deviation of the current.
It should be noted that all the results presented in Part 1, and those in figure 2, are the steady-state responses of all the various systems considered (unipolar, ideal bipolar and non-ideal bipolar). In order to calculate the steady-state response (shown in Part 1), we first calculated the equilibrium state (i.e. 
$\langle \bar{i}\rangle = V = 0$), and then at 
$t = 0$ we applied a non-zero voltage drop across the system 
$V \ne 0$. In Part 2, we shall focus on the time-transient responses, leading up to the more intuitive steady-state response of Part 1.
Since the steady-state response and the numerical simulations are already discussed in Part 1, in § 2, we only briefly review the theoretical model and numerical simulations (see Part 1 and its Supplementary Material for a detailed discussion regarding the methods). Thereafter, we divide the paper into three sections that discuss each of the three scenarios depicted in figure 2: unipolar, ideal bipolar and non-ideal bipolar. We start in § 3 with the time-dependent behaviour of the highly investigated unipolar system. While most of the results here have been presented in many past works (references are provided below), a thorough discussion and review are extremely useful as most of the phenomena observed in unipolar systems are observed, in one form or another, in bipolar systems. Thus, the unipolar section ‘sets the stage’ and provides all of the required background needed to understand the bipolar response. In § 4, we move to the other extremity of an ‘ideal’ bipolar response and consider the time-dependent behaviour for positive voltages 
$V > 0$. We shall show that, at initial times, two depletion regions are formed, and the EOI appears on both sides of the permselective material. Eventually, due to the symmetry of the system, the EOI decays completely and results in 
$\langle \bar{i}\rangle - V$ being unchanged relative to the scenario without EOF. Finally, in § 5, we discuss the transient response of a non-ideal bipolar system and relate this response to the two previous scenarios. For the sake of brevity, the negative voltage, 
$V < 0$, responses have been moved to the Supplementary Material (SM), available at https://doi.org/10.1017/flo.2024.24. In general, 
$V < 0$ is uninteresting, as double enrichment layers are formed and EOI cannot appear – such that the response with and without EOF are identical (figure 2 inset). Since we provide a preliminary summary for each of the scenarios, the concluding discussion in § 6 is relatively short.
2. Problem formulation
We solve the non-dimensional time-dependent equations that govern ion transport through a permselective medium. These are the Poisson–Nernst–Planck and the Stokes equations for a symmetric and binary electrolyte (
${z_ + } ={-} {z_ - } = 1$) with ions of equal diffusivities (
${\tilde{D}_ \pm } = \tilde{D}$) in a four-layer system, as shown in figure 1. Note that tilded notations are used for dimensional variables, whereas untilded variables are non-dimensional. Our control parameters are the non-dimensional Debye length, 
$\varepsilon $, the Péclet number, 
$Pe$, the non-dimensional voltage, V, and the non-dimensional volumetric excess counterion charge density, N,
\begin{equation}\varepsilon = \frac{{{{\tilde{\lambda }}_D}}}{{\tilde{L}}} = \frac{1}{{\tilde{L}}}\sqrt {\frac{{{{\tilde{\varepsilon }}_0}{\varepsilon _r}\tilde{\Re }\tilde{T}}}{{2{{\tilde{F}}^2}{{\tilde{c}}_0}}}} ,\quad Pe = \frac{{{{\tilde{\varepsilon }}_0}{\varepsilon _r}\tilde{\varphi }_{th}^2}}{{\tilde{\mu }\tilde{D}}},\quad V = \frac{{\tilde{V}}}{{{{\tilde{\varphi }}_{th}}}},\quad N = \frac{{\tilde{N}}}{{\tilde{F}{{\tilde{c}}_0}}}.\end{equation}
Here, 
$\tilde{\Re }$ is the universal gas constant, 
$\tilde{T}$ is the absolute temperature, 
$\tilde{F}$ is the Faraday constant, 
${\tilde{\varepsilon }_0}$ and 
${\varepsilon _r}$ are, respectively, the permittivity of the vacuum and the relative permittivity, 
${\tilde{c}_0}$ is the bulk concentration, 
$\tilde{\mu }$ is the dynamic viscosity of the fluid and 
${\tilde{\varphi }_{th}} = \tilde{\Re }\tilde{T}/\tilde{F}$ is the thermal potential. The spatial variables have been normalized by a characteristic length 
$\tilde{L}$ (which can be any of the lengths in the system). The excess counterion is related to the average volumetric space charge density needed to counterbalance the surface charge density. We remind the reader that time, 
$\tilde{t}$, has been normalized by the diffusion time 
${\tilde{L}^2}/\tilde{D}$, the ionic fluxes have been normalized by 
${\tilde{j}_0} = \tilde{D}{\tilde{c}_0}/\tilde{L}$, while the non-dimensional space charge density, 
${\rho _e}$ has been normalized by 
$\tilde{F}{\tilde{c}_0}$. The non-dimensional velocity vector 
$\boldsymbol{u} = u\hat{\boldsymbol{x}} + v\hat{\boldsymbol{y}}$ and pressure p have been normalized, respectively, by a typical velocity 
${\tilde{u}_0} = ({\tilde{\varepsilon }_0}{\varepsilon _r}\tilde{\varphi }_{th}^2)/(\tilde{\mu }\tilde{L})$ and pressure 
${\tilde{p}_0} = \tilde{\mu }{\tilde{u}_0}/\tilde{L}$.
Also, it is important to note that we have chosen to work in a non-dimensional formulation to reduce the number of parameters to a minimal number. In particular, in all of our simulations, we keep 
$\varepsilon$ and 
$Pe$ constant while we vary V and N (three tables of all simulation parameters are given in the SM of Part 1). Note that a constant non-dimensional number does not imply that all parameters are not varied but rather that a multiplication of all the relevant parameters is not varied. As such, performing the simulation in a non-dimensional manner is more robust than performing them in a dimensional manner.
Our two-dimensional (2-D) system consists of either three or four regions (always two diffusion layers, with one or two permselective regions) representing unipolar or bipolar systems, respectively. A detailed discussion on the boundary conditions is given in Part 1, and it will not be repeated here. We remind the reader that within the permselective regions, we do not account for hydrodynamic effects. Finally, we remind the reader that we will use the 
${\varDelta _k}$ notations to denote cumulative lengths within the system such that 
${\varDelta _1} = {L_1}$, 
${\varDelta _2} = {\Delta _1} + {L_2}$, 
${\varDelta _3} = {\varDelta _2} + {L_3}$, 
${\varDelta _4} = {\varDelta _3} + {L_4}$.
Throughout, we will use the following (1-D and 2-D) spatial (denoted with single and double overbars, respectively) and temporal (denoted with chevron brackets) averaging operators of any quantity f (e.g. 
${c_ \pm }$, 
$\varphi $, .
${\rho _e}$ and their fluxes, etc.) defined by
\begin{gather}\bar{f}(y,t) = \frac{1}{W}\int_0^W {f(x,y,t)\,\textrm{d}\kern 0.06em x} ,\quad {\bar{\bar{f}}_{k = 1,4}}(t) = \frac{1}{{{L_{k = 1,4}}}}\int_{{L_{k = 1,4}}} {\bar{f}(y,t)\,\textrm{d}y} ,\end{gather}
\begin{gather}\langle \,\bar{f}\rangle = \frac{1}{T}\int_{{t_0}}^{{t_0} + T} {\bar{f}\,\textrm{d}t} ,\quad \langle\, \bar{\bar{f}}\rangle = \frac{1}{T}\int_{{t_0}}^{{t_0} + T} {\bar{\bar{f}}\,\textrm{d}t} ,\end{gather}
where, 
${t_0}$ is the time at which the system reaches a state where the state is perturbed about a ‘steady state’ (i.e. a statistical steady state), and T is the time duration of this ‘steady state’ (typically the end of the simulations). Of particular importance are the averages for the non-dimensional normal electrical current density (normalized by 
$\tilde{F}\tilde{D}{\tilde{c}_0}/\tilde{L}$), 
${i_y}(t)$, at 
$y = 0$ given by 
$\bar{i}(t) = {W^{ - 1}}\smallint _0^W{i_y}(x,y = 0,t)\,\textrm{d}\kern 0.06em x$ and the non-dimensional kinetic energy density (normalized by 
$\tilde{\rho }\tilde{u}_0^2$) 
${E_k} = {\textstyle{1 \over 2}}|\boldsymbol{u}{|^2} = {\textstyle{1 \over 2}}({u^2} + {v^2})$. In the following, we will present an analysis of the time-dependent behaviour of the 1-D averages of the cation concentration, 
${\bar{c}_ + }$, the space charge density, 
${\bar{\rho }_e}$, and kinetic energy 
$\overline {{E_k}}$ in regions 1 and 4.
Similar to Part 1, which focuses on the steady-state 
$\langle \bar{i}\rangle - V$, here, we will consider the time-transient response of three scenarios: unipolar, non-ideal bipolar and ideal bipolar. All these systems are characterized by the ratio of the geometry and excess counterion charge density of both permselective regions (see Reference Green, Edri and YossifonGreen, Edri & Yossifon 2015 for the 2-D version of this equation)
\begin{equation}\eta =
\frac{{{L_3}}}{{{L_2}}} \times
\left|{\frac{{{N_3}}}{{{N_2}}}} \right| = \left\{
{\begin{array}{ * {20}{l}} {0,}&{\textrm{unipolar}}\\ {0
< \eta < 1,}&{\textrm{non-ideal}\,\textrm{bipolar}}.\\
{1,}&{\textrm{ideal}\ \textrm{bipolar}} \end{array}}
\right.\end{equation}
It is trivial to see that, when either 
${N_3} = 0$ or 
${L_3} = 0$, the second charged region does not exist, and the bipolar system reduces to the unipolar system. If, however, there are two charged regions, the ratio 
$\eta$ will determine the overall response of the system. If 
$\eta = 1$ (and 
${L_2} = {L_3}$), the total excess counterion charges in both regions are equal (but of opposite sign), and the response is what we now term ‘ideal’ bipolar. If 
$1 > \eta > 0$, we term the response non-ideal bipolar. The non-ideal bipolar system will exhibit time-dependent and steady-state characteristics of both unipolar and ideal bipolar.
3. Time-dependent response of a unipolar system
This section considers the scenario of a unipolar system with only one permselective region. The unipolar system has been investigated both in a one-layer system with a perfect permselective surface (Reference Rubinstein and ZaltzmanRubinstein & Zaltzman 2000; Reference Zaltzman and RubinsteinZaltzman & Rubinstein 2007; Reference Rubinstein, Manukyan, Staicu, Rubinstein, Zaltzman, Lammertink, Mugele and WesslingRubinstein et al. 2008; Reference Chang, Yossifon and DemekhinChang, Yossifon & Demekhin 2012; Reference Pham, Li, Lim, White and HanPham et al. 2012; Reference Demekhin, Nikitin and ShelistovDemekhin, Nikitin & Shelistov 2013; Reference Deng, Dydek, Han, Schlumpberger, Mani, Zaltzman and BazantDeng et al. 2013; Reference Druzgalski, Andersen and ManiDruzgalski, Andersen & Mani 2013; Reference de Valença, Wagterveld, Lammertink and Tsaide Valença et al. 2015; Reference Mani and WangMani & Wang 2020; Reference Sensale, Ramshani, Senapati and ChangSensale et al. 2021; Reference Pandey and BhattacharyyaPandey & Bhattacharyya 2022; Reference Zhang, Zhang, Luo, Yi and WuZhang et al. 2022; Reference Chen, Zhang, Zhang, Luo and YiChen et al. 2023) and a three-layer system with symmetric diffusion layers, 
${L_1} = {L_4}$ (Reference Kim, Wang, Lee, Jang and HanKim et al. 2007; Reference Yossifon and ChangYossifon & Chang 2008; Reference Rubinstein and ZaltzmanRubinstein & Zaltzman 2015; Reference Abu-Rjal, Rubinstein and ZaltzmanAbu-Rjal, Rubinstein & Zaltzman 2016; Reference Abu-Rjal, Prigozhin, Rubinstein and ZaltzmanAbu-Rjal et al. 2017; Reference Demekhin, Ganchenko and KalaydinDemekhin, Ganchenko & Kalaydin 2018). Without loss of generality, we set 
${N_3} = 0$ such that the four-layer system (figure 1) is reduced to a three-layer unipolar system in which a highly cation-permselective medium is flanked by two diffusion layers (figure 3a). Here, we will investigate a three-layer system with asymmetric diffusion layers for the particular case that 
$L_{4}=2 L_{1}$. It is essential to note that our definition of symmetry and asymmetry in the unipolar system is regarding the diffusion layers. In contrast, for bipolar systems, we deal with a completely different kind of ‘symmetry/asymmetry’ discussed in Part 1.
Figure 3. (a) Schematic describing a 2-D three-layered system comprising two asymmetric diffusion layers (regions 1 and 4) connected by a single highly permselective medium (region 2). (b) The electric current density versus time, 
$\bar{i}(t)$, response of the system with and without EOF. The inset of (b) shows a zoom-up on a semilog10 plot. (c–h) The spatially averaged time-dependent profiles of (c,d) the concentration 
${\bar{c}_ + }$, (e,f) the space charge density 
${\bar{\rho }_e}$ and (g,h) the kinetic energy 
$\overline {{E_k}}$ in regions 1 and 4 (based on the first equation of (2.2)). The colours for each line correspond to the markers shown in (b) and the legend is given in (f). The insets in (c), (e) and (g) show profiles for the last time point for 
${c_ + }$, 
${\rho _e}$ and 
${E_k}$, respectively, while the shaded grey regions denote the standard deviation of each variable respectively ( 
${\sigma _c}$, 
${\sigma _\rho }$ and 
${\sigma _E}$). The inset in (h) is the time evolution of the surface averaged kinetic energy, 
$\overline {\overline {{E_k}} } (t)$, in region 1 (based on the second equation of (2.2)). This figure uses 
$V = 60$.
Figures 3 and 4 show the time-transient response of a three-layer highly permselective system with asymmetric diffusion layers at a positive OLC regime. Figure 3(b) plots 
$\bar{i}(t)$ for two scenarios: without EOF (dashed lines) and with EOF (solid lines). Without EOF, it can be observed that the current density decreases monotonically (Reference Abu-Rjal, Leibowitz, Park, Zaltzman, Rubinstein and YossifonAbu-Rjal et al. 2019), this is also true for a ‘symmetric’ three-layered system without EOF (SM, figure S1). In contrast, for the scenario with EOF, the change is highly non-monotonic. This can be attributed to the appearance of the EOI (the time evolution is shown in figure 4 and Movies 1–3 in the SM), which is also responsible for increasing 
$\bar{i}(t)$ (relative to the scenario without EOF).
Figure 4. Time evolution of the positive concentrations (2-D colour plot), 
${c_ + }$, and velocity streamlines (white lines) in regions 1 and 4 for the simulation of figure 3. The coloured titles of each panel follow the colour notation of figure 3. The surface average of the kinetic energy, 
$\overline {\overline {{E_k}} }$ (based on the second equation of (2.2)), is given for each snapshot (subscripts of numbers denote regions).
At 
$t = 0$, the system starts at equilibrium; this is tantamount to setting the concentrations to unity, 
${c_ \pm }(t = 0) = 1$, and 
${\rho _e}(t = 0) = 0$ everywhere outside of the equilibrium electrical double layer (EDL). Upon application of a positive voltage, 
$V > 0$, it can be observed from figure 3(c) that outside the quasi-equilibrium EDL, the concentration in region 1 is always smaller than unity, while in region 4 (figure 3d), the concentration is larger than unity (see Movie 1 in the SM). Naturally, these regions are called the depleted and enriched regions, respectively, and are an essential component of concentration polarization (see Reference LevichLevich 1962; Reference Rubinstein and ZaltzmanRubinstein & Zaltzman 2000; Reference Zaltzman and RubinsteinZaltzman & Rubinstein 2007; Reference Chang, Yossifon and DemekhinChang et al. 2012; Reference Mani and WangMani & Wang 2020 for thorough discussions on concentration polarization).
In the enriched region (region 4), the large concentrations (figure 3d) ensure that 
${\bar{\rho }_e}$ (figure 3f) maintains its 1-D quasi-equilibrium structure such that EOI does not appear and the average kinetic energy is zero at all times (figure 3h). As will be discussed further below, it is essential to realize and remember that enriched regions cannot support non-equilibrium extended-space charges (ESCs) and that the ESC can only form on the depleted side when the concentration at the interface (
$y = {\varDelta _1}$) approaches zero.
Naturally, the dynamics on the depleted side is by far more complicated. Around 
$t = 0.011$, the concentration at the interface (
$y = {\varDelta _1}$) approaches zero (figure 3c). As a result, the space charge density moves from a quasi-equilibrium to a non-equilibrium structure, which is commonly known as the ESC (figure 3e) (Reference Rubinstein and ShtilmanRubinstein & Shtilman 1979; Reference Zaltzman and RubinsteinZaltzman & Rubinstein 2007; Reference YarivYariv 2009). This non-equilibrium structure is unstable to lateral perturbations (Reference Rubinstein and ZaltzmanRubinstein & Zaltzman 2000; Reference Kim, Wang, Lee, Jang and HanKim et al. 2007; Reference Zaltzman and RubinsteinZaltzman & Rubinstein 2007; Reference Rubinstein, Manukyan, Staicu, Rubinstein, Zaltzman, Lammertink, Mugele and WesslingRubinstein et al. 2008; Reference Yossifon and ChangYossifon & Chang 2008; Reference Pham, Li, Lim, White and HanPham et al. 2012; Reference Mani and WangMani & Wang 2020), and eventually, the EOI forms (Figure 4 and Movie 2). As a result, the average kinetic energy, 
$\overline {{E_k}}$, also increases such that the effects of the velocity span the entire depleted region (figure 3g) (Reference Demekhin, Nikitin and ShelistovDemekhin et al. 2013; Reference Druzgalski, Andersen and ManiDruzgalski et al. 2013).
3.1 Preliminary summary
In this section, we have revisited the highly investigated time-dependent response for a unipolar system. We have shown that, upon application of a positive voltage (negative voltages are discussed in the SM), a depleted layer and an enriched layer form in each of the diffusion layers. This is in contrast to the bipolar systems, which at early times will have two depleted layers and at later times have no depleted and enriched layers (ideal bipolar) or a depleted and an enriched layer (non-ideal bipolar). Above a critical voltage, an EOI forms only on the depleted sides. The EOI, which introduces and injects kinetic energy into the system, leads to OLCs. Importantly, so long as a supercritical voltage is applied, the EOI does not decay. This stands in contrast to the dynamics that will be observed for bipolar systems.
Finally, for the sake of completion and brevity, we note that the unipolar system differs from the bipolar systems in that the EOI can form in either diffusion layer if the positive voltage is changed to negative voltage. Also, we note that breaking the geometric symmetry of the diffusion layers leads to quantitatively different responses that are without any qualitative change. All these are shown in the SM.
4. Time-dependent response of an ideal bipolar system
This section focuses on an ideal bipolar system where the system satisfies 
$\eta = 1$ (i.e. 
${N_3} ={-} 25$). We will show that the ‘ideal’ bipolar system has completely different characteristics compared with those of the unipolar system (§ 3). The layout of figures 5 and 6 are similar to those of figures 3 and 4, respectively, save that we now consider an ideal bipolar four-layer set-up (figure 5a).
Figure 5. (a) Schematic describing a 2-D four-layered ideal bipolar system comprising two diffusion layers (regions 1 and 4) connected by two permselective mediums (regions 2 and 3). (b) The electric current density versus time, 
$\bar{i}(t)$, response of the system with and without EOF. The inset of (b) shows a zoom-up of early times on a semilog10 plot. (c–h) The spatially averaged time-dependent profiles of (c,d) the concentration 
${\bar{c}_ + }$, (e,f) the space charge density 
${\bar{\rho }_e}$ and (g,h) the kinetic energy 
$\overline {{E_k}}$ in regions 1 and 4. Note that, in (e), we present 
${\bar{\rho }_e}$, while, in (f), we present 
$- {\bar{\rho }_e}$. The colours for each line correspond to the markers shown in (b) and the legend is given in (f). The insets in (g) and (h) are the time evolution of the surface average of the kinetic energy, 
$\overline {\overline {{E_k}} } (t)$, in regions 1 and 4, respectively. The blue curved arrows indicate the direction of increasing time, showing the non-monotonic response. Here, we have used 
$V = 100$.
Figure 6. Time evolution of the positive concentrations (2-D colour plot), 
${c_ + }$, and velocity streamlines (white lines) in regions 1 and 4 for the simulation of figure 5. The coloured titles of each panel follow the colour notation of figure 5. The surface average of the kinetic energy, 
$\overline {\overline {{E_k}} }$, is given for each snapshot (subscripts of numbers denote regions).
Figure 5(b) plots 
$\bar{i}(t)$ for two scenarios: without EOF (dashed lines) and with EOF (solid lines). The scenario without EOF shows a sharp decrease and then a monotonic rise to the predicted steady-state current (Reference Green, Edri and YossifonGreen et al. (2015) and Reference Abu-Rjal and GreenAbu-Rjal & Green (2021) provide the expression for the convection-less 
$\langle \bar{i}\rangle - V$ – which extends the 
$\langle \bar{i}\rangle - V$ of Reference Vlassiouk, Smirnov and SiwyVlassiouk et al. (2008), from one dimension to two and accounts for the resistances of the diffusion layers). The scenario with EOF shows a similar drop and rise with several remarkable features. First, unsurprisingly, when accounting for EOF, the EOI appears (figure 6 and Movies 4–6). However, surprisingly, the EOI appears in both diffusion layers. Second, with the appearance of EOI, we see a relative increase of 
$\bar{i}(t)$ (relative to the scenario without EOF). This, too, is unsurprising since EOF is more efficient than diffusion in stirring the electrolyte. Third, since EOI is an unstable process, it is natural that the changes are non-monotonic and ‘noisy’ (inset of figure 5b). Fourth, surprisingly, the steady-state current with EOF equals the current of the no-EOF scenario. To understand this last point, it is essential to understand the behaviour of 
${\bar{c}_ \pm }$, and 
${\bar{\rho }_e}$ in regions 1 and 4.
At 
$t = 0$, the system starts at equilibrium. For the concentrations, this is tantamount to 
${c_ \pm }(t = 0) = 1$ and 
${\rho _e}(t = 0) = 0$ everywhere outside of the equilibrium EDL. Figure 5(c,d) shows the behaviour of 
${\bar{c}_ + }$ in regions 1 and 4, respectively, while figure 5(e,f) shows the behaviour of 
${\bar{\rho }_e}$ in regions 1 and 4, respectively. In contrast to the unipolar case considered in figure 3, where one of the diffusion layers is depleted and the other is enriched, here we observe a much more complicated dynamics – including double diffusion layers and much more.
Figure 5(c,d) shows that both diffusion layers become depleted. Figure 5(e,f) shows that in both layers, an ESC layer is formed (in the unipolar system, only one ESC was formed – in the single depletion layer). Equally remarkable, the dynamics of the concentration fields is symmetric (for the unipolar case, there were both depleted and enriched layers), while the space charge density fields are antisymmetric 
$({\bar{\rho }_{e,1}} ={-} {\bar{\rho }_{e,4}})$. At initial times, the effects of depletion increase such that the interfacial concentrations decrease, but the space charge still has an equilibrium structure (
$t < 0.004$). At intermediate times (
$t \sim 0.004$), a non-equilibrium space charge region is formed (pink lines figure 5c–f). Since the non-equilibrium space charge region is unstable, EOI appears in both regions! For times larger than 
$t > 3$, the effect of EOI appears to decay. The reason for this is explained below. Finally, at steady state after the EOI has completely decayed (best observed by the decay of the kinetic energy in figure 5g,h), 
${\bar{c}_ \pm }$ and 
${\bar{\rho }_e}$ return to their quiescent states.
Movie 4 in the SM plots the time-dependent behaviour of the distribution of 
${\bar{c}_ + }$ and 
${\bar{\rho }_e}$ for the two scenarios: without and with EOF. Except for the appearance of EOI and its eventual decay, the qualitative nature of the EOF results is almost qualitatively identical to the response of an ideal convection-less bipolar system (Reference Abu-Rjal and GreenAbu-Rjal & Green 2021; Reference Tepermeister and SilbersteinTepermeister & Silberstein 2023). We found this result somewhat surprising, leading to the question we alluded to in Part 1 (Reference Abu-Rjal and GreenAbu-Rjal & Green 2024) – why are the steady-state responses without and with EOF identical?
The answer is surprisingly simple and is best understood by considering the ideal convection-less bipolar system for the particular case 
$\eta = 1$. In our past work (Reference Green, Edri and YossifonGreen et al. 2015), we showed that one can derive an analytical 
$\langle \bar{i}\rangle - V$ for the case of 
$\eta = 1$ if one assumes that the salt current density, 
$j = {j_ + } + {j_ - }$, is zero (
$j = 0$). In this scenario, positive ions are transported via 
${j_ + }$ from region 1 to region 4 through the bipolar region. During the time they are transported to the interface of regions 3 and 4, the (negative) ESC has formed. The arrival of the positive ions stabilizes the ESC by reducing its magnitude. A similar process occurs for negative ions via 
${j_ - }$. By virtue of symmetry but opposite directions, 
${j_ + } ={-} {j_ - }$, the response is symmetric in both regions 1 and 4.
The time-dependent dynamics is a bit more complicated; however, the simplicity of the explanation remains. At 
$t = 0$, prior to the application of the voltage drop, V, the fluxes are zero such that 
$i = j = 0$. As stated above, at steady state, there is still a constraint that 
$j = 0$. However, in between 
$t = 0$ and 
$t \to \infty $, 
$j(t)$ varies locally (and globally) such that it increases before going back to zero. Since j and i are linked, this also leads to changes in i. The local behaviour can be observed in SM Movie 6, while figure S4 shows the correlation between 
$\bar{\bar{j}}(t)$ in region 1 and 
$\bar{i}(t)$ for both scenarios (without and with convection). It can be observed that the non-monotonic changes in 
$\bar{i}(t)$ are related to the non-monotonic changes in 
$\bar{\bar{j}}(t)$ due to the instability. A thorough discussion of the time-dependent dynamics of all the fluxes, without the EOF, is given in Reference Abu-Rjal and GreenAbu-Rjal & Green (2021).
For the current scenario of an ideal bipolar system with EOF, one must return to the problem definition in § 2 (and Part 1 § 3.2), where we assumed that the velocity within the permselective material is zero. As a result, the response of the bipolar membrane (regions 2 and 3) must remain unchanged relative to the response of the bipolar membrane in the convection-less scenario. Here, too, positive ions are transported via 
${j_ + }$ from region 1 to region 4 in such a manner as the ESC in region 4 is stabilized. Naturally, the exact same thing occurs for negative ions via 
$\,{j_ - }$ which stabilize the positive ESC in region 1. In this manner, both ESCs are stabilized and diminished such that the electric body force that drives the EOF instability is diminished, and the instability decays – this can be observed in the decay of the kinetic energy (figure 5g,h).
Before summarizing this section, three last comments are needed. First, we start with a technical comment. In general, for the ideal bipolar scenario, the response in both diffusion layers should have a mirror reflection (and possibly a sign difference) – this is what occurs for a bipolar system without EOI. However, in several time steps in figure 6, the 
$\overline {\overline {{E_k}} }$ in regions 1 and 4 are not identical. Since the spatial averages of the concentrations (figure 5c,d) and space charge densities (figure 5e,f) are the same, as are the averages for the velocities, we believe this few per cent error can be attributed to meshing and the fact that we are considering the square of the velocities (where any mismatch is drastically enhanced). Second, we point out that the appearance and decay of the EOI do not occur for a specific voltage. Rather, it occurs for all the voltages we considered. In the SM, we demonstrate this statement for several voltages. Third, for an ideal bipolar system subject to negative voltage drops, 
$V < 0$, the EOI does not appear at all. This is because for 
$V < 0$, double enrichment layers, that do not support an ESC, forms. This, too, is demonstrated in the SM (and discussed thoroughly in our past work (Reference Abu-Rjal and GreenAbu-Rjal & Green 2021)).
4.1 Preliminary summary
In summary, the response of the ideal bipolar system differs substantially from the unipolar system in that either double depletion (for 
$V > 0$) or double enrichment (for 
$V < 0$) layers are formed. When double depletion layers are formed, the EOI appears in both layers. However, in contrast to the unipolar, even when an over-critical voltage is applied, the EOI decays such that the steady-state OLC cannot be observed. Importantly, the time-dependent response of the system with EOI includes a clear signature that can be observed relative to the rather quiescent scenario without EOI (figure 5b).
5. Time-dependent response of a non-ideal bipolar system
This section considers the non-ideal bipolar case, 
$1 > \eta > 0$ (figure 7a). We will show that the response of this system displays some characteristics observed in the ideal bipolar case (§ 4) and characteristics observed in the unipolar system (§ 3). To highlight the main shared features and differences, figures 7 and 8 have the same layout as figures 3 and 4 in § 3 and figures 5 and 6 in § 4.
Figure 7. (a) Schematic describing a 2-D four-layered non-ideal bipolar system comprising two diffusion layers (regions 1 and 4) connected by two permselective mediums (regions 2 and 3). (b) The electric current density versus the time, 
$\bar{i}(t)$, response of the system with and without EOF. The inset of (b) shows a zoom-up on a semilog10 plot. (c–h) The spatially averaged time-dependent profiles of (c,d) the concentration 
${\bar{c}_ + }$, (e,f) the space charge density 
${\bar{\rho }_e}$ and (g,h) the kinetic energy 
$\overline {{E_k}}$ in regions 1 and 4. The colours for each line correspond to the markers shown in (b), and the legend is given in (f). The insets in (c), (e) and (g) show the profiles for the last time point for 
${c_ + }$, 
${\rho _e}$ and 
${E_k}$, respectively, while the shaded grey regions denote the standard deviation of each variable respectively (i.e. 
${\sigma _c}$, 
${\sigma _\rho }$ and 
${\sigma _E}$). The inset in (h) is the time evolution of the surface averaged kinetic energy, 
$\overline {\overline {{E_k}} } (t)$, in region 1. The blue curved arrows indicate the direction of increasing time. Here, we have used 
${N_3} ={-} 10$ and 
$V = 60$.
Figure 8. Time evolution of the positive concentrations (2-D colour plot), 
${c_ + }$, and velocity streamlines (white lines) in regions 1 and 4 for the simulation of figure 7. The coloured titles of each panel follow the colour notation of figure 7. The surface average of the kinetic energy, 
$\overline {\overline {{E_k}} }$, is given for each snapshot (subscripts of numbers denote regions).
Figure 7(b) plots 
$\bar{i}(t)$ for the two scenarios of without EOF (dashed lines) and with EOF (solid lines). The scenario without EOF shows a sharp decrease and then a monotonic rise to the steady-state current. Here, too (similar to the 
$\eta = 1$ scenario), the non-monotonic change of the current should also be attributed to the non-monotonic change of 
$\bar{\bar{j}}(t)$ (figure S7), which shows that 
$\bar{\bar{j}}(t = 0) = 0$ but changes until it reaches its steady-state value. Similarly, the scenario with EOF exhibits a non-monotonic change similar to the unipolar and ideal bipolar scenarios. This can be related to the appearance of EOI (figure 8 and Movies 7–9). Importantly, in contrast to the ideal bipolar scenario (figure 5b) where the steady-state currents without EOF and with EOF were identical, here we observe that the steady-state current without EOF is substantially smaller than that with EOF – as was observed in the unipolar scenario (figure 3b).
To better understand the change in the behaviour of 
$\bar{i}(t)$, we consider the behaviour of 
${\bar{c}_ + }$, and 
${\bar{\rho }_e}$ in regions 1 and 4. Once more, at 
$t = 0$, the system starts at equilibrium (
${c_ \pm }(t = 0) = 1$ and 
${\rho _e}(t = 0) = 0$ everywhere outside of the equilibrium EDL). The 
$t > 0$ results show a combination of both unipolar and ideal bipolar behaviours:
Concentrations. Figure 7(c) shows that the depletion layer in region 1 forms almost monotonically, and the concentration does not return to its equilibrium state 
$[{c_ + }(t = 0) = 1]$. This behaviour is reminiscent of the unipolar response. The behaviour of the concentration in region 4 (figure 7d), is by far more complicated. Similar to the ideal bipolar scenario, the interfacial concentration decreases to form a depletion layer. Once the depletion layer has achieved a minimal value, a reversal is observed. Then, in contrast to the ideal bipolar scenario and more reminiscent of the unipolar response, the concentration does not return to its equilibrium state 
$[{c_ + }(t = 0) = 1]$ and an enrichment layer is formed.
Space charge density. The behaviour of the space charge density (figure 7e,f) follows the behaviour of the concentrations in that the dynamics is by far more complicated. In particular, we observe that 
${\bar{\rho }_{e,1}} \ne - {\bar{\rho }_{e,4}}$. The ESC in region 1 (figure 7e), whose dynamics is slightly non-monotonic, does not decay at later times (this is reminiscent of the unipolar response). The space charge density in region 4 (figure 7f) mirrors the dynamics of the concentrations: when the layer is depleted, an ESC forms; when the layer becomes enriched, the ESC disappears.
Kinetic energy. Naturally, the behaviour of the velocity must follow that of the space charge density (figure 8). As in the ideal bipolar scenario, EOI appears in both regions 1 and 4. In contrast to the ideal bipolar scenario and similar to the unipolar scenario, the EOI in region 1 does not decay (figure 7g). Similar to the ideal bipolar scenario and the unipolar scenario, the EOI in region 4 decays (figure 7h). Importantly, the permanence of the EOI in region 1 leads to persistent OLC in the non-ideal bipolar case (figure 7b).
It should be noted that, since 
$1 > \eta > 0$, there is an inherent asymmetry between 
${j_ + }$ and 
${j_ - }$ fluxes. In particular, we note that the 
${j_ + }$ fluxes, leaving region 1 and arriving in region 4, deplete the ESC in region 1 quicker than the 
${j_ - }$ fluxes, leaving region 4 and arriving in region 1, can stabilize the ESC. In contrast, it is easy to observe that the 
${j_ + }$ fluxes stabilize the ESC in region 4 quicker than the depletion of 
${j_ - }$ fluxes. As a result, the EOF can be observed at steady state in region 1 while it is not observed in region 4. The difference in the behaviour of 
${j_ + }$ and 
${j_ - }$ in regions 1 and 4 also leads to differences in the behaviour of 
$\bar{\bar{j}}$ and 
$\bar{i}$ in these regions. Similar to figure S4 shown in § 2 in the SM, figure S7 presents the time dependence of these quantities, showing that differences in the behaviour of 
${\bar{\bar{j}}_1}(t)$ and 
${\bar{\bar{j}}_4}(t)$, which are also responsible for the non-monotonic change in 
$\bar{i}(t)$. See Movie 9 for more details.
To complete our analysis of non-ideal bipolar systems, we need two more things: more positive voltages and negative voltages. For the sake of completion, in the SM, we provide several 
$\bar{i}(t)$ curves for several additional positive voltages where we demonstrate that the overall characteristics remain the same with the sole difference that when the velocity field becomes chaotic, so does 
$\bar{i}(t)$. Here, similar to the ideal bipolar response, for negative applied voltages, both diffusion layers exhibit enrichment at early times. At later times, the small limiting currents, determined by the permselective regions (and by the diffusion layers, as in the case of the unipolar system), ensure that the concentration gradients are extremely small, and the space charge density is always in a quasi-equilibrium-like structure. Hence, instability is not observed.
5.1 Preliminary summary
Here, we have shown that the response of the non-ideal bipolar scenario includes many characteristics that can be attributed to either the unipolar response or the ideal bipolar response (e.g. two time-transient depletion layers, which can be found in the ideal bipolar system, but one depletion and one enrichment layer found in the steady-state response of the unipolar system). The reason that the depletion layers continue to be sustained in one of the regions – is because there is a symmetry breaking between 
${j_ + }$ and 
${j_ - }$, due to 
$\eta \ne 1$, leading to both 
$j \ne 0$ and at least one persistent depletion layer. Importantly, so long as there is a single depletion layer (that is completely depleted), there will also be a non-equilibrium space charge layer that can sustain the EOI (which leads to OLCs).
Admittedly, while our results have provided new insights into the behaviour of the EOI in bipolar systems, many questions remain unanswered. Notably, the conditions that lead to the instability and the various scaling laws observed in Part 1, are not understood. To that end, a thorough linear stability analysis is needed to determine the effects of the control parameters given in (2.1). However, this is left for future work.
6. Discussion and summary
This work focuses on elucidating the behaviour of the time-transient response of ideal and non-ideal bipolar systems subjected to EOI (§§ 4 and 5). To that end, we simulated several ideal and non-ideal bipolar systems for a wide range of V and 
$\eta $. To highlight the surprising and remarkable results uncovered in our simulations of the bipolar system, we first leveraged our understanding of the time-transient response of the more investigated unipolar system (§ 3). There it has been known, for quite some time, that the EOI is responsible for OLCs.
In ideal bipolar systems, regardless of whether EOF is included, double depletion layers are formed. Whether EOI forms or not, at steady-state, the system returns to its equilibrium state, such that OLCs cannot be observed. However, it leaves a clear and unique EOI signature: EOI can be observed in both diffusion layers and in 
$\tilde{i}(t)$ (figure 5b). The response of non-ideal bipolar systems comprises several characteristics of the unipolar and ideal bipolar responses. This includes the formation of EOI on both sides, but the eventual decay of EOI on one side and the appearance of OLCs (figure 2).
In Part 2, we have elucidated the time-transient response of unipolar and bipolar systems subject to over-critical voltages, which, in the right conditions, lead to OLCs. Nonetheless, whether EOI is desired or not and whether it can be leveraged is subjective based on the application (desalination, energy harvest, biosensing and more). Optimization of the effects of EOI and delineating its variations on the performance of each application require a lot more fine tuning associated with the particular application. Yet, this work focuses on more objective and robust issues. Namely, we show, based on the system parameters, whether or not EOI appears in bipolar systems – without this knowledge, an advanced and accurate experimental protocol cannot be formulated.
It is our hope that the systematic approach undertaken in this work will provide much-needed novel insights into the behaviour of these systems as well as how to (optimally) characterize them through their very unique time-dependent signatures that can be observed or measured. These time-dependent signatures, in conjunction with fluorescent microscopy, can be used to confirm or validate whether or not the EOI forms in experimental bipolar nanofluidics systems.
Supplementary Material and Movies
Supplementary material and movies are available at https://doi.org/10.1017/10.1017/jfm.2024.1006.
Acknowledgements
We thank O. Lavi for several in-depth discussions and very useful comments on how to streamline the results of this two-part work.
Funding
This work was supported by Israel Science Foundation grants 337/20 and 1953/20. We acknowledge the support of the Ilse Katz Institute for Nanoscale Science & Technology.
Competing interests
The authors declare no conflict of interest.
Data availability statement
Raw data will be made available upon request from the corresponding author (Y.G.).
Open access
