# Plasma broadening of autoionizing resonances

A K Pradhan1,2,3 1 Department of Astronomy, 2 Chemical Physics Program, 3
Biophysics Graduate Program, Ohio State University, Columbus, Ohio 43210, USA

###### Abstract

A general formulation is developed to demonstrate that atomic autoionizing
(AI) resonances are broadened and shifted significantly due to plasma effects
across bound-free continua. The theoretical and computational method presented
accounts for broadening mechanisms: electron collisional, ion microfields
(Stark), thermal Doppler, core excitations, and free-free transitions.
Extrinsic plasma broadening redistributes and shifts AI resonance strengths
while broadly preserving naturally intrinsic asymmetries of resonance
profiles. Integrated oscillator strengths are conserved as resonance
structures dissolve into continua with increasing electron density. As
exemplar, the plasma attenuation of photoionization cross sections computed
using the R-matrix method is studied in neon-like Fe XVII in a critical range
$N_{e}=10^{21-24}$cc along isotherms $T=1-2\times 10^{6}$K, and its impact on
Rosseland Mean opacities. The energy-temperature-density dependent cross
sections would elicit and introduce physical features in resonant processes in
photoionization, (e + ion) excitation and recombination. The method should be
generally applicable to atomic species in high-energy-density (HED) sources
such as fusion plasmas and stellar interiors.

  * September 2022

##  1 Introduction

Autoionization (AI) resonances are ubiquitous in atomic reactions. As bound-
quasibound-free transitions they manifest themselves in cross sections of
atomic processes in a wider variety of shapes, sizes and extent in energy than
line formation due to bound-bound transitions. As such, due to generally
higher intrinsic decay rates relative to radiative rates, it would be expected
that AI resonances broaden, smear out, and dissolve into the bound-free
continuum far more readily than lines when subjected to extrinsic high-energy-
density (HED) plasma environments. However, whereas line broadening is an
advanced field and is accounted for via elaborate treatments, a theoretical
method for AI resonance broadening is not available. Also, although line
broadening treatments are precise for hydrogenic and simple atomic systems,
several approximations need to be made using fitting formulae, Gaunt factors,
etc., for complex atomic systems involving large numbers of transitions for
practical applications in laboratory and astrophysical plasmas [2, 3, 4, 5, 6,
7, 8, 9, 10].

Whereas the main broadening mechanisms in AI broadening are physically similar
to line broadening, their theoretical and computational treatment is quite
different. Superimposed on intrinsic AI broadening in atomic cross sections
the extent of resonances owing to extrinsic plasma effects renders much of the
line broadening theory inapplicable, particularly for multi-electron systems.
The unbroadened AI resonances themselves vary by orders of magnitude in width,
shapes and heights, and incorporate two types: large features due to
photoexcitation-of-core (PEC) below thresholds corresponding to dipole core
transitions [11], and infinite Rydberg series of resonances converging on to
each excited core level of the (e + ion) system. The generally employed Voigt
line profiles obtained by convolution of a Lorentzian function for radiative
and collisional broadening, and a Gaussian function for Doppler or thermal
broadening, are found to be practically inapplicable for AI broadening.
Numerically, the Voigt kernel is ill-conditioned since the collisional-to-
Doppler width ratio $\Gamma_{c}$ /$\Gamma_{d}$ varies over a far wider range
for resonances than lines and therefore unconstrained a priori [12, 13].

##  2 Theoretical and computational method

The physical processes for broadening of AI resonances differ from lines
qualitatively and quantitatively. However, line broadening processes and
formulae may be generalized to develop a theoretical treatment and
computational algorithm outlined herein (details to be presented elsewhere).
The convolved bound-free photoionization cross section of level $i$ may be
written as:

$\sigma_{i}(\omega)=\int\tilde{\sigma}(\omega^{\prime})\phi(\omega^{\prime},\omega)d\omega^{\prime},$
(1)

where $\sigma$ and $\tilde{\sigma}$ are the cross sections with plasma-
broadened and unbroadened AI resonance structures, $\omega$ is the photon
energy (Rydberg atomic units are used throughout), and
$\phi(\omega^{\prime},\omega)$ is the normalized Lorentzian profile factor in
terms of the total width $\Gamma$ due to all AI broadening processes included:

$\phi(\omega^{\prime},\omega)=\frac{\Gamma(\omega)/\pi}{x^{2}+\Gamma^{2}},$
(2)

where $x\equiv\omega-\omega^{\prime}$. The crucial difference with line
broadening is that AI resonances in the (e + ion) system correspond to and are
due to quantum mechanical interference between discretized continua defined by
excited core ion levels in a multitude of channels. The coupled channel (CC)
approximation, such as implemented by the R-matrix (RM) method (viz. [2, 5,
6]), accounts for AI resonances in an (e + ion) system with generally
asymmetric profiles (unlike line profiles that are usually symmetric).
Schematically, CC-RM calculations are as shown in Fig. 1 for an (e + ion)
system with bound and continuum levels.

Figure 1: Schematic diagram of coupled channel calculations: photoionization
of bound (e + ion) levels into continua of core ion levels (solid lines)
$\hbar\omega+X^{z+}\rightarrow e+X^{z+1}$. Rydberg series of AI resonances
(dashed lines) converge on to excited core levels at energies
$E=E_{i}-z^{2}/\nu_{i}^{2}$, where $\nu_{i}$ is a continuous variable
corresponding to each threshold energy $E_{i}$. The schematics approximately
refers to Fe xvii cross sections exemplified in Figs. 2 and 3, with 218
n=2,3,4 energy levels of Fe xviii .

Given $N$ core ion levels corresponding to resonance structures,

$\sigma(\omega)=\sum_{i}^{N}\left[\int\tilde{\sigma}(\omega^{\prime})\left[\frac{\Gamma_{i}(\omega)/\pi}{x^{2}+\Gamma_{i}^{(}\omega)}\right]d\omega^{\prime}\right].$
(3)

With $x\equiv\omega^{\prime}-\omega$, the summation is over all excited
thresholds $E_{i}$ included in the $N$-level CC or RM wavefunction expansion,
and corresponding to total damping width $\Gamma_{i}$ due to all broadening
processes. The profile $\phi(\omega^{\prime},\omega)$ is centered at each
continuum energy $\omega$, convolved over the variable $\omega^{\prime}$ and
relative to each excited core ion threshold $i$ . In the present formulation
we associate the energy to the effective quantum number relative to each
threshold $\omega^{\prime}\rightarrow\nu_{i}$ to write the total width as:

$\displaystyle\Gamma_{i}(\omega,\nu,T,N_{e})$ $\displaystyle=$
$\displaystyle\Gamma_{c}(i,\nu,\nu_{c})+\Gamma_{s}(\nu_{i},\nu_{s}^{*})$
$\displaystyle+$
$\displaystyle\Gamma_{d}(A,\omega)+\Gamma_{f}(f-f;\nu_{i},\nu_{i}^{\prime}),$

pertaining to collisional $\Gamma_{c}$, Stark $\Gamma_{s}$, Doppler
$\Gamma_{d}$, and free-free transition $\Gamma_{f}$ widths respectively, with
additional parameters as defined below. Without loss of generality we assume a
Lorentizan profile factor that describes collisional-ion broadening which
dominates in HED plasmas. We assume this approximation to be valid since
collisional profile wings extend much wider as $x^{-2}$, compared to the
shorter range $exp(-x^{2})$ for thermal Doppler, and $x^{-5/2}$ for Stark
broadening (viz. [14]). In Eq. (4) the limits $\mp\infty$ are then replaced by
$\mp\Gamma_{i}/\sqrt{\delta}$; $\delta$ is chosen to ensure the Lorentzian
profile energy range for accurate normalization. Convolution by evaluation of
Eqs. (1-3) is carried out for each energy $\omega$ throughout the tabulated
mesh of energies used to delineate all AI resonance structures, for each cross
section, and each core ion threshold. We employ the following expressions for
computations:

$\Gamma_{c}(i,\nu)\ =\
5\left(\frac{\pi}{kT}\right)^{1/2}a_{o}^{3}N_{e}G(T,z,\nu_{i})(\nu_{i}^{4}/z^{2}),$
(5)

where T, $N_{e}$ , $z$, and $A$ are the temperature, electron density, ion
charge and atomic weight respectively, and $\nu_{i}$ is the effective quantum
number relative to each core ion threshold $i$ : $\omega\equiv
E=E_{i}-\nu_{i}^{2}/z^{2}$ is a continuous variable. The Gaunt factor
$G(T,z,\nu_{i})=\sqrt{3}/\pi[1/2+ln(\nu_{i}kT/z)]$ ([14, 8, 23, 20]). A factor
$(n_{x}/n_{g})^{4}$ is introduced for $\Gamma_{c}$ to allow for doubly excited
AI levels with excited core evels $n_{x}$ relative to the ground configuration
$n_{g}$ (e.g. for Fe xviii $n_{x}=3,4$ relative to the ground configuration
$n_{g}=2$). A treatment of the Stark effect for complex systems entails two
approaches, one where both electron and ion perturbations are combined (viz.
[9, 23]), or separately (viz. [5, 14]) employed herein. Excited Rydberg levels
are nearly hydrogenic and ion perturbations are the main broadening effect,
though collisional broadening competes significantly increasing with density
as well as $\nu_{i}^{4}$ (Eq. 5). The total Stark width of a given $n$
-complex is $\approx(3F/z)n^{2}$, where F is the plasma electric microfield.
Assuming the dominant ion perturbers to be protons and density equal to
electrons, $N_{e}$ =$N_{p}$, we take $F=[(4/3)\pi a_{o}^{3}N_{e})]^{2/3}$, as
employed in the Mihalas-Hummer-Däppen equation-of-state formulation [15].

$\Gamma_{s}(\nu_{i},\nu_{s}^{*})=[(4/3)\pi a_{o}^{3}N_{e}]^{2/3}\nu_{i}^{2}.$
(6)

In addition, in employing Eq. (6) a Stark ionization parameter
$\nu_{s}^{*}=1.2\times 10^{3}N_{e}^{-2/15}z^{3/5}$ is introduced such that AI
resonances may be considered fully dissolved into the continuum for
$\nu_{i}>\nu_{s}^{*}$ (analogous to the Inglis-Teller series limit [16, 15]).
Calculations are carried out with and without $\nu_{s}^{*}$ as shown later in
Table 1. The Doppler width is:

$\Gamma_{d}(A,T,\omega)=4.2858\times 10^{-7}\sqrt{(}T/A),$ (7)

where $\omega$ is not the usual line center but taken to be each AI resonance
energy. The last term $\Gamma_{f}$ in Eq. (5) accounts for free-free
transitions among autoionizing levels with $\nu_{i},\nu_{i}^{\prime}$ such
that

$X_{i}+e(E_{i},\nu_{i})\longrightarrow
X_{i}^{\prime}+e^{\prime}(E_{i}^{\prime},\nu_{i}^{\prime}).$ (8)

The large number of free-free transition probabilities for $+ve$ energy AI
levels $E_{i},E_{i}^{\prime}>0$ may be computed using RM or atomic structure
codes (viz. [17, 22]).

Whereas Eq.(3) has an analytical solution in terms of
$tan^{-1}(x/\Gamma)/\Gamma$ evaluated at limiting values of
$x\rightarrow\mp\Gamma/\sqrt{\delta}$, its evaluation for practical
applications entails piece-wise integration across multiple energy ranges
spanning many excited thresholds and different boundary conditions. For
example, the total width $\Gamma$ is very large at high densities and the
Lorentzian profile may be incomplete above the ionization threshold and
therefore not properly normalized. We obtain the necessary redward left-wing
correction for partial renormalization as

$\lim_{a\rightarrow-\Gamma/2\sqrt{\delta}}\int_{a}^{+\Gamma/\sqrt{\delta}}\phi(\omega,\omega^{\prime})d\omega^{\prime}=\left[\frac{1}{4}-\frac{tan^{-1}(\frac{a}{\Gamma/2\sqrt{\delta}})}{\pi}\right],$
(9)

where $a$ is the lower energy range up to the ionization threshold, reaching
the maximum value $-\Gamma/2\sqrt{\delta}$.

##  3 Results and discussion

The complexity and magnitude of computations is demonstrated for Fe xvii that
is of considerable importance in astrophysical and laboratory plasmas
described in a number of previous works ([19] and references therein), owing
to its neon-like ground configuration and many excited configurations and
$n$-complexes of levels and transitions. We utilize new results from an
extensive Breit-Pauli R-Matrix (BPRM) calculation with 218 fine structure
levels dominated by $n=2,3,4$ levels of the core ion Fe xviii (to be reported
elsewhere). The 587 Fe xvii bound levels ($E<0$) considered are dominated by
configurations $1s^{2}2s^{2}2p^{6}(^{1}S_{0}),1s^{2}2s^{p}2p^{q}n\ell,[SLJ]\
(p,q=0-2,\ n\leq 10,\ \ell\leq 9,\ J\leq 12$). The core Fe xvii levels
included in the CC calculation for the (e + Fe xviii ) $\rightarrow$Fe xvii
system
are:$1s^{2}2s^{2}2p^{5}(^{2}P^{o}_{1/2,3/2}),1s^{2}2s^{2}2p^{q},n\ell,[S_{i}L_{i}J_{i}]\
(p=4,5,\ n\leq 4,\ell\leq 3)$. The Rydberg series of AI resonances correspond
to $(S_{i}L_{i}J_{i})\ n\ell,\ n\leq 10,\ell\leq 9$, with effective quantum
number defined as a continuous variable $\nu_{i}=z/\sqrt{(}E_{i}-E)\ (E>0)$,
throughout the energy range up to the highest $218^{th}$ Fe xviii core level;
the $n=2,3,4$ core levels range from E=0-90.7 Ry ([18, 19]). The Fe xvii BPRM
calculations were carried out resolving the bound-free cross sections at
$\sim$40,000 energies for 454 bound levels with AI resonance structures. Given
217 excited core levels of Fe xviii , convolution is carried out at each
energy or approximately $10^{9}$ times for each (T,$N_{e}$) pair.

Figure 2: Plasma broadened photoionization cross sections for
$\hbar\omega+{\rm Fe~{}\sc xvii}\ \rightarrow e~{}+~{}{\rm Fe~{}\sc xviii}\ $
of the bound level $2s^{2}2p^{5}[^{2}P^{o}_{3/2}]3p(^{3}D_{2})$ (ionization
energy 37.707 Ry) along isotherm $T=2\times 10^{6}$K and electron densities
$N_{e}=10^{21,22,23}$cc: black — unbroadened, red — broadened, blue —
broadened with Stark ionization cut-off $\nu_{s}^{*}$ (Table 1). Rydberg
series of AI resonance complexes with $\nu_{i}\leq 10$ belonging to 217
excited Fe xviii levels $E_{i}$ broaden and shift with increasing density,
also resulting in continuum raising and threshold lowering.  Figure 3: Plasma
broadened photoionization cross sections of Fe xvii level
$2s^{2}2p^{5}[^{2}P^{o}_{3/2}]4d(^{1}F^{o}_{3})$ (ionization energy 17.626
Ry), along isotherm $T=10^{6}$K and electron densities
$N_{e}=10^{21.5,22.5,23.5}$cc, as in Fig. 2.

Fig. 2 displays detailed results for plasma broadened and unbroadened
photoionization cross section of one particular excited level
$2s^{2}2p^{5}[^{2}P^{o}_{3/2}]4d(^{1}F^{o}_{3})$ (ionization energy = 17.626
Ry) of Fe xvii at three representative densities (note the $\sim$10 orders of
magnitude variation in resonance heights along the Y-axis). The main feature
evident in the figure are as follows. (i) AI resonances begin to show
significant broadening and smearing of a multitude of overlapping Rydberg
series at $N_{e}=10^{21}$cc. The narrower high-$n$ $l$ resonances dissolve
into the continua but stronger low-$n$ $l$ resonance retain their asymmetric
shapes with attenuated heights and widths. (ii) As the density increases by
one to two order of magnitude, to $N_{e}=10^{22-23}$cc, resonance structures
not only broaden but their strengths shift and redistributed over a wide range
determined by total width $\Gamma(\omega,\nu_{i},T,N_{e})$ at each energy
$\hbar\omega$ (Eq. 4). (iii) Stark ionization cut-off (Table 1) results in
step-wise structures that represent the average due to complete dissolution
into continua. (iv) The total AI resonance strengths are conserved, and
integrated values generally do not deviate by more than 1-2%. For example, the
three cases in Fig. 2: unbroadened structure (black), and broadened without
(red) and with Stark cut-off (blue), the integrated numerical values are
59.11, 59.96, 59.94 respectively. This is also an important accuracy check on
numerical integration and the computational algorithm, as well as the choice
of the parameter $\delta$ that determines the energy range of the Lorentizan
profile at each T and $N_{e}$; in the present calculations it varies from
$\delta$ = 0.01-0.05 for $N_{e}$ =$10^{21-24}$cc.

Fig. 3 shows similar results ito Fig. 2 for another excited Fe xvii level
$2s^{2}2p^{5}[^{2}P^{o}_{3/2}]4d(^{1}F^{o}_{3})$ (ionization energy 17.626
Ry), along a lower temperature $10^{6}$K isotherm at different intermediate
densities. Both Figs. 2 and 3 show a redward shift of low-$n$ resonances and
dissolution of high-$n$ resonances. In addition, the background continuum is
raised owing to redistribution of resonance strengths, which merge into one
across high ying and overlapping thresholds. Free-free transitions are not
considered in the results in Figs. 2 and 3 but included in the resutls
discussed in Table 1, although it is found to be practically negligible.

Table 1 gives plasma parameters corresponding to Figs. 2 and 3. Their physical
significance is demonstrated by a representative sample tabulated temperature
T(K) and $N_{e}$ . The maximum width $\Gamma_{10}$ corresponding to
$\nu_{i}=10$ in Eqs. (4-7) is set by the CC-BPRM calculations which delineate
unbroadened AI resonance profiles up to $\nu\leq 10$, and employ an averaging
procedure up to each threshold $10<\nu<\infty$ using quantum defect (QD)
theory (viz. [21, 5, 6] and references therein). $\Gamma_{c}(10)$ and
$\Gamma_{s}(10)$ are the maximum collisional and Stark width components. The
Doppler width $\Gamma_{d}$ is much smaller, $1.18\times 10^{-3}$ and
$1.67\times 10^{-3}$ Ry at $10^{6}$K and $2\times 10^{6}$K respectively,
validating its inclusion in Eq. (4) in HED plasma sources but possibly not
when $\Gamma_{d}$ is comparable to $\Gamma_{c}$ or $\Gamma_{s}$. The
$\nu^{*}_{s}$ and $\nu_{D}$ are effective quantum numbers corresponding to
Stark ionization cut-off and the Debye radius respectively. We obtain
$\nu_{D}=\left[\frac{2}{5}\pi z^{2}\lambda_{D}^{2}\right]^{1/4}$, where the
Debye length $\lambda_{D}=(kT/8\pi N_{e})^{1/2}$. It is seen in Table 1 that
$\nu_{D}>\nu^{*}_{s}$ at the T, $N_{e}$ considered, justifying neglect of
plasma screening effects herein, but which may need to be accounted for at
even higher densities.

The aggregate effect of AI broadening for large-scale applications is
demonstrated in Table 1 by the ratio R of the Rosseland Mean Opacity using
broadened/unbroadened cross sections for 454 Fe xvii levels with AI resonances
(other higher bound levels have negligible resonances) [24, 19]. For any atom
or ion R is highly dependent on T and $N_{e}$ ; for Fe xvii R yields up to 58%
enhancement due to plasma broadening with increasing $N_{e}$ along the
$2\times 10^{6}$K isotherm, but decreasing to 6% along the $10^{6}$K isotherm.
Approximately 70,000 free-free transitions among +ve energy levels are
included in the calculation of R, but their contribution has no significant
broadening effect since they entail very high-lying levels with negligible
level populations. However, different plasma environments with intense
radiation fields, or a different equation-of-state than [15] employed here,
may lead to more discernible effect due to free-free transitions. AI
broadening in a plasma environment affects each level cross section
differently, and hence its contribution to opacities or rate equations for
atomic processes in general. A critical (T,$N_{e}$ ) range can therefore be
numerically ascertained where redistribution and shifts of atomic resonance
strengths would be significant and cross sections should be modified.

Table 1: Plasma parameters along isotherms in Fig. 2 and 3; $\nu_{D}$ corresponds to Debye radius; R is the ratio of Fe xvii Rosseland Mean Opacity with and without broadening [24]; $\Gamma_{10}$ is the maximum AI resonance width at $\nu=10$. T(K) | $N_{e}(cc)$ | $\Gamma_{10}(Ry)$ | $\Gamma_{c}(10)$ | $\Gamma_{s}(10)$ | $\nu_{s}^{*}$ | $\nu_{D}$ | R
---|---|---|---|---|---|---|---
|  | $\nu=10$ |  |  |  |  |
hline $2\times 10^{6}$ | $10^{21}$ | 3.42(-1) | 8.55(-2) | 2.57(-1) | 10.4 | 28.1 | 1.35
$2\times 10^{6}$ | $10^{22}$ | 2.05(0) | 8.55(-1) | 1.19(0) | 7.7 | 15.8 | 1.43
$2\times 10^{6}$ | $10^{23}$ | 1.41(1) | 8.55(0) | 5.53(0) | 5.6 | 8.9 | 1.55
$2\times 10^{6}$ | $10^{24}$ | 1.11(2) | 8.55(1) | 2.57(1) | 4.1 | 5.0 | 1.58
$10^{6}$ | $3.1\times 10^{21.5}$ | 8.17(-1) | 2.71(-1) | 5.46(-1) | 9.0 | 17.8 | 1.47
$10^{6}$ | $3.1\times 10^{22.5}$ | 5.25(0) | 2.71(0) | 2.53(0) | 6.6 | 10.0 | 1.13
$10^{6}$ | $3.1\times 10^{23.5}$ | 3.89(0) | 2.71(1) | 1.18(0) | 4.8 | 5.6 | 1.06

##  4 Conclusion

The main conclusions are: (I) The method described herein is generally
applicable to AI resonances in atomic processes in HED plasmas. (II) The cross
sections become energy-temperature-density dependent in a The cross sections
become energy-temperature-density dependent in a critical range leading to
broadening, shifting, and dissolving into continua. (III) Among the
approximations necessary to generalize the formalism is the assumption that
thermal Doppler widths are small compared to collisional and Stark widths as
herein, but given the intrinsic asymmetries of AI resonances it may not lead
to significant inaccuracies (although that needs to be verified in future
works). (IV) The treatment of Stark broadening and ionization cut-off is ad
hoc, albeit based on the equation-of-state formulation [15] and consistent
with previous works [5]. (V) Since it is negligiblhy small, the free-free
contribution is included post-facto in the computation of the ratio R in Table
1 and not in the cross sections and results shown in Figs. 2 and 3, but may be
important in special HED environments with intense radiation and should then
be incorporated in the main calculations of total AI width (Eq. 4). (VI) The
predicted redward shift of AI resonances as the plasma density increases
should be experimentally verifiable. (VII) Redistribution of AI resonance
strengths should particularly manifest itself in rate coeffcients for (e +
ion) excitiaton and recombination in plasma models and simulations, and for
photoabsorption in opacity calculations, using temperature-dependent
Maxwellian, Planck, or other particle distribution functions. (VIII) The
treatment of individual contributions to AI broadening may be improved, and
the theoretical formulation outlined here is predicated on the assumption that
external plasma effects are perturbations subsumed by and overlying the
intrinsic autoionization effect. (IX) The computational formalism is designed
to be amenable for practical applications and the computational algorithm and
a general-purpose program AUTOBRO are optimized for large-scale computations
of AI broadened atomic cross sections in HED plasma models. The CPU time
required depends mainly on the density which determines the total width
$\Gamma$; for example, in the reported calculations for Fe xvii at T=$2\times
10^{6}$K it is few minutes for $N_{e}$ = $10^{21}$ cc and $\sim$3 hours for
$N_{e}$ =$10^{24}$cc.

Acknowledgements

I would like to thank Sultana Nahar for Fe xvii atomic data and discussions.

References

## References

  * [1]
  * [2] Burke P G 2011 R-Matrix Theory of Atomic Collisions, Springer
  * [3] Drake R P 2006 High Energy Density Physics, Springer
  * [4] Hubeny Ivan and Mihalas Dimitri (2015) Theory of Stellar Atmospheres, Princeton University Press
  * [5] The Opacity Project Team, 1995 The Opacity Project, Institute of Physics Publishing, Vol 1
  * [6] Pradhan A K and Nahar S N 2011 Atomic Astrophysics and Spectroscopy Cambridge University Press
  * [7] Peach G 1981 Advances in Physics 30 367
  * [8] Dimitrijevic M S and Konjevic N 1981 Astron. Astrophys. 102 93
  * [9] Dimitrijevic M S and Konjevic N 1987 Astron. Astrophys. 172 345
  * [10] Konjevic N 1999 Phys. Repts. 316 339
  * [11] Yu Y and Seaton M J J. Phys. B 20 6409
  * [12] Hummer D 1965 Mon. Not. R. astr. Soc. 70 1
  * [13] Several formulations were employed to implement a Voigt profile but were numerically unstable owing to $\Gamma_{c}>>\Gamma_{d}$.
  * [14] Seaton M J 1990 J. Phys. B 23 3255; see also other papers in [5]
  * [15] Mihalas D, Hummer D G and Däppen W 1988 Astrophys. J. 331, 815
  * [16] Inglis D R and Teller E 1939 Astrophys. J. 90 439
  * [17] Seaton M J 2000 J. Phys. B 33 2677
  * [18] Nahar S N, Pradhan A K, Chen G.-X and Eissner W 2011 Phys. Rev. A , 83, 053417
  * [19] Nahar S N and Pradhan A K 2016 Phys. Rev. Lett. 116 235003; Ibid. Phys. Rev. Lett. , 117, 249502
  * [20] The Gaunt factor $G(T,z,\nu_{i}$ depends on T and z, and differs from the similar factor $G(z)=(8\pi/3)\sqrt{(}3)\overline{g}$, and approaches the fitting formula G(z)=6.3-5.9/(1+z) as $\nu_{i}\rightarrow 10$ used in [14, 5].
  * [21] The QD regions $10<\nu_{i}<\infty$ are generally small and overlapping for closely spaced excited core levels, as herein.
  * [22] Eissner W, Jones M and Nussbaumer H 1974 Comput. Phys. Commun. 8 270
  * [23] Xiang W, Zeng J, Fu Y and Gao C 2012 J. Mod. Phys., 3, 1670
  * [24] We eschew a discussion of opacities per se, as they also depend on accuracy and completeness of atomic data, equation-of-state of the plasma source and other quantities, and focus on AI resonance broadening in atomic processes.


# Theory of electrolyte solutions in a slit charged pore: effects of
structural interactions and specific adsorption of ions

Victoria A. Vasileva  School of Applied Mathematics, HSE University,
Tallinskaya st. 34, 123458 Moscow, Russia  Daria A. Mazur  School of Applied
Mathematics, HSE University, Tallinskaya st. 34, 123458 Moscow, Russia  Yury
A. Budkov<EMAIL_ADDRESS>School of Applied Mathematics, HSE University,
Tallinskaya st. 34, 123458 Moscow, Russia  Laboratory of Computational
Physics, HSE University, Tallinskaya st. 34, 123458 Moscow, Russia  Laboratory
of Multiscale Modeling of Molecular Systems,G.A. Krestov Institute of Solution
Chemistry of the Russian Academy of Sciences, 153045, Akademicheskaya st. 1,
Ivanovo, Russia

###### Abstract

In this paper, we present a continuation of our research on modeling
electrolyte solutions within charged slit pores. We make use of the model
developed by Blossey et al., which takes into account the structural
interactions between ions through a bilinear form over the gradients of local
ionic concentrations in the grand thermodynamic potential, as well as their
steric interactions through the lattice gas model. The structural interactions
may describe effects of the molecular structure of ions at a phenomenological
level. For example, these effects include steric effects due to non-spherical
shapes of ions, their conformation lability, and solvent effects. In addition,
we explore their specific interactions with the pore walls by incorporating
external attractive potentials. Our primary focus is on observing the behavior
of ionic concentration profiles and the disjoining pressure as the pore width
changes. By starting with the local mechanical equilibrium condition, we
derive a general expression for the disjoining pressure. Our findings indicate
that considering the structural interactions of ions leads to a pronounced
minimum on the disjoining pressure profiles at small pore widths. We attribute
this minimum to the formation of electric double layers on the electrified
surfaces of the pore. Additionally, our results demonstrate that inclusion of
the attractive interactions of ions with the pore walls enhances this minimum
and shifts it to smaller pore thicknesses. Our theoretical discoveries may be
useful for those involved in supercapacitor electrochemical engineering,
particularly when working with porous electrodes that have been infused with
concentrated electrolyte solutions.

##  I Introduction

The study of electrolyte solutions (ES) in confined geometries, such as
charged pores or slit-like channels, has gained significant importance due to
its inherent involvement in various scientific and technological applications,
ranging from energy storage devices [1, 2, 3], to water purification systems
[4, 5, 6, 7], and even to biological systems, including ion transport in
cellular and ion-exchanged membranes [8]. One of the critical challenges in
understanding these systems is the accurate modeling of structural
interactions of ions within the confinement, as well as the coupling between
ion-specific effects and the electrostatic interactions [9, 10, 11, 12, 13,
14, 15, 16]. To address this challenge, the self-consistent field theory
(SCFT) has emerged as a powerful and versatile tool in modeling the behavior
of ionic liquid-phase systems [17, 18, 19].

To effectively study electrolyte solutions confined within solid
nanostructures of varying shapes (e.g. nano-sized pores), it is essential to
calculate the mechanical stress using the stress tensor, along with ionic
concentrations and electrostatic potential profiles. Calculating the local
stress tensor, corresponding to certain SCF equations, we can determine
practically important physical properties like solvation pressure and shear
stresses [20, 21, 22, 23, 24, 25, 26]. Such properties can help estimate the
porous material deformation of a particular elastic modulus, which is crucial
for batteries and supercapacitors that use microporous electrodes impregnated
with liquid-phase electrolytes (see, for instance, ref. [27, 28, 29, 30, 25,
26]). Additionally, the stress tensor can be used to measure the macroscopic
force exerted on charged macroscopic conductors or dielectrics that are
immersed in ionic liquids. As such, a first-principles approach that enables
us to extract the stress tensor of inhomogeneous ionic liquids from the
thermodynamic potential would be valuable for practical purposes.

Some progress have been recently achieved in this area [31, 32, 33]. In ref.
[31], Budkov and Kolesnikov attempted to apply Noether’s theorem to the grand
thermodynamic potential of an ionic liquid as a functional of the
electrostatic potential. The authors established a conservation law,
$\partial_{i}\sigma_{ik}=0$, which represents the local mechanical equilibrium
condition in terms of the symmetric stress tensor, $\sigma_{ik}$. The stress
tensor obtained consisted of two components: the electrostatic (Maxwell)
stress tensor and the hydrostatic isotropic stress tensor. The former is
associated with the local electric field, while the latter, in general case
[31, 34, 35], is determined by both the local field and potential. The authors
generalized the local mechanical equilibrium condition for the cases when
external potential forces act on the ions. They then derived a general
analytical expression for the electrostatic disjoining pressure of an ionic
liquid confined in a charged nanopore slit, which extended the well-known DLVO
expression to different reference models of liquid. Budkov and Kalikin
presented a SCFT of macroscopic forces in inhomogeneous flexible chain
polyelectrolyte solutions in ref.[32]. The authors derived an analytical
expression for a stress tensor by subjecting the system to a small dilation
and considering self-consistent field equations resulting from the extremum
condition of the grand thermodynamic potential. This stress tensor, in
addition to the previously mentioned hydrostatic and Maxwell stress tensors,
also includes a conformational stress tensor generated by the conformational
(Lifshitz) entropy of flexible polymer chains. The authors applied their
theory to the investigation of polyelectrolyte solution constrained in a
conducting slit nanopore and noted anomalies in disjoining pressure and
electric differential capacitance at small pore thicknesses. Brandyshev and
Budkov [33] proposed the general covariant approach based on Noether’s second
theorem allowing them to derive the symmetric stress tensor from a grand
thermodynamic potential for an arbitrary model of inhomogeneous liquid. They
applied their approach to several models of inhomogeneous ionic liquids that
consider electrostatic correlations of ions or short-range correlations
related to packing effects. Specifically, they derived analytical expressions
for the symmetric stress tensors of the Cahn-Hilliard-like model [36], Bazant-
Storey-Kornyshev model [37], and Maggs-Podgornik-Blossey model [36].

In this paper, we model liquid electrolytes in charged pores, utilizing the
Cahn-Hilliard-like model [38, 18]. This model considers the structural and
steric interactions in the grand thermodynamic potential. Theory also takes
into account the short-range specific interactions between the ions and the
pore walls through attractive external potentials (specific adsorption).
Within this approach, we investigate the behavior of the disjoining pressure
in connection with the local ion concentrations.

##  II Theoretical background

We consider an electrolyte solution model with account for the so-called
structural interactions [38] within the Cahn-Hilliard [39] approach, which
takes into account molecular structure effects via the bilinear form over the
concentration gradients [39, 38]. Thus, the grand thermodynamic potential
(GTP) of the ES has the form

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

where

$F_{el}=\int
d\mathbb{r}\left(-\frac{\epsilon\epsilon_{0}\left(\nabla\psi\right)^{2}}{2}+\rho\psi\right)$
(2)

is the electrostatic free energy of the ES with the local electrostatic
potential $\psi(\mathbb{r})$; the local charge density of ions is
$\rho(\mathbb{r})=\sum_{\alpha}q_{\alpha}c_{\alpha}(\mathbb{r})$, $q_{\alpha}$
is the electrostatic charge of ion of the $\alpha$th type; $\epsilon$ is the
permittivity of solvent; we will model the latter as a uniform dielectric
medium; $\epsilon_{0}$ is the vacuum permittivity; $c_{\alpha}(\mathbb{r})$
are the local ionic concentrations. The GTP of the reference liquid system is

$\Omega_{liq}=\int
d\mathbb{r}\left(f(\\{c_{\alpha}\\})+\sum\limits_{\alpha}c_{\alpha}w_{\alpha}+\frac{1}{2}\sum\limits_{\alpha\gamma}\varkappa_{\alpha\gamma}\nabla
c_{\alpha}\cdot\nabla
c_{\gamma}-\sum\limits_{\alpha}\mu_{\alpha}c_{\alpha}\right),$ (3)

where the first term in the integrand is the free energy density of liquid as
the function of the local ionic concentrations, $c_{\alpha}$; the second term
is the potential energy density of ions in the external potential fields with
potential energies $w_{\alpha}(\mathbb{r})$; the third term is the
contribution of the so-called structural interactions of the ions within the
Cahn-Hilliard approach with the structural constants
$\varkappa_{\alpha\gamma}$ (see ref. [38]), $\mu_{\alpha}$ are the bulk ionic
chemical potentials. Note that in the general case, the structural
interactions may describe effects related to the molecular structure of ions
at a phenomenological level. Some examples of these effects include steric
effects due to non-spherical shapes of ions, their conformation lability, and
solvent effects. We do not specify here the physical nature of the structural
contribution to the GTP, as well as the nature of the external potentials,
wondering how it influences the local concentration profiles and disjoining
pressure (see below) relative to the regular modified Poisson-Boltzmann theory
[40, 41, 42].

Thus, the GTP can be rewritten as follows

$\Omega=\int d\mathbb{r}\omega(\mathbb{r}),$ (4)

where the GTP density is

$\omega=-\frac{\epsilon\epsilon_{0}\left(\nabla\psi\right)^{2}}{2}+\rho\psi+\frac{1}{2}\sum\limits_{\alpha\gamma}\varkappa_{\alpha\gamma}\nabla
c_{\alpha}\cdot\nabla
c_{\gamma}+f(\\{c_{\alpha}\\})+\sum\limits_{\alpha}c_{\alpha}w_{\alpha}-\sum\limits_{\alpha}\mu_{\alpha}c_{\alpha}.$
(5)

Note that summations in (5) are performed over all kinds of the ions in the
ES. Note also that integration in (4) is performed over the volume of the ES
bounded by the surfaces of conducting or dielectric macroscopic bodies. The
Euler-Lagrange equations,

$\frac{\delta\Omega}{\delta\psi}=\frac{\partial\omega}{\partial{\psi}}-\partial_{i}\left(\frac{\partial\omega}{\partial\psi_{,i}}\right)=0,\leavevmode\nobreak\
\frac{\delta\Omega}{\delta
c_{\alpha}}=\frac{\partial\omega}{\partial{c_{\alpha}}}-\partial_{i}\left(\frac{\partial\omega}{\partial
c_{\alpha,i}}\right)=0,$ (6)

yield

$\nabla^{2}\psi=-\frac{1}{\epsilon\epsilon_{0}}\sum\limits_{\alpha}q_{\alpha}c_{\alpha},\leavevmode\nobreak\
\sum\limits_{\gamma}\varkappa_{\alpha\gamma}\nabla^{2}c_{\gamma}=v_{\alpha},$
(7)

where $v_{\alpha}=q_{\alpha}\psi+w_{\alpha}+\bar{\mu}_{\alpha}-\mu_{\alpha}$,
where $\bar{\mu}_{\alpha}=\partial{f}/\partial{c_{\alpha}}$ are the intrinsic
chemical potential of the ions in the local density approximation. The first
equation is the standard Poisson equation for the electrostatic potential,
whereas the second one describes the thermodynamic equilibrium condition for
ions of the $\alpha^{th}$ kind. The system of these self-consistent field
equations can be solved with the appropriate boundary conditions for the
electrostatic potential and local ionic concentrations. Note that the boundary
conditions for the Poisson equation are determined by the nature of the
macroscopic bodies immersed in the ionic liquid. In what follows, we will
consider the case of conducting body with fixed surface potential, immersed
into the ES. For the local ionic concentrations we assume that at the surface
of the immersed bodies $c_{\alpha}=0$. The latter boundary condition is
related to the fact that the ions undergo strong repulsive forces in close
proximity to the surface. We assume that the specific attractive interactions
of the ions with the immersed bodies (specific adsorption) is included in the
potential energies $w_{\alpha}(\mathbb{r})$. Note that the short-range
specific interactions between ions [12, 43] are beyond the scope of this work.
We also do not explicitly consider the polar or polarizable solvent molecules
[44, 43, 31, 45, 46]. Furthermore, we neglect the orientation and static
polarizability of the ions [11, 10, 15, 14, 31]. We do not discuss the effect
of dielectric decrement [9, 47] and electrostatic correlations [37, 21]
either. Although these effects can be directly incorporated into the current
theory, they are irrelevant to the physical effects that are being discussed
below.

Then, with using eq. (7), we get

$\frac{\partial\omega}{\partial
x_{i}}=\frac{\partial}{\partial{x_{k}}}\left(\psi_{,i}\frac{\partial{\omega}}{\partial{\psi_{,k}}}\right)+\frac{\partial}{\partial
x_{k}}\left(\sum\limits_{\alpha}c_{\alpha,i}\frac{\partial{\omega}}{\partial
c_{\alpha,k}}\right)+\sum\limits_{\alpha}c_{\alpha}\frac{\partial{w}_{\alpha}}{\partial{x}_{i}}.$
(8)

or

$\frac{\partial\sigma_{ik}}{\partial{x_{k}}}-\sum\limits_{\alpha}c_{\alpha}\frac{\partial
w_{\alpha}}{\partial x_{i}}=0,$ (9)

where $c_{\alpha,i}=\partial_{i}c_{\alpha}$, $\psi_{,i}=\partial_{i}{\psi}$,
$\partial_{i}=\partial/\partial{x_{i}}$ and

$\sigma_{ik}=\omega\delta_{ik}-\psi_{,i}\frac{\partial{\omega}}{\partial{\psi_{,k}}}-\sum\limits_{\alpha}c_{\alpha,i}\frac{\partial{\omega}}{\partial
c_{\alpha,k}}$ (10)

is the total stress tensor [33, 31, 32]. Note that a summation over the
repeated coordinate indices in (8) and (9) is implied. With the use of the SCF
equations (7), the total stress tensor can be divided into three terms

$\sigma_{ik}=\sigma_{ik}^{(M)}+\sigma_{ik}^{(h)}+\sigma_{ik}^{(s)},$ (11)

where

$\sigma_{ik}^{(M)}=\epsilon\epsilon_{0}\bigg{(}E_{i}E_{k}-\frac{1}{2}\mathbb{E}^{2}\delta_{ik}\bigg{)}$
(12)

is the Maxwell electrostatic stress tensor with the electric field components,
$E_{i}=-\partial_{i}{\psi}$,

$\sigma_{ik}^{(h)}=-P\delta_{ik}$ (13)

is the hydrostatic stress tensor with the local osmotic pressure,
$P=\sum_{\alpha}c_{\alpha}\partial{f}/\partial{c_{\alpha}}-f$, of the ions,
whereas

$\sigma_{ik}^{(s)}=\sum\limits_{\alpha\gamma}\varkappa_{\alpha\gamma}\left[\left(\frac{1}{2}\nabla
c_{\alpha}\cdot\nabla
c_{\gamma}+c_{\alpha}\nabla^{2}c_{\gamma}\right)\delta_{ik}-\partial_{i}c_{\alpha}\partial_{k}c_{\gamma}\right]$
(14)

is the contribution of the structural interactions to the total stress tensor.
Note that the stress tensor (14) has been recently obtained by Brandyshev and
Budkov within a general covariant approach [33].

Eq. (9) is nothing but a local mechanical equilibrium condition of the ES
which will be used in the next section to derive the expression for the
disjoining pressure of the ES in a slit-like pore.

##  III Disjoining pressure

Now we would like to discuss how we can derive the general mean-field
approximation of the disjoining pressure for an ES in a slit charged nanopore.
Placing the origin of $z$ axis on one charged wall and another one at $z=H$,
we can calculate the disjoining pressure by the standard relation [24]

$\Pi=-\frac{\partial(\Omega/\mathcal{A})}{\partial H}-P_{b},$ (15)

where $P_{b}$ is the ES osmotic pressure in the bulk and

$\Omega/\mathcal{A}=\int\limits_{0}^{H}\omega(z)dz$ (16)

is the GTP per unit area of the pore walls, $\mathcal{A}$ is the total area of
the pore walls. Thus, we have

$\frac{1}{\mathcal{A}}\frac{\partial\Omega}{\partial
H}=\int\limits_{0}^{H}dz\left(\frac{\delta\Omega}{\delta\psi(z)}\frac{\partial\psi(z)}{\partial
H}+\sum\limits_{\alpha}\frac{\delta\Omega}{\delta c_{\alpha}(z)}\frac{\partial
c_{\alpha}(z)}{\partial
H}\right)+\sigma_{zz}(H)+\int\limits_{0}^{H}dz\sum\limits_{\alpha}\frac{\partial\omega}{\partial
w_{\alpha}}\frac{\partial w_{\alpha}}{\partial H}$ (17)
$=\int\limits_{0}^{H}dz\sum\limits_{\alpha}{c}_{\alpha}(z)\frac{\partial
w_{\alpha}(z)}{\partial H}+\sigma_{zz}(H),$ (18)

where we took into account that ${\partial\omega}/{\partial
w_{\alpha}}=c_{\alpha}$ and ${\delta\Omega}/{\delta\psi(z)}=0$,
${\delta\Omega}/{\delta c_{\alpha}(z)}=0$. Therefore, we eventually obtain

$\Pi=-\sum\limits_{\alpha}\int\limits_{0}^{H}dzc_{\alpha}(z)\frac{\partial
w_{\alpha}(z)}{\partial H}-\sigma_{zz}(H)-P_{b},$ (19)

where
$\sigma_{zz}(H)=(\omega-\psi^{\prime}{\partial\omega}/{\partial\psi^{\prime}}-\sum_{\alpha}c_{\alpha}^{\prime}{\partial{\omega}}/{\partial
c_{\alpha}^{\prime}})|_{z=H}$ is the normal stress at $z=H$.

Now let us consider the special practically important case, when the external
potentials are created by identical walls, i.e.

$w_{\alpha}(z)=u_{\alpha}(z)+u_{\alpha}(H-z),$ (20)

where $u_{\alpha}$ is the single wall potential. Then, taking into account
that ${c}_{\alpha}(z)={c}_{\alpha}(H-z)$, we arrive at

$\Pi=-\sum\limits_{\alpha}\int\limits_{0}^{H}dz{c}_{\alpha}(z)u_{\alpha}^{\prime}(z)-\sigma_{zz}(H)-P_{b}.$
(21)

Eq. (21) can be rewritten in a form that is more useful for practical
calculations. Using the local mechanical equilibrium condition

$\frac{d\sigma_{zz}(z)}{dz}-\sum\limits_{\alpha}{c}_{\alpha}(z)w_{\alpha}^{\prime}(z)=0,$
(22)

after the integration from $z=H/2$ to $z=H$ we obtain

$\sigma_{zz}(H)=\sigma_{zz}\left(\frac{H}{2}\right)+\sum\limits_{\alpha}\int\limits_{H/2}^{H}dz{c}_{\alpha}(z)u_{\alpha}^{\prime}(z)-\sum\limits_{\alpha}\int\limits_{0}^{H/2}dz{c}_{\alpha}(z)u_{\alpha}^{\prime}(z),$
(23)

where we accounted for eq. (20) and mentioned equality
$c_{\alpha}(z)=c_{\alpha}(H-z)$. Substituting expression (23) for
$\sigma_{zz}(H)$ in eq.(21), after some algebra, we obtain

$\Pi=P_{n}-P_{b}-2\sum\limits_{\alpha}\int\limits_{H/2}^{H}dzc_{\alpha}(z)u_{\alpha}^{\prime}(z),$
(24)

where

$P_{n}=-\sigma_{zz}\left(\frac{H}{2}\right)=P_{m}-\sum\limits_{\alpha\gamma}\varkappa_{\alpha\gamma}c_{\alpha}\left(\frac{H}{2}\right)c_{\gamma}^{\prime\prime}\left(\frac{H}{2}\right)$
(25)

is the normal osmotic pressure in the pore middle, where we took into account
that $\psi^{\prime}(H/2)=c_{\alpha}^{\prime}(H/2)=0$ and introduced the
osmotic pressure of the ions in the midpoint of the pore,
$P_{m}=P\left(H/2\right)$. Using the SCF equations (7), eq. (25) can be
simplified to

$P_{n}=P_{m}-\sum\limits_{\alpha}c_{\alpha,m}v_{\alpha,m},$ (26)

where $c_{\alpha,m}=c_{\alpha}\left(H/2\right)$,
$v_{\alpha,m}=v_{\alpha}\left(H/2\right)$. The second term in the right hand
side of eq. (26) describes contribution of the structural interactions to the
normal osmotic pressure of the ions in the slit pore with the identical
charged walls.

The main analytical result of this work, eq. (24) alongside with eq.(26), is a
generalization of the previously obtained disjoining pressure within the
modified Poisson-Boltzmann theory [31] for the case of an account of the
structural interactions of the ions. To calculate the disjoining pressure, it
is necessary to solve the SCF equations (7) with the appropriate boundary
conditions.

##  IV Numerical results and discussion

###  IV.1 Basic model

Before numerical calculations, we would like to specify the model of the ES
confined in a slit electrified pore with identical walls possessing the fixed
surface electrostatic potential, $\psi_{0}$. Following Maggs and Podgornik
[36], we use the asymmetric lattice gas model which allows us to account for
the ionic size asymmetry. The free energy density is [36]

$f=\frac{k_{B}T}{v}\left(\dfrac{\phi_{+}}{N_{+}}\ln\phi_{+}+\dfrac{\phi_{-}}{N_{-}}\ln\phi_{-}+(1-\phi_{+}-\phi_{-})\ln(1-\phi_{+}-\phi_{-})\right),$
(27)

where $v$ is a cell volume of the lattice gas, $\phi_{\pm}=N_{\pm}c_{\pm}v$
are the local volume fractions of the ions, $N_{\alpha}$ – number of cells
occupied by ion of the $\alpha$th kind ($\alpha=\pm$); $k_{B}$ is the
Boltzmann constant, $T$ is the temperature. The intrinsic chemical potentials
are

$\bar{\mu}_{\pm}=\dfrac{\partial f}{\partial
c_{\pm}}={k_{B}T}\left(\ln{\phi_{\pm}}+1-N_{\pm}\ln\left(1-\phi_{+}-\phi_{-}\right)-N_{\pm}\right),$
(28)

and the bulk chemical potentials

$\mu_{\pm}=k_{B}T\left(\ln\phi_{\pm}^{(b)}+1-N_{\pm}(\ln{(1-\phi_{+}^{(b)}-\phi_{-}^{(b)})}+1)\right),$
(29)

where $\phi^{(b)}_{\pm}=N_{\pm}c^{(b)}_{\pm}v$ are the bulk volume fractions
of the ions. Assuming that $q_{\pm}=\pm q$, we obtain from the bulk
electroneutrality condition, $q_{+}c^{(b)}_{+}+q_{-}c^{(b)}_{-}=0$, that
$c^{(b)}_{\pm}=c$. To describe interactions of the ions with pore wall, we use
the following Gaussian well short-range potentials

$u_{\pm}(z)=-\varepsilon_{\pm}\exp\left[-\frac{z^{2}}{2\sigma^{2}_{\pm}}\right],$
(30)

where $\varepsilon_{\pm}$ are the positive energetic parameters (depths) and
$\sigma_{\pm}$ are the effective widths of the Gaussian wells.

The SCF equations and the boundary conditions in the dimensionless units take
the form

$\begin{cases}\tilde{c}^{{}^{\prime\prime}}_{\alpha}(x)=\sum_{\gamma=\pm}\hat{\chi}_{\alpha\gamma}\tilde{v}_{\gamma}(x),\quad&\alpha=\pm\\\
u^{{}^{\prime\prime}}(x)=\frac{1}{2}(\tilde{c}_{-}(x)-\tilde{c}_{+}(x)),&x\in[0,\tilde{H}]\\\
\tilde{c}_{\alpha}(0)=\tilde{c}_{\alpha}(\tilde{H})=0,\\\
u(0)=u(\tilde{H})=u_{0},\end{cases}$ (31)

where we have introduced the following dimensionless variables:
$\tilde{c}_{\pm}={c_{\pm}}/{c}$, $c^{(b)}_{\pm}=c$, $u={q\psi}/{k_{B}T}$,
$u_{0}={q\psi_{0}}/{k_{B}T}$, $x={z}/{r_{D}}$, $\tilde{H}=H/{r_{D}}$,
$\chi_{\alpha\gamma}=(k_{B}Tr^{2}_{D}/c)^{-1}\varkappa_{\alpha\gamma}$,
$\hat{\chi}_{\alpha\gamma}=\chi^{-1}_{\alpha\gamma}$,
$\tilde{v}_{\alpha}={v_{\alpha}}/{k_{B}T}$;
$r_{D}=\sqrt{\epsilon\epsilon_{0}k_{B}T/2q^{2}c}$ is the Debye length. The
dimensionless disjoining pressure, $\tilde{\Pi}=\Pi/(ck_{B}T)$, can be written
as follows

$\displaystyle\tilde{\Pi}=-\dfrac{1}{\phi_{0}}\cdot\ln\dfrac{1-(N_{+}\tilde{c}_{+\,m}+N_{-}\tilde{c}_{-\,m})\phi_{0}}{1-(N_{+}+N_{-})\phi_{0}}-N_{+}(\tilde{c}_{+\,m}-1)-N_{-}(\tilde{c}_{-\,m}-1)$
(32)
$\displaystyle+\tilde{c}_{+\,m}+\tilde{c}_{-\,m}-2-\tilde{c}_{+\,m}\tilde{v}_{+\,m}-\tilde{c}_{-\,m}\tilde{v}_{-\,m}-2\sum\limits_{\alpha=\pm}\int\limits_{\tilde{H}/2}^{\tilde{H}}dx\tilde{c}_{\alpha}(x)\tilde{u}_{\alpha}^{\prime}(x),$

where the parameter $\phi_{0}=cv$ determines the packing of the ions in the
bulk; the subscript "$m$" denotes that corresponding variables to be
calculated at $x=\tilde{H}/2$; $\tilde{u}_{\alpha}=u_{\alpha}/k_{B}T$. For
numerical calculations, we use the following physical parameters: $T=300$ K,
$q=1.66\times 10^{-19}$ C, $\epsilon=40$, $c=1$ mol/l, $v^{1/3}=0.5$ nm. The
Debye radius in this case is $r_{D}\approx 0.2$ nm and packing parameter
$\phi_{0}\approx 0.08$. Thus, we consider the case of rather concentrated 1:1
ES with organic polar solvent like acetonitrile.

After numerical solution of the SCF equations (31), we return to the physical
units, in accordance with the above definitions. Numerical calculations were
performed using the $solve\\_bvp$ method from the Scipy library in Python[48].

Figure 1: Dimensionless ionic concentrations, $c_{\pm}/c$, corresponding to
different pore thicknesses plotted for different negative surface potentials,
$\psi_{0}$. Data are shown for $\chi_{++}=1$, $\chi_{--}=0.5$,
$\chi_{+-}=\chi_{-+}=0.1$, $\varepsilon_{\pm}=0$.

###  IV.2 Ionic concentrations

Fig. 1 illustrates how change in widths and surface potential affect the
concentration profiles of cations and anions. At narrow pore widths, both ions
concentrate predominantly at the middle of the pore, resulting in the absence
of an electric double layer (EDL) at the ES/wall boundary. However, as the
surface potential increases, the cation concentration increases and the anion
concentration decreases due to electrostatic attraction and repulsion,
respectively. When the pore width increases, electric double layers begin to
form near the walls. At a width of approximately 5 nm, the EDLs become
practically isolated, creating an electrically neutral bulk solution in the
middle of the pore. It is important to note that positive surface potentials
produce similar concentration behaviors.

Figure 2: Disjoining pressure dependences on the pore thickness plotted for
different surface potentials, $\psi_{0}$. Data are shown for $\chi_{++}=1$,
$\chi_{--}=0.5$, $\chi_{+-}=\chi_{-+}=0.1$, $\varepsilon_{\pm}=0$. Figure 3:
Ion concentration profiles at pore thicknesses up to ($H=1$ nm), at ($H=1.75$
nm), and after ($H=4$ nm) minimum disjoining pressure. Data are shown for
$\chi_{++}=1$, $\chi_{--}=0.5$, $\chi_{+-}=\chi_{-+}=0.1$,
$\varepsilon_{\pm}=0$, $\psi_{0}=-0.1$ V. Figure 4: Disjoining pressure
dependence on the pore thickness. Data are shown for $\chi_{++}=1$,
$\chi_{--}=0.5$, $\chi_{+-}=\chi_{-+}=0.1$, $\varepsilon_{\pm}=0$,
$\psi_{0}=-0.1$ V. Figure 5: Disjoining pressure dependences on the pore
thickness plotted for different specific adsorption parameter,
$\varepsilon_{+}$. Data are show for $\chi_{++}=1$, $\chi_{--}=0.5$,
$\chi_{+-}=\chi_{-+}=0.1$, $\varepsilon_{-}=0$, $\sigma_{+}/v^{1/3}=1$,
$\psi_{0}=-0.1$ V.

###  IV.3 Disjoining pressure

In this subsection, we will examine the behavior of disjoining pressure. Fig.
2 displays the disjoining pressure as a function of pore width at various
negative surface potentials. At sufficiently small pore thickness, the
disjoining pressure exhibits non-monotonic behavior, indicating effective
attraction between the walls ($\Pi<0$). This minimum is more pronounced with
an increase in surface potential absolute value. To comprehend the physical
nature of this behavior, we will analyse the ionic concentration profiles at
pore widths before ($H=1$ nm), at ($H=1.75$ nm), and after ($H=4$ nm) the
minimum (see Fig. 3). At a narrow pore width of 1 nm, unimodal concentration
profiles are present, with ions concentrated at the pore centre (see Fig. 4).
At $H=1.75$ nm, we observe maxima in cation concentration profiles near the
walls and minimum at the center of the pore. This suggests that a greater
number of cations migrate from the pore center to the walls, in an attempt to
form EDLs. In this configuration, there are three layers: two layers of
cations near the pore walls and a layer of mixed cations and anions in
between. Formation of such "sandwich" structure manifests itself in the
effective "structural" attraction between the walls. At pore width of $4$ nm,
as previously mentioned, we can observe isolated EDLs near the walls, with
bulk solution present in the middle of the pore. In this region, the
disjoining pressure exponentially damps (as shown in the upper inset in Fig.
4) in accordance with classical DLVO theory [49, 50].

We would like to state that the current model is inadequate for capturing the
layering effect that occur in both room temperature ionic liquids and highly
concentrated electrolyte solutions on the electrified interfaces [25, 22, 19]
which results in the strong fluctuations on the disjoining pressure. To
accurately describe these phenomena, it is necessary to take into account
higher order derivative terms in the GTP [51] or apply nonlocal theories [22,
33]. Nevertheless, present simple statistical theory could be used for
modeling of the ES of moderate concentrations, where the ion layering is
irrelevant.

Figure 6: Ion concentration profiles corresponding to different pore
thicknesses plotted for different specific adsorption parameter,
$\varepsilon_{+}$. Data are shown for $\chi_{++}=1$, $\chi_{--}=0.5$,
$\chi_{+-}=\chi_{-+}=0.1$, $\varepsilon_{-}=0$, $\sigma_{+}/v^{1/3}=1$,
$\psi_{0}=-0.1$ V. Figure 7: Disjoining pressure dependences on the pore
thickness plotted for different ionic size asymmetry $N_{+}=[1,2,3]$ and
$N_{-}=1$. Data are shown for $\chi_{++}=1$, $\chi_{--}=0.5$,
$\chi_{+-}=\chi_{-+}=0.1$, $\varepsilon_{\pm}=0$, $\sigma_{+}/v^{1/3}=1$,
$\psi_{0}=-0.1$ V. Figure 8: Ion concentration profiles corresponding to
different pore thicknesses plotted for different ionic size asymmetry
$N_{+}=[1,2,3]$ and $N_{-}=1$. Data are shown for $\chi_{++}=1$,
$\chi_{--}=0.5$, $\chi_{+-}=\chi_{-+}=0.1$, $\varepsilon_{\pm}=0$,
$\sigma_{+}/v^{1/3}=1$, $\psi_{0}=-0.1$ V. Figure 9: Disjoining pressure
dependences on the pore thickness plotted for different ionic size asymmetry
$N_{+}=[1,2,3]$ and $N_{-}=1$. Data are shown for $\chi_{++}=1$,
$\chi_{--}=0.5$, $\chi_{+-}=\chi_{-+}=0.1$, $\varepsilon_{\pm}=0$,
$\sigma_{+}/v^{1/3}=1$, $\psi_{0}=0.1$ V. Figure 10: Ion concentration
profiles corresponding to different pore thicknesses plotted for different
ionic size asymmetry $N_{+}=[1,2,3]$ and $N_{-}=1$. Data are shown for
$\chi_{++}=1$, $\chi_{--}=0.5$, $\chi_{+-}=\chi_{-+}=0.1$,
$\varepsilon_{\pm}=0$, $\sigma_{+}/v^{1/3}=1$, $\psi_{0}=0.1$ V.

###  IV.4 Effects of specific adsorption

It is informative to discuss how the specific adsorption of ions, described by
the Gaussian-well potentials (30), affects the local ionic concentrations and
disjoining pressure. Fig. 5 illustrates various ionic concentration profiles
for different pore thicknesses while altering the specific cation adsorption
energy parameter, $\varepsilon_{+}$. The results demonstrate that for narrow
pores ($H=1$ nm), an increase in $\varepsilon_{+}$ results in a predictable
increase in the cation local concentration. However, it also leads to a
substantial increase in anion local concentration due to cation-anion Coulomb
attraction. With larger $H$, an increase in $\varepsilon_{+}$ results in
sharper cation concentration profile maxima, leading to the formation of a
dense part of the EDLs (Stern layers). Furthermore, $\varepsilon_{+}$
enhancement intensifies the minimum on the disjoining pressure profiles and
shifts it to narrower pores. This behavior can be attributed to the fact that
cation specific adsorption promotes EDL formation at smaller $H$ than is
possible in the absence of specific adsorption. This is reflected in
pronounced cation concentration profile maxima in Fig. 6 at $H=1$ nm and
$\varepsilon_{+}=10k_{B}T$.

###  IV.5 Effects of ionic size asymmetry

Previously, we only discussed the case of equal ionic sizes, where
$N_{+}=N_{-}=1$. Now, we will briefly discuss the effects of ionic size
asymmetry on disjoining pressure behavior. Fig. 7 shows typical disjoining
pressure profiles for different cation sizes, specifically $N_{+}=1$,
$N_{+}=2$, and $N_{+}=3$, with a fixed anion size of $N_{-}=1$ for negatively
charged pore walls ($\psi_{0}=-0.1$ V). As is seen in the figure, an increase
in cation size relative to anion size practically eliminates the minimum on
the disjoining pressure curve. This is because a larger cation size increases
steric interactions and decreases the local cation concentration within the
pore volume (see Fig. 8), that in turn "compensates" for the structural
attraction of the walls. As in present case, according to observations made in
ref. [52], the electrostatic correlation of like-charge attraction between
surfaces can be neutralized by increasing the diameter of counterions.

When considering walls with a positive charge ($\psi_{0}=0.1$ V), the size of
cations has an opposite effect on the disjoining pressure - resulting in a
more pronounced minimum (see Fig. 9). This may seem counterintuitive, but it
can be attributed to the well-known phenomenon of Asakura-Oosawa depletion
forces [53]. As cations are placed in close proximity to positively charged
walls, a wall-cation depletion occurs [54, 55]. The latter effect results in
the expulsion of both cations and anions from the immediate vicinity of the
wall surface, leading to their accumulation behind the cations. Therefore, an
increase in cation effective size makes the depletion cation-wall attraction
stronger. This is evident from the concentration profiles depicted in Fig. 10
\- the larger the cation size, the higher the local cation and anion
concentration. We would like to mention that conventional modified Poisson-
Boltzmann equations [31, 44, 40, 41, 36] do not account for depletion effects.
However, incorporating structural interactions in the GTP enables us to
qualitatively describe the depletion effect. It should be noted that the
expulsion of counterions from the vicinity of a charged surface can be
predicted using mean-field modified Poisson-Boltzmann models, provided that
rotational ordering of finite water dipoles is taken into account as well [34,
44]. Specifically, an article [56] has demonstrated that for sufficiently high
surface charge densities, the traditional trend of increasing counterion
number density towards the charged surface may be reversed, resulting in a
decrease in the number of counterions near the surface. This phenomenon is
caused by the reduction of partial depletion of water molecules at high
surface charge densities, due to the competition between the counterion
Boltzmann factor and the rotationally averaged water Boltzmann factor in the
closely-bound region [56]. The simulation results presented in ref.[35] are
consistent with this conclusion. Therefore, the theoretical model must
consider both the orientational ordering of finite-sized water molecules and
the ionic size asymmetry, as demonstrated in [57].

##  V Concluding remarks

To summarize, we employed Cahn-Hilliard-like model of inhomogeneous
electrolyte solutions that accounted for the structural and steric
interactions between ions and incorporated external potentials to examine
their specific interactions with pore walls. Our main focus was to observe
changes in ionic concentration profiles and disjoining pressure as the pore
width varied. We derived a general expression for disjoining pressure starting
from the local mechanical equilibrium conditions. Our results showed that
considering the structural interactions of ions resulted in a pronounced
minimum on disjoining pressure profiles at small pore widths. We attributed
this minimum to the formation of electric double layers on the electrified
pore surfaces. Moreover, including the attractive interactions between ions
and pore walls enhanced and shifted the minimum to smaller pore thicknesses.
We believe that our theoretical findings could be of interest to
electrochemical engineers working on supercapacitors that utilize porous
electrodes impregnated with moderately concentrated electrolyte solutions.

Data availability statement. The data that supports the findings of this study
are available within the article.

Acknowledgements. This research is supported by the Russian Science Foundation
(No. 21-11-00031). The numerical calculations were partially performed on the
supercomputer facilities provided by NRU HSE.

## References

  * Hassan  _et al._ [2023] H. Hassan, M. W. Iqbal, S. Gouadria, A. M. Afzal, and H. Hegazy, Journal of Energy Storage 66, 107448 (2023).
  * Xia  _et al._ [2017] L. Xia, L. Yu, D. Hu, and G. Z. Chen, Materials Chemistry Frontiers 1, 584 (2017).
  * Zheng  _et al._ [2017] J. Zheng, J. A. Lochala, A. Kwok, Z. D. Deng, and J. Xiao, Advanced Science 4, 1700032 (2017).
  * Xu  _et al._ [2023] Y. Xu, Z. Mao, R. Qu, J. Wang, J. Yu, X. Luo, M. Shi, X. Mao, J. Ding, and B. Liu, Nature Water , 1 (2023).
  * Sharma, Borkute, and Gumfekar [2022] V. Sharma, G. Borkute, and S. P. Gumfekar, Chemical Engineering Journal 433, 133823 (2022).
  * Aluru  _et al._ [2023] N. R. Aluru, F. Aydin, M. Z. Bazant, D. Blankschtein, A. H. Brozena, J. P. de Souza, M. Elimelech, S. Faucher, J. T. Fourkas, V. B. Koman, _et al._ , Chemical reviews 123, 2737 (2023).
  * Shin, Guiver, and Lee [2017] D. W. Shin, M. D. Guiver, and Y. M. Lee, Chemical reviews 117, 4759 (2017).
  * Luo, Abdu, and Wessling [2018] T. Luo, S. Abdu, and M. Wessling, Journal of membrane science 555, 429 (2018).
  * Ben-Yaakov, Andelman, and Podgornik [2011] D. Ben-Yaakov, D. Andelman, and R. Podgornik, The Journal of chemical physics 134, 074705 (2011).
  * Frydel [2011] D. Frydel, The Journal of chemical physics 134, 234704 (2011).
  * Hatlo, Van Roij, and Lue [2012] M. M. Hatlo, R. Van Roij, and L. Lue, Europhysics Letters 97, 28010 (2012).
  * Goodwin, Feng, and Kornyshev [2017] Z. A. Goodwin, G. Feng, and A. A. Kornyshev, Electrochimica Acta 225, 190 (2017).
  * Uematsu, Netz, and Bonthuis [2018] Y. Uematsu, R. R. Netz, and D. J. Bonthuis, Journal of Physics: Condensed Matter 30, 064002 (2018).
  * Budkov, Zavarzin, and Kolesnikov [2021] Y. A. Budkov, S. V. Zavarzin, and A. L. Kolesnikov, The Journal of Physical Chemistry C 125, 21151 (2021).
  * Budkov  _et al._ [2020] Y. A. Budkov, A. V. Sergeev, S. V. Zavarzin, and A. L. Kolesnikov, The Journal of Physical Chemistry C 124, 16308 (2020).
  * Podgornik [2018] R. Podgornik, The Journal of chemical physics 149, 104701 (2018).
  * Naji  _et al._ [2013] A. Naji, M. Kanduč, J. Forsman, and R. Podgornik, The Journal of chemical physics 139, 150901 (2013).
  * Blossey [2023] R. Blossey, in  _The Poisson-Boltzmann Equation: An Introduction_ (Springer, 2023) pp. 53–96.
  * Budkov and Kolesnikov [2021] Y. A. Budkov and A. L. Kolesnikov, Current Opinion in Electrochemistry 33, 100931 (2021).
  * Misra  _et al._ [2019] R. P. Misra, J. P. de Souza, D. Blankschtein, and M. Z. Bazant, Langmuir 35, 11550 (2019).
  * de Souza and Bazant [2020] J. P. de Souza and M. Z. Bazant, The Journal of Physical Chemistry C 124, 11414 (2020).
  * de Souza  _et al._ [2020] J. P. de Souza, Z. A. Goodwin, M. McEldrew, A. A. Kornyshev, and M. Z. Bazant, Physical Review Letters 125, 116001 (2020).
  * Shi  _et al._ [2023] K. Shi, E. R. Smith, E. E. Santiso, and K. E. Gubbins, The Journal of Chemical Physics 158, 040901 (2023).
  * Kolesnikov, Budkov, and Gor [2021] A. Kolesnikov, Y. A. Budkov, and G. Gor, Journal of Physics: Condensed Matter 34, 063002 (2021).
  * Gurina  _et al._ [2022] D. Gurina, E. Odintsova, A. Kolesnikov, M. Kiselev, and Y. Budkov, Journal of Molecular Liquids 366, 120307 (2022).
  * Kolesnikov, Mazur, and Budkov [2022] A. L. Kolesnikov, D. A. Mazur, and Y. A. Budkov, Europhysics Letters 140, 16001 (2022).
  * Koczwara  _et al._ [2017] C. Koczwara, S. Rumswinkel, C. Prehal, N. Jackel, M. S. Elsasser, H. Amenitsch, V. Presser, N. Husing, and O. Paris, ACS applied materials & interfaces 9, 23319 (2017).
  * Chen  _et al._ [2022] Z. Chen, D. L. Danilov, R.-A. Eichel, and P. H. Notten, Advanced Energy Materials 12, 2201506 (2022).
  * Da Silva  _et al._ [2020] L. M. Da Silva, R. Cesar, C. M. Moreira, J. H. Santos, L. G. De Souza, B. M. Pires, R. Vicentini, W. Nunes, and H. Zanin, Energy storage materials 27, 555 (2020).
  * Li  _et al._ [2018] X. Li, J. Shao, S.-K. Kim, C. Yao, J. Wang, Y.-R. Miao, Q. Zheng, P. Sun, R. Zhang, and P. V. Braun, Nature communications 9, 2578 (2018).
  * Budkov and Kolesnikov [2022] Y. A. Budkov and A. L. Kolesnikov, Journal of Statistical Mechanics: Theory and Experiment 2022, 053205 (2022).
  * Budkov and Kalikin [2023] Y. A. Budkov and N. N. Kalikin, Physical Review E 107, 024503 (2023).
  * Brandyshev and Budkov [2023] P. E. Brandyshev and Y. A. Budkov, The Journal of chemical physics 158, 174114 (2023).
  * Gongadze  _et al._ [2014] E. Gongadze, A. Velikonja, Š. Perutkova, P. Kramar, A. Maček-Lebar, V. Kralj-Iglič, and A. Iglič, Electrochimica Acta 126, 42 (2014).
  * Marcovitz, Naftaly, and Levy [2015] A. Marcovitz, A. Naftaly, and Y. Levy, The Journal of Chemical Physics 142, 02B614_1 (2015).
  * Maggs and Podgornik [2016] A. Maggs and R. Podgornik, Soft matter 12, 1219 (2016).
  * Bazant, Storey, and Kornyshev [2011] M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, Physical review letters 106, 046102 (2011).
  * Blossey, Maggs, and Podgornik [2017] R. Blossey, A. Maggs, and R. Podgornik, Physical Review E 95, 060602 (2017).
  * Cahn [1965] J. W. Cahn, The Journal of chemical physics 42, 93 (1965).
  * Borukhov, Andelman, and Orland [1997] I. Borukhov, D. Andelman, and H. Orland, Physical review letters 79, 435 (1997).
  * Kornyshev [2007] A. A. Kornyshev, The Journal of Physical Chemistry B 111, 5545 (2007).
  * Kralj-Iglič and Iglič [1996] V. Kralj-Iglič and A. Iglič, Journal de Physique II 6, 477 (1996).
  * Budkov  _et al._ [2018] Y. A. Budkov, A. L. Kolesnikov, Z. A. Goodwin, M. G. Kiselev, and A. A. Kornyshev, Electrochimica Acta 284, 346 (2018).
  * Iglič, Gongadze, and Bohinc [2010] A. Iglič, E. Gongadze, and K. Bohinc, Bioelectrochemistry 79, 223 (2010).
  * Budkov, Kolesnikov, and Kiselev [2016] Y. A. Budkov, A. Kolesnikov, and M. Kiselev, The Journal of chemical physics 144, 184703 (2016).
  * Budkov, Kolesnikov, and Kiselev [2015] Y. A. Budkov, A. Kolesnikov, and M. Kiselev, Europhysics Letters 111, 28002 (2015).
  * Nakayama and Andelman [2015] Y. Nakayama and D. Andelman, The Journal of chemical physics 142, 044706 (2015).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, _et al._ , Nature methods 17, 261 (2020).
  * Derjaguin, Churaev, and Muller [1987] B. Derjaguin, N. Churaev, and V. Muller, in  _Surface Forces_ (Springer, 1987) pp. 293–310.
  * Israelachvili [2011] J. N. Israelachvili,  _Intermolecular and surface forces_ (Academic press, 2011).
  * Ciach [2018] A. Ciach, Journal of Molecular Liquids 270, 138 (2018).
  * Zelko  _et al._ [2010] J. Zelko, A. Iglič, V. Kralj-Iglič, and P. Kumar, The Journal of chemical physics 133 (2010).
  * Barrat and Hansen [2003] J.-L. Barrat and J.-P. Hansen,  _Basic concepts for simple and complex liquids_ (Cambridge University Press, 2003).
  * Taboada-Serrano, Yiacoumi, and Tsouris [2006] P. Taboada-Serrano, S. Yiacoumi, and C. Tsouris, The Journal of chemical physics 125, 054716 (2006).
  * Kewalramani  _et al._ [2016] S. Kewalramani, G. I. Guerrero-García, L. M. Moreau, J. W. Zwanikken, C. A. Mirkin, M. Olvera de la Cruz, and M. J. Bedzyk, ACS central science 2, 219 (2016).
  * Gongadze and Iglič [2013] E. Gongadze and A. Iglič, Gen Phys Biophys 21, 143 (2013).
  * Gongadze and Iglič [2015] E. Gongadze and A. Iglič, Electrochimica Acta 178, 541 (2015).