# Poisson-Nernst-Planck charging dynamics of an electric double layer capacitor: symmetric and asymmetric binary electrolytes

Ivan Palaia, Adelchi Asta, Megh Dutta, Patrick Warren, Benjamin Rotenberg,
Emmanuel Trizac

► To cite this version:

Ivan Palaia, Adelchi Asta, Megh Dutta, Patrick Warren, Benjamin Rotenberg, et al. Poisson-Nernst-Planck charging dynamics of an electric double layer capacitor: symmetric and asymmetric binary electrolytes. Physical Review E, 2025, 112 (3), pp.035417. 10.1103/p4dg-snqf. hal-04358062

HAL Id: hal-04358062

https://hal.science/hal-04358062v1

Submitted on 13 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Poisson-Nernst-Planck charging dynamics of an electric double-layer capacitor: Symmetric and asymmetric binary electrolytes

Ivan Palaia<sup>1,1,*</sup> Adelchi J. Asta,<sup>2</sup> Megh Dutta,<sup>3</sup> Patrick B. Warren,<sup>4</sup> Benjamin Rotenberg,<sup>5</sup> and Emmanuel Trizac<sup>6</sup>

<sup>1</sup>Department of Physics, King's College London, London WC2R 2LS, United Kingdom and Institute of Science and Technology Austria, 3400 Klosterneuburg, Austria

<sup>2</sup>Dessia, 92160 Antony, France

<sup>3</sup>Sorbonne Université, CNRS, Physico-chimie des Électrolytes et Nanosystèmes Interfaciaux, PHENIX, 75005 Paris, France

<sup>4</sup>Hartree Centre, Science and Technology Facilities Council (STFC), Sci-Tech Daresbury, Warrington WA4 4AD, United Kingdom

<sup>5</sup>Sorbonne Université, CNRS, Physico-chimie des Électrolytes et Nanosystèmes Interfaciaux, PHENIX, 75005 Paris, France and Réseau sur le Stockage Electrochimique de l'Energie (RS2E), FR CNRS 3459, 80039 Amiens Cedex, France

<sup>6</sup>Université Paris-Saclay, CNRS, LPTMS, 91405 Orsay, France and Ecole Normale Supérieure de Lyon, 69364 Lyon, France

<figure><img src="image_2.png" /></figure>

(Received 15 March 2023; revised 1 June 2025; accepted 20 August 2025; published 29 September 2025)

A parallel plate capacitor containing an electrolytic solution is the simplest model of a supercapacitor or electric double-layer capacitor. Using both analytical and numerical techniques, we solve the Poisson-Nernst-Planck equations for such a system, describing the mean-field charging dynamics of the capacitor, when a constant potential difference is abruptly applied to its plates. Working at constant total number of ions, we focus on the physical processes involved in the relaxation and, whenever possible, give its functional shape and exact time constants. We first review and study the case of a symmetric binary electrolyte, where we assume the two ionic species to have the same charges and diffusivities. We then relax these assumptions and present results for a generic strong (i.e., fully dissociated) binary electrolyte. At low electrolyte concentration, the relaxation is simple to understand, as the dynamics of positive and negative ions appear decoupled. At higher electrolyte concentration, we distinguish several regimes. In the linear regime (low voltages), relaxation is multiexponential, it starts by the buildup of the equilibrium charge profile and continues with neutral mass diffusion, and the relevant timescales feature both the average and the Nernst-Hartley diffusion coefficients. In the purely nonlinear regime (intermediate voltages), the initial relaxation is slowed down exponentially due to increased capacitance, while bulk effects become more and more evident. In the fully nonlinear regime (high voltages), the dynamics of charge and mass are completely entangled and, asymptotically, the relaxation is linear in time. We finally discuss nonideal behavior in real capacitors and provide conditions for which mean-field is expected to hold.

DOI: 10.1103/p4dg-snqf

## I. INTRODUCTION

In a wealth of physical systems, charged surfaces confine a liquid containing ions. If the surfaces are conductive, then a so-called electric double-layer capacitor (EDLC) is formed, owing its name to the two layers of opposite charge that build up at the interface between each surface and the electrolytic solution. Compared to standard capacitors with insulating dielectrics, the ability of EDLCs to store charge and energy is enhanced by the local rearrangement of the confined electrolyte and very large capacities per unit weight can be reached when the conductive electrodes are made of porous or fibrous materials, as this hugely increases the surface area in contact with the electrolyte [1–3]. EDLCs

of this kind can store more energy than common capacitors, can release it more rapidly than common batteries, and last longer than devices based on chemical reactions. These three characteristics, placing them precisely at the border between capacitors and batteries [4], explain their widespread usage as well as numerous potential applications: from consumer electronics and wearable devices [5] to energy production [6–8] and means of transportation, where large power exchanges are needed to propel and halt electric vehicles, harvest braking energy, or promptly activate emergency devices [1].

The need to control the dynamics of charge and discharge of EDLCs has motivated substantial interest in fundamental sciences. Impedance spectroscopy experiments are routinely used to measure the linear response of the system to oscillating fields [9,10]. Computer simulations have addressed paramount questions concerning *inter alia* the effects of electrode polarizability, ion and pore sizes, or electrostatic correlations [11]. In parallel, analytical studies have proved fundamental to inform experiments and simulations, emphasizing the role played by the molecular structure of the double layer (not captured by usual mean-field approaches) in deter-

<footnote>*Contact author: ivan.palaia@kcl.ac.uk</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Cartoon representing a charged, planar EDLC. Red cations and blue anions are treated at mean-field level, as suggested by the color of the solution, representing charge density. Within the mean-field approximation, for large-enough $L$ and sufficiently low applied voltage, the thickness of the electric double layer at equilibrium is $\lambda_D$.</figcaption></figure>

mining the EDLC capacitance [12,13]. Nonetheless, recent works have highlighted the need to better understand the dynamics of simple model systems, such as a planar capacitor in mean field, to clearly isolate the effect of geometry, hydrodynamics and nonidealities in real devices [14–20] or to develop strategies to speed up relaxation to equilibrium [21–23]. Namely, the relaxation dynamics of a planar EDLC subjected to a sudden change in the electric potential between plates has been described by the seminal works [24,25] for the simpler case of symmetric electrolytes and with particular attention for the linear (low-voltage) regime, of which Ref. [25] provided an exact solution.

The charging dynamics of an ideal EDLC, such as the one presented in Fig. 1, is an interesting problem *per se* from a physical perspective. Several length scales are involved: the Bjerrum length $l_B$, setting the distance between two ions at which thermal energy becomes comparable with electrostatic repulsion; the Debye screening length $\lambda_D$, setting the range of the electrostatic potential; the point along $z$ at which the electrostatic potential vanishes, coinciding with the geometrical center of the capacitor only for symmetric electrolytes; and, last, the size of the system $L$. In strongly polarized situations, a fifth length can in principle be relevant, the Gouy-Chapman length, a signature of the less efficient screening occurring when only counterions are in proximity of the plates. The relaxation to equilibrium observed when an external voltage is imposed can depend in principle on any of these lengths, as well as on the diffusion coefficients of the ionic species and on the applied voltage. These quantities can be combined in many ways to give plausible relaxation times: Understanding which of those are relevant for the dynamics, for different electrolytes and at different times, is a challenge. Moreover, the relaxation is not simply exponential, as in many model physical systems, and can be hard to solve and describe in closed form.

In this paper, and in its companion Letter [26], we present a detailed study of the relaxation dynamics of an ideal,

planar EDLC in mean field. In Sec. II we recapitulate the Poisson-Nernst-Planck formalism, within which we operate throughout the whole paper using the numerical method described in the Appendix. Then we present results concerning the case where the two ions species have the same valence and diffusion coefficients (fully symmetric case, Sec. III), the case where the two species have the same diffusion coefficient but different valences (partially asymmetric case, Sec. IV) and the case where the two species have different valences and diffusion coefficients (fully asymmetric case, Sec. V). For each of these cases we vary electrolyte concentration and applied voltage and we study the linear, purely nonlinear, and depleted nonlinear regimes. The relaxation involves in general more than one process, so that each regime features many relaxation timescales, each of which may become relevant at different times. We give a thorough characterization of this complex phenomenon, focusing with more care on the cases for which a simple analytical interpretation is possible. We conclude by discussing the range of validity of our mean-field results in Sec. VI.

## II. MODEL

Our ideal EDLC (Fig. 1) is treated within the Poisson-Nernst-Planck formalism [27]. We neglect hydrodynamic effects and suppose purely Coulombic interactions between ions and Coulombic and hard-core interactions between ions and walls. No absorption phenomena are considered, so the compact part of the electric double layer, usually called Stern layer, is neglected. The two plates, distant $2L$ from each other, are connected to a time-dependent ideal voltage source, imposing a potential difference $2V(t)$ between them.

The salt is a strong binary electrolyte and is therefore completely dissociated in the solvent: densities are denoted $n_+$ and $n_-$, diffusion coefficients $D_+$ and $D_-$, and integer ionic valences $\bar{q}_+ = +q_+$ and $\bar{q}_- = -q_-$ (for instance, for $MgCl_2$, $\bar{q}_+ = +2$, $\bar{q}_- = -1$, and $q_+ = 2$, $q_- = 1$). Electrodes are modelled as parallel, infinite, charged planes. Densities $n_\pm$ and potential $\phi$ are functions of the sole spatial coordinate $z$ in the direction perpendicular to the planes and of time $t$. Figure 1 shows a sketch of the system, with size $2L$ and $z = 0$ at the center of the capacitor.

We work in the canonical ensemble, i.e., with a fixed number of ion pairs per unit surface $2n_0L$, where $n_0$ is the uniform density of salt when the power source is off and the capacitor at equilibrium. The initial density for each species is $n_\pm^0$ and is such that $q_+n_+^0 = q_-n_-^0 = q_+q_-n_0$ by electroneutrality (the last equality holding when $q_+$ and $q_-$ are coprime). Choosing to work in the canonical ensemble allows to avoid the question of where exactly ions are injected into the system from the reservoir during the dynamics, which in a real nanocapacitor might depend on the size and the topology of the pores. In addition, the constant ion number ensemble is a good approximation for open systems where the electrodes have large lateral dimensions and it takes a long time to ensure the chemical equilibrium of the entire pore with the bath. In fact, studying relaxation in the canonical ensemble reveals precisely which timescales should be compared to chemical equilibration timescales to assess whether the system is effectively canonical or grand canonical. Finally, there are
regimes in which densities in the bulk solution (at the center of the capacitor) remain approximately constant, at least to first order, and the canonical ensemble is not quantitatively different from the grand-canonical ensemble.

The Poisson-Nernst-Planck equation describes the change in time and space of the ionic densities [27]: It consists of a drift-diffusion model for ionic currents, complemented with a continuity equation. The average current density $j_\alpha(z, t)$ for ions of type $\alpha$ ($\in \{-, +\}$), moving along $z$ under the action of a potential $\phi(z, t)$, is

$$j_\alpha(z, t) = -\beta D_\alpha n_\alpha(z, t) \bar{q}_\alpha e \frac{\partial \phi}{\partial z}(z, t) - D_\alpha \frac{\partial n_\alpha}{\partial z}(z, t), \quad (1)$$

where $\beta = 1/(k_B T)$ is the inverse temperature and $e$ the elementary charge. The first term represents a drift current, obtained as the product of a mobility $\beta D_\alpha$ (in agreement with Einstein's relation), a density and an electric force. The second term is a diffusion current, described by Fick's law. We assume a diagonal diffusivity tensor, meaning that the diffusion of a species is not influenced by the presence of other species. In general, the microscopic potential depends on the discrete positions of the ions. Within the Poisson-Nernst-Planck approximation, the potential is assumed to be a coarse-grained average of the microscopic potential, hence the mean-field nature of Eq. (1). In the static case, currents vanish and the Poisson-Boltzmann distribution is retrieved, with the average density proportional to the exponential of an average potential [27,28].

For currents $j_\alpha(z, t)$, the following exact continuity equation must hold, ensuring local ionic mass and charge conservation:

$$\frac{\partial n_\alpha}{\partial t}(z, t) = -\frac{\partial j_\alpha}{\partial z}(z, t). \quad (2)$$

Substituting Eq. (1) in Eq. (2) yields

$$\frac{\partial n_\alpha}{\partial t}(z, t) = D_\alpha \frac{\partial}{\partial z} \left( n_\alpha(z, t) \beta \bar{q}_\alpha e \frac{\partial \phi}{\partial z}(z, t) + \frac{\partial n_\alpha}{\partial z}(z, t) \right). \quad (3)$$

The potential is related to the density by the Poisson equation, also exact:

$$-\frac{\partial^2 \phi}{\partial z^2}(z, t) = \frac{\rho(z, t)}{\varepsilon_0 \varepsilon_r}, \quad (4)$$

where $\rho = q_+ e n_+ - q_- e n_-$ is the charge density, $\varepsilon_0$ the permittivity of vacuum, and $\varepsilon_r$ the relative permittivity. Equations (3) and (4) constitute the Poisson-Nernst-Planck theory.

The fact that ions are confined within the two slabs of the capacitor imposes the boundary condition of vanishing currents at $+L$ and $-L$. From the definition in Eq. (1):

$$\begin{aligned} -\beta n_\alpha(-L, t) \bar{q}_\alpha e \frac{\partial \phi}{\partial z}(-L, t) - \frac{\partial n_\alpha}{\partial z}(-L, t) &= 0, \\ -\beta n_\alpha(+L, t) \bar{q}_\alpha e \frac{\partial \phi}{\partial z}(+L, t) - \frac{\partial n_\alpha}{\partial z}(+L, t) &= 0. \end{aligned} \quad (5)$$

We will focus on the situation where the EDLC is at equilibrium, with zero applied potential for times $t < 0$, and is subject to a potential difference $2V_0$ for times $t > 0$. The electric potential must be continuous between slab and solution,

<table><caption>TABLE I. List of units of the nondimensional quantities used. Note that they do not form a coherent system of units: for instance, the units of distance, volumic density, and surface density are not related by simple powers.</caption><thead><tr><th>Observable</th><th>Symbol</th><th>Unit</th></tr></thead><tbody><tr><td>Time</td><td>t, τ<sub>i</sub></td><td>L/λ<sub>D</sub></td></tr><tr><td>Inverse time</td><td>s, s<sub>i</sub></td><td>D/λ<sub>D</sub></td></tr><tr><td>Distance</td><td>z</td><td>L</td></tr><tr><td>Volumic ion density</td><td>n<sub>±</sub></td><td>n<sub>0</sub></td></tr><tr><td>Volumic charge density</td><td>ρ</td><td>2en<sub>0</sub></td></tr><tr><td>Electric potential</td><td>φ, V</td><td>1/βe</td></tr><tr><td>Electric field</td><td>E</td><td>1/βeL</td></tr><tr><td>Surface density</td><td>σ</td><td>ε<sub>0</sub>ε<sub>r</sub>/βe<sup>2</sup>L = 1/4πl<sub>B</sub>L</td></tr></tbody></table>

so that

$$\phi(\pm L, t) = \pm V(t), \quad \text{with } V(t) = \begin{cases} 0 & \text{for } t < 0 \\ V_0 & \text{for } t \ge 0 \end{cases}. \quad (6)$$

The zero of the potential is arbitrary.

In the following, we solve the Poisson-Nernst-Planck equation analytically, when possible, and numerically, with a flux-conservative method whose details are given in the Appendix. The results of our numerical scheme were successfully compared to those of constant-potential lattice-Boltzmann electrokinetics [16], where lattice-Boltzmann is coupled with an iterative resolution of the Poisson equation.

### III. FULLY SYMMETRIC CASE: $q_+ = q_-$, $D_+ = D_-$

We focus first on the case of a symmetric binary electrolyte, whose species have valences $q_+ = q_- = 1$, and diffusivities $D_+ = D_- = D$. The problem of studying relaxation in this case has been previously addressed in Refs. [24] and [25], with particular focus on the linear regime. For the sake of completeness, we reobtain here some of the results from Ref. [25] following a simpler approach. In the linear regime, the relaxation is multiexponential. We then seek a clear characterization of the nonlinear regime, building on Ref. [24].

#### A. Linear regime

Linearizing Eqs. (3) around the initial densities $n_+^0 = n_-^0 = n_0$ and taking the difference of the two equations for the two species, yields, together with Eq. (4), the Debye-Falkenhagen equation for the charge density $\rho$:

$$\frac{\partial \rho}{\partial t} = D \left( \frac{\partial^2 \rho}{\partial z^2} - \lambda_D^{-2} \rho \right). \quad (7)$$

Here we defined the Debye length $\lambda_D$ as a constant in terms of the initial density $n_0$ by $\lambda_D^{-2} = 2n_0\beta e^2/(\varepsilon_r\varepsilon_0) = 8\pi l_B n_0$, where $l_B = \beta e^2/(4\pi\varepsilon_r\varepsilon_0)$ is the Bjerrum length.

We proceed by making the equations nondimensional, according to the mapping described in Table I. The system is completely described by the two dimensionless parameters,

$$\epsilon = \frac{\lambda_D}{L} \quad \text{and} \quad \nu = \beta e V_0, \quad (8)$$

that we will keep using throughout the rest of the paper. Note that $\epsilon$ can be, but need not be, a small quantity, while, as long
as we deal with the linear regime, $v$ has to be much smaller than unity.

In these units, Eq. (7) is readily rewritten as

$$ \epsilon \frac{\partial \rho}{\partial t} = \epsilon^2 \frac{\partial^2 \rho}{\partial z^2} - \rho \quad (9) $$

and the Poisson equation (4) reads

$$ -\epsilon^2 \frac{\partial^2 \phi}{\partial z^2} = \rho. \quad (10) $$

The boundary condition in Eq. (5) can be written in terms of the charge density upon linearization:

$$ -\frac{\partial \rho}{\partial z}(\pm 1, t) - \frac{\partial \phi}{\partial z}(\pm 1, t) = 0. \quad (11) $$

Finally, the boundary condition in Eq. (6) reads

$$ \phi(\pm 1, t) = \pm V(t). \quad (12) $$

We make the ansatz that the potential $\phi(z, t)$ and the electric charge density $\rho(z, t)$ relax to equilibrium as

$$ \phi(z, t) = v \frac{\sinh\left(\frac{z}{\epsilon}\right)}{\sinh\left(\frac{1}{\epsilon}\right)} + v \sum_{i=0}^{\infty} b_i(z) e^{s_i t} \quad (13) $$

$$ \rho(z, t) = -v \frac{\sinh\left(\frac{z}{\epsilon}\right)}{\sinh\left(\frac{1}{\epsilon}\right)} + v \sum_{i=0}^{\infty} B_i(z) e^{s_i t}, \quad (14) $$

for $t > 0$. The time-independent terms correspond to the only solutions of the steady-state Debye-Falkenhagen equation allowed by symmetry: Indeed, $\phi$ and $\rho$ must be odd with respect to $z$ and $v$. Additionally, by Eq. (10), we have $B_i(z) = -\epsilon^2 b_i''(z)$. We suppose $s_i < 0$, so that the characteristic relaxation times are $1/|s_i|$. Substituting Eq. (14) in Eq. (9) and enforcing an odd charge density gives

$$ B_i(z) = c_i \sinh \left( \frac{\sqrt{1 + \epsilon s_i}}{\epsilon} z \right). \quad (15) $$

This expression can be integrated to obtain the corresponding $b_i(z)$. Fixing the gauge $\phi(0) = 0$ and imposing boundary conditions Eq. (11), one finds

$$ b_i(z) = -c_i \left[ \frac{\sinh\left(\frac{\sqrt{1+\epsilon s_i}}{\epsilon} z\right)}{1+\epsilon s_i} + z \frac{s_i \cosh\left(\frac{\sqrt{1+\epsilon s_i}}{\epsilon}\right)}{\sqrt{1+\epsilon s_i}} \right]. \quad (16) $$

At this point, Eq. (12) gives, for any $i$, either $c_i = 0$ or

$$ 1 + s_i \sqrt{1 + \epsilon s_i} \coth \left( \frac{\sqrt{1 + \epsilon s_i}}{\epsilon} \right) = 0. \quad (17) $$

This equation is exactly equivalent to the one solved in Ref. [25] to obtain relaxation times. An alternative solution, consistent with our notation, was independently derived in Ref. [22]. The modes characterizing the linear response can also be considered in the frequency domain via the impedance of the cell (see, e.g., Refs. [29,30] for symmetric electrolytes).

For $\epsilon = \lambda_D/L \gg 1$, relaxation modes turn out to be of order $-\epsilon\pi^2(i + 1/2)^2$; the dominant (slowest) mode is $s_0 \sim -\epsilon\pi^2/4$, corresponding to a dimensional relaxation time $\tau_0 \sim 4L^2/(\pi^2 D)$. For $\epsilon = \lambda_D/L \ll 1$, the dominant relaxation rate results exactly from the solution of the equation

$$ \tanh \left( \frac{Z}{\epsilon} \right) = -\frac{Z}{\epsilon} (Z^2 - 1), \quad (18) $$

with $Z = \sqrt{1 + \epsilon s_0}$. In physical units (adopted from now on for characteristic timescales), this corresponds to a time

$$ \tau_0 \approx \frac{L \lambda_D}{D} - \frac{\lambda_D^2}{2D}. \quad (19) $$

The scaling $L\lambda_D/D$ was pointed out in Ref. [24] and previously reported in Refs. [31–33], while the second term represents the exact finite-double-layer correction as commented in Refs. [16,22,25]. The timescale from Eq. (19) is visible in Fig. 2(a), which shows numerical solutions of the ion density at contact with the electrode $\rho(-1, t)$ and of the electrode surface charge density $\sigma(t)$.

Note that assuming a monoexponential relaxation in Eqs. (13) and (14), for $\epsilon \ll 1$, still leads to a surprisingly good approximation for the dominant relaxation rate [22].

The exact relaxation profile can be obtained by an analysis of Eqs. (9)–(12) in the Laplace domain [22,25], allowing to identify $c_i$ and write

$$ \rho(z, t) = -v \frac{\sinh \frac{z}{\epsilon}}{\sinh \frac{1}{\epsilon}} + v \sum_{n=0}^{\infty} \frac{2(1 + \epsilon s_i)}{3 - s_i^2} \frac{\sinh \frac{z\sqrt{1+\epsilon s_i}}{\epsilon}}{\sinh \frac{\sqrt{1+\epsilon s_i}}{\epsilon}} e^{s_i t}. \quad (20) $$

A cumbersome computation shows that this equation, that we report for its perhaps more readable shape, is equivalent to Eq. (40) in Ref. [25].

As proposed in Ref. [24], when $\epsilon = \lambda_D/L \ll 1$, it is possible to establish an analogy between the EDLC and a simple RC circuit of the first order. We summarize this analogy here because it will prove useful in the following. Each double layer can be identified with a planar capacitor, opposing the charge of the electrode with an equal and opposite charge in the solution, distributed over a distance $\sim \lambda_D$. Its capacitance per unit surface reads

$$ C = \frac{\epsilon_0 \epsilon_r}{\lambda_D}. \quad (21) $$

Along the same spirit, the bulk, of length $2L - 2\lambda_D$, carries an electric resistance according to Drude’s model of conduction. For some applied potential difference $\Delta V$, the current density can be expressed as

$$ j = 2n_0 e^2 \beta D \frac{\Delta V}{2L - 2\lambda_D}. \quad (22) $$

The electric resistance of the bulk times unit surface (i.e., the inverse conductance per unit surface) is then

$$ R = \frac{\Delta V}{j} = \frac{L - \lambda_D}{n_0 e^2 \beta D}. \quad (23) $$

The circuit results as a series of a capacitance $C$, a resistance $R$, and a second capacitance $C$. Being the equivalent capacitance $C/2$, the characteristic time of the circuit turns out to be

$$ \tau = R \frac{C}{2} = \frac{L - \lambda_D}{2n_0 \frac{\beta e^2}{\epsilon_0 \epsilon_r} \lambda_D D} = \frac{L \lambda_D}{D} - \frac{\lambda_D^2}{D}. \quad (24) $$

This simple circuit analogy reproduces the dominant relaxation rate from Eq. (19)—the correction of order $\lambda_D^2/D$ is not exact but has the right scaling.
<figure><img src="image_4.png" /><figcaption>FIG. 2. Logarithmic plots of the quantities defined at the top left, where $\rho^{\pm L} = \rho(-L, t) = -\rho(L, t)$ is the charge density in the solution at contact with the electrode, $\sigma(t)$ is the surface charge density of the electrode, and $\rho_{eq}^{\pm L}$ and $\sigma_{eq}$ are their respective values at equilibrium. Here, time $t$ is in units of $L^2/D$, $\epsilon = 0.01$ and $v$ increases progressively, as indicated at the top right of each panel. When the curves shown are linear, their slopes correspond to exponential relaxation rates. $s_0$ corresponds to a relaxation time $\lambda_D L/D = 0.01 L^2/D$ and $s'$ of $L^2/(\pi^2 D) \simeq 0.10 L^2/D$. At $v = 0.001$ (a), only the double layer formation process is visible with its rate $s_0$ (unless too close to $t = 0$, where other rates $s_i$ can play a role). At $v = 0.1$ and 1 [(b) and (c)], the first slope corresponds to the buildup of the double layer, while the second slope reflects its reorganization as the bulk is depleted from ions. The system is in the purely nonlinear regime at $v = 1$ and the double layer builds up in a time $\tau_{PNL} \propto \cosh(v/2)$, as defined in Eq. (35) in Sec. III C. At $v = 10$ (d) the system is partially depleted at equilibrium, as discussed in Secs. III C and III D [see Fig. 3(c)] but continues to relax with a time $\tau_{PNL}$. At $v = 40$ and 100 [(e) and (f)], the system is fully depleted (Sec. III D): The early-time curves represent ion migration (for these parameters, the relaxation is not yet linear in time, but not perfectly exponential either); the late-time part represents relaxation of the counterionic double layers, in a time $\sim \mu_{nen}^2/D$. Insets show $|\rho^{\pm L} - \rho_{eq}^{\pm L}|/\rho_{eq}^{\pm L}$ and $|(\sigma - \sigma_{eq})/\sigma_{eq}|$ in linear scale.</figcaption></figure>

## B. Depletion

In the nonlinear regime, it is impossible to write a Debye-Falkenhagen equation for the charge density and it is harder to make analytical predictions on the dynamics. The first clearly visible nonlinear effect is depletion. We define depletion as the decrease in the bulk population from the initial value $n_0$ to some smaller final value. As $v$ increases, more ions are attracted to the oppositely charged electrode than are repelled by the like-charged one. Since the total number of ions (the integral of the ion densities) must be conserved, this calls for a decrease in the bulk density. Indeed, in the symmetric electrolyte case, Eq. (3) is not invariant under the transformation $(n_{\pm} - n_0 \rightarrow n_0 - n_{\pm}, z \rightarrow -z, \bar{q}_{\pm} \rightarrow \bar{q}_{\mp})$—only valid for asymptotically small voltages—meaning that exchanging the two ion species cannot be reduced to a simple flip of their excess densities about their initial values. More intuitively, if the applied potential is sufficiently high, then all positive ions condense on the negative electrode and all negative ions condense on the positive electrode: In this extreme case the bulk of the capacitor is basically empty and we have full depletion, as described in Sec. III D. In this section, we focus on depletion in the linear regime $v \ll 1$, as a weakly nonlinear effect. The phenomenon has a clear signature on the relaxation dynamics and it introduces a purely diffusive timescale,

that, at least in the $\epsilon \ll 1$ range, is slower than the dominant relaxation time $\tau_0$ from Eq. (19), related to the electric double layer formation.

To quantify depletion, we introduce a quantity, inspired by the Dukhin number [34,35], that quantifies the fraction of ions that remains in the bulk after relaxation:

$$ \mathrm{Du}_n = \frac{n_+(0, \infty)}{n_0} = \frac{n_-(0, \infty)}{n_0}. \qquad (25) $$

This quantity is 1 when there is no depletion and 0 when ions are completely absent at the center of the capacitor. We say that the system is weakly depleted when $0.9 < \mathrm{Du}_n < 1$, depleted when $\mathrm{Du}_n < 0.9$ and fully depleted when $\mathrm{Du}_n \simeq 0$.

We perform a weakly nonlinear analysis of the Poisson-Nernst-Planck equations. Indeed, in the linear regime, the magnitude of the corrections to $\rho$ and $\phi$ are orders of magnitude below the dominant process (the formation of the two double layers), but a new depletion-related timescale appears, slower than the previously mentioned ones. We expand ionic densities, charge density and potential—that we already know up to linear order in $v$—and consider small terms of order $v^2$ and $v^3$. Since $\rho$ and $\phi$ must be odd functions of $v$, we write

$$ \begin{aligned} n_{\pm}(z, t) = & n_{\pm}^{(0)} + v n_{\pm}^{(1)}(z, t) + v^2 n_{\pm}^{(2)}(z, t) \\ & + v^3 n_{\pm}^{(3)}(z, t) + O(v^4) \end{aligned} \qquad (26) $$
$$ \rho(z, t) = v\rho^{(1)}(z, t) + v^3\rho^{(3)}(z, t) + O(v^5) \quad (27) $$

$$ \phi(z, t) = v\phi^{(1)}(z, t) + v^3\phi^{(3)}(z, t) + O(v^5), \quad (28) $$

where all quantities are dimensionless and $\rho = \frac{n_+ - n_-}{2}$ according to the units in Table I.

Introducing the above expansions in the Poisson-Nernst-Planck equation (3) and collecting $O(v^0)$ terms, we find that at order $0, n_{\pm}^{(0)} = n_0$. From $O(v^1)$ terms, we find that at linear order, $\phi^{(1)}(z, t)$ is given by Eq. (13) and $\rho^{(1)}(z, t) = \pm n_{\pm}^{(1)}(z, t)$ by Eq. (14). Then, introducing these results into the equation for $O(v^2)$ terms, we obtain

$$ \frac{\partial n_{\pm}^{(2)}}{\partial t} = \epsilon \frac{\partial^2 n_{\pm}^{(2)}}{\partial z^2} - \frac{1}{\epsilon} \frac{\cosh \frac{2z}{\epsilon}}{\sinh^2 \frac{1}{\epsilon}} + \sum_{i=0}^{\infty} e^{s_i t} f_i(z) + \sum_{i=0}^{\infty} \sum_{j=0}^{\infty} e^{(s_j+s_i)t} g_{ij}(z) \quad (29) $$

where the exact form of the functions $f_i$ and $g_{ij}$, which can be expressed in terms of the $B_i$ and $b_i$ given in Eqs. (15) and (16), are in fact not relevant for the present analysis. Indeed, we look for relaxation modes slower than the purely linear one [$s_0$, as determined by Eq. (17)] for large $t$, so that the last two terms can be neglected. Imposing mass conservation in the form $\int_{-1}^{1} n_{\pm}^{(2)}(z, t) dz = 0$ and requiring $n_{\pm}^{(2)}$ to be even in $z$ to respect symmetry, one obtains

$$ n_{\pm}^{(2)}(z, t) = -\frac{\epsilon}{4} \coth \frac{1}{\epsilon} + \frac{\cosh \frac{2z}{\epsilon}}{4 \sinh^2 \frac{1}{\epsilon}} + A \cos(\pi z) e^{s't} + o(e^{s't}), \quad (30) $$

with A constant and $s' = -\pi^2\epsilon$. It can be checked that the first two terms in the right-hand side correspond to a particular solution of Eq. (29), while the last one is the standard solution for the corresponding homogeneous diffusion equation satisfying the symmetry and boundary conditions. In physical units, this corresponds to a diffusive timescale,

$$ \tau' = \frac{L^2}{\pi^2 D}. \quad (31) $$

This second-order correction to the ionic density represents depletion in the $\epsilon \ll 1$ regime (where $s'$ is indeed slower than $s_0$ and (30) makes sense). In this case, the equilibrium density in the bulk (i.e., at $|z| \ll 1 - \epsilon$), which in our units is nothing but $Du_n$, reads to second order

$$ n_{\pm}(z, \infty) \approx n_{\pm}(0, \infty) = Du_n = 1 - v^2 \frac{\epsilon}{4} \coth \frac{1}{\epsilon}. \quad (32) $$

Figure 3(a) compares numerical results with our analytical approximation, showing that Eq. (30) well predicts late-time density profiles as a function of $z$ and $t$. In this figure, similar for this symmetric case to Figs. 8(a), 9, and 10 (top right) in Ref. [24], the time-damped cosine shape is visible outside the double layer (light green curves). The depletion phenomenon emerges as a (neutral) mass diffusion of both species from the center of the EDLC toward the boundaries of the bulk region: Ions to constitute the double layer are initially recruited from the regions close to the electrode, leaving a nonuniform mass distribution in the bulk [light blue curves in Fig. 3(a)]. The

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 3. (a) Mass as a function of $z$ (in units of Table I), at different times (here in units of $L^2/D$), for $\epsilon = 0.01$ and $v = 2$. In these units, the electric double layer formation occurs on a timescale $\sim 0.02 L^2/D$. At much shorter times (dark blue) the system has not moved yet. At later times (lighter greens) the double layer has formed already and the mass diffusion process manifests itself, with a relaxation timescale $\tau' \simeq 0.1 L^2/D$: the sinusoid of Eq. (30) is visible. At times much larger than $\tau'$ (dark green) the system is at equilibrium at a new value of bulk density, predicted by Eq. (32). Dotted lines indicate the analytical predictions for later times as per Eq. (30), where the parameter A was set to 0.006 for all curves.</figcaption></figure>
corresponding potential across the capacitor is shown as a function of time in Fig. 3(b).

The mass imbalance then evens out with a purely diffusive process on a length scale $L$ rather than $2L$ (ions move from the center to the borders and not from one border to another). This explains the absence of the diffusive mode $4L^2/(\pi^2 D)$ for symmetric electrolytes, otherwise legitimate. As shown in Fig. 3(c), Eq. (32) predicts values of $Du_n$ in the weakly nonlinear regime.

Finally, the rate $s'$ can be proven to emerge naturally also in $\rho^{(3)}(z, t)$ and $\phi^{(3)}(z, t)$. This is why the relaxation rate $s'$ is visible in late-time profiles of $\rho(-1, t)$ and $\sigma(t)$, plotted in Figs. 2(b) and 2(c), for $\epsilon = 0.01$.

To summarize, the linear regime features two processes: double layer formation and, at higher order, depletion. The double layer forms at a dominant rate $s_0$. During this process the bulk stays electroneutral but becomes inhomogeneous as for mass density: The bulk region closer to the electrodes becomes less populated than the central region. A slower process, at least for $\epsilon < \pi^{-2} \simeq 0.1$, then onsets with rate $s'$: This is a mere diffusion of neutral excess mass within the bulk, with positive and negative ions moving together from the center toward the double layer boundaries. The diffusive depletion process is asymptotically absent at $v \to 0$, but its relative importance compared to the double layer formation process grows with $v$ (see the upwards shift of late-time lines in Figs. 2(b) and 2(c)).

A last word of caution concerns the quantity $\lambda_D$ defined after Eq. (7), and its dimensionless equivalent $\epsilon = \lambda_D/L$. This “reservoir” Debye length is defined in terms of the initial concentration $n_0$, but it is important to note that the physically relevant Debye length in the capacitor, defined for instance in terms of the depleted midpoint salt concentration, can be much larger than $\lambda_D$.

### C. Purely nonlinear regime

Figure 3(c) shows that for $\epsilon \gg 1$ depletion occurs as soon as $v \gtrsim 1$, so it is impossible to tell apart the purely nonlinear effects, the ones that would appear even in a grand-canonical formulation of the problem, from the strictly canonical effects of depletion. In this case, we say that no purely nonlinear regime exists: Upon increasing $v$, the system goes directly from the linear regime to a depleted nonlinear regime, which we will describe in Sec. III D. On the contrary, for $\epsilon \ll 1$, more ions are available and depletion is only observed at voltages significantly higher than 1: This makes purely nonlinear effects visible at intermediate voltages, starting from $v \simeq 1$. These consist in a clear asymmetry, for a given species,

between left and right double layer at equilibrium; however, the total number of ions is much larger than the number of ions involved in the double layers, so that the bulk population stays almost unaffected (no depletion).

In this purely nonlinear regime, the linear analysis of the relaxation times from Sec. III A is not valid and the double layer formation does not happen anymore at a rate $s_0$. However, it is possible to understand the change in the rate of formation of the electric double layer by using the Grahame equation [27,28]. In the units of Table I, this equation reads

$$|\sigma(\infty)| = \frac{2}{\epsilon} \sinh\left(\frac{v}{2}\right) \quad (33)$$

and is exact in the $\epsilon \to 0$ limit, even for $v > 1$. The following differential capacitance emerges:

$$\frac{\partial |\sigma(\infty)|}{\partial v} = \frac{1}{\epsilon} \cosh\left(\frac{v}{2}\right). \quad (34)$$

If this capacitance is used in the circuit model, replacing the one of Eq. (21), then the relaxation time in Eq. (24) becomes

$$\tau_{\text{PNL}} = \frac{\lambda_D L}{D} \cosh\left(\frac{v}{2}\right), \quad (35)$$

where PNL denotes purely nonlinear. The same scaling appears already in Refs. [24,36] and is in quantitative agreement with our numerical calculations (see Fig. 2(d) and Fig. 2 in Ref. [26], where timescales extracted from numerical simulations are summarized).

The Grahame Eq. (33) also allows us to estimate for what values of $v$ and $\epsilon$ depletion starts to become relevant, determining a boundary between the purely nonlinear and the depleted regimes for small $\epsilon$. We introduce the Dukhin number

$$Du = \frac{|\sigma(\infty)|}{2n_0 L}, \quad (36)$$

giving information on the maximum surface charge neutralisable by the ions [34,35]. We define a system as depleted when $Du \simeq 0.1$, in a manner conceptually equivalent to the previously employed $Du_n < 0.9$. Using Eq. (33), this condition reads

$$Du = 2\epsilon \sinh\left(\frac{v}{2}\right) \simeq 0.1. \quad (37)$$

This relation can be used to identify the limit between purely nonlinear and fully depleted nonlinear regime in a $(v, \epsilon)$ diagram, as we do in Fig. 3 of Ref. [26], showing that it is consistent with numerical results.

### D. Fully depleted nonlinear regimes

Upon increasing $v$ at $\epsilon \ll 1$, depletion becomes more and more important. The Dukhin number $Du_n$ correspondingly decreases [Fig. 3(c)] and the late-time line representing depletion in Figs. 2(b) and 2(c) gradually shifts upwards. At the same time, still for $\epsilon \ll 1$, the timescale related to charging increases as per Eq. (35). In terms of Fig. 2, the slope of the early-time curve gradually decreases as $v$ increases. Eventually the two processes (double layer charging and depletion) become indistinguishable and the system gradually enters the fully depleted nonlinear regime. A similar thing happens at

<footnote>(b) Potential as a function of $z$, at different times (here in units of $L^2/D$), for $\epsilon = 0.01$ and $v = 2$. The inset is a zoom close to the left electrode, while the arrow indicates the point of null potential (see Sec. IV B), which in the present case is $z_0 = 0$. (c) Points show $Du_n$ as a function of $v$, for different values of $\epsilon$, as extracted from numerical solutions of the Poisson-Nernst-Planck equation. Smaller $Du_n$ represent stronger depletion. Solid lines are a guide to the eye. The curves at $\epsilon = 100$ and 1000, not shown here, coincide with the $\epsilon = 10$ curve. In dotted lines, the prediction from Eq. (32).</footnote>
$\epsilon \gg 1$, where depletion coincides with nonlinearity, as discussed above. The transition to this new fully depleted regime is not abrupt and defines a band on the $(\epsilon, v)$ diagram where depletion is only partial (say, $0.1 < Du_n < 0.9$ or $Du \sim 1$). Such a band is relatively narrow: Figure 3(c) shows that the passage from $Du_n = 1$ to $Du_n = 0$ happens in no more than a decade of $v$, for all probed values of $\epsilon$ (see also Fig. 3 in Ref. [26]).

For parameters at which depletion is fully achieved, the “bulk” ionic concentration at equilibrium $n_\pm(0, \infty)$ can be orders of magnitude lower than the initial one $n_0$. In this case, all positive ions are concentrated in proximity of the negative electrodes and vice versa. At long times, the double layer contains only counterions, so that the equilibrium distribution is governed by a sort of Gouy-Chapman length for either subsystem made of one electrode and its counterions. Such subsystem is in general not electroneutral (it cannot be so if $Du > 1$): the Poisson equation shows that its Gouy-Chapman length reads, in physical units, $\mu_{nen} = (2\pi qL_B\sigma_{res})^{-1}$, where $\sigma_{res} = |\sigma(\infty)| - 2n_0L$ is the unscreened residual part of the electrodes’ surface charge (where “nen” denotes nonelectroneutral and “res” denotes residual). After the counterionic double layers have formed, it is reasonable to think that the last dynamic phenomenon to happen is a rearrangement of each double layer on the smallest length scale available, $\mu_{nen}$; this would correspond to a relaxation time $\mu_{nen}^2/D$. In the following we show how to predict $\sigma(\infty)$ and thus $\mu_{nen}$, and what dynamics leads to the formation of the counterionic double layers.

We take a look at the large voltage asymptotics to simplify the problem. We suppose $v$ to be so high that even when equilibrium is reached, ions have negligible impact on the potential profile between the electrodes: the electric field, at any time, is approximately $V_0/L$ everywhere, except possibly within a very small distance from the electrodes. We will verify later when this assumption is correct. As soon as the power source is switched on, ions, of valence $q = 1$, move by electric drag with a constant velocity $v = \beta DeV_0/L = vD/L$, directed toward the oppositely charged electrode. Once the ions reach the electrode, we suppose the latter to be so highly charged that a very thin double layer is formed, of negligible size compared to $L$. If the number of ions reaching the double layer is somehow proportional to the density at contact, then we expect $n_\pm(\mp L, t)$, $\rho(\pm L, t)$ and, more rigorously, the double layer charge $\int_{EDL} \rho(z, t) dz$ to grow linearly with time, approaching their final value in a time $t^*$. Such time is easily computed, in physical units, as the time needed for the furthermost ion to reach the oppositely charged electrode, a distance $2L$ away:

$$t^* = \frac{2L}{v} = \frac{2L^2}{Dv}. \qquad (38)$$

Figure 4 confirms that $\rho(-L, t)$ is linear in time from $t = 0$ until $t^*$.

This allows us to compute the time evolution of $\sigma$. By definition, integrating the electric field to get the electric potential across a half-capacitor (from $-L$ to $0$) for $0 < t < t^*$, one

<figure><img src="image_8.png" /><figcaption>FIG. 4. Ratios of $|\sigma(t)|$ (black) and $|\rho(\pm L, t)|$ (gray) versus their equilibrium values, as a function of time. Here $\epsilon = 0.1$ and $v = 200$. The blue dashed lines represent the ideal linear evolution of $\rho$ and the theoretical prediction for $\sigma$ from Eqs. (38)–(41).</figcaption></figure>

must obtain the potential difference $V_0$. In physical units:

$$V_0 = \frac{e|\sigma(t)|}{\epsilon_0\epsilon_r}L - \frac{e}{\epsilon_0\epsilon_r}\int_{-L}^{0} dz \int_{-L}^{z} dz' \times \underbrace{\left(2n_0L\frac{t}{t^*}\delta(z' + L) + n_0h(z', t)\right)}_{\rho(z',t)/e} \quad (39)$$

with

$$h(z', t) = \begin{cases} I_{(-L, -L+vt)}(z') & \text{if } 0 < t < \frac{t^*}{2} \\ I_{(-L, L-vt)}(z') & \text{if } \frac{t^*}{2} < t < t^*, \end{cases} \quad (40)$$

where $I_{(z_1, z_2)}(z')$ is the gate function, equal to 1 if $z' \in (z_1, z_2)$ and to 0 otherwise. In Eq. (39), the inner integral represents the electric field in $z$ due to the ions between $-L$ and $z$. The first term in parenthesis represents positive ions that have adhered to the wall in $z = -L$ at time $t$; we make the drastic choice of a $\delta$-function distribution, but this is not crucial, as long as ions stay confined within a length $\ll L$. The second term in parenthesis is the charge density at point $z'$ in the rest of the solution, whose functional shape Eq. (40) is determined by the following observation. For $t < t^*/2$, the nonzero charge density is due to negative ions leaving the region $(-L, 0)$ at velocity $v$ from left to right, while the concentration of positive ions stays constant in this half of the capacitor (though not at $-L$); for $t > t^*/2$ all negative ions have left the region and the last positive ions approach the negative electrode with velocity $v$ from right to left. Solving Eq. (39) for $0 < t < t^*$ gives, in dimensionless units,

$$|\sigma(t)| = v + \frac{1}{\epsilon^2}\left(\frac{2t}{t^*} - \frac{t^2}{t^{*2}}\right). \quad (41)$$

The electrodes’ charge when the double layer formation is concluded ($t = t^*$) is then $|\sigma(t^*)| = v + \epsilon^{-2}$. This is of straightforward interpretation, since, in our units of surface charge density, $2n_0L$ reads $\epsilon^{-2}$: The surface charge developed by the electrodes to maintain a field $v$ across the EDLC is indeed $v + \epsilon^{-2}$, because exactly $\epsilon^{-2}$ of its charge is screened by the counterions stuck at the electrode. Agreement of Eq. (41) with numerical data is shown in Fig. 4. This parabolic time dependence of the surface charge density is also consistent with
the linear current derived in Ref. [37] for the fully symmetric case considered in this section.

For times $t > t^*$, we expect to observe the fast relaxation mentioned at the beginning of this section, on a timescale $\mu_{\text{nen}}^2/D$. Now the residual surface charge to be used in $\mu_{\text{nen}}$ is nothing but $\sigma_{\text{res}} \simeq |\sigma(t^*)| - \epsilon^{-2} = v$. This corresponds to a relaxation time

$$ \tau = \frac{4L^2}{Dv^2}, \quad (42) $$

which is indeed observed at sufficiently high voltages, for all values of $\epsilon$ (see Fig. 2 of Ref. [26]). This exponential relaxation and corresponding timescale were not reported in Ref. [37]. We note that in the latter reference, the authors argued for a power-law current with exponent $-3/4$, which seem supported by the experimental results of Ref. [38]. Here we did not attempt to give scaling laws: The system transitions to the unscreened regime in a complex fashion, even more complex for asymmetric electrolytes (see Secs. IV and V), and, at least within the range of parameters that we can probe, it is hard to find scalings that work for the whole regime or to neatly separate sub-regimes.

We made so far the crucial assumption that ions do not affect the linear potential profile through almost all the capacitor. A necessary condition for this to happen is that the electrodes' surface charge be larger than the integrated density of ions, meaning $|\sigma(t)|/(2n_0L) > 1$. Imposing this to Eq. (41), we obtain

$$ v > \frac{1}{\epsilon^2}. \quad (43) $$

This defines the unscreened fully depleted nonlinear regime and was used to delimit such a region in Fig. 3 in Ref. [26].

There is a transition region from the purely nonlinear regime to the (unscreened) fully depleted regime we just described: The transition region is the area that comprises approximately between the curve (37) and the curve $v = 1/\epsilon^2$ (see Fig. 3 in Ref. [26]). We call this the partially screened fully depleted nonlinear regime. In such region, the relaxation process gradually changes from exponential, with a well-defined relaxation time given by Eq. (35), to linear in time, as in Eq. (38), due to depletion. At equilibrium, the electric field generated by the electrodes is reduced (partially screened) to a fraction of $v$ in the fully depleted zone at the center of the capacitor: double layers are indeed sufficiently populated to screen a non-negligible fraction of it. In addition, the electric field, as a function of $z$, changes with time, as more and more ions reach the electrodes. As a consequence, analytical examination is hard and the functional time dependence of relaxation processes had to be retrieved numerically. Numerical results [Figs. 2(e) and 2(f)] show anyway the presence of a late-time relaxation on a fast scale $\sim \mu_{\text{nen}}^2/D$, that, as $v$ increases, converges to the one (42) predicted for the unscreened fully depleted nonlinear regime (see Fig. 2 in Ref. [26]).

Results from this section are summarized in Fig. 3 in Ref. [26].

## IV. PARTIALLY ASYMMETRIC CASE: $q_+ \neq q_-$, $D_+ = D_-$

### A. Linear regime and depletion

We now tackle the case of asymmetric valences: $q_+ \neq q_-$. This requires a redefinition of our proxy for ion concentration

$\lambda_D$, so that in the linear regime this equals the Debye length [27]:

$$ \lambda_D^{-2} = 4\pi l_B (q_+^2 n_+^0 + q_-^2 n_-^0), \quad (44) $$

where $n_+^0$ and $n_-^0$ are the initial concentrations of positive and negative ions.

The linear regime exhibits the same dynamics as the symmetric valence case, as shown by a linearization of the Poisson-Nernst-Plank equation (3). For $v > 0$, depletion occurs analogously to Sec. III B, with one difference: With asymmetric valences, there is no reason to expect $n_+^{(2)}(z, t)$ to be even functions of $z$, nor to have any particular symmetry. This invalidates the reason why depletion-related diffusion took place on a length $L$ rather than $2L$ in a system with equal valences [Eq. (31)]. In an asymmetric-valence system, the relaxation mode corresponding to a time

$$ \tau' = \frac{(2L)^2}{\pi^2 D} \quad (45) $$

is then permitted. This characteristic time for depletion is observed numerically (see Fig. 5) at long times; for some parameter values (e.g., 1:10 case, $\lambda_D/L = 0.01$, $v = 0.1$), we observe, however, that after the end of the double layer relaxation and before the onset of this slower depletion mode, the faster depletion mode $L^2/(\pi^2 D)$ is still visible.

### B. Purely nonlinear regime

The nonlinear regime also deviates from the symmetric-valence case: This is shown again in Fig. 5 (with numerical results shown as symbols), which has to be compared to Fig. 2 in Ref. [26]. The exponential increase of the relaxation time with the applied voltage at small $\epsilon$, which in the symmetric-valence case is explained by Eq. (35), is not valid anymore in the general $q_+ : q_-$ case. This is particularly evident in the 1:10 case (Fig. 5, right panel), where the relaxation time appears to be nonmonotonic with $v$: an initial slight decrease, absent in the 1:1 case, is followed by an increase, steeper than in the 1:1 case. An analytical estimate for this curve can be obtained by a procedure analogous to the one leading to Eq. (35), making use of the simple RC circuit analogy. In the asymmetric valence situation, the Grahame equation (33), used to compute the capacitance, can be rewritten as per Ref. [39]. The equilibrium (infinite time) charge density of the negative electrode, in the limit where the two double layers are completely separated ($\epsilon \to 0$) reads:

$$ |\sigma(\infty)| = \sqrt{\frac{n_0}{2\pi l_B}} (q_-e^{-q_+\psi_-} + q_+e^{q_-\psi_-} - (q_- + q_+))^{1/2}, \quad (46) $$

where $\psi_- = \phi(-L) - \phi(0) < 0$ is the potential of the negative electrode compared to the neutral bulk, in units of $(\beta e)^{-1}$, and $n_0 = n_+^0/q_- = n_-^0/q_+$ is the salt concentration. If the potential on the right electrode is $\psi_+ = \psi_- + 2v > 0$, then the surface charge on the positive electrode reads:

$$ |\sigma(\infty)| = \sqrt{\frac{n_0}{2\pi l_B}} (q_+e^{q_-(\psi_-+2v)} + q_-e^{-q_+(\psi_-+2v)} - (q_+ + q_-))^{1/2}. \quad (47) $$

This corresponds to Eq. (46), with the changes $q_\pm \to q_\mp$ and $\psi_- \to -\psi_+$. The charge densities on the two planes must
<figure><img src="image_9.png" /><figcaption>FIG. 5. Relaxation times $\tau$ as a function of $v$, for the 1:2 case (left) and the 1:10 case (right); different colors correspond to different values of $\epsilon = \lambda_D/L$. Times are extracted from linear fits in logarithmic scale of $\sigma$, similarly to what shown in Fig. 2 for the symmetric case. When two exponential relaxation times are visible from numerical data, one at short times and one at long times, they are both represented here (if the same $v$ and $\epsilon$ correspond to two symbols of the same color, the solid one represents early-time and the empty one late-time relaxation). The dense dotted lines represent the purely nonlinear time obtained from Eqs. (51) and (52) and preceding ones; for comparison, we also plot the time given by Eq. (35), relevant for the symmetric case and shown with sparser dots for $\epsilon = 0.01$. The gray and black lines represent the times given by Eqs. (55) and (59), respectively. We only present numerical results for exponential relaxation times here, so no data is shown for the time-linear ion migration (gray line), which is described later in Fig. 8.</figcaption></figure>

be equal in absolute value because of global neutrality. This allows to equate Eqs. (46) and (47) to obtain $\psi_-$.

$$ \psi_{-} = \frac{1}{q_{+} + q_{-}} \left( -2vq_{-} + \ln \frac{q_{-} \sum_{m=0}^{q_{+}-1} e^{-2vm}}{q_{+} \sum_{m=0}^{q_{-}-1} e^{-2vm}} \right). \quad (48) $$

Note that $\psi_{-} = -v + O(v^2)$ as $v \to 0$, as it should, and $\psi_{-} = -\frac{q_{-}}{q_{+}+q_{-}}2v + O(1)$ as $v \to \infty$.

The differential capacitance per unit surface of the negative electrode is the derivative of Eq. (46) with respect to $|\psi_-|$, computed at the point determined by Eq. (48):

$$ C_{-} = \frac{\varepsilon_{0}\varepsilon_{r}}{\lambda_{D}} \sqrt{\frac{q_{+}q_{-}}{2(q_{+} + q_{-})}} \frac{-e^{q_{-}\psi_{-}} + e^{-q_{+}\psi_{-}}}{\sqrt{q_{-}e^{-q_{+}\psi_{-}} + q_{+}e^{q_{-}\psi_{-}} - (q_{+} + q_{-})}} \quad (49) $$

Analogously, for the positive electrode we have:

$$ C_{+} = \frac{\varepsilon_{0}\varepsilon_{r}}{\lambda_{D}} \sqrt{\frac{q_{+}q_{-}}{2(q_{+} + q_{-})}} \times \frac{e^{q_{-}(\psi_{-} + 2v)} - e^{-q_{+}(\psi_{-} + 2v)}}{\sqrt{q_{-}e^{-q_{+}(\psi_{-} + 2v)} + q_{+}e^{q_{-}(\psi_{-} + 2v)} - (q_{+} + q_{-})}} \quad (50) $$

Interestingly, for $q_- / q_+ \gtrsim 3.18$, $C_-$ exhibits a nonmonotonic behavior, reflected in the total capacitance per unit surface,

$$ C = \left( \frac{1}{C_{-}} + \frac{1}{C_{+}} \right)^{-1}, \quad (51) $$

which is represented in Fig. 6.

The resistance times unit surface of the equivalent RC circuit for the general $q_+ : q_-$ case, in analogy with Eq. (23), reads

$$ R = \frac{2L - 2\lambda_D}{\beta e^2 n_0 q_+ q_- (q_+ + q_-) D}, \quad (52) $$

which can be rewritten in terms of $\lambda_D$ using the fact that $q_+ n_+^0 = q_- n_-^0 = n_0 q_+ q_-$. Multiplying this resistance by the capacitance from Eqs. (51) and preceding ones, one obtains an equivalent-RC-circuit timescale that is of order $L\lambda_D/D$ at small $v$ and diverges exponentially with $v$. This time is represented by the densely dotted curves in Fig. 5 and well captures the numerical results, both in the nonmonotonic behavior and in the steeper ascent compared to the 1:1 case.

The loss of left/right symmetry in the composition of the electric double layers and in the potential profile comes with a shift in the zero of the potential, i.e., the point $z_0$ on the z axis where the solution is neutral. This point is at the center of the capacitor ($z_0 = 0$) when $v \to 0$; upon increasing $v$ it moves toward the electrode producing the smaller voltage drop, i.e., the electrode of the same sign as the least charged species. This is shown in Fig. 7. The shift appears linear in $v$ at small applied voltages. At high voltages (where depletion is present, though), it seems to saturate at $(q_- - q_+)/((q_- + q_+))$, thus

<figure><img src="image_10.png" /><figcaption>FIG. 6. Total capacitance per unit surface $C$, as a function of the applied voltage $v$, for $q_+ = 1$ and $q_- = 1, 2, 3, 5, 10$. $C$ is computed from Eqs. (48)–(51). The vertical scale is logarithmic and the slope of the oblique asymptotes, for $v \to \infty$, is $q_+q_-/(q_+ + q_-)$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 7. Point of null potential z₀, as a function of applied voltage, in logarithmic scale. z₀ is retrieved numerically as the point where the charge density is zero. For both the 1:2 and the 1:10 cases, different values of ε (10, 1, 0.1) are shown. Circles represent results from a numerical solution of the nonlinear Poisson-Boltzmann equation, for a grand-canonical system with Debye length ε<sub>GC</sub>L = 0.1L; the arrows highlight that depletion, absent in the grand-canonical ensemble, pushes z₀ toward the positive electrode. The dotted gray line is a guide with slope 1. Numerical errors on z₀/L are of order 10⁻³.</figcaption></figure>

dividing the cell into two parts, whose sizes are in the same ratio q<sub>-</sub>/q<sub>+</sub> as the potential drops across the two electrodes.

Finally, to determine the boundary between the purely non-linear and the fully depleted regime, we extend our definition of the Dukhin number in Eq. (36) to the asymmetric valences case, where to a salt concentration n₀ corresponds a maximum charge of q<sub>+</sub>n<sub>+</sub><sup>0</sup> = q<sub>-</sub>n<sub>-</sub><sup>0</sup> = q<sub>+</sub>q<sub>-</sub>n₀:

$$
\mathrm{Du} = \frac{|\sigma(\infty)|}{2Lq_{+}q_{-}n_{0}}. \qquad (53)
$$

Using Eqs. (46) or (47), the condition Du ≃ 0.1 reads

$$
\epsilon \sqrt{\frac{q_+ + q_-}{2q_+ q_-}} (q_- e^{-q_+ \psi_-} + q_+ e^{q_- \psi_-} - (q_+ + q_-))^{1/2} \approx 0.1, \qquad (54)
$$

to be complemented with Eq. (48). Note that Eq. (54) reduces
to Eq. (37) when q<sub>+</sub> = q<sub>-</sub>, as it should.

C. Fully depleted nonlinear regime

The unscreened fully depleted nonlinear regime, for
asymptotically high applied voltages, is now characterized by
the fact that the two ionic species have different velocities: A
constant electric field drags the negative species q-/q+ times
faster (or slower) than the positive one. The charge in the two
Debye layers grows linearly in time, like in the symmetric
case, but arrives at its final value at two different times. Com-
pared to Fig. 4, this means that the curve ρ(-L, t)/ρ(-L, ∞)
has a different slope than the curve ρ(+L, t)/ρ(+L, ∞),
which in turn translates to an abrupt change in curvature of
the σ(t) curve.

In the quantitative analysis of this regime, we proceed as in
the symmetric case (Sec. III D). As said, the velocities of the
two species are now different: v± = Dq±eβV₀/L, assuming
again very large electric field and neglecting the effect of
the electrolyte distribution on the velocities. If each species
travels at its (constant) velocity, then the times at which the
furthermost ion of each species has reached the oppositely
charged electrode are as follows:

$$
t_{\pm}^{*} = \frac{2L}{v_{\pm}} = \frac{2L^2}{Dq_{\pm}v}. \qquad (55)
$$

After a time t<sub>+</sub><sup>*</sup>, the total number of positive charges in the
system 2n'<sub>0</sub>L (we define n'<sub>0</sub> = q<sub>+</sub>n<sub>+</sub><sup>0</sup> = q<sub>-</sub>n<sub>-</sub><sup>0</sup> = q<sub>+</sub>q<sub>-</sub>n<sub>0</sub>) are
adsorbed on the negative electrode, and after a time t<sub>-</sub><sup>*</sup>, the
same number of negative charges are adsorbed on the positive
electrode. As in Sec. III D, we assume that the electrolyte
density at contact with the electrodes grows linearly in time
and is localized in a region much smaller than the system size.
For simplicity, we take q<sub>-</sub> > q<sub>+</sub>, i.e., t<sub>-</sub><sup>*</sup> < t<sub>+</sub><sup>*</sup>. The analogous
of Eq. (39) for the asymmetric case, for 0 < t < t<sub>+</sub><sup>*</sup>, is as
follows:

$$
2V_0 = \frac{2e|\sigma(t)|}{\epsilon_0 \epsilon_r} L - \frac{e}{\epsilon_0 \epsilon_r} \int_{-L}^{L} dz \int_{-L}^{z} dz' \underbrace{\left( 2n'_0 L \frac{t}{t_+^*} \delta(z' + L) - 2n'_0 L \min\left(\frac{t}{t_-^*}, 1\right) \delta(z' - L) + n'_0 h(z', t) \right)}_{\rho(z',t)/e}, \quad (56)
$$

with

$$
h(z', t) = \begin{cases} I_{(-L, -L+\nu_{-t})}(z') - I_{(L-\nu_{+t}, L)}(z') & \text{if } 0 < t < t_{\text{meet}} \\ I_{(-L, L-\nu_{+t})}(z') - I_{(-L+\nu_{-t}, L)}(z') & \text{if } t_{\text{meet}} < t < t_{-}^{*} \\ I_{(-L, L-\nu_{+t})}(z') & \text{if } t_{-}^{*} < t < t_{+}^{*} \end{cases} \tag{57}
$$

Again, I(z₁, z₂)(z') is the gate function. In Eq. (56), the first term
in parenthesis represents positive ions that have adhered to
the wall in z = -L at time t; the second term in parenthe-

sis represents negative ions that have adhered to the wall in
z = L at time t (we write it for clarity, but this term does not
contribute to the integral, as the presence of adsorbed ions is
already encoded in the charge neutrality condition that sets
the electrodes' field inside the capacitor); the third term in
parenthesis is the charge density at point z' in the rest of the
system. It helps to think of two trains of ions rigidly moving
toward the oppositely charged electrodes, where the head of
each train continuously brings new ions to each corresponding
double layer. Outside the double layer, on the z < 0 side
of the capacitor, the nonzero charge density is initially due
to the negative ions leaving altogether toward positive z, while
the concentration of positive ions stays constant. The same
<figure><img src="image_12.png" /><figcaption>FIG. 8. Ratios of $\sigma(t)$, $\rho^{+L} = \rho(+L, t)$ and $\rho^{-L} = \rho(-L, t)$ versus their equilibrium values, as a function of time. The electrode charge density is piecewise parabolic. The blue dotted lines represent the ideal linear evolution of $\rho(\pm L, t)$ according to times (55), while the blue dashed line represents Eq. (58). On the left: $v = 200$, $\epsilon = 0.1$ and $q_+ : q_- = 1 : 2$, so that $t^*_- = 0.005$ and $t^*_+ = 0.01$; the curve for the electrode charge density starts at $v/[v + 2\epsilon^{-2}(q_+ + q_-)^{-1}] = 0.75$, as predicted by Eq. (58). On the right: $v = 50$, $\epsilon = 0.1$, and $q_+ : q_- = 1 : 10$, so that $t^*_- = 0.004$ and $t^*_+ = 0.04$; the curve for the electrode charge density starts at 0.73.</figcaption></figure>

happens on the $z > 0$ side, with reversed roles. In other words, as the tails of the two trains of ions travel at speeds $v_+$ and $v_-$ toward the center of the electrode, they define three regions of positive charge (no negative ions present), zero charge (both species present) and negative charge (no positive ions present). The two tails meet at $t = t_{\text{meet}}$: At this time the neutral region has shrunk to a point. For $t_{\text{meet}} < t < t^*_-$, the neutral region re-expands, this time because no ions are present in the central region anymore. The two trains continue

to move until time $t = t^*$, when the fastest species (negative ions) has reached the oppositely charged electrode. After that, the right part of the bulk solution is neutral, whereas the left part is still populated by positive ions that have not reached the electrode yet. They do so at time $t = t^*_+$. After this time, all ions are adsorbed and the whole bulk is neutral.

Solving the integral in Eq. (56) and nondimensionalizing, one obtains the following surface charge as a function of time:

$$|\sigma(t)| = \begin{cases} v + \frac{2}{(q_+ + q_-)\epsilon^2} \left[t\left(\frac{1}{t^*_+} + \frac{1}{t^*_-}\right) - \frac{1}{2}t^2\left(\frac{1}{t^*_+} + \frac{1}{t^*_-}\right)\right] & \text{if } 0 < t < t^*_- \\ v + \frac{2}{(q_+ + q_-)\epsilon^2} \left[\frac{1}{2} + \frac{t}{t^*_+} - \frac{t^2}{2t^*_+^2}\right] & \text{if } t^*_- < t < t^*_+ \end{cases}. \quad (58)$$

At equilibrium, the electrodes' charge is $|\sigma(t^*_+)| = v + 2/[(q_+ + q_-)\epsilon^2]$, as it should be, the second term representing the final amount of charge adsorbed on either electrode ($2n_0^2L$ in physical units). The simple $t^2$ dependence observed in the symmetric case is here split into two parts, each relevant before or after the faster species has reached the electrode. The two parabolas in Eq. (58) have different curvatures and are centered at different times. Equation (58) agrees with numerical data, as shown in Fig. 8.

As for the symmetric case, for times $t > t^*_-$ we expect a fast relaxation over the Gouy-Chapman lengths of the two counterionic double layers. However, the relaxation of $\sigma$ is necessarily dominated by the relaxation of the slower double layer. Since the double layers are coupled through the electrodes' surface charge (equal and opposite on the two electrodes), the relaxation of the fast ions is limited by slow changes in $\sigma$ and therefore happens with the same timescale as the slow double layer. Hence the exponential relaxation for $t > t^*_+$ occurs in a time

$$\tau = \frac{4L^2}{Dv^2q_+^2}, \quad (59)$$

which is observed at sufficiently high voltages, for all values of $\epsilon$ (see Fig. 5).

Again, to obtain Eq. (58), it was crucial to assume that ions do not affect the linear potential profile through almost all the capacitor. This happens if the change in electric potential (or field) caused by their motion is small with respect to the applied voltage. The analog of condition Eq. (43), determining the phase-space boundary of the regime we just described, is then

$$v > \frac{2}{(q_+ + q_-)\epsilon^2}. \quad (60)$$

We also assumed that the distance an ion can travel by diffusion in a time $t^*_+$ is much smaller than $2L$, so that we could neglect diffusion currents and only consider drift: This amounts to $v > 1/q_-$. In summary, the regime we just discussed, as for the symmetric case, is defined by Eq. (60) at small $\epsilon$ and coincides with the nonlinear regime for large $\epsilon$.

In the latter large-$\epsilon$ case, it is worth mentioning that the asymmetric-valence phenomenology is richer than for symmetric valences: In the $\epsilon \gtrsim 1$ regime, ions of different valence are depleted in different proportions (Fig. 9). This can give rise to a hybrid behavior in terms of relaxation rates that, however, covers a substantial range of voltages $v$ only for large, nonrealistic valences. For small $\epsilon$, between the purely nonlinear regime at $v > 1/q_-$ (nondepleted) and the unscreened
<figure><img src="image_13.png" /><figcaption>FIG. 9. Dukhin number $Du_n^± = n_±(0, t = ∞)/n_±(0, t = 0)$, quantifying depletion, for $q_+ : q_− = 1 : 2$ (top row) and 1:10 (bottom row). At large $ϵ = λ_D/L$ (left column), depletion sets in at lower voltages for the large-valence negative ions (blue) than for the small-valence positive ones (red). In the 1:10 case, there exists a range of voltages at which one species is fully depleted while the other is almost nondepleted. This does not happen in the symmetric valences case [Fig. 3(c)].</figcaption></figure>

nonlinear regime described by Eq. (60) (fully depleted), we have what we called in Sec. IIID a partially screened fully depleted regime. At such intermediate $v$, the initial relaxation is neither exponential nor polynomial in time. However, the relaxation can be divided in two parts, similarly to what happens in the unscreened regime: a first part before migration and depletion are completed and a second part concerning relaxation inside the two counterion-only double layers. The latter process cannot but happen on the same timescale as in Eq. (59). This intermediate regime is in the asymmetric valence case even more complex, because at a given applied voltage the two species present two different levels of depletion.

## V. FULLY ASYMMETRIC CASE: $q_+ \neq q_-$, $D_+ \neq D_-$

### A. Linear regime

We finally consider the fully asymmetric case, with ions of different valences $q_+ \neq q_-$ and diffusivities $D_+ \neq D_-$. In the linear regime, we solve the Poisson-Nernst-Planck equations analytically in the Laplace domain, as done for the fully symmetric case [22,24,25]. Since Eqs. (3) and (4) are coupled PDEs, two successive nontrivial diagonalizations, one in the time domain and one in the Laplace domain are needed, leading to solvable uncoupled equations for linear combinations of the ion densities. In summary, the first transformation aims to diagonalize drift currents. After passing to the Laplace domain and diagonalizing a second time, the diffusion equation in space can be solved. Once time and space derivatives do not appear anymore, a cumbersome calculation leads back to

the original basis {$\hat{n}_+(z, s), \hat{n}_-(z, s)$} (Laplace transforms of the ion densities) and boundary conditions can be imposed (potential at the electrodes, no flux through the electrodes, electroneutrality). As in the symmetric case, the modes characterizing the linear response can also be considered in the frequency domain via the impedance of the cell (see, e.g., Refs. [40–43] for partially asymmetric electrolytes with opposite valences but unequal diffusivities).

The nonzero poles of these functions represent the relaxation rates of the system and can be analyzed numerically. In the following, we take as reference $\hat{\rho}(-L, s)$, the Laplace transform of the charge density at contact $\rho(-L, t)$. This quantity depends on $ϵ$, on the valences, and on the additional dimensionless parameter,

$$ \delta = \frac{D_+}{D_-}. \qquad (61) $$

We assume here that the positive species is the slower one, so that $δ \le 1$. A numerical analysis of $\hat{\rho}$ reveals, for any $δ < 1$, the presence of new poles that did not exist in the case $δ = 1$.

For large $ϵ$, $\hat{\rho}$ has twice as many poles (and zeros). The number of poles (and zeros) stays infinite, but each pole from the $δ = 1$ case splits into two poles and a zero as soon as $δ < 1$. The poles can be classified into two independent hierarchies of diffusive timescales: one, {$τ_{+,i}$}, for the positive species and one, {$τ_{-,i}$}, for the negative species. While their exact values can be retrieved numerically, these times are well described by the following equation:

$$ \tau_{\pm,i} \approx \frac{4L^2}{(2i + 1)^2 \pi^2 D_{\pm}}, \quad \text{with } i = 0, 1, 2, \dots, \qquad (62) $$

asymptotically exact for $ϵ \to \infty$. Still assuming that positive ions are slow compared to negative ones, the slowest timescale, i.e., the dominant one at large times, is $τ_{+,0} = 4L^2/(\pi^2 D_+)$. Note that the valences of the two species do not enter in the expression for the relaxation times $τ_{±,i}$; they do play a role, though, in determining the importance of each mode (e.g., the weight of positive modes increases with $q_+$).

For small $ϵ$, the picture changes and the slowest timescale approaches

$$ \tau_{\text{NH}} = \frac{4L^2}{\pi^2 D_{\text{NH}}}, \quad \text{with } D_{\text{NH}} = \frac{(q_+ + q_-)D_+ D_-}{q_+ D_+ + q_- D_-}, \qquad (63) $$

featuring the so-called Nernst-Hartley diffusion coefficient $D_{\text{NH}}$ [44]. Its presence reflects the fact that ions diffuse together under the action of an internal electric field that pulls the slowest species, while slowing down the fastest one [44]. The scaling of the relaxation time is shown in Fig. 10(a) for the equal-valence case. Figures 10(b) and 10(c) show instead how the slowest relaxation time $τ_0$ transitions from $τ_{+,0}$ at large $ϵ$ [Fig. 10(b)] to $τ_{\text{NH}}$ at small $ϵ$ [Fig. 10(c)], for the 1:2 and 2:1 cases, at any $δ$.

The pole corresponding to Eq. (63), the closest one to the origin, is, however, not the one with the largest residue. This is important because the residue of a pole is proportional to the weight of its corresponding mode in the time domain: a weight larger than the slowest mode’s weight allows a fast mode to be visible, at least at short times. A fast, but large-weight mode represents a relaxation phenomenon that is dominant for a certain time and eventually fades out giving way to slower
<figure><img src="image_14.png" /><figcaption>FIG. 10. (a) Slowest relaxation time $\tau_0$ extracted numerically from the analytical $\hat{\rho}(-L, s)$ in units of its asymptotic value from Eq. (63), as a function of $\delta$, for small values of $\epsilon$. The dotted horizontal line represents Eq. (63). $q_+$: $q_-$: 1: 1. [(b) and (c)] Same quantity, in units of its large-$\epsilon$ value $\tau_{+,0}$ from Eq. (62) (b), and in units of its small-$\epsilon$ value from Eq. (63) (c). Continuous lines are for a 1:2 electrolyte, dashed lines are for a 2:1 electrolyte.</figcaption></figure>

processes. Returning to the RC-circuit analogy, one can estimate what the time with the largest weight $\tau_w$ might look like, assuming that it represents double layer charging. Recomputing the electric resistance times unit surface of Eqs. (23) and (52) for a bulk with asymmetric diffusivities gives

$$R = \frac{2L - 2\lambda_D}{\beta e^2 n_0 q_+ q_- (q_+ D_+ + q_- D_-)} \quad (64)$$

Using a capacitance per unit surface $\epsilon_0\epsilon_r/\lambda_D$ for each electrode, one finds the characteristic time

$$\tau_{RC} = \frac{L\lambda_D - \frac{\lambda_D^2}{2}}{D_{ave}}, \quad \text{with} \quad D_{ave} = \frac{q_+D_+ + q_-D_-}{q_+ + q_-}, \quad (65)$$

where the factor 1/2 in front of the $\lambda_D^2$ term was added ad hoc so as to match the fully symmetric case [Eq. (19)].

This expression contains an average of the diffusion coefficients weighted by ion valences, $D_{ave}$, rather than the Nernst-Hartley diffusivity of Eq. (63). It explains results both from the numerical analysis of $\hat{\rho}$ [Fig. 11(a)] and from solutions of the Poisson-Nernst-Planck equations (Fig. 12, small $v$ and small $\lambda_D/L$). An analysis of the weights of the mode $\tau_w \simeq \tau_{RC}$ and of the slowest mode $\tau_0 \simeq \tau_{NH}$ is presented in Fig. 11(b); $\tau_w$ has a larger weight on the whole range of $\delta$.

For $\delta \to 1$, the weight of $\tau_0$ goes to 0—the time $\tau_{NH}$ is indeed absent from the linear analysis of the symmetric valence case.

In summary, at large $\epsilon$ the relaxation is fully described by Eq. (62), with positive and negative ions decoupled as in the fully symmetric case. At small $\epsilon$, the double layer charging occurs at early times with a timescale $\tau_w \simeq \tau_{RC}$, described by Eq. (65). Subsequently, the purely diffusive mode from Eq. (63) sets in, featuring the Nernst-Hartley coefficient and signaling that positive and negative ions diffuse together. An inspection of the curves for $n_+(z, t)$ highlights the origin of this two-step relaxation. Initially, the double layer builds up mostly thanks to the faster species: The equilibrium total charge density is reached through an immediate rearrangement of the fast species around the instantaneous local concentration of the slow species. Mass concentrations are therefore different from the equilibrium ones: There is a neutral excess of mass at the electrode with the same sign as the slow species (as the slow species has not had time to escape), and a defect of mass at the opposite side. This requires a mass relaxation process, which happens precisely with the diffusive timescale of Eq. (63). Note that even though the diffusive scaling might recall depletion, cf. Eq. (31), the phenomenon we just described is linear in all respects and occurs also for $v \to 0$.

<figure><img src="image_15.png" /><figcaption>FIG. 11. (a) Largest-weight relaxation time $\tau_w$ extracted numerically from $\hat{\rho}(-L, s)$, in units of its theoretical value from Eq. (65), as a function of $\delta = D_+/D_-$, for different values of $\epsilon = \lambda_D/L$. 1:1 electrolyte. (b) Ratio of the residues of $\hat{\rho}(-L, s)$ with respect to the pole corresponding to $\tau_0$ and to the pole corresponding to $\tau_w$, as a function of $\delta$; 1:1 electrolyte.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 12. Relaxation times $\tau$ as a function of $v$, for diffusivity ratios $\delta = D_+/D_- = 1/1.1$ or $1/10$ and valence ratios $q_+:q_- = 1:2$ or $2:1$. Times are in units of $L^2/D_+$ and different colors correspond to different values of $\epsilon = \lambda_D/L$. As in Fig. 5, times are extracted from exponential relaxation(s) of $\sigma$ and, when present, empty symbols distinguish late-time processes from earlier-time processes (full symbols). Two black arrows point at the values of $\tau_{0,+}$ and $\tau_{NH}$, respectively, from Eqs. (62) and (63). The colored dotted lines represent the purely nonlinear time obtained from Eqs. (64) and (51) and preceding ones (see Fig. 6): by construction, their small-$v$ value is $\tau_{RC}$ from Eq. (65). The dashed black line represents the time from Eq. (42) with $D = D_+$. The dashed gray lines are the two times from Eq. (66) observed at short times, where the relaxation is linear in time.</figcaption></figure>

## B. Nonlinear regimes

An analysis of the relaxation times observed in the nonlinear regime is presented in Fig. 12 for $\delta = 1/1.1$ and $1/10$ and for valences $1:2$ and $2:1$. Results obtained for the symmetric diffusivity cases seem to extend naturally to asymmetric diffusivities. In the purely nonlinear regimes, the capacitance computed in Eqs. (48)–(51) is still applicable and can be combined with the resistance from Eq. (64). Qualitatively, at least for small valences, this results in a simple vertical shift of the exponential curve representing increased nonlinear capacitance in Fig. 12, such that at small $\epsilon$ the curve coincides with the linear regime time from Eq. (65). The curve is most evident for the $\delta = 1/1.1$, $1:2$ case at $\epsilon = 0.01$ and is partially masked by high-weight diffusive timescales for other parameters. Its shape mostly depends on the valence ratio through the capacitance, as per Fig. 6.

In the unscreened ($v \to \infty$) fully depleted regime, the discussion for the symmetric and asymmetric valence cases still applies, with the *caveat* that asymmetry in the diffusivities will affect the initial linear drag phase. The times from

Eq. (55) must then be replaced by

$$t_{\pm}^* = \frac{2L}{v_{\pm}} = \frac{2L^2}{D_{\pm}q_{\pm}v}. \qquad (66)$$

Once each ion train has reached the electrode, each counterionic double layer will start its internal exponential relaxation with relaxation time $4L^2/(D_{\pm}v^2q_{\pm}^2)$. For the same reasons that led us to Eq. (59) in the partially asymmetric case, the last relaxation should feature the slower of the two relaxation times of the two counterionic double layers:

$$\tau = \frac{4L^2}{v^2 \min(D_+q_+^2, D_-q_-^2)}. \qquad (67)$$

This seems compatible with the numerical results presented in Fig. 12, that are, however, computationally more involved to obtain for significantly distinct values of valences and diffusivities.
Results from this section are summarized in Fig. 4 in Ref. [26].

## VI. NONIDEAL BEHAVIOR

In this section, we analyze in brief the applicability range of the Poisson-Nernst-Planck theory. Specifically, we discuss the impact of ionic correlations due to electrostatics and to ion packing by estimating the electrostatic coupling parameter and the maximum ion density in the device, across our regime diagram. Our model indeed provides a simple framework to benchmark further simulations, theories, or experiments, allowing for a more profound assessment of the effect of nonidealities (e.g., electrostatic correlations, ion size, polarizability, and electrode geometry).

### A. Deviation from mean field

Ionic correlations are an intrinsically discrete phenomenon, whose magnitude depends on local density and whose effect can be anisotropic and counterintuitive. A possible approach consists in introducing a constant effective correlation length scale within mean-field models. The effects of this effective correlation length on the small-voltage EDLC dynamics have been studied in Refs. [45,46], which reported relatively small deviations from the mean-field dynamics for reasonable values of correlation length. Nonetheless, the mapping between the effective correlation length and physical parameters is somewhat elusive. We take here a different approach, relying only on the electrostatic coupling parameter.

When only one kind of ions is present in the double layer, which is often the case in the nonlinear regimes, a useful tool to quantify the importance of ion-ion electrostatic correlations is the electrostatic coupling parameter [47–50]. It represents the squared ratio between the electrostatic energy between two ions of charge eq at a distance a and the thermal energy k_B T:

$$ \Xi = \left( \frac{q^2 l_B}{a} \right)^2. \qquad (68) $$

If ions adhere to a wall of charge density eΣ and the system is electroneutral, then the typical distance is fixed by the wall charge density (a = $\sqrt{q/\Sigma}$) and the coupling parameter assumes the familiar expression $\Xi = 2\pi l_B^2 q^3 \Sigma$ [47–50]. Poisson-Boltzmann is an excellent approximation when $\Xi \lesssim 1$, while for $\Xi \gtrsim 10$ the ion density noticeably deviates from the mean-field one [51–53]. In addition, since usual ion sizes are smaller than $l_B = 0.7$ nm in water at ambient temperature, one can reasonably assume that absence of electrostatic correlations implies absence of ion packing effects. $\Xi \lesssim 1$ is therefore a good proxy for the absence of ionic correlations.

### B. Fully depleted nonlinear regime

We analyze here the fully depleted nonlinear regime, described in Secs. III D, IV C, and V B, and represented by the light-shaded regions of Figs. 3(a) and 4(a) in Ref. [26]. In this regime, the vicinity of the electrodes is populated by counterions only. In addition, the system (electrode and counterions) is in general not electroneutral as the electrode charge Σ can be larger than the amount of ions present in the system $2n_0L$. Looking, for simplicity, at the symmetric 1:1 case, at equilibrium all ions are in the vicinity of the electrode.

In the worst-case scenario, they are adsorbed on the plate and the typical distance between them is $a \simeq 1/\sqrt{2n_0L}$, so that the coupling parameter (68) takes the form

$$ \Xi = \frac{l_B}{4\pi L \epsilon^2} = 2l_B^2 L n_0, \qquad (69) $$

which is independent of v.

For $v > 1$ and $\epsilon > 1$, we always have $\Xi < 1$, as realistically $l_B \ll L$ for any aqueous-electrolyte EDLC (where $l_B = 0.7$ nm), so that correlations are negligible. Referring to the regime diagram in Fig. 3(a) in Ref. [26], in the no-screening small-$\epsilon$ sector [$v > 1$ and $v^{-1/2} < \epsilon < 1$ as per Eq. (43)], we have $\Xi < v l_B / (4\pi L)$. For $v \lesssim 50$ for aqueous electrolytes, it is still $\Xi \lesssim 1$ for realistic EDLCs and correlations have little relevance. Finally, in the partially screened depleted sector [$v > 1$ and $1/\sinh(v/2) < \epsilon < v^{-1/2}$], $\Xi$ may exceed unity for some combinations of voltage, concentration, permittivity, pore size and temperature, and the validity of mean field should be checked on a case-by-case basis computing $\Xi$ through Eq. (69). For example, a 5-mM 1:1 aqueous electrolyte in a 3-mm-wide EDLC (experimental parameters from Ref. [54]) would result in $\Xi$ exceeding 4000. In contrast, a 1-µm cell in the same conditions would experience a fairly weak coupling, with $\Xi$ close to 1.

### C. Purely nonlinear regime

We now turn to the purely nonlinear regime, where depletion is absent. This is described in Secs. III C, IV B, and V B, represented in the central-bottom regions of Figs. 3 and 4 in Ref. [26], and identified, in the fully symmetric case, by the conditions $v > 1$ and $\epsilon < (20 \sinh(v/2))^{-1}$ as per Eq. (37). In this regime, the system (electrode and double-layer) is electroneutral. While a theory of electrostatic coupling for two-species double layers does not exist, for large voltages counterions are expelled from the EDL, so that one can still use $\Xi = 2\pi l_B^2 q^3 \Sigma$. Once equilibrium is reached, Eq. (33) for the electrode charge gives

$$ \Xi = \frac{l_B}{L\epsilon} \sinh\left(\frac{v}{2}\right) \qquad (70) $$

for monovalent ions. By Eq. (37) again, in this sector $\Xi \le \frac{l_B}{20L\epsilon^2}$, with the equality holding only at the boundary with the depleted region. Replacing for $\epsilon$, this means

$$ \Xi \le \frac{4\pi}{10} l_B^2 L n_0 \approx \frac{LX}{3 \text{ nm}}, \qquad (71) $$

where X is the numeric value of the molarity in moles per liter and we assumed an aqueous solution at ambient temperature. In practice, for nanopores with $L \approx 100$ nm and concentrations below 30 mM (i.e., $X < 0.03$), $\Xi$ stays below 1 in the whole sector. While Eq. (71) sets an upper bound, $\Xi$ decreases exponentially with decreasing v, rapidly moving away from this bound as per Eq. (70). As a result, for $v \sim 1$, our results stay relevant even at substantially higher concentrations. For the millimeter-scale EDLCs from Ref. [54], at concentrations between 5 and 400 mM and applied voltages between 0.05 and 0.6 V, Eq. (70) yields a $\Xi$ ranging from 0.08 to 230, corresponding to a very weak to moderately strong coupling.
### D. Early-time behavior ($t \to 0$)

Even when mean field is in principle bound to fail, it should still capture the early response, at least for initial concentrations low enough that correlations can be neglected at small times. At $t = 0$ it is $a = n_0^{-1/3}$, which translates to

$$ \Xi_{\text{bulk}} = (l_B n_0^{1/3})^2 \lesssim 1. \qquad (72) $$

This makes our analysis valid at least for short times for any electrolyte concentrations of order 1 M or below. In cells from Ref. [54], the early-time electrostatic coupling is completely negligible, with $\Xi$ between 0.01 and 0.2.

### E. Final comments

Overall, while in many regions of parameter space the mean-field assumption strictly holds, it is in general not possible to determine regions in the $(v, \epsilon)$ diagram where mean field surely fails, so that the value of the electrostatic coupling parameter $\Xi$ should be checked case by case. This is because electrostatic correlations can introduce new length scales that make it impossible to reduce the $\Xi \lesssim 1$ condition to a simple condition on $v$ and $\epsilon$, with $\Xi$ becoming in general a different function of the system size, ion concentration, voltage, solvent permittivity, and temperature. We also note that Eqs. (3) for the ion dynamics neglect the coupling with the solvent, not only in terms of electrokinetic effects (see, e.g., Refs. [16,55]) but also Maxwell-Stefan fluxes which may contribute as the salt concentration increases (see, e.g., Ref. [56]).

Finally, in the present work we have assumed that there are no electrochemical reactions (in particular involving the solvent) during the whole charging process. Considering the most used media, the electrochemical window of water is $v \simeq 50$, that of organic solvents $v \simeq 60$, while that of ionic liquids at most $v \simeq 120$ [57]. Our numerical results beyond these limits have the purpose of highlighting the dominant process, which is in general more evident far away from the transitions with other regimes.

## VII. CONCLUSION

Despite its applications in electrochemistry, the relaxation to equilibrium of EDLCs has been somewhat elusive, in particular for the most common case of electrolytes with different valences and/or diffusivities. This work, together with Ref. [26], tries to fill this gap by giving a new perspective on existing results for the linear fully symmetric case and by subsequently addressing the asymmetric cases and the nonlinear regimes. We characterize the relaxation behavior focusing on the dominant timescales, which are usually the slowest ones. We define different regimes in the parameter space spanned by applied voltage and ion concentration, providing analytical boundaries between regimes whenever possible. For low ion concentrations, the behaviors of positive and negative ions can be mutually decoupled, with a relaxation that is diffusive (exponential in time) in the small-voltage regime and dominated (linear in time) in the large-voltage one. Beyond this low concentration limit, the picture is different. At asymptotically small voltage, the relaxation is faster than diffusion, yet followed by a slower diffusion process of neutral mass in

the case of asymmetric diffusivities. As voltage increases, the relaxation time increases exponentially due to nonlinearity, until the potential difference is strong enough to deplete ions from the bulk. Then relaxation gradually goes from exponential to linear in time, through a regime where the electric field across the bulk is highly varying. Last, if the electrode charge is much larger than the charge that the confined electrolyte can possibly neutralize, then the picture from low salt concentration is recovered, where ions are dragged at constant velocity to the electrodes.

The results presented here rely on the validity of the mean-field approximation. We discussed in Sec. VI how exceedingly large ion concentrations, caused either by too-large initial concentrations or applied voltages, might bring about effects beyond mean field that can at least partially invalidate our results, in particular in the partially screened fully depleted regime. The significance of electrostatic correlations [48,58–60], but also of ion pair formation [61] and finite-size effects [18,62–67], may depend on the specific kind of ions used, the permittivity of the medium [50,68], the pore size [65,69], the temperature, and the nature of the electrode.

Our analysis contributes to the understanding of electrochemical devices and of confined charged materials, providing a clear mapping between any point in parameter space and mean-field relaxation times. This facilitates comparison among theories that include ion specificity, finite ion sizes, ionic correlations, hydrodynamics, or geometric effects [17–19,66,70,71]. Finally, our work paves the way to the design of optimization procedures for the charging of EDLCs [21–23], with promising applications to energy production and recovery.

## ACKNOWLEDGMENTS

This work has received funding from the European Union’s Horizon 2020 and Horizon Europe research and innovation programs under the Marie Skłodowska-Curie Grants No. 674979-NANOTRANS (I.P., P.B.W., B.R., and E.T.), No. 101034413 (I.P.), and No. 101119598-FLUXIONIC (M.D., B.R., and E.T.), as well as from the European Research Council under Grant No. 863473 (B.R.). B.R. acknowledges financial support from the French Agence Nationale de la Recherche (ANR) under Grant No. ANR-21-CE29-0021-02 (DIADEM). I.P. thanks Anděla Šarić for further support at ISTA.

## APPENDIX: NUMERICAL METHODS

The Poisson-Nernst-Planck equations are solved numerically using a flux-conservative finite-difference method. We discretize space into nodes, positioned at $z_{k+\frac{1}{2}}$ for $k = 0, \dots, N-1$, and edges, located at $z_k$ for $k = 0, \dots, N$. If the spacing between nodes and between edges is constant, then we have $z_{k+\frac{1}{2}} = -L + (k + \frac{1}{2})\Delta z$ and $z_k = -L + k\Delta z$, with $\Delta z = 2L/N$. The extension to a nonlinear spacing is straightforward and useful. For simplicity, we will describe the algorithm assuming constant spacing.

Ion densities $n_\pm$ and potential $\phi$ are defined on nodes $z_{k+\frac{1}{2}}$, while electric field and ionic currents are defined on edges $z_k$; this reduces the error associated with numerical derivation
or integration. Densities are initially set to their $t = 0$ value, $n_±(z, 0) = n_0^±$, for all nodes. Then the density and potential profiles are evolved in time in steps of $\Delta t$. More precisely, at every step $i$:

(1) The ionic contribution to the electric field is computed via Gauss' theorem, by numerically integrating the following charge density:

$$ \rho(z_{k+\frac{1}{2}}, i \Delta t) = \bar{q}_- n_-(z_{k+\frac{1}{2}}, i \Delta t) + \bar{q}_+ n_+(z_{k+\frac{1}{2}}, i \Delta t). $$

(2) The ionic contribution to the potential is computed by numerically integrating the ionic contribution of the electric field.

(3) The electrodes' (constant) contribution to the electric field is computed by imposing that the overall potential difference across the capacitor (ions' and electrodes' contributions) be $2V(i \Delta t) = 2V_0$ for $i > 0$. This corresponds to computing the surface charge density $\sigma$ of the electrodes at time $i \Delta t$, which is proportional to the electrodes' electric field through the dielectric permittivity.

(4) The overall electric field (ions' and electrodes' contributions) determines the ionic currents through the discrete analogous of Eq. (1).

(5) Ion densities at time $(i+1)\Delta t$ are computed through the discrete analogous of Eq. (2), starting from the just determined currents at time $i \Delta t$.

Iteration stops when the prescribed final time is reached.

The algorithm is flux conservative [72] because it conserves the numerical integral of the densities for both ionic species, i.e., the total number of ions, up to machine precision. This is a consequence of the fact that density updates are computed by the finite-difference equivalent of Eq. (2). At any time, the total number of ions is indeed

$$ \sum_k n_±(z_{k+\frac{1}{2}}, t + \Delta t) \Delta z = \sum_k \left( n_±(z_{k+\frac{1}{2}}, t) + \Delta t \frac{j_±(z_k, t) - j_±(z_{k+1}, t)}{\Delta z} \right) \Delta z $$

[1] P. Simon and Y. Gogotsi, Materials for electrochemical capacitors, Nat. Mater. 7, 845 (2008).

[2] S. Chen, L. Qiu, and H.-M. Cheng, Carbon-based fibers for advanced electrochemical energy storage devices, Chem. Rev. 120, 2811 (2020).

[3] P. Simon and Y. Gogotsi, Perspectives for electrochemical capacitors and related devices, Nat. Mater. 19, 1151 (2020).

[4] M. Winter and R. J. Brodd, What are batteries, fuel cells, and supercapacitors? Chem. Rev. 104, 4245 (2004).

[5] D. P. Dubal, N. R. Chodankar, D.-H. Kim, and P. Gomez-Romero, Towards flexible solid-state supercapacitors for smart and wearable electronics, Chem. Soc. Rev. 47, 2065 (2018).

$$ \begin{aligned} &= \left( \sum_k n_±(z_{k+\frac{1}{2}}, t) \Delta z \right) + \Delta t (j_±(z_0, t) - j_±(z_N, t)) \\ &= \sum_k n_±(z_{k+\frac{1}{2}}, t) \Delta z. \end{aligned} \quad (A1) $$

The last equality follows from the zero-current condition on the electrodes, placed at $z_0 = -L$ and $z_N = +L$. This property is exact for constant spacing $\Delta z$, as shown, but also for irregular spacings, provided that $\Delta z$ in Eq. (A1) is replaced everywhere by $(\Delta z)_{k+\frac{1}{2}} = z_{k+1} - z_k$. The usage of irregular spacings, with nodes and edges more dense close to the walls and less dense in the bulk, is essential to speed up the calculation when $L$ is orders of magnitude larger than the thickness of the double layer and it also favors stability in the nonlinear regime, as mentioned below.

A necessary condition for the stability of the algorithm is that the time step $\Delta t$ verify the Courant-Friedrichs-Lewy condition [72], requiring that $\Delta t < \Delta z^2/\max_\alpha\{D_\alpha\}$. Otherwise said, assuming purely diffusive motion, the typical distance traveled by the faster species in a time $\Delta t$ must be smaller than the lattice spacing. For the numerical calculations reported in the paper, the Courant factor $\Delta t/(\Delta z^2/\max_\alpha\{D_\alpha\})$ is between 0.05 and 0.9.

Last, it is necessary to have a sufficient number of nodes inside the double layer, for the results to be accurate and stable. This is particularly important when nonlinear effects set in, as densities and potential curves become much steeper in the double layer. For this reason, irregular node spacing is fundamental and was chosen in such a way that the distribution of nodes was linear in the bulk region and became exponentially dense closer to the electrodes. In addition, nonlinear effects favor electromigration with respect to diffusion and can elicit a rather violent response to the perturbation. This makes the Courant condition insufficient at strong voltages, calling for a further reduction of the time step.

In contrast to canonical simulations, the results for grand-canonical systems (shown in Fig. 7) are obtained by solving numerically the nonlinear Poisson-Boltzmann equation describing the final equilibrium state. In that case the Debye length sets the salt concentration in the reservoir instead of that in the system. In practice, we solve the Poisson-Boltzmann equation by an iterative procedure.

[6] D. Brogioli, Extracting renewable energy from a salinity difference using a capacitor, Phys. Rev. Lett. 103, 058501 (2009).

[7] N. Boon and R. van Roij, 'Blue energy' from ion adsorption and electrode charging in sea and river water, Mol. Phys. 109, 1229 (2011).

[8] M. Janssen, A. Härtel, and R. van Roij, Boosting capacitive blue-energy and desalination devices with waste heat, Phys. Rev. Lett. 113, 268501 (2014).

[9] V. Vivier and M. E. Orazem, Impedance analysis of electrochemical systems, Chem. Rev. 122, 11131 (2022).

[10] G. Barbero and L. R. Evangelista, Impedance spectroscopy of a cell: The role of the ions, in *Adsorption Phenomena and*
Anchoring Energy in Nematic Liquid Crystals (CRC Press, Boca Raton, FL, 2005).

[11] G. Jeanmairet, B. Rotenberg, and M. Salanne, Microscopic simulations of electrochemical double-layer capacitors, *Chem. Rev.* **122**, 10860 (2022).

[12] A. A. Kornyshev, Double-layer in ionic liquids: Paradigm change? *J. Phys. Chem. B* **111**, 5545 (2007).

[13] M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, Double layer in ionic liquids: Overscreening versus crowding, *Phys. Rev. Lett.* **106**, 046102 (2011).

[14] J. Wu, Understanding the electric double-layer structure, capacitance, and charging dynamics, *Chem. Rev.* **122**, 10821 (2022).

[15] M. Janssen, Curvature affects electrolyte relaxation: Studies of spherical and cylindrical electrodes, *Phys. Rev. E* **100**, 042602 (2019).

[16] A. J. Asta, I. Palaia, E. Trizac, M. Levesque, and B. Rotenberg, Lattice Boltzmann electrokinetics simulation of nanocapacitors, *J. Chem. Phys.* **151**, 114104 (2019).

[17] C. Lian, M. Janssen, H. Liu, and R. van Roij, Blessing and curse: How a supercapacitor's large capacitance causes its slow charging, *Phys. Rev. Lett.* **124**, 076001 (2020).

[18] K. Ma, M. Janssen, C. Lian, and R. van Roij, Dynamic density functional theory for the charging of electric double layer capacitors, *J. Chem. Phys.* **156**, 084101 (2022).

[19] J. Yang, M. Janssen, C. Lian, and R. van Roij, Simulating the charging of cylindrical electrolyte-filled pores with the modified Poisson–Nernst–Planck equations, *J. Chem. Phys.* **156**, 214105 (2022).

[20] L. Varela, S. Andraus, E. Trizac, and G. Téllez, Relaxation dynamics of two interacting electrical double-layers in a 1D Coulomb system, *J. Phys.: Condens. Matter* **33**, 394001 (2021).

[21] K. Breitsprecher, C. Holm, and S. Kondrat, Charge me slowly, I am in a hurry: Optimizing charge-discharge cycles in nanoporous supercapacitors, *ACS Nano* **12**, 9733 (2018).

[22] I. Palaia, Charged systems in, out of, and driven to equilibrium: From nanocapacitors to cement, Ph.D. thesis, Paris-Saclay University, 2019, https://cnrs.hal.science/tel-02926717/.

[23] K. Breitsprecher, M. Janssen, P. Srimuk, B. L. Mehdi, V. Presser, C. Holm, and S. Kondrat, How to speed up ion transport in nanopores, *Nat. Commun.* **11**, 6085 (2020).

[24] M. Z. Bazant, K. Thornton, and A. Ajdari, Diffuse-charge dynamics in electrochemical systems, *Phys. Rev. E* **70**, 021506 (2004).

[25] M. Janssen and M. Bier, Transient dynamics of electric double-layer capacitors: Exact expressions within the Debye-Falkenhagen approximation, *Phys. Rev. E* **97**, 052616 (2018).

[26] I. Palaia, A. J. Asta, M. Dutta, P. B. Warren, B. Rotenberg, and E. Trizac, companion paper, Charging dynamics of electric double-layer nanocapacitors in mean field, *Phys. Rev. Lett.* **135**, 148002 (2025).

[27] R. J. Hunter, *Foundations of Colloid Science*, 2nd ed. (Oxford University Press, Oxford, 2001).

[28] D. Andelman, Introduction to electrostatics in soft and biological matter, in *Soft Condensed Matter Physics in Molecular and Cell Biology*, edited by W. C. K. Poon and D. Andelman (CRC Press, Boca Raton, FL, 2006), pp. 97–122.

[29] G. Barbero and A. L. Alexe-Ionescu, Role of the diffuse layer of the ionic charge on the impedance spectroscopy of a cell of liquid, *Liq. Cryst.* **32**, 943 (2005).

[30] G. Barbero, F. Batalioto, and A. M. Figueiredo Neto, Theory of small-signal ac response of a dielectric liquid containing two groups of ions, *Appl. Phys. Lett.* **92**, 172908 (2008).

[31] J. R. Macdonald, Double layer capacitance and relaxation in electrolytes and solids, *Trans. Farad. Soc.* **66**, 943 (1970).

[32] A. A. Kornyshev and M. A. Vorotyntsev, Electric current across the metal-solid electrolyte interface II. Low-amplitude alternating current, *Phys. Stat. Sol. (a)* **39**, 573 (1977).

[33] A. Kornyshev and M. Vorotyntsev, Conductivity and space charge phenomena in solid electrolytes with one mobile charge carrier species, a review with original material, *Electrochim. Acta* **26**, 303 (1981).

[34] L. Bocquet and E. Charlaix, Nanofluidics, from bulk to interfaces, *Chem. Soc. Rev.* **39**, 1073 (2010).

[35] Electric double layers, in *Fundamentals of Interface and Colloid Science*, Vol. 2, edited by J. Lyklema (Academic Press, San Diego, CA, 1995).

[36] J. R. Macdonald, Theory of the differential capacitance of the double layer in unadsorbed electrolytes, *J. Chem. Phys.* **22**, 1857 (1954).

[37] F. Beunis, F. Strubbe, M. Marescaux, J. Beeckman, K. Neyts, and A. R. M. Verschueren, Dynamics of charge transport in planar devices, *Phys. Rev. E* **78**, 011502 (2008).

[38] F. Beunis, F. Strubbe, M. Marescaux, K. Neyts, and A. R. M. Verschueren, Diffuse double layer charging in nonpolar liquids, *Appl. Phys. Lett.* **91**, 182911 (2007).

[39] D. C. Grahame, Diffuse double layer theory for electrolytes of unsymmetrical valence types, *J. Chem. Phys.* **21**, 1054 (1953).

[40] I. Lelidis and G. Barbero, Effect of different anionic and cationic mobilities on the impedance spectroscopy measurements, *Phys. Lett. A* **343**, 440 (2005).

[41] G. Barbero and I. Lelidis, Evidence of the ambipolar diffusion in the impedance spectroscopy of an electrolytic cell, *Phys. Rev. E* **76**, 051501 (2007).

[42] C. Chassagne, E. Dubois, M. L. Jiménez, J. P. M. van der Ploeg, and J. van Turnhout, Compensating for electrode polarization in dielectric spectroscopy studies of colloidal suspensions: Theoretical assessment of existing methods, *Front. Chem.* **4**, 1 (2016).

[43] A. Antonova, G. Barbero, L. R. Evangelista, and P. Tilli, Ambipolar diffusion in the low frequency impedance response of electrolytic cells, *J. Stat. Mech.* (2020) 043202.

[44] R. A. Robinson and R. H. Stokes, *Electrolyte Solutions: Second Revised Edition* (Dover, Garden City, NY, 2002), pp. 286–292.

[45] H. Zhao, Diffuse-charge dynamics of ionic liquids in electrochemical systems, *Phys. Rev. E* **84**, 051504 (2011).

[46] D. Fertig and M. Janssen, Charging dynamics of electric double layer capacitors including beyond-mean-field electrostatic correlations, *Phys. Rev. E* **112**, 025414 (2025).

[47] R. R. Netz and H. Orland, Beyond Poisson-Boltzmann: Fluctuation effects and correlation functions, *Eur. Phys. J. E* **1**, 203 (2000).

[48] A. Naji, M. Kanduč, J. Forsman, and R. Podgornik, Perspective: Coulomb fluids–Weak coupling, strong coupling, in between and beyond, *J. Chem. Phys.* **139**, 150901 (2013).

[49] E. Trizac and L. Samaj, Like-charge colloidal attraction: A simple argument, in *Proceedings of the International School of Physics Enrico Fermi*, Vol. 184, edited by C. Bechinger, F.
Sciortino, and P. Ziherl (Società Italiana di Fisica, Bologna; IOS Press, Amsterdam, 2012), pp. 61–73.

[50] I. Palaia, A. Goyal, E. Del Gado, L. Šamaj, and E. Trizac, Like-charge attraction at the nanoscale: Ground-state correlations and water destructuring, *J. Phys. Chem. B* **126**, 3143 (2022).

[51] A. G. Moreira and R. R. Netz, Strong-coupling theory for counter-ion distributions, *Europhys. Lett.* **52**, 705 (2000).

[52] C. D. Santangelo, Computing counterion densities at intermediate coupling, *Phys. Rev. E* **73**, 041512 (2006).

[53] L. Šamaj, M. Trulsson, and E. Trizac, Strong-coupling theory of counterions between symmetrically charged walls: From crystal to fluid phases, *Soft Matter* **14**, 4040 (2018).

[54] C. Zhao, T. Yang, S. Jin, and B. Wu, Measurement of electric double layer charging dynamics on platinum electrodes in aqueous solutions of alkali sulfates and nitrates, *J. Phys. Chem. C* **128**, 5964 (2024).

[55] V. Lobaskin and R. R. Netz, Diffusive-convective transition in the non-equilibrium charging of an electric double layer, *Europhys. Lett.* **116**, 58001 (2016).

[56] B. Balu and A. S. Khair, Role of Stefan–Maxwell fluxes in the dynamics of concentrated electrolytes, *Soft Matter* **14**, 8267 (2018).

[57] M. Armand, F. Endres, D. R. MacFarlane, H. Ohno, and B. Scrosati, Ionic-liquid materials for the electrochemical challenges of the future, *Nat. Mater.* **8**, 621 (2009).

[58] B. D. Storey and M. Z. Bazant, Effects of electrostatic correlations on electrokinetic phenomena, *Phys. Rev. E* **86**, 056303 (2012).

[59] M. Kanduč, A. Naji, J. Forsman, and R. Podgornik, Dressed counterions: Strong electrostatic coupling in the presence of salt, *J. Chem. Phys.* **132**, 124701 (2010).

[60] H. Bonneau, V. Démery, and E. Raphaël, Temporal response of the conductivity of electrolytes, *J. Stat. Mech.* (2023) 073205.

[61] R. M. Adar, T. Markovich, and D. Andelman, Bjerrum pairs in ionic solutions: A Poisson-Boltzmann approach, *J. Chem. Phys.* **146**, 194904 (2017).

[62] I. Borukhov, D. Andelman, and H. Orland, Steric effects in electrolytes: A modified Poisson-Boltzmann equation, *Phys. Rev. Lett.* **79**, 435 (1997).

[63] D. Ben-Yaakov, D. Andelman, D. Harries, and R. Podgornik, Beyond standard Poisson-Boltzmann theory: Ion-specific interactions in aqueous solutions, *J. Phys.: Condens. Matter* **21**, 424106 (2009).

[64] J. Jiang, D. Cao, D. E. Jiang, and J. Wu, Time-dependent density functional theory for ion diffusion in electrochemical systems, *J. Phys.: Condens. Matter* **26**, 284102 (2014).

[65] A. Härtel, Structure of electric double layers in capacitive systems and to what extent (classical) density functional theory describes it, *J. Phys.: Condens. Matter* **29**, 423002 (2017).

[66] A. Gupta, A. Govind Rajan, E. A. Carter, and H. A. Stone, Ionic layering and overcharging in electrical double layers in a Poisson-Boltzmann model, *Phys. Rev. Lett.* **125**, 188004 (2020).

[67] M. Bültmann and A. Härtel, The primitive model in classical density functional theory: Beyond the standard mean-field approximation, *J. Phys.: Condens. Matter* **34**, 235101 (2022).

[68] A. Schlaich, A. P. Dos Santos, and R. R. Netz, Simulations of nanoseparated charged surfaces reveal charge-Induced water reorientation and nonadditivity of hydration and mean-field electrostatic repulsion, *Langmuir* **35**, 551 (2019).

[69] T. Mendez-Morales, M. Burbano, M. Haefele, B. Rotenberg, and M. Salanne, Ion-ion correlations across and between electrified graphene layers, *J. Chem. Phys.* **148**, 193812 (2018).

[70] B. L. Werkhoven, J. C. Everts, S. Samin, and R. van Roij, Flow-induced surface charge heterogeneity in electrokinetics due to Stern-Layer conductance coupled to reaction kinetics, *Phys. Rev. Lett.* **120**, 264502 (2018).

[71] M. Janssen and M. Bier, Transient response of an electrolyte to a thermal quench, *Phys. Rev. E* **99**, 042136 (2019).

[72] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes: The Art of Scientific Computing*, 3rd ed. (Cambridge University Press, New York, 2007).