# Graphene-based enhancement of near-field radiative-heat-transfer
rectification

Simon Landrieux  Université Paris-Saclay, Institut d’Optique Graduate School,
CNRS, Laboratoire Charles Fabry, 91127, Palaiseau, France  Philippe Ben-
Abdallah  Université Paris-Saclay, Institut d’Optique Graduate School, CNRS,
Laboratoire Charles Fabry, 91127, Palaiseau, France  Riccardo Messina
<EMAIL_ADDRESS>Université Paris-Saclay, Institut
d’Optique Graduate School, CNRS, Laboratoire Charles Fabry, 91127, Palaiseau,
France

(August 28, 2024)

###### Abstract

We present a thermal device based on the near-field interaction between two
substrates made of a polar and a metal-insulator-transition material. As a
result of the temperature dependence of the optical properties, this device
acts as a thermal rectifier, implying a strong asymmetry in the heat flux when
reversing the two temperatures. By covering both substrates with a graphene
sheet we show a significant enhancement of rectification coefficient. By
investigating the flux spectral properties along with its distance dependence,
we prove that this enhancement is associated to a change in the power-law
dependence of heat flux with respect to the separation distance in the
electrostatic regime due to the presence of graphene sheets. Our results
highlight the promising role of graphene-based hybrid structures in the domain
of nanoscale thermal management.

Two bodies at different temperatures exchange energy under the form of thermal
radiation even when separated by vacuum. This radiative heat transfer (RHT) is
well described at large separation distances by Stefan-Boltzmann’s law,
setting an upper bound to this energy exchange, reached only by two
blackbodies. During the second half of the 20th century, the works of Rytov
[1], Polder and van Hove [2] and the development of fluctuational
electrodynamics showed that this limit does not apply in the near field, i.e.
for separation distances smaller than the thermal wavelength
$\lambda_{\text{th}}=\hbar c/k_{B}T$ (some microns at ambient temperature). In
this regime, the flux can overcome even by several orders of magnitude the
blackbody limit, especially in the presence of resonant surface modes of the
electromagnetic field, such as surface phonon-polaritons for polar materials
[3]. Since then, these theoretical predictions have been confirmed by a large
number of experiments, stimulating in turn further theoretical investigations
(see Refs. 4, 5, 6 and references therein).

These developments have paved the way to the exploration of the role played by
near-field RHT in several applications, including heat-assisted data recording
[7, 8], infrared spectroscopy [9, 10], energy-conversion techniques [11, 12,
13, 47] and thermotronics [15, 16], i.e. the design of thermal equivalent of
circuit elements. A significant attention has been devoted, both theoretically
[17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44] and experimentally [45, 46, 48, 47], to
the rectification of RHT in the near field, consisting in a configuration
implying a strong asymmetry of heat flux when exchanging the temperatures of
the two bodies. In a two-terminal system, this effect is due to a temperature
dependence of the optical properties of the materials, and has been
extensively studied both for planar and non-planar systems [22, 36, 38], from
cryogenic to ambient and high temperatures. This dependence is not required in
systems made of more than two objects [44], in which rectification can result
from purely many-body effects. The sensitivity of the rectification
coefficient to a variety of aspects, such as film thickness [27, 29, 39],
doping [18], dielectric coating [20], gratings [28, 32, 33, 42, 43] has been
explored. Concerning the choice of materials, while phase-transition materials
[19, 45, 24, 23, 26, 28, 46, 31, 32, 47, 42, 39, 43] such as vanadium dioxide
(VO2) have gathered a remarkable attention, also superconductors [30, 40] and
graphene-based systems [31, 35, 34] have been investigated. While it is often
difficult to compare different near-field rectification schemes due to
strongly different geometrical and physical features (see Ref. 36 for a
summary of recent results), the most recent works [36, 43, 41, 39, 42, 40]
propose scenarios with an efficiency, defined as
$\eta=\frac{\Phi_{\text{max}}-\Phi_{\text{min}}}{\Phi_{\text{max}}}$
($\Phi_{\text{max}}$ and $\Phi_{\text{min}}$ being the maximum and minimum
fluxes, respectively), going beyond 90%, and limited mainly by the intrinsic
losses of the materials and the mismatch between their surface resonances.

In this work we address the role played by graphene in modulating and
enhancing the efficiency of a near-field radiative thermal rectifier based on
two parallel substrates made of VO2 and silica (SiO2). Graphene has recently
proved to be promising for near-field RHT, since it supports a delocalized
surface plasmon [49] which can be used to improve the coupling between two
different materials and whose properties can be modulated by acting on the
graphene chemical potential. Based on this behavior and on the coupling
between graphene plasmons and surface modes of an underlying substrate [50],
the possibility of tuning the heat transfer and increasing the performances of
energy-conversion devices has been investigated [51, 52, 53, 54, 55, 56, 57].
In this Letter, we propose a system in which a graphene sheet is deposed both
on the VO2 and the SiO2 films. By optimizing for each distance the chemical
potentials of the two sheets, we show a modulation of the efficiency of the
rectification mechanism and an increase by a maximum of 14% for a vacuum gap
of 100 nm. By studying the dependence on the separation distance and the
graphene chemical potential along with the flux spectral features, we
highlight the crucial role played by graphene plasmons in the increased
coupling between the two films. The main difference between our results and
previous works [31, 35, 34] focusing on graphene is both the simultaneous use
of two graphene sheets and the optimization of the chemical potentials for any
separation distance, allowing us to maximize the coupling between dissimilar
materials exchanging heat in the near field. Moreover, we highlight a non-
monotonic behavior of the efficiency enchancement ascribed to a transition in
the power-law behavior of RHT. These results pave the way to the active
manipulation of heat flux at the nanoscale.

The radiative thermal diode (RTD) we consider, depicted in the inset of Fig.
1, consists of two semi-infinite substrates separated by a vacuum gap of width
$d$, one made of vanadium dioxyde (VO2) and the other made of fused silica
(SiO2), both covered with a monolayer graphene sheet at a chemical potential
$\mu_{\text{VO}_{2}}$ and $\mu_{\text{SiO}_{2}}$, respectively. For a given
temperature gradient $\Delta T=T_{\text{high}}-T_{\text{low}}$, $\Phi_{F}$ is
the heat flux received by VO2 at $T_{\text{low}}$ from SiO2 at
$T_{\text{high}}$, whereas $\Phi_{B}$ is the one received by SiO2 at
$T_{\text{low}}$ from VO2 at $T_{\text{high}}$. The rectification coefficient
$\eta=\frac{\Phi_{F}-\Phi_{B}}{\Phi_{F}}$ is used to assess the efficiency of
the system as a thermal diode.

Figure 1: Rectification efficiency $\eta$ for the thermal diode without
graphene (dashed black line) or with the graphene sheets tuned to the chemical
potentials that maximize $\eta$ (red solid line) as a function of the
separation distance $d$. The calculations were done for
$T_{\text{high/low}}=T_{c}\pm 100\,$K.

It is consequently necessary to have $\Phi_{F}\gg\Phi_{B}$ in order to have
$\eta\simeq 1$. This is achieved (with or without graphene) by exploiting the
first-order Mott transition of VO2 [58] from a low-temperature insulating
phase to a high-temperature metallic one close to room temperature
($T_{c}=340\,$K). Indeed, if $T_{c}\in[T_{\text{low}},T_{\text{high}}]$, VO2
is dielectric at $T_{\text{low}}$ and metallic at $T_{\text{high}}$,
considerably reducing the transmission coefficient in the latter case due to
the absence of a resonant mode in the infrared range. This results in a
$\Phi_{B}$ significantly lower than $\Phi_{F}$, hence in a high $\eta$.
However, the use of a second material different from VO2 to break the symmetry
of the system implies a mismatch between the frequencies of the surface modes
of the two materials. Due to the crucial role they play in the RHT at the
nano-scale, this has a negative effect on $\Phi_{F}$ and consequently on
$\eta$.

The way to enhance the rectification coefficient $\eta$ of the RTD we present
consists in increasing $\Phi_{F}$ by adjusting the chemical potentials
$\mu_{\text{VO}_{2}}$ and $\mu_{\text{SiO}_{2}}$ of the two graphene sheets to
tune the surface modes of the two sides of the system. An upper limit of
$\mu=0.5\,$eV has been set for the calculations to be coherent with
experimentally attainable values. We compute the heat fluxes
$\Phi_{\text{F},\text{B}}$ between the two bodies by means of the following
Landauer-like decomposition [59, 60]:

$\Phi=\int_{0}^{\infty}\frac{d\omega}{2\pi}\Delta\Theta(\omega)\sum_{p}\int\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}\mathcal{T}_{p}(\omega,\mathbf{k}),$
(1)

where
$\Delta\Theta(\omega)=\Theta(\omega,T_{\text{high}})-\Theta(\omega,T_{\text{low}})$
is the difference in mean energy of a harmonic oscillator
$\Theta(\omega,T)=\hbar\omega/[\exp(\hbar\omega/k_{B}T)-1]$ at the two
temperatures, $\mathbf{k}$ is the projection of the wave vector on the surface
of the semi-spaces and $\mathcal{T}_{p}$ is the Landauer transmission
coefficient for the polarization $p$ [transverse electric (TE) or magnetic
(TM)], taking values between 0 and 1, giving the proportion of the energy
discrepancy $\Delta\Theta$ between the modes of the two sides that is actually
transferred from the hot to the cold side. For two parallel semi-infinite
planar bodies separated by a vacuum gap of width $d$, it is given by

$\mathcal{T}_{p}(\omega,\mathbf{k})=\begin{cases}\displaystyle\frac{(1-|r_{1,p}|^{2})(1-|r_{2,p}|^{2})}{|1-r_{1,p}r_{2,p}e^{2ik_{z}d}|^{2}},&k<\omega/c\\\
\displaystyle\frac{4\
\text{Im}(r_{1,p})\text{Im}(r_{2,p})e^{-2\text{Im}(k_{z})d}}{|1-r_{1,p}r_{2,p}e^{-2\text{Im}(k_{z})d}|^{2}},&k>\omega/c\end{cases}$
(2)

where $k_{z}=\sqrt{(\omega/c)^{2}-\mathbf{k}^{2}}$ is the component of the
wave vector normal to the surfaces, $r_{i,p}$ is the reflection coefficient
for medium $m$ at polarization $p$, which implicitly depends on $\omega$ and
$k$. In the case of a graphene-covered surface, the reflection coefficients
read [54]:

$\begin{split}r_{i,\text{TE}}&=\displaystyle\frac{k_{z}-k_{z,m}-\mu_{0}\omega\sigma(\omega)}{k_{z}+k_{z,m}+\mu_{0}\omega\sigma(\omega)},\\\
r_{i,\text{TM}}&=\displaystyle\frac{\epsilon_{m}(\omega)k_{z}-k_{z,m}+\frac{1}{\epsilon_{0}\omega}\sigma(\omega)k_{z}k_{z,m}}{\epsilon_{m}(\omega)k_{z}+k_{z,m}+\frac{1}{\epsilon_{0}\omega}\sigma(\omega)k_{z}k_{z,m}},\end{split}$
(3)

where $\epsilon_{m}(\omega)$ is the dielectric permittivity of medium $m$,
$k_{z,m}=\sqrt{(\omega/c)^{2}\epsilon_{m}(\omega)-\mathbf{k}^{2}}$ is the
normal component of the wave vector inside medium $m$ and $\sigma(\omega)$ is
the conductivity of graphene. The latter is typically expressed as the sum of
an intraband (Drude) contribution $\sigma_{D}(\omega)$ and an interband
contribution $\sigma_{I}(\omega)$, which are given by the following
expressions [61]:

$\begin{split}&\sigma_{D}(\omega)=\frac{i}{\omega+i/\tau}\frac{2e^{2}k_{B}T}{\pi\hbar^{2}}\log\Big{[}2\cosh\Big{(}\frac{\mu}{2k_{B}T}\Big{)}\Big{]}\text{,}\\\
&\sigma_{I}(\omega)=\frac{e^{2}}{4\hbar}\Big{[}G\Big{(}\frac{\hbar\omega}{2}\Big{)}+i\frac{4\hbar\omega}{\pi}\int_{0}^{\infty}\frac{G(\xi)-G(\frac{\hbar\omega}{2})}{(\hbar\omega)^{2}-4\xi^{2}}d\xi\Big{]}\text{,}\end{split}$
(4)

where $\mu$ is the chemical potential of the graphene, $\tau$ is the
relaxation time, for which we take a value of $10^{-13}\,$s (see Refs. 52,
54), and
$G(x)={\sinh\frac{x}{k_{B}T}}/(\cosh\frac{x}{k_{B}T}+\cosh\frac{\mu}{k_{B}T})$.
As clear from Eq. (4), the chemical potential $\mu$ has an influence on the
conductivity of graphene, resulting in the possibility to tune, by controlling
(chemically or electrically) $\mu$, the Landauer coefficient $\mathcal{T}_{p}$
and thus the RHT between the two bodies.

In order to numerically highlight the enhancement of rectification efficiency
offered by the presence of graphene, we have calculated the forward and
backward fluxes for $T_{\text{high/low}}=T_{c}\pm 100\,$K in two
configurations, namely in the absence and in the presence of graphene. In the
latter scenario, for each distance $d$ between 10 nm and 10 $\mu$m, we have
found the values of $\mu_{\text{VO}_{2}}$ and $\mu_{\text{SiO}_{2}}$
maximizing the efficiency $\eta$. Concerning the permittivities, for SiO2 we
calculate it by interpolating experimental data [63], while for VO2, the
dielectric phase is approximated by an isotropic medium of dielectric constant
$\epsilon_{\text{VO}_{2},\text{d}}(\omega)$ computed by a classical oscillator
model and the metallic phase is described by a Drude-Lorentz model [62]. The
enhancement of the efficiency $\eta$ resulting from the presence of graphene
is illustrated in Fig. 1. We clearly observe, for values of $d$ lower than
$4\,\mu$m, an increase of the efficiency, becoming negligible for larger
distances. More specifically, the graphene-enhanced RTD has a maximal
rectification coefficient of $\eta=85\%$ for $d\simeq 100\ \text{nm}$, whereas
the one without graphene has $\eta=71\%$ for the same $d$. Moreover, the
adjunction of graphene enables to attain values of $\eta$ $6\%$ higher than
the absolute maximum ($\eta\simeq 78\%$) of the efficiency without graphene,
values which would be hard to achieve due to the small gap widths $d$. On the
other hand, at closer separation distance, we observe that the increase of
rectification coefficient becomes less pronounced. As discussed below, this
result can be ascribed to the behavior of graphene in the electrostatic
regime.

It is well known [54] that the plasmons of graphene play a relevant role for
relatively high values of the parallel wave-vector $k$ [54]. This explains the
impossibility to enhance the efficiency at distances above 1 $\mu$m. On the
contrary, in order to get further insight into the mechanism behind the
efficiency amplification taking place at smaller separation distances, we
focus on the configuration $d=100\,$nm, both without graphene and with
graphene at the chemical potentials $\mu_{\text{VO}_{2}}=0.4\,$eV and
$\mu_{\text{SiO}_{2}}=0.2\,$eV for VO2 and SiO2, respectively, corresponding
to the configuration maximizing $\eta$ at $d=100\,$nm.

Figure 2: Transmission coefficient
$\mathcal{T}_{\text{TM}}(\omega,\mathbf{k})$ in the $(k,\omega)$ plane for the
backward temperature gradient [(a), (b)] and for the forward temperature
gradient [(c), (d)]. The left column [(a), (c)] correponds to the RTD without
graphene, while the right column [(b), (d)] correponds to the RTD with
graphene at the optimal chemical potentials for $d\simeq 100\,$nm. The
horizontal dashed [(a)] and dotted (a) and (c)] lines correspond to the
surface resonances of SiO2 and VO2 (dielectric phase), respectively.

The participation of the surface resonant modes to the RHT in the two
configurations are represented in Fig. 2, where the transmission coefficient
in TM polarization $\mathcal{T}_{\text{TM}}(\omega,\mathbf{k})$, the one
contributing the most to the RHT, are represented in the $(k,\omega)$ plane.
In the first row (corresponding to the backward temperature gradient), the
main resonant modes observable are the ones of SiO2 at frequencies
$\omega_{\text{SiO}_{2},1}=0.9\times 10^{14}\,$rad/s and
$\omega_{\text{SiO}_{2},2}=2.3\times 10^{14}\,$rad/s, both without (a) and
with graphene (b). In the optimal configuration with graphene (b), a third
branch of low intensity for low $\omega$ is barely visible, and corresponds to
the surface plasmon of graphene alone. We conclude that the adjunction of
graphene has little effect on the transmission coefficients for the backward
temperature gradient due to the absence of resonant modes in the Planck window
for VO2 in its metallic phase. In the second row (forward temperature
gradient), a multitude of resonant modes are observable for the RTD without
graphene [panel (c)] in the $[0.8,1.1]\times 10^{14}\,$rad/s range, which can
be linked to resonant surface modes of VO2 in its dielectric phase. For the
optimal configuration with graphene [panel (d)], while the signature of the
modes of graphene alone can still be seen for low $\omega$, the most
noticeable effect of graphene is its coupling with the modes of the two
substrates, which results in the broadening of their frequency domain and in
an increased participation of modes with higher $k$.

These different effects on $\mathcal{T}_{\text{TM}}(\omega,\mathbf{k})$ have a
noticeable impact on the spectral features of the RHT, shown in Fig. 3. In the
first row, the mild effect of graphene on the RHT can be confirmed, as its
adjunction only modifies slightly the shape of the two peaks
$\omega_{\text{SiO}_{2},1}$ and $\omega_{\text{SiO}_{2},2}$ of SiO2 and adds a
small contribution for low $\omega$. Consequently, the effect of the graphene
on $\Phi_{\text{B}}$ is not substantial. Unlike $\Phi_{\text{B}}$, the
addition of graphene at the optimal chemical potentials has a considerable
effect on $\Phi_{\text{F}}$. As the broadening in $\omega$ of the resonant
modes observed with $\mathcal{T}_{\text{TM}}(\omega,\mathbf{k})$ in Fig. 2
suggested, the spectral flux in the forward scenario is significantly higher
on a wide range of $\omega$ (namely from $0.8\times 10^{14}$ rad/s to
$1.7\times 10^{14}$ rad/s). This broadening of the RHT spectrum leads to a
significative enhancement of $\Phi_{\text{F}}$ compared to the RTD without
graphene. This explains why the adjunction of graphene has such an effect on
$\eta$.

Figure 3: Spectral flux $\Phi_{\omega}$ for the backward temperature gradient
[(a), (b)] and for the forward temperature gradient [(c), (d)]. The left
column [(a), (c)] correponds to the RTD without graphene, while the right
column [(b), (d)] correponds to the RTD with graphene at the optimal chemical
potentials for $d\simeq 100$ nm.

It is now interesting to investigate the origin of the non-monotonic behavior
of the efficiency increase observed in Fig. 1. As a matter of fact, we notice
that for $d<30$ nm the optimal $\eta$ decreases when $d$ gets smaller. An
explanation to this phenomenon can be found by looking at Fig. 4, where the
evolution of $\Phi_{\text{F}}$ and $\Phi_{\text{B}}$ with respect to $d$ is
represented for three different configurations: in the absence of graphene
(dashed black lines), in the optimal configuration corresponding to $d=10$ nm
(red curves), and in the optimal configuration corresponding to $d=100$ nm
(blue curves). We stress that the optimal parameters for the two distances are
different, and more specifically we have $\mu_{\text{VO}_{2}}=0.16\,$eV and
$\mu_{\text{SiO}_{2}}=0.09\,$eV for $d=10$ nm, whereas
$\mu_{\text{VO}_{2}}=0.4\,$eV and $\mu_{\text{SiO}_{2}}=0.2\,$eV for $d=100$
nm.

Figure 4: Evolution of the heat fluxes $\Phi_{\text{B}}$ (three lower curves)
and $\Phi_{\text{F}}$ (three upper curves) for different configurations. The
black dashed lines correspond to the RTD without graphene, the red solid lines
to the optimal configuration for $d\simeq 10$ nm and the blue solid lines to
the optimal configuration for $d\simeq 100$ nm. The two gray dotted lines are
visual references for evolutions proportional to $1/d$ and $1/d^{2}$ as
indicated in the plot.

We start by observing that in the absence of graphene both $\Phi_{\text{F}}$
and $\Phi_{\text{B}}$ present the well-known $1/d^{2}$ diverging behavior [3]
of RHT in the near field. On the contrary, in the optimal configuration for
$d=100$ nm (blue curves), there is a clear transition below 35 nm towards a
$1/d$ behavior, due to the behavior of graphene in the electrostatic regime
[64]. At the same small distances, it is possible to guess the beginning of
the same transition in $\Phi_{\text{B}}$, which nevertheless still behaves
approximately as $1/d^{2}$ in the entire range of distances. As a consequence
of this transition in power-law behavior, we have an intermediate distance
range, containing $d=100\,$nm, where the curve corresponding the forward flux
$\Phi_{\text{F}}$ gets farther from the one associated with $\Phi_{\text{B}}$,
resulting in the observed increased efficiency. This analysis shows that this
transition in the power-law behavior plays indeed a key role in the
optimization of the RTD effect. If we keep focusing on the optimal
configuration for $d=100\,$nm, we also observe that for smaller distances the
transition to a $1/d$ behavior of $\Phi_{\text{F}}$ brings the two fluxes
closer to each other, producing then a reduction of efficiency.

If we know switch to the optimal configuration corresponding to $d=10\,$nm, we
clearly observe that a deviation from a $1/d^{2}$ behavior of both fluxes is
barely visible. We deduce that the characteristic distance at which this
power-law transition takes place depends, not surprisingly, on the values of
the two chemical potentials, and more specifically this distance decreases
with respect to both chemical potentials. As a consequence, the reduction of
the optimal values of $\mu_{\text{SiO}_{2}}$ and $\mu_{\text{VO}_{2}}$ for
$d=10\,$nm, compared to $d=100\,$nm, is expected in the sense that a
transition at smaller distances would be beneficial for the efficiency when
reducing $d$. Nevertheless, one should bear in mind that reducing these values
is also producing a detuning between the surface modes of the two materials.
While this compromise works at intermediate distances such as $d=100\,$nm, it
becomes less efficient when reducing the separation distance, explaining the
non-monotonic behavior of the optimized $\eta$ in Fig. 1 as a function of the
distance.

In conclusion, we have highlighted an enhancement of the efficiency of the
radiative thermal rectification in near field using tunable graphene sheets to
mediate the heat transfer between a polar and a phase-change material. We have
shown that this enhancement results from a better coupling between the surface
modes supported by the the two interacting solids thanks to the presence of
graphene sheets, and connected this effect to a transition in the power-law
dependence of heat flux from a $1/d^{2}$ to a $1/d$ behavior. Our findings
shed further light on the remarkable potential of graphene in the domain of
active thermal management.

## Author declarations

### Conflict of interest

The authors have no conflicts to disclose.

## Data availability

The data that support the findings of this study are available from the
corresponding authors upon reasonable request.

## References

  * [1] S. M. Rytov, Y. A. Kravtsov, and V. I. Tatarskii, _Principles of Statistical Radiophysics_ (Springer, New York, 1989).
  * [2] D. Polder and M. van Hove, Phys. Rev. B 4, 3303 (1971).
  * [3] K. Joulain, J.-P. Mulet, F. Marquier, R. Carminati, and J.-J. Greffet, Surf. Sci. Rep. 57, 59 (2005).
  * [4] J. C. Cuevas and F. J. García-Vidal, ACS Photonics 5, 3896 (2018).
  * [5] B. Song, A. Fiorino, E. Meyhofer, and P. Reddy, AIP Adv. 5, 053503 (2015).
  * [6] S.-A. Biehs, R. Messina, P. S. Venkataram, A. W. Rodriguez, J. C. Cuevas, and P. Ben-Abdallah, Rev. Mod. Phys. 93, 025009 (2021).
  * [7] W. A. Challener, C. Peng, A. V. Itagi et al., Nat. Photon. 3, 220 (2009).
  * [8] B. C. Stipe et al., Nat. Photon. 4, 484 (2010).
  * [9] Y. De Wilde, F. Formanek, R. Carminati, B. Gralak, P.-A. Lemoine, K. Joulain, J.-P. Mulet, Y. Chen, and J.-J. Greffet, Nature 444, 740 (2006).
  * [10] A. C. Jones and M. B. Raschke, Nano Lett. 12, 1475 (2012).
  * [11] R. S. DiMatteo, P. Greiff, S. L. Finberg, K. A. Young-Waithe, H. K. Choy, M. M. Masaki, and C. G. Fonstad, Appl. Phys. Lett. 79, 1894 (2001).
  * [12] A. Narayanaswamy and G. Chen, Appl. Phys. Lett. 82, 3544 (2003).
  * [13] S. Basu, Y.-B. Chen, and Z. M. Zhang, Int. J. Energy Res. 31, 689 (2007)
  * [14] A. Fiorino, L. Zhu, D. Thompson, R. Mittapally, P. Reddy, and E. Meyhofer, Nat. Nanotech. 13, 806 (2018).
  * [15] P. Ben-Abdallah and S.-A. Biehs, Phys. Rev. Lett. 112, 044301 (2013).
  * [16] P. Ben-Abdallah and S.-A. Biehs, AIP Adv. 5, 053502 (2015).
  * [17] C. R. Otey, W. T. Lau, and S. Fan, Phys. Rev. Lett. 104, 154301 (2010).
  * [18] S. Basu and M. Francoeur, Appl. Phys. Lett. 98, 113106 (2011).
  * [19] P. J. van Zwol, K. Joulain, P. Ben-Abdallah, and J. Chevrier, Phys. Rev. B 84, 161413(R) (2011).
  * [20] H. Iizuka and S. Fan, J. Appl. Phys. 112, 3303 (2012).
  * [21] L. P. Wang and Z. M. Zhang, Nanosc. Microsc. Therm. 17, 337 (2013).
  * [22] L. Zhu, C. R. Otey, and S. Fan, Phys. Rev. B 88, 184301 (2013).
  * [23] Y. Yang, S. Basu, and L. Wang, Appl. Phys. Lett. 103, 163101 (2013).
  * [24] J. Huang, Q. Li, Z. Zheng, and Y. M. Xuan, Int. J. Heat Mass Transfer 67, 575 (2013).
  * [25] K. Joulain, Y. Ezzahri, J. Drevillon, B. Rousseau, and D. Domingos, Opt. Express 23, A1388 (2015).
  * [26] W. Gu, G. H. Tang, and W. Q. Tao, Int. J. Heat Mass Tran. 82, 429 (2015).
  * [27] Y. Yang, S. Basu, and L. P. Wang, J. Quant. Spectrosc. Radiat. Transfer 158, 69 (2015).
  * [28] A. Ghanekar, J. Ji, and Y. Zheng, Appl. Phys. Lett. 109, 123106 (2016).
  * [29] L. Tang and M. Francoeur, Opt. Express 25, A1043 (2017).
  * [30] J. Ordoñez-Miranda, K. Joulain, M. D. De Sousa, Y. Ezzahri, and J. Drevillon, J. Appl. Phys. 122, 093105 (2017).
  * [31] Z. Zheng, X. Liu, A. Wang, and Y. M. Xuan, Int. J. Heat Mass Transfer 109, 63 (2017).
  * [32] A. Ghanekar, Y. Tian, M. Ricci, S. Zhang, O. Gregory, and Y. Zheng, Opt. Express 26, A209 (2018).
  * [33] J. Shen, X. Liu, H. He, W. Wu, and B. Liu, J. Quant. Spectrosc. Radiat. Transfer 211, 1 (2018).
  * [34] G. Xu, J. Sun, H. Mao, and T. Pan, J. Appl. Phys. 124, 183104 (2018).
  * [35] G. Xu, J. Sun, H. Mao, and T. Pan, J. Quant. Spectrosc. Radiat. Transfer 220, 140 (2018).
  * [36] S. Wen, X. Liu, S. Cheng, Z. Wang, S. Zhang, and C. Dang, J. Quant. Spectrosc. Radiat. Transfer 234, 1 (2019).
  * [37] A. Ott, R. Messina, P. Ben-Abdallah, and S.-A. Biehs, Appl. Phys. Lett. 114, 163105 (2019).
  * [38] O. R. Toyin, W. Ge, and L. Gao, Chin. Phys. Lett. 38, 016801 (2021).
  * [39] Q. Li, H. He, Q. Chen, and B. Song, Phys. Rev. Applied 16, 014069 (2021).
  * [40] E. Moncada-Villa and J. C. Cuevas, Phys. Rev. Applied 15, 024036 (2021).
  * [41] D. Feng, S. K. Yee, and Z. M. Zhang, Appl. Phys. Lett. 119, 181111 (2021).
  * [42] Y. Liu, Y. Tian, F. Chen, A. Caratenuto, X. Liu, M. Antezza, and Y. Zheng, Appl. Phys. Lett. 119, 123101 (2021).
  * [43] F. Chen, X. Liu, Y. Tian, and Y. Zheng, Adv. Eng. Mater. 23, 2000825 (2021).
  * [44] I. Latella, P. Ben-Abdallah, and M. Nikbakht, Phys. Rev. B 104, 045410 (2021).
  * [45] P. J. van Zwol, L. Ranno, and J. Chevrier, Phys. Rev. Lett. 108, 234301 (2012).
  * [46] K. Ito, K. Nishikawa, A. Miura, H. Toshiyoshi, and H. Iizuka, Nano Lett. 17, 4347 (2017).
  * [47] A. Fiorino, D. Thompson, L. Zhu, R. Mittapally, S.-A. Biehs, O. Bezencenet, N. El-Bondry, S. Bansropun, P. Ben-Abdallah, E. Meyhofer, and P. Reddy, ACS Nano 12, 5774 (2018).
  * [48] M. Elzouka and S. Ndao, Sci. Rep. 7, 44901 (2017).
  * [49] O. Ilic, M. Jablan, J. D. Joannopoulos, I. Celanovic, H. Buljan, and M. Soljačić, Phys. Rev. B 85, 155422 (2012).
  * [50] R. Messina, J.-P. Hugonin, J.-J. Greffet, F. Marquier, Y. De Wilde, A. Belarouci, L. Frechette, Y. Cordier, and P. Ben-Abdallah, Phys. Rev. B 87, 085421 (2013).
  * [51] V. B. Svetovoy, P. J. van Zwol, and J. Chevrier, Phys. Rev. B 85, 155418 (2012).
  * [52] O. Ilic, M. Jablan, J. D. Joannopoulos, I. Celanovic, and M. Soljačić, Opt. Express 20, A366 (2012).
  * [53] M. Lim, S. S. Lee, and B. J. Lee, Opt. Express 21, 22173 (2013).
  * [54] R. Messina and P. Ben-Abdallah, Sci. Rep. 3, 1383 (2013).
  * [55] V. B. Svetovoy and G. Palasantzas, Phys. Rev. Applied 2, 034006 (2014).
  * [56] R. Messina, P. Ben-Abdallah, B. Guizal, and M. Antezza, Phys. Rev. B 96, 045402 (2017).
  * [57] G. T. Papadakis, B. Zhao, S. Buddhiraju, and S. Fan, ACS Photonics 6, 709 (2019).
  * [58] M. M. Qazilbash, M. Brehm, B. G. Chae, P.-C. Ho, G. O. Andreev, B. J. Kim, S. J. Yun, A. V. Balatsky, M. B. Maple, F. Keilmann, H. T. Kim, and D. N. Basov, Science 318, 1750 (2007).
  * [59] P. Ben-Abdallah and K. Joulain, Phys. Rev. B 82, 121419(R) (2010).
  * [60] S.-A. Biehs, E. Rousseau, and J.-J. Greffet, Phys. Rev. Lett. 105, 234301 (2010).
  * [61] L. A. Falkovsky, J. Phys. Conf. Ser. 129, 012004 (2008).
  * [62] A. S. Barker, H. W. Verleur, and H. J. Guggenheim, Phys. Rev. Lett. 17, 1286 (1966).
  * [63] _Handbook of Optical Constants of Solids_ , edited by E. Palik (Academic Press, New York, 1998).
  * [64] P. Rodriguez-López, W.-K. Tse, and D. A. R. Dalvit, J. Phys.: Condens. Matter 27, 214019 (2015).


# Electrochemistry meets polymer physics: polymerized ionic liquids on an
electrified electrode

Yury A. Budkov<EMAIL_ADDRESS>School of Applied Mathematics, HSE University,
Tallinskaya st. 34, 123458 Moscow, Russia  G.A. Krestov Institute of Solution
Chemistry of the Russian Academy of Sciences, 153045, Akademicheskaya st. 1,
Ivanovo, Russia  Nikolai N. Kalikin  G.A. Krestov Institute of Solution
Chemistry of the Russian Academy of Sciences, 153045, Akademicheskaya st. 1,
Ivanovo, Russia  Andrei L. Kolesnikov  Institut für Nichtklassische Chemie
e.V., Permoserstr. 15, 04318 Leipzig, Germany

###### Abstract

Polymeric ionic liquids are emerging polyelectrolyte materials for modern
electrochemical applications. In this paper, we propose a self-consistent
field theory of the polymeric ionic liquid on a charged conductive electrode.
Taking into account the conformation entropy of rather long polymerized
cations within the Lifshitz theory and electrostatic and excluded volume
interactions of ionic species within the mean-field approximation, we obtain a
system of self-consistent field equations for the local electrostatic
potential and average concentrations of monomeric units and counterions. We
solve these equations in the linear approximation for the cases of a point-
like charge and a flat infinite uniformly charged electrode immersed in a
polymeric ionic liquid and derive analytical expressions for local ionic
concentrations and electrostatic potential. We also derive an analytical
expression for the linear differential capacitance of the electric double
layer. Solving numerically the self-consistent field equations, we establish a
range of the potential drop, when the linear theory is valid and analyze the
differential capacitance as a function of applied voltage for a pure polymeric
ionic liquid and a polymeric ionic liquid dissolved in a highly polar organic
solvent. For the first case, we obtain a slightly asymmetric bell-shaped
differential capacitance profile, whereas for the second one – strongly
asymmetric camel-shaped profile. We elucidate the nature of this asymmetry,
analyzing the behavior of the volume fractions of monomeric units and
counterions on the electrode at different voltages. Based on the numerical
solution of the self-consistent field equations, we show that at a rather
large potential drop in both considered cases, the differential capacitance
has the same asymptotic behavior as in regular "monomeric" ionic liquids,
predicted within Kornyshev’s analytical self-consistent field theory. Our
findings can be used in different electrochemical applications, such as
supercapacitors, batteries, fuel cells, electrodeposition, etc.

##  I Introduction

Ionic liquids (ILs) at charged electrodes have already been quite thoroughly
studied experimentally and theoretically due to their wide range of
applications as "designer solvents" and "filling" for supercapacitors Fedorov
and Kornyshev (2014). Their polymeric forms (PIL) are a relatively new-found
Ohno and Ito (1998); Ohno, Yoshizawa, and Ogihara (2004); Ohno (2007) class of
materials with a shallow theoretical background, although demonstrating a wide
spectrum of interesting properties of both ionic liquids and polymeric
systems, which can lead to their potential application as ionic and solid-
state electrolytes, fuel cell membranes, dispersants, ordered nano- and
mesostructures, sorbents, carbon precursors, gating material in field-effect
transistors, systems for water treatment, osmotic energy conversion, solar
modulation, capillary electrophoresis, electrochromatography, among others
Yuan and Antonietti (2011); Mecerreyes (2011); Eftekhari and Saito (2017);
Fadeeva  _et al._ ; Chen  _et al._ (2018); Peltekoff  _et al._ (2019); Gupta,
Liang, and Hu (2019); Yuki Kohno  _et al._ (2015); Tang  _et al._ (2014); Gao
_et al._ (2021).

The theoretical investigations on the charged polymers near charged interfaces
are practically limited to the cases of polyelectrolyte adsorption modeling
(see reviews Dobrynin and Rubinstein (2005); Netz and Andelman (2003) and
references therein). It is instructive to discuss the existing theoretical
models of polyelectrolyte solutions on a charged surface. One way to
theoretically consider such systems is multi-Stern layer theory Van Der Schee
and Lyklema (1984); Papenhuijzen, Van Der Schee, and Fleer (1985); Evers  _et
al._ (1986) \- a lattice model extending classical works Roe (1974);
Scheutjens and Fleer (1979, 1980); Borisov, Zhulina, and Birshtein (1994);
Friedsam, Gaub, and Netz (2005) on uncharged polymer adsorption modeling by
incorporating the electrostatic free energy in the partition function of the
system or modeling of a single polyelectrolyte chain adsorption on a charged
surface. Another approach treats the electrostatic potential and the
concentration of the polyelectrolyte as continuous functions Varoqui, Johner,
and Elaissari (1991); Varoqui (1993); Borukhov, Andelman, and Orland (1995);
Châtellier and Joanny (1996); Borukhov, Andelman, and Orland (1998); Joanny
(1999); Podgornik and Jönsson (1993); Brilliantov, Budkov, and Seidel (2016).
Taking into account the free energy contributions of the monomeric unit
excluded volume via the second virial coefficient, the entropy of small ions
and electrostatic interactions and minimizing it with respect to the
concentrations and the potential, the authors obtained a system of two
differential equations. The first of them is a generalized Poisson-Boltzmann
equation and the second one – a generalization of the self-consistent field
equation of neutral polymers. In work Borukhov, Andelman, and Orland (1995)
the authors theoretically investigated the behavior of a polyelectrolyte
solution confined between two charged surfaces for several cases of monomeric
charge distributions. The authors described the system within the mean-field
approach, introducing a system of coupled differential equations, consisting
of a modified Poisson-Boltzmann equation for the electrostatic potential and a
self-consistent field equation for the polymer order parameter (square root of
monomeric units concentration). Previously, similar results for the smeared
case of the monomeric charge distribution without accounting for the monomeric
unit excluded volume were obtained in Varoqui, Johner, and Elaissari (1991);
Varoqui (1993); Bohmer, Evers, and Scheutjens (1990). A numerical calculation
of the obtained combined equations allowed the authors to compare the monomer
concentration profiles with different types of charge distributions and also
to conclude that the surface charge drastically influences the amount of
polymer confined between the plates. In order to continue their theoretical
investigation of the polyelectrolyte behavior near a single charged surface,
based on the obtained numerical results, the authors Borukhov, Andelman, and
Orland (1998) also derived expressions for the amount of the adsorbed polymers
for high- and low-salt regimes. In the high-salt case, the authors obtained
two limiting regimes dependent on the value of the charged monomeric unit
fraction, i.e. polyelectrolyte strength.

A study of polyelectrolyte solution adsorption at the ideal charged interface
within the mean-field theory Châtellier and Joanny (1996) in the linear
approximation led to the analytical results for the polymer concentration
profiles, which were shown to oscillate in the case of low salt
concentrations. Further development of the self-consistent field theory
approach Joanny (1999), based on the ground state dominance approximation for
the polymer order parameter of the Gaussian chains and consideration of the
small surface charge, which justifies the linearization of the self-consistent
field equation for the electrostatic potential, allowed the author to describe
the polyelectrolyte adsorption in two limiting regimes of low and high ionic
strength. The model successfully described the surface charge overcompensation
and the following inversion of the surface charge due to the polyelectrolyte
adsorption.

A nonlocal density functional theory of the polyelectrolyte solution near a
charged planar surface was developed in Li and Wu (2006a, b). The account of
the hard-sphere contribution within the fundamental measure theory, the chain
connectivity via the thermodynamic perturbation theory, and electrostatic
correlations allowed the authors to reach good agreement with the results of
simulations for ion and polyion segment density distribution as well as
polyionic surface excess at varying properties. They also studied the changes
in the ion valence, polyion chain length and difference of the ion and polyion
sizes on the microscopic structure and mean electrostatic potential, while
accurately capturing the charge inversion and layering effects.

Two theoretical papers Kumar  _et al._ (2017, 2020) have been published, where
the investigation of the electrode-PIL interface by broadband dielectric
spectroscopy with theoretical modeling, based on the Rayleighian dissipation
function formalism, allowed the authors to conclude that at zero potential
there is a preabsorbed layer that influences the results of the impedance and
capacitance measuring. Interestingly, the steady-state capacitance curves do
not show pronounced asymmetry, and the overall shape is qualitatively the same
as that of the corresponding IL.

To sum up the above, there are a decent number of studies devoted to
theoretical modeling of charged polymer adsorption onto a charged surface.
However, to the best of our knowledge, there are no thermodynamic models of
EDL on a PIL-conductive electrode interface similar to those formulated for
the conventional monomeric ILs Fedorov and Kornyshev (2014); Kornyshev (2007);
Budkov  _et al._ (2018); Wu  _et al._ (2011); Jiang, Meng, and Wu (2011);
Forsman, Woodward, and Trulsson (2011); Goodwin, Feng, and Kornyshev (2017);
Maggs and Podgornik (2016); May (2019); Cruz  _et al._ (2019); Bazant, Storey,
and Kornyshev (2011); Avni, Adar, and Andelman (2020); Budkov, Zavarzin, and
Kolesnikov (2021), describing differential capacitance as a function of
applied voltage. Moreover, most of the models discussed above deal with dilute
or semi-dilute polyelectrolyte solutions, accounting for the excluded volume
interactions via the second virial term.

In this paper, we propose a self-consistent field theory of the PIL condensed
phase on a charged conductive electrode. Taking into account the conformation
entropy of rather long polymerized cations, electrostatic and excluded volume
interactions of ionic species, we will derive a system of self-consistent
field equations for the local electrostatic potential and local concentrations
of the monomeric units and counterions. We will analyze the behavior of the
EDL differential capacitance as a function of applied voltage and compare it
with the differential capacitance behavior predicted in Kornyshev’s analytical
mean-field theory of low-molecular-weight ILs with the same physical
parameters.

##  II Self-consistent field theory for polymeric ionic liquids

In this section, we will formulate a self-consistent field (SCF) theory of PIL
electrostatically disturbed by some external charges fixed in space with the
density $\rho_{ext}(\mathbb{r})$. Let us consider a PIL consisting of
polymerized cations, whose monomeric units carry a charge $q>0$ and low-
molecular-weight anions with a charge $-q$. We assume that a polymerization
degree of the cations is very large ($N\gg 1$) and we can neglect the
translation entropy of the mass center of the polymer chains. Thus, the grand
thermodynamic potential of the PIL is

$\Omega=F_{el}+\Omega_{ref},$ (1)

where

$F_{el}=\frac{\varepsilon\varepsilon_{0}}{2}\int
d\mathbb{r}\psi(\mathbb{r})\Delta\psi(\mathbb{r})+\int
d\mathbb{r}\rho(\mathbb{r})\psi(\mathbb{r})$ (2)

is the electrostatic energy in the mean-field approximation with the local
charge density

$\rho(\mathbb{r})=\rho_{ext}(\mathbb{r})+q\left(n_{p}(\mathbb{r})-n_{c}(\mathbb{r})\right),$
(3)

where the first term is the local charge density of the external charges and
the second term is the local charge density of the monomer segments and
counterions; $\varepsilon$ is the dielectric permittivity of polymeric ionic
liquid taking into account polarization of the medium, $\varepsilon_{0}$ is
the empty space permittivity, $\Delta$ is the Laplace operator. The grand
thermodynamic potential of the reference polymer system can be written as
follows

$\Omega_{ref}=F_{conf}+F_{vol}-\int
d\mathbb{r}\left(\mu_{p}n_{p}(\mathbb{r})+\mu_{c}n_{c}(\mathbb{r})\right),$
(4)

where

$F_{conf}[n_{p}]=-\frac{k_{B}Tb^{2}}{6}\int
d\mathbb{r}n_{p}^{1/2}(\mathbb{r})\Delta n_{p}^{1/2}(\mathbb{r})$ (5)

is the Lifshitz conformation entropy Lifshitz (1969); Grosberg and Khokhlov
(1994) of the flexible polymer chains with a bond length $b$; $k_{B}$ is the
Boltzmann constant, $T$ is the temperature. The contribution of the volume
interactions of monomeric units and counterions is

$F_{vol}[n_{p},n_{c}]=\int
d\mathbb{r}f(T,n_{p}(\mathbb{r}),n_{c}(\mathbb{r}))$ (6)

where $f=f(T,n_{p}(\mathbb{r}),n_{c}(\mathbb{r}))$ is the free energy density
as a function of local concentrations, $n_{p,c}(\mathbb{r})$, of monomeric
units and counterions, which can be estimated within the asymmetric lattice
gas model Sanchez and Lacombe (1978); Maggs and Podgornik (2016)

$f=\frac{k_{B}T}{v}\left(\phi_{c}\ln\phi_{c}+\left(1-\phi_{c}-\phi_{p}\right)\ln\left(1-\phi_{c}-\phi_{p}\right)\right),$
(7)

where $\phi_{c}=n_{c}v$ and $\phi_{p}=n_{p}v$ are the local volume fractions
of the counterions and monomeric units, respectively, $v$ is the elementary
cell volume related to the bond length by the condition natural for a lattice
model, $v=b^{3}$; $\mu_{p}$ and $\mu_{c}$ are the bulk chemical potentials of
the monomeric units and counterions, respectively. The Euler-Lagrange
equations, $\delta\Omega/\delta n_{p,c}(\mathbb{r})=0$ and
$\delta\Omega/\delta\psi(\mathbb{r})=0$, have the following form

$\displaystyle\bar{\mu}_{c}(\mathbb{r})-q\psi(\mathbb{r})=\mu_{c}$
$\displaystyle\bar{\mu}_{p}(\mathbb{r})-\frac{k_{B}Tb^{2}}{6n_{p}^{1/2}(\mathbb{r})}\Delta
n_{p}^{1/2}(\mathbb{r})+q\psi(\mathbb{r})=\mu_{p}$ (8)
$\displaystyle-\varepsilon\varepsilon_{0}\Delta\psi(\mathbb{r})=\rho_{ext}(\mathbb{r})+q\left(n_{p}(\mathbb{r})-n_{c}(\mathbb{r})\right),$

where

$\bar{\mu}_{c}=\frac{\partial f}{\partial
n_{c}}=k_{B}T\ln\left(\frac{\phi_{c}}{1-\phi_{c}-\phi_{p}}\right),$ (9)
$\bar{\mu}_{p}=\frac{\partial f}{\partial
n_{p}}=-k_{B}T\ln\left(1-\phi_{c}-\phi_{p}\right)$ (10)

are the intrinsic chemical potentials of the monomeric units and counterions,
respectively. Taking into account that in the bulk solution, where $\psi=0$,
the local electroneutrality condition, $n_{p}=n_{c}=n_{0}$, is fulfilled, we
obtain the following expressions for the bulk chemical potentials of the
species

$\mu_{c}=k_{B}T\ln\left(\frac{\phi_{0}}{1-2\phi_{0}}\right),$ (11)
$\mu_{p}=-k_{B}T\ln\left(1-2\phi_{0}\right),$ (12)

where $\phi_{0}=n_{0}v$ is the bulk volume fraction of the monomeric units and
counterions.

##  III Linear theory

Unfortunately, system (II) of SCF equations obtained above for the local
concentrations of monomeric units and counterions, $n_{p,c}(\mathbb{r})$, and
electrostatic potential, $\psi(\mathbb{r})$, cannot be, in general,
analytically solved. However, in a similar way as in the classical Debye-
Hueckel theory of simple electrolyte solutions, we can analyze these equations
in the limit of weak electrostatic interactions, where the local
concentrations of monomeric units and counterions only slightly deviate from
the bulk value, $n_{0}$, i.e.

$n_{p}=n_{0}+\bar{n}_{p},\leavevmode\nobreak\ n_{c}=n_{0}+\bar{n}_{c}.$ (13)

where we have introduced small perturbations of concentrations
($|\bar{n}_{p,c}|\ll n_{0}$). Thus, introducing the small dimensionless
variables $x_{p}={\bar{n}_{p}}/{n_{0}}$, $\bar{\psi}={q\psi}/{k_{B}T}$,
$x_{c}={\bar{n}_{c}}/{n_{0}}$, we can linearize SCF equations (II) and, after
some algebra, obtain the following system of linear equations

$\displaystyle\Delta^{2}x_{p}-\lambda^{2}\Delta
x_{p}+\frac{\chi^{4}}{4}x_{p}=-\frac{12q}{\varepsilon\varepsilon_{0}k_{B}Tb^{2}(1-\phi_{0})}\rho_{ext}$
$\displaystyle\bar{\psi}=\frac{b^{2}}{12}(1-\phi_{0})\Delta
x_{p}-\phi_{0}x_{p}$ (14) $\displaystyle
x_{c}=\frac{(1-2\phi_{0})\bar{\psi}-\phi_{0}x_{p}}{1-\phi_{0}},$

where we have also introduced the auxiliary parameters

$\lambda^{2}=\frac{12\phi_{0}}{b^{2}(1-\phi_{0})}+\frac{1-2\phi_{0}}{r_{D}^{2}(1-\phi_{0})},\leavevmode\nobreak\
\chi^{4}=\frac{48}{r_{D}^{2}b^{2}}\frac{1+2\phi_{0}}{1-\phi_{0}},$ (15)

with the Debye radius,
$r_{D}=\left({q^{2}n_{0}}/{\varepsilon\varepsilon_{0}k_{B}T}\right)^{-1/2}$,
which comes from the presence of the mobile counterions. In what follows, we
will analyze the solution of the system of linear equations (III) for two
physically relevant cases: (A) a point-like charge and (B) a flat charged
electrode, immersed in a PIL.

###  III.1 Case of point-like charge

Here we would like to solve eqs. (III) for the case, when the polymeric ionic
liquid is disturbed by the fixed point-like test charge $q_{0}$ placed at the
origin. Taking into account that the charge density of the test point-like
charge is

$\rho_{ext}(\mathbf{r})=q_{0}\delta(\mathbf{r}),$ (16)

where $\delta(\mathbb{r})$ is the Dirac delta-function, we obtain

$\Delta^{2}x_{p}-\lambda^{2}\Delta
x_{p}+\frac{\chi^{4}}{4}x_{p}=-\frac{12qq_{0}}{\varepsilon\varepsilon_{0}k_{B}Tb^{2}(1-\phi_{0})}\delta(\mathbf{r}).$
(17)

Further, using the integral Fourier transformation

$x_{p}(\mathbf{r})=\int\frac{d\mathbf{k}}{(2\pi)^{3}}e^{i\mathbf{k}\mathbf{r}}\tilde{x}_{p}(\mathbf{k}),$
(18)

we arrive at the following expression for the Fourier-image

$\tilde{x}_{p}(\mathbf{k})=-\frac{12qq_{0}}{\varepsilon\varepsilon_{0}k_{B}Tb^{2}(1-\phi_{0})}\frac{1}{k^{4}+\lambda^{2}k^{2}+\frac{\chi^{4}}{4}},$
(19)

We only consider here the physically relevant case, $\lambda<\chi$, which can
be reformulated as the following condition for the Debye radius and bond
length of the chain:

$\frac{\sqrt{(1+2\phi_{0})(1-\phi_{0})}-1}{2\sqrt{3}\phi_{0}}<\frac{r_{D}}{b}<\frac{\sqrt{(1+2\phi_{0})(1-\phi_{0})}+1}{2\sqrt{3}\phi_{0}}.$
(20)

Condition (20) covers the cases of pure PILs and PILs dissolved in a polar
organic solvent like acetonitrile (see the physical parameters below). In this
case, the inverse Fourier-transform yields

$x_{p}=a\frac{\exp\left[-\varkappa r\right]}{r}\sin{kr},$ (21)

where we have introduced the constants
$\varkappa={\sqrt{\chi^{2}+\lambda^{2}}}/{2}$,
$k={\sqrt{\chi^{2}-\lambda^{2}}}/{2}$ (see the discussion of their physical
sense below) and coefficient

$a=-\frac{3qq_{0}}{2\pi\varepsilon\varepsilon_{0}k_{B}Tb^{2}\varkappa
k(1-\phi_{0})}.$ (22)

The dimensionless electrostatic potential has the following form

$\bar{\psi}=\frac{\exp\left[-\varkappa
r\right]}{r}(b_{1}\sin{kr}+b_{2}\cos{kr}),$ (23)

where

$b_{1}=a\left(\theta-\phi_{0}\right),\leavevmode\nobreak\ b_{2}=-a\eta,$ (24)

where $\theta=\chi^{2}b^{2}(1-\phi_{0})/24$ and $\eta=\varkappa
kb^{2}(1-\phi_{0})/6$.

The perturbation of the cation concentration is determined by

$x_{c}=\frac{\exp\left[-\varkappa r\right]}{r}(c_{1}\sin{kr}+c_{2}\cos{kr})$
(25)

with the coefficients

$c_{1}=\frac{(1-2\phi_{0})b_{1}-\phi_{0}a}{1-\phi_{0}},\leavevmode\nobreak\
c_{2}=\frac{1-2\phi_{0}}{1-\phi_{0}}b_{2}.$ (26)

It is interesting to note that for $\phi_{0}=0.5$ (there are no voids in the
bulk liquid structure) we get $x_{p}=-x_{c}$ or $n_{p}+n_{c}=2n_{0}$. To the
best of our knowledge, eqs. (25) and (21) are new analytical results. We would
like to note that expression, similar to eq. (23), was obtained within the
field theories of polyelectrolyte solutions in papers Borue and Erukhimovich
(1988); Muthukumar (1996).

We would like to note that the length $\varkappa^{-1}$ characterizes the
spatial decay of the perturbations of the concentrations of monomeric units
and counterions, and electrostatic potential relative to the bulk values, and
is analogous to the Debye radius in simple electrolyte solutions. The length
$k^{-1}$ can be considered as a measure of the $"$periodicity$"$ of positive
and negative charge-rich domains, which is the result of the fine interplay
between the conformation entropy of the polymer chains, electrostatic
interactions, and excluded volume interactions of monomeric units. We would
like to mention that similar behavior of the electrostatic potential of the
point-like charge has been recently predicted in salt solutions of quadrupolar
molecules Budkov (2020, 2019); Slavchov (2014) and ionic liquids with an
account of the quadrupolar ionic clusters Avni, Adar, and Andelman (2020).

###  III.2 Case of a flat charged electrode

Now, we would like to consider the above linear response theory for the case
of a flat infinite uniformly charged electrode with the surface charge density
$\sigma$. The electrode charge density is

$\rho_{ext}(z)=\sigma\delta(z),$ (27)

so that the equation for the perturbation of monomeric unit concentration has
the form

$x_{p}^{(IV)}(z)-\lambda^{2}x_{p}^{\prime\prime}(z)+\frac{\chi^{4}}{4}x_{p}(z)=-\frac{12q\sigma}{\varepsilon\varepsilon_{0}k_{B}Tb^{2}(1-\phi_{0})}\delta(z),$
(28)

where the superscript ${(IV)}$ denotes the fourth derivative. The boundary
conditions can be written as follows

$x_{p}^{\prime}(0)=0,\leavevmode\nobreak\
x_{p}^{\prime\prime\prime}(0)=-\frac{12q\sigma}{\varepsilon\varepsilon_{0}k_{B}Tb^{2}(1-\phi_{0})},\leavevmode\nobreak\
x_{p}(\infty)=x_{p}^{\prime}(\infty)=0,$ (29)

where we adopt that there is no physical adsorption of polymer chains on the
electrode surface, so that the boundary condition $n_{p}^{\prime}(0)=0$ is
fulfilled Netz and Andelman (2003) and take into account the delta-function
term in the right-hand side of eq. (28). The generalization of the theory for
the case of physical adsorption of monomeric units on the electrode is
straightforward. The solution of eq. (28) with the help of the above boundary
condition results in

$x_{p}(z)=e^{-\varkappa z}\left(f_{1}\cos(kz)+f_{2}\sin(kz)\right),$ (30)

with the following coefficients

$f_{1}=-\frac{6q\sigma}{\varkappa(\varkappa^{2}+k^{2})\varepsilon\varepsilon_{0}k_{B}Tb^{2}(1-\phi_{0})},\leavevmode\nobreak\
f_{2}=-\frac{6q\sigma}{k(\varkappa^{2}+k^{2})\varepsilon\varepsilon_{0}k_{B}Tb^{2}(1-\phi_{0})}.$
(31)

The dimensionless electrostatic potential takes the following form

$\bar{\psi}=e^{-\varkappa z}\left(g_{1}\cos(kz)+g_{2}\sin(kz)\right),$ (32)

where

$g_{1}=\left(\theta-\phi_{0}\right)f_{1}-\eta f_{2},$ (33)
$g_{2}=\left(\theta-\phi_{0}\right)f_{2}+\eta f_{1}.$ (34)

The perturbation of the counterion concentration has the form

$x_{c}=e^{-\varkappa z}\left(h_{1}\cos(kz)+h_{2}\sin(kz)\right),$ (35)

with the following coefficients linearly expressed via the ones written above

$h_{1}=\frac{(1-2\phi_{0})g_{1}-\phi_{0}f_{1}}{1-\phi_{0}},\leavevmode\nobreak\
h_{2}=\frac{(1-2\phi_{0})g_{2}-\phi_{0}f_{2}}{1-\phi_{0}}.$ (36)

As in the above case of the point-like test charge, it is interesting to note
that for $\phi_{0}=0.5$ we obtain $x_{p}=-x_{c}$ or
$n_{c}(z)+n_{p}(z)=2n_{0}$.

Introducing the potential drop on the electrode $\psi_{0}=\psi(0)$, using eq.
(23), after some algebraic calculations, we obtain the following expression
for the surface charge density

$\sigma=C_{0}\psi_{0},$ (37)

with the linear differential capacitance of the EDL

$C_{0}=\frac{\varepsilon\varepsilon_{0}}{d},$ (38)

where

$d=\left(1+\frac{\phi_{0}}{\theta}\right)\left(\lambda^{2}+\chi^{2}\right)^{-1/2}$
(39)

can be interpreted as effective thickness of the diffuse EDL. We would like to
note that the thickness is proportional to determined above the effective
screening length, $\varkappa^{-1}$. Eqs. (38) and (39) are the main analytical
result of this paper. It is interesting to note that in contrast to the case
of low-molecular-weight ionic liquids described within the symmetric lattice
gas model Kornyshev (2007), the differential capacitance, $C_{0}$, of the
polymeric ionic liquids in the linear theory depends on the molecular
parameters of the ionic species, namely, the bond length, $b$, and the
effective excluded volume of particles, $v$. In the limit of negligible
excluded volume interactions ($v\to 0$) at constant bond length, $b$, we
arrive at the expression

$C_{0}=C_{GC}\sqrt{1+\frac{4\sqrt{3}r_{D}}{b}},$ (40)

which is an analog of the well-known linear Gouy-Chapman or Debye linear
capacitance Kornyshev (2007), $C_{GC}=\varepsilon\varepsilon_{0}/r_{D}$, for
low-molecular-weight electrolyte solutions and monomeric ionic liquids. Note
that eq. (40) should be valid for rather dilute polyelectrolyte solutions on a
charged electrode in the theta-solvent regime, when the polymer chains behave
as the ideal Gaussian coils Grosberg and Khokhlov (1994). In the region
$r_{D}\ll b$, eq. (40) transforms into the standard Gouy-Chapman expression,
i.e. $d\simeq r_{D}$. However, when $r_{D}\gg b$, the electric double layer
thickness is determined by $d\simeq 48^{-1/4}(r_{D}b)^{1/2}$. Note that a
decrease in the bond length, $b$, as is seen from eq. (40), increases the
linear differential capacitance. Such behavior can be explained by the fact
that the decrease in the bond length leads to higher linear charge density of
the polymer chains and, thus, higher accumulated charge in the EDL.

##  IV Polymeric ionic liquid on a charged electrode: nonlinear regime

###  IV.1 Statement of problem

Now we would like to consider the case of a PIL on a flat charged electrode
with a constant surface potential, $\psi_{0}$. The self-consistent field
equations take the following simplified form

$\displaystyle\bar{\mu}_{p}-\frac{k_{B}Tb^{2}}{6\phi_{p}^{1/2}}(\phi_{p}^{1/2})^{\prime\prime}+q\psi=\mu_{p}$
$\displaystyle\psi^{\prime\prime}=-\frac{q}{\varepsilon\varepsilon_{0}v}\left(\phi_{p}-\phi_{c}\right)$
(41)
$\displaystyle\phi_{c}=(1-\phi_{p})\frac{e^{\frac{\mu_{c}+q\psi}{k_{B}T}}}{1+e^{\frac{\mu_{c}+q\psi}{k_{B}T}}}$

with the boundary conditions

$\psi(0)=\psi_{0},\leavevmode\nobreak\
\psi^{\prime}(\infty)=0,\leavevmode\nobreak\
\phi_{p}^{\prime}(0)=0,\leavevmode\nobreak\ \phi_{p}(\infty)=\phi_{0}.$ (42)

Solving eqs. (IV.1) at different potential drops, $\psi_{0}$, we can calculate
the surface charge density as
$\sigma=-\varepsilon\varepsilon_{0}\psi^{\prime}(0)$ and then the differential
capacitance, $C=d\sigma/d\psi_{0}$.

###  IV.2 Numerical results and discussions

Solving numerically the self-consistent field equations (IV.1), we can
calculate the differential capacitance as a function of applied voltage. Fig.
1 shows a comparison of the differential capacitance profile of a PIL
numerically calculated within the present theory and differential capacitance
for monomeric IL determined by the following analytical expression derived by
Kornyshev in ref. Kornyshev (2007)

$C=\frac{\left({2\beta\varepsilon\varepsilon_{0}c^{2}vq^{2}}\right)^{1/2}|\sinh{\beta
q\psi_{0}}|}{\left(1+4\phi_{0}\sinh^{2}\left(\frac{\beta
q\psi_{0}}{2}\right)\right)\sqrt{\ln\left(1+4\phi_{0}\sinh^{2}\left(\frac{\beta
q\psi_{0}}{2}\right)\right)}}.$ (43)

Figure 1: Differential capacitance profile obtained numerically by present
theory for pure PIL and differential capacitance predicted by the analytical
mean-field theory for a monomeric IL (see eq. (43)). The data are shown for
$\phi_{0}=0.4$, $b=v^{1/3}=0.5\leavevmode\nobreak\ nm$,
$T=298\leavevmode\nobreak\ K$, $\varepsilon=5$, $q=1.6\times
10^{-19}\leavevmode\nobreak\ C$. Figure 2: Local volume fractions of monomeric
units and counterions on the electrode as functions of the potential drop for
the cases of a pure polymeric ionic liquid and and a polymeric ionic liquid
solution. The parameters for the pure polymeric ionic liquid (blue lines) are
$\phi_{0}=0.4$, $b=v^{1/3}=0.5\leavevmode\nobreak\ nm$,
$T=298\leavevmode\nobreak\ K$, $\varepsilon=5$. The parameters for the
polymeric ionic liquid solution (black lines) are $\phi_{0}=0.1$,
$b=v^{1/3}=0.5\leavevmode\nobreak\ nm$, $T=298\leavevmode\nobreak\ K$,
$\varepsilon=40$,$q=1.6\times 10^{-19}\leavevmode\nobreak\ C$.

As is seen, as in case of a monomeric IL, the PIL differential capacitance
profile has a quite symmetric bell shape (see also the discussion below). It
is interesting to note that at sufficiently large potential drops, as in case
of a monomeric IL, the differential capacitance monotonically decreases as
$C\sim|\psi_{0}|^{-0.5}$. Moreover, at rather large positive potential drops,
the differential capacitance profile coincides with the one predicted by
Kornyshev’s theory. At sufficiently large positive surface potentials, the
monomeric units are fully expelled from the electrode surface by the
counterions (see Fig. 2). Thus, the EDL has the same local structure as in the
case of monomeric ILs, which explains the correspondence to Kornyshev’s
theory. At rather large negative surface potentials, the differential
capacitance values are also close to those predicted by the Kornyshev mean-
field theory. This can be explained by the fact that at rather large local
volume fractions of the monomeric units near the electrode ($\phi_{p}\sim 1$),
occurring at sufficiently large potential drops, there is no difference
between the tied charged monomeric units and freely moving ions from the
entropic standpoint – the translational entropy in both cases is negligible.

Figure 3: Differential capacitance profiles predicted by the present theory
for a PIL solution and predicted by Kornyshev’s mean-field theory for a simple
electrolyte solution (eq. (43)). The data are shown for $\phi_{0}=0.1$,
$b=v^{1/3}=0.5\leavevmode\nobreak\ nm$, $T=298\leavevmode\nobreak\ K$,
$\varepsilon=40$,$q=1.6\times 10^{-19}\leavevmode\nobreak\ C$.

Now it is instructive to discuss the differential capacitance profile realized
for PILs dissolved in a highly polar organic solvent with the dielectric
permittivity $\varepsilon\approx 40$. We assume the bulk volume fraction of
the monomeric units to be $\phi_{0}=0.1$. Fig. 3 demonstrates a comparison
between the differential capacitance profile predicted by the present theory
and Kornyshev’s mean-field theory for the cases of a PIL solution and a simple
electrolyte solution, respectively. Unlike the case of a pure PIL discussed
above, where we obtained quite a symmetric bell-shaped capacitance profile,
the differential capacitance profile for the PIL solution case has two
strongly asymmetric peaks. The higher one is located in the region of negative
voltages, whereas the less pronounced peak is in the region of positive
voltages. As is seen, such capacitance behavior drastically differs from the
behavior observed in simple electrolyte solutions, predicted by Kornyshev’s
theory for the same physical parameters. In particular, at negative voltages,
when there are more positively charged monomeric units than counterions in the
vicinity of the electrode surface, we obtain much higher differential
capacitance values than those realized for the low-molecular-weight
electrolyte solution. Such abnormal differential capacitance behavior can be
explained by the fact that an increase in the negative voltage value leads to
a rather dramatic increase in the local volume fraction of monomeric units on
the electrode, $\phi_{p}(0)$, and the same decrease in the local volume
fraction of the counterions, $\phi_{c}(0)$ (see Fig. 2). However, an increase
in the positive voltage leads to a slower growth in the local volume fraction
of the counterions and decrease in the local volume fraction of monomeric
units. Such difference is associated with the fact that in the case of very
large polymerization degree the long polymer chains have zero translation
entropy and, therefore, can be more easily attracted to the negatively charged
electrode than the mobile counterions to the positively charged electrode. We
would like to note that in case of a PIL solution, as in case of a pure PIL,
the differential capacitance at sufficiently high potential drop values is
described by the same power law, $C\sim|\psi_{0}|^{-0.5}$.

Now we turn to the discussion of how a variation in the bond length,
$b=v^{1/3}$, influences the differential capacitance. As is seen from Fig. 4,
an increase in the bond length results in a monotonic decrease in the
differential capacitance at all potential drops. Such differential capacitance
behavior can be explained by the fact that the EDL accumulates a smaller
charge with the increase in the bond length and excluded volume of particles
due to stronger steric interactions of the ionic species in the EDL. It is
also interesting to note that the differential capacitance maximum (see also
Fig. 1) is slightly shifted to the negative potential drop values. Thus the
predicted bell-shaped capacitance profile can be considered a symmetric one
with good accuracy. This is in agreement with the conclusion of paper Kumar
_et al._ (2020) discussed in the Introduction. Note that the same behavior
with an increase in the bond length is observed in the PIL solution discussed
above – differential capacitance increases with the bond length.

Figure 4: Differential capacitance profiles of a pure PIL on a charged
electrode plotted for different bond lengths, $b=v^{1/3}$. The horizontal
dashed lines determine the linear capacitance (see eq. (38)). The data are
shown for $\phi_{0}=0.4$, $T=298\leavevmode\nobreak\ K$,
$\varepsilon=5$,$q=1.6\times 10^{-19}\leavevmode\nobreak\ C$.

Figure 5: Electrostatic potential (top) and concentration (bottom) profiles
calculated within the linear analytical theory and nonlinear theory for
positive potential drops. The data are shown for $\phi_{0}=0.4$,
$b=v^{1/3}=0.5\leavevmode\nobreak\ nm$, $T=298\leavevmode\nobreak\ K$,
$\varepsilon=5$, $q=1.6\times 10^{-19}\leavevmode\nobreak\ C$. The local
concentrations are determined via $x_{p,c}(z)=(n_{p,c}(z)-n_{0})/n_{0}$.

Now we would like to discuss the electrostatic potential profiles and local
ionic concentrations following from the numerical solution of the SCF
equations (IV.1) and compare them with those predicted by the linear theory.
As is seen from Fig. 5, at positive voltages the potential profile has a
minimum which becomes more pronounced with a voltage increase. Such potential
profile behavior is determined by an interplay of electrostatics and
conformation entropy of the macromolecules similar to the case of the
electrostatic potential of a point-like charge immersed in a PIL (see eq.
(23)) or a polyelectrolyte solution Borue and Erukhimovich (1988); Muthukumar
(1996). Note that the potential profiles predicted by the linear theory
practically coincide with the numerical solutions at the voltages in the
interval $|\psi_{0}|\leq 0.05\leavevmode\nobreak\ V$. At the voltages
$|\psi_{0}|>0.05\leavevmode\nobreak\ V$ the linear theory already gives
significant discrepancy for the local electrostatic potential values. However,
the interval of potential drops, where the linear theory gives a satisfactory
approximation for the concentration profiles is even narrower. Indeed, already
at $\psi_{0}\approx\pm 0.05\leavevmode\nobreak\ V$ the local concentrations of
ionic species predicted by the linear theory slightly deviate from the
numerical solution near the electrode surface ($z=0$). The deviation increases
at larger potential drops. We would like to note that eqs. (30), (32), and
(35), obtained within the linear theory, determine the asymptotic behavior of
the potential profile and local ionic concentration at sufficiently large
distances from the electrode at any applied voltages. It is interesting to
note that the oscillation effect for the concentration profiles is very small,
so that they are qualitatively the same as those predicted by Kornyshev’s
mean-field theory for regular ILs.

Figure 6: Electrostatic potential (top) and concentration (bottom) profiles
calculated within the linear analytical theory and nonlinear theory for
negative potential drops. The data are shown for $\phi_{0}=0.4$,
$b=v^{1/3}=0.5\leavevmode\nobreak\ nm$, $T=298\leavevmode\nobreak\ K$,
$\varepsilon=5$, $q=1.6\times 10^{-19}\leavevmode\nobreak\ C$. The local
concentrations are determined via $x_{p,c}(z)=(n_{p,c}(z)-n_{0})/n_{0}$.

Figure 7: Electrostatic potential (top) and concentration (bottom) profiles
numerically calculated within the nonlinear theory for sufficiently high
voltages $\psi_{0}=\pm 0.2\leavevmode\nobreak\ V$. The data are shown for
$\phi_{0}=0.4$, $b=v^{1/3}=0.5\leavevmode\nobreak\ nm$,
$T=298\leavevmode\nobreak\ K$, $\varepsilon=5$, $q=1.6\times
10^{-19}\leavevmode\nobreak\ C$.

Finally, it is interesting to discuss the behavior of the potential and
concentration profiles at a sufficiently large potential drop, at which the
linear theory does not work. Fig. 7 shows the electrostatic potential (top)
and local concentrations (bottom) of species for the voltages $\psi_{0}=\pm
0.2\leavevmode\nobreak\ V$. As is seen, in contrast to the salt-free
polyelectrolyte solution case, where the electrostatic potential or
concentration decrease (increase) is accompanied by rather strong oscillations
Châtellier and Joanny (1996), these quantities in the case of a pure PIL do
not show pronounced oscillation behavior. Instead, there are quite pronounced
maxima or minima at a certain distance from the electrode. It is also
interesting to note that the potential profiles at $\psi_{0}=\pm
0.2\leavevmode\nobreak\ V$ are almost symmetric. Such symmetry can be
explained by the fact that at a sufficiently high voltage the local volume
fraction of the monomer units (for a negatively charged electrode) or
counterions (for a positively charged electrode) near the electrode is close
to unity so that it does not matter whether the particles are tied in a chain
or freely move – the translation and conformation contributions to the entropy
are equal to zero in both cases. We would like to underline once again that
the almost symmetric form of the differential capacitance discussed above (see
Fig. 2) has the same causes.

##  V Conclusions

We have proposed a self-consistent field theory of the polymeric ionic liquid
on a charged metal electrode. Taking into account the conformation entropy of
polymerized cations within the Lifshitz theory and electrostatic and excluded
volume interactions of ionic species within the asymmetric lattice gas model
at the mean-field level of approximation, we have obtained a system of self-
consistent field equations for the local electrostatic potential and average
concentrations of monomeric units and counterions. In the linear approximation
of the self-consistent field equations, for the cases of point-like charge and
a flat infinite uniformly charged electrode immersed in polymeric ionic
liquid, we have analyzed the behavior of local ionic concentrations and
electrostatic potential. We have derived an analytical expression for the
linear differential capacitance of the electric double layer. We have shown
that in contrast to the linear capacitance predicted by Kornyshev’s self-
consistent field theory of regular ionic liquids on a charged electrode,
linear capacitance for the case of polymeric ionic liquids strongly depends on
the molecular parameters of polymerized cations: effective excluded volume of
the monomeric unit and bond length of the polymeric cation. We have analyzed
the differential capacitance behavior as a potential drop function. We have
shown that at moderate voltages slightly asymmetric bell-shaped differential
capacitance profile is lower and wider than that predicted in Kornyshev’s
mean-field theory for monomeric ionic liquids. However, at a rather high
voltage, the differential capacitance monotonically decreases in accordance
with the same power law, $C\sim|\psi_{0}|^{-0.5}$, as for monomeric ionic
liquids, predicted within Kornyshev’s analytical self-consistent field theory
and more sophisticated calculations based on the non-local density functional
theory formalism Wu  _et al._ (2011). For the case of a polymeric ionic liquid
solution on a charged electrode, we have predicted a strongly asymmetric
camel-shaped differential capacitance profile. We have obtained a rather high
peak on the differential capacitance profile in the region of negative
voltages and an unpronounced local maximum – at positive potential drops. We
have attributed the higher peak to the dramatic increase in the local volume
fraction of the monomeric units on the electrode when moving towards higher
negative voltage values, while the small peak – to a rather smooth increase in
the local volume fraction of the counterions when increasing the positive
voltage. We have analyzed the electrostatic potential and concentration
profiles for a pure polymeric ionic liquid demonstrating that in contrast to
sufficiently dilute salt-free polyelectrolyte solutions Châtellier and Joanny
(1996), these quantities have no pronounced oscillation behavior.

We believe that our findings can be of interesting for electrochemical
applications, such as supercapacitors, batteries, fuel cells, catalysis,
electrodeposition, etc.

## Acknowledgments

The development of the self-consistent field theory presented in section 2 is
supported by the grant of the President of the Russian Federation (project No.
MD-341.2021.1.3). The results presented in sections 3 and 4 have been obtained
with supporting of the Russian Science Foundation (Grant No. 21-11-00031).

## References

  * Fedorov and Kornyshev (2014) M. V. Fedorov and A. A. Kornyshev, ‘‘Ionic liquids at electrified interfaces,’’ Chemical reviews 114, 2978–3036 (2014).
  * Ohno and Ito (1998) H. Ohno and K. Ito, ‘‘Room-temperature molten salt polymers as a matrix for fast ion conduction,’’  _Chemistry Letters_ ,  , 751–752 (1998).
  * Ohno, Yoshizawa, and Ogihara (2004) H. Ohno, M. Yoshizawa, and W. Ogihara, ‘‘Development of new class of ion conductive polymers based on ionic liquids,’’ Electrochimica Acta 50, 255–261 (2004).
  * Ohno (2007) H. Ohno, ‘‘Design of ion conductive polymers based on ionic liquids,’’ Macromolecular Symposia 249-250, 551–556 (2007).
  * Yuan and Antonietti (2011) J. Yuan and M. Antonietti, ‘‘Poly(ionic liquid)s: Polymers expanding classical property profiles,’’ Polymer 52, 1469–1482 (2011).
  * Mecerreyes (2011) D. Mecerreyes, ‘‘Polymeric ionic liquids: Broadening the properties and applications of polyelectrolytes,’’ Progress in Polymer Science 36, 1629–1648 (2011).
  * Eftekhari and Saito (2017) A. Eftekhari and T. Saito, ‘‘Synthesis and properties of polymerized ionic liquids,’’ European Polymer Journal 90, 245–272 (2017).
  * (8) Y. A. Fadeeva, L. E. Shmukler, M. S. Gruzdev, and L. P. Safonova, ‘‘Imidazolium zwitterion based protic ionic liquids: from monomers to polymer membranes,’’ Polymer International https://doi.org/10.1002/pi.6238.
  * Chen  _et al._ (2018) S. Chen, F. Frenzel, B. Cui, F. Gao, A. Campanella, A. Funtan, F. Kremer, S. S. Parkin, and W. H. Binder, ‘‘Gating effects of conductive polymeric ionic liquids,’’ Journal of Materials Chemistry C 6, 8242–8250 (2018).
  * Peltekoff  _et al._ (2019) A. J. Peltekoff, V. E. Hiller, G. P. Lopinski, O. A. Melville, and B. H. Lessard, ‘‘Unipolar Polymerized Ionic Liquid Copolymers as High-Capacitance Electrolyte Gates for n-Type Transistors,’’ ACS Applied Polymer Materials 1, 3210–3221 (2019).
  * Gupta, Liang, and Hu (2019) N. Gupta, Y. N. Liang, and X. Hu, ‘‘Thermally responsive ionic liquids and polymeric ionic liquids: emerging trends and possibilities,’’ Current Opinion in Chemical Engineering 25, 43–50 (2019).
  * Yuki Kohno  _et al._ (2015) Yuki Kohno, Shohei Saita, Yongjun Men, Jiayin Yuan, and Hiroyuki Ohno, ‘‘Thermoresponsive polyelectrolytes derived from ionic liquids,’’ Polymer Chemistry 6, 2163–2178 (2015).
  * Tang  _et al._ (2014) S. Tang, S. Liu, Y. Guo, X. Liu, and S. Jiang, ‘‘Recent advances of ionic liquids and polymeric ionic liquids in capillary electrophoresis and capillary electrochromatography,’’ Journal of Chromatography A 1357, 147–157 (2014).
  * Gao  _et al._ (2021) T. Gao, J. Itliong, S. P. Kumar, Z. Hjorth, and I. Nakamura, ‘‘Polarization of ionic liquid and polymer and its implications for polymerized ionic liquids: An overview towards a new theory and simulation,’’ Journal of Polymer Science (2021).
  * Dobrynin and Rubinstein (2005) A. V. Dobrynin and M. Rubinstein, ‘‘Theory of polyelectrolytes in solutions and at surfaces,’’ Progress in Polymer Science (Oxford) 30, 1049–1118 (2005).
  * Netz and Andelman (2003) R. R. Netz and D. Andelman, ‘‘Neutral and charged polymers at interfaces,’’ Physics Reports 380, 1–95 (2003), arXiv:0203364 [cond-mat] .
  * Van Der Schee and Lyklema (1984) H. A. Van Der Schee and J. Lyklema, ‘‘A lattice theory of polyelectrolyte adsorption,’’ Journal of Physical Chemistry 88, 6661–6667 (1984).
  * Papenhuijzen, Van Der Schee, and Fleer (1985) J. Papenhuijzen, H. A. Van Der Schee, and G. J. Fleer, ‘‘Polyelectrolyte adsorption. I. A new lattice theory,’’ Journal of Colloid And Interface Science 104, 540–552 (1985).
  * Evers  _et al._ (1986) O. A. Evers, G. J. Fleer, J. M. Scheutjens, and J. Lyklema, ‘‘Adsorption of weak polyelectrolytes from aqueous solution,’’ Journal of Colloid And Interface Science 111, 446–454 (1986).
  * Roe (1974) R. J. Roe, ‘‘Multilayer theory of adsorption from a polymer solution,’’ The Journal of Chemical Physics 60, 4192–4207 (1974).
  * Scheutjens and Fleer (1979) J. Scheutjens and G. Fleer, ‘‘Statistical theory of the adsorption of interacting chain molecules. 1. partition function, segment density distribution, and adsorption isotherms,’’ Journal of Physical Chemistry 83, 1619–1635 (1979).
  * Scheutjens and Fleer (1980) J. M. Scheutjens and G. J. Fleer, ‘‘Statistical theory of the adsorption of interacting chain molecules. 2. Train, loop, and tail size distribution,’’ Journal of Physical Chemistry 84, 178–190 (1980).
  * Borisov, Zhulina, and Birshtein (1994) O. Borisov, E. Zhulina, and T. Birshtein, ‘‘Polyelectrolyte molecule conformation near a charged surface,’’ Journal de Physique II 4, 913–929 (1994).
  * Friedsam, Gaub, and Netz (2005) C. Friedsam, H. Gaub, and R. Netz, ‘‘Adsorption energies of single charged polymers,’’ EPL (Europhysics Letters) 72, 844 (2005).
  * Varoqui, Johner, and Elaissari (1991) R. Varoqui, A. Johner, and A. Elaissari, ‘‘Conformation of weakly charged polyelectrolytes at a solid–liquid interface,’’ The Journal of chemical physics 94, 6873–6878 (1991).
  * Varoqui (1993) R. Varoqui, ‘‘Structure of weakly charged polyelectrolytes at a solid-liquid interface,’’ Journal de Physique II 3, 1097–1108 (1993).
  * Borukhov, Andelman, and Orland (1995) I. Borukhov, D. Andelman, and H. Orland, ‘‘Polyelectrolyte solutions between charged surfaces,’’ Epl 32, 499–504 (1995).
  * Châtellier and Joanny (1996) X. Châtellier and J.-F. Joanny, ‘‘Adsorption of polyelectrolyte solutions on surfaces: a debye-hückel theory,’’ Journal de Physique II 6, 1669–1686 (1996).
  * Borukhov, Andelman, and Orland (1998) I. Borukhov, D. Andelman, and H. Orland, ‘‘Scaling laws of polyelectrolyte adsorption,’’ Macromolecules 31, 1665–1671 (1998).
  * Joanny (1999) J. F. Joanny, ‘‘Polyelectrolyte adsorption and charge inversion,’’ European Physical Journal B 9, 117–122 (1999).
  * Podgornik and Jönsson (1993) R. Podgornik and B. Jönsson, ‘‘Stretching of polyelectrolyte chains by oppositely charged aggregates,’’ EPL (Europhysics Letters) 24, 501 (1993).
  * Brilliantov, Budkov, and Seidel (2016) N. V. Brilliantov, Y. A. Budkov, and C. Seidel, ‘‘Generation of mechanical force by grafted polyelectrolytes in an electric field,’’ Physical Review E 93, 032505 (2016).
  * Bohmer, Evers, and Scheutjens (1990) M. R. Bohmer, O. A. Evers, and J. M. Scheutjens, ‘‘Weak polyelectrolytes between two surfaces: adsorption and stabilization,’’ Macromolecules 23, 2288–2301 (1990).
  * Li and Wu (2006a) Z. Li and J. Wu, ‘‘Density functional theory for planar electric double layers: Closing the gap between simple and polyelectrolytes,’’ Journal of Physical Chemistry B 110, 7473–7484 (2006a).
  * Li and Wu (2006b) Z. Li and J. Wu, ‘‘Density functional theory for polyelectrolytes near oppositely charged surfaces,’’ Physical Review Letters 96, 1–4 (2006b).
  * Kumar  _et al._ (2017) R. Kumar, J. P. Mahalik, V. Bocharova, E. W. Stacy, C. Gainaru, T. Saito, M. P. Gobet, S. Greenbaum, B. G. Sumpter, and A. P. Sokolov, ‘‘A Rayleighian approach for modeling kinetics of ionic transport in polymeric media,’’ Journal of Chemical Physics 146 (2017), 10.1063/1.4975309.
  * Kumar  _et al._ (2020) R. Kumar, J. P. Mahalik, K. S. Silmore, Z. Wojnarowska, A. Erwin, J. F. Ankner, A. P. Sokolov, B. G. Sumpter, and V. Bocharova, ‘‘Capacitance of thin films containing polymerized ionic liquids,’’ Science Advances 6 (2020), 10.1126/sciadv.aba7952.
  * Kornyshev (2007) A. A. Kornyshev, ‘‘Double-layer in ionic liquids: paradigm change?’’ The Journal of Physical Chemistry B 111, 5545–5557 (2007).
  * Budkov  _et al._ (2018) Y. A. Budkov, A. L. Kolesnikov, Z. A. Goodwin, M. G. Kiselev, and A. A. Kornyshev, ‘‘Theory of electrosorption of water from ionic liquids,’’ Electrochimica Acta 284, 346–354 (2018).
  * Wu  _et al._ (2011) J. Wu, T. Jiang, D. Jiang, Z. Jin, and D. Henderson, ‘‘A classical density functional theory for interfacial layering of ionic liquids,’’ Soft Matter 7, 11222–11231 (2011).
  * Jiang, Meng, and Wu (2011) D. Jiang, D. Meng, and J. Wu, ‘‘Density functional theory for differential capacitance of planar electric double layers in ionic liquids,’’ Chemical Physics Letters 504, 153–158 (2011).
  * Forsman, Woodward, and Trulsson (2011) J. Forsman, C. E. Woodward, and M. Trulsson, ‘‘A classical density functional theory of ionic liquids,’’ The Journal of Physical Chemistry B 115, 4606–4612 (2011).
  * Goodwin, Feng, and Kornyshev (2017) Z. A. Goodwin, G. Feng, and A. A. Kornyshev, ‘‘Mean-field theory of electrical double layer in ionic liquids with account of short-range correlations,’’ Electrochimica Acta 225, 190–197 (2017).
  * Maggs and Podgornik (2016) A. Maggs and R. Podgornik, ‘‘General theory of asymmetric steric interactions in electrostatic double layers,’’ Soft matter 12, 1219–1229 (2016).
  * May (2019) S. May, ‘‘Differential capacitance of the electric double layer: mean-field modeling approaches,’’ Current Opinion in Electrochemistry 13, 125–131 (2019).
  * Cruz  _et al._ (2019) C. Cruz, S. Kondrat, E. Lomba, and A. Ciach, ‘‘Effect of proximity to ionic liquid-solvent demixing on electrical double layers,’’ Journal of Molecular Liquids 294, 111368 (2019).
  * Bazant, Storey, and Kornyshev (2011) M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, ‘‘Double layer in ionic liquids: Overscreening versus crowding,’’ Physical review letters 106, 046102 (2011).
  * Avni, Adar, and Andelman (2020) Y. Avni, R. M. Adar, and D. Andelman, ‘‘Charge oscillations in ionic liquids: A microscopic cluster model,’’ Physical Review E 101, 010601 (2020).
  * Budkov, Zavarzin, and Kolesnikov (2021) Y. A. Budkov, S. V. Zavarzin, and A. L. Kolesnikov, ‘‘Theory of ionic liquids with polarizable ions on a charged electrode,’’ J. Phys. Chem. C (2021), https://doi.org/10.1021/acs.jpcc.1c05548.
  * Lifshitz (1969) I. Lifshitz, ‘‘Some problems of the statistical theory of biopolymers,’’ Sov. Phys. JETP 28, 1280–1286 (1969).
  * Grosberg and Khokhlov (1994) A. Y. Grosberg and A. R. Khokhlov,  _Statistical physics of macromolecules_ (Amer Inst of Physics, 1994).
  * Sanchez and Lacombe (1978) I. C. Sanchez and R. H. Lacombe, ‘‘Statistical thermodynamics of polymer solutions,’’ Macromolecules 11, 1145–1156 (1978).
  * Borue and Erukhimovich (1988) V. Y. Borue and I. Y. Erukhimovich, ‘‘A Statistical Theory of Weakly Charged Polyelectrolytes: Fluctuations, Equation of State and Microphase Separation,’’ Macromolecules 21, 3240–3249 (1988).
  * Muthukumar (1996) M. Muthukumar, ‘‘Double screening in polyelectrolyte solutions: Limiting laws and crossover formulas,’’ The Journal of chemical physics 105, 5183–5199 (1996).
  * Budkov (2020) Y. A. Budkov, ‘‘Statistical field theory of ion-molecular solutions,’’ Physical Chemistry Chemical Physics (2020).
  * Budkov (2019) Y. A. Budkov, ‘‘A statistical field theory of salt solutions of "hairy" dielectric particles,’’ Journal of Physics: Condensed Matter 32, 055101 (2019).
  * Slavchov (2014) R. I. Slavchov, ‘‘Quadrupole terms in the maxwell equations: Debye-hückel theory in quadrupolarizable solvent and self-salting-out of electrolytes,’’ The Journal of Chemical Physics 140, 074503 (2014).