Breakdown of electroneutrality in polyelectrolyte gels

Abstract

Mathematical models of polyelectrolyte gels are often simplified by assuming the gel is electrically neutral. The rationale behind this assumption is that the thickness of the electric double layer (EDL) at the free surface of the gel is small compared to the size of the gel. Hence, the thin-EDL limit is taken, in which the thickness of the EDL is set to zero. Despite the widespread use of the thin-EDL limit, the solutions in the EDL are rarely computed and shown to match to the solutions for the electrically neutral bulk. The aims of this paper are to study the structure of the EDL and establish the validity of the thin-EDL limit. The model for the gel accounts for phase separation, which gives rise to diffuse interfaces with a thickness described by the Kuhn length. We show that the solutions in the EDL can only be asymptotically matched to the solutions for an electrically neutral bulk, in general, when the Debye length is much smaller than the Kuhn length. If the Debye length is similar to or larger than the Kuhn length, then phase separation can be initiated in the EDL. This phase separation spreads into the bulk of the gel and gives rise to electrically charged layers with different degrees of swelling. Thus, the thin-EDL limit and the assumption of electroneutrality only generally apply when the Debye length is much smaller than the Kuhn length.

1. Introduction

Polyelectrolyte gels are soft, electroactive materials that are used in a wealth of applications including smart materials [8, 26], fuel cells [17], gel diodes [32], regenerative medicine [19] and drug-delivery systems [20]. A polyelectrolyte gel consists of a deformable network of polymers that is swollen with fluid. The polymers carry a fixed electric charge and can therefore electrostatically interact with ions that are dissolved in the imbibing fluid. Typically, polyelectrolyte gels are surrounded by a salt bath composed of a solvent, such as water, and dissolved ions. Solvent and ion exchange across the gel–bath interface occurs until an equilibrium is established. The degree of swelling and hence the gel volume are set by this equilibrium, which can be controlled through a number of factors such as temperature and electric fields, as well as the pH and salt content of the surrounding bath [1]. Slight alterations in these environmental parameters can trigger enormous changes in the gel volume. In some cases, the volume of the gel will undergo a discontinuous change, a phenomenon that is called a volume phase transition [7, 22]. Environmental stimuli can also induce phase separation, whereby a homogeneous gel spontaneously separates into co-existing phases with different compositions [18, 27]. Phase separation has been proposed as a facile means of self-assembling nanostructures in polyelectrolyte gels [30, 31].

When a polyelectrolyte gel is surrounded by a salt solution, ions from the solution will migrate to the free surface of the gel in order to screen the electric charges on the polymers. The accumulation of ions leads to a diffuse layer of electric charge that is known as the electric double layer (EDL). The thickness of the EDL is described by the Debye length and is often on the order of tens of nanometres. An interesting feature of polyelectrolyte gels is that the EDL is diffuse on both sides of the gel–bath interface due to the mobile ions in the gel migrating to counter the accumulation of charge in the surrounding bath. Similar doubly-diffuse EDLs also appear in immiscible liquid drops that are suspended in a different ion-carrying liquid [24]. Outside of the EDL, the gel and the bath are electrically neutral to a very good approximation.

Despite the intricate structure of the EDL, it is generally believed to play a passive role in the gel dynamics. Moreover, due to the smallness of the Debye length (tens of nanometres) relative to the typical dimensions of a polyelectrolyte gel (microns to centimetres), any impact of the EDL on the gel dynamics is assumed to be confined to an extremely thin region near the free surface. Thus, the EDL is often neglected in studies that involve modelling polyelectrolyte gels [9, 10, 14, 15, 23, 34, 35]. The few exceptions include the works by Hong et al. [13, 29], who compute solutions in the EDL for a limited range of parameters. Hill [12] accounts for internal EDLs that form in polyelectrolyte gels with liquid inclusions, neglecting the elasticity of the polymer network and solvent–polymer interactions.

Neglecting the EDL due to its thinness is called the thin-EDL limit. The singular nature of the thin-EDL limit means there are two components to it. The first involves computing an outer solution that is valid in regions away from the EDL where the bath and gel are approximately electrically neutral. The second involves computing an inner solution that is valid within the EDL. Although the thin-EDL limit is used extensively when modelling polyelectrolyte gels, very little attention is paid to computing the inner solution and checking that it can, in fact, be asymptotically matched to the outer solution. For instance, Mori et al. [21] used matched asymptotic expansions to derive a model for an electrically neutral polyelectrolyte gel, but did not compute solutions in the inner and outer regions nor did they explore the structure of the EDL in detail.

The aims of this paper are to use matched asymptotic expansions to: (i) revisit the assumption that the EDL plays a passive role in the dynamics of polyelectrolyte gels and (ii) ascertain the validity of the thin-EDL limit. In particular, we explore when the outer solutions which govern the electrically neutral bulk can be asymptotically matched to the inner solutions in the EDL.

The main result of our work is that asymptotic matching of solutions cannot always be carried out because the EDL can trigger a mode of phase separation that leads to a breakdown of electroneutrality across the entire gel. In this case, the charge density in the gel oscillates in space, corresponding to the formation of alternating layers of positive and negative charge. Similar oscillations have been observed in ionic liquids in contact with a charged surface, where they are attributed to the finite size of ions and their short-range interactions [6, 11].

Our asymptotic analysis of the EDL builds on that of Yariv [33] by accounting for the nonlinear electro-chemo-mechanics of the gel. A crucial feature of our analysis is that it is based on a thermodynamically consistent phase-field model of a polyelectrolyte gel. The phase-field model introduces an additional length scale into the problem, the Kuhn length, which is proportional to the thickness of the diffuse, internal interfaces that form within the gel if it undergoes phase separation. The governing equations are fourth order in space and capture the energy cost of mixing ions with finite volume, which is similar to continuum models for the layering of ionic liquids [2]. Most models for polyelectrolyte gels do not account for phase separation and take the Kuhn length to be zero. However, we find that the thin-EDL limit is only asymptotically consistent, in general, when the Kuhn length greatly exceeds the Debye length. Thus, we argue that particular care must be taken when using the thin-EDL limit to describe the behaviour of polyelectrolyte gels.

The paper is organised as follows. In Sec. 2, the governing equations for a cylindrical polyelectrolyte gel that is in equilibrium with a salt solution are presented. In Sec. 3, we carry out the asymptotic analysis of the EDL assuming the Kuhn length is much larger than the Debye length. In Sec. 4, we discuss how the analysis differs if the Kuhn length is zero, which is more typical across the literature. The asymptotic framework is then used to investigate the structure of the EDL in Sec. 5. The paper concludes in Sec. 6.

2. Mathematical model

We consider a cylindrical polyelectrolyte gel that is in equilibrium with a stationary salt bath, as shown in Figure 1. Motivated by the experiments of Horkay et al. [14], we assume the gel can freely swell in the radial and orthoradial directions but is confined in the axial direction. The gel is composed of a deformable network of polymers that carry electric charges of the same sign. The bath consists of a solvent and a dissolved binary salt such as NaCl or CaCl $_2$ . We assume that the system remains axisymmetric and that the gel remains cylindrical; that is, we do not allow for instabilities along the axial or orthoradial directions.

Figure 1. A swollen polyelectrolyte gel surrounded by a bath. The bath consists of a solvent and a dissolved binary salt. The polymers of the gel carry an electric charge, which is assumed to be positive. An electric double layer of thickness $O(\varepsilon )$ forms near the gel–bath interface, located at $r = a$ , where charge neutrality is violated. The non-dimensional Debye length $\varepsilon$ is defined in (2.1).

A thermodynamically consistent model of a polyelectrolyte gel surrounded by a viscous bath has been derived by Celora et al. [5]. We employ this model here but specialise it to a steady, cylindrical configuration. For brevity, we only present the non-dimensional form of the governing equations in the main text; however, the dimensional model is provided in Appendix A. In the equations below, the subscript $m$ is used to represent quantities associated with the solvent ( $s$ ), cation ( $+$ ) or the anion ( $-$ ). The set $\mathbb{M} = \{s, +, -\}$ contains all of the mobile species that move into and out of the polymer network. We let $\mathbb{I} = \{+,-\}$ denote the ionic species.

In non-dimensionalising the model, spatial variables are scaled with $a_0$ , the radius of the gel in its unswollen (dry) state. The chemical potentials of the mobile species, $\mu _m$ , are written as $\mu _m = \mu _m^0 + k_B T \mu'_{\!\!m}$ , where $\mu _m^0$ is a reference chemical potential, $k_B$ is Boltzmann’s constant and $T$ is the absolute temperature. Primes are used to denote dimensionless quantities. The electric potential in the bath and the gel is scaled with the thermal voltage and written as $\Phi = (k_B T/ e) \Phi '$ , where $e$ is the elementary charge. The stresses and pressure in the gel are non-dimensionalised using the shear modulus of the polymer network, $G$ , as a scale. In the bath, the pressure is non-dimensionalised as $p = [\epsilon ^{\text{bath}} (k_B T/ e)^2/ a_0^2] p'$ , with $\epsilon ^{\text{bath}}$ denoting the electrical permittivity of the bath. The pressure scaling for the bath can be motivated by the condition of mechanical equilibrium for a motionless fluid, which demands that the fluid pressure balances the Maxwell stress, as these are the only two forces at play. Due to the diluteness of the ions, the electric permittivity of the gel and the bath, $\epsilon ^{\text{gel}}$ and $\epsilon ^{\text{bath}}$ , respectively, are treated as constants. However, the composition dependence of the electric permittivity has been shown to impact the swelling of polyelectrolyte gels [16].

This scaling introduces three key dimensionless parameters given by

(2.1) \begin{align}{\mathcal{G}} = \frac{\nu G}{k_B T}, \quad \omega = \frac{L_K}{a_0}, \quad \varepsilon = \frac{L_D}{a_0}, \end{align}

where $\nu$ is the volume of solvent molecule, $L_K$ is the Kuhn length of a polymer chain and $L_D = (\nu \epsilon ^{\text{gel}} k_B T)^{1/2}/e$ is the Debye length. The parameter $\mathcal{G}$ characterises the energetic cost of elastically deforming the gel relative to the energy that is released upon insertion of a solvent molecule into the polymer network. The parameters $\omega$ and $\varepsilon$ describe the thickness of diffuse internal interfaces and the EDL relative to the size of the gel, respectively. Alternatively, $\omega$ can be related to the energetic cost of gradients in the solvent concentration; see Celora et al. [5] for details. The magnitudes of these numbers will be estimated in Sec. 2.4.

2.1 Governing equations for the gel

The equations for the gel are formulated in terms of Eulerian coordinates $\boldsymbol{x} = r \boldsymbol{e}_r(\theta ) + z \boldsymbol{e}_z$ associated with the current state of the system, where $r$ , $\theta$ and $z$ denote the radial, angular and axial coordinates, respectively, and $\boldsymbol{e}_r$ , $\boldsymbol{e}_\theta$ and $\boldsymbol{e}_z$ are the corresponding cylindrical basis vectors. An Eulerian coordinate system enables the equations to be written in a physically intuitive way and it facilitates coupling the gel and bath models via boundary conditions. A detailed account of Eulerian-based hydrogel modelling is provided by Bertrand et al. [3]. We let $\boldsymbol{X} = R(r) \boldsymbol{E}_{R}(\theta ) + Z(z) \boldsymbol{E}_Z$ denote the Lagrangian coordinates associated with the stress-free reference configuration, which we assume is a dry gel. Here, $R$ , $\Theta = \theta$ and $Z$ represent the radial, angular, and axial Lagrangian coordinates, and $\boldsymbol{E}_R$ , $\boldsymbol{E}_{\Theta }$ and $\boldsymbol{E}_Z$ are the basis vectors.

The deformation gradient tensor $\boldsymbol{\mathsf{F}}$ describes the distortion of material elements relative to the dry state of the gel. For an axisymmetric geometry which remains cylindrical, the deformation gradient tensor can be written as $\boldsymbol{\mathsf{F}} = \lambda _r\, \boldsymbol{e}_r\otimes \boldsymbol{E}_R + \lambda _\theta \, \boldsymbol{e}_\theta \otimes \boldsymbol{E}_\Theta + \lambda _z \boldsymbol{e}_z \otimes \boldsymbol{E}_Z$ , where

(2.2) \begin{align} \lambda _r = \left (\frac{{\textrm{d}} R}{{\textrm{d}} r} \right )^{-1}, \quad \lambda _\theta = \frac{r}{R}, \quad \lambda _z = \left (\frac{{\textrm{d}} Z}{{\textrm{d}} z}\right )^{-1} \end{align}

denote the radial, orthoradial and axial stretches, respectively. The axial stretch $\lambda _z$ is imposed, whereas the radial and orthoradial stretches $\lambda _r$ and $\lambda _\theta$ are unknown and must be solved for. The determinant $J = \det \boldsymbol{\mathsf{F}} = \lambda _r \lambda _\theta \lambda _z$ characterises volumetric changes in material elements. Both the polymers and the imbibed salt solution are assumed to be incompressible. As a result, any volumetric change in a solid element must be due to a variation in the amount of fluid contained within that element. This leads to the so-called molecular incompressibility condition

(2.3) \begin{align} J = \left (1 - \sum _{m\in \mathbb{M}} \phi _m\right )^{-1}, \end{align}

where $\phi _k$ represent the volume fraction of species $k$ . Since $J$ describes the volume of swollen material elements relative to their dry volume, we also refer to it as the swelling ratio. It will be convenient to formulate the incompressibility condition as

(2.4) \begin{align} R \frac{{\textrm{d}} R}{{\textrm{d}} r} = \frac{\lambda _z r}{J}, \end{align}

where $J$ is given by (2.3).

The chemical potentials of the mobile species, $\mu _m$ , can be written as

(2.5a) \begin{align} \mu _s &= \Pi _s +{\mathcal{G}} p - \frac{\omega ^2}{r}\frac{{\textrm{d}} }{{\textrm{d}} r}\left (r \frac{{\textrm{d}} \phi _s}{{\textrm{d}} r}\right ), \end{align}

(2.5b) \begin{align} \mu _{\pm } &= \Pi _{\pm } +{\mathcal{G}} p + z_{\pm } \Phi, \end{align}

where $z_{\pm }$ is the valence of the ions, $p$ is the mechanical pressure in the gel and $\Pi _m$ are osmotic pressures defined as

(2.6a) \begin{align} \Pi _s &= \log \phi _s + \chi \,J^{-1}(1 - \phi _s) + J^{-1}, \end{align}

(2.6b) \begin{align} \Pi _\pm &= \log \phi _\pm + J^{-1}(1 - \chi \phi _s). \end{align}

Here, $\chi$ is the Flory interaction parameter, which describes (unfavourable) enthalpic interactions between the solvent molecules and the polymers. Due to our assumption that the system is in equilibrium, the chemical potentials are spatially uniform. Hence, $\mu _m$ are constants that will be specified below.

The electric potential satisfies Poisson’s equation

(2.7) \begin{align} -\frac{\varepsilon ^2}{r}\frac{{\textrm{d}} }{{\textrm{d}} r}\left (r \frac{{\textrm{d}} \Phi }{{\textrm{d}} r}\right ) = z_+\phi _+ + z_{-} \phi _- + \frac{z_f \varphi _f}{J}, \end{align}

where $\varphi _f$ is the nominal volume fraction of fixed charges on the polymer network and $z_f$ denotes the valence of these charges. We will focus on cationic gels with positive fixed charges, $z_f \gt 0$ .

The conservation of linear momentum in the gel leads to

(2.8) \begin{align} \frac{{\textrm{d}} \mathsf{T}_{rr}}{{\textrm{d}} r} + \frac{\mathsf{T}_{rr} - \mathsf{T}_{\theta \theta }}{r} = 0, \end{align}

where $\mathsf{T}_{rr}$ and $\mathsf{T}_{\theta \theta }$ are the radial and orthoradial components of the Cauchy stress tensor. These stresses can be expressed as

(2.9a) \begin{align} \mathsf{T}_{rr} &= \mathsf{T}_{e,rr} + \mathsf{T}_{K,rr} + \mathsf{T}_{M,rr} - p, \end{align}

(2.9b) \begin{align} \mathsf{T}_{\theta \theta } &= \mathsf{T}_{e,\theta \theta } + \mathsf{T}_{K,\theta \theta } + \mathsf{T}_{M,\theta \theta } - p. \end{align}

The first contributions, $\mathsf{T}_{e,rr}$ and $\mathsf{T}_{e,\theta \theta }$ , represent elastic stresses, which are calculated by assuming the polymer network behaves as a neo-Hookean material. This leads to

(2.10a) \begin{align} \mathsf{T}_{e,rr} &= J^{-1}(\lambda _r^2 - 1), \end{align}

(2.10b) \begin{align} \mathsf{T}_{e,\theta \theta } &= J^{-1}(\lambda _\theta ^2 - 1). \end{align}

The second and third contributions to the Cauchy stresses in (2.9) correspond to Korteweg ( $\mathsf{T}_K$ ) and Maxwell ( $\mathsf{T}_M$ ) stresses, respectively, which capture the forces generated within the bulk of the gel due to internal interfaces and electric fields. The radial and orthoradial components of these stresses are given by

(2.11a) \begin{align} \mathsf{T}_{K,rr} &={\mathcal{G}}^{-1} \omega ^2 \left [\frac{\phi _s}{r}\frac{{\textrm{d}} }{{\textrm{d}} r}\left (r \frac{{\textrm{d}} \phi _s}{{\textrm{d}} r}\right ) - \frac{1}{2}\left (\frac{{\textrm{d}} \phi _s}{{\textrm{d}} r}\right )^2\right ], \end{align}

(2.11b) \begin{align} \mathsf{T}_{K,\theta \theta } &={\mathcal{G}}^{-1} \omega ^2 \left [\frac{\phi _s}{r}\frac{{\textrm{d}} }{{\textrm{d}} r}\left (r \frac{{\textrm{d}} \phi _s}{{\textrm{d}} r}\right ) + \frac{1}{2}\left (\frac{{\textrm{d}} \phi _s}{{\textrm{d}} r}\right )^2\right ], \end{align}

(2.11c) \begin{align} \mathsf{T}_{M,rr} &= \frac{1}{2}{\mathcal{G}}^{-1} \varepsilon ^2\left (\frac{{\textrm{d}} \Phi }{{\textrm{d}} r}\right )^2, \end{align}

(2.11d) \begin{align} \mathsf{T}_{M,\theta \theta } &= -\frac{1}{2}{\mathcal{G}}^{-1} \varepsilon ^2\left (\frac{{\textrm{d}} \Phi }{{\textrm{d}} r}\right )^2. \end{align}

The final contribution to the Cauchy stresses represents the stress induced by the fluid pressure.

2.2 Governing equations for the bath

The spatially uniform chemical potentials of the solvent and ions are

(2.12a) \begin{align} \mu _s &= \log \phi _s + \epsilon _r \varepsilon ^2 p, \end{align}

(2.12b) \begin{align} \mu _{\pm } &= \log \phi _{\pm } + \epsilon _r \varepsilon ^2 p + z_{\pm } \Phi, \end{align}

where $\epsilon _r = \epsilon ^{\text{bath}}/ \epsilon ^{\text{gel}}$ . The electric potential satisfies

(2.13) \begin{align} -\frac{\epsilon _r \varepsilon ^2}{r}\frac{{\textrm{d}} }{{\textrm{d}} r}\left (r \frac{{\textrm{d}} \Phi }{{\textrm{d}} r}\right ) = z_+\phi _+ + z_{-} \phi _-. \end{align}

Conservation of linear momentum in the bath implies that

(2.14) \begin{align} \frac{{\textrm{d}} \mathsf{T}_{rr}}{{\textrm{d}} r} + \frac{\mathsf{T}_{rr} - \mathsf{T}_{\theta \theta }}{r} = 0, \end{align}

where the components of the Cauchy stress tensor are

(2.15a) \begin{align} \mathsf{T}_{rr} &= \mathsf{T}_{M,rr} - p, \end{align}

(2.15b) \begin{align} \mathsf{T}_{\theta \theta } &= \mathsf{T}_{M,\theta \theta } - p. \end{align}

The Maxwell stresses are given by

(2.16a) \begin{align} \mathsf{T}_{M,rr} &= \frac{1}{2}\left (\frac{{\textrm{d}} \Phi }{{\textrm{d}} r}\right )^2, \end{align}

(2.16b) \begin{align} \mathsf{T}_{M,\theta \theta } &= -\frac{1}{2}\left (\frac{{\textrm{d}} \Phi }{{\textrm{d}} r}\right )^2. \end{align}

2.3 Boundary conditions

At the centre of the gel, we impose

(2.17) \begin{align} R|_{r=0} = 0, \quad \left .\frac{{\textrm{d}} \Phi }{{\textrm{d}} r}\right |_{r=0} = 0, \quad \left .\frac{{\textrm{d}} \phi _s}{{\textrm{d}} r}\right |_{r=0} = 0. \end{align}

The first condition ensures that the origin in Lagrangian coordinates is mapped to the origin in Eulerian coordinates. The second and third can be viewed as symmetry conditions. The boundary condition on $\phi _s$ is needed due to the presence of a second derivative in the expression for the chemical potential of solvent in the gel (2.5a).

Far from the bath, $r \to \infty$ , we set the electric potential to $\Phi ^{\text{bath}}$ and the pressure to zero. In addition, we assume that the volume fraction of cations has been fixed to $\phi ^{\text{bath}}_{+}$ . Assuming electroneutrality then requires that the volume fraction of anions is given by $\phi ^{\text{bath}}_{-} = |z_{+}/ z_{-}| \phi ^{\text{bath}}_{+}$ . The far-field solvent fraction is then fixed at $\phi ^{\text{bath}}_s = 1 - (1 + |z_{+}/z_{-}|) \phi ^{\text{bath}}_{+}$ . Consequently, the far-field chemical potentials are given by

(2.18a) \begin{align} \mu ^{\text{bath}}_{s} = \log \phi ^{\text{bath}}_{s}, \quad \mu ^{\text{bath}}_{\pm} = \log \phi ^{\text{bath}}_{\pm} + z_{\pm } \Phi ^{\text{bath}}. \end{align}

The radius of the deformed gel, and hence the position of the gel–bath interface, is denoted by $a$ . We use the notation $r \to a^{\pm }$ to describe approaching the interface from the bath ( $+$ ) and gel ( $-$ ). Due to the formulation of the model in terms of Eulerian coordinates, the deformed radius of the gel, $a$ , is unknown. Since the undeformed radius of the gel has been scaled to unity, the deformed radius is implicitly determined by the equation

(2.19) \begin{align} R|_{r= a^{-}} = 1. \end{align}

At the gel–bath interface, $r = a$ , thermodynamic equilibrium demands that the chemical potentials are continuous. Moreover, since the chemical potentials must also be spatially uniform, we have that

(2.20) \begin{align} \mu _m = \mu ^{\text{bath}}_m \end{align}

in both the gel and the bath. We also impose the condition

(2.21) \begin{align} \left .\frac{{\textrm{d}} \phi _s}{{\textrm{d}} r}\right |_{r = a^{-}} = 0 \end{align}

at the gel surface, which is needed due to the second derivative in (2.5a). From a physical point of view, (2.21) implies that the solvent does not preferentially wet or dewet the interface, both of which would lead to a localised gradient in the solvent composition. The balance of radial stresses at the interface leads to

(2.22) \begin{align} \left .{\mathcal{G}} \mathsf{T}_{rr} \right |_{r = a^{-}} = \left .\epsilon _r \varepsilon ^2 \mathsf{T}_{rr}\right |_{r = a^{+}}. \end{align}

We assume there are no surface charges on the interface and therefore impose continuity of the electric potential and electric displacement:

(2.23a) \begin{align} \left .\Phi \right |_{r = a^{-}} &= \left .\Phi \right |_{r = a^{+}}, \end{align}

(2.23b) \begin{align} \left .\frac{{\textrm{d}} \Phi }{{\textrm{d}} r}\right |_{r= a^{-}} &= \epsilon _r\left .\frac{{\textrm{d}} \Phi }{{\textrm{d}} r}\right |_{r = a^{+}}. \end{align}

2.4 Parameter estimation

We assume that the molecular volume is $\nu \sim 10^{-28}$ m $^{3}$ [34], the system is held at a temperature of $T = 300$ K and dry radius of the gel is $a_0\sim 1$ cm. Horkay et al. [14] measured the shear moduli of polyelectrolyte gels to be around $G \sim 10$ kPa, which leads to ${\mathcal{G}} \sim 10^{-4}$ . Yu et al. [34] reported values of ${\mathcal{G}} \sim 10^{-3}$ . The Flory interaction parameter $\chi$ is generally a function of the gel composition and temperature. However, we treat $\chi$ as a constant, which is a common simplification in the literature. Yu et al. [34] use constant values of $\chi$ that range from 0.1 to 1.6.

We assume that the electrical permittivity of the gel and the bath are approximately the same as water due to the ions being dilute. Thus, we set $\epsilon ^{\text{gel}} \simeq \epsilon ^{\text{bath}} \simeq 80\,\epsilon _0$ , where $\epsilon _0$ is the permittivity of free space. Hence, $\epsilon _r = \epsilon ^{\text{gel}}/ \epsilon ^{\text{bath}} \simeq 1$ . The non-dimensional width of the EDL is then $\varepsilon \sim 10^{-8}$ , corresponding to a dimensional value of $0.1$ nm. However, we will show in Sec. 5 that the value of $\varepsilon$ underestimates the width of the EDL because it is based on an imprecise estimate of the ionic volume fractions.

The dimensionless parameter $\omega$ is difficult to estimate due to uncertainties in the values of the Kuhn length. Hua et al. [15] set $L_K = 0.9$ nm in their modelling study. Similarly, Wu et al. [31] take $L_K = 1$ nm. Both values lead to an estimate of $\omega \sim 10^{-7}$ .

3. Asymptotic analysis for large Kuhn lengths

Matched asymptotic expansions in the limit $\varepsilon \to 0$ will now be used to formulate the governing equations away from and within the EDL at the gel–bath interface. The analysis in this section will focus on the case when the Kuhn length is much larger than the Debye length, $\varepsilon \ll \omega$ . Thus, we will consider the limit $\varepsilon \to 0$ with $\omega$ fixed. Although our estimates suggest that $\omega$ and $\varepsilon$ are similar in magnitude and hence the limit $\varepsilon \to 0$ with $\omega = O(\varepsilon )$ may be more physically accurate, we will show that the asymptotic solutions cannot generally be matched in this case. Analysing the case when $\varepsilon \to 0$ with $\omega$ fixed, i.e., $\varepsilon \ll \omega$ , provides mathematical and physical insights into why the matching fails.

The asymptotic analysis is split into two parts. In Sec. 3.1, we derive the model in the outer region away from the gel–bath interface. In Sec. 3.2, we formulate the problem in the inner region near the gel–bath interface to resolve the EDL.

3.1 The outer problem

We now consider the limit $\varepsilon \to 0$ with $r = O(1)$ . The outer variables are expanded as $f(r) = f^{(0)}(r) + O(\varepsilon )$ , where $f$ is an arbitrary quantity.

3.1.1 Electroneutral equations for the bath

Taking $\varepsilon \to 0$ in (2.13) leads to the electroneutrality condition for the bath

(3.1) \begin{align} z_+ \phi _+^{(0)} + z_- \phi _-^{(0)} = 0. \end{align}

The chemical potentials (2.12) reduce to

(3.2a) \begin{align} \mu ^{\text{bath}}_s &= \log \phi _s^{(0)}, \end{align}

(3.2b) \begin{align} \mu ^{\text{bath}}_{\pm } &= \log \phi _{\pm }^{(0)} + z_{\pm } \Phi ^{(0)}. \end{align}

Solving these four equations shows that the outer solution in the bath corresponds to a homogeneous mixture with the same composition and voltage as the far field:

(3.3) \begin{align} \phi ^{(0)}_{m}(r) = \phi ^{\text{bath}}_m, \quad \Phi ^{(0)}(r) = \Phi ^{\text{bath}}. \end{align}

The radial stress balance (2.14) reduces to ${\textrm{d}} p^{(0)}/{\textrm{d}} \xi = 0$ . Solving and matching to the far field implies that

(3.4) \begin{align} p^{(0)}(r) = 0, \quad \mathsf{T}_{rr}^{(0)}(r) = 0. \end{align}

Thus, the bath is stress free to leading order.

3.1.2 Electroneutral equations for the gel

Motivated by the constant outer solutions for the bath, as well as past works in the literature [9, 10, 14, 15, 23, 34, 35], we assume that the outer solution for the gel represents a homogeneous state. We therefore write the leading-order contributions to the outer solution as $f^{(0)}(r) = f^{\text{gel}}$ , where $f^{\text{gel}}$ is a constant, for all variables except the Lagrangian radius $R^{(0)}$ , which retains a dependence on $r$ .

The $O(1)$ contributions to (2.7) lead to the electroneutrality condition for the gel,

(3.5) \begin{align} z_+ \phi ^{\text{gel}}_+ + z_- \phi ^{\text{gel}}_- = - z_f \varphi _f/ J^{\text{gel}}. \end{align}

The solvent and ionic chemical potentials (2.5) are given by

(3.6a) \begin{align} \Pi ^{\text{gel}}_s +{\mathcal{G}} p^{\text{gel}} = \mu ^{\text{bath}}_s, \end{align}

(3.6b) \begin{align} \Pi ^{\text{gel}}_{\pm } +{\mathcal{G}} p^{\text{gel}} + z_{\pm } \Phi ^{\text{gel}} = \mu ^{\text{bath}}_{\pm }. \end{align}

For a homogeneous gel that is free to swell in the radial and orthoradial directions, the radial and orthoradial stretches are equal; thus, $\lambda ^{\text{gel}}_r = \lambda ^{\text{gel}}_{\theta } = (J^{\text{gel}}/\lambda _z)^{1/2}$ . The leading-order contribution to the Lagrangian coordinate can then be obtained from (2.4) as $R^{(0)}(r) = (\lambda _z/ J^{\text{gel}})^{1/2} r$ . We define

(3.7) \begin{align} R^{\text{gel}} \equiv R^{(0)}(a^{-}) = (\lambda _z/ J^{\text{gel}})^{1/2} a, \end{align}

where $R^{\text{gel}}$ must be determined by matching to the inner solution. The final quantity to determine is the mechanical pressure in the gel. The radial stress balance (2.8) implies that the radial Cauchy stress $\mathsf{T}_{rr}^{\text{gel}}$ must be a constant. Using asymptotic matching, we will show that $\mathsf{T}_{rr}^{\text{gel}} = 0$ . Hence, from (2.9a), we have that $p^{\text{gel}} = 1/ \lambda _z - 1/ J^{\text{gel}}$ .

3.2 The inner problem

The inner problem is formulated by introducing the change of variable $r = a + \varepsilon \xi$ , where $\xi$ is a radial coordinate that is localised to the free surface of the gel. By definition, $\xi \gt 0$ corresponds to the regions in the bath whereas $\xi \lt 0$ corresponds to regions in the gel. Tildes are used to denote dependent variables in the inner region, which are generally expanded as $\tilde{f} = \tilde{f}^{(0)} + \varepsilon \tilde{f}^{(1)} + O(\varepsilon ^2)$ , where $\tilde{f}$ is an arbitrary quantity. However, additional rescaling is required in some cases; this will be made explicit in the proceeding discussion. Near the interface, the outer solutions for the bath and gel are expanded in terms of inner variables, respectively, as

(3.8a) \begin{align} &f(a + \varepsilon \xi ) = f^{(0)}(a^{+}) + O(\varepsilon ) = f^{\text{bath}} + O(\varepsilon ), \end{align}

(3.8b) \begin{align} &f(a + \varepsilon \xi ) = f^{(0)}(a^{-}) + O(\varepsilon ) = f^{\text{gel}} + O(\varepsilon ), \end{align}

which will be used for asymptotic matching.

3.2.1 Inner problem for the bath

3.2.1.1 Mechanics

After changing variables, we anticipate that the Maxwell stresses (2.16) scale as $\mathsf{T}_M = O(\varepsilon ^{-2})$ because the electric potential $\Phi$ should remain $O(1)$ in size across the EDL. Moreover, we expect that the pressure will scale like the Maxwell stress so that $p = O(\varepsilon ^{-2})$ ; the rationale behind this choice is discussed in the introduction of Sec. 2. Therefore, we write $\mathsf{T}_M = \varepsilon ^{-2} \tilde{\mathsf{T}}_M$ and $p = \varepsilon ^{-2} \tilde{p}$ . Consequently, the Cauchy stresses must also be scaled as $\mathsf{T} = \varepsilon ^{-2} \tilde{\mathsf{T}}$ . Expanding $\tilde{\mathsf{T}}_{rr}$ and $\tilde{\mathsf{T}}_{M,rr}$ in powers of $\varepsilon$ and matching to the far field leads to the stress-free condition $\tilde{\mathsf{T}}_{rr}^{(0)}\to 0$ as $\xi \to \infty$ . The leading-order contribution to the radial stress balance in the bath (2.14) leads to ${\textrm{d}} \tilde{\mathsf{T}}^{(0)}_{rr}/{\textrm{d}} \xi = 0$ . Integrating and imposing the far-field condition leads to the conclusion that $\tilde{\mathsf{T}}_{rr}^{(0)} = 0$ . In terms of the original scaling, this means that, in the EDL, the total stress in the bath must be $O(\varepsilon ^{-1})$ in size. The pressure can be obtained from (2.15a) and is found to be

(3.9) \begin{align} \tilde{p}^{(0)} = \frac{1}{2}\left (\frac{\partial \tilde{\Phi }^{(0)}}{\partial \xi }\right )^2 = \tilde{\mathsf{T}}_{M,rr}^{(0)}. \end{align}

Thus, as expected, the pressure in the bath balances the Maxwell stresses.

3.2.1.2 Chemical equilibrium

Expanding the chemical potentials (2.12) gives, to leading order,

(3.10a) \begin{align} \log \tilde{\phi }_s^{(0)} + \epsilon _r \tilde{p}^{(0)} &= \mu ^{\text{bath}}_s, \end{align}

(3.10b) \begin{align} \log \tilde{\phi }_{\pm }^{(0)} + \epsilon _r \tilde{p}^{(0)} + z_{\pm } \tilde{\Phi }^{(0)} &= \mu ^{\text{bath}}_{\pm }. \end{align}

Equating (3.10) with (2.18a) and using (3.9) to eliminate the pressure leads to an expression for the ion fractions in the EDL,

(3.11) \begin{align} \tilde{\phi }_{\pm }^{(0)} = \phi ^{\text{bath}}_{\pm } \exp \left [z_{\pm }(\Phi ^{\text{bath}} - \tilde{\Phi }^{(0)}) - \frac{\epsilon _r}{2}\left (\frac{{\textrm{d}} \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi }\right )^2\right ]. \end{align}

This is a modification of the Boltzmann distribution for the ions, which arises from accounting for the pressure dependence of the ionic chemical potentials.

3.2.1.3 Electrostatics

The leading-order electrical problem is obtained by combining (2.13) with the ionic volume fractions (3.11) to obtain a modified Poisson–Boltzmann equation given by

(3.12) \begin{align} -\epsilon _r \frac{{\textrm{d}}^2 \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi ^2} = \exp \left [-\frac{\epsilon _r}{2} \left (\frac{{\textrm{d}} \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi }\right )^2\right ]\sum _{i\in \mathbb{I}}z_i \phi ^{\text{bath}}_i\exp \left (z_{i}(\Phi ^{\text{bath}} - \tilde{\Phi }^{(0)})\right ). \end{align}

Equation (3.12) can be integrated once and the conditions ${\textrm{d}} \tilde{\Phi }^{(0)}/{\textrm{d}} \xi \to 0$ and $\tilde{\Phi }^{(0)} \to \Phi ^{\text{bath}}$ as $\xi \to \infty$ used to obtain

(3.13) \begin{align} \frac{{\textrm{d}} \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi } = \mp \sqrt{\frac{2}{\epsilon _r}\log \left \{ 1 + \sum _{i\in \mathbb{I}}\phi ^{\text{bath}}_i\left [\exp \left (z_i(\Phi ^{\text{bath}} - \tilde{\Phi }^{(0)})\right ) -1 \right ]\right \}}. \end{align}

The minus sign is taken if $\Phi ^{\text{gel}} - \Phi ^{\text{bath}} \gt 0$ , which will generally be the case if the fixed charges on the polymer chains are positive, as assumed here.

3.2.2 Inner problem for the gel

3.2.1.4 Chemical equilibrium

The chemical potential of solvent (2.5a) can be expanded as

(3.14) \begin{align} \mu ^{\text{bath}}_s = \tilde{\Pi }_s^{(0)} +{\mathcal{G}} \tilde{p}^{(0)} - \varepsilon ^{-2}\omega ^2\Bigg [&\frac{{\textrm{d}}^2 }{{\textrm{d}} \xi ^2}\left ( \tilde{\phi }_s^{(0)} + \varepsilon \tilde{\phi }_s^{(1)} + \varepsilon ^2 \tilde{\phi }_s^{(2)}\right ) + \frac{\varepsilon }{a}\frac{{\textrm{d}} }{{\textrm{d}} \xi }\left (\tilde{\phi }_s^{(0)} + \varepsilon \tilde{\phi }_s^{(1)}\right ) - \varepsilon ^2 \frac{\xi }{a^2} \frac{{\textrm{d}} \tilde{\phi }_s^{(0)}}{{\textrm{d}} \xi }\Bigg ] + O(\varepsilon ). \end{align}

Similarly, the boundary condition at the gel–bath interface (2.21) can be expanded to give ${\textrm{d}} \tilde{\phi }_s^{(n)}/{\textrm{d}} \xi = 0$ at $\xi = 0$ for $n = 0, 1, 2$ . The $O(\varepsilon ^{-2})$ and $O(\varepsilon ^{-1})$ contributions to (3.14) along with the boundary and matching conditions show that the solvent concentration is uniform to leading and next order,

(3.15) \begin{align} \tilde{\phi }_s^{(0)}(\xi ) = \phi _s^{\text{gel}}, \quad \tilde{\phi }_s^{(1)}(\xi ) = \phi _s^{(1)}, \end{align}

which is a distinguishing feature of the asymptotic limit in which $\varepsilon \to 0$ with $\omega$ fixed. Physically, this result is a consequence of gradients in the solvent concentration having a high energy cost when the Kuhn length is large. The ionic chemical potentials (2.5b) can be expanded as

(3.16) \begin{align} \mu ^{\text{bath}}_{\pm } = \log \tilde{\phi }_{\pm }^{(0)} + \frac{1}{\tilde{J}^{(0)}}(1 - \chi \phi ^{\text{gel}}_s) +{\mathcal{G}} \tilde{p}^{(0)} + z_{\pm } \tilde{\Phi }^{(0)}. \end{align}

By combining (3.16) and (2.18a) and using (3.15), we find that

(3.17) \begin{align} \tilde{\phi }_{\pm }^{(0)} = \phi ^{\text{bath}}_{\pm }\exp \left [z_\pm (\Phi ^{\text{bath}}-\tilde{\Phi }^{(0)}) -{\mathcal{G}}\tilde{p}^{(0)} - \frac{1}{\tilde{J}^{(0)}}(1 - \chi \phi ^{\text{gel}}_s)\right ]. \end{align}

Although (3.17) appears to be a closed-form expression for the volume fraction of ions, it is important to recall that the swelling ratio $\tilde{J}^{(0)}$ also depends on the these quantities; see (2.3).

3.2.1.5 Kinematics and incompressibility

The $O(\varepsilon ^{-1})$ part of the incompressibility condition (2.4) implies that ${\textrm{d}} \tilde{R}^{(0)}/{\textrm{d}} \xi = 0$ . Imposing the boundary condition $\tilde{R}^{(0)}(0) = 1$ leads to $\tilde{R}^{(0)} = 1$ . Matching the inner and outer solutions leads to the condition $R^{\text{gel}} = 1$ , which can be used in combination with (3.7) to find that the radius of the deformed gel is given by

(3.18) \begin{align} a = (J^{\text{gel}}/ \lambda _z)^{1/2}. \end{align}

The $O(1)$ part of (2.4) gives ${\textrm{d}} \tilde{R}^{(1)}/{\textrm{d}} \xi = \lambda _z a/ \tilde{J}^{(0)}$ . The leading-order radial stretch can then be found by expanding (2.2) to find that $\tilde{\lambda }_r^{(0)} = (\lambda _z J^{\text{gel}})^{-1/2}\tilde{J}^{(0)}$ . Similarly, the leading-order orthoradial stretch is given by $\tilde{\lambda }_{\theta }^{(0)} = (J^{\text{gel}}/\lambda _z)^{1/2}$ .

3.2.1.6 Mechanics

The leading-order part of the stress balance in the gel (2.8) is ${\textrm{d}} \tilde{\mathsf{T}}^{(0)}_{rr}/{\textrm{d}} \xi = 0$ . Hence, the total stress in the gel is constant across the EDL. Imposing stress continuity at the interface (2.22) and using the fact that the stress in the bath is $O(\varepsilon ^{-1})$ shows that the gel is stress free to leading order, $\tilde{\mathsf{T}}_{rr}^{(0)} = 0$ . Matching to the outer solution leads to $\mathsf{T}_{rr}^{\text{gel}} = 0$ , as previously claimed.

The mechanical pressure can be obtained from (2.9a) as

(3.19) \begin{align} \tilde{p}^{(0)} = \tilde{\mathsf{T}}_{e,rr}^{(0)} + \tilde{\mathsf{T}}_{K,rr}^{(0)} + \tilde{\mathsf{T}}_{M,rr}^{(0)}. \end{align}

The leading-order radial components of the elastic, Korteweg, and Maxwell stresses can be obtained from (2.10a), (2.11a), and (2.11c) as

(3.20a) \begin{align} \tilde{\mathsf{T}}_{e,rr}^{(0)} &= \frac{1}{\lambda _z} \frac{\tilde{J}^{(0)}}{J^{\text{gel}}} - \frac{1}{\tilde{J}^{(0)}}, \end{align}

(3.20b) \begin{align} \tilde{\mathsf{T}}_{K,rr}^{(0)} &= \frac{\omega ^2}{{\mathcal{G}}} \phi ^{\text{gel}}_s\frac{{\textrm{d}}^2 \tilde{\phi }_s^{(2)}}{{\textrm{d}} \xi ^2}, \end{align}

(3.20c) \begin{align} \tilde{T}_{M,rr}^{(0)} &= \frac{1}{2{\mathcal{G}}}\left (\frac{{\textrm{d}} \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi }\right )^2. \end{align}

In order to evaluate the Korteweg stresses without explicitly solving for $\tilde{\phi }_s^{(2)}$ , the $O(1)$ contributions to the solvent chemical potential (3.14) can be used in (3.20b) to obtain

(3.21) \begin{align} \tilde{\mathsf{T}}_{K,rr}^{(0)} = \frac{\phi ^{\text{gel}}_s}{{\mathcal{G}}}\left (\tilde{\Pi }_s^{(0)} +{\mathcal{G}} \tilde{p}^{(0)} - \mu ^{\text{bath}}_s\right ). \end{align}

Note that setting $\omega = 0$ leads to $\tilde{\Pi }_s^{(0)} +{\mathcal{G}} \tilde{p}^{(0)} - \mu ^{\text{bath}}_s = 0$ from (3.14) and hence $\tilde{\mathsf{T}}_{K,rr}^{(0)} = 0$ . Substitution of (3.21) into (3.19) gives an algebraic relation for the pressure $\tilde{p}^{(0)}$ .

3.2.1.7 Electrostatics

The leading-order electrical problem in the gel is given by

(3.22) \begin{align} -\frac{{\textrm{d}}^2 \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi ^2} = z_+ \tilde{\phi }_{+}^{(0)} + z_-\tilde{\phi }_-^{(0)} + \frac{z_f\varphi _f}{\tilde{J}^{(0)}}, \end{align}

which is coupled to the algebraic equations for the volume fractions of ions (3.17) and the mechanical pressure (3.19). The electrical problems for the bath and gel can be decoupled by combining the first integral for the electric potential in the bath (3.13) with the electrostatic boundary conditions (2.23) to obtain

(3.23a) \begin{align} \left .\frac{{\textrm{d}} \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi }\right |_{\xi = 0^{-}} = \left .\mp \sqrt{2\epsilon _r \log \left \{ 1 + \sum _{i\in \mathbb{I}}\phi ^{\text{bath}}_i\left [\exp \left (z_i(\Phi ^{\text{bath}} - \tilde{\Phi }^{(0)})\right ) -1 \right ]\right \}}\right |_{\xi = 0^{-}}, \end{align}

which acts as a boundary condition for (3.22). The electrical problem in the gel is closed by imposing the matching condition

(3.23b) \begin{align} \tilde{\Phi }^{(0)} \to \Phi ^{\text{gel}}, \quad \xi \to -\infty. \end{align}

3.3 Summary

The asymptotic analysis has produced a closed system of algebraic equations that determines the outer solution in the gel. The inner problems for the gel and the bath decouple. In the case of the gel, the inner problem can be condensed into a nonlinear system of differential-algebraic equations; this system will be presented in Sec. 5. The inner problem for the bath has been solved in terms of the electric potential; this can be obtained by integrating (3.13) once the electric potential at the gel surface has been determined by solving the inner problem for the gel.

4. Asymptotic analysis for Kuhn lengths of zero

We now briefly mention how the asymptotic analysis of the inner region is carried out when the dimensionless Kuhn length, $\omega$ , is naively set to zero. Mathematical models of polyelectrolyte gels in the zero-Kuhn-length limit ( $\omega = 0$ ) are common throughout the literature. Analysing the inner region when $\omega = 0$ will serve as a useful point of comparison. However, the zero-Kuhn-length limit is only valid when phase separation does not occur. The terms that are proportional to $\omega ^2$ in the governing equations provide a regularisation that ensures the problem remains well posed when the system undergoes phase separation. As we will show in Sec. 5, the EDL can trigger phase separation which then spreads into the bulk of the gel. Hence, setting $\omega = 0$ is not trivial and may render the model ill posed.

Assuming that phase separation does not occur, we can set $\omega = 0$ in (3.14), in which case the leading-order contribution (in $\varepsilon$ ) to the solvent chemical potential in the EDL becomes

(4.1a) \begin{align} \tilde{\Pi }_s^{(0)} +{\mathcal{G}} \tilde{p}^{(0)} = \mu ^{\text{bath}}_s, \end{align}

The osmotic pressure $\tilde{\Pi }_s^{(0)}$ is given by (2.6a). The mechanical pressure $\tilde{p}^{(0)}$ can be calculated directly from (3.19) after neglecting the Korteweg stresses. Thus, (4.1a) can be interpreted as a nonlinear algebraic equation for the solvent fraction $\tilde{\phi }^{(0)}_s$ , which can now be a function of $\xi$ and hence vary across the EDL. The corresponding ion fractions are given by

(4.1b) \begin{align} \tilde{\phi }_{\pm }^{(0)} = \phi ^{\text{bath}}_{\pm }\exp \left [z_\pm (\Phi ^{\text{bath}}-\tilde{\Phi }^{(0)}) -{\mathcal{G}}\tilde{p}^{(0)} - \frac{1}{\tilde{J}^{(0)}}\left (1 - \chi \tilde{\phi }^{(0)}_s\right )\right ]. \end{align}

5. Case studies

We now use our formulation to study the structure of the EDL by computing numerical solutions to the inner and outer problems. We restrict our attention to monovalent salts with $z_{\pm } = \pm 1$ . The cation fraction in the bath, $\phi ^{\text{bath}}_+$ , is treated as a control parameter.

In Sec. 5.1, the outer problem in the gel is formulated. This consists of a system of nonlinear algebraic equations for homogeneously swollen states that are in equilibrium with the bath. In Sec. 5.2, the corresponding inner problems are formulated when $\omega = 0$ and $\omega \gg \varepsilon$ . The inner solutions are used to explore the structure of the EDL in Sec. 5.3.

5.1 Solution of the outer problem for the gel

In the outer region of the gel, the volume fraction of solvent and ions, as well as the electric potential, are determined from (3.5) to (3.6). This nonlinear algebraic system can be formulated as

(5.1a) \begin{align} \log \phi ^{\text{gel}}_s + \frac{1}{J^{\text{gel}}} + \frac{\chi (1 - \phi ^{\text{gel}}_s)}{J^{\text{gel}}} +{\mathcal{G}}\left (\frac{1}{\lambda _z} - \frac{1}{J^{\text{gel}}}\right ) = \log (1 - 2 \phi ^{\text{bath}}_+), \end{align}

(5.1b) \begin{align} \phi ^{\text{gel}}_{\pm } = \phi ^{\text{bath}}_+ \exp \left [\pm (\Phi ^{\text{bath}} - \Phi ^{\text{gel}}) -{\mathcal{G}}\left (\frac{1}{\lambda _z}-\frac{1}{J^{\text{gel}}}\right ) - \frac{1}{J^{\text{gel}}}\left (1 - \chi \phi ^{\text{gel}}_s\right )\right ], \end{align}

(5.1c) \begin{align} 2 \phi ^{\text{bath}}_+ \sinh (\Phi ^{\text{bath}} - \Phi ^{\text{gel}}) = -\frac{z_f \varphi _f}{J^{\text{gel}}}\exp \left [{\mathcal{G}}\left (\frac{1}{\lambda _z}-\frac{1}{J^{\text{gel}}}\right ) + \frac{1}{J^{\text{gel}}}\left (1 - \chi \phi ^{\text{gel}}_s\right )\right ], \end{align}

where $J^{\text{gel}}$ is given by (2.3). When the cation fraction in the bath is small, $\phi ^{\text{bath}}_{+} \ll 1$ , the nonlinear system (5.1) can be reduced to a single equation, as described in Appendix B.

We numerically solve (5.1) over a range of values of $\phi ^{\text{bath}}_{+}$ using pseudo-arclength continuation. Three values of $\lambda _z \leq 1$ are considered, corresponding to gels in axial compression. The equilibrium swelling ratios $J^{\text{gel}}$ computed from the solutions of (5.1) are shown as solid curves in Figure 2. The dashed black lines represent numerical solutions to the reduced model derived in Appendix B. The figure shows that for each value of $\lambda _z$ there are three branches of solutions, two of which intersect and then vanish as the salt fraction in the bath $\phi ^{\text{bath}}_{+}$ increases beyond a critical value. The loss of equilibrium solutions at this critical point indicates that a volume phase transition can occur, as the gel volume will undergo a discontinuous decrease as the salt content of the bath is increased. For a given value of $\lambda _z$ , the branch of solutions corresponding to the largest and smallest values of $J^{\text{gel}}$ are referred to as the swollen and collapsed branches, respectively. The branch of solutions corresponding to intermediate values of $J^{\text{gel}}$ is unstable [4] and will not be considered further.

Figure 2. (a) Equilibrium swelling ratio $J^{\text{gel}}$ as a function of cation fraction in the bath $\phi ^{\text{bath}}_+$ showing swollen and collapsed branches. (b) The swelling ratio along the collapsed branch. Solid lines correspond to solutions of (5.1). Dashed lines represent solutions to the reduced equation (B.4) for a dilute concentration of cations. The parameter values are ${\mathcal{G}} = 5 \times 10^{-4}$ , $\chi = 1.2$ , $\varphi _f = 0.05$ $z_\pm = \pm 1$ , $z_f = 1$ .

For a fixed value of the salt fraction, increasing the axial compression reduces the degree of swelling that occurs. Moreover, increasing the axial compression also decreases the critical salt fraction at which the volume phase transition occurs. Both of these findings are in agreement with experimental observations [14]. Due to the incompressibility of the gel, imposing an axial compression results in a radial stretch. The elastic energy cost of inserting a molecule into a pre-stretched gel is greater than for an unstretched gel. Hence, the balance between the mixing and elastic energies is established at smaller concentrations, resulting in the equilibrium swelling ratio $J^{\text{gel}}$ decreasing with the axial stretch $\lambda _z$ .

5.2 Formulation of the inner problems

The inner problem for the gel can now be constructed using the results from the previous sections. In particular, if $\omega = 0$ , then the governing equations for the gel can be condensed into

(5.2a) \begin{align} \tilde{\Pi }_s^{(0)}+{\mathcal{G}}\tilde{p}^{(0)} = \log (1 - 2 \phi ^{\text{bath}}_+), \end{align}

(5.2b) \begin{align} \tilde{\phi }^{(0)}_{\pm } = \phi ^{\text{bath}}_+ \exp \left [\pm (\Phi ^{\text{bath}} - \tilde{\Phi }^{(0)}) -{\mathcal{G}}\tilde{p}^{(0)} - \frac{1}{\tilde{J}^{(0)}}\left (1 - \chi \tilde{\phi }^{(0)}_s\right )\right ], \end{align}

(5.2c) \begin{align} -\frac{{\textrm{d}} \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi } = \tilde{\phi }_{+}^{(0)} - \tilde{\phi }_-^{(0)} + \frac{z_f \varphi _f}{\tilde{J}^{(0)}}, \end{align}

(5.2d) \begin{align} \tilde{p}^{(0)} = \frac{1}{\lambda _z} \frac{\tilde{J}^{(0)}}{J^{\text{gel}}} - \frac{1}{\tilde{J}^{(0)}} + \frac{1}{2{\mathcal{G}}}\left (\frac{{\textrm{d}} \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi }\right )^2, \end{align}

(5.2e) \begin{align} \tilde{J}^{(0)} = (1 - \tilde{\phi }^{(0)}_s - \tilde{\phi }_{+}^{(0)} - \tilde{\phi }_{-}^{(0)})^{-1}, \end{align}

(5.2f) \begin{align} \tilde{\Pi }_s^{(0)} = \log \tilde{\phi }^{(0)}_s + \frac{\chi (1 - \tilde{\phi }_s^{(0)})}{\tilde{J}^{(0)}} + \frac{1}{\tilde{J}^{(0)}}. \end{align}

In the case $\omega \gg \varepsilon$ , equation (5.2a) is replaced with $\tilde{\phi }_s^{(0)} = \phi ^{\text{gel}}_s$ , resulting in the system

(5.3a) \begin{align} \tilde{\phi }^{(0)}_{\pm } = \phi ^{\text{bath}}_+ \exp \left [\pm (\Phi ^{\text{bath}} - \tilde{\Phi }^{(0)}) -{\mathcal{G}}\tilde{p}^{(0)} - \frac{1}{\tilde{J}^{(0)}}\left (1 - \chi \phi ^{\text{gel}}_s\right )\right ], \end{align}

(5.3b) \begin{align} -\frac{{\textrm{d}}^2 \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi ^2} = \tilde{\phi }_{+}^{(0)} - \tilde{\phi }_-^{(0)} + \frac{z_f \varphi _f}{\tilde{J}^{(0)}}, \end{align}

(5.3c) \begin{align}{\mathcal{G}} (1 - \phi ^{\text{gel}}_s) \tilde{p}^{(0)} ={\mathcal{G}}\left (\frac{1}{\lambda _z} \frac{\tilde{J}^{(0)}}{J^{\text{gel}}} - \frac{1}{\tilde{J}^{(0)}}\right ) + \phi ^{\text{gel}}_s \left (\tilde{\Pi }_s^{(0)} - \mu ^{\text{bath}}_s\right ) + \frac{1}{2}\left (\frac{{\textrm{d}} \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi }\right )^2, \end{align}

(5.3d) \begin{align} \tilde{J}^{(0)} = (1 - \phi ^{\text{gel}}_s - \tilde{\phi }_{+}^{(0)} - \tilde{\phi }_{-}^{(0)})^{-1}, \end{align}

(5.3e) \begin{align} \tilde{\Pi }_s^{(0)} = \log \phi ^{\text{gel}}_s + \frac{\chi (1 - \phi ^{\text{gel}}_s)}{\tilde{J}^{(0)}} + \frac{1}{\tilde{J}^{(0)}}. \end{align}

In both cases, the boundary conditions for the electrical potential are given by (3.23a). The expression for the hoop stress in the gel is the same in both cases as well:

(5.4) \begin{align} \tilde{\mathsf{T}}^{(0)}_{\theta \theta } = \frac{1}{\lambda _z}\left (\frac{J^{\text{gel}}}{\tilde{J}^{(0)}} - \frac{\tilde{J}^{(0)}}{J^{\text{gel}}}\right ) -{\mathcal{G}}^{-1}\left (\frac{{\textrm{d}} \tilde{\Phi }^{(0)}}{{\textrm{d}} \xi }\right )^2. \end{align}

The first term represents the elastic contribution to the total hoop stress, which can be compressive or tensile. The second term captures the contribution from the Maxwell stresses, which is always compressive.

Figure 3. Numerical solutions of the inner problems with far-field conditions corresponding to the collapsed state. The regions of gel and bath are defined by $\xi \lt 0$ and $\xi \gt 0$ , respectively. Dashed lines represent the solution to (5.2) when $\omega = 0$ . Solid lines represent the solution to (5.3) when $\omega \gg \varepsilon$ . The parameter values are $\chi = 1.2$ , ${\mathcal{G}} = 5 \times 10^{-4}$ , $\varphi _f = 0.05$ , $\phi ^{\text{bath}}_+ = 10^{-5}$ , $\lambda _z = 1$ , $\epsilon _r = 1$ , $z_{\pm } = \pm 1$ , and $z_f = 1$ .

Figure 4. Numerical solution of the inner problem with far-field conditions corresponding to the swollen state when $\omega \gg \varepsilon$ . The regions of gel and bath are defined by $\xi \lt 0$ and $\xi \gt 0$ , respectively. The curves are obtained by solving (5.2). Parameter values are the same as in Figure 3: $\chi = 1.2$ , ${\mathcal{G}} = 5 \times 10^{-4}$ , $\varphi _f = 0.05$ , $\phi ^{\text{bath}}_+ = 10^{-5}$ , $\lambda _z = 1$ , $\epsilon _r = 1$ , $z_{\pm } = \pm 1$ , and $z_f = 1$ .

Figure 5. Phase separation in the inner region. (a)–(c) The swelling ratio and (d)–(f) total electric charge density for different salt fractions in the bath $\phi ^{\text{bath}}_{+}$ . The curves are obtained by numerically solving the inner problem in the intermediate asymptotic limit $\omega = O(\varepsilon )$ as $\varepsilon \to 0$ ; see (5.5). We have taken $\omega = \Omega \varepsilon$ with $\Omega = 10^{-1}$ . The remaining parameters are $\chi = 0.7$ , ${\mathcal{G}} = 4\times 10^{-3}$ , $\varphi _f = 0.04$ , $z_{\pm } = \pm 1$ , $z_f = 1$ , $\epsilon _r = 1$ , and $\lambda _z = 1$ .

5.3 The structure of the electric double layer

The systems (5.2) and (5.3) are discretised using finite differences and solved using Newton’s method. The asymptotic solution is validated against numerical solutions of the full problem in Appendix C.

We consider the case where the axial stretch and salt content in the bath are set to $\lambda _z = 1$ and $\phi ^{\text{bath}}_{+} = 10^{-5}$ , with the remaining parameters being the same as those in Figure 2. There are three possible solutions to the outer problem. We are only concerned with two of these, which correspond to the collapsed state ( $J^{\text{gel}} \simeq 1.447$ ) and the swollen state ( $J^{\text{gel}} \simeq 82$ ).

In Figure 3, we plot the inner solution when the outer solution corresponds to the collapsed state ( $J^{\text{gel}} \simeq 1.447$ ). The solid and dashed lines correspond to the cases $\omega \gg \varepsilon$ and $\omega = 0$ , respectively. In both cases, we see that our non-dimensionalisation underestimates the width of the double layer, which is about $10 \varepsilon$ in the gel (or 1 nm) and $1000 \varepsilon$ in the bath (or 100 nm). For this parameter set, the value of $\omega$ does not lead to noticeable changes in the electric potential and ion fractions; see Figure 3 (a)–(b). However, substantial differences arise in the gel pressure and the solvent fraction; see Figure 3 (c)–(d). When $\omega = 0$ , the gel pressure balances a large Maxwell stress. This large pressure causes a local decrease in the solvent fraction and a slight volumetric contraction of the gel (Figure 3 (d)), which can be rationalised in terms of equation (4.1a). At equilibrium, the osmotic pressure $\tilde{\Pi }_s$ must balance the mechanical pressure $\tilde{p}$ . To compensate for the increase in mechanical pressure that arises from the Maxwell stress, the osmotic pressure must decrease, which drives solvent out of the gel and causes it to shrink. When $\omega \gg \varepsilon$ , gradients in the solvent fraction are energetically penalised; thus, the solvent fraction remains uniform across the EDL. From a mechanical perspective, this penalisation occurs through the development of a large Korteweg stress, which counters the effect of the Maxwell stress in order to maintain a uniform solvent fraction. The mechanical contribution from the Korteweg stress manifests as an increase in the gel pressure compared to the $\omega = 0$ case, as seen in Figure 3 (c). Although the solvent fraction is constant across the EDL when $\omega \gg \varepsilon$ , the swelling ratio still decreases relative to the bulk value (Figure 3 (d)) due to the variation in ionic content (Figure 3 (b)).

The inset of Figure 3 (c) shows the total hoop stress in the gel, which is the same in both models owing to the strong similarities in the electric potential. Due to the large Maxwell stress, the gel experiences a substantial compressive hoop stress, which leads to the intriguing possibility of mechanical instabilities in the EDL.

In Figure 4, we show the numerical solution of the inner problem with $\omega \gg \varepsilon$ when the outer solution corresponds to the swollen state ( $J^{\text{gel}} \simeq 82$ ). The qualitative features of the inner solution are similar to those obtained when the outer solution corresponds to the collapsed state (Figure 3). However, an important difference is that the volume fraction of anions has decreased by more than a factor of ten. This decrease is driven by the reduction in the volume fraction of fixed charges on the polymers that occurs when the gel is highly swollen; see Figure 4 (b). Consequently, the EDL in the gel has increased in thickness by a roughly factor of ten to approximately $250 \varepsilon$ (or 25 nm). The gradient in the electric potential in the gel is therefore ten times weaker, resulting in a 100-fold reduction in the Maxwell stress and the total hoop stress; see Figure 4 (c). Despite these decreases, the pressure in the gel remains large because of the Korteweg stress. Due to convergence issues, it was not possible to compute the corresponding inner solution when $\omega = 0$ .

To understand the origin of these numerical difficulties, we consider an intermediate asymptotic limit where $\omega = \Omega \varepsilon$ , with $\Omega = O(1)$ as $\varepsilon \to 0$ . By assuming that the gel remains in a homogeneous, swollen state away from the EDL, the outer problem in this limit is still given by (5.1). The corresponding inner problem can be formulated by changing (3.14) or (4.1a) to

(5.5a) \begin{align} \tilde{\Pi }_s^{(0)} +{\mathcal{G}} \tilde{p}^{(0)} + \Omega ^2 \frac{{\textrm{d}}^2 \tilde{\phi }_s^{(0)}}{{\textrm{d}} \xi ^2} = \mu ^{\text{bath}}_s. \end{align}

The pressure (3.19) can be evaluated using a Korteweg stress given by

(5.5b) \begin{align} \tilde{\mathsf{T}}_{K,rr}^{(0)} ={\mathcal{G}}^{-1} \Omega ^2 \left [\tilde{\phi }_s^{(0)}\frac{\textrm{d} ^2 \tilde{\phi }_s^{(0)}}{\textrm{d} \xi ^2} - \frac{1}{2}\left (\frac{\textrm{d} \tilde{\phi }_s^{(0)}}{\textrm{d} \xi }\right )^2\right ]. \end{align}

We solve the intermediate asymptotic model based on (5.5) by imposing the boundary conditions ${\textrm{d}} \tilde{\phi }_s^{(0)}/{\textrm{d}} \xi = 0$ as $\xi \to -\infty$ and $\xi \to 0^{-}$ . A second parameter set is used to reduce the degree of swelling that occurs in the gel. This parameter set leads to the outer problem having single branch of equilibrium solutions that does not fold back on itself. Thus, the gel monotonically and continuously decreases in volume as the salt fraction in the bath $\phi ^{\text{bath}}_+$ increases.

Figure 6. Phase separation drives the breakdown of charge neutrality in the gel when the Debye length is comparable to the Kuhn length. (a) The swelling ratio and (b) the total electric charge density computed from the full steady problem in cylindrical coordinates. The gel self-organises into a highly swollen, negatively charged core (dark blue); a moderately swollen interior with alternating electric charge (blue), and a weakly swollen, positively charged shell (light blue). We set $\varepsilon = 10^{-2}$ , $\omega = 10^{-3}$ . The remaining parameter values correspond to Figure 5 (b) and (e) and are $\phi ^{\text{bath}}_+ = 6.6 \times 10^{-4}$ , $\chi = 0.7$ , ${\mathcal{G}} = 4\times 10^{-3}$ , $\varphi _f = 0.04$ , $z_{\pm } = \pm 1$ , $z_f = 1$ , $\epsilon _r = 1$ , and $\lambda _z = 1$ .

It is important to point out that the intermediate asymptotic model based on (5.5) is only fourth order in space. Therefore, it can be solved without explicitly imposing that its solutions tend to the homogeneous and electrically neutral outer solutions determined by (5.1). However, if the inner solutions tend to constants in the far field, then these constants must satisfy (5.1) and hence a match with the outer solution will be obtained.

The inner problem in the intermediate asymptotic limit is solved at three specific values of $\phi ^{\text{bath}}_{+}$ with $\Omega = 0.1$ . The swelling ratio $\tilde{J}^{(0)}$ and total charge density $\tilde{Q}^{(0)} = \tilde{\phi }^{(0)}_{+} - \tilde{\phi }_{-}^{(0)} + z_f \varphi _f/ \tilde{J}^{(0)}$ are computed and plotted in Figure 5. In this case, decreasing the salt fraction in the bath from $\phi ^{\text{bath}}_+ = 10^{-3}$ triggers the onset of phase separation, which gives rise to a periodic array of electrically charged phases that spans the entire inner region. Charge neutrality is not recovered in the far field, even if the domain used to numerically solve the inner problem is increased. Moreover, enforcing the boundary condition $\tilde{\phi }_s^{(0)} \to \phi ^{\text{gel}}_s$ as $\xi \to -\infty$ leads to convergence issues. Hence, the inner solution cannot be matched to the homogeneous and electroneutral outer solutions computed from (5.1). We therefore posit that homogeneous outer solutions do not always exist in the thin-EDL limit $\varepsilon \to 0$ if $\omega = O(\varepsilon )$ or $\omega = 0$ .

To explore the hypothesis that the bulk of the gel may not be homogeneous and electrically neutral at equilibrium, we solved the full steady problem in the regime $\omega = O(\varepsilon )$ by setting $\varepsilon = 10^{-2}$ and $\omega = 10^{-3}$ . The salt fraction in the bath was set to $\phi ^{\text{bath}}_+ = 6.6\times 10^{-4}$ , corresponding to the parameters in Figure 5 (b) and (e). The swelling ratio $J$ and the total charge density $Q$ , which are shown in Figure 6, reveal that phase separation occurs throughout the entire gel and gives rise to a periodic arrangement of electrically charged domains. Using numerical integration, we find that the total amount of electric charge contained within a pair of adjacent domains is on the order of $10^{-7}$ . Thus, the gel effectively separates into three distinct regions consisting of an electrically negative, highly swollen core ( $0 \lt r \lt 0.073$ ); a moderately swollen interior that is electrically neutral on average ( $0.073 \lt r \lt 2.0$ ); and a positively charged, collapsed shell ( $2.0 \lt r \lt 2.1$ ). Overall, the gel carries a net positive charge which exactly balances the net negative charge in the bath to ensure that charge neutrality holds on a global scale. The pointwise breakdown of charge neutrality across the gel indicates that it is not always appropriate to decompose the problem into inner and outer regions that are characterised by the local charge density of the gel.

6. Discussion and conclusion

Asymptotic and numerical methods are used to study the EDL that forms at the interface between a polyelectrolyte gel and a salt bath. The gel is described using a phase-field model, which introduces an additional length scale, the Kuhn length, into the problem. The Kuhn length measures the thickness of diffuse internal interfaces that can form due to phase separation within the gel. The ratio of the non-dimensional Kuhn and Debye lengths, $\omega$ and $\varepsilon$ , has a profound influence on the structure of equilibrium solutions.

When $\omega \gg \varepsilon$ , there is a high energy cost associated with forming gradients in the solvent volume fraction in the gel. Therefore, the solvent volume fraction is spatially uniform across the EDL and phase separation is suppressed. In this case, it is always possible to match the inner solutions to electrically neutral, homogeneous outer solutions. In contrast, when $\omega = 0$ , there is no energy penalty associated with forming gradients in the solvent fraction. In this case, the solvent fraction, which is set by a balance between the osmotic and mechanical pressures, can vary across the EDL. However, it is not always possible to compute a numerical solution to the inner problem when $\omega = 0$ .

Our preliminary investigation of the intermediate asymptotic limit where $\varepsilon \to 0$ with $\omega = O(\varepsilon )$ reveals that phase separation can result in heterogeneous gels consisting of repeating pairs of positively and negatively charged domains. The breakdown of charge neutrality means that the inner region effectively spans the entire gel. The difficulties in numerically computing inner solutions with $\omega = 0$ are attributed to the gel undergoing phase separation and the loss of homogeneous, electrically neutral outer solutions.

The breakdown of electroneutrality due to phase separation can be rationalised as follows. Phase separation leads to the formation of diffuse interfaces that separate domains with distinct compositions and electric potentials. The gradient in the electric potential across the diffuse interface generates an electric field. When the Kuhn and Debye lengths are commensurate, the electric field near the diffuse interface will be of sufficient magnitude to trigger the formation of an EDL within the gel. If the Kuhn length greatly exceeds the Debye length, then the electric field is too weak to generate an internal EDL and hence the gel remains electrically neutral.

In Celora et al. [5], we used numerical continuation to track solutions of the full steady problem as the salt fraction in the bath is varied in the regime when $\omega$ and $\varepsilon$ are comparable. We found that the breakdown of charge neutrality in the gel occurs via a cascade of saddle-node bifurcations associated with a spatially localised mode of phase separation that originates from the EDL. A more in-depth analysis of the asymptotic limit $\varepsilon \to 0$ with $\omega = O(\varepsilon )$ can shed light on how phase separation is triggered near the free surface of the gel and spreads into the bulk. Setting $\varepsilon \to 0$ with $\omega = O(\varepsilon )$ is expected to be mathematically interesting as it requires relaxing the assumption that the outer solutions are homogeneous and it involves taking the limit in which the thickness of the EDL and the thickness of diffuse internal interfaces simultaneously tend to zero.

Models of polyelectrolyte gels usually do not account for phase separation and thus implicitly set $\omega = 0$ . Homogeneous and hence electrically neutral solutions that neglect the EDL are often sought and compared against experimental data. However, our results show that these homogeneous ‘solutions’ may be asymptotically inconsistent because there is no inner solution that can be matched to them. Importantly, when $\omega = 0$ , the bulk behaviour of the gel can be strongly coupled to the behaviour in the EDL and thus the latter must be explicitly considered when constructing model solutions. The extensive use of homogeneous, electroneutral solutions to characterise the response of highly swollen polyelectrolyte gels is more consistent with the assumption that $\omega \gg \varepsilon$ , as this limit enables the successful matching of inner and outer solutions and prohibits the breakdown of electroneutrality in the bulk of the gel.

A useful extension of this work is to carry out the asymptotic analysis for a general geometry to produce an asymptotically consistent electroneutral model in three dimensions. Such a model would be a useful tool for designing polyelectrolyte gels that undergo programmable shape changes driven by mechanical instabilities [28], phase transitions [4] or imposed electric fields [25].