Asymmetric rectified electric and concentration fields in multicomponent electrolytes with surface reactions

Abstract

Recent experimental studies have utilized AC electric fields and electrochemical reactions in multicomponent electrolyte solutions to control colloidal assembly. However, theoretical investigations have thus far been limited to binary electrolytes and have overlooked the impact of electrochemical reactions. In this study, we address these limitations by analyzing a system with multicomponent electrolytes, while also relaxing the assumption of ideally blocking electrodes to capture the effect of surface electrochemical reactions. Through a regular perturbation analysis in the low-applied-potential regime, we solve the Poisson–Nernst–Planck equations and obtain effective equations for electrical potential and ion concentrations. By employing a combination of numerical and analytical calculations, our analysis reveals a significant finding: electrochemical reactions alone can generate asymmetric rectified electric fields (AREFs), i.e., time-averaged, long-range electric fields, even when the diffusivities of the ionic species are equal. This finding expands our understanding beyond the conventional notion that AREFs arise solely from diffusivity contrast. Furthermore, we demonstrate that AREFs induced by electrochemical reactions can be stronger than those resulting from asymmetric diffusivities. Additionally, we report the emergence of asymmetric rectified concentration fields (ARCFs), i.e., time-averaged, long-range concentration fields, which supports the electrodiffusiophoresis mechanism of colloidal assembly observed in experiments. We also derive analytical expressions for AREFs and ARCFs, emphasizing the role of imbalances in ionic strength and charge density, respectively, as the driving forces behind their formation. The results presented in this article advance the field of colloidal assembly and also have implications for improved understanding of electrolyte transport in electrochemical devices.

---

1 Introduction

Colloidal particles immersed in an electrolyte aggregate along a plane near an electrode when an AC field is applied^1–8 due to attractive electrohydrodynamic flows between the particles.^2,9 Intriguingly, Woehl et al.¹⁰ and Bukosky and Ristenpart¹¹ reported that the planar height at which colloids aggregate exhibits a bifurcation that depends on the electrolyte type and the frequency of the AC field. This bifurcation was particularly surprising because the particles levitate several diameters away from the electrode.¹⁰

Since this discovery, Hashemi et al.¹² demonstrated through direct numerical simulations of the Poisson–Nernst–Planck (PNP) equations for a binary electrolyte that diffusivity contrast between anions and cations induces a long-range steady electric field, also referred to as an asymmetric rectified electric field (AREF). The strength of the AREF is dependent on the diffusivity contrast and frequency,¹³ which consequently determines the electrophoretic force on the particle and the bifurcation height. This mechanism was experimentally validated in a study by Bukosky et al.¹⁴ Recently, Hashemi et al.¹⁵ have also argued that AREFs can directly impact or even dominate flows from mechanisms such as induced-charged electrophoresis.

While AREFs are able to recapitulate experimental observations, their direct numerical simulation requires high-order adaptive meshing,^12,16 which poses its own challenge. To this end, for a binary electrolyte, Hashemi et al.¹⁷ performed a regular perturbation expansion on the PNP equations in the low-applied-potential limit and showed that AREFs appear at the second order in applied potential. Balu and Khair,¹⁸ in contrast, performed a singular perturbation expansion in the thin-double-layer limit and demonstrated that AREFs are recovered at the second order in the ratio between double layer and cell length. While both of these studies have furthered our understanding of AREFs, they are limited to binary electrolytes. Moreover, the aforementioned analyses on AREFs rely on the ideally blocking electrode approximation, which can be a limiting factor.^19–21 Wang et al.²⁰ note that “the AREF theory assumes no flux for all ions at the electrodes; essentially, it does not account for faradaic reactions (electrochemistry), which will take place at frequencies below 1 kHz in water”. In their work, in addition to the AC electric field, Wang et al.²⁰ employed water splitting reactions and experimentally demonstrated that colloidal aggregation occurs at the location where the pH of the solution is at its maximum. A similar influence of pH on colloidal aggregation was also reported by Rath et al.,²² where the authors employed the electroreduction of para-benzoquinone while also including a variable steady component of the applied potential.

Given the increasing interest in utilizing electrochemical reactions for manipulating colloidal assembly, we generalize the regular perturbation analysis of Hashemi et al.¹⁵ in the small potential limit to multicomponent electrolytes while also relaxing the ideally blocking electrode assumption. We find that AREFs can also be induced solely through electrochemical reactions, even for symmetric diffusivities, and can be stronger than AREFs created by diffusivity contrast alone. This demonstrates that AREFs may be present in a wider parameter space than previously anticipated. In addition to AREFs, we report the formation of asymmetric rectified concentration fields (ARCFs). While AREFs induce an electrophoretic force, ARCFs induce a diffusiophoretic (or osmotic, depending on the definition) force.^23–32 We discover that ARCFs are primarily observed in systems with diffusivity contrast and that electrochemical reactions do not produce ARCFs, but can enhance ARCFs caused by diffusivity asymmetry.

The simultaneous inclusion of AREFs and ARCFs could rationalize recent experimental findings. For instance, the colloidal aggregation reported in Wang et al.²⁰ closely resembles the diffusiophoretic focusing for an acid–base reaction reported in Shi et al.³³ and Banerjee et al.,²⁷ but as their system also includes an imposed electric field, the authors invoked the phenomena of electrodiffusiophoresis. Our findings also provide mechanistic insights into the formation of AREFs and ARCFs. Specifically, we highlight that the imbalances in ionic strength and charge density lead to AREFs and ARCFs, respectively. We also provide convenient analytical expressions for AREFs and ARCFs, which although valid only in the limiting case of small applied potentials and thin electrical double layers, provide a good starting point for estimating their spatial variations.

We provide a simplified model for ease of understanding in Section 2. We outline the problem formulation in Section 3. We perform a regular perturbation in the low-applied-potential limit and derive effective equations for AREFs and ARCFs; see Sections 4 and 5. Next, in Section 6, we validate our numerical results with analytical calculations and report the dependencies of AREFs on various parameters. We briefly discuss the factors which control the strength of ARCFs. Finally, we describe the limitations of our analysis, outline potential future research directions, and discuss the implications of our findings on colloidal assembly and electrochemical devices.

2 Toy model

Before delving into the details of the electrokinetic equations, inspired by Hashemi et al.,^12,15 we propose a toy model that is able to capture the effect of surface reactions. At the outset, we would like to clarify that the toy model described here ignores several complexities that are present in a real system. However, it provides a convenient choice to grasp the dominant physics within the system.

We consider a system with a cation and an anion. Both the ions can move in the z-direction due to an applied sinusoidal electric field E = E₀ cos(ωt), where E₀ is the amplitude, ω is the frequency and t is time. In addition to the electric field, the cation also moves in the z-direction due to the presence of a redox flux field, denoted here as N = N₀ cos(ωt) = g₁E₀ cos(ωt); where g₁ is a constant. This assumption assumes that the redox flux is proportional to the electric field, which is valid for small amplitude oscillations.⁶ Note that the redox flux is not consuming/producing the cation but is rather inducing a velocity on the cation; see Fig. 1. The two ions have valences of +1 and −1 and their locations are denoted z₊ and z₋, respectively. It is assumed that both the ions are at the location z_± = 0 at t = 0. The cation and anion have diffusivities of D₊ and D₋, respectively.

Fig. 1 Schematic of the toy model. (a) A cation and an anion move in the z-direction as a response to an electric field E. The cation also moves in response to a surface redox flux N. Both E and N are sinusoidal in time. N is proportional to E. The velocity induced by the electric field is dependent on the charge and the diffusivity of the ion. The velocity induced by surface reactions on the cation is proportional to the strength of the surface reaction. (b) The position of the cation (z₊) and anion (z₋) as a function of time t. When the diffusivities of each ion are equal and there are no reactions, the movement of ions has equal amplitude and thus no AREF is formed. If the diffusivities are unequal or surface reactions are present, the amplitudes are no longer equal and an AREF is formed.

The applied electric field is known to create an electromigrative flux, and the induced velocities for the cation and anion are given by where e is the charge on an electron, k_B is Boltzmann's constant, and T is the temperature. The velocity induced on the cation due to the reactive flux is (obtained by equating the convective flux to the redox flux), where C is the concentration scale of the cations. This implies that one can write

which upon integration yield where and . By assuming the ions are point charges, the net electric field induced, E_induced, by the ions at a location z can be calculated by applying Coulomb's law. For |z| ≫ z_± and time averaging, it is straightforward to obtain

〈E_induced〉 ∝ E₀²(m₊² − m₋²), (1c)

where 〈 〉 represents time averaging and where we have ignored the higher order terms beyond the second order in .

Clearly, as per eqn (1c), if g₁ = 0 and D₊ = D₋, the induced electric field vanishes. If g₁ = 0 and D₊ ≠ D₋, the induced electric field is an AREF and falls under the scenario described by Hashemi et al.^12,15 However, even if D₊ = D₋, an AREF is also possible when g₁ ≠ 0; see Fig. 1. This is the key finding that is explored in this paper, as a surface redox flux can also produce AREFs without the requirement of asymmetric diffusivities by enhancing the effective mobility of one ion. We reiterate that the toy model has its limitations, as it is not able to capture the subtle features that we discuss in the remainder of this manuscript. Before proceeding further, we would like to emphasize that the mathematical definitions of symbols in Section 2 are not applicable for the remainder of the article. From the next section, the symbols will be defined as they are introduced and will remain valid for the rest of the article.

3 Problem setup

3.1 Dimensional problem

We study a one-dimensional electrochemical system with an arbitrary number of ions and two electrodes separated by a distance 2L; see Fig. 2. An AC field of frequency Ω is applied to the electrochemical cell. X = 0 is at the centerline of the cell, [small script l] is the length of the concentration boundary layer (conc. BL) and λ is a measure of the length of the electrical double layer (EDL). Here, we investigate the formation of AREFs and ARCFs in the presence of surface reactions without imposing any restrictions on ionic diffusivities.

Fig. 2 Schematic of the model problem. (a) We consider a cell of length 2L with an arbitrary number of ions. X is the spatial coordinate and J₀sin(Ωτ) is the charge flux due to surface reactions. (b) Zoomed-in schematic of the dashed box. The cell consists of three spatial regions: the electrical double layer (EDL), a concentration boundary layer (conc. BL), and the bulk. The thickness of the EDL is denoted by λ and the thickness of the conc. BL is denoted by [small script l]. is the diffusivity of the ith ion, is a characteristic diffusivity corresponding to the conc. BL length, and Ω is the frequency of the applied field. We show that asymmetry and surface reactions can cause AREFs, and both attributes do so due to an imbalance in ionic strength. We also show that asymmetry can create an ARCF due to an imbalance in charge, and surface reactions can further enhance them.

We seek to describe the spatial and temporal variations of ionic concentrations and potential in the system to subsequently determine the AREF and ARCF. Therefore, we write Poisson's equation^34–36

with ε being the electrical permittivity of the solvent, Φ being the potential, X being the spatial coordinate, and Q_e being the volumetric charge density. where e is the charge on an electron and z_i and C_i are the valence and concentration of the ith ion, respectively.

Ion transport is modeled using the Nernst–Planck equations,^34–36i.e.,

with N_i as the flux of ion i, given bywhere τ is time, is the diffusivity of the ith ion, k_B is Boltzmann's constant, and T is temperature. We note that eqn (2b) ignores any volumetric reactions. The charge flux (or current per unit area) is evaluated as .

Eqn (2) are subjected to the following boundary and initial conditions. The sinusoidal potential boundary conditions are given as

Φ(±L,τ) = ±Φ_D sin(Ωτ). (3a)

The above equation ignores any native zeta potential of the electrodes, similar to Hashemi et al.¹⁷ We consider a surface reactive flux condition (i.e., non-ideally blocking) at the two electrodes

N_i(±L,τ) = N_i0 sin(Ωτ). (3b)

We note that the flux amplitude N_i0 may not be identical at the two electrodes,³⁷ but is assumed to be the same and time-independent for simplicity. Typically, N_i0 is dependent on applied potential and ion concentrations. We would like to clarify that that the applied flux is dependent on frequency and has the same sinusoidal dependency as the potential, i.e., it is assumed that the potential and fluxes are in-phase. The dependency of the amplitude N_i0 on Φ_D is discussed in Sections 4 and 6.4. We also define current amplitude see Fig. 2.

At τ = 0, the electrical potential is

Φ(X,0) = 0, (3c)

and the concentrations are given by

C_i(X,0) = C_i0. (3d)

The initial concentrations are required to maintain electroneutrality, i.e., .

3.2 Dimensionless problem

We write dimensionless concentration diffusivity time potential spatial coordinate charge density charge flux species fluxes and frequency where C* is a reference concentration and is a reference diffusivity. We also define as a representative measure of double-layer thickness and . We would like to clarify that λ is not the true Debye length, as it is based on a reference value C* and not the ionic strength. We make this choice on purpose to decouple the effects of ionic valences and κ. As we will discuss later, the true Debye length is given by a combination of λ and ionic strength.

In dimensionless variables, eqn (2) take the form

and eqn (3) become

ϕ(±1,t) = ±ϕ_Dsin(ωt), (5a)

n_i(±1,t) = n_i0sin(ωt), (5b)

ϕ(x,0) = 0, (5c)

c_i(x,0) = c_i0, (5d)

where , , , and , while we require to maintain electroneutrality.

4 Asymptotic solution for small applied potentials

In this work, we employ a regular perturbation expansion in the low-applied-potential limit, i.e. ϕ_D ≪ 1. While this limit is not directly observed in experiments (which generally tend to operate in moderate to large potential limits), it is able to capture the essential physics of the electrokinetic problem, albeit qualitatively.^7,11,14 This limit is also a common choice for theoretical developments.^17,38–43

The perturbation expansions in powers of ϕ_D are ϕ = ϕ⁽⁰⁾ + ϕ_Dϕ⁽¹⁾ + ϕ_D²ϕ⁽²⁾ + O(ϕ_D³), c_i = c⁽⁰⁾_i + ϕ_Dc⁽¹⁾_i + ϕ_D²c⁽²⁾_i + O(ϕ_D³), n_i = n⁽⁰⁾_i + ϕ_Dn⁽¹⁾_i + ϕ_D²n⁽²⁾_i + O(ϕ_D³), j = j⁽⁰⁾ + ϕ_Dj⁽¹⁾ + ϕ_D²j⁽²⁾ + O(ϕ_D³), and ρ_e = ρ⁽⁰⁾_e + ϕ_Dρ⁽¹⁾_e + ϕ_D²ρ⁽²⁾_e + O(ϕ_D³), where the superscripts (0), (1), and (2) refer to the leading-, first-, and second-order terms, respectively.

In the small potential limit, we assume that the equilibrium cell potential is 0 and invoke the linearized Butler–Volmer kinetic equation to write n_i0 = ϕ_Dn⁽¹⁾_i0. Thus, neglecting the equilibrium potential implies that the applied potential is the overpotential and the fluxes (and consequently the current) are proportional to the overpotential. This relationship was systematically derived by Prieve et al.⁶ We acknowledge that this assumption ignores the impact of equilibrium cell potential^35,36 and also neglects higher-order effects. These effects can become important in experimental systems^20–22 where the voltage amplitude or the steady bias in the applied potential impacts the reaction rate, boundary conditions, and leading order solution. Therefore, the analysis presented here will need to be adjusted to incorporate these effects. We outline the modifications required to incorporate these effects in Section 6.4.

4.1 Leading order

The leading-order limit refers to the condition of no applied potential, and thus the leading-order solutions are set by the initial conditions

ϕ⁽⁰⁾(x,t) = 0, (6a)

c⁽⁰⁾_i(x,t) = c_i0. (6b)

One can also verify that the solution above satisfies the governing equations and boundary conditions at the leading order.

4.2 First order

At the first order, eqn (4a) takes the formFor the ith ion, eqn (4b) reduces toIn order to separate temporal and spatial variables, we consider solutions to eqn (7a) and (7b) of the forms and c⁽¹⁾_i(x,t) = Im[e^iωtĉ⁽¹⁾_i(x)]. The governing equations with only spatial dependency becomeFollowing a similar process for the boundary conditions listed in eqn (5), we find the boundary conditions for eqn (8) areEqn (8) and (9) enable the determination of ĉ⁽¹⁾_i and . Since the variables are periodic at this order, i.e., the average of e^iωt is 0, neither an AREF nor an ARCF are observed. Thus, we examine the second order.

From the results of eqn (8), we write salt concentration or salt charge density ionic strength and electric field . These variables are employed at the second order.

4.3 Second order

We time average the governing equations at the second order over one period of the applied potential such that eqn (4a) reads We add eqn (4b) for all ions and time average to getWe multiply eqn (4b) by z_i, time average, and sum the equations to obtainVariables with bar are complex conjugates of variables with hat, and 〈 〉 corresponds to time-averaged variables. Note that eqn (10) are sufficient to determine the presence and forms of the AREF and ARCF. Further, we determine that diffusivity has no explicit effect on the second order time-averaged results, though D_i indirectly influences the first-order variables.

The boundary conditions given in eqn (5) at the second order become

〈ϕ⁽²⁾〉|_x=±1 = 0. (11c)

We integrate eqn (10b) and (10c) with boundary conditions in eqn (11a) and (11b) to writewhere we have also utilized eqn (10a). Note that we define the ARCF as the salt gradient . To integrate 〈s⁽²⁾〉 using eqn (12a), since the flux of salt is zero at both boundaries at the second order, can be used as a boundary condition.

Similarly, eqn (12b) is a third-order equation in 〈ϕ⁽²⁾〉, but we only have two boundary conditions in eqn (11c). To find the third boundary condition,³⁷ we note that since the flux of charge is zero at both boundaries at the second order, . By substituting eqn (10a) in the aforementioned condition, we obtain

Eqn (12b) can thus be solved with eqn (11c) and (13) as boundary conditions.

Experimentally,^{8,10,11,14,20,22} the relevant limit is thin EDLs, or κ ≫ 1. In this limit, a singular perturbation is required where solutions are divided into the EDL regions and the region outside the EDLs.^18,37,44 However, since AREFs occur outside the EDLs, one can simplify eqn (12b) in the limit κ ≫ 1 to directly write

where . We note that eqn (14) is valid only outside the EDL regions, and thus directly predicts the AREF. Eqn (14) highlights that the presence of an AREF is only dependent on the first-order ionic strength and the first-order electric field. Even further, it is known that as Ē⁽¹⁾ ≠ 0 in the concentration boundary layer due to the asymmetric boundary conditions (see eqn (9) and the discussion in Section 6.1.1), 〈E⁽²⁾〉 will be nonzero when Î⁽¹⁾ ≠ 0. Physically, this implies that an imbalance in ionic strength outside the EDL regions creates an AREF. This is a crucial physical insight that our analysis reveals. We would like to emphasize that this requirement is true for an arbitrary number of ions without any restriction on valences and diffusivities, albeit within the limits of thin double layers and small potentials.

We provide a brief physical explanation regarding the requirement of Î⁽¹⁾ ≠ 0 to create an AREF. Physically, EDL charging and/or redox reactions produce electric currents which subsequently induce a net electric field in the regions outside the EDLs. Still, electroneutrality is required to hold in the regions excluding the EDLs. Therefore, outside the EDLs, the only possible charge flux is an electromigrative flux, which has a magnitude that is directly dependent on the local ionic strength. This induced imbalance in ionic strength produces an asymmetry in the local conductivity of the electrolyte, which results in a time-averaged charge flux. An AREF forms to balance this charge flux.

Eqn (12a) shows that there can be steady salt gradients at the second order. Since salt gradients give rise to diffusiophoresis,^{23–25,27–29,32} the salt gradient or ARCF can simply be evaluated by utilizing eqn (12a). We emphasize that ARCFs are also a phenomenon that occurs outside EDLs and are induced by an imbalance in the first-order charge density; see eqn (12b). At this point, we note that diffusiophoretic mobility may vary for different combinations of ions,²⁸ and calculating a steady concentration field for each ion may be important in some scenarios. These calculations are straightforward, but we do not detail them for brevity.

We summarize the procedure for estimating AREFs and ARCFs in the form of a flowchart in Fig. 3. Ion valences, diffusivities, charge fluxes due to surface reactions, relative double-layer thickness, and frequency are the required inputs of the first-order solution. At the first order, the Poisson eqn (4a) and the Nernst–Planck (NP) equation eqn (4b) are solved together to obtain ion concentrations and potential. They are subsequently used to determine the first-order charge density and ionic strength. Once the first-order charge density, potential (and by extension electric field), and ionic strength are determined, the AREF is evaluated using eqn (14). A nonzero Î⁽¹⁾ outside of the EDLs indicates a nonzero 〈E⁽²⁾〉, and thus indicates the presence of an AREF. Similarly, first-order charge density and potential are used to determine the ARCF using eqn (12a). A nonzero [small rho, Greek, circumflex]⁽¹⁾_e outside the EDLs is a requirement for an ARCF.

Fig. 3 Methodology for calculating the AREF and ARCF analytically in low-potential and thin-double-layer limits for multi-ion electrolytes. The key finding is that a nonzero first-order ionic strength is a sufficient criterion to observe an AREF and a nonzero first-order charge density outside EDLs is a sufficient condition to observe an ARCF.

4.4 Numerical solution procedure

We first solve eqn (8) with boundary conditions in eqn (9) using the bvp4c functionality in MATLAB. D_i, z_i, c_i0, n⁽¹⁾_i0, κ and ω are used as inputs to our code. bvp4c deploys an adaptive grid meshing in the spatial dimension. We benchmark these numerical results with the results from Hashemi et al.¹⁷ for a binary electrolyte with no reactions and asymmetric diffusivities.

Next, we also use bvp4c to solve eqn (12b) with boundary conditions set by eqn (11c) and (13). These equations require the additional input of Î⁽¹⁾Ē⁽¹⁾, which we calculate from the first-order equations. Similarly, we solve eqn (12a) with the salt conservation boundary condition. This calculation requires the input of [small rho, Greek, circumflex]⁽¹⁾_eĒ⁽¹⁾, which we obtain from the first-order equations. We benchmark our results with the binary electrolyte results from Hashemi et al.¹⁷ Comparisons between numerical and analytical results are shown in Fig. 4 and 6, and our adaptation of the analytics from Hashemi et al.¹⁷ are given in the Appendix.

Fig. 4 Analysis of the AREF for a binary electrolyte. (a) 〈E⁽²⁾〉 vs. x for surface reactions (pink, left y-axis) and diffusivity asymmetry (blue, right y-axis) cases. Dependency of the maximum value of the AREF, i.e., 〈E⁽²⁾_max〉 = AREF_max on (b) κ with ω = 100 and (c) ω with κ = 100. The solid lines represent analytical solutions and the orbs are numerical calculations. The pink color represents surface reactions, and the blue color represents diffusivity asymmetry. All cases have ions with valences of z₁ = 1 and z₂ = −1, respectively. Surface reactions have j⁽¹⁾₀ = −0.5, D₁ = 1, and D₂ = 1, while diffusivity asymmetry has j⁽¹⁾₀ = 0, D₁ = 2, and D₂ = 1.

5 Analytical solution for symmetric diffusivities and thin double layers

To make analytical progress, we assume equal diffusivities of all ions, i.e., D_i = 1, and the thin-double-layer limit, i.e., κ ≫ 1. Eqn (4) at first order in ϕ_D for the variables [small rho, Greek, circumflex]⁽¹⁾_e, Î⁽¹⁾, and becomeThe corresponding boundary conditions for these equations are given byWe require asymmetric solutions to these equations since has asymmetric boundary conditions. We begin by finding an expression for [small rho, Greek, circumflex]⁽¹⁾_e. As a shorthand notation, let and . Then, by substituting eqn (15a) into eqn (15b) and integrating, we determine

[small rho, Greek, circumflex]⁽¹⁾_e = B sinh(λ₂x), (17)

where B is an unknown constant. Eqn (17) is utilized in eqn (15a) and integrated twice to findwhere we have already employed the boundary conditions in eqn (16a). With functional forms for both and [small rho, Greek, circumflex]⁽¹⁾_e, we use eqn (16b) to showwhere . With the full forms of and [small rho, Greek, circumflex]⁽¹⁾_e determined, we solve for Î⁽¹⁾. Based upon eqn (15c), we observe that we will require a particular and a homogeneous solution to determine Î⁽¹⁾ due to the inhomogeneity brought about by the electromigrative term. Rewriting eqn (15a) in the form of an operator on the left side, we showNote that from eqn (15a) and (15b), . This means that we can apply the operator to both sides of eqn (19) to reach a homogeneous equation. Taking only the asymmetric solutions for Î⁽¹⁾, we arrive at

Î⁽¹⁾ = F sinh(λ₁x) + G sinh(λ₂x), (20a)

where F and G are constants that need to be determined. We take G such that it cancels out the inhomogeneity in eqn (19), which results in

Next, we apply the boundary conditions in eqn (16c) and determine

With [small rho, Greek, circumflex]⁽¹⁾_e and Î⁽¹⁾ fully determined, we invoke eqn (14) to obtain an expression for the AREF. We also invoke eqn (12a) to calculate the ARCF.

The analytical expressions obtained shed light on the physics of the EDLs and concentration boundary layers. In the experimentally relevant limit of thin double layers relative to the length of the concentration boundary layer, i.e., ω/(I⁽⁰⁾κ²) ≪ 1, the charge density obtained from eqn (17) is entirely screened over the EDL regions. However, there are salt dynamics introduced by the surface reactions, resulting in an ionic strength imbalance over the concentration BL regions, as seen from the homogeneous solution in eqn (20a). AREFs, therefore, result from the simultaneous presence of first-order electric field and ionic strength. Local maxima of the components of ionic strength are found over the concentration boundary layer dimensionless length scale ω^−1/2; see Fig. 2(b). Simultaneously, surface reactions, an AC field, or both effects can produce a residual electric field in the concentration boundary layers. In fact, it can be shown by integrating eqn (15a) and (15b) over the EDLs that in order for the combination of diffusive and electromigrative fluxes to match the surface charge flux and the rate of change of accumulated charge in the EDLs, there must be a homogeneous residual electric field in the concentration BLs. As seen from eqn (18a), this residual field is given by which leads to the conclusion that only a nonzero Î⁽¹⁾ is requisite to lead to an AREF. Additionally, we note that the charge density is fully screened over the EDLs in this scenario, meaning no ARCF will develop (see also Section 6.2).

6 Results and discussion

6.1 Asymmetric rectified electric fields (AREFs)

6.1.1 Binary electrolyte. We first analyze the formation of an AREF due to surface reactions and compare it with the previously known requirement of diffusivity asymmetry to produce AREFs.^12,17,18 We focus on the limit of κ ≫ 1 such that eqn (14) is valid. We begin by comparing the two mechanisms for a monovalent binary electrolyte. Fig. 4 displays a comparison between numerical (orbs) and analytical (solid lines) results between the two mechanisms for AREF formation, i.e. diffusivity asymmetry (blue) and surface reactions (pink). Numerical results are calculated as per the procedure outlined in Section 4.4. Analytical results are determined by the approach described in Section 5. Results by Hashemi et al.¹⁷ are employed for the diffusivity asymmetry case to obtain first-order results which are subsequently utilized in eqn (14) to obtain the AREF; see the Appendix for details. The surface reaction case is presented for j⁽¹⁾₀ = −0.5 and D₁ = D₂ = 1, while the diffusivity asymmetry case is presented for j⁽¹⁾₀ = 0, D₁ = 2 and D₂ = 1. We ignore the spatial regions located between x = −1 to x = −0.9 and x = 0.9 to x = 1 to focus on the regions outside the EDLs. As evident from Fig. 4, we obtain excellent quantitative agreement between analytical and numerical results in all scenarios and for either mechanism.

First, we note that an AREF is present even with symmetric diffusivities due to the presence of surface reactions, indicating a wider parameter space that can give rise to AREFs than was previously anticipated. We observe that both the shape and the magnitude of the AREF are different for the two formation mechanisms employed; see Fig. 4(a). Specifically, for the parameters chosen, reaction-driven AREFs display a maximum near the left electrode, while AREFs produced by diffusivity asymmetry possess a maximum near the right electrode. The magnitude of the AREF maxima with surface reactions is greater than that of the AREF maxima with asymmetric diffusivities by roughly a factor of 5.

We discuss the dependency of the maximum value of the AREF, denoted here by AREF_max, with κ in Fig. 4(b). For both mechanisms, AREF_max decreases in magnitude as κ increases. However, the decrease observed with surface reactions is significantly lower than the decrease observed with diffusivity asymmetry; we find that AREFs with diffusivity asymmetry decay as κ⁻², consistent with Balu and Khair.¹⁸ Physically, AREF formation due to diffusivity asymmetry is driven by the currents arising from the EDLs. An increase in κ reduces the volume of charge (and the current) in the EDLs. Therefore, AREF_max also decreases. On the other hand, the current produced by surface reactions is not directly related to κ, leading to a weaker dependence of AREF_max on κ. The variation of AREF_max with ω is given in Fig. 4(c). The surface reactions AREF_max is insensitive to the change in ω, while the diffusivity asymmetry AREF_max increases with an increase in ω. For ω = 100, AREF_max with reactions is still greater than with D_i asymmetry by roughly a factor of 5.

The results outlined in Fig. 4(b and c) demonstrate that if surface reactions are present, AREFs could be stronger than the AREFs created by diffusivity asymmetry alone, at least for the parameters chosen. To better understand the dependencies of AREF_max on κ and ω, we employ numerical results to expand our parameter sweep to κ = 10²–10⁴ and ω = 10–10³ in Fig. 5. Fig. 5(a) shows the results for AREFs driven by surface reactions. The largest values of AREF_max occur around κ = 100 and ω = 100. Additionally, the different contours shown are all within an order of magnitude of one another, indicating that AREF_max is weakly dependent on κ and ω. This weak dependency for a constant reactive charge flux can be understood through the dependency of first-order ionic strength. A nonzero Î⁽¹⁾ with surface reactions is impacted strongly by the boundary conditions n⁽¹⁾_i0 and weakly due to the effects of κ and ω, as shown analytically in eqn (18b) and (20c). Furthermore, our analysis shows that Ē⁽¹⁾ also has a strong dependency on the surface reactive flux and a weak dependency on the frequency of the applied field. These coupled dependencies directly lead to the non-monotonic behavior observed in Fig. 5(a), and also explain why AREF_max values are within an order of magnitude of one another for a wide range of κ and ω values.

Fig. 5 Contour plot of the maximum AREF for a binary electrolyte for (a) surface reactions only and (b) diffusivity asymmetry only. Panel (a) is simulated with j⁽¹⁾₀ = −0.5, D₁ = 1, and D₂ = 1, while panel (b) is simulated with D₁ = 2, j⁽¹⁾₀ = 0, and D₂ = 1. Both cases have two ions with valences of z₁ = 1 and z₂ = −1, respectively. Variation of AREF_max for the surface reactions case alone is less sensitive to changes in κ and ω than the diffusivity asymmetry case alone.

In contrast, Fig. 5(b) shows the contours of AREF_max with κ and ω values for AREFs caused by diffusivity asymmetry. Here, a monotonic increase in AREF_max with an increase in ω and a decrease in κ are observed. Unlike the surface reactions case, the differences in AREF_max are on the scale of several orders of magnitude between different contours. This indicates a strong dependency of the AREF on both κ and ω.

The results outlined in Fig. 4 and 5 emphasize that the behavior of AREFs with surface reactions is different from the behavior of AREFs with asymmetric diffusivities alone. Specifically, we find that the magnitude of an AREF tends to be larger with surface reactions for the range of parameters explored. In fact, the maximum values are also insensitive to κ and ω, meaning surface reactions are an important mechanism to tune AREFs. Next, we discuss AREFs in the presence of more than two ions.

6.1.2 Three ions. To move beyond the comparisons outlined in the prior subsection, we investigate the parameter dependencies of AREF_max and its location, denoted here by x_c, for a solution with three ions. We note that x_c is a crucial physical parameter that has been previously observed in experiments.^{7,8,11,14,20,22}

First, we compare the values of AREF_max obtained from both numerical simulations and analytical calculations for a three-ion cell; see Fig. 6(a). To focus on the effect of surface reactions, we keep diffusivities constant, or D₁ = D₂ = D₃ = 1. We perform a comprehensive sweep of parameters including ω, κ, z_i and j⁽¹⁾₀. Excellent quantitative agreement is observed between simulations and analytical solutions obtained over the entire space of parameters. Since the dependency of AREF_max on κ and ω was discussed previously, we focus this discussion on the dependency of AREF_max on the other listed parameters. We find that AREF_max increases with an increase in j⁽¹⁾₀. This trend is expected since increasing the value of j⁽¹⁾₀ leads to a larger Î⁽¹⁾, which dictates the strength of the AREF. The changes in valence, in contrast, result in a non-monotonic behavior. The non-monotonic behavior of z_i can be explained from the trend of κ. We recall that κ⁻¹ is not the true Debye length, but instead a measure of Debye length in the system. This choice was made out of mathematical convenience; see Section 3.2. As such, the green arrow in Fig. 6(a) starts at z₁ = −3 and ends at z₁ = 3. Therefore, the true Debye length first increases and then decreases. In effect, the trend should follow a decrease in κ (i.e., opposite to the blue arrow) and then an increase in κ (i.e., in the direction of the blue arrow). This is consistent with the observed behavior of changes in ionic valence.

Fig. 6 Comparison of data for analytical calculations and simulations for an electrolyte solution with three ions for (a) the maximum value of AREF, i.e., AREF_max and (b) the location of the maximum AREF near the left electrode, i.e., x_c. In both panels, j⁽¹⁾₀ = −0.01 to −0.85 (orange), κ = 50 to 400 (blue), ω = 10 to 100 (pink), and both z₁ = −3 to 3 and z₂ = −2, −3 (green). When not varied, κ = 100, ω = 100, j⁽¹⁾₀ = −0.5, z₁ = 1, z₂ = −1, and z₃ = 1. In all cases, we find that the results from analytics and numerical simulations collapse onto the diagonal, indicating strong agreement between the results of the two methods.

We now focus on x_c closest to the left electrode; see Fig. 6(b). We observe strong quantitative agreement between simulations and analytical calculations. An increase in j⁽¹⁾₀ magnitude moves the location of the maximum further from the electrode. In contrast, an increase in ω moves the location of the maximum closer to that of the electrode. Surprisingly, the location of the maximum yields a non-monotonic behavior with κ. We investigate this in more detail in Fig. 7.

Fig. 7 Dependency of the maximum AREF location on κ and j⁽¹⁾₀ for a three-ion electrolyte. κ vs. x_c for j⁽¹⁾₀ = −0.005, −0.025, −0.01, −0.25, −0.5, −0.75. For all cases, ω = 100, z₁ = 1, z₂ = −1, z₃ = 1, and D₁ = D₂ = D₃ = 1.

To probe this non-monotonic behavior, Fig. 7 shows the location of the maximum with κ for six different values of j⁽¹⁾₀ ranging from −0.005 to −0.75. For small j⁽¹⁾₀ values, we observe that the location moves towards the electrode with an increase in κ. However, as j⁽¹⁾₀ grows, this behavior becomes more non-monotonic, suggesting that competing effects are present in the system. In fact, when κ values are small, the different curves of j⁽¹⁾₀ appear to converge. We believe the non-monotonic behavior present is due to κ driving movement of x_c towards the electrode whereas j⁽¹⁾₀ is driving x_c away from the electrode.

The results described in Fig. 6 demonstrate that the analytical procedure outlined in this manuscript is able to uncover the complex dependencies of AREFs on system parameters even for more than 2 ions. However, our analytical methodology is limited to the scenario of equal diffusivities. Asymmetric diffusivities can be included in the analytical formulation, but the relevant physics are not clearly visible with asymmetric diffusivities. As such, we solve the case with both asymmetric diffusivities and surface reactions numerically for convenience. With this in mind, we next focus on a 5-ion case with both surface reactions and diffusivity asymmetries.

6.1.3 Five-ion case. We investigate a five-ion problem (Na⁺, H⁺, Ca²⁺, Al³⁺, and Cl⁻, where the diffusivity values for each ion are taken from literature⁴⁵) such that only one of the cations is reactive; a schematic of the problem is provided in Fig. 8(a). We would like to emphasize that the setup described here is hypothetical and these surface reactions may not be necessarily observable in a real electrochemical cell. The intent of this exercise is to study the simultaneous impact of diffusivity asymmetry and surface reactions, such as the ones present in experiments.^20–22

Fig. 8 Impact of the reactive ion on the AREF in an electrolyte solution with 5 ions. (a) Schematic of the model problem with H⁺, Na⁺, Ca²⁺, Al³⁺, and Cl⁻, and a reactive charge flux j⁽¹⁾₀sin(ωt). (b) 〈E⁽²⁾〉 vs. x with H⁺ (pink), Na⁺ (blue), Ca²⁺ (black), and Al³⁺ (orange) reacting, respectively. j⁽¹⁾₀ = −0.5, κ = 100, ω = 100, c_Na0 = 1, c_H0 = 1, c_Ca0 = 1, c_Al0 = 1, and c_Cl0 = 7 in all cases. Results are obtained through numerical simulations.

We focus on a constant charge flux magnitude due to surface reactions j⁽¹⁾₀. The resulting AREFs given in Fig. 8(b) show the cases where each of the cations is reactive with j⁽¹⁾₀ = −0.5. We note that the maximum AREF values are smaller than previously discussed scenarios due to the inverse relationship between AREF_max and I⁽⁰⁾; see eqn (14).

We find that the peak AREF location and magnitude is weakly dependent upon which ion is reacting. Even though the diffusivity of the H⁺ ion is larger than the other reacting ions by approximately one order of magnitude, no appreciable change in the AREF is observed. This indicates that if ionic diffusivities are on the same order of magnitude and j⁽¹⁾₀ = O(1), estimating AREFs by assuming equal diffusivities could serve as a good starting point. This underscores the utility of the analytical procedure outlined in Section 5.

6.2 Asymmetric rectified concentration fields (ARCFs)

We investigate the formation of ARCFs in a binary electrolyte. We employ eqn (12a) to estimate the ARCF based on first-order results from (i) numerical calculations described in Section 4.4 and (ii) analytical calculations described in Section 5 and the Appendix.

We note that ARCFs form when [small rho, Greek, circumflex]⁽¹⁾_e ≠ 0 beyond the EDLs. To this end, we find that surface reactions alone are unable to produce ARCFs in the limit which is valid in experiments.^10,14,20,22 In this limit, for symmetric diffusivities, eqn (17) reduces to which suggests that the charge is only accumulated in the EDLs and ARCFs do not form. We note that Bazant et al.⁴⁶ demonstrated that a region outside the EDLs could accumulate charge for large potentials even for symmetric diffusivities, but this effect cannot be captured in our analysis and could be a potential avenue for future research.

Now, we shift our focus to the case of asymmetric diffusivities only. We show that [small rho, Greek, circumflex]⁽¹⁾_e ≠ 0 outside EDLs for D₁ ≠ D₂; see Fig. 9(a). As expected, the magnitude of [small rho, Greek, circumflex]⁽¹⁾_e increases with an increase in D₁. This charge imbalance occurs because the equations for the charge and salt are coupled for asymmetric diffusivities;^15,38 also see the Appendix. The induced ARCFs due to these charge imbalances are shown in Fig. 9(a). The ARCF is also a long-range steady field, indicating that it could be important for relatively large distances away from the electrodes. Interestingly, the spatial dependency of the ARCF remains identical even if the values of D₁ and D₂ are interchanged (results not shown), in contrast to AREFs.^12,15,18 We note that the magnitude of the ARCF is smaller compared to the AREF; see Fig. 4. While this might suggest that the role of the ARCF is minor, we note that the shape of the ARCF is different than the AREF and could thus influence the regions where the AREF is smaller. We also anticipate that ARCFs will become stronger in the nonlinear applied potential limit.⁴⁶ Furthermore, the relative importance of AREFs and ARCFs will depend on the interactions between the ionic species and particles.²⁴ This is particularly important for the phenomena of electrodiffusiophoresis, where both electrophoresis and diffusiophoresis are present.^20,21 We emphasize that while surface reactions cannot induce ARCFs on their own, they can enhance ARCFs caused by diffusivity asymmetry; see Fig. 9(b), which shows the linear dependency of ARCF_max on j⁽¹⁾₀.

Fig. 9 Presence of ARCFs. (a) ARCF vs. x for diffusivity asymmetry for D₁ = 2 (blue), D₁ = 3 (pink), D₁ = 4 (black), and D₁ = 8 (orange). D₂ = 1 in all cases. (b) ARCF_maxvs. −j⁽¹⁾₀ with both diffusivity asymmetry and surface reactions (only cation reactive). D₁ = 1 and D₂ = 4. In all cases, analytical results are solid lines and numerical simulations are orbs. κ = 100, ω = 1000, z₁ = 1, and z₂ = − 1.

6.3 Validity of the proposed framework

Since our work assumes ϕ_D ≪ 1, the results are most accurate up to an applied potential of ±25 mV. However, as shown by Balu and Khair,¹⁸ the peak value of AREFs is relatively linear up to ϕ_D = 10. This means that the results shown here may be extra-polated until applied potentials of ±250 mV with a low to moderate loss in accuracy. With even larger voltages, this framework is only suitable for a qualitative analysis. The surface reactive flux boundary conditions used in this work break down when reaction rate is no longer linear with potential; the potential value at which they break down will depend on the transfer coefficients and Stern layer thickness,⁴⁷ and we thus refrain from making a quantitative estimate of this limit.

6.4 Limitations of the proposed framework

The proposed framework demonstrates that electrochemical reactions can produce AREFs and that ARCFs are also present in the system. Further, our work elucidates that imbalances in ionic strength and charge density outside the EDLs produce AREFs and ARCFs, respectively. Nonetheless, our work has two primary limitations.

The first limitation of this work is that it assumes a small applied potential. In experiments, the voltage applied is significantly larger,^10,14,20,21 and thus the nonlinear PNP equations would need to be solved. In this scenario, typically, the thin EDL limit is invoked and a singular perturbation expansion is performed.^{18,37,44,46,48} The analysis for a binary electrolyte with asymmetric diffusivities and no electrochemical reactions has been investigated by Balu and Khair.¹⁸ Based on the trends observed in Fig. 4(b), we anticipate that the AREFs due to chemical reactions will appear at leading- and first-order expansion terms in the singular perturbation expansion, unlike asymmetric diffusivities where they are observed at the second-order,¹⁸ though a complete analysis is required to be certain. For a system with electrochemical reactions and symmetric diffusivities, the framework proposed in our prior work³⁷ can be extended.

The second limitation of our work is that it ignores the effect of equilibrium cell potential. Clearly, the surface reactions become significant beyond a certain cell voltage. To capture such an effect, Frumkin–Butler–Volmer kinetics could be invoked.^47,49 In this approach, forward and reverse reactions are written in terms of overpotential, defined as the difference between the applied potential and the equilibrium cell potential. Additionally, a Stern layer needs to be accounted for in the system to write the Frumkin–Butler–Volmer kinetic equation. There are two primary differences from the approach described in this paper. First, the leading order solution requires the formation of equilibrium double layers. Second, the boundary conditions for potential will now be written as a linear potential drop across each Stern layer. We note that our prior work³⁷ provides a methodology to capture these effects, while also simultaneously capturing the effect of EDLs and the region outside EDLs, albeit in the limit of symmetric diffusivities. To capture these effects for asymmetric diffusivities, the frameworks laid out by Balu and Khair¹⁸ and Jarvey et al.³⁷ would need to be combined, but the analysis is likely to be limited to a binary electrolyte. For multicomponent electrolytes, numerical calculations would be required.

We emphasize that while the opportunities outlined above will improve the accuracy of the current results, qualitatively, the small potential calculations are able to capture the essential physics of the system.

7 Conclusion and outlook

In this article, we performed a regular perturbation expansion in the small-applied-potential limit on the Poisson–Nernst–Planck equations for multicomponent electrolytes for arbitrary diffusivities and valences, while also including the effect of electrochemical surface reactions. Our results highlight that surface electrochemical reactions are an additional mechanism for AREF formation. We show that an imbalance in ionic strength is a prerequisite for a nonzero AREF. Further, we find that ARCFs may also be present in electrochemical cells and could induce a diffusiophoretic force on the particles. We show analytical expressions for AREFs and ARCFs are possible by further invoking the thin-double-layer limit. AREFs caused by surface reactions are less sensitive to parameters as compared to AREFs from diffusivity contrast. Lastly, we find that ARCFs appear primarily due to diffusivity contrast, though electrochemical reactions can enhance them.

Our contributions are directly applicable for the directed assembly of colloids using electric fields, including ellipsoids,⁵⁰ colloidal dumbbells,⁵¹ colloidal dimers,⁵² dicolloids,⁵³ Janus particles,⁵⁴ patchy anisotropic microparticles,⁵⁵ and chiral clusters,⁵⁶ among others. This assembly method holds great promise in creating materials with superior optical and electrical properties,^52,55–57 such as photonic crystals, microactuators, and colloidal robots. The work presented here can help estimate the assembly of colloids in an electrochemical cell, where reactions provide an additional knob to tune colloidal assembly.^19–22

We highlight that the implications of our findings extend beyond the realm of colloids. Specifically, multicomponent electrolytes and surface reactions are often used in batteries,^58,59 faradaic desalination,^60,61 carbon dioxide reduction,⁶² hybrid capacitors,⁶³ and reversible metal electrodeposition windows.⁶⁴ While multicomponent electrolytes and surface electrochemical reactions are prevalent in experimental literature, theoretical understanding of ion transport in these systems still remains elusive. In our previous work,³⁷ we analyzed the coupled effects of electrical double layers and surface electrochemical reactions on ion transport in multicomponent electrolytes for a DC potential. This work furthers the literature on AC fields, which have important implications in electrochemical impedance spectroscopy⁴⁴ and transport in porous materials.^38–40,65 The results outlined here are also a potential avenue for improving upon the modified Donnan potential approach^47,60 for modeling faradaic capacitive deionization.⁶⁶

Conflicts of interest

There are no conflicts to declare.

Appendix: adapting results from Hashemi et al.¹⁷

At first order in ϕ_D, we directly apply the analytical results of Hashemi et al.¹⁷ Note that their derivation is only valid for binary electrolytes, and as such is used exclusively to compare with results for binary electrolytes. To match the formulation used in their work, we definewhere ion 1 is the cation and ion 2 is the anion. We then write the additional parameters in their analytical solutions such that

Δ = 1 − 4βω_h(iγ + βω_h), (22a)

z = 0.5(z₁ − z₂), (22c)

Continuing to build to the forms of the first-order solutions the authors arrive at, we writeFinally, we write the first order variables

To determine the AREF and ARCF from these first-order results, we employ eqn (14) and (12a). Note that the remainder of the analytics in Section 5 do not hold for asymmetric diffusivities, but eqn (14) and (12a) are valid for asymmetric diffusivities.

Acknowledgements

A. G. thanks the National Science Foundation (CBET – 2238412) CAREER award for financial support. Acknowledgement is made to the donors of the American Chemical Society Petroleum Research Fund for partial support of this research. N. J. thanks the ARCS Foundation Scholarship and GAANN fellowship in Soft Materials for financial support. F. H. acknowledges the Ryland Family Graduate Fellowship for financial assistance.

References

1. P. J. Sides, Langmuir, 2001, 17, 5791–5800.
2. W. Ristenpart, I. A. Aksay and D. Saville, J. Fluid Mech., 2007, 575, 83–109.
3. C. L. Wirth, R. M. Rock, P. J. Sides and D. C. Prieve, Langmuir, 2011, 27, 9781–9791.
4. C. L. Wirth, P. J. Sides and D. C. Prieve, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 032302.
5. J. D. Hoggard, P. J. Sides and D. C. Prieve, Langmuir, 2008, 24, 2977–2982.
6. D. C. Prieve, P. J. Sides and C. L. Wirth, Curr. Opin. Colloid Interface Sci., 2010, 15, 160–174.
7. T. J. Woehl, K. L. Heatley, C. S. Dutcher, N. H. Talken and W. D. Ristenpart, Langmuir, 2014, 30, 4887–4894.
8. S. Saini, S. C. Bukosky and W. D. Ristenpart, Langmuir, 2016, 32, 4210–4216.
9. W. Ristenpart, I. A. Aksay and D. Saville, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 021405.
10. T. Woehl, B. Chen, K. Heatley, N. Talken, S. Bukosky, C. Dutcher and W. Ristenpart, Phys. Rev. X, 2015, 5, 011023.
11. S. C. Bukosky and W. D. Ristenpart, Langmuir, 2015, 31, 9742–9747.
12. A. Hashemi, S. C. Bukosky, S. P. Rader, W. D. Ristenpart and G. H. Miller, Phys. Rev. Lett., 2018, 121, 185504.
13. A. Hashemi, G. H. Miller and W. D. Ristenpart, Phys. Rev. E, 2019, 99, 062603.
14. S. C. Bukosky, A. Hashemi, S. P. Rader, J. Mora, G. H. Miller and W. D. Ristenpart, Langmuir, 2019, 35, 6971–6980.
15. A. Hashemi, G. H. Miller and W. D. Ristenpart, Phys. Rev. Fluids, 2020, 5, 013702.
16. M. Mirzadeh and F. Gibou, J. Comput. Phys., 2014, 274, 633–653.
17. A. Hashemi, G. H. Miller, K. J. Bishop and W. D. Ristenpart, Soft Matter, 2020, 16, 7052–7062.
18. B. Balu and A. S. Khair, J. Eng. Math., 2021, 129, 4.
19. C. A. Silvera Batista, H. Rezvantalab, R. G. Larson and M. J. Solomon, Langmuir, 2017, 33, 10861–10867.
20. K. Wang, S. Leville, B. Behdani and C. A. S. Batista, Soft Matter, 2022, 18, 5949–5959.
21. K. Wang, B. Behdani and C. A. Silvera Batista, Langmuir, 2022, 38, 5663–5673.
22. M. Rath, J. Weaver, M. Wang and T. Woehl, Langmuir, 2021, 37, 9346–9355.
23. D. Prieve, J. Anderson, J. Ebel and M. Lowell, J. Fluid Mech., 1984, 148, 247–269.
24. J. L. Anderson, Ann. Rev. Fluid Mech., 1989, 21, 61–99.
25. S. Shin, E. Um, B. Sabass, J. T. Ault, M. Rahimi, P. B. Warren and H. A. Stone, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 257–261.
26. D. Velegol, A. Garg, R. Guha, A. Kar and M. Kumar, Soft Matter, 2016, 12, 4686–4703.
27. A. Banerjee and T. M. Squires, Sci. Adv., 2019, 5, eaax1893.
28. A. Gupta, B. Rallabandi and H. A. Stone, Phys. Rev. Fluids, 2019, 4, 043702.
29. A. Gupta, S. Shim and H. A. Stone, Soft Matter, 2020, 16, 6975–6984.
30. H. C. Chu, S. Garoff, R. D. Tilton and A. S. Khair, Soft Matter, 2022, 18, 1896–1910.
31. A. Ganguly and A. Gupta, Phys. Rev. Fluids, 2023, 8, 014103.
32. R. R. Raj, C. W. Shields and A. Gupta, Soft Matter, 2023, 16, 892–904.
33. N. Shi, R. Nery-Azevedo, A. I. Abdel-Fattah and T. M. Squires, Phys. Rev. Lett., 2016, 117, 258001.
34. W. Deen, Analysis of Transport Phenomena, Oxford University Press, 2012.
35. A. Bard and L. Faulkner, Electrochemical Methods: Fundamentals and Applications, John Wiley & Sons, Incorporated, 2nd edn, 2000.
36. J. Newman and K. Thomas-Alyea, Electrochemical Systems, Wiley, 2012.
37. N. Jarvey, F. Henrique and A. Gupta, J. Electrochem. Soc., 2022, 169, 093506.
38. F. Henrique, P. J. Zuk and A. Gupta, Soft Matter, 2022, 18, 198–213.
39. F. Henrique, P. J. Zuk and A. Gupta, Electrochim. Acta, 2022, 433, 141220.
40. R. de Levie, Electrochim. Acta, 1963, 8, 751–780.
41. H. J. Keh and Y. K. Wei, Langmuir, 2000, 16, 5289–5294.
42. A. Hollingsworth and D. Saville, J. Colloid Interface Sci., 2003, 257, 65–76.
43. J. Krozel and D. Saville, J. Colloid Interface Sci., 1992, 150, 365–373.
44. B. Balu and A. S. Khair, J. Electroanal. Chem., 2022, 911, 116222.
45. W. M. Haynes, CRC Handbook of Chemistry and Physics, CRC press, 2016.
46. M. Z. Bazant, K. Thornton and A. Ajdari, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 021506.
47. P. Biesheuvel, Y. Fu and M. Z. Bazant, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 061507.
48. B. Balu and A. S. Khair, Soft Matter, 2018, 14, 8267–8275.
49. M. Z. Bazant, K. T. Chu and B. J. Bayly, SIAM J. Appl. Math., 2005, 65, 1463–1484.
50. M. Mittal and E. M. Furst, Adv. Funct. Mater., 2009, 19, 3271–3278.
51. A. F. Demirörs, P. M. Johnson, C. M. van Kats, A. van Blaaderen and A. Imhof, Langmuir, 2010, 26, 14466–14471.
52. F. Ma, S. Wang, L. Smith and N. Wu, Adv. Funct. Mater., 2012, 22, 4334–4343.
53. M. M. Panczyk, J.-G. Park, N. J. Wagner and E. M. Furst, Langmuir, 2013, 29, 75–81.
54. L. Zhang and Y. Zhu, Langmuir, 2012, 28, 13201–13207.
55. C. W. Shields IV, S. Zhu, Y. Yang, B. Bharti, J. Liu, B. B. Yellen, O. D. Velev and G. P. López, Soft Matter, 2013, 9, 9219–9229.
56. F. Ma, S. Wang, D. T. Wu and N. Wu, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 6307–6312.
57. O. Kruglova, P.-J. Demeyer, K. Zhong, Y. Zhou and K. Clays, Soft Matter, 2013, 9, 9072–9087.
58. M. Li, J. Lu, Z. Chen and K. Amine, Adv. Mater., 2018, 30, 1800561.
59. S. S. Zhang, J. Power Sources, 2006, 162, 1379–1394.
60. P. M. Biesheuvel, Y. Fu and M. Z. Bazant, Russ. J. Electrochem., 2012, 48, 580–592.
61. E. T. Sayed, M. Al Radi, A. Ahmad, M. A. Abdelkareem, H. Alawadhi, M. A. Atieh and A. G. Olabi, Chemosphere, 2021, 275, 130001.
62. J. Wu, Y. Huang, W. Ye and Y. Li, Adv. Sci., 2017, 4, 1700194.
63. P. Simon and Y. Gogotsi, Nat. Mater., 2020, 19, 1151–1163.
64. X. Tao, D. Liu, J. Yu and H. Cheng, Adv. Opt. Mater., 2021, 9, 2001847.
65. R. de Levie, Electrochim. Acta, 1964, 9, 1231–1245.
66. F. He, P. Biesheuvel, M. Z. Bazant and T. A. Hatton, Water Res., 2018, 132, 282–291.

---