1. Introduction
Polymeric networks hosting fixed and mobile ionic charges have important technological functions, e.g. as energy convertors, sensors, signal processors and actuators (Tepermeister et al. Reference Tepermeister, Bosnjak, Dai, Zhang, Kielar, Wang, Tian, Suntivich and Silberstein2022). The electrical conductivity and capacitance of these ‘soft ionics’ are often reported from impedance spectroscopy as frequency spectra or their Nyquist counterparts. These present distinctive characteristics that reflect, for example, charge concentration and mobility. Despite a long history in many scientific fields, analysis of impedance spectroscopy data often rests on empirical equivalent-circuit modelling approaches, which are relatively simple to apply, but not necessarily straightforward to interpret (von Hauff Reference von Hauff2019). Physically grounded models have highlighted complex macroscale dynamics as emerging from materials with simple, uniform microstructures, e.g. electrolytes with frequency independent conductivity and dielectric permittivity (Macdonald Reference Macdonald1953; Ho, Raistrick & Huggins Reference Ho, Raistrick and Huggins1980). None have considered the frequency dependence that fluid dynamics, viscoelasticity and microstructural heterogeneity bring to such measurements. Accordingly, the present study targets the microscale heterogeneity of soft ionics, and how this promotes the coupled fluid, elastic and charge dynamics that shape their intrinsic frequency dependent conductivity and dielectric permittivity. This is similar to the manner in which the intrinsic dielectric response of colloidal dispersions, which are soft, heterogeneous ionic conductors, has been elucidated on the basis of what is known as the ‘standard electrokinetic model’, namely electric-flux conservation, ion conservation with Nernst–Planck ion fluxes and low-Reynolds-number hydrodynamics (Delacey & White Reference Delacey and White1981).
Chang & Jaffé (Reference Chang and Jaffé1952) are attributed with the first physically based model furnishing the complex-valued current in an ion-conducting sample subjected to a small-amplitude oscillatory electric field, soon after generalised by Macdonald (Reference Macdonald1953). These models consider a uniform continuum in which dynamics, in the absence of microhydrodynamic and rheological physics, arise from electromigration and diffusion of charge carriers with a frequency-independent mobility. As illustrated succinctly by Coelho (Reference Coelho1991), the dynamics in a parallel electrode configuration (electrode separation 
$L$) furnish a frequency response that depends on 
$\kappa L$, which is the electrode spacing 
$L$ scaled with the Debye length 
$\kappa ^{-1}$, and on 
$\omega \tau$, which is the angular frequency 
$\omega$ scaled with a diffusion time 
$\tau = \kappa ^{-2} / D$ (
$D$ is the charge-carrier diffusion coefficient). At low frequencies, charge accumulation at the electrodes polarises the sample, manifesting as a decrease in the apparent conductivity and an increase in the apparent dielectric constant. Coelho (Reference Coelho1991) shows that the response in this frequency range is a Debye relaxation (von Hauff Reference von Hauff2019).
In the context of measuring the impedance of electrolyte solutions, Hollingsworth & Saville (Reference Hollingsworth and Saville2003) derived a formula, based of the standard electrokinetic model, that can be demonstrated to be equivalent to that of Coelho (Reference Coelho1991). The analogy highlights the standard electrokinetic model as being a generalisation, to accommodate fluid advection, of the charge-transport model of Chang & Jaffé (Reference Chang and Jaffé1952). However, in the one-dimensional context of a uniform electrolyte bounded by parallel, rigid, impermeable, electrodes, incompressible fluid dynamics were completely arrested. This simplification is relaxed, in part, by the present work, which accounts for heterogeneity of a compressible polyelectrolyte network. This permits a uniform translation of an incompressible electrolyte relative to planar deformation (compression and rarefaction) of a permeable, viscoelastic, charged polymer skeleton.
Hollingsworth & Saville (Reference Hollingsworth and Saville2003) highlight that intrinsic material properties from impedance spectroscopy can be contaminated or completely obscured by electrode polarisation. For example, the apparent dielectric constant (low-frequency limit of the real part of the apparent dielectric permittivity) of a simple electrolyte approaches a factor that is 
$\kappa L \gg 1$ relative to its intrinsic dielectric constant. The model also predicts an apparent conductivity that vanishes at low frequency and, perhaps paradoxically, tends to a high-frequency plateau that equals the intrinsic d.c. electrolyte conductivity. Hollingsworth & Saville (Reference Hollingsworth and Saville2004) demonstrate how the electrode-polarisation model may be applied to extract from impedance spectra the intrinsic polarisability spectrum of dispersed colloidal microspheres: a soft, heterogeneous ionic conductor in which fluid dynamics play an important role (Delacey & White Reference Delacey and White1981). These ideas have been applied by Klein et al. (Reference Klein, Zhang, Dou, Jones, Colby and Runt2006) to interpret the dielectric spectra of single-ion polymer electrolytes. Despite such polymers having mobile ions at high (
${\sim }1$ M) concentrations, the absence of a solvent furnishes low conductivities (
${\lesssim }0.1\ {\rm mS}\ {\rm cm}^{-1}$), which Klein et al. (Reference Klein, Zhang, Dou, Jones, Colby and Runt2006) attribute to low ion mobilities (diffusion coefficients 
$\sim 0.1$–
$1 \times 10^{-11}\ {\rm m}^2{\rm s}^{-1}$ at 
$20\,^\circ$C) and a small fraction of these ions being mobile (
$\sim 0.1$–
$1 \times 10^{-4}$ at 
$20\,^\circ$C).
Complex microscale dynamics are evident from impedance spectroscopy applied to Nafion membranes (Kusoglu & Weber Reference Kusoglu and Weber2017). Despite decades of study (Mauritz & Fu Reference Mauritz and Fu1988), many features of their dielectric permittivity spectra are still not well understood (Matos Reference Matos2020). The polymers have hydrophobic (perfluorovinyl ether) groups and ionic (sulfonate) side groups, which together produce nanoscale heterogeneity by phase separation when hydrated. High-resolution imaging reveals an interconnected channel network with domain spacing 
${\approx }5$ nm (Allen et al. Reference Allen, Comolli, Kusoglu, Modestino, Minor and Weber2015). Such heterogeneity is also a prominent feature of the conducting polymer hydrogels of Pan et al. (Reference Pan2012), which have pore sizes from the angstrom to micrometre scales. Heterogeneity in the fixed charge density, which is achieved while maintaining a uniform polymer-segment density (Decher Reference Decher1997), is evident from layered cationic and anionic polyelectrolytes (Schönhoff Reference Schönhoff2003). As demonstrated by Schwarz & Schönhoff (Reference Schwarz and Schönhoff2002), the (hydrated) thicknesses for poly(allylamine hydrochloride) (cationic) and poly(sodium 4-styrenesulfate) (anionic) deposited on microspheres is 
${\approx }3$ nm with the surface 
$\zeta$-potential alternating between 
$\pm \approx 45$ mV. Much thicker layers have been achieved for multilamellar films. For example, Farhat & Hammond (Reference Farhat and Hammond2005) demonstrated layering of poly(ethylene oxide) and poly(acrylic acid) for fuel cells. These furnished conductivities up to 
$0.1\ {\rm mS}\ {\rm cm}^{-1}$ for films under ambient conditions (
$50$ % relative humidity) comprising some 
$40$ bilayers with a total thickness 
${\sim }9\ \mathrm {\mu }$m (
${\approx }225$ nm per bilayer).
Drawing specifically on impedance spectroscopy, Durstock & Rubner (Reference Durstock and Rubner2001) synthesised sequential adsorbed layers of poly(allylamine hydrochloride) with poly(acrylic acid) and sulfonated polystyrene. The effects of NaCl in the dipping solutions, film hydration (wet vs dry) and temperature on the dielectric spectra were explored. Large increases in the conductivity and real part of the dielectric permittivity with respect to increasing temperature and moisture content were interpreted on the basis of the mobility and effective concentration of small ions, determined, in part, by the solution pH and NaCl concentration. Dielectric spectra were modelled as a parallel bulk resistance and capacitance in series with an interfacial capacitance. From the circuit and geometrical parameters (Durstock & Rubner Reference Durstock and Rubner2001, figure 2), a bulk conductivity 
$\sim 2.1 \times 10^{-6}\ {\rm mS}\ {\rm cm}^{-1}$ can be identified with a bulk specific capacitance 
$0.075\ \mathrm {\mu }{\rm F}\ {\rm cm}^{-2}$ and interfacial specific capacitance 
$0.83\ \mathrm {\mu }{\rm F}\ {\rm cm}^{-2}$. Note that a specific (also termed ‘areal’) capacitance of 
$10\ \mathrm {\mu }{\rm F}\ {\rm cm}^{-2}$ corresponds to a charge density 
${\approx }0.02 e$ per nm
$^2$ at a voltage of 
$k_B T / e \approx 25$ mV (
$k_B T$ is the thermal energy and 
$e$ the fundamental charge); for a parallel-plate capacitor, the specific capacitance can be expressed 
$\epsilon _r \epsilon _0 / L$, where 
$L$ is the electrode/plate spacing and 
$\epsilon _r \epsilon _0$ is the dielectric permittivity (relative permittivity 
$\epsilon _r$, vacuum permittivity 
$\epsilon _0$). However, it is not clear how the conductivity and capacitance are related to the lamellar microstructure, since, for example, interfacial capacitance can be ascribed to the sample-electrode interfaces and/or the lamellar polyelectrolyte domains. Indeed, the authors acknowledged the possibility of ion accumulation at internal interfaces.
Abstractions of the foregoing classes of (cavity-doped and lamellar) microstructures are depicted in figure 1. These impart various degrees of heterogeneity in the fixed charge, ion concentrations, hydrodynamic permeability and elasticity, thus shaping coupled microscale relaxation processes. In the present study, these properties are prescribed continuous spatial variations (on periodic domains), thus avoiding the need to implement internal boundary conditions. The models draw on a limited number of empirical parameters, which are bulk property values and microstructural lengths. In the anionic cavity-doped microstructure, the cationic counter-ions electromigrate in the direction of electroosmotic flow (in the absence of a mean pressure gradient), opposite to the direction in which electrical forces act on the negatively charged (anionic) skeleton. In the lamellar microstructure, however, the overall fixed-charge density is zero, and so with a uniform hydrodynamic permeability, there is zero electroosmotic flow, even though the positively and negatively charged domains experience electrical forces that drive them in opposite directions, opposite to their respective counter-ion electromigration fluxes.
Figure 1. Schematics (one-dimensional abstractions) of structured polyelectrolyte microstructures: anionic polyelectrolyte with cavities (top); and lamellar polyelectrolyte with anionic (
$-$) and cationic (
$+$) enriched domains (bottom). Arrows identify the direction of electrical forces exerted on mobile and fixed (bound to the skeleton) ions under a positive (left to right) electric field. Lines under each schematic depict fixed (prescribed, green) and mobile (calculated, black) equilibrium charge-density spatial modulation.
The present work explores the roles of electrical polarisation and electroosmosis, based on a continuum model that couples ion electromigration, diffusion and advection upon perturbing a nonlinear equilibrium state by the application of a weak uniform electric field and/or pressure gradient. This is in the spirit of the models of solid-state ionic conductors pioneered by Chang & Jaffé (Reference Chang and Jaffé1952) and Macdonald (Reference Macdonald1953), but bringing coupled fluid and elastic dynamics, which will be demonstrated to play a significant, and notably under-appreciated, role in shaping dynamic conductivity spectra, particularly at low frequencies where extremely high dielectric constants have been identified.
The theory in § 2 presents the model, first in a dimensional non-reduced form that identifies the key physical processes accounted for (§ 2.1). Next, adopting scaled variables in the one-dimensional context, the nonlinear equilibrium problem is identified (§ 2.2), proceeding to address the linearised dynamics under small-amplitude oscillatory forcing (§ 2.3). This furnishes a coupled set of ordinary differential equations, which are solved numerically, from which the current density is obtained and averaged within the periodic unit cell (§ 2.4). The results in § 3 begin with preliminary observations under considerably simplified conditions, namely with a single charge carrier in a uniform medium without advection (§ 3.1); steady-state, uniform fixed charge with electroosmosis (§ 3.2); and layered microstructures without electroosmosis (§ 3.3). The main results are presented in the following subsections, for cavity-doped networks (§ 3.4) and then lamellar networks (§ 3.5). These examine the respective equilibrium states, and then turn to examining the dynamics, comparing spectra (and their Nyquist representations) for elastic networks with their rigid counterparts. The paper concludes with a summary in § 4.
2. Theory
The theory is applied to model ionic transport in polyelectrolyte networks that are illustrated schematically in figure 1. These are one-dimensional and spatially periodic, expediting a calculation of the averaged current density in response to a small-amplitude oscillatory electric field. The ratio of the current density to the strength of the electric field is the effective conductivity, which has real and imaginary parts that may be converted into complex-valued impedance and dielectric permittivity spectra. The microstructures of the polyelectrolytes are principally distinguished by the spatial distribution of the fixed charge density, and the accompanying hydrodynamic permeability and elastic stiffness. The equilibrium and dynamics of the mobile anions and cations in the electrolyte solution that bathes the network are captured by what is often termed (in the colloids literature) the ‘standard electrokinetic model’ (O'Brien & White Reference O'Brien and White1978; Delacey & White Reference Delacey and White1981, Reference Delacey and White1982; Hollingsworth & Saville Reference Hollingsworth and Saville2003), but augmented with Brinkman fluid dynamics, i.e. averaged porous-medium flow that is subject to viscous tractions and hydrodynamic-drag body forces (Brinkman Reference Brinkman1947; Durlofsky & Brady Reference Durlofsky and Brady1987). Note the that fluid dynamics in the standard electrokinetic model are of low-Reynolds-number viscous flow with electrical body forces and, possibly, temporal inertia. In the present work, a polyelectrolyte network (charged Brinkman porous medium) is endowed with elasticity, so the charge that is fixed to its skeleton may contribute to the electrical current.
2.1. Model
The one-dimensional model described here is derived from the following differential conservation relationships governing the electrostatic potential 
$\psi (\boldsymbol {x},t)$, 
$N$ mobile ion concentrations 
$c_i(\boldsymbol {x},t)$, fluid velocity 
$\boldsymbol {u}(\boldsymbol {x},t)$ and pressure 
$p(\boldsymbol {x},t)$ and network displacement 
$\boldsymbol {v}(\boldsymbol {x},t)$.
The Poisson equation
\begin{equation} -\epsilon_s \epsilon_0 \nabla^2 \psi = \rho_m + \rho_f, \end{equation}
is the differential form of Gauss's law, which may be considered an electric-flux conservation relationship with 
$\rho _m(\boldsymbol {x},t)$ and 
$\rho _f(\boldsymbol {x},t)$ the mobile and fixed charge densities. The dielectric constant of the medium 
$\epsilon _s \epsilon _0$ (solvent dielectric permittivity 
$\epsilon _s$, vacuum permittivity 
$\epsilon _0$) is taken to be uniform, equal to the value of the solvent. Note that the mobile charge density
\begin{equation} \rho_m = \sum_{i=1}^N z_i e c_i ,\end{equation}
depends on 
$N$ mobile ion concentrations in the fluid, whereas the fixed-charge density 
$\rho _f$ depends on the skeleton dynamics and a fixed-charge conservation relationship, which are presented in the following.
The mobile-ion conservation relationships are
\begin{equation} \frac{\partial c_i}{\partial t} ={-} \boldsymbol{\nabla} \boldsymbol{\cdot}\left(c_i \boldsymbol{u} - D_i \boldsymbol{\nabla} c_i - z_i e \frac{D_i}{k_B T} c_i \boldsymbol{\nabla}\psi \right) \quad (i = 1,\ldots, N), \end{equation}
which have Nernst–Planck fluxes for ion advection, diffusion and electromigration; the 
$N$ independent ion diffusion coefficients 
$D_i = k_B T / \gamma _i$ are the ratio of the thermal energy 
$k_B T$ (Boltzmann constant 
$k_B$, absolute temperature 
$T$) to the ion friction coefficients 
$\gamma _i$; and 
$z_i e$ are the ion charges (valence 
$z_i$, fundamental charge 
$e$). Note that ions may have significantly contrasting mobilities, so, in the present computational framework, the model has been developed for any 
$N$. Examples are presented with 
$N = 2$ ion species (Na
$^+$ and OH
$^{-}$), which have mobilities (in water) varying by a factor of approximately four. It is assumed in the present work that the ion mobilities are uniform and, therefore, unhindered by the network.
Fluid momentum and mass conservation relationships are
\begin{equation} \rho \frac{\partial \boldsymbol{u}}{\partial t} = \boldsymbol{\nabla} \boldsymbol{\cdot} [ \eta (\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla}{\boldsymbol{u}}^T) - p \boldsymbol{I} ] - \frac{\eta}{\ell^{2}} \left(\boldsymbol{u} - \frac{\partial \boldsymbol{v}}{\partial t} \right) - \rho_m \boldsymbol{\nabla} \psi, \quad \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u} = 0. \end{equation}
The fluid is assumed incompressible with temporal inertia (manifesting at high frequency, but with Reynolds number 
$\ll 1$ due to a small microstructural length scale), Newtonian stress (shear viscosity 
$\eta$) and hydrodynamic and electric body forces. Note that the second term on the right-hand side of the momentum equation is a Darcy drag force; this models the hydrodynamic drag that the polymer skeleton exerts on the fluid due to their relative velocity 
$\boldsymbol {u} - \partial \boldsymbol {v} / \partial t$ with prefactor 
$\eta / \ell ^2$ the Darcy permeability. In this paper, the Brinkman screening length 
$\ell (\boldsymbol {x},t)$ (Brinkman Reference Brinkman1947; Durlofsky & Brady Reference Durlofsky and Brady1987) is prescribed to be either spatially modulated (with an accompanying non-uniform charge density) or uniform (with a spatially modulated charge density). When the model is expanded to linear order in the perturbations, it is sufficient to take 
$\ell (\boldsymbol {x},t) \rightarrow \ell (\boldsymbol {x})$.
The polymer-skeleton dynamics are governed by an inertialess momentum conservation relationship with linearly elastic stress:
\begin{equation} 0 = \boldsymbol{\nabla} \boldsymbol{\cdot} [\mu (\boldsymbol{\nabla} \boldsymbol{v} + \boldsymbol{\nabla}{\boldsymbol{v}}^T) + \lambda (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}) \boldsymbol{I}] + \frac{\eta}{\ell^{2}} \left(\boldsymbol{u} - \frac{\partial \boldsymbol{v}}{\partial t} \right) - \rho_f \boldsymbol{\nabla} \psi. \end{equation}
The hydrodynamic (Darcy) body force is equal in magnitude, but with opposite direction, to that on the fluid, whereas the electrical body force depends on 
$\rho _f(\boldsymbol {x},t)$, which is subject to skeleton compressibility. The elastic stress arises from shear and dilation with shear modulus 
$\mu (\boldsymbol {x})$ and Lamé parameter 
$\lambda (\boldsymbol {x})$ (bulk modulus 
$K = \lambda + 2 \mu / 3$) (Landau & Lifshitz Reference Landau and Lifshitz1986). Accordingly, the fixed-charge density is linked to the divergence of the skeleton displacement, which itself is a perturbation from a non-uniform equilibrium state with displacement 
$\boldsymbol {v}^0(\boldsymbol {x})$ and prescribed fixed-charge density 
$\rho _f^0(\boldsymbol {x})$. Note that the model is independent of 
$\boldsymbol {v}^0(\boldsymbol {x})$, since 
$\rho _f^0(\boldsymbol {x})$ is prescribed according to specific network architectures (e.g. cavity-doped or lamellar, as illustrated in figure 1) to ascertain how the equilibrium charge heterogeneity affects charge dynamics.
Assuming a fixed-charge flux 
$\rho _f^0 \partial \boldsymbol {v} / \partial t$ (fixed charge attached to the skeleton), conservation demands
\begin{equation} \frac{\partial \rho_f}{\partial t} ={-} \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \rho_f^0 \frac{\partial\boldsymbol{v}}{\partial t} \right), \end{equation}thus furnishing
\begin{equation} \rho_f = \rho_f^0 - \boldsymbol{\nabla} \boldsymbol{\cdot}[ \rho_f^0 (\boldsymbol{v} - \boldsymbol{v}^0)], \end{equation}which couples the network dynamics/displacement to the fluid dynamics/velocity and electric field.
The solution of the foregoing model equations is set out below in a context where all the dependent variables are functions of position 
$x$ with harmonic time dependence, i.e. having oscillatory temporal dynamics at the same frequency as the forcing. This is accomplished by solving an equilibrium problem (with no external forcing), treating the dynamics (under forcing) as perturbations from equilibrium. The harmonic time dependence (Delacey & White Reference Delacey and White1982) is captured implicitly with prefactors 
${\rm e}^{-{\rm i} \omega t}$, where 
$\omega$ is the angular frequency and 
$i^2 = -1$. In the following sections, the variables will be scaled and used without additional notation (to distinguish between dimensional and non-dimensional counterparts). The characteristic scales are identified (with subscripts 
$c$) as they are introduced, and summarised in table 1.
Table 1. Characteristic scales, representative dimensional and dimensionless parameters, and conversion factors (e.g. with 
$f = 10$ kHz, 
$I = 10$ mM and 
$\rho _f^0 / e = -10$ mM).
$^{\dagger}$The Reynolds number is taken to be zero.
2.2. Equilibrium
In the absence of external forcing, (2.1)–(2.7) furnish an equilibrium electrostatic potential 
$\psi ^0$ that satisfies the Poisson–Boltzmann equation (Russel, Saville & Showalter Reference Russel, Saville and Showalter1989)
\begin{equation} - \psi^0_{xx} = (\kappa a)^2 \left( \sum_{i=1}^N z_i c^0_i + \rho_f^0 \right), \end{equation}
where the fixed-charge density 
$\rho _f^0$ is a prescribed phenomenological function of 
$x$ to model, e.g. cavity-doped and lamellar microstructures (see figure 1). The accompanying equilibrium ion concentrations and pressure gradient are
\begin{equation} c^0_i = c^\infty_i \, {\rm e}^{- z_i \psi^0} \quad \mbox{and}\quad p^0_x ={-} \sum_{i=1}^N z_i c^0_i \psi^0_x, \end{equation}
where the ion concentrations in an external bath 
$c^\infty _i$ satisfy electroneutrality, 
$\sum _{i=1}^N z_i c^\infty _i = 0$. Note that the electrostatic potential 
$\psi$ is now scaled with 
$\psi _c = k_B T / e$, position 
$x$ with a microstructural length scale 
$x_c = a$ (a model parameter that is later defined according to specific network architectures), concentrations 
$c^0_i$ and 
$c^\infty _i$ with 
$c_c = 2 I$ and pressure 
$p$ with 
$p_c = 2 I k_B T$. Here, 
$I = 0.5 \sum _{i=1}^N z_i^2 c^\infty _i$ (dimensional 
$c_i^\infty$) is the ionic strength of the bath, and 
$\kappa ^{-1} = \sqrt {k_B T \epsilon _s \epsilon _0 /(2 I e^2)}$ is a Debye length (based on the electrolyte-bath concentration 
$I$).
2.3. Linearised dynamics
When subjected to a uniform electric field 
$E \boldsymbol {e}_x$, uniform pressure gradient 
$P \boldsymbol {e}_x$ or ion-concentration perturbation of the form 
$C_i + x B_i \boldsymbol {e}_x$ (
$C_i$ is a constant and 
$B_i$ is a uniform gradient), the spatially periodic (hatted) perturbations to the electrostatic potential 
$\hat {\psi }(x) X$, ion concentrations 
$\hat {c}_i(x) X$, fluid velocity 
$\hat {u}(x) X$ and skeleton displacement 
$\hat {v}(x) X$ satisfy the following linear ordinary differential equations emerging from (2.1)–(2.7) (discarding quadratic non-linearities):
$$\begin{gather} -\psi'_{xx} = (\kappa a)^2 \sum_{i=1}^N z_i c_i' - (\kappa a)^2 \rho^0_f v'_x - (\kappa a)^2 \rho^0_{f,x} v', \end{gather}$$
$$\begin{gather}- {\rm i} \varOmega \frac{\nu}{D_i} c'_i ={-} {Pe}_i u' c^0_{i,x} + c'_{i,xx} + z_i (c^0_{i,x} \psi'_x + c'_{i,x} \psi^0_x + c^0_i \psi'_{xx} + c'_i \psi^0_{xx}), \quad i = 1,\ldots,N, \end{gather}$$
$$\begin{gather}u' [- {\rm i} \varOmega + (a / \ell)^2] + (a/ \ell)^2 i (\omega a/u_c) v' ={-}p'_x - \rho^0_m \psi'_x - \sum_{i=1}^N z_i c'_i \psi^0_x, \end{gather}$$
$$\begin{gather}u'_x = 0, \end{gather}$$
$$\begin{gather}(a / \ell)^2 [u' + i (\omega a/u_c) v' ] + \frac{M}{2 I k_B T} v'_{xx} + \frac{M_x}{2 I k_B T} v'_{x} - \rho_f^0 \psi'_x + (\rho_f^0 v'_x + v' \rho_{f,x}^0) \psi^0_x = 0 , \end{gather}$$where the primed variables (perturbations from equilibrium) are
\begin{equation} \left.\begin{array}{c} \psi' = \psi - \psi^0 = \hat{\psi}(x) X - x E,\\ c'_i = c_i - c_i^0 = \hat{c}_i(x) X + C_i + x B_i,\\ p' = p - p^0 = \hat{p}(x) X + x P,\\ u' = \hat{u}(x) X,\\ v' = v- v^0 = \hat{v}(x)X. \end{array}\right\} \end{equation}
This dynamic model is solved with 
$X \in \{E, P, B_i,C_i \}$, with only one of these forcing variables being non-zero (conveniently 
$X = 1$) at a time. The solution for any linear combination of such responses may then be constructed by superposition.
Equations (2.10)–(2.14) are the one-dimensional (scaled, linearised/perturbed) Poisson equation, ion-conservation equations (with Nernst–Planck fluxes) and fluid and network momentum equations. In the one-dimensional context, the incompressible fluid has a velocity 
$\hat {u} X$ that must be uniform. The electrostatic potential, pressure and ion-concentration perturbations are scaled as their equilibrium counterparts above; 
$t$ with 
$t_c = \omega ^{-1}$; and the fluid velocity 
$u$ and network displacement 
$v$ with 
$u_c = 2 I k_B T a / \eta$ and 
$v_c = a$, respectively. The applied electric field 
$E$ is scaled with 
$\psi _c / x_c = k_B T / (e a)$; concentration 
$C_i$ and concentration gradients 
$B_i$ with 
$c_c = 2 I$ and 
$c_c / x_c = 2 I / a$, respectively; and the pressure gradient 
$P$ with 
$p_c / x_c = 2 I k_B T / a$.
Note that the fluid momentum balance includes all the temporal inertia, whereas the network dynamics are assumed to be inertialess due to the characteristically low polymer volume fraction in swollen hydrogels. The fluid and network dynamics are coupled by hydrodynamic friction that is linear in the relative velocity of the two phases (skeleton velocity 
$- \textrm {i} \omega v$, dimensional). The fluid is subject to the electrical body force of the mobile free charge and the network is subject to the electrical body force of the fixed charge and elastic stresses (shear modulus 
$\mu$ and Lamé parameter 
$\lambda$). In the one-dimensional context, the elastic moduli combine into a single function of position 
$M(x) = 2 \mu (x) + \lambda (x)$, which is termed the 
$p$-wave modulus (appearing as a ratio with the osmotic pressure 
$2 I k_B T$). Other dimensionless groups include 
$\varOmega = \omega a^2 / \nu$ (squared Womersley number, fluid kinematic viscosity 
$\nu = \eta / \rho$) and 
${Pe}_i = 2 I k_B T a^2 / (\eta D_i)$ (ion Péclet numbers). The Brinkman screening length, similarly to the fixed-charge density and 
$p$-wave modulus, is a prescribed function of position 
$\ell (x)$ (according to the specific network architecture).
2.4. Solution and ion current
The nonlinear equilibrium (2.8) is solved for 
$\boldsymbol {y}^0(x) = [\psi _x^0, \psi ^0]$ on 
$-L/2 \le x \le L/2$ as two coupled first-order ordinary differential equations (using Matlab's bvp4 or bvp5 functions) with periodic boundary conditions 
$\psi _x^0(-L/2) = \psi _x^0(L/2)$ and 
$\psi ^0(-L/2) = \psi ^0(L/2)$. With these, the linear (2.10)–(2.14) for 
$N = 2$ mobile ion species are solved as 
$10$ coupled first-order ordinary differential equations (also using Matlab's bvp4 or bvp5 functions) for the hatted perturbations:
\begin{equation} \frac{\partial \hat{\boldsymbol{y}}}{\partial x} = \boldsymbol{M}^{{-}1} \boldsymbol{\cdot} (\boldsymbol{f} - \boldsymbol{A} \boldsymbol{\cdot} \hat{\boldsymbol{y}}), \end{equation}
where the components of the matrices 
$\boldsymbol {A}$ and 
$\boldsymbol {M}$ and vector 
$\boldsymbol {f}$ are functions of the components of 
$\boldsymbol {y}^0(x)$ with 
$\hat {\boldsymbol {y}}(x) =[\hat {\psi }_x, \hat {c}_{1,x}, \hat {c}_{2,x}, \hat {v}_x, \hat {\psi }, \hat {c}_1, \hat {c}_2, \hat {p}, \hat {v}, \hat {u}]$ and boundary conditions:
\begin{align} \left.\begin{array}{c} \hat{\psi}({-}L/2) = 0, \quad \hat{\psi}({-}L/2) = \hat{\psi}(L/2), \quad \hat{p}({-}L/2) = 0, \quad \hat{p}({-}L/2) = \hat{p}(L/2),\\ \hat{v}_{x}({-}L/2) = \hat{v}_{x}(L/2), \quad \hat{v}({-}L/2) = \hat{v}(L/2),\\ \hat{c}_{1,x}({-}L/2) = \hat{c}_{1,x}(L/2), \quad \hat{c}_{1}({-}L/2) = \hat{c}_{1}(L/2),\\ \hat{c}_{2,x}({-}L/2) = \hat{c}_{2,x}(L/2), \quad \hat{c}_{2}({-}L/2) = \hat{c}_{2}(L/2). \end{array}\right\} \end{align}
Note that 
$\hat {\psi }(-L/2)$ and 
$\hat {p}(-L/2)$ are arbitrarily set to zero, since the governing equations depend only on the first and second derivatives of these functions.
The (dimensional) volume-averaged, linearised current density is
\begin{equation} \langle \boldsymbol{i} \rangle = \sum_{i = 1}^N z_i e \langle\kern 1.5pt \boldsymbol{j}'_i \rangle - {\rm i} \omega \langle \rho_f^0 \boldsymbol{v}' \rangle + {\rm i} \omega \epsilon_s \epsilon_0 \langle \boldsymbol{\nabla}{\psi}' \rangle, \end{equation}where the advective, diffusive and electromigrative fluxes for the mobile ions are
\begin{equation} \boldsymbol{j}'_i = c^0_i \boldsymbol{u} - D_i \left(\boldsymbol{\nabla} c'_i + c'_i \frac{z_i e}{k_B T} \boldsymbol{\nabla} \psi^0 + c^0_i \frac{z_i e}{k_B T} \boldsymbol{\nabla} \psi' \right), \end{equation}
$- \textrm {i} \omega \rho _f^0 \boldsymbol {v}'$ is the current density from fixed-charge translation, and 
$\textrm {i} \omega \epsilon _s \epsilon _0 \boldsymbol {\nabla }{\psi }'$ is termed the displacement current from dielectric polarisation of the solvent (Delacey & White Reference Delacey and White1982).
From the solutions 
$\boldsymbol {y}^0(x)$ and 
$\hat {\boldsymbol {y}}(x)$, the volume-averaged current density, scaled with 
$i_c = (2 I)^2 k_B T e a / \eta$ and linearised with respect to the perturbations relative to the displacement current is
\begin{align} \langle i \rangle &={-} \frac{{\rm i} \omega \eta}{2 I k_B T} \langle \rho_f^0 \hat{v} \rangle X + \sum_{i=1}^N [ - z_i {Pe}_i^{{-}1} B_i - z_i^2 {Pe}_i^{{-}1} \langle \psi^0_x \rangle (C_i + B_i x) + z_i^2 {Pe}_i^{{-}1} \langle c^0_i \rangle E ] \nonumber\\ &\quad +\sum_{i=1}^N [ z_i \langle c_i^0 \hat{u} \rangle - z_i {Pe}_i^{{-}1} \langle \hat{c}_{i,x} \rangle - z_i^2 {Pe}_i^{{-}1} (\langle \hat{c}_{i} \psi^0_x \rangle + \langle c^0_i \hat{\psi}_x \rangle) ] X, \end{align}
where angled brackets denote volume averages 
$\langle {\cdot } \rangle = L^{-1}\int {\cdot } \mbox {d}\kern 0.06em x$. The volume-averaged (dimensionless) displacement current 
$-\textrm {i} \omega \epsilon _s \epsilon _0 (E - X \langle \hat {\psi }_x \rangle ) \eta /(2 I e a)^2$ is dominated by the term having dimensional counterpart 
$- \textrm {i} \omega \epsilon _s \epsilon _0 E$ (Delacey & White Reference Delacey and White1982); both have been discarded in (2.20). Note that 
$\langle i \rangle$ denotes the averaged current density, otherwise 
$i$ (without angled brackets) denotes the imaginary unit.
The spatial periodicity of 
$\psi ^0$ and 
$\hat {c}_i$ in (2.20) means that 
$\langle \psi ^0_x \rangle = \langle \hat {c}_{i,x} \rangle = 0$. The first term in (2.20) arises from the motion of the fixed charge (dimensional skeleton velocity 
$-\textrm {i} \omega \hat {v} X$). The first sum collects the diffusion and electromigration fluxes due to uniform ion concentrations 
$C_i$, ion-concentration gradients 
$B_i$, and electric field 
$E$ acting on the equilibrium fields and the second sum collects the coupled advection, diffusion and electromigration perturbation currents generated by the coupling of the equilibrium and hatted perturbations. Recall that the calculations are undertaken with only one non-zero 
$X \in \{E, P, B_i,C_i \}$ at a time. The response to macroscale concentration perturbations (
$B_i$, 
$C_i$) are beyond the scope of the present study, and so may be investigated elsewhere. Accordingly, the response to a bulk electric field or pressure-gradient perturbation may be superposed to meet macroscale boundary conditions for which 
$\langle u \rangle = \langle \hat {u}^P \rangle P + \langle \hat {u}^E \rangle E = 0$. Linear superposition of electric-field- and pressure-driven current and volume fluxes gives
\begin{equation} \langle i \rangle = \alpha_{IE} E - \alpha_{IP} P, \quad \langle u \rangle = \alpha_{UE} E - \alpha_{UP} P, \end{equation}
so under conditions where the average velocity 
$\langle u \rangle = 0$, we have
\begin{equation} P / E = \alpha_{UE} / \alpha_{UP} \quad \text{and}\quad \langle i \rangle / E = \alpha_{IE} - \alpha_{IP} \alpha_{UE} / \alpha_{UP}. \end{equation}
According to the Onsager principle (Thovert, Shapiro & Adler Reference Thovert, Shapiro and Adler1996; Doi Reference Doi2013), 
$\alpha _{UE} = \alpha _{IP}$, giving
\begin{equation} \langle i \rangle / E \rightarrow \alpha_{IE} - \alpha_{UE}^2 / \alpha_{UP}. \end{equation}
Since the dimensional hydrodynamic permeability 
$\alpha _{UP} \sim \ell ^2 /\eta > 0$ and 
$\alpha _{UE}^2 / \alpha _{UP} > 0$, the effective conductivity 
$\langle i \rangle / E$ when 
$\langle u \rangle = 0$ is less than the effective conductivity 
$\alpha _{IE}$ when 
$P = 0$. However, Onsager reciprocity furnishing (2.23) appears to hold only for networks with zero net displacement, e.g. perfectly rigid networks, and compliant networks bearing zero net charge and uniform hydrodynamic permeability. For this reason, the current densities with 
$\langle u \rangle = 0$ below are reported according to the superposition (2.21a,b) rather than (2.23).
In addition to the complex conductivity/current density, we examine the complex-dielectric permittivity, defined as
\begin{equation} \epsilon = \epsilon' + i \epsilon'' = \frac{\sigma' + {\rm i} \sigma''}{- {\rm i} \omega \epsilon_0}, \end{equation}
so 
$\epsilon ' = - \sigma '' / (\omega \epsilon _0)$ and 
$\epsilon '' = \sigma ' / (\omega \epsilon _0)$, where 
$\sigma = \sigma ' + \textrm {i} \sigma ''$ is the (dimensional) complex conductivity.
In the results section, a dimensionless dielectric permittivity is computed by dividing the scaled conductivity with a scaled angular frequency 
$\omega \kappa ^{-2} / D_1$ (e.g. diffusivity/mobility of the Na
$^+$ counter-ion). This is equivalent to scaling 
$\epsilon$ with 
$\epsilon _c = i_c e a / (k_B T \epsilon _0 \kappa ^2 D_1) = 2 I a^2 k_B T \epsilon _s / (\eta D_1) = 12 {\rm \pi}I a^2 a_1 \epsilon _s$, where 
$a_1$ is defined by the Stokes–Einstein formula 
$D_1 = k_B T / (6 {\rm \pi}\eta a_1)$ and 
$I a^2 a_1$ may be considered the number of electrolyte ions in a volume of the order 
$a^2 a_i$. Moreover, the dimensionless current density 
$\langle i \rangle$ and dielectric permittivity 
$\epsilon$ are explicitly multiplied by the harmonic-mean ion Péclet number 
${Pe}_h = 2 / ({Pe}_1^{-1} + {Pe}_2^{-1})$. Note that 
${Pe}_h \langle i \rangle / E$ with dimensionless 
$E = 1$ is precisely the reduced conductivity 
$\sigma / \sigma _\infty$, where 
$\sigma _\infty$ is the bulk electrolyte conductivity.
3. Results
Calculations have been undertaken with a variety of prescribed functions 
$\rho _f^0(x)$, 
$\ell (x)$ and 
$M(x)$ to model network microstructures having cavities, alternating and interpenetrating layers of cationic and anionic polyelectrolytes, and networks with spatially modulated/segregated anionic fixed charge. Before explicitly examining cavity-doped and lamellar microstructures (figure 1) in §§ 3.4 and 3.5, respectively, some preliminary observations are highlighted in §§ 3.1–3.3 under three simplified conditions. First, the model is applied for a single mobile charge carrier (counter-cation) in a uniform polyelectrolyte, without advection; this demonstrates compatibility with already established theory. Second, the d.c. conductivity for a uniform polyelectrolyte (with added salt) is addressed, including steady advection (zero pressure gradient). Third, the d.c. conductivity for uniform layers (alternating cavity and uniform polyelectrolyte domains with cavity fraction 
$\phi$) is addressed, albeit without advection.
The results are presented in terms of dimensionless/scaled quantities; for example, a scaled current density (per unit of electric field strength) that is equivalent to a scaled complex conductivity. This helps to identify universal, qualitative characteristics of the spectra and Nyquist representations. Formulae for the factors that convert the scaled current density 
${Pe}_h \langle i \rangle$ to, for example, a dimensional conductivity, or the scaled dielectric permittivity 
${Pe}_h \epsilon$ to a relative dielectric permittivity are provided in the text and recapitulated, for convenience, in table 1.
Recall that the scaled current density, computed with unit scaled electric field, i.e. 
${Pe}_h \langle i \rangle / E$ with 
$E = 1$ is equivalent to the medium conductivity 
$\sigma$ divided by the value for the bulk electrolyte 
$\sigma _\infty$, i.e. 
${Pe}_h \langle i \rangle \equiv \sigma / \sigma _\infty$, as adopted in §§ 3.4 and 3.5. Moreover, the real or imaginary part of this complex quantity will be explicitly identified in the figure captions. A scaled (complex) dielectric permittivity is computed as 
$\epsilon = (\langle i \rangle / E) / (\omega \kappa ^{-2} / D_1)$ with 
$E = 1$, and so may be converted into a dielectric permittivity by multiplying by a factor 
$\epsilon _c$. However, to be consistent with the foregoing scaling of the current density/conductivity, the scaled dielectric permittivity (real part) is plotted as 
${Pe}_h \epsilon ^\prime$, and so may be converted into a relative dielectric permittivity by multiplying by a factor 
$\epsilon _c / {Pe}_h$.
To ensure that the parameters are in, or not too far outside of, physically relevant ranges, dimensionless parameters are mostly varied by systematically changing dimensional variables, similarly to an (ideal) experimental parametric study. Of course, varying one such dimensional variable may affect several dimensionless model parameters.
Table 1 summarises the dimensional scales and variables, and the primary dimensionless parameters arising from the model equations, as presented in (2.8) and (2.10)–(2.14). Other dependent dimensionless parameters will be highlighted in the text as needed. Note that the characteristic fluid velocity emerges from a balance of an 
$O[2 I e k_B T / (e a)]$ electrical body force, itself based on the product of a characteristic charge density 
$2 I e$ and electric field strength 
$k_B T / (e a)$, and an 
$O(\eta u_c / a^2)$ viscous drag (per unit volume). As indicated by the table, the Reynolds number 
$Re = u_c a / \nu = 2 I k_B T a^2 / (\eta \nu ) \ll 1$, thus justifying having neglected the nonlinear inertial term in the momentum equation. Of the four dimensionless parameters depending on the frequency, the (squared) Womersley number is very small at subgigahertz frequencies. The others, which reflect ion-concentration and, more notably, Darcy-dissipation dynamics, are suggested by their values in table 1 to be relevant at at-least megahertz frequencies, accessible with impedance and dielectric spectrometers. Several other (dependent) relaxation time scales are identified in table 1 (bottom section), since these may manifest through the coupling of fluid momentum, ion, charge and elastic dynamics: two (involving 
$M$) depend on the network stiffness, and one (involving 
$\ell$) depends on the network hydrodynamic permeability, both of which may vary in space according to the specific microstructure under investigation.
3.1. Charge dynamics, uniform medium, advection free
If the equilibrium state is uniform and there is no advection, the one-dimensional Poisson equation and ion-conservation equations reduce to
\begin{equation} - {\rm i} \varOmega \frac{\nu} {D_i} (\hat{c}_i + C_i / X) = \hat{c}_{i,xx} - z_i c^0_i (\kappa a)^2 \sum_{j=1}^N z_j (\hat{c}_j + C_j / X). \end{equation}
If we further prescribe only one mobile charge carrier (counter-ions, so 
$c^0_1 = 1$), then its concentration disturbance satisfies (combining 
$\hat {c}_1 X + C_1$ into one dependent variable)
\begin{equation} [z_1^2 (\kappa a)^2 - {\rm i} \varOmega \nu / D_1] \hat{c}_1 = \hat{c}_{1,xx}, \end{equation}
which upon rescaling 
$x$ as 
$y = x \sqrt {z_1^2 (\kappa a)^2 - \textrm {i} \omega a^2 / D_1}$ may be written
\begin{equation} \hat{c}_1 = \hat{c}_{1,yy}. \end{equation}
Then, defining 
$a$ as a macroscopic sample size 
$L$, and identifying 
$D_1 = k_B T \mu _1$ with 
$\mu _1$ the ion/electron mobility, this is equivalent to the charge transport model of Macdonald (Reference Macdonald1953). Accordingly, writing in terms of the dimensional position 
$x$ gives 
$y = x / L = \kappa x \sqrt {1 - \textrm {i} \omega \tau }$, where
\begin{equation} \tau = \frac{\kappa^{{-}2}}{z_1^2 D_1} = \frac{\epsilon_s \epsilon_0}{2 I e^2 z_1^2 \mu_1}, \end{equation}
which has been termed the dielectric relaxation time (Coelho Reference Coelho1991). Here, this resurfaces as the characteristic time for charge diffusion on the Debye length scale 
$\kappa ^{-1}$ (Hollingsworth & Saville Reference Hollingsworth and Saville2003).
3.2. Direct current conductivity, uniform medium
For a uniform network (with 
$z$–
$z$ electrolyte, counter-ion valence 
$z_1$) subjected to a steady electric field, the scaled d.c. conductivity is the sum of electromigrative and electroosmotic-advective terms:
\begin{align} \sigma / \sigma_\infty &= {Pe}_h \langle i \rangle / E = 0.5 z^2 {Pe}_h [{Pe}_1^{{-}1} \, {\rm e}^{-\mbox{asinh}(\rho_f^0 / z_1)} + {Pe}_2^{{-}1} \, {\rm e}^{\mbox{asinh}(\rho_f^0 / z_1)}] \nonumber\\ &\quad - 0.5 z_1 \rho_f^0 {Pe}_h (\ell / a)^2 [{\rm e}^{-\mbox{asinh}(\rho_f^0 / z_1)} - \, {\rm e}^{\mbox{asinh}(\rho_f^0 / z_1)}], \end{align}
where the bulk-electrolyte conductivity 
$\sigma _\infty = \sigma _c / {Pe}_h = I e^2 (z_1^2 D_1 + z_2^2 D_2) / (k_B T)$. Note that 
$\exp {[\pm \mbox {asinh}(x)]} =\sqrt {1 + x^2} \pm x$, and so (3.5) may be written (now in terms of dimensional variables)
\begin{equation} \sigma = \frac{\langle i \rangle}{E} = \frac{I (z e)^2}{k_B T} [D_1\, {\rm e}^{-\mbox{asinh}({\rho_f^0}/{2 I e z_1})} +D_2 \, {\rm e}^{\mbox{asinh}({\rho_f^0}/{2 I e z_1})} ] + (\rho_f^0)^2 \ell^2 / \eta. \end{equation}
Of course, this is independent of the length scale 
$x_c = a$ embedded in the Péclet numbers, also highlighting the advective contribution from the counter-ions as being proportional to 
$(\rho _f^0)^2$, whereas the electromigrative contribution is more generally proportional to 
$\rho _f^0$.
Equation (3.5) is plotted in figures 2(a) and 2(b) for networks that are in equilibrium with high- and low-ionic-strength (NaOH) electrolytes, respectively. The range of 
$-\rho _f^0 / (2 I e)$ in each panel ensures that the range of 
$\rho _f^0$ is physically motivated for the respective fixed value of 
$I$. Blue lines are the limit 
$\ell \rightarrow 0$ in which hydrodynamic friction completely arrests electroosmotic flow; here the conductivity reflects the composition of the electrolyte in the network, as prescribed by Donnan equilibrium with the bath. Increasing 
$\ell$ promotes electroosmotic flow that increases 
$\sigma$. When the fixed-charge density is low, the current density is that for electrolyte in the network with the same composition as the bulk electrolyte, thus furnishing 
$\sigma / \sigma _\infty = 1$. As expected, increasing the network charge tends to increase 
$\sigma$, due to the added counter-ion (Na
$^+$). However, because the counter-ion (Na
$^+$) has a lower mobility than the co-ion (OH
$^-$), this manifests in 
$\sigma / \sigma _\infty < 1$ when 
$|\rho _f^0| / (2 I e) \lesssim 2$.
Figure 2. Scaled current density/conductivity vs fixed-charge density (scaled with 
$-2 I e$) in a uniform medium with asymmetric 
$1$–
$1$ electrolyte having 
${Pe}_1 / {Pe}_2 = D_2 / D_1 \approx 3.97$ and Brinkman length 
$\ell = 0$ (blue), 
$0.5$ (red), 
$1$ (yellow), 
$2$ (violet) and 
$10$ nm (green) (increasing upward): (a) 
$I = 100$ mM; (b) 
$I = 1$ mM. Blue lines (bottom) are the advection-free limit.
3.3. Direct current conductivity, advection-free, layered medium
Now consider a medium that comprises alternating layers of uniform polyelectrolyte and pure electrolyte, the latter representing cavities with volume fraction 
$\phi = 2 a / L$. Approximating the ionic conductivity inside and outside the cavity phase as being uniform with values 
$\sigma _i$ and 
$\sigma _o$, respectively, a continuous current density requires 
$\sigma _i (E + E_i) = \sigma _o (E + E_o)$, where 
$E_i = -{\rm \Delta} \psi / a$ and 
$E_o = {\rm \Delta} \psi / (L - 2 a)$ are the periodic disturbances to the applied electric field 
$E$ inside and outside the cavities. With a continuous electrostatic potential, these furnish
\begin{equation} {\rm \Delta} \psi / a = \frac{\sigma_i / \sigma_o -1}{ (a/L) / (1 - 2 a / L) + \sigma_i / \sigma_o} E, \end{equation}and so the effective conductivity of the medium is
\begin{equation} \sigma = \sigma_i (E + E_i) / E = \sigma_i - \frac{\sigma_i {\rm \Delta} \psi}{a E} = \sigma_i \frac{\phi + 2 (1 - \phi)}{ \phi + 2 (1 - \phi)\sigma_i / \sigma_o}. \end{equation}
If we assume that the conductivities reflect Donnan equilibrium of a 
$z$–
$z$ electrolyte partitioned between the cavity (also bulk electrolyte) and polyelectrolyte phases, also neglecting electroosmotic flow, then this model interpolates between (dimensional variables)
\begin{align} \left.\begin{aligned}\sigma_i &=
\sigma_\infty = \frac{I (z e)^2}{k_B T} (D_1 + D_2) \quad
\mbox{and}\nonumber\\
\sigma_o &= \frac{I (z e)^2}{k_B T}
[D_1 \, {\rm e}^{-\mbox{asinh}({\rho_f^0}/{2 I e
z_1})} + D_2 \, {\rm e}^{\mbox{asinh}({\rho_f^0}/{2 I e
z_1})}].\end{aligned}\right\}
\end{align}
The approximation requires 
$\kappa a \gg 1$ and 
$\kappa (L - 2 a) \gg 1$, so that the cavity and polyelectrolyte domains do indeed have uniform ion concentrations, thus furnishing an effective conductivity 
$\sigma$ that is independent of 
$a$ with fixed 
$\phi = 2 a / L$. Equation (3.8), which has a small-
$\phi$ expansion
\begin{equation} \sigma = \sigma_o \left[ 1 - \frac{(1 - \sigma_i / \sigma_o)}{2 \sigma_i / \sigma_o} \phi + O(\phi^2) \right], \end{equation}
is plotted in figure 3 vs the cavity volume fraction with 
$\sigma _i = \sigma _\infty$, ion mobilities furnishing 
${Pe}_1 / {Pe}_2 = D_2 / D_1 \approx 3.97$ and various fixed charge densities. In figure 3(a), the bulk ionic strength is high (
$I = 100$ mM) and the effective conductivity reflects the changing composition of the electrolyte that is in equilibrium with the polyelectrolyte. In figure 3(b), the bulk ionic strength is low (
$I = 1$ mM), so the effective conductivity is dominated by the counter-ion concentration and mobility. Note that the values of 
$\rho _f^0 / (2 I e)$ adopted in each panel ensure that the underlying values of 
$\rho _f^0$ are physically acceptable. The non-monotonic variation of 
$\sigma$ with 
$\rho _f^0$ in figure 3(a) is consistent with expectations from figure 2(a) with 
$\ell = 0$, as is the monotonic variation in figure 3(b) with figure 2(b) (for 
$-\rho _f^0 / (2 I) \gg 1$).
Figure 3. Scaled current density/conductivity vs cavity volume fraction in an hydrodynamically impermeable, lamellar polyelectrolyte with asymmetric 
$1$–
$1$ electrolyte having 
${Pe}_1 / {Pe}_2 = D_2 / D_1 \approx 3.97$: (a) 
$I = 100$ mM; 
$- \rho _f^0 / (2 I e) = 1$ (blue), 
$0.1$ (red), 
$0.2$ (yellow), 
$0.4$ (violet) and 
$0.6$ (green) (solid, increasing downward); 
$- \rho _f^0 / (2 I e) = 1$ (blue), 
$1.25$ (red), 
$1.5$ (yellow), 
$1.75$ (violet) and 
$2$ (green) (dashed, increasing upward);(b) 
$I = 1$ mM; 
$- \rho _f^0 / (2 I e) = 0$ (blue), 
$50$ (red), 
$100$ (yellow), 
$150$ (violet) and 
$200$ (green) (increasing upward).
3.4. Cavity-doped networks
For an anionic network containing cavities, the fixed-charge density is prescribed
\begin{equation} \rho_f^0(x) = \rho_f^0 \{1 - \text{erfc}[(x / a - 1) / \delta]/2 + \text{erfc}[(x / a +1)/ \delta)]/2 \}, \end{equation}
where 
$-L/2 \le x \le L/2$ and 
$\delta \sim \ell / a \lesssim 1$ (a parameter controlling the sharpness of the transition) so that the cavity volume fraction 
$\phi \approx 2 a / L$. Unless noted otherwise, calculations were undertaken with 
$\delta = \ell / a$, so the diffuseness of the cavity wall scales with 
$\ell$. The Brinkman length 
$\ell (x)$ is set inversely proportional to 
$|\rho _f^0(x)|^{1/2}$, plateauing to a constant 
$\ell$ outside the cavity, and the 
$p$-wave modulus 
$M(x)$ is set proportional to 
$| \rho _f^0(x) |$, plateauing to a constant 
$M$ outside the cavity. These, reflect a hydrodynamic/Darcy drag and elastic moduli that scale with the network charge density. Primary dimensional parameters here are the cavity width 
$2 a$, bulk ionic strength 
$I$ and (anionic) fixed-charge density 
$\rho _f^0 < 0$. The solvent is taken to be water at 
$T = 298$ K with cations and anions (Na
$^+$ and OH
$^-$) having diffusion coefficients that are ascertained from their limiting molar conductivities (
${\approx }50.1$ and 
$199\ \textrm {S}\ \textrm {cm}^{-2}\ \textrm {mol}^{-1}$, respectively). For a highly charged network with dimensional 
$|\rho _f^0| \gtrsim 2 I e$, the intervening electrolyte is dominated by the Na
$^+$ counter-cation. This arises from the Donnan equilibrium between the network and bath, as captured by the Poisson–Boltzmann equation. In a weakly charged network with dimensional 
$|\rho _f^0| \lesssim 2 I e$, the intervening electrolyte is more uniform with equal concentrations of cations and anions.
An equilibrium charge density profile according to (3.11) is plotted as the green line in figure 4, with accompanying electrostatic potential, ion concentrations, mobile and total charge density. Note that the net charge density (black line) is predominantly non-zero at the cavity–polyelectrolyte interfaces (
$x / a = \pm 1$), with positively and negatively charged Debye layers on the inside and outside surfaces, respectively. The electrolyte inside the cavity has equal concentrations of cations and anions (red and yellow lines), whereas the polyelectrolyte is dominated by an excess of counter-cations (red line). Note that the mobility of the counter-cation is lower than of the co-anions, the ratio being 
$D_1 / D_2 = 1.3 / 5.3 \approx 0.25$. Here, Donnan equilibrium between the cavity and polyelectrolyte domains gives rise to a scaled Donnan potential that is of the order of the bulk scaled fixed-charge density, corresponding to 
${\approx }-25$ mV. This proportionality holds only in the Debye–Hückel regime (dimensional 
$|\psi | \lesssim k_B T / e \approx 25$ mV); otherwise the nonlinearity of the Poisson–Boltzmann equation (counter-ion screening) furnishes a scaled Donnan potential that is much smaller than of the scaled fixed-charge density when 
$|\rho _f^0| / (2 I e) \gtrsim 1$ (e.g. 
$I \sim 1$ mM, 
$|\rho _f^0| / e \sim 100$ mM).
Figure 4. Scaled equilibrium electrostatic potential (
$\psi ^0$, blue), ion concentrations (
$c_i^0$, red and yellow), mobile (
$\rho _m^0$, violet), fixed (
$\rho _f^0$, green) and total (
$\rho _m^0 + \rho _f^0$, black) charge density for a cavity-filled polyelectrolyte [fixed-charge density according to (3.11)]: 
$a = 5$ nm, 
$L = 5 a$, 
$I = 100$ mM (
$\kappa a \approx 5.1, {Pe}_1 \approx 9.5, {Pe}_2 \approx 2.3$) and 
$\rho _f^0 / (2 I e) = -1$.
Perturbations to the electrostatic potential gradient, ion-concentration gradients, and the skeleton displacement and velocity when subjected to an electric field 
$E$ (with 
$P = 0$) at a frequency 
$f = 100$ kHz are shown in figure 5. In this example, the fixed-charge density is high (
$- \rho _f^0 / e = 100$ mM) and the bulk electrolyte concentration is low (
$I = 1$ mM), with the bulk 
$p$-wave modulus set to 
$M / (2 I k_B T) = 10$, corresponding to 
$M \approx 50$ kPa. Here, the dominant ion-concentration perturbation is that of the counter-cation (red lines). This reflects its large equilibrium concentration, as required to satisfy electroneutrality. It is this perturbation that dominates the free charge and accompanying electrostatic potential perturbations, which reflect an induced dipole. In one spatial dimension, electrostatic interaction of the cavities is very strong, reflecting a linear variation in the electrostatic potential with respect to position. As shown in figure 5(a), the spatially periodic electric-field perturbation is negative inside the cavity and positive outside. The network displacement profiles in figure 5(c) reflect dynamics of the soft/compliant interfacial region responding to electrical, hydrodynamic and elastic stresses. The accompanying network velocity in figure 5(d) reflects the skeleton's (negative) electrophoretic mobility, which takes the largest values in regions where the fixed-charge density and elastic moduli are low.
Figure 5. Real (a) and imaginary (b) parts of the perturbation gradients of electrostatic potential (blue), counter-cation concentration (red) and co-anion concentration (yellow) for compressible (
$M / (2 I k_B T) = 10$) cavity-filled polyelectrolytes. (c,d) Real (solid) and imaginary (dashed) parts of the scaled skeleton displacement (c) and velocity (d). All the hatted quantities are dimensionless, calculated with dimensionless 
$E = 1$. Other parameters: 
$f = 100$ kHz, 
$a = 5$ nm, 
$\ell / a = 1/5$, 
$L /a =5$, 
$I = 1$ mM (
$\kappa a \approx 0.51$, 
${Pe}_1 \approx 0.093$, 
${Pe}_2 \approx 0.023$), 
$\rho _f^0 / (2 I e) = -100$ and 
$P = 0$.
Note that a crude estimate of the skeleton's in-phase (real) part of the electrophoretic mobility (dimensional 
$-\textrm {i} \omega v / E$) may be estimated by balancing an 
$O(\rho _f^0 E)$ body force with an 
$O(\textrm {i} \omega \hat {v} E \eta / \ell ^2)$ Darcy drag (
$- \textrm {i} \omega \hat {v} E$ is the skeleton velocity), thus giving a mobility 
$\sim \textrm {i} \omega \hat {v} \sim \rho _f^0 \ell ^2 / \eta$. In terms of the scaled ordinate in figure 5(d), this has a magnitude 
$\sim (\ell / a)^2 \rho _f^0 / (2 I e) \approx -4$. Similarly, balancing an 
$O(\hat {v} E M / a^2)$ elastic force with the foregoing Darcy drag furnishes an out-of-phase (imaginary) part of the mobility 
$- \textrm {i} \omega \hat {v} \sim M \ell ^2 / (\eta a E)$, which, in terms of the scaled ordinate, has a magnitude 
$\sim (\ell / a)^2 M / (2 I k_B T) \approx 0.4$. Thus, the electric field imparts a negatively signed velocity (on account of the negatively signed network charge) when the field is in its positive phase, relaxing under the elastic stresses with a positive velocity when the electric field is changing its sign, as captured quantitatively by the real (solid) and imaginary (dashed) parts of the scaled velocity in figure 5(d).
The real part of the averaged current density and dielectric permittivity for perfectly rigid cavity-filled networks with three cavity volume fractions are shown in figure 6. Note that the dashed lines are under electric-field forcing (
$E = 1$ with 
$P = 0$), which permits a non-zero electroosmotic advection, whereas the solid lines are constructed so that the average fluid velocity vanishes (
$\langle u \rangle = 0$). The real part of the current density increases from a low-frequency plateau that reflects quasi-steady electrostatic polarisation. With increasing frequency, this polarisation is suppressed by the diminishing time available for polarisation, thus causing the current density to approach high-frequency plateaus. These plateaus reflect the overall composition of the electrolyte, which is increasingly dominated by the Na
$^+$-rich polyelectrolyte phase (lower ion mobility) with decreasing cavity volume fraction. Recall, this was captured by the approximate d.c. conductivity model furnishing (3.8). That the low- and high-frequency plateaus both vary in proportion to 
$L / a = 2 / \phi$ suggests that the low-frequency polarisation effect is approximately independent of the cavity volume fraction. This is consistent with the cavities being well separated with respect to the Debye length, i.e. 
$\kappa (2 a - L) \gg 1$.
Figure 6. Real part of the scaled current density (a) and scaled dielectric permittivity (b) vs scaled angular frequency: 
$L / a = 2.5$ (blue), 
$5$ (red), 
$10$ (yellow) for rigid cavity-filled polyelectrolytes. Solid and dashed lines correspond to 
$\langle u \rangle = 0$ and 
$P = 0$, respectively. Other parameters: 
$a = 10$ nm, 
$\ell = 1$ nm, 
$I = 1$ mM and 
$\rho _f^0 / (2 I e) = -100$. (c) and (d) are Nyquist plots of the scaled current density (as a measure of the complex admittance) and its reciprocal (as a measure of the complex impedance), respectively. Frequency in (c) and (d) increases from left to right, and from right to left, respectively, spanning the range 
$f = 1$–
$10^8$ kHz. The negative imaginary parts at high frequency in (c) and (d) with 
$P = 0$ (dashed lines) (absent when 
$\langle u \rangle = 0$, solid lines) reflect temporal fluid inertia.
For rigid networks, the characteristic relaxation time reflects ion diffusion. For the highly charged networks (in low-ionic-strength electrolyte) in figure 6(a), the diffusion length scale is set by a Debye length that depends on the counter-cation concentration 
$| \rho _f | / e$, not the ionic strength of the electrolyte bath 
$I$. This is demonstrated by the spectra in figure 7(a) for three cavity sizes 
$a$ (with fixed 
$L / a$). These exhibit dispersions that are reasonably well centred on 
$\omega \kappa ^{-2} / D_1 \approx 60$, thus identifying an effective Debye length for the medium that is greater than 
$\kappa ^{-1}$ (for the bath electrolyte). In distinct contrast are the weakly charged counterparts (in high-ionic-strength electrolyte) in figure 7(b), for which the dispersions are centred on 
$\omega a^2 / D_1 \approx 6$, thus identifying the cavity size 
$a = L / 4$ as the characteristic diffusion length. Note that designating these networks as being highly or weakly charged is from the perspective of the scaled charge density, since the dimensional charge densities are the same (
$\rho _f^0 / e = - 200$ mM).
Figure 7. Imaginary part of the scaled current density vs scaled angular frequency (note the different relaxations with respect to frequency) for rigid cavity-filled polyelectrolytes: 
$a = 2.5$ (blue), 
$5$ (red), 
$10$ nm (yellow) with 
$L / a = 4$. (a) ‘Highly charged’ network, low-ionic-strength electrolyte: 
$\rho _f^0 / (2 I e) = -100$, 
$I = 1$ mM. (b) ‘Weakly charged’ network, high-ionic-strength electrolyte: 
$\rho _f^0 / (2 I e) = -1$, 
$I = 100$ mM. Other parameters: 
$\ell = 1$ nm and 
$P = 0$.
The responses of the dielectric permittivity spectra to network elasticity with three 
$p$-wave moduli 
$M$ are shown in figure 8(b), here with 
$a = 5$ nm. Note that the real part of the current density in figure 8(a) is similar to its rigid-network counterpart in figure 6(a) (red), and is practically independent of 
$M$. However, the real part of the dielectric permittivity now reflects, in addition to the foregoing diffusion relaxation time scale, slower relaxations that depend on 
$M$. With an assumption that 
$a \sim L$, a poroelastic draining time (e.g. balancing the 
$\sim M v / a^2$ elastic and 
$\sim \eta u / \ell ^2$ hydrodynamic friction forces with 
$u \sim v / \tau _d$) is
\begin{equation} \tau_d = (\eta / M)(a / \ell)^2, \end{equation}
where 
$\tau _e = \eta / M$ is the viscoelastic time. For a cavity size 
$a = 5$ nm, 
$\ell = 1$ nm, and 
$M / (2 I k_B T) = 10$ giving 
$M \approx 50$ kPa with 
$I = 1$ mM, we find 
$\tau _d^{-1} \kappa ^{-2} / D_1 \approx 0.14$ and 
$\tau _e^{-1} \kappa ^{-2} / D_1 \approx 3.5$, which enclose the relaxation centred at 
$\sim \omega \kappa ^{-2} / D_1 = 1$ in figure 8(b) (red). Whereas the Nyquist plots for rigid networks (figures 6c and 6d) are semicircular, the poroelastic and viscoelastic relaxations distort these in their respective low-frequency regions.
Figure 8. Real part of the scaled current density (a) and dielectric permittivity (b) vs scaled angular frequency, and Nyquist representations of the scaled admittance (c) and scaled impedance (d) for compressible cavity-filled polyelectrolytes: 
$M / (2 I k_B T) = 3$ (blue), 
$10$ (red), 
$30$ (yellow). Other parameters: 
$a = 5$ nm, 
$\ell / a = 1/5$, 
$L / a = 5$, 
$\delta = \ell / a$, 
$I = 1$ mM (
$\kappa a \approx 0.51$, 
${Pe}_1 \approx 0.093$, 
${Pe}_2 \approx 0.023$) and 
$\rho _f^0 / (2 I e) = -100$. Solid and dashed lines correspond to 
$\langle u \rangle = 0$ and 
$P = 0$, respectively. Frequency in (a–c) increases from left to right [right to left in (d)], spanning the range 
$f = 10^2$–
$10^6$ kHz.
The factor by which the scaled dielectric permittivity in figure 8(b) may be converted into a relative dielectric permittivity is 
$\epsilon _c / {Pe}_h \approx 200$. Accordingly, the low-frequency plateaus for 
$P = 0$ (dashed lines) with magnitude 
$\sim 0.6$ correspond to a relative dielectric permittivity 
$\epsilon _r \sim 120$, which is somewhat larger than the value 
$\epsilon _s \approx 80$ for the solvent (water). A relative dielectric permittivity of 
$120$ furnishes an areal capacitance 
$C = \epsilon _r \epsilon _0 / L \approx 430\ \mathrm {\mu }\textrm {F}\ \textrm {cm}^{-2}$. At lower frequencies, however, 
${Pe}_h \epsilon '$ for
$\langle u \rangle = 0$ (solid lines) diverges, as highlighted by the double logarithmic plots in figure 9.
Figure 9. Same as figure 8, but with double logarithmic axes to highlight (i) the low-frequency power-law divergences and plateaus, and (ii) comparison of the solid blue line (softest sample) in panel (b) with its 
$\epsilon '$ counterpart for Nafion sulfonate as measured by Mauritz & Fu (Reference Mauritz and Fu1988, figure 4a).
For the softest network in figure 8(b) (blue, solid line) 
${Pe}_h \epsilon '$ diverges in the intermediate frequency range as 
$\omega ^{-n}$ with 
$0< n \lesssim 1$. Similar divergences are reported from experiments conducted on Nafion sulfonate films equilibrated in (concentrated 
$> 7$ M) NaOH electrolytes by Mauritz & Fu (Reference Mauritz and Fu1988). In their figure 4(a), for example, the relative dielectric permittivities increase from 
$\epsilon _r \sim 10$ at 
$f \sim 1$ MHz to 
$\epsilon _r \sim 10^5$ at 
$f \sim 10$ Hz, with power-law exponents 
$n \lesssim 1.2$. Such large dielectric permittivities have been a source of long-standing controversy pertaining to the polarisation mechanism, a cross-over from a.c. to d.c. conductivity, and roles of bulk and electrode polarisation (Matos Reference Matos2020). Note that 
${Pe}_h \epsilon ''$ (not shown) diverges as 
$\omega ^{-1}$ (without a plateau, independent of 
$M$), as do the counterparts for the Nafion sulfonate films reported by Mauritz & Fu (Reference Mauritz and Fu1988) and, more recently, by Matos (Reference Matos2020, figure 8g,h).
The present model clearly identifies network compliance as a necessary factor that couples with microscale charge polarisation to shape the bulk dielectric permittivity. In some sense, this unifies the ‘segmental motion’ and ‘interfacial polarisation’ mechanisms discussed in this literature for several decades (Matos Reference Matos2020). While the dielectric permittivity spectra calculated here bear resemblance to those measured by Mauritz & Fu (Reference Mauritz and Fu1988) and reported in more recent studies (Matos Reference Matos2020), it should be cautioned that Nafion membranes have a connected pore network arising from a hydrophobic condensed phase from which mobile ions are excluded (Kusoglu et al. Reference Kusoglu, Vezzu, Hegde, Nawn, Motz, Sarode and Herring2020). The present model suggests that it is the heterogeneity of the network compliance and fixed charge that principally shape the diverging dielectric permittivity, not the specific pore geometry or microscale contrast in dielectric constant.
At low enough frequencies, figure 9(a) shows how 
${Pe}_h \epsilon '$ increases as 
$\omega ^{-2}$ in an intermediate frequency range, eventually reaching plateaus where 
$\epsilon '$ may be considered a dielectric constant. The 
$\omega ^{-2}$ scaling is the same as captured by the analytical theory for electrode polarisation in simple electrolytes by Hollingsworth & Saville (Reference Hollingsworth and Saville2003). Here, however, the polarisation arises from the bulk medium rather than the electrodes.
Further insight into the distinctive differences between the low-frequency behaviour of 
$\epsilon '$ for 
$P = 0$ (dashed lines) and 
$\langle u \rangle = 0$ (solid lines) can be gleaned from a model of a single charge 
$q$ that is tethered to a linear spring (spring constant 
$k$, friction coefficient 
$\gamma$), subjected to a flow with velocity 
$\langle u \rangle = - (\rho _f^0 E + P) \ell ^2 / \eta$. The contribution of this charge (tether number density 
$n$) to the current density is
\begin{equation} -{\rm i} \omega \hat{v} n q ={-} n q ({\rm i} \omega \gamma / k) \frac{1 + {\rm i} \omega \gamma / k}{1 +(\omega \gamma/ k)^2} (q / \gamma + \langle u \rangle / E) E. \end{equation}
Thus, when the electroosmotic flow is arrested by a pressure gradient (
$\langle u \rangle = 0$), there is a positive low-frequency plateau contributing to the real part of the dielectric permittivity 
$\epsilon ' = n q^2 / (k \epsilon _0)$. Similarly to the full model, this is inversely related to 
$k$, and is accompanied by an 
$\sim \omega ^{-2}$ scaling when 
$\omega \gtrsim k / \gamma$. Moreover, in the absence of a pressure gradient (
$P = 0$), (3.13) identifies the possibility of 
$\epsilon '$ being small (or even negative) due to the opposite signs of the electrophoretic mobility of the tethered charge 
$q / \gamma$ and the electroosmotic mobility of the fluid 
$- \rho _f^0 \ell ^2 / \eta$. How the model parameters (e.g. the ratio 
$k / \gamma$) should be applied to real polyelectrolyte networks is beyond the present approximation that 
$M(x)$ is proportional to 
$\ell ^{-2}(x)$. Moreover, real networks have distributions of relaxation times that reflect a variety of chain-relaxation mechanisms, as captured, for example, by the empirical generalisation of a Debye relaxation by Havriliak & Negami (Reference Havriliak and Negami1967), e.g. as fitted to 
$\epsilon ''$ spectra for Nafion membranes by Matos (Reference Matos2020, figure 7). In the present model, such a spectrum is captured, in part, by the elastic modulus decaying (from its bulk value 
$M$) to zero through the cavity interface.
Note that the transitions to the low-frequency plateaus in figure 9 occur at frequencies well below those for the foregoing poroelastic and ion-diffusion relaxations. This suggests relaxations occurring on a larger length scale, e.g. 
$L > a$, and/or arising from relaxations in regions, such as the cavity interface, where 
$M$ is locally small, as captured in the low-dimensional model furnishing (3.13) when 
$k$ is small. One may then question why 
$\ell ^2$ appearing in (3.12) should not also be large in such regions, thus limiting 
$\tau _d$ when 
$M$ is locally small. It must be noted, however, that 
$\ell ^2$ in (3.12) should be restricted to the small value that is representative of the bulk. This is because fluid draining is limited by those parts of the microstructure that present the greatest resistance to flow.
In summary, a large low-frequency dielectric constant is favoured by highly compliant, charged microdomains, which are present in this model within the diffuse cavity interface. Note that calculations (not shown) undertaken with a less diffuse (sharper) cavity interface, achieved by setting 
$\delta = 0.2 \ell / a$, furnished dielectric permittivity spectra with much lower low-frequency plateaus (
${Pe}_h \epsilon ' \sim 10$–
$100$). Moreover, these were reached within a much higher frequency range (
$\omega \kappa ^{-2} / D_1 \sim 0.01$–
$0.1$). Interestingly, with 
$M / (2 I k_B T) = 3$, the low-frequency power-law divergence in figure 9(a) was suppressed completely, the plateau being reached via the 
$\omega ^{-1}$ scaling for 
$M / (2 I k_B T) = 3$ in figure 9(b) when 
$\omega \kappa ^{-2} / D_1 \sim 1$.
3.5. Lamellar networks
In this section, networks with uniform hydrodynamic permeability 
$\ell ^2$ and 
$p$-wave modulus 
$M$ are endowed with a spatially modulated (oscillatory) fixed-charge density
\begin{equation} \rho_f^0(x) = \rho_f^0 \{\cos^2{[{\rm \pi} (x / a -1/4)]} - \sin^2{[{\rm \pi} (x / a -1/4)]}\} = \rho_f^0 \sin{(2 {\rm \pi}x / a)}, \end{equation}
where 
$a$ and 
$\rho _f^0$ are now the period and amplitude of the fluctuating fixed-charge density, respectively. As shown by the example in figure 10, the fixed charge (green line) alternates between positive and negative values. Note that the sum of the periodic cationic and anionic number densities is uniform, consistent with a uniform permeability and 
$p$-wave modulus. This mimics polyelectrolyte networks synthesised using layer-by-layer deposition (Durstock & Rubner Reference Durstock and Rubner2001; Schönhoff Reference Schönhoff2003), assuming, for simplicity, that the layers (having equal thickness) interpenetrate, maintaining an approximately uniform segment and (effective) cross-linking densities (Decher Reference Decher1997). In contrast to the cavity-doped media in the previous section, with zero net fixed charge, uniform hydrodynamic permeability and an incompressible fluid, the network dynamics (compression and rarefaction) occur with zero electroosmotic flow (
$\langle u \rangle = 0$) and zero pressure gradient (
$P = 0$). Accordingly, all the computations in this section were undertaken with 
$E = 1$ and 
$P = 0$, furnishing 
$\langle u \rangle \approx 0$.
Figure 10. Scaled equilibrium electrostatic potential (
$\psi ^0$, blue), mobile-ion concentrations (
$c_i^0$, red and yellow), mobile-charge density (
$\rho _m^0$, violet), fixed-charge density (
$\rho _f^0$, green) and total charge density (
$\rho _m^0 + \rho _f^0$, black) in lamellar polyelectrolytes [fixed-charge density according to (3.14)]: 
$L = a = 5$ nm, 
$I = 100$ mM (
$\kappa a \approx 5.1$) and 
$\rho _f^0 /(2 I e) = -1$.
The results in figure 11, which are for a highly charged network (
$\rho _f^0 /(2 I e) = -100$, zero net charge) in equilibrium with a low-ionic-strength electrolyte (
$I = 1$ mM), explore how 
$M$ influences the conductivity (panels a,b), permittivity spectra (panels c,d) and Nyquist representations of the scaled complex conductivity (panel e, admittance) and its reciprocal (panel f, impedance). The stiffest network with 
$M / (2 I k_B T) = 10\ 000$ (green) produces spectra that are similar to those above for rigid, cavity-filled anionic networks. This reflects the electrical polarisation and relaxation dynamics of the mobile-ion populations. However, with a sufficiently compliant network, low-frequency dispersions reflect poroelastic network dynamics. Here, these impart a pair of semicircular arcs to the Nyquist impedance for compliant networks, morphing into a single arc for increasingly rigid networks. With the prevailing 
$a = 10$ nm and 
$\ell = 0.5$ nm and 
$M / (2 I k_B T) = 300$ (yellow) giving 
$M \approx 1.5$ MPa with 
$I = 1$ mM, for example, we find 
$\tau _d^{-1} \kappa ^{-2} / D_1 \approx 0.3$ and 
$\tau _e^{-1} \kappa ^{-2} / D_1 \approx 100$. These span the frequency range in which the relaxation (varying with 
$M$) in figure 11 falls, shifting to the left and right with smaller and larger values of 
$M$, respectively. Note that the magnitude and shape of such plots is sensitive to the period 
$a$ and hydrodynamic permeability 
$\ell ^2$.
Figure 11. (a) Real and (b) imaginary parts of the scaled current density vs scaled angular frequency (
$\,f = 10^2$–
$10^6$ kHz) for compressible, lamellar polyelectrolytes [fixed-charge density according to (3.14)]: 
$M / (2 I k_B T) = 30$ (blue), 
$100$ (red), 
$300$ (yellow), 
$1000$ (violet) and 10 000 (green). Other parameters: 
$I = 1$ mM, 
$\rho _f^0 /(2 I e) = -100$, 
$L = a = 10$ nm, 
$\ell = 0.5$ nm and 
$P = 0$. (c) Imaginary and (d) real parts of the scaled dielectric permittivity vs scaled angular frequency. (e,f) Nyquist plots of the scaled current density (as a measure of the complex admittance) and its reciprocal (as a measure of the complex impedance). Frequency in (e,f) increases from left to right, and from right to left, respectively, spanning the range 
$f = 1$–
$10^8$ kHz.
Figure 12 details the inner region of the Nyquist plots of admittance and impedance for networks that are the same as in figure 11, but with a larger period 
$a = 20$ nm. Nyquist representations of the impedance for polyelectrolyte multilayers reported by Durstock & Rubner (Reference Durstock and Rubner2001), albeit with systematic variation in the temperature, bear similar resemblance to those in figure 12(b). From their experiments, a systematic increase in temperature furnished a qualitatively similar effect as decreasing the network stiffness, as might be expected by a thermal softening.
Figure 12. Nyquist plots (in the high-frequency regime) of the scaled current density (a, as a measure of the complex admittance) and its reciprocal (b, as a measure of the complex impedance) for compressible, lamellar polyelectrolytes [fixed-charge density according to (3.14)]: 
$M / (2 I k_B T) = 30$ (blue), 
$100$ (red), 
$300$ (yellow), 
$1000$ (violet) and 10 000 (green). Other parameters: 
$I = 1$ mM, 
$\rho _f^0 /(2 I e) = -100$, 
$L = a = 20$ nm, 
$\ell = 0.5$ nm and 
$P = 0$.
Qualitatively similar classes of impedance spectra have also been reported from neural implants (Williams et al. Reference Williams, Hippensteel, Dilgen, Shain and Kipke2007) and deep brain stimulation (Lempka et al. Reference Lempka, Miocinovic, Johnson, Vitek and McIntyre2009). The lower frequencies at which relaxations occur in these studies may reflect larger geometrical features of the microstructure than for the model in figure 12. For example, Williams et al. (Reference Williams, Hippensteel, Dilgen, Shain and Kipke2007) demonstrated temporal transitions in the shape of Nyquist plots of the impedance as arising from the evolving structure of neural tissue, associating electrical circuit elements with cells, extracellular spaces and electrodes. Their Nyquist plots resemble parts of those in figure 12(b). The Nyquist plot of the impedance reported by Lempka et al. (Reference Lempka, Miocinovic, Johnson, Vitek and McIntyre2009) features a semicircle and line, as captured by the blue and red lines in figure 12(b). These authors attributed the semicircular arcs to neural cells accumulating around the electrode, captured empirically using a circuit-element model (with a constant-phase element). Here, such features may also be ascribed to changing characteristics of the microstructure, such as length scales, permeability and stiffness.
Nyquist plots of admittance and impedance in figure 13 highlight the significance of hydrodynamic coupling and elastic compliance. Here the hydrodynamic coupling has been enhanced by a smaller Brinkman length (
$\ell = 0.2$ nm with 
$L = a = 10$ nm). The Nyquist impedance bears resemblance to those reported by Farhat & Hammond (Reference Farhat and Hammond2005, figure 4) for lamellar polyelectrolyte fuel-cell membranes under dry and increasing humidity (frequencies 
$1$–
$10^7$ Hz). The authors modelled these data using two Randles cells (parallel capacitor and resistor) in series: one was attributed polarisation of the membrane pores and the other to the membrane–electrode interface, the latter including a Warberg constant-phase element (with capacitance and fractional exponent). Thus, six fitting parameters were adopted for each level of humidity. However, according to the present model, the two semicircular arcs can be firmly attributed ion and poroelastic relaxations of the membrane, since electrode polarisation is absent. According to the present model, electrode polarisation must be attributed to the experimental low-frequency slanted line.
Figure 13. Nyquist plots (
$\,f = 10^2$–
$10^{8}$ kHz) of the scaled current density (a, as a measure of the complex admittance) and its reciprocal (b, as a measure of the complex impedance) for compressible, lamellar polyelectrolytes [fixed-charge density according to (3.14)]: 
$M / (2 I k_B T) = 100$ (blue), 
$300$ (red), 
$1000$ (yellow), 
$3000$ (violet) and 10 000 (green). Other parameters: 
$I = 1$ mM, 
$\rho _f^0 /(2 I e) = -100$, 
$L = a = 10$ nm, 
$\ell = 0.2$ nm and 
$P = 0$.
Figure 14(b) shows the Nyquist impedance with a smaller 
$\ell = 0.2$ nm than in figure 13(b), all other model parameters are the same. Increasing the hydrodynamic friction attenuates the skeleton/fixed-charge contribution to the current, thus increasing the real part of the impedance where the semicircular arcs meet, leaving the low- and high-frequency impedances unchanged, since these reflect mobile-ion concentrations and their low- and high-frequency mobilities. The changes in these limits, as registered by the experiments of Farhat & Hammond (Reference Farhat and Hammond2005) by varying humidity, may be captured by adjusting the counter-ion concentration and/or mobility. In figure 14(b), the fixed-charge concentration as been halved (
$\rho _f^0 /(2 I e) = -50$) and the bilayer thickness doubled (
$a = 20$ nm) to mimic swelling, albeit keeping 
$\ell$ and 
$M$ fixed. Interestingly, the high-frequency impedance decreases and the low-frequency increases.
Figure 14. Same Nyquist impedance plots as in figure 13(b), but with adjusted parameters: (a) 
$\ell = 0.15$ nm, 
$L = a = 10$ nm (
$\,f = 10^2$–
$10^{8}$ kHz), 
$M / (2 I k_B T) = 100$ (blue), 
$300$ (red), 
$1000$ (yellow), 
$3000$ (violet) and 10 000 (green); and (b) 
$\ell = 0.15$ nm, 
$L = a = 20$ nm, 
$\rho _f^0 /(2 I e) = -50$ (
$\,f = 1$–
$10^{8}$ kHz), 
$M / (2 I k_B T) = 30$ (blue), 
$100$ (red), 
$300$ (yellow), 
$1000$ (violet) and 
$3000$ (green).
The low- and high-frequency conductivities (reciprocal to the impedance) are naively expected to increase in proportion to the fixed charge density, since this increases the mobile counter-ion concentration. While this is superficially the case at low frequency, the high-frequency limit is, perhaps surprisingly, to the contrary, since this limit it not subject to microscale polarisation. The sequence of equilibrium ion concentrations and charge densities shown in figure 15 reveal that increasing the bilayer thickness from 
$10$ to 
$80$ nm transforms the microstructure from one that is in a regime where 
$\kappa a \sim 1$ with the interpenetrating anionic and cationic polyelectrolyte layers electrostatically neutralising each other, releasing a large portion of their counter-ions to the electrolyte bath. Increasing the bilayer thickness brings the layers into a regime where 
$\kappa a \gtrsim 1$, for which the layers are increasingly (screened) by their respective Na
$^+$ or OH
$^-$ counter-ions.
Figure 15. Equilibrium profiles [electrostatic potential (
$\psi ^0$, blue), mobile-ion concentrations (
$c_i^0$, red and yellow), mobile-charge density (
$\rho _m^0$, violet), fixed-charge density (
$\rho _f^0$, green) and total charge density (
$\rho _m^0 + \rho _f^0$, black)] for lamellar microstructures with systematically varying fixed-charge density (amplitude) and bilayer thickness: (a) 
$L = a = 10$ nm (
$\kappa a \approx 1.03$), 
$\rho _f^0 /(2 I e) = -100$; (b) 
$L = a = 20$ nm (
$\kappa a \approx 2.06$), 
$\rho _f^0 /(2 I e) = -50$; (c) 
$L = a = 40$ nm (
$\kappa a \approx 4.12$) , 
$\rho _f^0 /(2 I e) = -25$; (d) 
$L = a = 80$ nm (
$\kappa a \approx 8.23$), 
$\rho _f^0 /(2 I e) = -12.5$. These demonstrate thin layers (with respect to the Debye length) releasing mobile counter-ions to the bath, presenting a lower concentration of mobile counter-ions than thicker layers with lower fixed charge density. Other parameters: 
$I = 1$ mM (
$\kappa ^{-1} \approx 9.72$ nm).
The low- and high-frequency conductivities accompanying the microstructures in figure 15 are listed in table 2. As expected, based on the foregoing interpretation of the equilibrium states, the high-frequency conductivities increase with decreasing fixed-charge density until the layers are sufficiently thick (
$\kappa a \gg 1$) to furnish a conductivity that is commensurate with the mobile ion concentration. A solution of the Poisson–Boltzmann equation with a sinusoidal fixed charge density (
$z$–
$z$ electrolyte, equal ion mobilities, see Appendix A) for 
$\kappa a \rightarrow \infty$ is (all dimensional variables)
\begin{equation} \frac{\sigma}{\sigma_\infty} = \left\langle \cosh{\left\{ |z| \mbox{asinh}\left[\frac{\rho_f^0}{2 I e|z|} \sin{(2 {\rm \pi}x / a)}\right] \right\}} \right\rangle, \end{equation}
where 
$\langle \cdot \rangle$ is the average over the period of 
$\sin {(2 {\rm \pi}x / a)}$. For the microstructure with 
$\kappa a \approx 8.23$ and 
$-\rho _f^0 /(2 I e) = 12.5$, this predicts 
$\sigma / \sigma _\infty \approx 8.07$, thus overestimating the high-frequency model calculation (
$\sigma / \sigma _\infty \approx 6.52$). The discrepancy seems to reflect the finite value of 
$\kappa a$. Indeed, solving the linearised Poisson–Boltzmann (LPB) equation for arbitrary 
$\kappa a$ (albeit for 
$|\rho _f| / (2 I e) \ll 1$, Appendix A) furnishes (all dimensional variables)
\begin{equation} \frac{\sigma}{\sigma_\infty} = 1 + \frac{1}{4} |z|^2 \left[ \frac{\rho_f^0}{(2 {\rm \pi})^2 k_B T \epsilon \epsilon_0 / (e a^2)+ 2 I |z|^2 e } \right]^2 + \cdots . \end{equation}
As highlighted by a comparison with its nonlinear counterpart in figure 17, the Debye–Hückel regime furnishes 
$\sigma / \sigma _\infty$ with 
$\kappa a = 8.23$ that are 
$\approx 2.4$ times larger than in the limit 
$\kappa a \rightarrow \infty$.
Table 2. Scaled conductivity for the microstructures in figure 15. These computations undertaken for a NaOH electrolyte with 
$I = 1$ mM, 
$\ell = 0.15$ nm and 
$M / (2 I k_B T) = 10\ 000$ (rigid skeleton).
Next, the low-frequency conductivities in table 2 are noted to be an order of magnitude smaller that their high-frequency counterparts. These vary in proportion to the fixed-charge density when 
$\kappa a \sim 1$ until 
$\kappa a \gtrsim 1$. As already seen for cavity-doped microstructures in the previous subsection, low-frequency (including the steady-state limit) imparts strong microscale polarisation due to counter-ion electromigration and diffusion, generally increasing with the magnitude of the equilibrium electrostatic potential 
$\psi ^0$. As indicated by table 2, the most highly charged interpenetrating polyelectrolytes furnish the lowest equilibrium electrostatic potential (max. 
$|\psi ^0| \approx 64$ mV). This reflects the small value of 
$\kappa a$ and accompanying co-neutralisation.
The network displacement and velocity for a lamellar microstructure with the real parts of the gradients of electrostatic potential, mobile concentrations, etc. are shown in figure 16. The odd symmetries 
$\hat {v} = -\hat {v}(-x)$ (figure 16a) and 
$\hat {\rho }_m(x) = -\hat {\rho }_m(-x)$ (figure 16d) are consistent with an accumulation of positive and negative fixed charge to the immediate left and right of 
$x = 0$. Accordingly, the region between the oppositely charged (and therefore oppositely migrating) domains (centred on 
$x = 0$) is stationary (figure 16a and b) with zero perturbation in the fixed and mobile charge at 
$x = 0$ (figure 16d). Such nodes correspond to even symmetry 
$\hat {\psi }_x(x) =\hat {\psi }_x(-x)$ (figure 16c) identifying deformation-induced charge dipoles. In this example, there is an ostensibly large (positive) perturbation to the electric field (
$- \hat {\psi }_x$) at 
$x = 0$ (figure 16c). Note that the same electric-field perturbation (
$- \hat {\psi }_x$) is induced at the other nodes of the fixed-charge density perturbation (
$x /a = \pm 0.5$). Here, however, the fixed-charge perturbation reflects local depletions of the fixed charge populations, since the neighbouring fixed charge-density populations migrate away from each other. These dynamics explain the perhaps unexpected symmetries 
$\hat {c}_{1,x}(x) = -\hat {c}_{2,x}(-x)$ (figure 16c) and 
$\hat {c}_1(x) = \hat {c}_2(-x)$ (figure 16d).
Figure 16. (a,b) Real (solid) and imaginary (dashed) parts of the scaled network displacement (a) and velocity (b) for soft, compressible, lamellar polyelectrolytes [fixed-charged according to (3.14)]: 
$M / (2 I k_B T) = 10$, 
$I = 1$ mM (
$\kappa a \approx 0.51$), 
$\rho _f^0 /(2 I e) = -100$. Other parameters: 
$f = 100$ kHz, 
$L = a = 5$ nm, 
$\ell = 0.5$ nm and 
$P = 0$. (c) Real parts of the (scaled) perturbation gradients of electrostatic potential (
$\hat {\psi }_x$, blue), cation concentration (
$\hat {c}_{1,x}$, red) and anion concentration (
$\hat {c}_{2,x}$, yellow). (d) Real parts of the perturbed mobile-charge density (
$\hat {\rho }_m$, blue), cation concentration (
$\hat {c}_{1}$, red) and anion concentration (
$\hat {c}_{2}$, yellow). All the hatted quantities are dimensionless, calculated with dimensionless 
$E = 1$.
4. Summary
Hydrodynamic and elastic physics have been integrated into an ion transport model for structured polyelectrolytes. The equations were solved numerically on one-dimensional, periodic unit cells, for anionic cavity-doped microstructures and lamellar media with alternating layers of positive and negative fixed charge, but uniform hydrodynamic permeability and elastic moduli.
As a benchmark, the dielectric permittivity spectra of rigid polyelectrolytes were established as presenting a well-defined, low-frequency dielectric constant. The single relaxation time is an ion-diffusion time that is based on a microstructural length, which is the cavity scale 
$a$ when the ionic strength in the polyelectrolyte is dominated by added electrolyte; otherwise, it becomes the Debye length based on the counter-ion concentration.
For compliant polyelectrolytes, a linearly elastic, fluid-permeable skeleton produced additional relaxations with poroelastic and viscoelastic time scales. These furnished a significant increase in the real part of the dielectric permittivity with respect to frequency at low frequencies, breaking the symmetry of otherwise semicircular Nyquist plots of the impedance.
For microstructures comprising alternating layers of positive and negative charge, but uniform hydrodynamic friction and network elasticity, much more intricate impedance spectra emerged. Examples were presented demonstrating a notable sensitivity to the coupling of electrical, elastic and hydrodynamic drag forces. These produced a variety of Nyquist impedance and admittance plots, featuring pairs of semicircular arcs, merging into one when increasing the network stiffness. Whereas empirical (equivalent circuit-based) modelling in the literature has attributed secondary arcs to electrode polarisation, the present theory provides a new interpretation.
Also noteworthy was the ability of thin, oppositely charged, interpenetrating polyelectrolyte layers to electrostatically neutralise each other, thus releasing a significant portion of mobile counter-ions, decreasing (decreasing) the conductivity (impedance). Between the low- and high-frequency limits, hydrodynamic coupling of the skeleton to the electrolyte, and polymer network stiffness, shape the Nyquist impedance (and admittance), as exemplified in experimental literature by changes in temperature and humidity.
In general, the real part of the scaled current density/conductivity transitions from a low-frequency plateau that reflects electrostatic polarisation of the microstructure, significantly slowing ion transport, manifesting experimentally as hindered ion mobility or diminished ion concentration. With increasing frequency, the low-frequency conductivity transits to a high-frequency plateau that, perhaps paradoxically, reflects the d.c. mobilities (and composition) of ions within the polyelectrolyte network. As already highlighted, sufficiently thin anionic and cationic bilayered microstructures can release their mobile counter-ions to an electrolyte bath, diminishing the expected counter-ion concentration and effective conductivity.
The electrostatic interactions that polarise the (one-dimensional) microstructures in this study are very strong. It will therefore be pertinent, and much more computationally challenging, to apply the model to more realistic two- and three-dimensional microstructures. Such calculations will be appropriate for quantifying how the cavity volume fraction, for example, affects impedance spectra: an aspect of the present one-dimensional model that was largely ignored. Predictions of averaged ion diffusion coefficients, which may be necessary for modelling low-frequency dynamics at the macroscale may be necessary to disentangle bulk- and electrode-proximate ion transport.
Funding
Financial support from an NSERC Discovery Grant is gratefully acknowledged.
Declaration of interests
The author reports no conflict of interest.
Appendix A. Analysis of the equilibrium state for the lamellar microstructure
For the lamellar microstructure in the main text, here simplified for a 
$z$–
$z$ electrolyte, the scaled Poisson–Boltzmann equation (
$\psi ^0$ scaled with 
$k_B T / e$, 
$x$ scaled with 
$a$ and 
$\rho _f^0$ scaled with 
$2 I e$) is
\begin{equation} \psi^0_{xx} = (\kappa a)^2 [|z| \sinh{(|z| \psi^0)} - \rho_f^0 \sin{(2 {\rm \pi}x)}]. \end{equation}
An approximate solution for 
$\kappa a \gg 1$ is the electrostatic potential for local Donnan equilibrium with the electrolyte bath:
\begin{equation} \psi^0(x) = \mbox{asinh}[(\rho_f^0 / |z|) \sin{(2 {\rm \pi}x)}]. \end{equation}The accompanying concentrations of mobile ions (assuming equal ion mobilities) furnish an averaged conductivity (scaled with the bulk electrolyte conductivity, high-frequency limit or d.c. without polarisation effects)
\begin{align} \frac{\sigma}{\sigma_\infty} &= \langle \cosh{[|z| \psi^0(x)]}\rangle \nonumber\\ &= \langle \cosh{\{ |z| \mbox{asinh}[(\rho_f^0 / |z|) \sin{(2 {\rm \pi}x)}] \}} \rangle = \frac{2}{\rm \pi} E(-{\rho_f^0}^2) {\quad \mbox{for } |z| = 1}, \end{align}
where 
$\langle \cdot \rangle = \int _{-1/2}^{1/2} \ \cdot\ \mbox {d}\kern 0.06em x$ and 
$E(x)$ is the complete elliptic integral of the second kind. With 
$\rho _f^0 = -12.5$ (table 2 in the main text with 
$\kappa a \approx 8.23$, 
$|z| = 1$), this furnishes 
$\sigma / \sigma _\infty \approx 8.07$, whereas the full model furnishes 
$6.52$ (at high frequency).
Now for arbitrary 
$\kappa a$, but under the Debye–Hückel approximation (
$|\psi ^0| \ll 1$), the LPB equation
\begin{equation} \psi^0_{xx} = (\kappa a)^2 [|z|^2 \psi^0 - \rho_f^0 \sin{(2 {\rm \pi}x)}], \end{equation}
has a periodic solution (on 
$|x| \le 1/2$)
\begin{equation} \psi^0(x) = \frac{(\kappa a)^2 \rho_f^0}{(2{\rm \pi})^2 + (\kappa a)^2 |z|^2} \sin{(2 {\rm \pi}x)}. \end{equation}The scaled concentrations of mobile ions are
\begin{equation} c^0_i(x) = \tfrac{1}{2} [1 \pm |z| \psi^0(x) + \tfrac{1}{2} |z|^2 {\psi^0(x)}^2 + \cdots], \end{equation}giving averaged conductivity
\begin{equation} \frac{\sigma}{\sigma_\infty} = 1 + \frac{1}{2} |z|^2 \langle {\psi^0}(x)^2 \rangle +\cdots = 1 + \frac{1}{4}|z|^2 \left[ \frac{(\kappa a)^2 \rho_f^0}{(2 {\rm \pi})^2 + (\kappa a)^2 |z|^2} \right]^2 +\cdots.\end{equation}
Figure 17 compares the solution of the nonlinear Poisson–Boltzmann equation (A3), albeit for 
$\kappa a \rightarrow \infty$, with its linearised counterpart (A7) for 
$\kappa a =10$ and 
$100$.
Figure 17. Reduced conductivity for a lamellar polyelectrolyte with sinusoidal fixed charge density (period 
$a$, amplitude 
$\rho _f^0$), as furnished by (A3) (from the nonlinear Poisson–Boltzmann equation for 
$\kappa a \rightarrow \infty$, solid line) and its linearised counterpart (A7) for 
$\kappa a = 10$ (dash-dotted) and 
$100$ (dashed).
Open access
