holds, we can show that $\mathcalbf{C}_{-}$ is stable in certain limits.

First, we notice that $\mathcalbf{C}_{-}$ is stable for sufficiently large
$\beta$. This follows from the expansion

$\displaystyle
c_{1}^{-}c_{2}^{-}-c_{0}^{-}=\dfrac{\zeta^{3}}{\eta^{2}}\beta^{2}+O(\beta),$
(87a) obtained using Mathematica to assist with manipulating complicated
algebraic expressions.

Next, $\upDelta=\smash{\left|\beta\zeta\kappa-V\delta\eta/\zeta\right|}+O(U)$
from its definition. If ${V<0}$, then $\beta\zeta^{2}\kappa>V\delta\eta$ is
clearly true. This inequality also holds if ${V>0}$ and $\mathcalbf{C}_{-}$ is
feasible, from Eqs. (75). Hence
${\upDelta=(\beta\zeta\kappa-V\delta\eta/\zeta)+O(U)}$, and we find

$\displaystyle
c_{1}^{-}c_{2}^{-}-c_{0}^{-}=\dfrac{\zeta}{\eta^{2}}u(u+\zeta\eta)+O(U)>0,$
(87b)

where $u=\beta\zeta+(\gamma+\delta)\eta$. Hence $\mathcalbf{C}_{-}$ is stable
for small enough $|U|$ if it is feasible, i.e. if $U,V<0$ or conditions (75)
hold.

Moreover, as $\beta\rightarrow 0$ at $\kappa\not=0$,
$\mathcalbf{C}_{-}\rightarrow\mathcalbf{C}$, so, by continuity,
$\mathcalbf{C}_{-}$ is stable for small enough $\beta$, since $\mathcalbf{C}$
is, at least under the conditions discussed above and in particular if the
persister scalings (7) are satisfied.

Finally, we discuss the stability of $\mathcalbf{C}_{-}$ on, and hence by
continuity near, the feasibility boundary $\beta=\beta_{\ast}$ defined in Eqs.
(75). On this boundary, $\upDelta=0$ and hence $c_{0}^{-}=0$. Since
$c_{2}^{-}>0$, stability is equivalent with $c_{1}^{-}>0$ there, as in the
analysis of Eqs. (26). Direct computation yields

$\displaystyle(2\beta\kappa)^{2}\eta\vartheta c_{1}^{-}$
$\displaystyle=\beta\kappa(\beta\zeta\kappa-\delta
W)\left[\beta\zeta\kappa+\delta W+2\eta(\delta+\gamma\kappa)\right]$
$\displaystyle\quad+(\beta\zeta\kappa+\delta
W)\left[(\alpha\vartheta-\eta)(2\beta\delta\vartheta-\beta\zeta\kappa+\delta
W)\right.$ $\displaystyle\hskip
76.82243pt+\left.2\beta\delta\kappa\eta\vartheta\right].$ (88a) From Eqs.
(75), feasibility requires $V>0$, and so, from definition (73a),
$\vartheta>\zeta$. Moreover, we expect, from the persister scalings (7), that
$\delta>\kappa$. Thus
$2\beta\delta\vartheta-\beta\zeta\kappa>\beta\zeta\kappa$. Since
${\beta\zeta\kappa>\delta|W|}$, Eq. (88a) then shows that $c_{1}^{-}>0$
provided that $\alpha\vartheta>\eta$. Conversely, if $\alpha\vartheta<\eta$ or
$\kappa<\delta$, then $c_{1}^{-}<0$ is possible (not shown), but we are not
aware of simple conditions that ensure the presence or absence of instability
in that case. However, after inserting the explicit expressions for
$\beta=\beta_{\ast}$ and $W$ in terms of the other model parameters into Eq.
(88a) and expanding using Mathematica, $\displaystyle
c_{1}^{-}=\dfrac{64\delta^{3}\eta^{3}\vartheta^{3}}{\zeta^{4}}+O\bigl{(}\zeta^{-3}\bigr{)},$
(88b) and hence $\mathcalbf{C}_{-}$ is stable on the feasibility boundary
$\beta=\beta_{\ast}$ for sufficiently small $\zeta$. We stress that this does
not follow from Eq. (87a), since the “sufficiently large” there depends on the
other model parameters and, in particular, on $\zeta$. Again expressing
$\beta=\beta_{\ast}$ and $W$ in terms of the other model parameters, we also
find $\displaystyle c_{1}^{-}$
$\displaystyle=\dfrac{4\delta\eta}{\vartheta}(\delta+\gamma\kappa)^{2}W^{2}$
$\displaystyle\text{if
}\zeta=\vartheta\left(1+\dfrac{\gamma\kappa}{\delta}\right)^{-1},$ (88c)
$\displaystyle=\dfrac{4\alpha\delta^{3}\vartheta}{\eta\kappa}(\alpha\vartheta-\eta+\eta\kappa)W^{2}$
$\displaystyle\text{if }\zeta=\dfrac{\eta}{\alpha}.$ (88d)

If $\eta/\alpha<\vartheta/(1+\gamma\kappa/\delta)<\vartheta$, then
$\alpha\vartheta-\eta>0$. This therefore shows that $c_{1}^{-}>0$ at
$\zeta=\min{\\{\eta/\alpha,\vartheta/(1+\gamma\kappa/\delta)\\}}$, and hence,
on referring to Eq. (76), that $\mathcalbf{C}_{-}$ is stable at the “endpoint”
of the feasibility boundary under discussion. For $\zeta=\eta/\alpha$, this
result confirms a particular case of the conclusion that we drew from Eq.
(87b).

####  B.3.3 Stability of the trivial steady states and permanence

Finally, we analyze the stability of the trivial steady states to obtain
conditions for permanent coexistence. As usual, $\mathcalbf{O}$ is always
unstable, with the Jacobian of Eqs. (28) at this fixed point being

$\displaystyle\left(\begin{array}[]{ccc}1-\gamma&\delta&0\\\
\gamma&-\delta&0\\\ 0&0&\zeta\end{array}\right),$ (89d) with an unstable
eigenvalue $\zeta>0$. Again omitting entries that do not affect the stability,
the Jacobian at $\mathcalbf{A}$ is
$\displaystyle\left(\begin{array}[]{ccc}-U/\eta-(\gamma+\beta\zeta/\eta)&\delta&0\\\
\gamma+\beta\zeta/\eta&-\delta&0\\\ &&-\zeta\end{array}\right).$ (89h) The
Jacobians at $\mathcalbf{O}$ and $\mathcalbf{A}$ are thus actually equal to
the Jacobians of the corresponding steady states for model (26) given by Eqs.
(51d) and (51h). Moreover, the Jacobian evaluated at $\mathcalbf{B}$ is
$\displaystyle\left(\begin{array}[]{ccc}-\gamma-J^{-1}&\delta-\kappa
J^{-1}&\\\ \gamma&-\delta&\\\ 0&0&-VJ^{-1}\\!\\!\end{array}\right),$ (89l)

with $J=1+\gamma\kappa/\delta$ again. On referring to the results for Eqs.
(26) in the first part of this Appendix, we find that $\mathcalbf{A}$ is
stable if and only if $U>0$. Moreover, on computing the eigenvalues of the
matrix in Eq. (89l), we obtain again that $\mathcalbf{B}$ is stable if and
only if $V>0$.

What is more, the eigenvectors of Jacobians in Eqs. (89) are easily seen to
have the same geometric properties as those of the corresponding Jacobians of
model (26) analyzed in the first part of this Appendix. We can therefore
conclude, as we have done there, that the two species coexist permanently,
irrespective of the initial conditions, if $\mathcalbf{A}$ and $\mathcalbf{B}$
are unstable, i.e. if $U,V<0$.

####  B.3.4 Stability diagrams of Eqs. (28)

We assemble all of these analytical results into the stability diagrams shown
in Fig. 13 for $\kappa>0$. Again, we will not discuss the singular case
$\kappa=0$. These exact diagrams establish some of the features that we have
already seen in the numerical stability diagrams in Figs. 6(e),(f).

Figure 13: Feasibility and stability of the coexistence states
$\mathcalbf{C}_{\pm}$ of the simplified model (28) in the $(\zeta,\beta)$
diagram for $\kappa>0$, in the cases (a)
$\eta/\alpha>\vartheta/(1+\kappa\gamma/\delta)$ and (b)
$\eta/\alpha<\vartheta/(1+\kappa\gamma/\delta)$, assuming that the persister
scalings (7) hold. Coexistence is feasible in the region bounded by the thick
solid black lines. Both $\mathcalbf{C}_{+}$ and $\mathcalbf{C}_{-}$ are
unstable in the region marked “unstable”, but $\mathcalbf{C}_{-}$ is stable
(and $\mathcalbf{C}_{+}$ is unstable) in the region marked “stable”. Whether
$\mathcalbf{C}_{-}$ is unstable in a subregion of the “stable” region must be
determined numerically. In the diagonally hatched region of parameter space in
panel (a), the trivial steady states $\mathcalbf{A}$ and $\mathcalbf{B}$ are
unstable, and the two species coexist permanently. If the “unstable” region
intersects this region, all steady states of Eqs. (28) are unstable within
this intersection. The horizontally hatched region is the region of stable
steady-state coexistence for averaged models of Eqs. (28) without phenotypic
variation and with stochastic switching only; coexistence may be stable in the
average with respect to $\mathcalbf{C}_{-}$ and in that with respect to
$\mathcalbf{C}_{+}$ (Appendix C).

Once again, the combination of responsive switching and bacteria-persister
interactions ($\beta\kappa>0$) leads to expanded regions of feasibility and
stability, compared to $\beta\kappa=0$, if $\beta>\beta_{\ast}^{\min}$, with
$\beta_{\ast}^{\min}$ now given by Eq. (77). As for models (26) and (27),
$\beta_{\ast}^{\min}=O(\varepsilon)$ is possible even if
$\kappa=O(\varepsilon)$ provided that $\eta-\alpha\vartheta\lesssim
O\bigl{(}\varepsilon^{2}\bigr{)}$. This is the now familiar requirement that
the intra-species competitions of bacteria and competitors be sufficiently
close to their inter-species competitions (Table 1).

We conclude our discussion by comparing Eqs. (28) to their averages without
phenotypic variation and with stochastic switching only. As shown in Appendix
C, the conditions for stability of coexistence if $\beta=0$, i.e. for
stochastic switching only, is precisely that of the averaged model without
phenotypic variation. This equivalence assumes that
${\delta[\delta+\vartheta(\gamma+\delta)]+\gamma\kappa(\delta-\zeta)>0}$,
expected to hold given the persister scalings (7). Thus stochastic switching
on its own has again no effect on stability. The comparison for $\beta>0$,
shown in Fig. 13, stresses again the importance of careful model averaging:
Although the region of parameter space in which coexistence is stable widens
with increasing $\beta$, the corresponding variation of the model parameter
means that coexistence is stable in the averaged model if and only if it is
feasible. Thus, if regions of instability as sketched in Fig. 13 exist,
steady-state coexistence is destabilized there by the combination of
responsive switching and competition between bacteria and persisters
($\beta\kappa>0$). Again, this argument does not preclude non-steady-state
coexistence there, unless there is a region, as shown in Fig. 13(a) in which
all steady states are unstable and hence coexistence is permanent there.

Anyway, numerical experiments (not shown) indicate that, while such regions of
instability do exist, they are rare for parameter values consistent with the
persister scalings (7). We have no analytical understanding of this
observation, however. Given that signatures of the possible regions of
instability for model (28) sketched in Fig. 13 are absent from the numerical
stablity diagrams in Figs. 6(e),(f), too, this suggests that the destabilizing
effect of the combination of responsive switching and bacteria-persister
competition ($\beta\kappa>0$) is weaker than the stabilizing effect of
responsive switching with persister-competitor interactions ($\beta\iota>0$)
and persister growth ($\beta\mu>0$) revealed respectively by the analyses of
models (26) and (27) and also seen in the numerical results in Figs. 6(a)–(c).

##  APPENDIX C COEXISTENCE WITHOUT PHENOTYPIC VARIATION

In this Appendix, we briefly rederive the classical results for the stability
of two-species Lotka–Volterra competition models [27]. We then obtain the
stability conditions for averaged models of this Lotka–Volterra form and that
correspond to the simplified models (26), (27), (28). We combine these results
with the stability calculations in Appendix B to compare the stability of
coexistence in these simplified models and in their averages without
phenotypic variation and with stochastic switching only.

###  C.1 Coexistence in a two-species Lotka–Volterra model without phenotypic
variation

Coexistence in a two-species Lotka–Volterra competition model without
phenotypic variation is a classical problem, discussed, for example, in Ref.
[27]. With the aim in mind of comparing this model to the models with
phenotypic variation considered in the main text, it will be useful to
rederive the results briefly. We consider the competition between species
$A^{\prime},B^{\prime}$ described by the differential equations

$\displaystyle\dot{B}^{\prime}=B^{\prime}\bigl{(}\omega^{\prime}-\alpha^{\prime}A^{\prime}-\chi^{\prime}B^{\prime}\bigr{)},$
$\displaystyle\dot{A}^{\prime}=A^{\prime}\bigl{(}\zeta^{\prime}-\eta^{\prime}A^{\prime}-\vartheta^{\prime}B^{\prime}\bigr{)},$
(90)

wherein
$\alpha^{\prime},\zeta^{\prime},\eta^{\prime},\vartheta^{\prime},\chi^{\prime},\omega^{\prime}\geqslant
0$ are parameters. As in our derivation of Eqs. (6) and (26), we could have
scaled time and $A^{\prime},B^{\prime}$ to set some parameters equal to $1$,
but we have not done so in order to be able to relate this model to a model
with phenotypic variation in the next subsection [37]. Equations (90) have a
single coexistence state $\mathcalbf{C}^{\prime}=(b^{\prime},a^{\prime})$,
where $b^{\prime}=u^{\prime}/w^{\prime}$ and
$a^{\prime}=v^{\prime}/w^{\prime}$, with

$\displaystyle u^{\prime}$
$\displaystyle=\alpha^{\prime}\zeta^{\prime}-\eta^{\prime}\omega^{\prime},\quad
v^{\prime}=\vartheta^{\prime}\omega^{\prime}-\zeta^{\prime}\chi^{\prime},$
(91a) $\displaystyle w^{\prime}$
$\displaystyle=\alpha^{\prime}\vartheta^{\prime}-\eta^{\prime}\chi^{\prime}=\dfrac{\chi^{\prime}u^{\prime}+\alpha^{\prime}v^{\prime}}{\omega^{\prime}}.$
(91b)

Hence the coexistence state is feasible if and only if $u^{\prime},v^{\prime}$
have the same sign. The Jacobian matrix evaluated at this steady state is

$\displaystyle\left(\begin{array}[]{cc}-\chi^{\prime}b^{\prime}&-\alpha^{\prime}b^{\prime}\\\
-\vartheta^{\prime}a^{\prime}&-\eta^{\prime}a^{\prime}\end{array}\right).$
(94)

Since $\mathrm{tr}=-\chi^{\prime}b-\eta^{\prime}a^{\prime}<0$, classical
stability results [27] imply that coexistence is stable if and only if
$0<\det=-u^{\prime}v^{\prime}/w^{\prime}$. Hence coexistence state is feasible
and stable if and only if $u^{\prime},v^{\prime}<0$, or equivalently, if and
only if

$\displaystyle\dfrac{\vartheta^{\prime}\omega^{\prime}}{\chi^{\prime}}<\zeta^{\prime}<\dfrac{\eta^{\prime}\omega^{\prime}}{\alpha^{\prime}}.$
(95)

Equations (90) have two additional nonzero steady states, namely
$\mathcalbf{A}^{\prime}=(0,\zeta^{\prime}/\eta^{\prime})$ and
$\mathcalbf{B}^{\prime}=(\omega^{\prime}/\chi^{\prime},0)$, at which the
Jacobian matrix evaluates to

$\displaystyle\left(\begin{array}[]{cc}-u^{\prime}/\eta^{\prime}&0\\\
-\vartheta^{\prime}\zeta^{\prime}/\eta^{\prime}&-\zeta^{\prime}\end{array}\right)\quad\text{and}\quad\left(\begin{array}[]{cc}-\omega^{\prime}&-\alpha^{\prime}\omega^{\prime}/\chi^{\prime}\\\
0&-v^{\prime}/\chi^{\prime}\end{array}\right),$ (100)

respectively. Hence $\mathcalbf{A}^{\prime}$ and $\mathcalbf{B}^{\prime}$ are
stable if $u^{\prime}>0$ and $v^{\prime}>0$, respectively. Comparing these
parameter ranges for stability, it follows that either exactly one of
$\mathcalbf{A}^{\prime},\mathcalbf{B^{\prime}},\mathcalbf{C}^{\prime}$ is
stable, or $\mathcalbf{A}^{\prime},\mathcalbf{B}^{\prime}$ are both stable.
From arbitrary initial conditions, Eqs. (90) converge to the stable steady
state if is unique; if $\mathcalbf{A^{\prime}},\mathcalbf{B}^{\prime}$ are
both stable, a separatrix through $\mathcalbf{C}^{\prime}$ separates initial
conditions converging to $\mathcalbf{A}^{\prime}$ from those converging to
$\mathcalbf{B}^{\prime}$ [27].

###  C.2 Averages without phenotypic variation of the two-species model (6)
and of the simplified models (26), (27), (28)

We seek to describe the populations $B,P,A$ that evolve according to Eqs. (6)
by an averaged model without phenotypic variation and two populations
$B^{\prime}$, which corresponds to $B$ and $P$, and $A^{\prime}$, which
corresponds to $A$ [Fig. 1(b)]. We have introduced this averaging in Ref.
[15], and we have motivated it again here, in Section II.

####  C.2.1 Derivation of the averaged model

A coexistence equilibrium $\mathcalbf{C}=(b,p,a)$ of Eqs. (6) is consistent
with the equilibrium $\mathcalbf{C}^{\prime}=(b^{\prime},a^{\prime})$ of Eqs.
(90) if and only if the populations, the births, and the competition are equal
at equilibrium, i.e. if and only if

$\displaystyle b^{\prime}$ $\displaystyle=b+p,$ $\displaystyle a^{\prime}$
$\displaystyle=a,$ (101a) $\displaystyle\omega^{\prime}b^{\prime}$
$\displaystyle=b+\mu p,$ $\displaystyle\zeta^{\prime}a^{\prime}$
$\displaystyle=\zeta a,$ (101b)
$\displaystyle\alpha^{\prime}a^{\prime}b^{\prime}$ $\displaystyle=a(\alpha
b+\xi p),$ $\displaystyle\vartheta^{\prime}a^{\prime}b^{\prime}$
$\displaystyle=a(\vartheta b+\iota p),$ (101c)
$\displaystyle\chi^{\prime}b^{\prime 2}$
$\displaystyle=b^{2}+(\kappa+\varpi)bp+\varsigma p^{2},$
$\displaystyle\eta^{\prime}a^{\prime 2}$ $\displaystyle=\eta a^{2}.$ (101d)

These conditions are to Eqs. (6) what conditions (25) are to Eqs. (1). We let
$q=p/b$, so that, on introducing $Q=(1+q)^{-1}$, they reduce to

$\displaystyle\alpha^{\prime}$ $\displaystyle=Q(\alpha+\xi q),$
$\displaystyle\zeta^{\prime}$ $\displaystyle=\zeta,$
$\displaystyle\omega^{\prime}$ $\displaystyle=Q(1+\mu q),$ (102a)
$\displaystyle\vartheta^{\prime}$ $\displaystyle=Q(\vartheta+\iota q),$
$\displaystyle\eta^{\prime}$ $\displaystyle=\eta,$
$\displaystyle\chi^{\prime}$
$\displaystyle=Q^{2}\left[1+(\kappa+\varpi)q+\varsigma q^{2}\right].$ (102b)

Hence Eq. (95), the stability condition for the averaged model, becomes

$\displaystyle\dfrac{\vartheta+\iota q}{1+(\kappa+\varpi)q+\varsigma
q^{2}}<\dfrac{\zeta}{1+\mu q}<\dfrac{\eta}{\alpha+\xi q}.$ (103)

####  C.2.2 Equivalence of the simplified models (26), (27), (28) with
$\beta=0$
to a model without phenotypic variation

If $\xi=\varpi=\varsigma=0$ and $\beta=0$, then $q=\gamma/\tau$, where
$\tau=\delta-\mu$. Feasibility requires $\tau>0$. In the cases $\mu=\kappa=0$,
$\iota=\kappa=0$, $\mu=\iota=0$, the interval in Eq. (103) thus reduces to

$\displaystyle\vartheta+\dfrac{\gamma\iota}{\delta}<\zeta<\dfrac{\eta}{\alpha},\;\vartheta<\dfrac{\zeta}{1+\mu\gamma/\tau}<\dfrac{\eta}{\alpha},\;\dfrac{\vartheta}{1+\kappa\gamma/\delta}<\zeta<\dfrac{\eta}{\alpha},$
(104)

respectively, which, from the calculations in Appendix B, are precisely the
stability conditions for Eqs. (26), (27), and (28), respectively, with
$\beta=0$. This establishes the one-to-one correspondence between these
simplified models and the corresponding averaged models claimed in the main
text. We note that, for model (28), this requires an additional condition on
the model parameters discussed in the analysis of that model; this condition
follows from the persister scalings (7).

####  C.2.3 Stability of the averages of the simplified models (26), (27),
(28)
for $\beta\not=0$

If now $\xi=\varpi=\varsigma=0$, but $\beta\not=0$, then $q=(\gamma+\beta
a)/\tau$, where, again, $\tau=\delta-\mu>0$. We discuss the three simplified
models (26), (27), (28) severally, simplifying the interval (103) using the
explicit expressions for $a$ derived in Appendix B.

For model (26), $\mu=\kappa=0$ and $a=a_{\pm}$, defined in Eq. (30a) and
corresponding to the equilibria $\mathcalbf{C}_{\pm}$. The stability
conditions (103) of the averaged model become

$\displaystyle U<0,\;V+\dfrac{\iota\beta}{\delta}a_{\pm}<0,$
$\displaystyle\mbox{with
}U=\alpha\zeta-\eta,\;V=\vartheta-\zeta+\dfrac{\gamma\iota}{\delta},$ (105)

as in definitions (31a). Since $a_{\pm}>0$ by feasibility, a necessary
condition for stability is $U,V<0$, but this is inconsistent with feasibility
of $\mathcalbf{C}_{+}$ (Appendix B), which is therefore unstable in the
averaged model. Now, using the expression for $a_{-}$ in Eq. (30a), we find

$\displaystyle V+\dfrac{\iota\beta}{\delta}a_{-}$
$\displaystyle=\dfrac{1}{2\alpha\delta}\left(Z-\sqrt{Z^{2}+\alpha\delta^{2}UV}\right),$
(106)

where $Z=\delta(\alpha V-U)+\iota\beta$, and infer that the second condition
above is satisfied if $U,V<0$. Since this condition is necessary for
stability, $\mathcalbf{C}_{-}$ is stable in the average of model (26) if and
only if $U,V<0$.

For model (27), $\iota=\kappa=0$ and $a$ is given in Eqs. (53) and corresponds
to the single coexistence equilibrium $\mathcalbf{C}$. The interval (103)
yields the inequalities

$\displaystyle V\left(1+\dfrac{\mu\vartheta\beta}{\tau W}\right)<0,\quad
U\left(1+\dfrac{\mu\vartheta\beta}{\tau W}\right)<0,$ (107)

wherein

$\displaystyle U$
$\displaystyle=\alpha\zeta-\eta-\dfrac{\mu}{\tau}\bigl{(}\gamma\eta+\beta\zeta\bigr{)},$
$\displaystyle V=\vartheta\left(1+\dfrac{\mu\gamma}{\tau}\right)-\zeta,$ (108)

are as in Eqs. (54a), and $W=(\theta U+\eta V)/\zeta$. Feasibility of
$\mathcalbf{C}$ requires $U,V$, and hence $W$ to be of the same sign (Appendix
B). Clearly, $U,V,W>0$ is not compatible with inequalities (107). If
$U,V,W<0$, then they require $1+\mu\beta\vartheta/\tau W<0$, which reduces to
$\vartheta<\eta/\alpha$. Hence coexistence in the averaged model corresponding
to Eqs. (27) is stable (and feasible) if and only if $\vartheta<\eta/\alpha$
and $U,V<0$.

Finally, for model (28), $\iota=\mu=0$ and $a=a_{\pm}$, where $a_{\pm}$ are
defined in Eq. (72a) and correspond to the equilibria $\mathcalbf{C}_{\pm}$.
The stability interval (103) reduces to

$\displaystyle U<0,\;V<\dfrac{\beta\kappa\zeta}{\delta}a_{\pm},\;\text{with
}U=\alpha\zeta-\eta,\;V=\vartheta-\left(1+\dfrac{\gamma\kappa}{\delta}\right)\zeta,$
(109)

as defined in Eqs. (73a), too. Since $a_{\pm}>0$ by feasibility, the second
condition is clearly satisfied if $V<0$. We may therefore suppose that $U<0$
and $V>0$. Now, on letting $Y=\delta(\vartheta U+\eta
V)/\zeta+\beta\kappa\zeta$ as in Eqs. (73b) and using the explicit form of
$a_{\pm}$ given in Eq. (72a), the second condition in Eqs. (109) becomes

$\displaystyle\mp\sqrt{Y^{2}-4\beta\delta\eta\kappa
V}>\dfrac{2\delta\eta}{\zeta}V-Y.$ (110)

We observe that

$\displaystyle\dfrac{2\delta\eta}{\zeta}V-Y<0\quad\Longleftrightarrow\quad\beta>\dfrac{\delta}{\kappa\zeta^{2}}\left(\eta
V-\vartheta U\right).$ (111)

Since we assume that $U<0,V>0$, feasibility of $\mathcalbf{C}_{\pm}$ requires,
from Eq. (75),

$\displaystyle\beta>\dfrac{\delta}{\kappa\zeta^{2}}\left(\sqrt{\eta
V}+\sqrt{-\vartheta U}\right)^{2}>\dfrac{\delta}{\kappa\zeta^{2}}\left(\eta
V-\vartheta U\right),$ (112)

where the second inequality holds since $(x+y)^{2}>x^{2}+y^{2}$ for $x,y>0$.
This shows that inequalities (111) hold true. Moreover,

$\displaystyle\left(Y^{2}-4\beta\delta\eta\kappa
V\right)-\left(\dfrac{2\delta\eta}{\zeta}V-Y\right)^{2}=\dfrac{4UV\delta^{2}\eta\vartheta}{\zeta^{2}}<0,$
(113)

since $U<0$, $V>0$. It follows from inequalities (111) and (113) that
condition (110) holds true. On comparing with the feasibility conditions
derived in Appendix B, we infer that $\mathcalbf{C}_{-}$ is stable in the
averaged model if and only if it is feasible. Moreover, $\mathcalbf{C}_{+}$ is
stable (and feasible) in its averaged model (which is in general different
from that for $\mathcalbf{C}_{+}$) if and only if conditions (75) hold.
Appendix B shows that $\mathcalbf{C}_{+}$ is an unstable equilibrium of Eqs.
(28). We therefore emphasize that, while responsive switching ($\beta>0$) thus
destabilizes the coexistence equilibrium $\mathcalbf{C}_{+}$, coexistence in
the unaveraged Eqs. (28) may still be stable at the other coexistence
equilibrium $\mathcalbf{C}_{-}$.

These results also enable us to compare the simplified models (26), (27), (28)
with responsive switching to averaged models with stochastic switching only:
The latter have an effective switching rate $\gamma^{\prime}=\gamma+\beta a$.
This follows similarly to the correspondence of Eqs. (1) and (2) established
in Sec. II and implies that $q^{\prime}=\gamma^{\prime}/\tau=(\gamma+\beta
a)/\tau$ as above, i.e. conditions (104), with $q$ replaced by this
$q^{\prime}$, are precisely the conditions that we have just analyzed. This
means that the averaged model with stochastic switching only is stable if and
only if the averaged model without phenotypic variation is stable. We stress
however that there is no reason to expect this “complete” correspondence to
hold for the full two-species model (6).

##  APPENDIX D STABILITY OF THE TRIVIAL STEADY STATES OF EQS. (6) AND
PERMANENT COEXISTENCE

While we did not obtain any analytical results bearing on the stability (or
indeed the feasibility) of the coexistence equilibria of Eqs. (6), more
meaningful progress can be made as far as the stability of the trivial steady
states is concerned. These results extend our results for the simplified
models by proving that coexistence is permanent for model (6) whenever these
trivial steady states are all unstable.

We begin by noting that the trivial steady state $\mathcalbf{O}=(0,0,0)$ of
model (6) is always unstable. Next, Eqs. (6) have a trivial steady state
$\mathcalbf{A}=(0,0,\zeta/\eta)$, which is feasible for all parameter values.
The Jacobian of Eqs. (6) evaluated at $\mathcalbf{A}$ is

$\displaystyle\left(\begin{array}[]{ccc}1-(\alpha+\beta)\zeta/\eta-\gamma&\delta&0\\\
\gamma+\beta\zeta/\eta&-\tau-\xi\zeta/\eta&0\\\ &&-\zeta\end{array}\right),$
(114d) wherein, again, $\tau=\delta-\mu$ and entries left blank clearly do not
affect stability; we shall assume that $\tau>0$, consistently with the
persister scalings (7). This Jacobian has one eigenvalue $-\zeta<0$, which
results from the trivial connection $\mathcalbf{O}\rightarrow\mathcalbf{A}$.
Classical stability results [27] imply that $\mathcalbf{A}$ is stable if and
only if the sub-Jacobian
$\displaystyle\left(\begin{array}[]{cc}1-(\alpha+\beta)\zeta/\eta-\gamma&\delta\\\
\gamma+\beta\zeta/\eta&-\tau-\xi\zeta/\eta\\\ \end{array}\right),$ (114g)

has $\text{tr}<0$ and $\det>0$. These conditions are at most quadratic in
$\zeta$ and can therefore be solved to show that

$\displaystyle\text{tr}<0\Longleftrightarrow\zeta>\zeta_{0},$
$\displaystyle\det>0\Longleftrightarrow\zeta<\zeta_{-}\text{ or
}\zeta>\zeta_{+},$ (115)

where

$\displaystyle\zeta_{0}$
$\displaystyle=\dfrac{(1-\gamma-\tau)\eta}{\alpha+\beta+\xi},$ (116a)
$\displaystyle\zeta_{\pm}$
$\displaystyle=\dfrac{\eta}{2(\alpha+\beta)\xi}\left(-z\pm\sqrt{z^{2}+4(\alpha+\beta)\xi(\tau+\gamma\mu)}\right),$
(116b)

with $z=\alpha\tau-\beta\mu-(1-\gamma)\xi$. In particular, $\zeta_{-}<0$.
Also, if $\zeta=\zeta_{0}$, then
$\det=-[1-\gamma-(\alpha+\beta)\zeta_{0}/\eta]^{2}-\delta(\gamma+\beta\zeta_{0}/\eta)<0$,
so $\zeta_{-}<\zeta_{0}<\zeta_{+}$ using the second of Eqs. (115). Hence the
necessary and sufficient condition for stability of $\mathcalbf{A}$ is
$\zeta>\zeta_{+}$.

More importantly, if the sub-Jacobian in Eq. (114g) is unstable, its
eigenvectors are

$\displaystyle\left(-\phi\pm\sqrt{\phi^{2}+4\delta\eta(\beta\zeta+\gamma\eta)},2(\beta\zeta+\gamma\eta)\right),$
(117)

with $\phi=(\alpha+\beta)\zeta-(1+\tau-\gamma)\eta-\zeta\xi$, and so one of
the eigendirections is not feasible. This shows that the geometric properties
of $\mathcalbf{A}$ match those of the corresponding trivial steady state of
Eqs. (26) analyzed in Appendix B.

If any more feasible trivial steady states exist, they are of the form
$\mathcalbf{B}=(b,p,0)$, where $b,p>0$ satisfy the simultaneous equations

$\displaystyle b(1-\gamma-b-\kappa p)+\delta p=p(-\tau-\varpi b-\varsigma
p)+\gamma b=0.$ (118)

The Jacobian of Eqs. (6) evaluated at $\mathcalbf{B}$ is, on simplification
using Eqs. (118),

$\displaystyle\left(\begin{array}[]{ccc}-\delta q-b&\delta-\kappa b&\\\
\gamma-\varpi p&-\gamma/q-\varsigma p&\\\ 0&0&\zeta-\vartheta b-\iota
p\end{array}\right),$ (119d) where $q=p/b$ and, once again, entries left blank
do not affect stability. This Jacobian has one eigenvalue $\zeta-\vartheta
b-\iota p$, associated with perturbations out of the plane $A=0$, which may be
of either sign. We note that the sub-Jacobian
$\displaystyle\left(\begin{array}[]{cc}-\delta q-b&\delta-\kappa b\\\
\gamma-\varpi p&-\gamma/q-\varsigma p\end{array}\right),$ (119g)

has trace and determinant

tr $\displaystyle=-\delta q-b-\gamma/q-\varsigma p<0,$ (120a)
$\displaystyle\det$ $\displaystyle=\delta p(\varsigma q+\varpi)+\gamma
b(1/q+\kappa)+(\varsigma-\kappa\varpi)pb>0,$ (120b)

respectively. We have assumed, in the final line and consistently with the
persister scalings (7), that $\varsigma>\kappa\varpi$. [It is because of the
need for this additional assumption in this argument that we have provided
separate proofs of our permanence result for the simplified models (26), (27),
(28), which do not require such an additional assumption.] It follows that the
remaining eigenvalues at $\mathcalbf{B}$ have negative real parts, and that
the stability of $\mathcalbf{B}$ is determined by the sign of $\zeta-\vartheta
b-\iota p$; determining this sign requires solving Eqs. (118). For our
purposes, it suffices to note that this implies that if $\mathcalbf{B}$ is
unstable, then the direction transverse to $A=0$ is unstable. Hence the
geometric properties of $\mathcalbf{B}$ match those of the corresponding
steady state of Eqs. (26) in Appendix B.

Compared to the discussion in Appendix B, there remains however one more case
to be discussed before we can conclude that coexistence in Eqs. (6) is
permanent if all trivial steady states are unstable by the results of Appendix
B. Indeed, there could be multiple steady states of the form
$\mathcalbf{B}=(b,p,0)$, and so, using the notation introduced in Appendix B,
the limit set $\mathcalbf{L}$ could also be a connection of several such
states. This is however impossible because each of them is stable in the plane
$A=0$. We can therefore conclude that, also in the full model (6), coexistence
is permanent if all trivial steady states are unstable.

## References

  * May [1972] R. M. May, Will a large complex system be stable?, Nature (London) 238, 413 (1972).
  * Serván  _et al._ [2018] C. A. Serván, J. A. Capitán, J. Grilli, K. E. Morrison, and S. Allesina, Coexistence of many species in random ecosystems, Nat. Ecol. Evol. 2, 1237 (2018).
  * Allesina and Tang [2012] S. Allesina and S. Tang, Stability criteria for complex ecosystems, Nature (London) 483, 205 (2012).
  * Mougi and Kondoh [2012] A. Mougi and M. Kondoh, Diversity of interaction types and ecological community stability, Science 337, 349 (2012).
  * Coyte  _et al._ [2015] K. Z. Coyte, J. Schluter, and K. R. Foster, The ecology of the microbiome: Networks, competition, and stability, Science 350, 663 (2015).
  * Grilli  _et al._ [2016] J. Grilli, T. Rogers, and S. Allesina, Modularity and stability in ecological communities, Nat. Commun. 7, 12031 (2016).
  * Grilli  _et al._ [2017a] J. Grilli, G. Barabás, M. J. Michalska-Smith, and S. Allesina, Higher-order competitive interactions stabilize dynamics in competitive network models, Nature (London) 548, 210 (2017a).
  * Butler and O’Dwyer [2018] S. Butler and J. P. O’Dwyer, Stability criteria for complex microbial communities, Nat. Commun. 9, 2970 (2018).
  * Roberts [1974] A. Roberts, The stability of a feasible random ecosystem, Nature (London) 251, 607 (1974).
  * Grilli  _et al._ [2017b] J. Grilli, M. Adorisio, S. Suweis, G. Barabás, J. R. Banavar, S. Allesina, and A. Maritan, Feasibility and coexistence of large ecological communities, Nat. Commun. 8, 14389 (2017b).
  * Gibbs  _et al._ [2018] T. Gibbs, J. Grilli, T. Rogers, and S. Allesina, Effect of population abundances on the stability of large random ecosystems, Phys. Rev. E 98, 022410 (2018).
  * Stone [2018] L. Stone, The feasibility and stability of large complex biological networks: a random matrix approach, Sci. Rep. 8, 8246 (2018).
  * Barron and Galla [2020] J. W. Barron and T. Galla, Dispersal-induced instability in complex ecosystems, Nat. Commun. 11, 6032 (2020).
  * Maynard  _et al._ [2019] D. S. Maynard, C. A. Serván, J. A. Capitán, and S. Allesina, Phenotypic variability promotes diversity and stability in competitive communities, Ecol. Lett. 22, 1776 (2019).
  * Haas  _et al._ [2020] P. A. Haas, N. M. Oliveira, and R. E. Goldstein, Subpopulations and stability in microbial communities, Phys. Rev. Research 2, 022036(R) (2020).
  * Maisonneuve and Gerdes [2014] E. Maisonneuve and K. Gerdes, Molecular mechanisms underlying bacterial persisters, Cell 157, 539 (2014).
  * Harms  _et al._ [2016] A. Harms, E. Maisonneuve, and K. Gerdes, Mechanisms of bacterial persistence during stress and antibiotic exposure, Science 354, aaf4268 (2016).
  * Radzikowski  _et al._ [2017] J. L. Radzikowski, H. Schramke, and M. Heinemann, Bacterial persistence from a system-level perspective, Curr. Opin. Biotechnol. 46, 98 (2017).
  * Kussell and Leibler [2005] E. Kussell and S. Leibler, Phenotypic diversity, population growth, and information in fluctuating environments, Science 309, 2075 (2005).
  * Andersson and Hughes [2014] D. I. Andersson and D. Hughes, Microbiological effects of sublethal levels of antibiotics, Nat. Rev. Microbiol. 12, 465 (2014).
  * Balaban  _et al._ [2004] N. Q. Balaban, J. Merrin, R. Chait, L. Kowalik, and S. Leibler, Bacterial persistence as a phenotypic switch, Science 305, 1622 (2004).
  * Kerr  _et al._ [2002] B. Kerr, M. A. Riley, M. W. Feldman, and B. J. M. Bohannan, Local dispersal promotes biodiversity in a real-life game of rock–paper–scissors, Nature (London) 418, 171 (2002).
  * Czárán  _et al._ [2002] T. L. Czárán, R. F. Hoekstra, and L. Pagie, Chemical warfare between microbes promotes biodiversity, Proc. Natl. Acad. Sci. USA 99, 786 (2002).
  * Kirkup and Riley [2004] B. C. Kirkup and M. A. Riley, Antibiotic-mediated antagonism leads to a bacterial game of rock–paper–scissors _in vivo_ , Nature (London) 428, 412 (2004).
  * Reichenbach  _et al._ [2007] T. Reichenbach, M. Mobilia, and E. Frey, Mobility promotes and jeopardizes biodiversity in rock–paper–scissor games, Nature (London) 448, 1046 (2007).
  * Oliveira  _et al._ [2021] N. M. Oliveira, J. H. R. Wheeler, C. Deroy, S. C. Booth, E. J. Walsh, W. M. Durham, and K. R. Foster, Suicidal chemotaxis in bacteria, bioRxiv:2021.12.21.473623v3 (2021).
  * Murray [2002] J. D. Murray, in  _Mathematical Biology_ , Vol. I (Springer, Berlin, Germany, 2002) 3rd ed.
  * Note [1] We imply, in Eqs. (1) and throughout Sec. II and Appendix A, elementwise multiplication of vectors and rows or columns of matrices by writing the corresponding symbols next to each other, and reserve dots to denote matrix multiplication.
  * Kelsic  _et al._ [2015] E. D. Kelsic, J. Zhao, K. Vetsigian, and R. Kishony, Counteraction of antibiotic production and degradation stabilizes microbial communities, Nature (London) 521, 516 (2015).
  * Hofsteenge  _et al._ [2013] N. Hofsteenge, E. van Nimwegen, and O. K. Silander, Quantitative analysis of persister fractions suggests different mechanisms of formation among environmental isolates of _E. coli_ , BMC Microbiol. 13, 25 (2013).
  * Lewis [2008] K. Lewis, Multidrug tolerance of biofilms and persister cells, in  _Bacterial Biofilms_, Current Topics in Microbiology and Immunology, Vol. 322, edited by T. Romeo (Springer, Berlin, Germany, 2008) pp. 107–131.
  * Brown  _et al._ [2001] L. D. Brown, T. T. Cai, and A. DasGupta, Interval estimation for a binomial proportion, Stat. Sci. 16, 101 (2001).
  * Nam and Blackwelder [2002] J.-m. Nam and W. C. Blackwelder, Analysis of the ratio of marginal probabilities in a matched-pair setting, Stat. Med. 21, 689 (2002).
  * Butler and Waltman [1986] G. Butler and P. Waltman, Persistence in dynamical systems, J. Differ. Equations 63, 255 (1986).
  * Hutson and Vickers [1983] V. Hutson and G. Vickers, A criterion for permanent coexistence of species, with an application to a two-prey one-predator system, Math. Biosci. 63, 253 (1983).
  * Note [2] The name “bacteria”, which we use to have a simple way of referring to phenotype B, is of course not strictly biologically accurate, because the persisters (phenotype P) and competitors (phenotype A) are, biologically, bacteria, too.
  * Note [3] The variables $B$ and $P$ in Eqs. (6) must be nondimensionalized using the same dimensional scalings lest the switching terms in Eqs. (6a) and (6b) become unbalanced. We could have chosen a different scaling for $A$ in Eqs. (6) from that of $B$ and $P$ to set e.g. $\alpha=1$, but equal scalings make comparing competition strengths easier. In the same spirit, we have reused the dimensional scalings from Eqs. (6) to nondimensionalize $B^{\prime}$, $A^{\prime}$ in the averaged model described by Eqs. (90) to avoid having to unravel different dimensional scalings when writing down the consistency conditions expressed by Eqs. (101).
  * Cox  _et al._ [2005] D. A. Cox, J. Little, and D. O’Shea, in  _Using Algebraic Geometry_ , Graduate Texts in Mathematics, Vol. 185 (Springer, Berlin, Germany, 2005) Chap. 2, pp. 26–76, 2nd ed.
  * Yap [2000] C.-K. Yap, in  _Fundamental Problems in Algorithmic Algebra_ (Oxford University Press, Oxford, England, 2000) Chap. 6 and Chap. 7, pp. 141–218.
  * Wood  _et al._ [2013] T. K. Wood, S. J. Knabel, and B. W. Kwan, Bacterial persister cell formation and dormancy, Appl. Environ. Microbiol. 79, 7116 (2013).
  * Zou  _et al._ [2022] J. Zou, B. Peng, J. Qu, and J. Zheng, Are bacterial persisters dormant cells only?, Front. Microbiol. 12, 708580 (2022).
  * Rivoire [2016] O. Rivoire, Informations in models of evolutionary dynamics, J. Stat. Phys. 162, 1324 (2016).
  * Zhang  _et al._ [2020] Z. Zhang, S. Guan, and H. Shi, Information geometry in the population dynamics of bacteria, J. Stat. Mech.: Theory Exp. 2020 (7), 073501.
  * Raynaud and Nunan [2014] X. Raynaud and N. Nunan, Spatial ecology of bacteria at the microscale in soil, PLoS ONE 9, e87217 (2014).
  * Stubbendieck  _et al._ [2016] R. M. Stubbendieck, C. Vargas-Bautista, and P. D. Straight, Bacterial communities: Interactions to scale, Front. Microbiol. 7, 1234 (2016).
  * Parker and Chua [1989] T. S. Parker and L. O. Chua, in  _Practical Numerical Algorithms for Chaotic Systems_ (Springer, New York, NY, 1989) Chap. 3.2, Chap. 5.5, and App. B, pp. 61–66, 118–131, and 305–306.
  * Freedman and Waltman [1985] H. I. Freedman and P. Waltman, Persistence in a model of three competitive populations, Math. Biosci. 73, 89 (1985).
  * Butler and Waltman [1988] G. J. Butler and P. Waltman, Persistence in three-dimensional Lotka–Volterra systems, Math. Comput. Model. 10, 13 (1988).
  * Hofbauer and Sigmund [1987] J. Hofbauer and K. Sigmund, Permanence for replicator equations, in  _Dynamical Systems_, Lecture Notes in Economics and Mathematical Systems, Vol. 287, edited by A. B. Kurzhanski and K. Sigmund (Springer, Berlin, Germany, 1987) pp. 70–92.
  * Glendinning [1994] P. Glendinning, in  _Stability, Instability and Chaos: An Introduction to the Theory of Nonlinear Differential Equations_ (Cambridge University Press, Cambridge, UK, 1994) Chap. 5, pp. 102–144.


Present address: ]Theory and Computation of Energy Materials (IEK-13),
Institute of Energy and Climate Research, Forschungszentrum Jülich GmbH, 52425
Jülich, Germany

# Accounting for the Quantum Capacitance of Graphite in Constant Potential
Molecular Dynamics Simulations

Kateryna Goloviznina  Sorbonne Université, CNRS, Physicochimie des
Électrolytes et Nanosystèmes Interfaciaux, F-75005 Paris, France  Réseau sur
le Stockage Electrochimique de l’Energie (RS2E), FR CNRS 3459, 80039 Amiens
Cedex, France  Johann Fleischhaker  Institute of Polymers and Composites,
Hamburg University of Technology, 21073 Hamburg, Germany  Sorbonne Université,
CNRS, Physicochimie des Électrolytes et Nanosystèmes Interfaciaux, F-75005
Paris, France  Tobias Binninger11footnotemark: 1 ICGM, Univ Montpellier, CNRS,
ENSCM, Montpellier, France  [  Benjamin Rotenberg  Sorbonne Université, CNRS,
Physicochimie des Électrolytes et Nanosystèmes Interfaciaux, F-75005 Paris,
France  Réseau sur le Stockage Electrochimique de l’Energie (RS2E), FR CNRS
3459, 80039 Amiens Cedex, France  Heigo Ers  University of Tartu, Ravila 14a,
51004 Tartu, Estonia  Vladislav Ivanistsev  University of Tartu, Ravila 14a,
51004 Tartu, Estonia  Robert Meissner  Institute of Polymers and Composites,
Hamburg University of Technology, 21073 Hamburg, Germany  Institute of Surface
Science, Helmholtz-Zentrum Hereon, 21502 Geesthacht, Germany  Alessandra Serva
Sorbonne Université, CNRS, Physicochimie des Électrolytes et Nanosystèmes
Interfaciaux, F-75005 Paris, France  Réseau sur le Stockage Electrochimique de
l’Energie (RS2E), FR CNRS 3459, 80039 Amiens Cedex, France  Mathieu Salanne
<EMAIL_ADDRESS>Sorbonne Université, CNRS,
Physicochimie des Électrolytes et Nanosystèmes Interfaciaux, F-75005 Paris,
France  Réseau sur le Stockage Electrochimique de l’Energie (RS2E), FR CNRS
3459, 80039 Amiens Cedex, France  Institut Universitaire de France (IUF),
75231 Paris, France

###### Abstract

Molecular dynamics simulations at a constant electric potential are an
essential tool to study electrochemical processes, providing microscopic
information on the structural, thermodynamic, and dynamical properties.
Despite the numerous advances in the simulation of electrodes, they fail to
accurately represent the electronic structure of materials such as graphite.
In this work, we introduce a simple parameterization method that allows to
tune the metallicity of the electrode based on a quantum chemistry calculation
of the density of states. As a first illustration, we study the interface
between graphite electrodes and two different liquid electrolytes, an aqueous
solution of NaCl and a pure ionic liquid, at different applied potentials. We
show that the simulations reproduce qualitatively the experimentally-measured
capacitance; in particular, they yield a minimum of capacitance at the point
of zero charge, which is due to the quantum capacitance contribution. An
analysis of the structure of the adsorbed liquids allows to understand why the
ionic liquid displays a lower capacitance despite its large ionic
concentration. In addition to its relevance for the important class of
carbonaceous electrodes, this method can be applied to any electrode materials
(e.g. 2D materials, conducting polymers, etc), thus enabling molecular
simulation studies of complex electrochemical devices in the future.

##  I Introduction

The electrochemical double-layer (EDL) is the region where the electrolyte
contacts the electrode in an electrochemical device. Despite its importance in
many fields, ranging from energy storage, [1] electrocatalysis, [2] sensing,
[3] water desalination, nanogenerators [4], etc., much remains to be known
about its structure and properties. In situ experiments provide many useful
details, [5] but they are often hindered by their spatial resolution since the
processes at play occur in a liquid region of a few nanometers only. From the
theoretical point of view, although textbooks still present the
Gouy–Chapman–Stern model as a realistic picture, its relevance for all these
applications has largely been questioned due to its drastic approximations.
[6, 7] Based on these observations, microscopic simulations appear as the
ideal tool for advancing our knowledge about the EDL. [8]

Yet, the application of microscopic simulations to such complicated
interfacial systems is not straightforward. Researchers are faced with the
conundrum of simulating large-scale systems, over relatively long timescales,
while keeping the accuracy of the calculation to the highest level. For
example, the characteristic adsorption lifetimes of ionic species inside
electrified nanoporous electrodes are of several nanoseconds. [9] Recent
developments have made density functional theory (DFT) accessible to calculate
the interactions in molecular dynamics (MD) simulation, but the associated
computational costs are very large. [10] This method was thus so far limited
to the important case of metal-water interfaces, [11] with a special focus on
the reactivity of water molecules. [12] In recent years, several works have
intended to account for the presence of ions by including them directly within
the electrochemical double-layer. This approach yielded accurate estimates of
the interfacial capacitance at the interfaces between gold, [13] copper or
graphene [14] and simple aqueous solutions, but the limited system size
prevents one from controlling the concentration of ions. A similar setup was
used in the framework of electrocatalysis applications to simulate the effect
of various cations on the CO2 reduction reaction mechanism. [15, 16] The main
example of non-aqueous system simulated using DFT-based MD is the study of the
potential drop at the interface formed between a pure ionic liquid and a
graphene electrode, [17] but it is worth noting that the computational cost of
the simulations was of the order of 1 million CPU hours to obtain trajectories
of $\approx$1 ps. For applications in which the effects of ion adsorption are
important and an efficient sampling of the liquid configurations needs to be
performed, the use of less costly alternative is therefore necessary.

Consequently, the main method currently employed for simulating
electrochemical double-layers is classical MD, in which the interactions are
represented using a classical force field. [8] In current electrode models,
the polarization by the electrolyte is generally accounted for by allowing the
partial atomic charges of the electrode atoms to fluctuate over time, while
maintaining a constant potential difference between the negative and the
positive electrodes. [18, 19] However, this approach did not include the so-
called quantum capacitance (QC) effects. As it has been known since the
pioneering work of Gerischer on graphite, [20] the density of states (DOS) of
the electrode may vary near the Fermi level, which reflects on the electrode
ability to respond to the polarization by the electrolyte. A simple solution
to this problem was found by exploiting the fact that the EDL and quantum
capacitance contributions are in series, and could thus be obtained by two
different calculations. [21, 22] The QC can then be readily obtained at the
DFT level from a static calculation; a recent example of this approach is the
determination of the capacitance of graphdiynes, a new class of 2D porous
materials that display a complex electronic structure. [23] While it gives
accurate capacitances, this approach cannot yield a good structural
description of the EDL at a given potential since the interaction of the
electrolyte with the electrode is not modified accordingly in the
computational workflow. Recently, several attempts to account for the
metallicity of the electrode were reported, in particular in the framework of
the Thomas–Fermi model. [24, 25] The constant-potential version of this model
includes an additional parameter, the Thomas–Fermi length, that accounts for
the screening of the potential within the electrode. It was used to
systematically study the effect of the electrode metallicity on properties
such as the interfacial free energies. [26] In principle, the Thomas–Fermi
model includes intrinsically the QC effects, yet the use of a single parameter
that is material-dependent does not allow to investigate the impact of the
variations of the number of charge carriers in the electrode with the applied
potential.

In this work, we show how to account explicitly for the DOS of an electrode
material in classical molecular dynamics simulations. To this end, we selected
graphite as a representative material due to its ubiquitous importance and to
its peculiar DOS in the vicinity of the Fermi level. The parameterization
requires a simple DFT calculation, and the model can be used in the current
implementations of the constant potential method. We show its application by
computing the full EDL capacitance, which now contains the quantum
contribution, for two different electrolytes, aqueous NaCl (with a
concentration of 1 M) and the 1-ethyl-3-methylimidazolium tetrafluoroborate
ionic liquid ([EMIM][BF4] IL). We show that the agreement with experiments is
almost quantitative, in particular our simulations reproduce the minimum of
capacitance at the point of zero charge (PZC). This then allows us to discuss
the microscopic details of the systems; in particular we show that despite its
much lower ionic concentration the aqueous system displays a larger
capacitance due to a strong contribution from the water molecules. In the
future, this approach will enable the study of electrode materials with more
complex electronic structure, in particular in the context of capacitive
energy storage.[27, 23, 28]

##  II Model

###  II.1 The electrode model and its parameterization

Figure 1:  (a-c) Illustration of the impact of the Gaussian width ($\sigma$)
on the charge density of a coronene molecule. For visualization purposes, two-
dimensional charge distributions centered on the carbon atoms were used in
this plot. The values in the colorbar are given in $e$/Å2. The represented
charge density corresponds to a case where all the partial charges of the
carbon atoms are equal; in the constant potential model, they fluctuate in
response to the electrolyte. (d) Density of states of graphite obtained from
DFT calculations

The constant potential electrode model introduced by Siepmann and Sprik [18]
is based on a representation of the electrode charge distribution, $\rho_{\rm
elec}$, as a collection of atom-centered Gaussians:

$\rho_{\rm elec}(\textbf{r})=\sum_{i=1}^{N_{\rm
elec}}q_{i}\dfrac{1}{(2\pi\sigma^{2})^{3/2}}\exp\Bigg{(}-\dfrac{\mid\textbf{r}-\textbf{r}_{i}\mid^{2}}{2\sigma^{2}}\Bigg{)},$
(1)

where $N_{\rm elec}$ is the number of atoms (with positions
${\textbf{r}_{i}}$) in the electrodes and $\sigma$ is the width of the
Gaussian. The partial charges $q_{i}$ of the electrode atoms behave as
additional degrees of freedom of the simulation, contrary to the electrolyte
charges represented using point charges of fixed magnitude. The fluctuating
charges should take values such that the potential is uniform within the whole
electrode at each timestep. In practice, we define two electrodes with
potentials $\Psi_{+}$ and $\Psi_{-}$ (and with the corresponding number of
atoms $N_{+}$ and $N_{-}$). The determination of the fluctuating charges is
done within the Born–Oppenheimer approximation, by minimizing the total
electrostatic potential energy of the system which includes all the Coulomb
interactions, as well as the electric work which is necessary to charge the
electrodes to the desired potential. The corresponding expressions are
provided in Reference 29.

In the original model, this potential energy term includes only one parameter,
which is the width of the Gaussian (taken the same for all the electrode
atoms). Tuning its value has two effects. Firstly, larger Gaussians lead to an
increased delocalization of the electrode charge. This effect can be
appreciated in the example of two-dimensional charge distributions for a
coronene molecule in Figure 1a-c, in the specific case for which all the atoms
would carry the same partial charge. For example, for
$\sigma=$0.207\text{\,}\mathrm{\text{Å}}$$, it is highly localized around the
atoms. Some overlap starts to be seen for a value of
$0.486\text{\,}\mathrm{\text{Å}}$, as a non-negligible charge density starts
to be observed along the C-C bonds. Finally, for a larger value of
$0.764\text{\,}\mathrm{\text{Å}}$, the charge density is now almost equally
distributed within the carbon rings. The latter representation is in agreement
with pictures extracted from single-molecule imaging of aromatic carbonaceous
molecules using scanning tunneling microscopy or atomic force microscopy. [30]
Secondly, larger Gaussian widths result in an enhanced ability to charge the
atom. Indeed, the Gaussian charges self-interact through the term:

$U_{\rm
self}=\frac{1}{4\pi\epsilon_{0}}\left(\frac{1}{\sigma\sqrt{4\pi}}\sum_{i\in(N_{+}+N_{-})}q_{i}^{2}\right),$
(2)

so that the magnitude of this self-interaction corresponds to an energy
penalty for charging the electrode.

The extension of the model to include Thomas–Fermi effects simply introduced
an additional energy penalty proportional to the Thomas–Fermi screening length
$l_{\rm TF}$:

$U_{\rm TF}=\frac{1}{4\pi\epsilon_{0}}\left(2\pi l_{\rm
TF}^{2}\rho\sum_{i\in(N_{+}+N_{-})}q_{i}^{2}\right)$ (3)

where $\rho$ is the number density of electrode atoms within the material.
Note that this term does not impact the interactions between neighbor atoms
and the delocalization of the charge within the plane, which remained
controlled by the Gaussian width. Quantities such as the interfacial
capacitance or the solid-liquid interfacial free energies were shown to vary
systematically with $l_{\rm TF}$ in previous works. [24, 26]

Here our objective is to go further and use the DOS of the material
($\mathcal{D}(E)$) to parameterize the electrode model. The DOS is defined as
a measure of the number of allowed electronic states whose energies lie at a
value $E$. A DOS plot of graphite, calculated at the DFT level (details given
in the SI), is shown in Figure 1d. It displays a characteristic minimum around
the Fermi level $E_{\mathrm{F}}$, which often leads graphite to be classified
as a “semi-metal”. To accurately describe the quantum capacitance effect, the
residual value of the DOS at the Fermi level must be correctly quantified,
requiring the use of a very dense $K$-point mesh and a consistent treatment of
thermal broadening in the calculation of the DOS (see SI for details) [31].
The Thomas–Fermi length can be expressed in terms of the DOS at the Fermi
level: $l_{\rm
TF}=\Bigg{(}\dfrac{e^{2}}{\epsilon_{0}}\mathcal{D}(E_{F})\Bigg{)}^{-1/2}$
which provides an easy parameterization of the model through Equation 3. This
approach was successfully used in a theoretical model to explain the
difference in freezing properties behavior for an ionic liquid in contact with
electrodes of varying metallicities. [32]

Here we extend this approach by using the fact that in the fixed band
approximation, the Fermi level is shifted by the bias potential to
$E_{F}+e\Psi$, which leads to an increase in the number of available states.
We use this correspondence to determine the DOS at a given electric potential,
which in turn gives a new value for the Thomas–Fermi screening length.
However, a difficulty arises in choosing the bias potential since in the
constant potential electrode model; the fixed quantity is the applied
potential ${\Delta\Psi}$ between the electrodes. A constraint of global
electroneutrality is imposed, [33] so that for ${\Delta\Psi}=0$ the two
electrodes are equivalent. Their charge must fluctuate around zero, and we can
define their potential as the PZC ($\Psi_{+}=\Psi_{-}=\Psi_{\rm PZC}$). For
non-null applied voltages, there exists a potential shift that depends on the
instantaneous configurations and fluctuates around an average value, so that
$\Psi_{+}=\frac{\Delta\Psi}{2}+\Psi_{\rm shift}$ and
$\Psi_{-}=-\frac{\Delta\Psi}{2}+\Psi_{\rm shift}$. In the following, based on
the symmetric shape of the experimental capacitance, we arbitrarily neglect
the variation of the electrode potentials ($\Psi_{\rm shift}=0$ V, but in
future work, it should be possible to allow the Thomas–Fermi length to vary
dynamically to account for it).

In practice, we chose to use the Gaussian width as the main parameter to
control the DOS of the simulation instead of the $l_{\rm TF}$. This leads to
the following parameterization for $\sigma$ depending on the potential:

$\sigma(\Psi)=\frac{1}{4\pi\epsilon_{0}}\left(\frac{e^{2}\mathcal{D}(E_{F}+e\Psi)}{\rho\sqrt{\pi}}\right)$
(4)

The obtained values of Gaussian width are given in Supplementary Table S1. At
the PZC, the Gaussian width takes a minimal value of 0.207 Å which corresponds
to the situation shown in Figure 1a for coronene. Since the DOS slightly
differs depending on the sign of the potential, we considered non-identical
Gaussian width for the negative and positive electrodes. However, for the
largest applied potential of 2 V, the corresponding Gaussian width would be
larger than 1.0 Å. This value corresponds to a situation where the overlap
between neighboring Gaussians is too large, and the problem becomes ill-
defined mathematically. [34] In practice, we observed that the chosen $\sigma$
cannot be larger than $0.764\text{\,}\mathrm{\text{Å}}$, which corresponds to
the case shown on Figure 1c for coronene. This value was thus selected for the
2 V potential simulation for both the positive and negative electrode atoms.

###  II.2 Non-Coulombic electrode/electrolyte interactions

Additional parameters need also to be defined on the electrolyte side.
Although, for the aqueous solution, a well-established non-polarizable force
field is employed, it is not the case for the [EMIM][BF4] ionic liquid, for
which we use a polarizable interaction potential. The induced dipoles are
calculated at each timestep together with the fluctuating charges, using a
preconditioned conjugate gradient method as described in Reference 35. In such
a case, the parameters of the conventional OPLS force field cannot be used in
order to avoid double-counting of some interactions as explained briefly in
the following. We therefore use the CL&Pol approach,[36, 37] which consists in
scaling down the non-bonded attractive term (in the case of the Lennard Jones
potential, the $\epsilon$ parameter) by the factor $k_{ij}$:

$k_{ij}=\frac{E_{\mathrm{disp}}}{E_{\mathrm{disp}}+E_{\mathrm{ind}}}$ (5)

so that the potential accounts only for dispersion interactions, while the
implicit induction contribution is removed. The scaling factor can be
evaluated from the symmetry-adapted perturbation theory (SAPT) calculation,
which allows isolating the dispersion and induction terms from the total
interaction energy of two fragments. IL ions are treated as entire fragments,
while a coronene molecule is considered as representative of the carbon
nanomaterial. The geometries of corresponding dimers (cation–coronene,
anion–coronene, and coronene–coronene) are optimized, and then potential
energy scans are performed at different SAPT levels (details given in SI). An
example of the SAPT scan is given in Supplementary Figure S2, and the final
$k_{ij}$ values are reported in Supplementary Table S2.

Figure 2: Normalised density profiles of selected atoms of [EMIM][BF4] as a
function of their distance from the electrode. Rose and blue lines correspond
to the trajectories of the same length, using the same initial configuration
($\sim$ $4\text{\,}\mathrm{ps}$), while the orange line corresponds to a long
$2\text{\,}\mathrm{ns}$ run.

In order to validate the proposed model, we compared the atomic forces (given
in Supplementary Figure S3) and density profiles to the ones of a DFT-based MD
simulation. We re-used a model (Ref. 17), which represents a single graphite
electrode with $3$ layers and $200$ ion pairs of [EMIM][BF4]. As discussed
above, due to the high computational cost associated with the DFT-MD method
(see methods description in SI), the total simulation time was only around 5
ps. Nevertheless, it allowed us to check whether the equilibrium distances
between the ions and the electrodes are well reproduced. An example of density
profiles is given in Figure 2. Both in the cases of the N atom of the EMIM+
cation and the B atoms of the BF4- anion, we observe an excellent agreement
between two methods in the peak positions and intensities up to
$\sim$$15\text{\,}\mathrm{\text{Å}}$ from the surface. The positions of the
two first cationic and anionic layers are also preserved when a long
$2\text{\,}\mathrm{ns}$ run is carried out. The intensities of these peaks
then decrease as expected from the better sampling of the configuration space.

###  II.3 Results and discussion

As a first application of our method, we focus on the canonical example of the
graphite electrode. By subtracting the EDL capacitance (estimated using
analytical contributions associated with the Helmholtz and the diffuse layer)
from the experimentally measured differential capacitance, in 1985, Gerischer
determined the QC, and hence the DOS around the Fermi level, in qualitative
agreement with the theoretical calculations available at the time. [20] Here
we use an inverse approach, in which the DOS is directly fed into the model,
with the aim to determine the differential capacitance without any empirical
input (other than the one used to build the force fields).

When using the constant potential method, a direct output of the simulation is
the total charge ($Q_{\rm elec}=\sum_{i\in N_{\rm elec}}q_{i}$) on the
electrodes. An electroneutrality constraint is imposed so that
$Q_{+}=-Q_{-}=Q$ at each timestep. In order to determine the capacitance, it
is necessary to extract the average value of $Q$ once the equilibrium has been
reached at a given potential. As shown in Supplementary Figure S4, this regime
is obtained well before a time of 1 nanosecond of simulation for the aqueous
solution, but several nanoseconds are necessary in the case of the ionic
liquid due to its much larger viscosity.

Figure 3: (a) Average surface charge on the graphite with respect to the
electrode potential (symbol: simulation data, dashed line: fit with a third-
order polynomial function). The blue stars are the surface charges obtained
using a constant Gaussian width of 0.207 Å, i.e., the Gaussian value obtained
from the DOS at $0\text{\,}\mathrm{V}$. (b) Differential capacitance
dependence on the electrode potential. The experimental values are extracted
from references 38 and 39.

The differential capacitance is then calculated by deriving the variation of
$<Q>$ with respect to the electrode potential, $\Psi_{+}$ or $\Psi_{-}$. In
practice, due to the limited number of points, a fit of the $Q=f(\Psi)$ plot
is necessary prior to the derivation. As can be seen in Figure 3a, a third-
order polynomial function yields a very good fit of the data. Comparison with
the experimental results obtained for the same electrolytes [38, 39], shown in
Figure 3b as dotted lines, gives qualitative agreement. In particular, our
simulations reproduce the minimum of capacitance at the point of zero charge
(which corresponds by construction to an applied potential of 0 V in our
case). Still, they overestimate the value of this minimum by approximately 2
µF cm-2. This discrepancy may be attributed to either some errors due to the
employed force field or to the fact that we simulate a perfect, pristine
graphite surface while the experimental materials may not be perfectly planar
despite the use of highly ordered pyrolytic graphite. Despite this difference,
the variation of the differential capacitance with potential is very similar
between simulations and experiments, which shows that the response to applied
voltage is well taken into account. In order to show further the impact of the
QC in the simulation, additional simulations of the aqueous NaCl solution were
performed with a fixed Gaussian width of 0.207 Å (that corresponds to the
value taken at rest in our scheme). The accumulated charge is then much
smaller, as can be seen in Figure 3a (blue stars), and its variation is better
fitted with a linear function, which corresponds to a constant interfacial
capacitance.

Notably, as in experiments, when accounting for the QC we observe that the
capacitance of the aqueous solution is larger than the one of the ionic
liquid. This result may be surprising since the latter is entirely composed of
ions, and could thus be expected to induce a larger charge at the surface of
the electrode. It is, therefore, interesting to compare the local structure of
the two systems and try to understand the microscopic mechanism at play. The
structure of the two interfaces was largely documented in various simulation
studies,[40, 41, 42, 43] but in the absence of explicit QC effect and no
systematic comparison was performed so far. On the one hand, in the ionic
liquid, at zero applied electric potential, a “checkerboard-like” structure is
obtained in which the first liquid layer is made of well-organized adsorbed
ions. In particular, as shown in Figure 4, the imidazolium rings tend to lie
flat on the surface[44] and to arrange collectively along lines. This is
consistent with previous observations based on atomic force microscopy
experiments as well as previous MD simulations on similar systems. [45, 46]
Upon application of a voltage, due to steric hindrance, the charge can only
accumulate via an exchange mechanism. In the case of the negative electrode
and for an applied voltage of $2\text{\,}\mathrm{V}$, the cations also tend to
reorient and become perpendicular to the electrode (Supplementary Figure S5),
thus allowing a larger ionic density on the surface.

Figure 4: Representative snapshots of the adsorbed layer for [EMIM][BF4] IL-
based simulations: EMIM^+ is given is red, BF4- in blue. The surface atoms are
colored according to their charge (white: neutral, blue: negative, red:
positive).
Figure 5: Representative snapshots of the adsorbed layer in the case of
$1\text{\,}\mathrm{M}$ aqueous NaCl (black and white: water molecules, red
spheres: Na+ ions, blue sphere: Cl- ions). The surface atoms are colored
according to their charge (white: neutral, red: positive, blue: negative). The
represented Na+ and Cl- ions are located beyond the plane of water molecules,
and they do not induce any noticeable surface charge on the electrode.

On the contrary, in the case of the aqueous sodium chloride solution, we do
not observe any adsorption of the Na+ and Cl- in the first liquid layer at low
potential (see Supplementary Figure S6). Some ions are present in the second
layer, but they do not lead to any noticeable surface charge on the electrode
surface, as shown in Figure 5. This is consistent with previous simulations
performed on gold surfaces, in which it was shown that the degree of
metallicity of the gold had to be very high to promote their adsorption. [47]
Similarly, no specific adsorption of Na+ and Cl- ions was observed in a
previous work when the Thomas–Fermi length was systematically varying.[24] The
polarization of the surface is thus almost entirely due to the behavior of the
liquid water molecules. At zero voltage, they adopt a parallel orientation
with respect to the electrode surface, and they form a 2-dimensional hydrogen
bond network inside the plane. Up to one volt, the charging of the electrode
does not imply any variation on the composition of the EDL, but it is rather
accompanied by a change of the orientation of the water molecules, as shown in
Supplementary Figure S7. In line with our results, a recent simulation study
of aqueous NaCl at different concentrations between Pt electrodes has shown
that the electrostatic potential of the EDL is almost independent on the ionic
concentration due to the molecular influence of water molecules.[48]

Based on the different charging mechanisms, it is now possible to interpret
the differences in capacitance between the two systems. From the energy point
of view, charging the ionic liquid EDL requires breaking some strong Coulombic
interactions between ions of opposite signs, while in the aqueous NaCl case,
it only needs to rearrange the hydrogen bond network, which is substantially
easier. Interestingly, this interpretation holds in the case of a planar
surface, but it cannot be applied to more complex cases such as porous carbons
for which the confinement effects lead to (i) an increased screening of the
Coulombic interaction, which was termed as a “superionic” effect by Kondrat
and Kornyshev [49, 50] and confirmed by in situ X-ray diffraction experiments;
[51] and (ii) the impossibility to induce a macroscopic polarization through
the water dipoles only, for symmetry reasons. Consequently, in the latter
case, the polarization occurs via the adsorption/desorption of the ions, and
the capacitances are not larger anymore than for ionic liquids. [52]

It is worth noting that the mechanisms described above are valid over the
whole investigated potential range, and so there is interference from the QC
on the overall behavior: It changes the onset potential at which the ion
exchange/water reorientations occur. This validates the approaches based on a
decoupling of the EDL and QC effect proposed in previous works, [21, 22, 23]
but our approach has the advantage of specifying the liquid structure at each
potential without any ambiguity.

##  III Conclusions & perspectives

In conclusion, in this work we have shown how the constant potential method
can be extended to explicitly account for the DOS of the material in MD
simulations of electrochemical double layers. This extension can readily be
used in any existing implementation, since it only requires fine-tuning of the
Gaussian width of the electrode charge depending on the applied potential. The
parameterization is easily made from a calculation of the DOS of the material
at rest, using first-principles methods. Our new method, therefore, accounts
directly for the QC contribution to the overall capacitance, contrarily to
previous works in which it was simply added to the double-layer contribution.
The advantage is that it is possible to extract the correct structure of the
liquid at the interface for a given applied potential.

As a first application of the method, we focused on the graphite electrode.
Despite its importance in energy applications and its apparent simplicity, its
peculiar electronic structure is characterized by a minimum of the DOS at the
Fermi level. Consequently, a correct description of the double-layer forming
at its interface with electrolytes remained lacking so far. In this work, we
obtained a qualitative agreement with experimental works for two important
electrolytes, an aqueous NaCl solution and a pure ionic liquid. The
simulations predict overall too large capacitances (by 1 to 2 µF cm-2) for the
two systems, but the variation with potential is well reproduced in both
cases. The larger capacitance obtained for the aqueous system is interpreted
based on the structure of the adsorbed liquid and its effect on the surface-
induced charges.

In the future, this work could expand the scope of constant potential MD
simulations to a variety of new materials for energy applications. In
particular, many porous materials, such as carbon-based materials or metal-
organic-frameworks, covalent organic frameworks, conducting polymers, [27]
etc, display a complex electronic structure, that cannot be approximated to
the one of a perfect conductor, and could benefit from such a representation.
Additional work is however required to discriminate the presence of different
chemical elements (with eventually different environments) in such materials.
This could be tackled by coupling the present model with machine-learning-
based approaches such as the PiNNwall interface [53], which was recently
applied to the cases of chemically doped graphene and graphene oxide
electrodes.

###### Acknowledgements.

This work was funded by the European Research Council (ERC) under the European
Union Horizon 2020 research and innovation program (grant agreements 771294
and 863473). It was also supported by the Estonian–French cooperation program
Parrot, funded by the Estonian Research Council and Campus France
Computational tasks were performed using resources from GENCI-IDRIS (Grant
A0140910463). DFT-based MD simulations were run at the HPC resources provided
by the Partnership for Advanced Computing in Europe (PRACE).

##  IV References

## References

  * Salanne  _et al._ [2016] M. Salanne, B. Rotenberg, K. Naoi, K. Kaneko, P.-L. Taberna, C. P. Grey, B. Dunn, and P. Simon, “Efficient Storage Mechanisms for Building Better Supercapacitors,” Nat. Energy 1, 16070 (2016).
  * Ledezma-Yanez  _et al._ [2017] I. Ledezma-Yanez, W. D. Z. Wallace, P. Sebastián-Pascual, V. Climent, J. M. Feliu, and M. T. Koper, “Interfacial water reorganization as a ph-dependent descriptor of the hydrogen evolution rate on platinum electrodes,” Nat. Energy 2, 17031 (2017).
  * Chipangura  _et al._ [2023] Y. E. Chipangura, B. D. Spindler, P. Bühlmann, and A. Stein, “Design criteria for nanostructured carbon materials as solid contacts for ion-selective sensors,” Adv. Mater. 36, 2309778 (2023).
  * Pace  _et al._ [2023] G. Pace, A. Esau del Rio Castillo, A. Lamperti, S. Lauciello, and F. Bonaccorso, “2d materials-based electrochemical triboelectric nanogenerators,” Adv. Mater. 35, 2211037 (2023).
  * Forse  _et al._ [2016] A. C. Forse, C. Merlet, J. M. Griffin, and C. P. Grey, “New perspectives on the charging mechanisms of supercapacitors,” J. Am. Chem. Soc. 138, 5731–5744 (2016).
  * Kornyshev [2007] A. A. Kornyshev, “Double-layer in ionic liquids: paradigm change?” J. Phys. Chem. B 111, 5545–5557 (2007).
  * Wu [2022] J. Wu, “Understanding the electric double-layer structure, capacitance, and charging dynamics,” Chem. Rev. 122, 10821–10859 (2022).
  * Jeanmairet, Rotenberg, and Salanne [2022] G. Jeanmairet, B. Rotenberg, and M. Salanne, “Microscopic simulations of electrochemical double-layer capacitors,” Chem. Rev. 122, 10860–10898 (2022).
  * Pean  _et al._ [2015] C. Pean, B. Daffos, B. Rotenberg, P. Levitz, M. Haefele, P.-L. Taberna, P. Simon, and M. Salanne, “Confinement, desolvation, and electrosorption effects on the diffusion of ions in nanoporous carbon electrodes,” J. Am. Chem. Soc. 137, 12627 (2015).
  * Groß [2023] A. Groß, “Challenges for ab initio molecular dynamics simulations of electrochemical interfaces,” Curr. Opin. Electrochem. 40, 101345 (2023).
  * Groß and Sakong [2022] A. Groß and S. Sakong, “Ab initio simulations of water/metal interfaces,” Chem. Rev. 122, 10746–10776 (2022).
  * Surendralal  _et al._ [2018] S. Surendralal, M. Rodorova, M. W. Finnis, and J. Neugebauer, “First-principles approach to model electrochemical reactions: Understanding the fundamental mechanisms behind mg corrosion,” Phys. Rev. Lett. 120, 246801 (2018).
  * Li  _et al._ [2019] C.-Y. Li, J.-B. Le, Y.-H. Wang, S. Chen, Z.-L. Yang, J.-F. Li, J. Cheng, and Z.-Q. Tian, “In situ probing electrified interfacial water structures at atomically flat surfaces,” Nat. Mater. 18, 697–701 (2019).
  * Li  _et al._ [2023] X.-Y. Li, X.-F. Jin, X.-H. Yang, X. Wang, J.-B. Le, and J. Cheng, “Molecular understanding of the helmholtz capacitance difference between cu(100) and graphene electrodes,” J. Chem. Phys. 158, 084701 (2023).
  * Monteiro  _et al._ [2021] M. C. O. Monteiro, F. Dattila, B. Hagedoorn, R. García-Muelas, N. López, and M. T. M. Koper, “Absence of CO2 electroreduction on copper, gold and silver electrodes without metal cations in solution,” Nat. Catal. 4, 654–662 (2021).
  * Qin, Vegge, and Hansen [2023] X. Qin, T. Vegge, and H. A. Hansen, “Cation-coordinated inner-sphere CO2 electroreduction at au–water interfaces,” J. Am. Chem. Soc. 145, 1897–1905 (2023).
  * Ers  _et al._ [2020] H. Ers, M. Lembinen, M. Mišin, A. P. Seitsonen, M. V. Fedorov, and V. B. Ivaništšev, “Graphene–ionic liquid interfacial potential drop from density functional theory-based molecular dynamics simulations,” J. Phys. Chem. C 124, 19548–19555 (2020).
  * Siepmann and Sprik [1995] J. I. Siepmann and M. Sprik, “Influence of Surface-Topology and Electrostatic Potential on Water Electrode Systems,” J. Chem. Phys. 102, 511–524 (1995).
  * Reed, Lanning, and Madden [2007] S. K. Reed, O. J. Lanning, and P. A. Madden, “Electrochemical Interface Between an Ionic Liquid and a Model Metallic Electrode,” J. Chem. Phys. 126, 084704 (2007).
  * Gerischer [1985] H. Gerischer, “An interpretation of the double layer capacity of graphite electrodes in relation to the density of states at the Fermi level,” J. Phys. Chem. 89, 4249–4251 (1985).
  * Kornyshev, Luque, and Schmickler [2014] A. A. Kornyshev, N. B. Luque, and W. Schmickler, “Differential capacitance of ionic liquid interface with graphite: the story of two double layers,” J. Solid State Electrochem. 18, 1345–1349 (2014).
  * Pak, Paek, and Hwang [2013] A. J. Pak, E. Paek, and G. S. Hwang, “Relative contributions of quantum and double layer capacitance toward the supercapacitor performance of carbon nanotubes in an ionic liquid,” Phys. Chem. Chem. Phys. 15, 19741–19747 (2013).
  * Mo  _et al._ [2023] T. Mo, Z. Wang, L. Zeng, M. Chen, A. A. Kornyshev, M. Zhang, Y. Zhao, and G. Feng, “Energy storage mechanism in supercapacitors with porous graphdiynes: Effects of pore topology and electrodemetallicity,” Adv. Mater. 35, 2301118 (2023).
  * Scalfi  _et al._ [2020a] L. Scalfi, T. Dufils, K. G. Reeves, B. Rotenberg, and M. Salanne, “A semiclassical Thomas–Fermi model to tune the metallicity of electrodes in molecular simulations,” J. Chem. Phys. 153, 174704 (2020a).
  * Schlaich  _et al._ [2022] A. Schlaich, D. Jin, L. Bocquet, and B. Coasne, “Electronic screening using a virtual thomas–fermi fluid for predicting wetting and phase transitions of ionic liquids at metal surfaces,” Nat. Mater. 21, 237–245 (2022).
  * Scalfi and Rotenberg [2021] L. Scalfi and B. Rotenberg, “Microscopic origin of the effect of substrate metallicity on interfacial free energies,” Proc. Natl. Acad. Sci. U.S.A. 118, e210879118 (2021).
  * Wang  _et al._ [2017] J. Wang, J. Wang, Z. Kong, K. Lv, C. Teng, and Y. Zhu, “Conducting-polymer-based materials for electrochemical energy conversion and storage,” Adv. Mater. 29, 1703044 (2017).
  * Zhu  _et al._ [2016] J. Zhu, A. S. Childress, M. Karakaya, S. Dandeliya, A. Srivastava, Y. Lin, A. M. Rao, and R. Podila, “Defect-engineered graphene for high-energy- and high- power-density supercapacitor devices,” Adv. Mater. 28, 7185–7192 (2016).
  * Coretti  _et al._ [2022a] A. Coretti, C. Bacon, R. Berthin, A. Serva, L. Scalfi, I. Chubak, K. Goloviznina, M. Haefele, A. Marin-Laflèche, B. Rotenberg, S. Bonella, and M. Salanne, “MetalWalls: simulating electrochemical interfaces between polarizable electrolytes and metallic electrodes,” J. Chem. Phys. 157, 184801 (2022a).
  * Gross  _et al._ [2009] L. Gross, F. Mohn, N. Moll, P. Liljeroth, and G. Meyer, “The chemical structure of a molecule resolved by atomic force microscopy,” Science 325, 1110–1114 (2009).
  * Binninger [2021] T. Binninger, “Piecewise nonlinearity and capacitance in the joint density functional theory of extended interfaces,” Phys. Rev. B 103, L161403 (2021).
  * Comtet  _et al._ [2017] J. Comtet, A. Niguès, V. Kaiser, B. Coasne, L. Bocquet, and A. Siria, “Nanoscale capillary freezing of ionic liquids confined between metallic interfaces and the role of electronic screening,” Nat. Mater. 16, 634–639 (2017).
  * Scalfi  _et al._ [2020b] L. Scalfi, D. T. Limmer, A. Coretti, S. Bonella, P. A. Madden, M. Salanne, and B. Rotenberg, “Charge fluctuations from molecular simulations in the constant-potential ensemble,” Phys. Chem. Chem. Phys. 22, 10480–10489 (2020b).
  * Gingrich and Wilson [2010] T. R. Gingrich and M. Wilson, “On the ewald summation of gaussian charges for the simulation of metallic surfaces,” Chem. Phys. Lett. 500, 178–183 (2010).
  * Coretti  _et al._ [2022b] A. Coretti, C. Bacon, R. Berthin, A. Serva, L. Scalfi, I. Chubak, K. Goloviznina, M. Haefele, A. Marin-Laflèche, B. Rotenberg, S. Bonella, and M. Salanne, “Metalwalls: Simulating electrochemical interfaces between polarizable electrolytes and metallic electrodes,” J. Chem. Phys. 157, 184801 (2022b).
  * Pádua [2017] A. A. H. Pádua, “Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids,” J. Chem. Phys. 146, 204501 (2017).
  * Goloviznina  _et al._ [2019] K. Goloviznina, J. N. Canongia Lopes, M. Costa Gomes, and A. A. H. Pádua, “Transferable, Polarizable Force Field for Ionic Liquids,” J. Chem. Theory Comput. 15, 5858–5871 (2019).
  * Oll  _et al._ [2017] O. Oll, T. Romann, C. Siimenson, and E. Lust, “Influence of chemical composition of electrode material on the differential capacitance characteristics of the ionic liquid$|$electrode interface,” Electrochem. commun 82, 39–42 (2017).
  * Iamprasertkun  _et al._ [2019] P. Iamprasertkun, W. Hirunpinyopas, A. Keerthi, B. Wang, B. Radha, M. A. Bissett, and R. A. W. Dryfe, “Capacitance of basal plane and edge-oriented highly ordered pyrolytic graphite: Specific ion effects,” J. Phys. Chem. Lett. 10, 617–623 (2019).
  * Haskins and Lawson [2016] J. B. Haskins and J. W. Lawson, “Evaluation of molecular dynamics simulation methods for ionic liquid electric double layers,” J. Chem. Phys. 144, 184707 (2016).
  * Ferreira, Sampaio, and Siqueira [2022] P. H. L. Ferreira, A. M. Sampaio, and L. J. A. Siqueira, “Energy and power performances of binary mixtures of ionic liquids in planar and porous electrodes by molecular dynamics simulations,” Electrochimica Acta 410, 139982 (2022).
  * He  _et al._ [2024] H. He, J. Zhou, L. Yang, C. Liang, S. Xu, M. Chen, and T. Mo, “Accurately simulating electrical double layers structure and formation using all-atom scaled-charge force fields,” J. Mater. Chem. A , – (2024).
  * Hu  _et al._ [2013] Z. Hu, J. Vatamanu, O. Borodin, and D. Bedrov, “A molecular dynamics simulation study of the electric double layer and capacitance of $[$bmim$]$$[$pf6$]$ and $[$bmim$]$$[$bf4$]$ room temperature ionic liquids near charged surfaces,” Phys. Chem. Chem. Phys. 15, 14234–14247 (2013).
  * Carstens  _et al._ [2014] T. Carstens, R. Gustus, O. Höfft, N. Borisenko, F. Endres, H. Li, R. J. Wood, A. J. Page, and R. Atkin, “Combined stm, afm, and dft study of the highly ordered pyrolytic graphite/1-octyl-3-methyl-imidazolium bis(trifluoromethylsulfonyl)imide interface,” J. Phys. Chem. C 118, 10833–10843 (2014).
  * Hayes, Warr, and Atkin [2010] E. Hayes, G. G. Warr, and R. Atkin, “At the interface: solvation and designing ionic liquids,” Phys. Chem. Chem. Phys. 12, 1709–1723 (2010).
  * Black  _et al._ [2013] J. M. Black, D. Walters, A. Labuda, G. Feng, P. C. Hillesheim, S. Dai, P. T. Cummings, S. V. Kalinin, R. Proksch, and N. Balke, “Bias-dependent molecular-level structure of electrical double layer in ionic liquid on graphite,” Nano Lett. 13, 5954–5960 (2013).
  * Serva  _et al._ [2021] A. Serva, L. Scalfi, B. Rotenberg, and M. Salanne, “Effect of the metallicity on the capacitance of gold–aqueous sodium chloride interfaces,” J. Chem. Phys. 155, 044703 (2021).
  * Limaye, Suvlu, and Willard [2024] A. Limaye, D. Suvlu, and A. P. Willard, “Water molecules mute the dependence of the double-layer potential profile on ionic strength,” Faraday Discuss. 249, 267–288 (2024).
  * Kondrat and Kornyshev [2011] S. Kondrat and A. A. Kornyshev, “Superionic State in Double-Layer Capacitors with Nanoporous Electrodes,” J. Phys.: Condens. Matter 23, 022201 (2011).
  * Kondrat  _et al._ [2023] S. Kondrat, G. Feng, F. Bresme, M. Urbakh, and A. A. Kornyshev, “Theory and simulations of ionic liquids in nanoconfinement,” Chem. Rev. 123, 6668–6715 (2023).
  * Futamura  _et al._ [2017] R. Futamura, T. Iiyama, Y. Takasaki, Y. Gogotsi, M. J. Biggs, M. Salanne, J. Ségalini, P. Simon, and K. Kaneko, “Partial breaking of the coulombic ordering of ionic liquids confined in carbon nanopores,” Nat. Mater. 16, 1225–1232 (2017).
  * Simoncelli  _et al._ [2018] M. Simoncelli, N. Ganfoud, A. Sene, M. Haefele, B. Daffos, P.-L. Taberna, M. Salanne, P. Simon, and B. Rotenberg, “Blue energy and desalination with nanoporous carbon electrodes: Capacitance from molecular simulations to continuous models,” Phys. Rev. X 8, 021024 (2018).
  * Dufils  _et al._ [2023] T. Dufils, L. Knijff, Y. Shao, and C. Zhang, “PiNNwall: Heterogeneous electrode models from integrating machine learning and atomistic simulation,” J. Chem. Theory Comput. 19, 5199–5209 (2023).