#  Quantum-aided Multi-Objective Routing Optimization Using Back-Tracing-Aided
Dynamic Programming

Dimitrios Alanis, Panagiotis Botsinis, Zunaira Babar, Hung Viet Nguyen, Daryus
Chandra, Soon Xin Ng, and Lajos Hanzo The authors are with the School of
Electronics and Computer Science, University of Southampton, Southampton, SO17
1BJ, UK (email: {da1d16, pb1y14, zb2g10, hvn08r, dc2n14, sxn,
lh}@ecs.soton.ac.uk).The financial support of the EPSRC under the grant
EP/L018659/1, that of the European Research Council, Advanced Fellow Grant and
that of the Royal Society’s Wolfson Research Merit Award is gratefully
acknowledged. Additionally, the authors acknowledge the use of the IRIDIS High
Performance Computing Facility, and associated support services at the
University of Southampton, in the completion of this work. he research data of
this paper can be found at https://doi.org/10.5258/SOTON/D0479.

###### Abstract

Pareto optimality is capable of striking the optimal trade-off amongst the
diverse conflicting QoS requirements of routing in wireless multihop networks.
However, this comes at the cost of increased complexity owing to searching
through the extended multi-objective search-space. We will demonstrate that
the powerful quantum-assisted dynamic programming optimization framework is
capable of circumventing this problem. In this context, the so-called
Evolutionary Quantum Pareto Optimization (EQPO) algorithm has been proposed,
which is capable of identifying most of the optimal routes at a near-
polynomial complexity versus the number of nodes. As a benefit, we improve
both the the EQPO algorithm by introducing a back-tracing process. We also
demonstrate that the improved algorithm, namely the Back-Tracing-Aided EQPO
(BTA-EQPO) algorithm, imposes a negligible complexity overhead, while
substantially improving our performance metrics, namely the relative frequency
of finding all Pareto-optimal solutions and the probability that the Pareto-
optimal solutions are indeed part of the optimal Pareto front.

###### Index Terms:

Quantum Computing, QoS, Dynamic Programming, Pareto Optimality, Routing,
Multi-objective Optimization.

##  I Introduction

Routing optimization in _Wireless Multihop Networks_ (WMHN) has to strike a
trade-off among diverse and often conflicting _Quality-of-Service_ (QoS)
requirements[1]. For this reason several metrics have been advocated, such as
the _Network Lifetime_ (NL) [2] or the _Network Utility_ (NU) [3], which are
_single-objective_ aggregate functions of multiple QoS requirements. However,
these single-objective metrics may not be giving justice to all design
objectives. This problem can be circumvented by employing the concept of
Pareto optimality [4, 5]. This comes at the cost of increased complexity
imposed by the extended search-space, which can be in turn circumvented by
utilizing the powerful optimization framework of quantum computing [6].

In this context, several contributions on quantum-aided multi-objective
routing exist in the literature [7, 8, 9, 10]. To elaborate further, the so-
called _Non-Dominated Quantum Optimization_ (NDQO) and the _Non-Dominated
Quantum Iterative Optimization_ (NDQO) algorithms have been proposed in [7]
and [8], respectively, relying on full-search-based database exploration. As
an intermediate step, the so-called _Non-Dominated Quantum Optimization_
(MODQO) algorithm of [9] exploited the database correlations emerging from the
formation of Pareto-optimal route-combinations for efficiently reducing the
database size, thus achieving a further complexity reduction. The database
correlation has been exploited in [10], where the _Evolutionary Quantum Pareto
Optimization_ (EQPO) algorithm has been introduced. More explicitly, the EQPO
algorithm, which is a feed-forward-style algorithm, achieved a further
complexity reduction by exploiting the potential correlations among the
individual links constituting Pareto-optimal routes. Nevertheless, this
complexity reduction comes at the price of reduced heuristic accuracy. Against
this background our contributions are summarized as follows:

  1. _1)_

We propose an improved version of the EQPO, namely the Back-Tracing-Aided EQPO
(BTA-EQPO) algorithm, by introducing novel Back-Tracing Processes (BTPs) by
extending the quantum-aided dynamic programming framework of [10].

  2. _2)_

We demonstrate that the BPTs impose an insignificant complexity overhead, when
compared to the complexity imposed by the EQPO algorithm, hence the BTA-EQPO
imposes the same order of complexity as its predecessor, namely the EQPO.

  3. 3)

We also demonstrate that the BTA-EQPO algorithm’s resultant error floor is an
order of magnitude below that of its predecessor.

The rest of this paper is organized as follows. In Section II, we will present
the network topology considered. In Section III, we will elaborate on the
novel back-tracing process of the BTE-EQPO algorithm. We will then evaluate
its performance versus complexity in Section IV.

##  II Network Specifications

We have adopted the WMHN model considered in [7, 8, 10], where the _Source
Node_ (SN) and the _Destination Node_ (DN) are located at the opposite corners
of a $(100\times 100)$ m2 square block. By contrast, the _Relay Nodes_ (RNs)
are mobile, having locations that are uniformly distributed within this square
block. We also assume that the DN acts as a cluster-head, which has access to
a universal quantum computer. Each node experiences random interference power,
relying on a normal distribution with its mean set to -90 dBm and its standard
deviation to 10 dB. An example of the network topology consisting of
$N_{\text{nodes}}=5$ nodes is shown in Fig. 1.

Fig. 1: Exemplified WMHN topology associated with $N_{\text{nodes}}=5$ nodes
[10].

As for our optimization metrics, we have jointly considered the routes’ end-
to-end delay $D$, their total _Bit Error Ratio_ (BER) $P_{e}$ as well as their
total power dissipation $L$ in a similar fashion to [7, 8, 10]. More
specifically, we have considered _Quadrature Phase Shift Keying_ (QPSK)
transmissions in an uncorrelated Rayleigh fading environment, where the packet
forwarding has been carried out using the _Decode-and-Forward_ (DF) scheme
[11]. Consequently, the route’s overall BER $P_{e}(x)$ can be calculated using
the following recursive formula [7]:

$P_{e,tot}=P_{e,1}+P_{e,2}-2P_{e,1}P_{e,2},$ (1)

where $P_{e,tot}$ corresponds to the output BER of a two-stage _Binary
Symmetric Channel_ (BSC) [7] with $P_{e,1}$ and $P_{e,2}$ representing the
individual BER of the first and the second stage, respectively. Additionally,
the route’s end-to-end delay $D$ is quantified in terms of the number of hops
composing the route, while the total power dissipation $L$ is determined by
the sum of the path-losses of each individual link $L_{ij}$ of the route.
Explicitly, each link between the $i$-th and the $j$-th nodes exhibits path-
losses quantified in dB as follows [8]:

$L^{dB}_{ij}=10\alpha\log_{10}\left(\frac{4\pi d_{ij}}{\lambda_{c}}\right),$
(2)

where $\alpha$ corresponds to the path-loss exponent, which is set to
$\alpha=3$, $d_{ij}$ denotes the Euclidean distance between the $i$-th and the
$j$-th, while $\lambda_{c}$ is the carrier’s wavelength set to
$\lambda_{c}=0.125$ m. Therefore, the _Utility Vector_ (UV) $\mathbf{f}(x)$ of
the $x$-th route can be expressed as follows:

$\mathbf{f}(x)=\left[P_{e}(x),L(x),D(x)\right].$ (3)

The concept of Pareto optimality [5] has been adopted for evaluating the
fitness of the UVs. In a nutshell, a specific route $x_{1}$ dominates another
route $x_{2}$, i.e. we have $\mathbf{f}(x_{1})\succ\mathbf{f}(x_{2})$, if all
the individual metrics of $\mathbf{f}(x_{1})$ are lower than the respective
components of $\mathbf{f}(x_{2})$. Based on this principle, a route is
considered to be Pareto optimal, if there are no other routes dominating it.
Note that our ultimate goal is to identify the entire set of Pareto optimal
routes, which jointly constitute the so-called _Optimal Pareto Front_ (OPF)
[7].

##  III Back-Tracing-Aided Quantum Pareto Optimization

The BTA-EQPO algorithm, which is presented in Alg. 1, is constituted by three
distinct parts: a stage of the single-objective optimization followed by the
so-called _Single-Objective Back-Tracing Process_ (SO-BTP), a stage of the
multi-objective optimization process in a similar fashion to the EQPO
algorithm and a stage invoking a _Multi-Objective Back-Tracing Process_ (MO-
BTP).

Algorithm 1 Back-Tracing-Aided Evolutionary Quantum Pareto Optimization (BTA-
EQPO) Algorithm.

1: Set $S^{\text{OPF}}_{(i)}\leftarrow\varnothing$ and
$S^{\text{surv}}_{(i)}\leftarrow\varnothing$ $\forall
i\in\\{0,..,N_{\text{nodes}}-1\\}$.

2: Determine the optimal routes $S^{\text{opt}}$ based on each individual
objective based on the optimal framework presented in [10, Sec. III] and store
accordingly the optimal routes visited in $S^{\text{OPF}}_{(i)}$, where $i$ is
the number of RNs constituting the visited route.

3: For each route in $S^{\text{opt}}$ perform SO-BTP based on Fig. 2 and store
accordingly the surviving routes visited to the set $S^{\text{surv}}_{(i)}$,
where $i$ is the number of RNs constituting the visited route.

4: Set $S^{\text{gen}}_{(0)}\leftarrow\\{SN\rightarrow DN\\}$,
$S^{\text{OPF}}_{(0)}\leftarrow S^{\text{gen}}_{(0)}$,
$S^{\text{surv}}_{(0)}\leftarrow S^{\text{gen}}_{(0)}$, $i\leftarrow 0$.

5: repeat

6: Set $i\leftarrow i+1$.

7: Generate the set of routes $S^{\text{gen}}_{(i)}$ from the set
$S^{\text{surv}}_{(i-1)}$ by appropriately inserting a single RN between two
intermediate nodes.

8: Set $S^{\text{gen}}_{(i)}\leftarrow S^{\text{gen}}_{(i)}\cup
S^{\text{OPF}}_{(i-1)}$.

9: Invoke the P-NDQIO algorithm of [10, Alg. 2] in the set
$S^{\text{gen}}_{(i)}$ and initialize the identified OPF to
$S^{\text{OPF}}_{(i)}\leftarrow S^{\text{OPF}}_{(i-1)}\cup
S^{\text{OPF}}_{(i)}$.

10: Set
$S^{\text{surv}}_{(i)}\leftarrow\left(S^{\text{OPF}}_{(i)}-S^{\text{OPF}}_{(i-1)}\right)\cup
S^{\text{surv}}_{(i)}$.

11: until $\left|S^{\text{surv}}_{(i)}\right|=0$ or $i=N_{\text{nodes}}-2$

12: Set $i\leftarrow i+1$.

13: For each route in $S^{\text{OPF}}_{(i)}$ perform MO-BTP for $n$ trellis-
stages based on Fig. 2 and store the surviving routes visited in
$S^{\text{gen}}_{(i)}$.

14: Invoke the P-NDQIO algorithm of [10, Alg. 2] in the set
$S^{\text{gen}}_{(i)}\leftarrow S^{\text{gen}}_{(i)}\cup
S^{\text{OPF}}_{(i-1)}$ and initialize the identified OPF to
$S^{\text{OPF}}_{(i)}\leftarrow S^{\text{OPF}}_{(i-1)}$.

15: Export the OPF $S^{\text{OPF}}_{(i)}$ and terminate.

As far as the first stage is concerned, we first invoke in Step 2 of Alg. 1
single-objective dynamic programming based optimization utilizing the optimal
dynamic programming framework of [10, Sec. III] for the sake of identifying
the optimal routes $S^{\text{opt}}$ in terms of each individual objective.
These routes will also be Pareto-optimal [5], when jointly optimizing the UV
of Eq. (3). Therefore, we will appropriately initialize the set
$\\{S^{\text{OPF}}_{(i)}\\}_{i=0}^{N_{\text{nodes}}-2}$ of Pareto-optimal
routes to the set $S^{\text{opt}}$ based on the trellis-stage index $i$,
during which they were identified. For instance, the optimal route
$1\rightarrow 2\rightarrow 3\rightarrow 4\rightarrow 5$ will be appended to
$S^{\text{OPF}}_{(3)}$, since it consists of 3 RNs and thus it was identified
at the second trellis-stage. We have opted for this optimal framework, since
it guarantees the detection of these globally optimal routes, while it imposes
a complexity111We quantify the complexity in terms of the number of dominance
comparisons; a single dominance comparison is defined as a single _Cost
Function Evaluation_ (CFE). We further distinguish the complexity into two
domains: the _parallel complexity_[10], which takes into account the
beneficial hardware parallelism exploited by the NDQIO-based algorithms, and
the _sequential complexity_[10], which neglects the benefits of hardware
parallelism and it is simply quantified in terms of the number of Pareto-
dominance comparisons. In our application we have utilized quantum Pareto-
dominance comparison operators that are identical to those of [8].
Consequently, assuming a total of $a$ reference routes and $k$ optimization
objectives, a single activation of this quantum dominance operator results in
a parallel and a sequential complexity of $1/k$ and $a$ _Cost Function
Evaluations_ (CFEs), respectively. on the order of $O(N_{\text{nodes}}^{3})$.
Explicitly, there exist precisely $N_{\text{nodes}}$ surviving routes at each
trellis-stage, thus a total of $N_{\text{nodes}}^{2}$ comparisons are required
per trellis-stage, while a total of $O(N_{\text{nodes}})$ trellis stages are
processed.

Subsequently, the SO-BTP is activated in Step 3 of Alg. 1 for each of the
globally optimal routes identified by the optimization process in Step 2 of
Alg. 1. During this process, starting from a single optimal route we
successively trace back to the direct route by removing the last RN of the
route, as portrayed in the upper sub-figure of Fig. 2. We conceived utilized
this specific strategy, since the routes of a specific trellis stage are
generated by appropriately inserting an RN between the last RN and the DN at
each of the surviving routes of the previous trellis-stage. Additionally, the
surviving routes w.r.t. an individual objective will be also classified as
surviving [5], when we jointly optimize the entire set of objectives, since
their sub-routes will remain non-dominated by any other route or sub-route.
Using this observation, we will appropriately initialize the set
$\\{S^{\text{surv}}_{(i)}\\}_{i=0}^{N_{\text{nodes}}-2}$ of surviving routes
to the specific routes visited during each of the SO-BTPs in a similar fashion
to the initialization of
$\\{S^{\text{OPF}}_{(i)}\\}_{i=0}^{N_{\text{nodes}}-2}$.

Fig. 2: Single- versus multi-objective back-tracing.

After the initialization of the surviving routes, a multi-objective
optimization process similar to that of the EQPO algorithm of [10, Alg. 1] is
activated in Steps 5-11 of Alg 1. Their main difference is that both the set
of surviving and Pareto optimal routes have been initialized by the SO-BTP, as
highlighted in Steps 9 and 10 of Alg. 1. Naturally, the initialization of the
surviving routes expands the search-space, hence rendering the BTA-EQPO
capable of identifying a more diverse set of Pareto optimal routes. This
search-overhead imposed by the additionally generated routes is on the order
of $O(N_{\text{nodes}})$ extra cost-function evaluations, when using a similar
approach to that of [10]. Since the number of generated routes excluding this
overhead at the $i$-th trellis-stage is on the order of
$O(N_{\text{OPF}}N_{\text{nodes}}i)$ with $N_{\text{OPF}}$ representing the
number of Pareto optimal routes, we may deem this overhead to be low.
Quantitatively, the second step imposes the same order of complexity as the
EQPO algorithm, whose parallel and sequential complexity were shown to be on
the order of $O(N_{\text{OPF}}^{3/2}N_{\text{nodes}}^{2})$ and
$O(N_{\text{OPF}}^{5/2}N_{\text{nodes}}^{2})$, respectively. Naturally, the
complexity order of the fist stage can be considered as negligible compared to
that of the second stage.

Finally, the third stage in Steps 13-14 of Alg. 1 is activated, which invokes
the MO-BTP for $n$ trellis stages and it is invoked for each of the hitherto
identified OPF routes. To further aid its exposition, its employment is
visually portrayed in the bottom sub-figure of Fig. 2. During the MO-BTP, the
inverse of Step 7 of Alg. 1 is carried out, i.e. we move to the previous
trellis stage by removing a single RN from the route examined. For instance,
observe in Fig. 2 that invoking the MO-BTP for the Pareto optimal route
$1\rightarrow 2\rightarrow 3\rightarrow 4\rightarrow 5$ results in visiting
the routes $1\rightarrow 2\rightarrow 3\rightarrow 5$, $1\rightarrow
2\rightarrow 4\rightarrow 5$ and $1\rightarrow 3\rightarrow 4\rightarrow 5$,
when back-tracing for $n=1$ trellis stage, and the routes $1\rightarrow
2\rightarrow 5$, $1\rightarrow 3\rightarrow 5$ as well as $1\rightarrow
4\rightarrow 5$, when back-tracing for $n=2$ trellis stages. During this
process, we keep track of the visited routes of the MO-BTP, storing them while
we reach the final set $S^{\text{gen}}_{(i)}$ of generated routes. We then
invoke the _Preinitialized NDQIO_ (P-NDQIO) algorithm [10, Alg. 2] with its
OPF initialized to the hitherto identified OPF emanating from the second stage
for the sake of finding any further Pareto optimal routes. The complexity
order of the P-NDQIO algorithm is proportional to $O(\sqrt{N})$ [10]. We have
chosen to optimize the routes over the entire database, since offers a
beneficial complexity reduction against performing the optimization for each
backward trellis transition, since we have
$\sqrt{\sum_{i}n_{i}}<\sum_{i}\sqrt{n_{i}}$.

Last but not least, let us quantify the extra complexity imposed by this
process. The total number of generated routes as a function of the number $n$
of backward trellis transitions can be readily shown to be on the order of
$O(N_{\text{OPF}}N_{\text{nodes}}^{n})$. Consequently, the parallel and
sequential complexities imposed by the P-NDQIO algorithms of the MO-BTP may be
shown to be on the orders of $O(N^{3/2}_{\text{OPF}}N_{\text{nodes}}^{n/2})$
and $O(N^{5/2}_{\text{OPF}}N_{\text{nodes}}^{n/2})$, respectively. Hence, the
total complexity imposed by the BTE-EQPO can be shown to be:

$\displaystyle L^{P}_{\text{BTA-EQPO}}$
$\displaystyle=O\left[N^{3/2}_{\text{OPF}}\left(N_{\text{nodes}}^{2}+N_{\text{nodes}}^{n/2}\right)\right],$
(4) $\displaystyle L^{S}_{\text{BTA-EQPO}}$
$\displaystyle=O\left[N^{5/2}_{\text{OPF}}\left(N_{\text{nodes}}^{2}+N_{\text{nodes}}^{n/2}\right)\right].$
(5)

Hence, the MO-BTP will dominate the complexity orders, when having more than
$n=4$ backward-trellis steps. Let us now proceed by examining the performance
versus complexity trade-off of the BTA-EQPO algorithm.

(a)  (b)

Fig. 3: Parallel (a) and sequential (b) complexities of the BTE-EQPO algorithm
compared to those of the EQPO, NDQIO and NDQIO algorithms. The results have
been averaged over 108 runs.

##  IV Performance versus Complexity

In this section we will provide some further insights concerning BTA-EQPO
algorithm’s performance versus complexity and compare it to the existing
quantum assisted algorithms, namely the EQPO [10], NDQIO [8] and NDQO [7]
algorithms. We will first examine the average complexity imposed by the
aforementioned algorithms as a function of the number $N_{\text{nodes}}$ of
nodes constituting the WMHN. In addition to the aforementioned algorithms we
investigate a hybrid algorithm, which uses the first two stages of the BTA-
EQPO, while the third stage is replaced by a full-database search carried out
by the P-NDQIO algorithm [10]. The latter will be referred to as “BTA-EQPO
with P-NDQIO” and it is used as the upper bound of the complexity imposed by
MO-BTP, when we have $n=N_{\text{nodes}}-1$.

The average parallel and sequential complexities are shown in Figs.
LABEL:fig:complexity-p and LABEL:fig:complexity-s. In these figures we vary
the number $n$ of backward trellis stages in the range of $\\{0,1,2,4\\}$.
Note that for $n=0$ only the SO-BTP is active, while for $n=4$ the MO-BTP
complexity orders match those of the BTA-EQPO algorithm’s second stage.
Observe in both figures that both the parallel and the sequential complexity
imposed by the BTA-EQPO algorithm approach that of the EQPO algorithm, hence
verifying our theoretical analysis of Sec. III, where we proved that the extra
complexity imposed both by the SO-BTP and by the MO-BTP is significantly lower
than the complexity of BTA-EQPO algorithm’s second stage. Furhtermore, observe
in Fig. LABEL:fig:complexity-p that the BTA-EQPO algorithm imposes almost the
same parallel complexity as BTA-EQPO with P-NDQIO algorithm. However, a a
factor of two sequential complexity increase is observed in Fig.
LABEL:fig:complexity-s for 9-node WMHNs. This is because the square root of
the total number of routes is close to that of the routes created by the MO-
BTP for the WMHN sizes we investigated; however, for larger WMHNs we expect
much higher complexity reduction for our BTA-EQPO algorithm. Additionally,
both a parallel and a sequential complexity reduction is achieved against the
NDQIO algorithm, which almost is a high as an order of magnitude for 9-node
WMHNs.

(a)  (b)  (c)  (d)

Fig. 4: Performance versus complexity in terms of the average Pareto distance
$E[P_{d}]$ (a,c) and the average Pareto completion $E[C]$ for 8-node WMHNs. In
(b) and (d) the value $1-E[C]$ is portrayed facilitating observation of the
differences among the algorithms. The results have been averaged over 108
runs.

Continuing with the BTA-EQPO algorithm’s performance evaluation, we will
utilize two metrics: the average _Pareto distance_ $E[P_{d}]$ [7], which is
defined as the probability of a route identified as Pareto optimal being truly
Pareto optimal, and the average _Pareto completion_ $E[C]$ [7], defined as the
average fraction of the true OPF being identified by a heuristic method.
Naturally, for $E[P_{d}]=0$ the identified OPF exclusively consists of true
Pareto optimal routes, while for $E[C]=1$ the entire true OPF has been
identified. The average Pareto distance $E[P_{d}]$ is shown in Figs.
LABEL:fig:pd-p and LABEL:fig:pd-s as a function of the parallel and sequential
complexity invested, respectively. Observe in these figures that the BTA-EQPO
algorithm associated with $n=0$, i.e. with the particular case where SO-BTP is
active, has a similar performance to that of the EQPO algorithm [10]. However,
the beneficial effects of MO-BTP are visible even for $n=1$, where $E[P_{d}]$
is reduced by a factor of 5 after 2,200 and 28,000 CFEs in the parallel and
sequential complexity domains, respectively, when compared to the EQPO
algorithm, where the latter is portrayed with the aid of the gray solid lines.
This improvement is further enhanced for $n=2$ and $n=4$, where $E[P_{d}]$ is
improved by an order of magnitude compared to that of the EQPO algorithm.
Additionally, observe in Figs. LABEL:fig:pd-p and LABEL:fig:pd-s that beyond
$n=2$ the BTP-EQPO algorithm exhibits an error floor formation, hence
rendering the application of further backward-trellis steps redundant. As for
the full-search-based methods, observe in Fig. LABEL:fig:pd-p that the BTA-
EQPO with P-NDQIO algorithm becomes more efficient than both the NDQO and the
NDQIO algorithms beyond a parallel complexity of 3,000 CFEs, while its
$E[P_{d}]$ decays to infinitesimally low levels beyond 3,500 CFEs. This trend
is also present in Fig. LABEL:fig:pd-s; however, observe that the NDQIO
algorithm is more efficient than the BTA-EQPO with P-NDQIO algorithm. However,
we expect this trend to change following that of Fig. LABEL:fig:pd-p as the
number of nodes increases, where the BTA-EQPO with P-NDQIO algorithm offers a
substantial sequential complexity reduction compared to the NDQIO algorithm.

As far as the average Pareto completion is concerned, observe in Figs.
LABEL:fig:c-p and LABEL:fig:c-s that the BTA-EQPO algorithm associated with
$n=0$ succeeds in identifying a larger fraction of the OPF by improving the
complementary Pareto completion metric by a factor of 3. This happens at a
parallel and a sequential complexity of 3,500 and 49,000 CFEs, respectively,
thus explicitly demonstrating the benefit of the SO-BTP. When the MO-BTP is
activated, this metric is further reduced, exhibiting of an order of magnitude
total improvement over EQPO algorithm. Additionally, we can observe that this
metric is slightly improved, as the number $n$ of backward-trellis steps
increases. Explicitly, the Pareto Completion error floor exhibited stems from
the BTA-EQPO and EQPO algorithms’ property of terminating the trellis stages,
when no Pareto optimal routes are detected. Thus, they are incapable of even
examining potential Pareto-optimal routes located at later trellis stages.
This limitation is partially mitigated by the SO-BTP, which rectifies the
deficiency, where a globally optimal route may be located several stages apart
from the rest of the OPF. Despite this inability, the BTA-EQPO algorithm’s
performance is near-optimal, identifying the Pareto optimal routes with 0.1%
probability of misdetection, while being able to detect 99.97% of the time the
true OPF.

##  V Conclusions

We have further developed the quantum-assisted multi-objective dynamic
programming framework of [10] by introducing the SO-BTP and the MO-BTP for the
sake of enhancing the heuristic accuracy attained. We have shown that the SO-
BTP enables the algorithm to detect almost all of the Pareto optimal
solutions, while the activation of MO-BTP also increases our confidence in
detecting only the true Pareto-optimal routes. Finally, we have proven that
the SO-BTP’s extra complexity is insignificant. Furthermore, we have
demonstrated for the MO-BTP that its extra complexity is insignificant as long
as we employ less than 5 backward-trellis steps. Finally, we have demonstrated
that with the above proviso the BTA-EQPO algorithm outperforms the EQPO and
exhibits a near-optimal accuracy.

## References

  * [1] Z. Fei, B. Li, S. Yang, C. Xing, H. Chen, and L. Hanzo, “A survey of multi-objective optimization in wireless sensor networks: Metrics, algorithms and open problems,” IEEE Communications Surveys Tutorials, vol. 19, no. 1, pp. 550–586, 2017.
  * [2] F. Tashtarian, M. H. Y. Moghaddam, K. Sohraby, and S. Effati, “On maximizing the lifetime of wireless sensor networks in event-driven applications with mobile sinks,” IEEE Transactions on Vehicular Technology, vol. 64, pp. 3177–3189, July 2015.
  * [3] A. Zhou, M. Liu, Z. Li, and E. Dutkiewicz, “Joint traffic splitting, rate control, routing, and scheduling algorithm for maximizing network utility in wireless mesh networks,” IEEE Transactions on Vehicular Technology, vol. 65, pp. 2688–2702, April 2016.
  * [4] S. Liao and Q. Zhang, “A multiutility framework with application for studying tradeoff between utility and lifetime in wireless sensor networks,” IEEE Transactions on Vehicular Technology, vol. 64, pp. 4701–4711, Oct 2015\.
  * [5] K. Deb, “Multi-objective optimization,” in Search Methodologies (E. K. Burke and G. Kendall, eds.), pp. 273–316, Springer US, 2005.
  * [6] M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information.  Cambridge university press, 2010.
  * [7] D. Alanis, P. Botsinis, S. X. Ng, and L. Hanzo, “Quantum-assisted routing optimization for self-organizing networks,” IEEE Access, vol. 2, pp. 614–632, 2014.
  * [8] D. Alanis, P. Botsinis, Z. Babar, S. X. Ng, and L. Hanzo, “Non-dominated quantum iterative routing optimization for wireless multihop networks,” IEEE Access, vol. 3, pp. 1704–1728, 2015.
  * [9] D. Alanis, J. Hu, P. Botsinis, Z. Babar, S. X. Ng, and L. Hanzo, “Quantum-assisted joint multi-objective routing and load balancing for socially-aware networks,” IEEE Access, vol. 4, pp. 9993–10028, 2016.
  * [10] D. Alanis, P. Botsinis, Z. Babar, H. V. Nguyen, D. Chandra, S. X. Ng, and L. Hanzo, “A quantum-search-aided dynamic programming framework for pareto optimal routing in wireless multihop networks,” IEEE Transactions on Communications, vol. PP, no. 99, pp. 1–1, 2018.
  * [11] S. Yang, X. Xu, D. Alanis, S. X. Ng, and L. Hanzo, “Is the low-complexity mobile relay aided FFR-DAS capable of outperforming the high-complexity CoMP?,” IEEE Transactions on Vehicular Technology, vol. 65, pp. 2154–2169, April 2016.


# Dynamics and energetics of ion adsorption at the interface between a pure
ionic liquid and carbon electrodes

Nicolas Gaudy  CIRIMAT, Université de Toulouse, CNRS, Université Toulouse 3 -
Paul Sabatier, 118 Route de Narbonne, 31062 Toulouse cedex 9, France  Mathieu
Salanne  Sorbonne Université, CNRS, Physicochimie des Électrolytes et
Nanosystèmes Interfaciaux, F-75005 Paris, France  Céline Merlet  CIRIMAT,
Université de Toulouse, CNRS, Université Toulouse 3 - Paul Sabatier, 118 Route
de Narbonne, 31062 Toulouse cedex 9, France<EMAIL_ADDRESS>
###### Abstract

Molecular dynamics simulations have been used extensively to determine
equilibrium properties of the electrode-electrolyte interface in
supercapacitors held at various potentials. While such studies are essential
to understand and optimize the performance of such energy storage systems,
investigations of the dynamics of adsorption during the charge of the
supercapacitors is also necessary. Dynamical properties are especially
important to get an insight into the power density of supercapacitors, one of
their main assets. In this work, we propose a new method to analyze the
trajectories of adsorbing ions. We focus on pure 1-ethyl-3-methylimidazolium
bis(trifluoromethylsulfonyl)imide in contact with planar carbon electrodes. We
characterize the evolution of the ion orientation and ion-electrode distance
during adsorption and show that ions reorientate as they adsorb. We then
determine the forces experienced by the adsorbing ions and demonstrate that
Coulomb forces dominate at long range while van der Waals forces dominate at
short range. We also show that there is an almost equal contribution from the
two forces at an intermediate distance, explaining the peak of ion density
close to the electrode surface.

Réseau sur le Stockage Électrochimique de l’Énergie (RS2E), Fédération de
Recherche CNRS 3459, HUB de l’Énergie, Rue Baudelocque, 80039 Amiens, France
Réseau sur le Stockage Électrochimique de l’Énergie (RS2E), Fédération de
Recherche CNRS 3459, HUB de l’Énergie, Rue Baudelocque, 80039 Amiens, France
Réseau sur le Stockage Électrochimique de l’Énergie (RS2E), Fédération de
Recherche CNRS 3459, HUB de l’Énergie, Rue Baudelocque, 80039 Amiens, France

##  1 Introduction

Supercapacitors, also known as electrochemical capacitors, have emerged as
promising energy storage devices due to their exceptional power delivery, long
cycle life, and fast charge/discharge capabilities1. They store energy through
the formation of electrical double layers (EDLs) at the electrode-electrolyte
interface, leading to the compensation of electrode charges through reversible
adsorption of ions2. Similarly to the bulk electrolyte transport properties
(diffusion coefficients, electrical conductivity, viscosity), the dynamics of
ion adsorption is influenced by various factors, such as electrostatic
interactions, ion size and coordination. It has a direct effect on the power
characteristics of capacitive systems,3 hence it is essential to describe and
understand interfacial properties dynamically at the molecular level in order
to improve the performance of supercapacitors.

Classical molecular dynamics (MD) is a powerful computational tool for
investigating the behaviour of supercapacitors. By performing MD simulations,
it is possible to gain insights into the intricate processes occurring at the
molecular level within these energy storage devices. Due to their capability
to obtain trajectories for large systems ($>$ 10,000 atoms) and for relatively
extended duration (several nanoseconds)4, 5, 6, these simulations provide
valuable information on the mechanisms underlying charge storage, ion
migration and other crucial factors that influence the efficiency and
reliability of supercapacitors. In particular, MD simulations enable the
exploration of various aspects of supercapacitor operation, such as ion
diffusion, charge/discharge dynamics, and interfacial phenomena (EDL
structure, EDL dynamic, ion orientation with respect to the electrode etc). In
recent years, significant progress has been made in understanding the complex
structure of ionic layers in EDLs through MD simulation studies. The effects
of overscreening and ion packing on the structure of the EDL and capacitance
value have been investigated7, 8, 9, 10, 11. However, most of these studies
have focused on the equilibrium properties of supercapacitors such as the ion
density, the angular distribution of adsorbed ions, and more recently the
energy aspects 12, 13. Although some MD simulation studies have investigated
the charging dynamics of overall devices14, 15, 16, 3, a dynamic analysis of
the individual ions adsorption process has yet to be carried out.

To address this, in this work we develop a set of methods allowing to provide
a comprehensive and dynamic depiction of the ion adsorption process for a
system composed of a neat IL (1-ethyl-3-methylimidazolium
bis(trifluoromethylsulfonyl)imide, EMIM-TFSI) in contact with planar carbon
electrodes. In particular, we focus on examining the dynamic orientation of
adsorbing ions, as well as comparing the electrostatic and van der Waals (vdW)
forces exerted on these ions throughout the charging process of the
supercapacitor. By analysing the evolution of these quantities over time, we
show that ions adsorb at the same time as they reorient. Once there is a
vacancy in the EDL, adsorption is the result of the competition between the
electrostatic and van der Waals forces.

##  2 Methods

The MetalWalls software17, 18 was used to perform MD simulations of neat EMIM-
TFSI in contact with two graphite electrodes, between which a constant
potential difference was maintained. Figure 1 provides a representative
snapshot of the system.

Figure 1: Supercapacitor composed of two graphite electrodes (in grey) in
contact with a pure ionic liquid electrolyte EMIM-TFSI. A constant potential
difference is applied between the two electrodes.

The electrolyte comprises 322 ion pairs, and the electrodes consist of three
parallel graphene sheets that are maintained rigid and separated from each
other by 3.4 Å. The simulation box has dimensions of 34 $\times$ 37 $\times$
112 $\textrm{{\AA}}^{3}$, with two-dimensional periodic boundary conditions
applied in the $x$ and $y$ directions. No periodic boundary conditions are
applied in the $z$ direction, perpendicular to the electrodes.

In the constant potential method, each carbon atom $i$ of the electrode
carries a Gaussian charge distribution, with a width
$\eta_{i}^{-1}=0.40\,\textrm{{\AA}}$ in this work. Thus, the charge
distribution of an atom in the electrode will be of amplitude $Q_{i}$ and
width $\eta_{i}^{-1}$ , such that

$\rho_{i}\left(\mathbf{r}\right)=Q_{i}\eta_{i}^{3}\pi^{-3/2}{\rm
e}^{-\eta_{i}^{2}\left|\mathbf{r}-\mathbf{r_{i}}\right|^{2}}$ (1)

where $\mathbf{r}=(x,y,z)$ is the position considered in the simulation box
and $\mathbf{r_{i}}=(x_{i},y_{i},z_{i})$ is the position of the $i-th$
electrode atom. The electrode atom charges fluctuate in order to satisfy the
constant potential condition. 19

On the electrolyte side, the charges carried by the atoms are considered to be
point charges. Thus, the charge distribution of an atom $i$ of charge $Q_{i}$
is given by

$\rho_{i}\left(\mathbf{r}\right)=Q_{i}\updelta\left(\mathbf{r}-\mathbf{r_{i}}\right)$
(2)

To represent EMIM-TFSI, the all-atom CL&P force field20, 21 was employed. This
force field includes electrostatic and van der Waals interactions. Atomic
charges are scaled by a factor of $0.8$ to mimic polarizability effects, a
technique that has been suggested in previous studies of bulk liquids22, 23,
24, 25. For van der Waals interactions, a Lennard-Jones potential with a
cutoff of 14.8 Å is implemented, while the Ewald summation method is used to
calculate long-range interactions, with a real space cutoff of 17 Å. The
initial configuration of the system was generated using the PACKMOL
software26. The liquid density was adjusted by allowing the electrodes to move
in the $z$ direction, for a few nanoseconds, using a piston and applying a
pressure of 0 bar on both sides. Once the positions of the electrodes became
stable, a bulk density of $1.6\,\mathrm{g.cm^{-3}}$ was attained, which is
very close to the experimental value at ambient temperature27.

The system undergoes charging by applying potential differences of
$\Delta\Psi$ = 1 V, 2 V and 4 V between the two electrodes. The simulations
are conducted in the NVT ensemble with the temperature set to 298 K using a
Nosé-Hoover thermostat chain of length 5 with a relaxation time of 500 fs. A
constant applied potential is maintained between the two electrodes. 19
Initially, an equilibration phase is performed at 0 V. Following this
equilibration, a production run of 63 ns is carried out, using a timestep of 2
fs. Figure S3 shows the evolution of the electrode charge along the
simulation.

##  3 Results and discussion

In this work, we aim to investigate both structural changes and forces
individually experienced by the ions during adsorption. More specifically, we
focus on the adsorption of EMIM cations on the negative electrodes and TFSI
anions on the positive electrodes during electrode charging with 1 V, 2 V and
4 V potential differences applied. To emphasize specific phenomena occurring
during adsorption, we have developed a method for dynamically monitoring the
system and statistically studying adsorption.

###  3.1 Interfacial structure and identification of the adsorbed ions

A common way to characterize the liquid structure at the interface with an
electrode is to plot ionic densities. Figure 2 shows the density of EMIM,
$\rho$, as a function of the distance to the electrode, $d_{electrode}$, at
different times during the charge of the supercapacitor.

Figure 2: Ionic density of EMIM as a function of $d_{electrode}$ for
$\Delta\Psi=1\,\mathrm{V}$. The first peak increases when the total electrode
charge, $Q$, increases. The black vertical line represents the boundary
between the first and second adsorption layers, noted $z_{ads}$.

This distance is defined here as $z-z_{elec}$ where $z$ is the $z$-component
of the position considered in the box and $z_{elec}$ is the position of the
graphene sheet in contact with the electrolyte. For quantities where ionic
positions are considered, $d_{electrode}$ is defined as $z_{CM}-z_{elec}$
where $z_{CM}$ is the $z$-component of the position of the center of mass of
the considered ion.

There is a first intense adsorption peak for $d_{electrode}\sim$ 3.6 Å that
increases over time. Indeed, during the charge of the supercapacitor,
counterions adsorb on the electrodes and compensate their charge, which leads
to an increase in the ion concentration in the first adsorbed layer. According
to this ionic density profile, ions are considered to be adsorbed if $z_{CM}$
is located between the graphene sheet in contact with the electrolyte and the
first minimum of the ion density, noted $z_{ads}.$ For the remainder of this
article, analyzes are conducted only on ions that are not in the first
adsorption layer at $t=0$ and that are there once the supercapacitor is
charged, i.e. only ions which adsorb during charging.

####  3.1.1 Evolution of quantities during the adsorption process

In order to generalize the study of the adsorption process, a method for
statistically studying the dynamics of adsorption was developed. During the
charging of the supercapacitor, the ions that adsorb are identified and their
orientation and the forces acting on them over time are computed.
Subsequently, each ion is then assigned a time $t_{ads}$ such that
$d_{electrode}\left(t_{ads}\right)=z_{ads}$, the time after which the ion is
considered to be adsorbed. The orientations and the forces undergone by the
ions are then realigned with respect to $t_{ads}$ ($t^{\prime}=t-t_{ads})$.
This ensures that for all the ions, the new origin $t^{\prime}=0$ corresponds
to the instant when the ions passes from the non-adsorbed state to the
adsorbed state. In other words, for $t^{\prime}<0$ the ion is considered not
to be adsorbed and for $t^{\prime}>0$ the ion is considered to be adsorbed. By
aligning all quantities around the time of adsorption, averages of various
properties can be computed and standard deviations extracted. To correlate
local properties during adsorption, large vectors containing i) all the ion-
electrode distances and ii) the observed quantities during adsorption are
generated. It is then possible to plot any property as a function of the
electrode distance and conduct subsequent analyzes. More details about this
way of analyzing the trajectories done are available in SI.

###  3.2 Structural description

In order to characterize the orientation of the ions with respect to the
electrode on which they are adsorbed, a unit vector normal to the electrode,
labelled $\mathbf{n}$, and a vector associated with each ion, noted
$\mathbf{u}$, are defined. For EMIM, $\mathbf{u}$ is defined as a vector
normal to the plane of the imidazolium ring. For TFSI, $\mathbf{u}$ is defined
as the vector linking the two carbon atoms (Figure 3). The angle
$\theta=\widehat{\mathbf{n\cdot}\mathbf{u}}$ is used to determine the general
orientation of the ions. For EMIM ions, if $\theta$ is between 0∘ and 45∘, the
ion is considered parallel to the electrode while it is considered
perpendicular to the electrode if $\theta$ is between 45∘ and 90∘. For TFSI
ions, if $\theta$ is between 0∘ and 45∘, the ion is considered perpendicular
to the electrode while it is considered parallel to the electrode if $\theta$
is between 45∘ and 90∘. Although stricter definitions could be used, it would
not change the conclusions obtained.

Figure 3: Scheme explaining the definition of $\mathbf{u}$ vectors for each
ion

Figure 4 shows the average distance and orientation of EMIM ions as a function
of time, during adsorption, for $\Delta\Psi=1\,\mathrm{V}$. The insert shows
the variation of the angle as a function of the distance from the electrode.

Figure 4: In blue, average distance between the electrode and the centre of
mass of the EMIM ions that adsorb over time. In orange, average orientation of
the EMIM ions that adsorb over time. Shaded regions indicate the range of
statistical uncertainty. The applied potential is equal to
$\Delta\Psi=1\,\mathrm{V}$. The blue dashed line corresponds to $z_{ads}$, the
orange dashed line corresponds to $\theta=45{{}^{\circ}}$ and the vertical
black line corresponds to $t=0$ ns. Insert: average angle as function of the
distance.

It can be seen that before adsorption all the ions approach the electrode in a
perpendicular orientation. Indeed, for $d_{electrode}$ between 5Å and 6Å,
there is no variation in the angle $\theta$. When adsorption occurs, the ions
reorient themselves parallel to the electrode at the same time as they enter
the first adsorbed layer. Once the ions are adsorbed, the angle varies
linearly with the variation in distance from the electrode. As the ions
reorient, their centre of mass moves and so their distance from the electrode
decreases. Figure S4 reports the raw data for distances and angles, before
averaging, showing the extent of the deviations observed for $d_{electrode}$
and $\theta$. The same behaviour was observed for all potential differences
and for TFSI anions.

###  3.3 Forces description

Following the analysis of the interfacial structure of the ions, the forces
felt by the ions during the adsorption process are determined. Specifically,
we investigate electrostatic interactions, which are characterized by the
Coulomb force, as well as van der Waals interactions, which are derived from
the Lennard-Jones potential. The force derived from a potential exerted by a
particle $j$ on a particle $i$ is given by

$\mathbf{F_{ij}}=-\nabla\mathcal{U}\left(\mathbf{r_{ij}}\right)=\frac{\partial\mathcal{U}\left(\mathbf{r_{ij}}\right)}{\partial
r_{ij}}\times\frac{1}{r_{ij}}\left(\mathbf{r_{j}-\mathbf{r_{i}}}\right)$ (3)

where $\mathbf{r_{i}}$ and $\mathbf{r_{j}}$ are the positions of particles $i$
and $j$, $r_{ij}$ the distance between the two particles and $\mathcal{U}$ the
potential energy.

In this analysis, ions are considered as single molecular entities (and not by
their individual atoms) so the interactions between two ions are defined in
terms of the interactions involving their respective centers of mass. Then,
the position of a molecule is determined by the position of its centre of
mass, while its charge is calculated as the sum of all its partial atomic
charges. Since adsorption takes place along the $z$ direction and in order to
account for an ion “field of vision”, we define the neighborhood of an ion as
an infinite cone oriented along $z$ (Figure 5). Only the centres of mass and
carbon atoms in this neighborhood of the adsorbing ion are taken into account
when calculating the contributions of the forces. The dimensions of the cone
have been chosen so that at 9 Å (corresponding to the first minimum of the
EMIM-TFSI pair distribution function, see Figure S5) the diameter of the cone
is equal to 10 Å (approximate size of an electrolyte ion).

Figure 5: Scheme describing the definition of the neighborhood of an ion in
this work. The two colors represent anions and cations.

Using the expression for the Coulomb and Lennard-Jones potential energies
given in SI, the equations for the forces are derived. The Coulomb force
exerted by the particle $j$ of charge $Q_{j}$ on the particle $i$ of charge
$Q_{i}$ is given by

$\mathbf{F_{ij}^{Coulomb}}=-\frac{1}{4\pi\varepsilon_{0}}\frac{Q_{i}Q_{j}}{r_{ij}^{3}}\left(\mathbf{r_{j}-\mathbf{r_{i}}}\right)$
(4)

where $\varepsilon_{0}$ denotes the permittivity of vacuum. The Coulomb force
between an ion in the electrolyte and an atom in the electrode is given by

$\mathbf{F_{ij}^{Coulomb}}=-\frac{1}{4\pi\varepsilon_{0}}\frac{Q_{i}Q_{j}}{r_{ij}^{3}}\left(\mathrm{erf\left(\textrm{$\eta\,r_{ij}$}\right)}-\frac{2}{\sqrt{\pi}}\eta\,r_{ij}{\rm
e}^{-\left(\eta\,r_{ij}\right)^{2}}\right)\left(\mathbf{r_{j}-\mathbf{r_{i}}}\right)$
(5)

where $\mathrm{erf}$ is defined as the primitive of a Gaussian function. The
force derived from the Lennard-Jones potential, describing the van der Waals
interactions, is given by:

$\mathbf{F_{ij}^{vdW}}=24\varepsilon_{ij}\left(\frac{\sigma_{ij}^{12}}{r_{ij}^{14}}-2\frac{\sigma_{ij}^{6}}{r_{ij}^{8}}\right)\left(\mathbf{r_{j}-\mathbf{r_{i}}}\right)$
(6)

To work in the centre of mass representation and in order to compute
$\mathbf{F_{ij}^{vdW}}$, it is necessary to define a new set of parameters
$\varepsilon_{ii}$ and $\sigma_{ii}$ corresponding to anions and cations.
Following the method described in SI, the obtained values are
$\sigma_{\mathrm{EMIM}}=4.02\,\textrm{{\AA}}$,
$\sigma_{\mathrm{TFSI}}=5.45\,\textrm{{\AA}}$,
$\varepsilon_{\mathrm{EMIM}}=0.71~{}\mathrm{kJ~{}mol^{-1}}$ and
$\varepsilon_{\mathrm{TFSI}}=1.05~{}\mathrm{kJ~{}mol^{-1}}$.

The $z$ component of the total force exerted on the centre of mass of the ions
adsorbing during the charging of the supercapacitor is then given by

$\mathbf{F_{ij}^{tot}.e_{z}}\left(t^{\prime}\right)=\underset{j}{\sum}\left(\mathbf{F_{ij}^{Coulomb}}\left(t^{\prime}\right)+\mathbf{F_{ij}^{vdW}}(t^{\prime})\right)\mathbf{.e_{z}}$
(7)

where $j$ is the $j-th$ particle contained in the cone (Figure 5) and
$\mathbf{e_{z}}$ is a unit vector along $z$.

Figure 6: $z$-component of the mean total force felt by the adsorbing EMIM (in
magenta) and number of neighboring EMIM and TFSI ions (in red) as a function
of the distance from the electrode for $\Delta\Psi=$1V. The horizontal dashed
lines highlight i) where the total force is zero and ii) the frontier where
total number of ions in the cone is exactly 1. The solid black line is the
cation density profile (arbitrary units) as a guide for the eyes to know which
layers correspond to given distance $d_{electrode}$.

Figure 6 shows the mean total force felt by the EMIM during adsorption for
$\Delta\Psi=1\mathrm{\,V}$ as a function of the distance from the electrode.
For the negative electrode, the EMIM cations adsorb in the direction opposite
to the axis $z$. So, if the force is directed toward the electrode, then
$\mathbf{F_{ij}^{tot}.e_{z}}<0$ and then the force is attractive. Conversely,
if $\mathbf{F_{ij}^{tot}.e_{z}}>0$ the force is repulsive.

Different regimes can be identified in Figure 6. The force is i) attractive
and constant for $d_{electrode}$ between 6.1 and 12 Å, ii) increasing linearly
for $d_{electrode}$ between 4.1 and 6.1 Å, iii) decreasing and becoming
strongly repulsive for $d_{electrode}$ smaller than 4.1 Å. The distance from
the electrode at which a maximum ion density is observed corresponds to the
distance at which the total force cancels out. At this distance, the force is
neither repulsive nor attractive which explains this maximum of density. In
addition, it should be noted that the regime change observed at
$d_{electrode}\sim 4.1\,\textrm{{\AA}}$, coincides with the time at which the
ions have, on average, less than one molecule in their field of vision. This
shows that a vacant space must be available to facilitate the adsorption
process. Additional plots showing the results of the same analysis with TFSI
and the statistical uncertainty from the average curves are available in
Figures S6, S7 and S8.

The kinetic energy of an ion is
$E_{K}=\frac{1}{2}m\left(v_{x}^{2}+v_{y}^{2}+v_{z}^{2}\right)$ where
$v_{\alpha}$ correspond to the different components of the ion velocity. For
$d_{electrode}\gtrsim 4.1\,\textrm{{\AA}}$ and $d_{electrode}\lesssim
8\,\textrm{{\AA}}$, $v_{x}^{2}$ and $v_{y}^{2}$ are close to zero, so the
kinetic energy behaves like $v_{z}^{2}$. Figure S9 shows that the component of
the total force along the $x$-axis and the $y$-axis is zero, which is an
additional argument for looking only along the $z$-axis.

Figure 7: $z$-component of the mean total force felt by the adsorbing EMIM (in
magenta) and a quantity proportional to the kinetic energy
($\frac{1}{2}\left\langle v_{z}\right\rangle^{2}$, in black) as a function of
the distance from the electrode for $\Delta\Psi=$1V. The black dashed line
correspond to the distance from the electrode where
$\mathbf{F_{ij}^{tot}.e_{z}}=0$.

For this reason, on Figure 7 we compare the total force felt by the ion with
$v_{z}^{2}$ as a function of the ion-electrode distance. Figure 7 shows that
the kinetic energy reaches its maximum at the distance for which the force is
the most attractive, i.e. for $d_{electrode}\approx z_{ads}.$ In the
neighborhood of the distance for which $\mathbf{F_{ij}^{tot}.e_{z}}=0$, the
force does not modify the kinetic energy, the work of the force is zero: the
ions diffuse in the plane of the electrode.

Figure 8: $z$-component of the van der Waals (in blue) and Coulomb forces (in
red) as a function of the distance from the electrode for $\Delta\Psi=$1V.

Figure 8 shows the breakdown of the contributions of the van der Waals and
Coulomb forces to the adsorption process. The Coulomb force is attractive over
the whole set of distances, which is due to the presence of an opposite charge
on the surface of the electrode. At long distances this attractive Coulomb
force is the major contribution to the total force and it drives the
adsorption of the ion. At short distances a larger contribution from the van
der Waals force is observed, which is due to the Pauli repulsion between the
electronic clouds, and prevent the ion to approach too close to the electrode.
Finally, it is worth noting that in the intermediate region, for
$d_{electrode}$ between 4.1 and 6.1 Å the two contributions are positive and
almost equal.

Figure 9 shows the van der Waals force felt by the adsorbing ions as a
function of time for $\Delta\Psi=1,\,2,\,4~{}\mathrm{V}$. It can be seen that
before being adsorbed, the forces experienced by the ions do not exhibit any
variation with respect to the applied potential and remain consistently
attractive. This is expected since this force, by definition, is not
contingent upon charge.

Figure 9: Average van der Waals force felt by EMIM (a) and TFSI (b) ions
during adsorption over time for $\Delta\Psi=1,\,2,\,4~{}\mathrm{V}$. Note that
the y-axis scales are different.

Once the ions are adsorbed, the force is essentially repulsive but the
standard deviations of the curves are such that it is difficult to draw
conclusions. Indeed, as seen in Figure 8, once the ions are adsorbed, for
small variations in distance there is a large variation in force.
Nevertheless, if one only considers the mean value, it seems that the
magnitude of the repulsive forces increases with the applied potential
difference. This is due to the ion-electrode distance decreasing with an
increasing applied potential, which consequently increases the force.
Furthermore, the fact that $\varepsilon_{\mathrm{EMIM}}$ is smaller than
$\varepsilon_{\mathrm{TFSI}}$ can explain why
$\left|\mathbf{F_{EMIM}^{LJ}.e_{z}}\right|$ is smaller than
$\left|\mathbf{F_{TFSI}^{LJ}.e_{z}}\right|$.

Figure 10 shows the Coulomb force felt by the adsorbing ions for
$\Delta\Psi=1,\,2,\,4~{}\mathrm{V}$. For all applied potentials the forces are
attractive and increase at the time of adsorption. Furthermore, the
attractiveness of the force increases with the potential. This observation is
consistent with the idea that when the potential increases, the charge carried
by the electrodes also increases. Consequently, this increased charge
amplifies the Coulomb force experienced by the ions, intensifying the
attraction between the ions and the electrodes.

Figure 10: Average Coulomb force felt by EMIM (a) and TFSI (b) ions during
adsorption over time for potentials
$\Delta\Psi=1,2~{}\mathrm{and}~{}4~{}\mathrm{V}$.

The force felt by the cations is greater compared to the force felt by the
anions. This can be explained by the fact that after adsorption, the centres
of mass of the cations are on average closer to the electrode than for the
anions (see Figure S10). It is worth noting that the quantitative results
presented here may depend on the nature of the neighborhood considered.

The comparison of the average absolute values of forces shown in Figure 11
reveals that, at equilibrium, the van der Waals forces exhibit greater
magnitude than the Coulomb forces. This is consistent with the recent findings
reported by de Araujo Chagas et al with graphene/graphyne electrodes and
[EMIM][BF4] electrolyte13.

Figure 11: Average of the absolute van der Waals force felt by EMIM and TFSI
(a) and average of the absolute Coulombic force felt by EMIM and TFSI (b)
after adsorption for $\Delta\Psi=1,\,2,\,4~{}\mathrm{V}$. Averages were done
after ions are adsorbed, for times $t^{\prime}\geq 1$..

This figure also shows clearly the dependance of the forces on the applied
potential.

##  4 Conclusion

The literature concerning the description, through molecular dynamics
simulations, of the EDLs formed by ionic liquids in supercapacitors focused
mostly on equilibrium results so far. Most of the studies present the ion
densities and ion orientations for different electrodes and electrolytes in
order to compare the capacitive properties of the simulated devices. However,
an understanding of the adsorption and charging dynamics of supercapacitors is
needed to improve the performance of such systems. So we have developed
methods for generalising and dynamically describing the ion adsorption
process.

In this work, we focused on a typical supercapacitor, made of two graphite
electrodes separated by pure ionic liquid [EMIM][TFSI] for potential
differences of $\Delta\Psi=1,\,2,\,4~{}\mathrm{V}$. Using molecular dynamics,
we studied the dynamics of EDL creation by describing the orientation of
adsorbing ions in relation to their distances from the electrode, and by
analyzing the forces acting on them.

Firstly, we have seen that EMIM and TFSI ions reorientate as they adsorb,
regardless of the potential applied. The closer the ions are to the electrode,
the more they tend to be parallel to it. We have seen that this tendency
varies linearly with the distance from the electrode, i.e. reorientation
occurs simultaneously to the adsorption.

Secondly, we described and detailed the contributions of the van der Waals and
Coulomb forces acting on the ions during adsorption. We have seen that Coulomb
forces dominate at long range while van der Waals forces dominate at short
range. In addition, we saw that there was an almost equal contribution from
the two forces at an intermediate distance. Interestingly, at short range,
both contributions can depend on the potential applied, as a consequence of
the reduced ion-electrode distance.

In the future, the method developed here could be applied to more complex
materials, such as MoS2, 28 MXenes, 29 etc, for which the presence of a
variety of elements may change the balance between the main interactions,
leading to different adsorption processes. The method could further be used to
explore the case of nanoporous materials, but it would be necessary to account
for their complex geometries by adapting dedicated tools 30 in particular when
detailing the forces experienced by the ions during adsorption.

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation program (grant
agreement no. 714581). This work was granted access to the HPC resources of
TGCC under the allocations AD010911061R1 and A0110911061 made by GENCI. The
authors acknowledge Camille Bacon for initial simulation files.

Supplementary Information available: Additional details on the methods for
defining adsorption events and finding a common time scale, additional details
on the definitions of energy and coarse-grained Lennard-Jones parameters,
evolution of the total charge as a function of time, raw data for the
orientation of ions and ion-electrode distance, pair distribution functions
between ions, forces for TFSI ions and additional data on forces, average ion-
electrode distances.

## Data availability

The data corresponding to the plots reported in this paper, as well as example
input files for the MD simulations, are available in the Zenodo repository
with the identifier 10.5281/zenodo.10695593.

## References

  * Miller and Simon 2008 Miller, J. R.; Simon, P. Electrochemical Capacitors for Energy Management. _Science_ 2008, _321_ , 651–652
  * Griffin et al. 2015 Griffin, J.; Forse, A.; Tsai, W.-Y.; Taberna, P.-L.; Simon, P.; Grey, C. In Situ NMR and Electrochemical Quartz Crystal Microbalance Reveal the Structure of the Electric Double-Layer in Supercapacitor Electrodes. _Nature Mater._ 2015,
  * Péan et al. 2014 Péan, C.; Merlet, C.; Rotenberg, B.; Madden, P. A.; Taberna, P.-L.; Daffos, B.; Salanne, M.; Simon, P. On the Dynamics of Charging in Nanoporous Carbon-Based Supercapacitors. _ACS Nano_ 2014, _8_ , 1576–1583
  * Xu et al. 2020 Xu, K.; Shao, H.; Lin, Z.; Merlet, C.; Feng, G.; Zhu, J.; Simon, P. Computational Insights into Charge Storage Mechanisms of Supercapacitors. _Ener. Environ. Mater._ 2020, _3_ , 235–246
  * Park and McDaniel 2020 Park, S.; McDaniel, J. G. Interference of Electrical Double Layers: Confinement Effects on Structure, Dynamics, and Screening of Ionic Liquids. _J. Chem. Phys._ 2020, _152_ , 074709
  * Jeanmairet et al. 2022 Jeanmairet, G.; Rotenberg, B.; Salanne, M. Microscopic Simulations of Electrochemical Double-Layer Capacitors. _Chem. Rev._ 2022, _122_ , 10860–10898
  * Fedorov and Kornyshev 2008 Fedorov, M. V.; Kornyshev, A. A. Towards Understanding the Structure and Capacitance of Electrical Double Layer in Ionic Liquids. _Electrochim. Acta_ 2008, _53_ , 6835–6840
  * Fedorov et al. 2010 Fedorov, M. V.; Georgi, N.; Kornyshev, A. A. Double Layer in Ionic Liquids: The Nature of the Camel Shape of Capacitance. _Electrochem. Commun._ 2010, _12_ , 296–299
  * Kirchner et al. 2013 Kirchner, K.; Kirchner, T.; Ivaništšev, V.; Fedorov, M. V. Electrical Double Layer in Ionic Liquids: Structural Transitions from Multilayer to Monolayer Structure at the Interface. _Electrochim. Acta_ 2013, _110_ , 762–771
  * Merlet et al. 2014 Merlet, C.; Limmer, D. T.; Salanne, M.; van Roij, R.; Madden, P. A.; Chandler, D.; Rotenberg, B. The Electric Double Layer Has a Life of Its Own. _J. Phys. Chem. C_ 2014, _118_ , 18291–18298
  * Vatamanu et al. 2012 Vatamanu, J.; Borodin, O.; Bedrov, D.; Smith, G. D. Molecular Dynamics Simulation Study of the Interfacial Structure and Differential Capacitance of Alkylimidazolium Bis(trifluoromethanesulfonyl)imide [Cnmim][TFSI] Ionic Liquids at Graphite Electrodes. _J. Phys. Chem. C_ 2012, _116_ , 7940–7951
  * de Araujo Chagas et al. 2022 de Araujo Chagas, H.; Fileti, E. E.; Colherinhas, G. A Molecular Dynamics Study of Graphyne-Based Electrode and Biocompatible Ionic Liquid for Supercapacitor Applications. _J. Mol. Liq._ 2022, _360_ , 119494
  * de Araujo Chagas et al. 2023 de Araujo Chagas, H.; Fileti, E. E.; Colherinhas, G. Comparing Supercapacitors with Graphene/Graphyne Electrodes and [Bmim][PF6], [Emim][BF4], [Ch][Gly] and [Pyr][Tfsi] Ionic Liquids Using Molecular Dynamics. _J. Mol. Liq._ 2023, _379_ , 121703
  * Jiang et al. 2014 Jiang, X.; Huang, J.; Zhao, H.; Sumpter, B. G.; Qiao, R. Dynamics of Electrical Double Layer Formation in Room-Temperature Ionic Liquids Under Constant-Current Charging Conditions. _J. Phys.: Condens. Matter_ 2014, _26_ , 284109
  * Kong et al. 2015 Kong, X.; Lu, D.; Liu, Z.; Wu, J. Molecular Dynamics for the Charging Behavior of Nanostructured Electric Double Layer Capacitors Containing Room Temperature Ionic Liquids. _Nano Res._ 2015, _8_ , 931–940
  * Noh and Jung 2019 Noh, C.; Jung, Y. Understanding the Charging Dynamics of an Ionic Liquid Electric Double Layer Capacitor via Molecular Dynamics Simulations. _Phys. Chem. Chem. Phys._ 2019, _21_ , 6790–6800
  * Marin-Laflèche et al. 2020 Marin-Laflèche, A.; Haefele, M.; Scalfi, L.; Coretti, A.; Dufils, T.; Jeanmairet, G.; Reed, S. K.; Serva, A.; Berthin, R.; Bacon, C. et al. MetalWalls: A Classical Molecular Dynamics Software Dedicated to the Simulation of Electrochemical Systems. _J. Open Source Softw._ 2020, _5_ , 2373
  * Coretti et al. 2022 Coretti, A.; Bacon, C.; Berthin, R.; Serva, A.; Scalfi, L.; Chubak, I.; Goloviznina, K.; Haefele, M.; Marin-Laflèche, A.; Rotenberg, B. et al. MetalWalls: Simulating Electrochemical Interfaces Between Polarizable Electrolytes and Metallic Electrodes. _J. Chem. Phys._ 2022, _157_ , 184801
  * Reed et al. 2007 Reed, S. K.; Lanning, O. J.; Madden, P. A. Electrochemical Interface Between an Ionic Liquid and a Model Metallic Electrode. _J. Chem. Phys._ 2007, _126_ , 084704
  * Canongia Lopes et al. 2004 Canongia Lopes, J. N.; Deschamps, J.; Pádua, A. Modeling Ionic Liquids Using a Systematic All-Atom Force Field. _J. Phys. Chem. B_ 2004, _108_ , 11250–11250
  * Canongia Lopes and Pádua 2004 Canongia Lopes, J. N.; Pádua, A. Molecular Force Field for Ionic Liquids Composed of Triflate or Bistriflylimide Anions. _J. Phys. Chem. B_ 2004, _108_ , 16893–16898
  * Bhargava and Balasubramanian 2007 Bhargava, B. L.; Balasubramanian, S. Refined Potential Model for Atomistic Simulations of Ionic Liquid [bmim][PF6]. _J. Chem. Phys._ 2007, _127_ , 114510
  * Chaban 2011 Chaban, V. Polarizability Versus Mobility: Atomistic Force Field for Ionic Liquids. _Phys. Chem. Chem. Phys._ 2011, _13_ , 16055–16062
  * Zhang and Maginn 2012 Zhang, Y.; Maginn, E. J. A Simple AIMD Approach to Derive Atomic Charges for Condensed Phase Simulation of Ionic Liquids. _J. Phys. Chem. B_ 2012, _116_ , 10036–10048
  * Schröder 2012 Schröder, C. Comparing Reduced Partial Charge Models with Polarizable Simulations of Ionic Liquids. _Phys. Chem. Chem. Phys._ 2012, _14_ , 3089–3102
  * Martínez et al. 2009 Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. _J. Comput. Chem._ 2009, _30_ , 2157–2164
  * Trenzado et al. 2021 Trenzado, J. L.; Rodríguez, Y.; Gutiérrez, A.; Cincotti, A.; Aparicio, S. Experimental and Molecular Modeling Study on the Binary Mixtures of [EMIM][BF4] and [EMIM][TFSI] Ionic Liquids. _J. Mol. Liq._ 2021, _334_ , 116049
  * Bi and Salanne 2022 Bi, S.; Salanne, M. Co-Ion Desorption as the Main Charging Mechanism in Metallic 1T-MoS2 Supercapacitors. _ACS Nano_ 2022, _16_ , 18658–18666
  * Xu et al. 2020 Xu, K.; Merlet, C.; Lin, Z.; Shao, H.; Taberna, P.-L.; Miao, L.; Jiang, J.; Zhu, J.; Simon, P. Effects of Functional Groups and Anion Size on the Charging Mechanisms in Layered Electrode Materials. _Energy Storage Mater._ 2020, _33_ , 460–469
  * Dick and Kirchner 2023 Dick, L.; Kirchner, B. CONAN – Novel Tool to Create and Analyze Liquids in Confined Space. _J. Chem. Inform. Model._ 2023, _63_ , 6706–6716