As shown schematically in Figure 1A, association between cations and anions can be written in terms of a reversible reaction, i.e., one anion reacts with one cation to form a neutral particle labeled as m. At equilibrium, the fraction of neutral particles is determined by association constant K through the law of mass action. In the presence of an electric field, the equilibrium distribution of ionic species depends on the local electrical potential, various forms of non-electrostatic interactions as well as the equilibrium constant.

Figure 1B sketches a cartoon representation of the ionic system considered in this work. For simplicity, we assume that all ionic species, including the neutral ion pairs, can be represented by hard spheres of equal size. While ILs confined between two parallel plates of opposite charge do not reflect the configuration of any specific electrochemical device, the symmetric one-dimensional setup allows us to capture the essential features of charging dynamics for electrolytes near an electrode surface, in particular in terms of the evolutions of the local charge density and ionic density profiles in response to an applied voltage. The relaxation process is coupled with an evolving electric field, particle transport, and association between cations and anions. A realistic model for IL charging should capture all such effects in addition to thermodynamic non-ideality arising from other forms of inter-particle interactions.

Our reaction-coupled modified PNP (RC-MPNP) equations are similar to those proposed by Suzuki who studied charge transport in weak electrolytes (Suzuki and Seki, 2018). The model ionic liquid consists of three kinds of ionic species: cations, anions, and neutral particles generated by cation-anion association. For simplicity, all these ionic species are represented by spherical particles of the same size. To describe the charging dynamics, we start with the continuity and diffusion-reaction equations ∂ ρ i ∂ t = − ∇ ⋅ J i + R i , (2) J i = − D i k B T ρ i ∇ μ i , (3) where ρ i stands for the local number density of specie i, J _i is the mass flux, μ _i is the local electrochemical potential, k _B is the Boltzmann constant, T is the absolute temperature, D _i is the diffusion constant of specie i. The source term, R _i arises from ionic association R ± = k d ρ m − k a ρ + ρ − , (4) and R m = − k d ρ m + k a ρ + ρ − , (5)

Equations 4, 5 follow the kinetics of ion association, i.e., the change in the local ion concentration reflects a balance of the forward and backward reactions. As mentioned above, the equilibrium constant is related to the rate coefficients K = k a / k d = ρ m / ( ρ + ρ − ) , (6)

We assume that the intrinsic chemical potential can be obtained by minimizing the grand potential (Detail derivations are given in Supplementary Material). Accordingly, the electrochemical potentials of the ionic species are μ ± = ± e ψ + k B T [ ln ( a 3 ρ ± ) − ln ( 1 − a 3 ρ + − a 3 ρ − − a 3 ρ m ) ] (7) μ m = k B T [ ln ( a 3 ρ m ) − ln ( 1 − a 3 ρ + − a 3 ρ − − a 3 ρ m ) ] (8) where e is the unit charge, a is the diameter of hard spheres, and ψ is the local electrical potential. The latter is solved from the Poisson equation: ∇ 2 ψ = − e ε 0 ε r ∑ i z i ρ i (9) where ε r the relative permittivity, ε ₀ = 8.854 × 10^–12 F/m is the absolute permittivity, z _i the valence of specie i. Substituting Eqs 7, 8 into Eqs. 3 leads to J ± = − D [ ± ρ ± e k B T ∇ ψ + ∇ ρ ± + a 3 ρ ± 1 − a 3 ρ + − a 3 ρ − − a 3 ρ m ∇ ( ρ + + ρ − + ρ m ) ] (10) J m = − D [ ∇ ρ m + a 3 ρ m 1 − a 3 ρ + − a 3 ρ − − a 3 ρ m ∇ ( ρ + + ρ − + ρ m ) ] . (11)

For simplicity, we assume that the diffusion coefficient and the particle size for ion pairs the same as those for ionic species.

In integration of the continuity equations, we need to solve the following equations simultaneously: ∂ 2 ϕ ∂ x ∗ 2 = − 1 2 ς 2 ( ρ + ∗ − ρ − ∗ ) ∂ ρ ± ∗ ∂ t ∗ = ς ( ∂ 2 ρ ± ∗ ∂ x ∗ 2 ± ∂ ∂ x ∗ ( ρ ± ∗ ∂ ϕ ∂ x ∗ ) + ∂ ∂ x ∗ [ η ρ ± ∗ 1 − η ( ρ + ∗ + ρ − ∗ + ρ m ∗ ) ∂ ∂ x ∗ ( ρ + ∗ + ρ − ∗ + ρ m ∗ ) ] ) + k d ∗ ρ m ∗ − k a ∗ ρ + ∗ ρ − ∗ ∂ ρ m ∗ ∂ t ∗ = ς ( ∂ 2 ρ m ∗ ∂ x ∗ 2 + ∂ ∂ x ∗ [ η ρ m ∗ 1 − η ( ρ + ∗ + ρ − ∗ + ρ m ∗ ) ∂ ∂ x ∗ ( ρ + ∗ + ρ − ∗ + ρ m ∗ ) ] ) − k d ∗ ρ m ∗ + k a ∗ ρ + ∗ ρ − ∗ . (12) where x ∗ = x / L , ϕ = ( e ψ ) / ( k B T ) , ρ i ∗ ( i = + , − , m ) = ρ i / ρ s with ρ s being the number density of cations or anions in the bulk, t ∗ = t / τ R C with τ R C = λ D L / D , λ D = ε 0 ε r k B T / 2 e 2 ρ s is the Debye screening length, η = a 3 ρ s , K ∗ = ρ s K , k d ∗ = τ R C k d and k a ∗ = τ R C ρ s k a . For later discussions, we also introduced a dimensionless scale factor, ς = λ D / L , and set k d ∗ = 1 as a unit time. With these dimensionless variables, the numerical results depend on only three basic parameters, i.e., η , ς and k a ∗ . Physically, controlling these parameters amounts to adjusting ionic density ρ s , the system size L, and association strength K.

In describing the charging dynamics, the RC-MPNP equations entail three types of parameters that are affiliated with the system characteristics, boundary conditions and initial values. In this work, we assume that cations, anions and ion pairs have the same size of a = 0.5 nm, i.e., they all occupy a single lattice site. Both cations and anions are monovalent, | z i | = 1 , and a single diffusivity is used for all particles D = 4.3 × 10^–11 m²/s. Similar parameters were used in previous studies on ILs and are in reasonable agreement with experimental results (Largeot et al., 2008; Lian et al., 2016; Yang et al., 2020a; Yang et al., 2020b; Ma et al., 2020). As ionic species are explicitly considered in this work, a fictitious dielectric constant ε _r = 2 was adopted to account for dispersion interactions and ionic polarizability of ILs not included in the coarse-grained model (Lian et al., 2018a; Lian et al., 2018b; Yang et al., 2020b). Temperature T = 293.15 K, and cell width L = 100 λ _D . Assuming that an ionic liquid has a typical bulk concentration of 3.3 mol/L, we have λ D ≈ 0.025 nm , L ≈ 2.5 nm . To investigate the ionic association effects, we consider four association constants, K* = 0.1, 1, 10, 148, representing the cases of weak, moderate, strong, and extreme association strength, respectively.

In studying the charging dynamics, we assume that a biased potential of 2 ϕ 0 is initially applied to the cathode and anode (at t = 0). Based on the one-dimensional model illustrated in Figure 1, we have the following boundary conditions: ϕ ( ± L , t > 0 ) = ± ϕ 0 ∂ ρ i ( i = + , − , m ) ∂ x ( ± L , t ) = 0 (13) where the initial particle densities are assumed identical to the corresponding bulk values, ρ + , b , ρ − , b , ρ m , b . With the assumption that the bulk ionic liquid can be represented by an ideal solution, these densities can be determined from the association constant K ρ m , b ρ + , b ρ − , b = K ρ + , b + ρ − , b + 2 ρ m , b 2 = ρ s . (14)

Because the bulk concentrations are fixed during the charging process, an explicit consideration of the thermodynamic non-ideality is equivalent to changing the bulk concentration. In the presence of an electric field, the RC-MPNP equations allow us to describe the relaxation of the ionic density profiles until a new equilibrium is established. The time range used in the calculation is from 10^–3 τ _RC to 10³ τ _RC with τ R C = λ D L / D .

The application of a biased potential introduces a local electric field driving ionic motions thereby breaking down the local ionic association equilibrium. Therefore, we expect that the evolution of the ionic density profiles is different from that corresponding to a non-reactive system. Figure 2 shows the variations of the total ionic density ( ρ + ∗ + ρ − ∗ ), the local charge density ( ρ + ∗ − ρ − ∗ ), and the density of neutral particles ( ρ m ∗ ) near the surface of the negative electrode. At the initial state, the whole system is assumed homogeneous as shown in uniform particle and charge distributions. Upon the application of an electrical potential, the counterions (cations) gradually accumulate at the electrode surface, concomitant with the receding of cations (anions). Interestingly, the total ionic density ρ ± ∗ first falls slightly and rises monotonically to the equilibrium value. Both the total ion density ( ρ + ∗ + ρ − ∗ ) and the local charge ( ρ + ∗ − ρ − ∗ ) exhibit a U-shaped charging dynamic as indicated in Figures 2J,K. Due to the volume effect caused by the accumulation of anions and anions near the electrode, the neutral particles are enriched in the middle region of the cell (Figure 2L), leading to a density profile with a convex shape across the slit as the charging continues.

Furthermore, Figure 2 shows the influence of ionic association on density distribution. It can be found that with the increase of association strength (i.e., a larger K ^*), the number densities of cations and anions are dramatically reduced while the number density of neutral particles increases. The ionic association affects the local concentrations of the three types of ionic species (cations, anions and ion pairs) hence the kinetic behavior and other physicochemical properties, such as timescales, capacitance and resistance.

In order to understand the structural relaxation under the influence of ionic association, we need to identify the decay length from the density distributions. The decay length, λ _S , may be extracted by taking the logarithm of the charge distribution (log|q*|), which yields approximately a straight line at each position as shown in Figure 3. A decay length can be obtained by averaging over the slopes of the straight lines at different positions.

Figure 3E presents the evolutions of the decay length λ _S for ionic liquids with different association constants. The decay length reaches a stable state at τ _RC approximately independent of K. Note that the dimensionless association constant is defined as K ∗ = ρ s K . The invariance of the relaxation time with the association strength conforms to the prediction of the RC transmission line theory (Lian et al., 2020). In Figure 3F, we show the decay length λ _S as a function of the association strength K* after the equilibrium is established. A larger association strength leads to a larger decay length, which can be intuitively understood in terms of the concentrations of free ions. According to Figure 3E, when K* = 0, λ _S = λ _D , the decay length λ _S in the figure is the Debye length λ _D , which characterizes the charge shielding effect. However, the decay length λ _S gradually deviates from the Debye length λ _D when K* > 0 due to the influence of association.

Alternatively, the decay length can be defined in terms of the density profiles of ionic species or neutral particles. Qualitatively, the conclusions are similar. As shown in Supplementary Figure S1, the stronger association that takes place between cations and anions, the longer decay length that we extract according to the evolution of the particle concentrations. While there is only one decay length if K = 0, The two slopes of the density profile allow us to identify two decay lengths at low association strength.

It has been well established that the PNP equations yield a quantitative relation between the dynamic timescale and the Debye length, τ R C = λ D L / D . In the presence of ionic association, the Debye length does not reflect the association effect thus the relation is no longer valid. To explore the relation between timescale τ _RC and decay length λ _S , we assume that the RC time corresponds to the charging constant of a capacitor, which describes the speed of charging process. As the local charge profile relaxes, counterions are adsorbed on the electrode surface while the coions are depleted. The Gauss law, σ ( t ) = − 2 ρ s λ D 2 ( ∂ x ϕ ) , can be used to calculate the surface charge density of the electrode. Figure 4A shows the normalized surface charge density on the electrode surface, σ ( t ) / σ e q ( σ e q ≡ σ ( t / τ R C → ∞ ) ) . As expected, the increase of the association strength expands the relaxation time. The charge relaxation 1 − σ ( t ) / σ e q (shown in Supplementary Figure S1) shows two slopes, suggesting that the relaxation process can be divided into two stages: a faster initial response followed by a slower relaxation. A time-dependent function τ(t) can be used to characterize the timescales of early and late stages of the dynamic process: τ ( t ) = − [ d ⁡ ln ( 1 − σ ( t ) / σ e q ) d t ] − 1 . (15)

Figure 4B plots τ(t) for different association strengths. In the early stage (t/τ _RC < 1), τ(t) is controlled by the EDL formation, and the characteristic time is represented by τ ₁. In the case of no association, the RC time is given by τ 1 = τ R C = λ D L / D . The timescale for the early stage of charging grows with the increase of the association strength because the ionic association leads to the reduction of the concentrations of anions and cations. As a result, association slows down the formation of EDL. Due to the ionic association, the decay length λ _S , not the Debye length λ _D , should be used in the RC model, i.e., τ 1 = λ S L / D .

Figure 4C shows a quantitative relation between the decay length λ _S and the Debye length at a fixed association strength. Empirically, the relation can be fitted by using the following analytical form (R ² > 0.99): λ S = [ 0.341 ( K ∗ ) 0.408 + 1 ] λ D . (16)

Substituting the above formula into τ 1 = λ S L / D gives: τ 1 = λ S L D = [ a ( K ∗ ) b + 1 ] λ D L D (17) where a = 0.341 and b = 0.408 are obtained by fitting the numerical results. By comparing λ S L / D and τ ₁ for different values of K*, we assert that τ 1 = λ S L / D reveals the quantitative relation between the RC time, the decay length and association strength.

In the later stage of charging (t/τ _RC > 1), τ(t) shows a step-like behavior with different plateaus τ _ad for different values of K*. In this case, τ(t) is controlled by the particle diffusion, and the dynamics is much slower than electrical motion in the early stage. The ion diffusion can be characterized by the timescale, τ 2 = τ a d = α L 2 / D . If there is no association, the coefficient is α = 0.05 . A different value should be used due to association. As shown in Figure 4D, we may derive a quantitative relation between τ ₂ and association constant K* (R ² > 0.99): τ 2 = { c ⁡ ln [ ln ( K ∗ + 1 ) ] + d } α L 2 D . (18) where c = 0.109 and d = 1.318. Eqs 16–18 provide convenient estimations of the timescales for different stages of charging under the influence of ionic association.

Qualitatively, the equilibrium and dynamic properties of ionic liquids discussed above are relatively insensitive to the bulk concentration of ionic species (Supplementary Figures S3, S4). While both the bulk concentration and association strength affect the timescales of charging and particle relaxation, it is important to recognize that ionic association is mainly responsible for the increase of decay length and relaxation time.