# The intrinsic electrostatic dielectric behaviour of graphite anodes in Li-
ion batteries – across the entire functional range of charge

Simon Anniés  Christoph Scheurer  Chiara Panosetti

###### Abstract

Lithium-graphite intercalation compounds (Li-GICs) are by far the most common
anode material for modern Li-ion batteries. However, the dielectric response
of this material in the electrostatic limit (and its variation depending on
the state of charge (SOC)) has not been investigated to a satisfactory degree
– neither by means of theory nor by experiment – and especially not for the
higher range of SOC. Nevertheless, said dielectric behaviour is a highly
desired property, particularly as an input parameter for charged kinetic Monte
Carlo simulations – one of the most promising modeling techniques for energy
materials. In this work, we make use of our recently published DFTB
parametrization for Li-GICs based on a machine-learned repulsive potential in
order to overcome the computational hurdles of sampling the long-ranged
Coulomb interactions within this material – as experienced by the charge
carriers within. This approach is rather novel due to computational cost, but
best suited for investigating our specific property of interest. For the first
time, we discover a mostly linear dependency of the relative permittivity
$\epsilon_{r}$ on the SOC, from $\approx 7$ at SOC 0% to $\approx 25$ at SOC
100%. In doing so, we also present a straightforward approach that can be used
in future research for other intercalation compounds – once sufficiently fast
and long-ranged computational methods – such as linear-scaling DFT, a good
DFTB parametrization, or atomic potentials with inbuilt electrostatics –
become available. However, while the presented qualitative behaviour is robust
and our results compare favourably with the very few experimental studies
available, we do stress that the quantitative results are strongly dependent
on our estimation of the partial charge transfer from the intercalated Li-ion
to the carbon host structure, and need to be verified by further experiments
and other calculations. Yet, our research shows that in principle, two
measurements – one at low and one at high SOC – should suffice for that
purpose.

###### keywords:

relative permittivity , dielectric response , graphite anodes , lithium
intercalated graphite , multiscale modeling

###### PACS:

0000 , 1111

###### MSC:

0000 , 1111

††journal: Electrochimica Acta

[inst1]organization=Fritz-Haber-Institute of the Max-Planck-
Society,addressline=Faradayweg 4-6, postcode=14195, city=Berlin,
country=Germany

[inst2]organization=Chair for Theoretical Chemistry, TU
Munich,addressline=Lichtenbergstr. 4, postcode=85747, city=Garching b.
München, country=Germany

The relative permittivity of lithium intercalated graphite anodes is linearly
dependent on the state of charge, with approximate values of around 7 at 0%
and around 25 at 100%.

Sampling the Coulomb interaction between two charge carriers within the
material is enabled by our recently published DFTB parametrization for lithium
intercalated graphite.

The dielectric screening in graphite is significantly larger in the direction
perpendicular to the graphene-sheets, than in the direction parallel to them.
This is important when comparing experimental results for graphite powder and
for perfect graphite crystals.

##  1 INTRODUCTION:

The relative permittivity (RP) is one of the defining properties of many
materials, as it describes the degree to which the Coulomb interactions
between charge carriers are screened within the material compared to vacuum.
Especially in older literature, this property is also known as ‘‘dielectric
constant’’, even though it is far from being constant, but dependent on
temperature, the frequency of probing electric fields and even the underlying
mechanisms in terms of polarizability, conductivity and others.

Traditionally, the RP has primarily been of interest for insulators, but in
recent times it has also been increasingly investigated for conducting
materials [1], where it stems from the interaction between a small fraction of
the charge carriers and the atoms. The fact that these charge carriers are
mobile in conducting materials fundamentally changes the way the property and
its dependencies need to be understood within these materials, compared with
insulators.

One group of materials of great interest are lithium-graphite intercalation
compounds (Li-GICs), which constitute by far the most common anode material in
modern lithium ion batteries. During the charging and discharging cycles of
the battery, Li-ions are stored between the layers of the graphite host
structures, up to a stoichiometry of LiC6, which traditionally translates to a
state of charge (SOC) of $100\%$ and corresponds to one Li-ion above every
third C6 ring of the hexagonal base lattice – even though recent studies have
shown that overlithiation beyond that point is possible at ambient conditions
[2]. The distribution of the Li-ions for intermediate SOCs is not uniform, but
ordered as shown in fig. 1, as first explained in [3].

Figure 1: Illustration of lithium-ions (purple spheres) intercalated into a
graphite host structure (grey lines) at $50\%$ state of charge, based on the
Daumas-Herold domain model [3]. Li-ions tend to fill up every second layer
completely (staging, left), before starting to intercalate into the other half
of the layers. However, this behaviour is not global, but occurs in finite-
sized domains (right) and is not expected to be perfect in real systems.
During the intercalation process, the distance between the graphene sheets
(interlayer distance) is increased by around $10\%$.

Within the scope of this specific material, we define the RP we are
investigating as ‘‘the damping of the electrostatic interaction between two
Li-ions or Li-ion vacancies embedded in the material caused by the surrounding
charge-carrier density’’. As of note, this is the electrostatic, low-frequency
RP, as opposed to what is measured in many experiments, which make use of
alternating AC-fields at a vast variety of frequencies [1, 4]. Furthermore, we
point out that this specific property is directionally separated – its
contribution in the $xy$-plane (in this work defined as parallel to the
graphene sheets) is expected to be different from the contribution in
$z$-direction (orthogonal to the graphene sheets). This makes direct
comparison with experiments performed on graphite powder as opposed to a
‘‘perfect’’ crystal non-trivial – a problem we address in section 3.4.

One of the primary motivations for investigating the RP of Li-GICs is the fact
that it is a required input parameter for including charge in kinetic Monte
Carlo simulations (kMCs), which in turn are a crucial method for studying
charge carrier dynamics in functional energy materials [5, 6, 7, 8]. Long-
ranged Coulomb interactions are a necessity when performing kMC on systems,
which include charged or partially charged particles. For example, Casalegno
et al. [9] have shown that not including such interactions (as it would be the
case when using e.g. force field approaches with finite-size descriptors)
causes an error of 14% in the protonic diffusion coefficients in doped
perovskites. The situation becomes significantly more complicated when looking
at anode and cathode materials, but also electrolytes, perovskites in solar
cells, and any other type of functional materials that involves changes in the
density and/or local ordering of charge carriers as part of their intended
function. This is due to the fact that the local relative permittivity then is
not only influenced by the ‘‘host’’ material, but also by the charge carriers
close by. Therefore, the RP of Li-GICs changes significantly depending on the
SOC (vide infra).

In this work, we put forward a systematic approach to investigating this
crucial property. Based on this, we determine for the first the relative
permittivity of Li-GICs, as a function of the SOC, for the entire functional
range of the material.

Beyond the previously outlined interest for charge-kMC, we believe there to be
many more valuable applications for the relative permittivity of Li-GICs also
at higher levels of the multiscale simulation hierarchy: one such motivation
is understanding charge gradients, as they occur during the fast charging of
modern batteries in electric vehicles, and the chemical pressure which leads
to plating and dendrite formation inside the batteries under certain
operational conditions. The latter phenomenon is typically investigated by
means of continuum simulations like e.g. by Hein et al. [10], which also rely
at least implicitly on knowledge of the dielectric response. Furthermore, a
simple model of the charge carrier electrostatics could be used as a physical
baseline for otherwise short-ranged machine learning models or cluster
expansions. Another related field is the development of functional materials
based on doped graphite [11].

The dielectric behaviour of Li-GICs (and solid materials in general) is
significantly more complicated than the expression ‘‘dielectric constant’’
would suggest, and is governed by vastly different physics at different
frequencies of a probing field. In the static limit (which this work aims to
investigate – the probing field is essentially the electrostatic field of the
intercalated ion itself), no periodic movement (beyond thermal fluctuations)
of the electrons is induced. However, this picture changes in the kHZ range,
where the entire charge carrier density oscillates with the probing field,
causing large polarization and large dielectric screening. For example, Chung
et al. [1] measured an RP of $\epsilon_{r}=2100$ for highly oriented pyrolytic
graphite (HOPG) and even higher ones for other carbon structures, at 2-10KHz.
Moving on to the GHz regime, a balance is reached where the field oscillations
are too fast for macroscopic bulk currents to build up, and a situation occurs
that is arguably similar to the static limit and may serve as comparison for
our research. Hotta et al. [4] put forward a dielectric constant for graphite
powder of $\epsilon_{r}\approx 15$, at 6GHz. Finally, at even higher
frequencies beyond THz, the electric field becomes high enough in energy to
excite a significant number of electrons, again creating a physically
different situation with much lower dielectric screening, which converges to
transparency in the $\omega\rightarrow\infty$ limit. A study by Jellison et
al. [12] in the frequency regime of visible light finds an RP in $xy$-plane
for HOPG of $\epsilon_{r}=4.21$, which, for the previously mentioned reasons,
cannot be used as comparison either and is expected to serve as a lower bound
in the following.

It is apparent that there is a glaring lack of studies investigating the exact
property of interest to us, which is – again – the electrostatic dielectric
response of a perfect graphite crystal in $xy$-plane, i.e. parallel to the
graphene sheets, as experienced by some internal charge carriers (in this case
Li-ions and vacancies). There are some studies available on graphene, either
on some substrate or quasi-freestanding, with results ranging from
$\epsilon_{r}=2.2-5.0$ by Elias et al. [13] to $\epsilon_{r}=15.4$ by Reed et
al. [11], and another study by Bostwick et al. finding $\epsilon_{r}\approx
4.4$ [14], none of which can serve as direct comparison to our research
either. However, there is a study on bilayer graphene (which according to our
calculations can be compared with graphite quite well) by Bessler et al. [15],
putting forward an RP of $\epsilon_{r}=6\pm 2$. This is likely the most
reliable direct experimental comparison currently available to us.

In terms of theoretical approaches to determining the dielectric response of
materials, substantial work has been done on water [16, 17, 18, 19]. There is
also some promising work by Gigli et al. [20] in the development of an
integrated machine learning model predicting the dielectric response of
$BaTiO_{3}$. However, all these methods are reliant on the presence of
polarizable dipoles within the system, which is not the case for Li-GICs.

##  2 METHODOLOGY:

###  2.1 Computational details:

For this study, we used self-consistent-charge Density Functional Tight
Binding (SCC-DFTB [21]) as implemented in in DFTB+ [22], with the
parametrization developed in our group. The corresponding Slater-Koster files
are available upon request and have been parametrized and tested as described
in [23, 24]. The repulsion potential is machine-learned by means of Gaussian
Process Regression (GPR) [25].

Geometries have been constructed and analysed by means of the Atomic
Simulation Environment (ASE [26]) which we also used as a base framework for
all force- and energy-calculations, structure relaxations (specifically using
the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm as an optimizer [27]),
and transition state calculations. For the latter, we employed the Nudged
Elastic Band (NEB, [28]) algorithm with the BFGS optimizer and climbing image
[29, 30] switched on.

For all DFTB calculations, we used a well converged k-point density of at
least 0.1/Å for the $z$-component of the unit cell. The $xy$ size of the cell
is large enough to sample at the Gamma point. The SCC-tolerance is $10^{-6}$.
We employed Fermi filling with a Fermi temperature of 0.001 Kelvin, as well as
a Broyden mixer [31] for convergence acceleration with a mixing parameter of
0.5. All of these settings have been tested with regard to convergence for the
whole range of SOC. As described in [23], our parametrization is meant to be
used with the Lennard-Jones dispersion correction [32] switched on.

In terms of supercell size convergence, satisfactory convergence of the
extracted slopes of the energies relative to the inverse distances (see
section 2 for details) is reached at distances of around $15\text{\r{A}}$
between the periodic images of the sampled areas of the investigated layers
(see fig. 2 and fig. 4), a value we reach or exceed with all used supercells.
We point out that convergence of total energies is still not reached at those
distances due to the long-ranged nature of the Coulomb interactions, however
our property of interest – the previously mentioned slope – is rather robust
to total energy shifts of the whole sampled area. Slight further improvements
would still be probable with even larger supercells, but due to computation
time constraints, and due to the fact that other effects introduce much larger
errors to the final results, we chose the supercells described in the
following sections.

###  2.2 Derivation of the relative permittivity $\epsilon_{r}$ from the
Coulomb law and discussion of the partial charge transfer:

In order to determine the relative permittivity within the system, our
approach is to sample the electrostatic interactions between two charge
carriers inside the system, placed at varying distances from each other. We
begin our considerations with the electrostatic energy $E_{Coul}$ of two
charge densities $\rho_{1}$ and $\rho_{2}$, governed by the Coulomb law

$\displaystyle
E_{Coul}=\frac{1}{4\pi\epsilon_{0}\epsilon_{r}}\cdot\int\int\frac{\rho_{1}\cdot\rho_{2}}{|\vec{r}_{1}-\vec{r}_{2}|}dr_{1}dr_{2}$
(1)

where $\epsilon_{0}$ is the vacuum permittivity. In this work, we chose to
approximate the sampled charge carriers as point charges, an approximation
that is unproblematic for Li-ions, but requires a bit more consideration for
Li-ion vacancies, which may be more diffuse in shape. In this approximation,
the expression becomes

$\displaystyle E_{Coul}=\frac{1}{4\pi\epsilon_{0}}\cdot\frac{Q_{1}\cdot
Q_{2}}{\epsilon_{r}}\cdot\frac{1}{|\vec{r}_{1}-\vec{r}_{2}|}$ (2)

which has been rearranged to clearly separate the continuum electrostatic from
the geometric quantities. With the two point charges being Li-ions or Li-ion
vacancies intercalated into the graphite host structure, we get

$\displaystyle
E_{Coul}=\frac{e^{2}}{4\pi\epsilon_{0}}\cdot\frac{Z^{2}_{Li}}{\epsilon_{r}}\cdot\frac{1}{R}$
(3)

where $e$ is the electron charge, $Z_{Li}$ the partial charge on the Li-ion
and $R$ the distance between the two charge carriers. In order to link this
Coulomb energy to the potential energy $E_{pot}$ of an entire supercell (as we
obtain from our DFTB calculations), an appropriate reference energy $E_{0}$
needs to be introduced, which is the host structure energy in the limit of
$R\rightarrow\infty$ at that specific stoichiometry. With this we finally get:

$\displaystyle
E_{pot}=\frac{e^{2}}{4\pi\epsilon_{0}}\cdot\frac{Z^{2}_{Li}}{\epsilon_{r}}\cdot\frac{1}{R}+E_{0}$
(4)

The first two terms of this expression can be accessed as the slope of a
linear regression, when plotting $E_{pot}$ over $1/R$. In order to relate said
slope to $\epsilon_{r}$, we approximate $E_{Coul}\approx E_{pot}-E_{0}$. While
this obviously holds in the macroscopic limit for a homogeneous medium, here
we are attempting to measure the Coulombic repulsion felt by two particles on
a potential energy surface (PES) resolved at the atomistic level. The latter
is not infinitely smooth, but presents minima, maxima and saddle points due to
the carbon host structure. It however shows a sufficiently regular pattern (as
shown in [24]) to assume that analogous points (i.e., comparing minima with
minima etc.) exhibit analogous local shape. Consistently, the DFTB total
energy is indeed expressed (within the formalism of the method [33]) as a sum
of three contributions: the so-called band structure energy, the Coulomb
energy (the expression of which depends explicitly on the partial charges on
the individual atoms, vide infra), and the repulsive energy. Ideally, we can
isolate the Coulomb energy between two particles in the system by considering
‘‘the rest of the DFTB total energy’’ as a background to subtract pointwise,
assuming that it will be similar enough in similar local environments.
Ideally, this would correspond to the limit for infinitely dilute, unperturbed
PES, where the Coulombic interaction between charge carriers vanishes. In this
light, for each sampled point we may identify a suitable (pointwise) reference
$E_{0}$ as the corresponding point on such unperturbed PES, e.g., in the case
of minima, the total energy of a minimum at the center of any C6 ring of the
graphitic host sufficiently far away from the other charge carrier.

We point out that this approximation neglects distortions of the local
structure and electron density that may be caused by two charge carriers being
close together – an approximation that holds well for the low SOC regime, but
not quite as well for the high SOC regime, as will be shown. The reference
energy $E_{0}$ does not impact that slope, but only causes an up or down shift
along the $y$-axis. Therefore, it is not immediately relevant for the
extraction of the $\epsilon_{r}$, but it is needed for consistent plotting of
the Coulomb energy. For this reason, we estimate $E_{0}$ by interpolating the
slope to the $1/R\rightarrow 0$ limit (see the following chapters for more
details).

The only other variable that is left, then, is the partial charge of the
intercalated lithium $Z_{Li}$. It is known that, upon intercalation, the
electron density of the Li-atom is partially transferred to the carbon host
structure, leaving the intercalant as something usually denominated as
$Li^{+}$, but the exact magnitude of said charge transfer is hard to pinpoint
– and not an observable. The self-consistent-charge cycles of our DFTB+
calculations output a local electron population of around $0.21$ to $0.26$ for
the Li-intercalants depending on the local environment of the respective Li-
ion, which corresponds to a partial charge of $+0.79e$ to $+0.74e$. However,
literature reports a variety of different values. Valencia et al. [34] pointed
out the large dependency of the partial charge transfer on the method of
analysis and presented values of $0.43e$ (Mulliken charge analysis), $0.47e$
(Voronoy), $0.6e$ (Löwdin) and $1.0e$ (Bader). Krishnan et al. [35] found
$0.86e$, also by Bader analysis. Song et al. [36] determined $0.68e$ by means
of quantum mechanical calculations and comparison with experimental layer
spacing, but also pointed out that this value may change with SOC. Finally,
Rakotomahevitra et al. [37] calculated $0.517e$, using an extra-orbital model.
To complicate things further, $\epsilon_{r}$ depends on $Z_{Li}$ quadratically
(see eq. 4) and is therefore very sensitive to it. For consistency, we will
move forward in this work assuming $Z_{Li}=0.765\pm 0.05$ – the median partial
charge to which the SCC converges. The latter indeed directly enters the
expression for the Coulombic contribution to the DFTB total energy, as
mentioned before, and therefore represents a natural choice within the
framework (and within the approximations) of DFTB. We add a generous, but
arbitrary measure of uncertainty, and we stress again that all absolute
numbers for the relative permittivity $\epsilon_{r}$ presented in this work
need to be understood with this assumption in mind. However, the same is not
true for the electrostatic screening (captured by the slopes of the linear
regression) – this can be taken at face value and may be, depending on the
model or method employed, the more valuable property to take into account for
future research.

##  3 RESULTS & DISCUSSION:

###  3.1 Lithium-intercalant pairwise interaction screening ($xy$-plane):

In order to sample the pairwise Li-ion interactions within the material, we
construct an investigated layer (called $inv(Li)$, consisting of 300 carbon
atoms and two Li-ions, one of them fixed in the corner, the other one sampling
a number of different positions, as schematically illustrated in fig. 2. We
then combine this investigated layer with varying stackings of empty and
filled layers (as they predominantly appear in the system according to the
staging model), perform full structure relaxations on the resulting supercells
(with two or three layers and 600-900 carbon atoms) for each of the 21
possible Li-ion positions (per data point on the SOC axis) and extract the
potential energies of the cells. We also perform nudged elastic band (NEB)
calculations between all combinations of two relaxed structures, which have
directly adjacent occupied Li-positions, thus acquiring 41 transition states
(per data point) of diffusive next-neighbour jumps. This is especially
relevant for kinetic Monte Carlo applications, where the dynamics is governed
by the rates (and therefore, the activation barriers) corresponding to
elementary processes bringing the system from one state to another. As our
results show, the Coulombic behaviour we get is identical (within method
accuracy) for both cases (fig. 2). This justifies the application of
simplified rescaling rules solely based on Coulombic interactions to include
the effect of next-neighbour occupations on the elementary barriers, allowing
to set up a charged kMC model with one, simple, elementary barrier (the jump
of one Li-ion from one site to another in the infinitely dilute limit), and
correcting the latter depending on the number and distance of nearest
neighbours as well as on the direction of the jump.

Figure 2: Illustration of the investigated layer for pairwise Coulombic Li-
ion-interactions with one fixed Li-ion ($purple$). Of the locations marked
with an ‘‘X’’, one is occupied by, the other Li-ion, while all others are
empty. This layer will be referred to as $inv(Li)$ throughout this work.

Figure 3: Approximate Coulomb energies $E_{Coul}\approx E_{pot}-E_{0}$ of minima and transition states relative to $E_{0}$ (at the respective SOC) in dependence of the distance (left) and inverse distance (right) between the investigated Li-ions at different SOCs. The slopes have been extracted by means of linear regression. The global minima $E_{0}$ in the $1/R\rightarrow 0$ limit are extrapolated from the linear regression of the slopes of the minima and set to $0$ in the plot. The DFT reference is taken from an Ising model based on BEEF-vdW DFT [38]. The data points represent the pair-interactions for next-neighbour, second, third and fourth next neighbour lithium positions. The SOC is not clearly defined without giving the size of the otherwise empty supercell (a limitation of the Ising model), but can be understood as ‘‘low’’ and taken as comparison for our blue data-points. Structure | Stoichiometry | SOC | Slope [eVÅ] | rel. permittivity
---|---|---|---|---
$inv(Li)-empty$ |  Li2C_600 | 2.0% | $1.168\pm 0.013$ | $7.23\pm 0.86$
$bilayer\ graphene$ |  Li2C_600 | 2.0% | $1.111\pm 0.018$ | $7.55\pm 1.04$
$inv(Li)-empty-full$ |  Li_52C_900 | 34.7% | $0.759\pm 0.018$ | $11.32\pm 1.74$
$inv(Li)-full$ |  Li_52C_600 | 52.0% | $0.407\pm 0.026$ | $21.04\pm 4.05$
Table 1: Results for the data points constructed with an otherwise empty
investigated layer with sampled lithium positions inv(Li). The slopes and
their RMSEs stem from linear regression of all data points (minima and
transition states) – shifted to the same baseline. The corresponding relative
permittivities have been calculated via the Coulomb law in eq. 4.

By plotting the estimated Coulomb energies $E_{Coul}=E_{pot}-E_{0}$ as
functions of the (inverse) distance between the two Li-ions in the
investigated layer, we clearly illustrate how our model captures the Coulombic
nature of the interaction close to perfectly in the case where no other Li-
ions are present in the adjacent layers (fig. 3, blue) that could distort the
electron density. The very minor scatter in this case likely stems from the
small distortions in the carbon structure close to the Li-intercalants and the
fact that our full cell relaxations cannot be converged to infinity, but have
to be stopped at some threshold forces. In the other two cases, there is some
more scatter present (caused by the slightly deformed charge density, due to
the filled adjacent layers (one in case of SOC 34.7%, two in case of SOC 52%))
but the overall behaviour is still predominantly Coulombic.

In a next step, we extract the slopes from the $1/R$-plots by means of linear
regression (see table 1). We achieve this by shifting the transition state
energy levels down to the ground states and then fit all data points at once.
Based on that, the relative permittivity $\epsilon_{r}$ can be extracted from
the Coulomb law in eq. 4. For the sake of comparison with experiment, we also
perform the same procedure with freestanding bilayer graphene, consisting of
one $inv(Li)$ layer between 2 graphene sheets, which according to our model,
has a very similar $\epsilon_{r}$ as (periodic) graphite at the same
stoichiometry.

As pointed out in fig. 3, there is one DFT study available for comparison by
Pande et al. [38], and the slope they find at low SOC agrees well with ours.
However, due to the computational cost of DFT, they were only able to provide
4 data points, and only at the low end of the range of charge, which is the
cheapest to compute.

Furthermore, Bessler et al. measured a relative permittivity of
$\epsilon_{r}=6\pm 2$ for bilayer graphene, which agrees quite well with our
result of $\epsilon_{r}=7.55\pm 1.04$.

###  3.2 Vacancy pairwise interaction screening ($xy$-plane):

It is clear that it is not possible to investigate SOCs close to 100%
following the same approach, i.e. using the same investigated layer described
in the previous chapter. Furthermore, for the higher SOCs, the diffusion
mechanism transitions towards vacancy hopping instead of Li-intercalant
hopping. Therefore, we introduce a second type of investigated layer (fig. 4),
which samples the pairwise Coulomb interactions between Li-vacancies within a
filled layer instead (called $inv(Vac)$). A key difference here is the fact
that, in a filled layer, only every third C6 ring is occupied by a Li-ion, so
the investigated vacancies cannot be placed on each $C_{6}$-ring, but only on
every third one. Due to this, a slightly larger cell is needed (384 carbon
atoms) in order to sample a reasonable number of data points (15 per SOC),
while also keeping the separation from the periodic image large enough to be
well converged.

Figure 4: Illustration of the investigated layer for pairwise Coulombic
vacancy-interactions with fixed Li-ions (purple) and fixed vacancy #1 (bottom
left, dashed circle). Of the positions marked with an ‘‘X’’, the sampled
vacancy #2 is located on one, while all others are occupied by Li-ions. This
layer will be referred to as $inv(Vac)$ throughout this work.

Another consequence of this is the fact that the diffusion path from one
vacancy location to another is not clearly defined – it may be a straight line
or pass through one of the two next-neighbour minima in between. Because of
this and because we already showed before that the slopes extracted from the
minima and from the diffusion path intermediates are very similar, we limit
our calculations to the minima in this chapter. Other than that, the procedure
is the same as in section 3.1.

Figure 5: Approximate Coulomb energies $E_{Coul}\approx E_{pot}-E_{0}$ of local minima relative to $E_{0}$ (at the respective SOC) in dependence of the distance (left) and inverse distance (right) between the investigated Li-ion vacancies at different SOCs. Slopes have been extracted by means of linear regression. The global minima $E_{0}$ in the $1/R\rightarrow 0$ limit are extrapolated from the linear regression of the slope and set to $0$ in the plot. Structure | Stoichiometry | SOC | Slope [eVÅ] | rel. permittivity
---|---|---|---|---
$inv(Vac)-empty$ |  Li_62C_768 | 48.4% | $0.602\pm 0.047$ | $14.28\pm 2.95$
$inv(Vac)-empty-full$ |  Li_126C_1152 | 65.6% | $0.448\pm 0.041$ | $19.28\pm 4.22$
$inv(Vac)-full$ |  Li_126C_768 | 98.4% | $0.350\pm 0.034$ | $24.70\pm 5.55$
Table 2: Results for the data points constructed with a filled investigated
layer with sampled vacancy positions inv(Vac). The slopes and their RMSEs stem
from linear regression of the data points (minima). The corresponding relative
permittivities have been calculated via the Coulomb law in eq. 4.

As can be seen in fig. 5, there is more distortion to the Coulombic behaviour
present in this case, some of which is even qualitative and cannot just be
attributed to the uneven electron density in the adjacent layers. Especially,
the data points at ‘‘medium’’ distance ($\approx 8\text{\r{A}}$ in the left
plot or $\approx 0.125\ 1/\text{\r{A}}$ in the right plot of fig. 5) have
lower energy than expected if it was purely governed by the Coulomb law. As
pointed out before, this is to be expected within our approximations, since a
combination of multiple effects are at play here, that can not be addressed
without explicit treatment of the overlapping local environments of the charge
carriers: firstly, for these compositions, there is not just vacuum present in
the space between the investigated sites (as it is the case in all
configurations in section 3.1, as well as the closest possible vacancy
positions (see fig. 4). This additional charge density between the sampled
sites seems to have an additional stabilizing effect in our DFTB calculations.
For vacancy pairs at distances larger than $\approx 10\text{\r{A}}$, this
effect probably becomes negligible again, since there is hardly any
interaction energy left to be screened at those distances. Secondly, the
point-charge approximation may be less accurate for a Li-vacancy, than it is
for a Li-ion, since the next-neighbour local minima within the host structure
are not occupied either. And thirdly, upon inspection of the fully relaxed
structures, one realizes that the Li-ions adjacent to the investigated
vacancies are not located in the middle of their respective C6 ring, but
slightly displaced towards the vacancy. Due to these additional effects,
fitting the potential energies found for these configurations to the
unperturbed Coulomb law ($E_{Coul}=E_{pot}-E_{0}$) is a more severe
approximation in the context of the $inv(Vac)$ investigated layer than it was
for the $inv(Li)$ investigated layer. Finally, the slopes and relative
permittivities $\epsilon_{r}$ are summarized in table 2, analogous to section
3.1.

###  3.3 Dilute adjacent layers:

One of the complications when investigating Li-GICs is the fact that the same
stoichiometry can be realized in a variety of different ways, and while it is
known that staged configurations are favoured in equilibrium and in perfect
crystals, other – dilute – configurations may still play a role, when the
system is under the effects of defects, grain boundaries or nonequilibrium
states caused by fast charging. Because of this, we choose to investigate some
of these configurations, as well. For this purpose, we make use of the same
investigated layers as before, but combine them with an adjacent layer, that
is dilutely filled at 50% capacity. In order to realize the periodicity of
said layer, a slightly larger cell with 360 carbon atoms is necessary in the
case of the Li-ion investigated layer.

Figure 6: Approximate Coulomb energies $E_{Coul}\approx E_{pot}-E_{0}$ of local minima relative to $E_{0}$ (at the respective SOC) in dependence of the inverse distance between the investigated sites in the $inv(Li)-dilute$ (left) and $inv(Vac)-dilute$ (right) configurations. The global minima $E_{0}$ in the $1/R\rightarrow 0$ limit are extrapolated from the linear regression of the slope and set to $0$ in the plot. Structure | Stoichiometry | SOC | Slope [eVÅ] | rel. permittivity
---|---|---|---|---
$inv(Li)-dilute$ |  Li_32C_720 | 26.7% | $0.594\pm 0.078$ | $14.03\pm 3.12$
$inv(Vac)-dilute$ |  Li_94C_768 | 73.4% | $0.463\pm 0.152$ | $21.45\pm 9.52$
Table 3: Results for the data points constructed with dilutely filled adjacent
layers. The slopes and their RMSEs stem from linear regression of the data
points (minima). The corresponding relative permittivities have been
calculated via the Coulomb law in eq. 4.

As expected, this dilutely filled adjacent layer causes a larger scatter of
the potential energies of the individual configurations, which is due to the
less even electron density of such a layer compared to a completely empty or
full one. In the $inv(Vac)-dilute$ case, this effect is also combined with the
already larger scatter we observed and explained previously when sampling this
type of investigated layer. We tried to mitigate this by also sampling over
two different orientations of the dilute adjacent layers – averaging out the
uneven electron density – but the scatter remains very large, which leads to a
much larger uncertainty of the resulting slope and relative permittivity
$\epsilon_{r}$ (see table 3).

Nevertheless, the Coulombic nature of the interaction is still clearly visible
despite the scatter at least for the $inv(Li)-dilute$ case (fig. 6, left), and
is expected to be captured by the (admittedly rough) linear fit in the other
case as well – yet of course with a larger uncertainty. Indeed, the resulting
slope is still in line with the overall behaviour we find throughout the whole
system.

###  3.4 Dielectric screening in the $xy$-plane vs. in $z$-direction:

Thus far, we have focused on the dielectric behaviour of Li-GICs parallel to
the graphene sheets ($xy$-plane), because that is what governs the diffusion
and intercalation of the Li-ions. However, real graphite anodes and
experimental samples often come in powder form, with graphite nanoparticles in
random spatial orientation. Experimental measurements of the dielectric
behaviour of such samples cannot be directly compared to the results we have
presented so far. Therefore, we also investigate the dielectric behaviour in
$z$-direction (perpendicular to the graphene sheets).

Figure 7: Approximate Coulomb energies $E_{Coul}\approx E_{pot}-E_{0}$ of
local minima relative to $E_{0}$ (at the respective SOC) in dependence of the
distance (left) and inverse distance (right) between the investigated Li-ions
in $z$-direction.

In order to do so, we construct a 15-layered supercell with a stoichiometry of
Li2C_720 which corresponds to an SOC of 1.7%. We sample 5 of the 15 layers (to
ensure no self-interaction with the periodic image) with one pair of Li-ions,
placed in positions both directly above each other and shifted by one next-
neighbour position (those latter interactions are not perfectly in
$z$-direction, but still ‘‘through’’ the graphene sheets). As can be seen, the
resulting electrostatic behaviour is not purely Coulombic (see fig. 7).

Similar to the results for the investigated layer $inv(Vac)$, the presence of
charge carriers in between the investigated Li positions seems to cause some
extra, nonlinear screening. Holding on to the assumption of a linear trend, a
linear fit finds a slope of $0.411\pm 0.054$, which translates to a relative
permittivity of $\epsilon_{z}=21.35\pm 5.50$ – significantly larger than the
in-plane relative permittivity $\epsilon_{xy}$ we find at similar SOC. Based
on this, we can average over the 3 spatial dimensions in order to provide a
rough estimate of the relative permittivity $\epsilon_{r}$ of randomly
oriented graphite powder at very low SOC, simply by applying:

$\displaystyle\epsilon_{r}=(2\epsilon_{xy}+\epsilon_{z})/3=11.93$ (5)

This result is comparable to the $\epsilon_{r}\approx 15$ found by Hotta et
al. [4] in the GHz range, but again, it is not entirely clear that this
comparison is physically valid. Once a more reliable consensus is reached in
experiment, our estimation of the partial charge can be further validated. It
is also necessary to point out that grain-boundary and grain-size effects may
play a significant role, especially the smaller the particles in the powder
get, but these effects are not considered at all in the numbers we present
here – those are simply intended to help with comparison between our results
and experiments.

###  3.5 Relative permittivity $\epsilon_{r}$ as a function of the SOC:

Putting the previous results together (see fig. 8), we find an approximately
linear dependency of the relative permittivity $\epsilon_{r}$ (in the
$xy$-plane) on the SOC of the Li-GIC. Between 20% and 80% SOC, we observe some
additional effects of the local ordering – dilute or staged – and of whether
Li-ions in the otherwise empty layers or Li-ion vacancies in the otherwise
full layers are sampled. These deviations are largest at around 50% SOC, but
the average $\epsilon_{r}$ of configurations at roughly equal SOC are still
close to the weighted linear regression fit, so the more macroscopic a
viewpoint is taken, the less these local phenomena matter.

Figure 8: Relative permittivity $\epsilon_{r}$ in $xy$-plane (parallel to the
graphene sheets), found for different configurations of the Li-GIC material as
a function of the respective SOC . The weighted linear fit outputs the
expression $\epsilon_{r}(SOC)=6.8+18.2\cdot SOC$.

The final uncertainty intervals are translated through eq. 4 and stem from a
combination of our (arbitrarily chosen) uncertainty in the partial charge
$Z_{Li}$ (as described in the methodology section) and the root mean square
error (RMSE) of the slope, which is returned by the linear regression and
caused by scatter around the Coulombic behaviour, which in turn is due to
variations in the electron densities of adjacent layers, local structural
perturbations and similar effects not captured by the approximations we made
when introducing $E_{Coul}\approx E_{pot}-E_{0}$ (as described throughout the
previous sections), as well as the intrinsic limited accuracy of DFTB.

##  4 CONCLUSION:

With this work, we present the first rigorous investigation of the dielectric
behaviour of lithium graphite intercalation compounds (Li-GICs) for the entire
functional range of charge during the application as anode of a modern Li-ion
battery. In doing so, we provide a straightforward approach for investigating
the intrinsic relative permittivity of materials with mobile charge carriers
in the electrostatic limit, that can be applied to other materials in the
future – given that sufficiently fast and long-ranged computational methods
are available. Thanks to our recently published DFTB parametrization [23, 24],
we are able to sample the long-ranged Coulomb interactions between two
intercalated charge carriers (Li-ions or vacancies) in a variety of
stoichiometries and configurational realizations thereof. By examining the
approximations we made during this process, we additionally outline ways to
further improve this methodology in the future.

The primary finding of this work is the mostly linear dependency (from
$\epsilon_{r}\approx 7$ at SOC 0% to $\epsilon_{r}\approx 25$ at SOC 100%) of
the relative permittivity on the state of charge, which we put forward for the
first time. With this, we make valuable contributions to the future modeling
of functional materials by means of charged kinetic Monte Carlo and continuum
simulations and to the general understanding of Li-GICs. Our results hold for
qualitatively different realizations of the intermediate stoichiometries
(staged or dilute), as well as both possible types of diffusion (Li-ion or
vacancy). The few available experimental studies agree reasonably well with
our results, but more investigation is necessary to really pinpoint the
quantitative dielectric response, especially at higher SOC. Thanks to our
results, it is now clear that only two measurements – one at low and one at
high SOC – would suffice for that purpose.

We find that an approximation neglecting the local distortions of the
structure and of the electron density caused by the sampled charge carriers
holds very well if those charge carriers are Li-ions, and slightly less well,
but still within reason, in the vacancy case. Future improvements could be
achieved by correcting with some local descriptor based machine learning model
that is specifically trained to pick up those local effects.

Additionally, we provide a rough estimate (grain-boundary effects are
neglected, but should be uniform with space direction and therefore should not
have any qualitative impact on our results) on how the relative permittivity
of a perfect crystal and a powder can be compared, or relatedly, of how an
experimental result for (intercalated) graphite powder can be translated to
the internal $xy$-plane relative permittivity the Li-ions feel locally (which
is the one that is actually relevant for their diffusion behaviour).

##  5 Declaration of competing interests

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

##  6 Acknowledgements

This work was funded by the German Federal Ministry of Education and Research
(BMBF) as part of the research cluster ‘‘AQua’’ within the project InOPlaBat
(grant number 03XP0352). The authors gratefully acknowledge the computational
and data resources provided by the Leibniz Supercomputing Centre (LRZ). The
authors jointly thank Sebastian Matera, Jakob Filser, David Egger, Cristina
Grosu and Julian Holland for fruitful discussions.

## References

  * Chung and Xi [2021] D. Chung, X. Xi,  Factors that govern the electric permittivity of carbon materials in the graphite allotrope family,  Carbon 184 (2021) 245–252.
  * Grosu et al. [2021] C. Grosu, C. Panosetti, S. Merz, P. Jakes, S. Matera, R.-A. Eichel, J. Granwehr, C. Scheurer,  Revisiting the storage capacity limit of graphite battery anodes: spontaneous lithium overintercalation at ambient pressure,  arXiv preprint arXiv:2107.11137 (2021).
  * Daumas and Herold [1969] N. Daumas, A. Herold,  Relations between phase concept and reaction mechanics in graphite insertion compounds,  C. R. Acad. Sci. C 268 (1969) 373.
  * Hotta et al. [2011] M. Hotta, M. Hayashi, M. T. Lanagan, D. K. Agrawal, K. Nagata,  Complex permittivity of graphite, carbon black and coal powders in the ranges of X-band frequencies (8.2 to 12.4 GHz) and between 1 and 10 GHz,  ISIJ int. 51 (2011) 1766–1772.
  * Gavilán-Arriazu et al. [2020] E. M. Gavilán-Arriazu, O. A. Pinto, B. L. De Mishima, D. Barraco, O. A. Oviedo, E. P. M. Leiva,  Kinetic Monte Carlo applied to the electrochemical study of the Li-ion graphite system,  Electrochim. Acta 331 (2020) 135439.
  * Gavilán-Arriazu et al. [2021] E. M. Gavilán-Arriazu, M. Mercer, D. E. Barraco, H. Hoster, E. Leiva,  Kinetic Monte Carlo simulations applied to Li-ion and post Li-ion batteries: a key link in the multi-scale chain,  Prog. Energy (2021).
  * Methekar et al. [2011] R. N. Methekar, P. W. Northrop, K. Chen, R. D. Braatz, V. R. Subramanian,  Kinetic Monte Carlo simulation of surface heterogeneity in graphite anodes for lithium-ion batteries: passive layer formation,  J. Electrochem. Soc. 158 (2011) A363.
  * Dean et al. [2021] J. M. Dean, S. W. Coles, W. R. Saunders, A. R. McCluskey, M. J. Wolf, A. B. Walker, B. J. Morgan,  Overscreening and underscreening in solid-electrolyte grain boundary space-charge layers,  Phys. Rev. Lett. 127 (2021) 135502.
  * Casalegno et al. [2010] M. Casalegno, G. Raos, R. Po,  Methodological assessment of kinetic Monte Carlo simulations of organic photovoltaic devices: The treatment of electrostatic interactions,  J.Chem. Phys. 132 (2010) 094705.
  * Hein et al. [2020] S. Hein, T. Danner, A. Latz,  An electrochemical model of lithium plating and stripping in lithium ion batteries,  ACS Appl. Energy Mater. 3 (2020) 8519–8531.
  * Reed et al. [2010] J. P. Reed, B. Uchoa, Y. I. Joe, Y. Gan, D. Casa, E. Fradkin, P. Abbamonte,  The effective fine-structure constant of freestanding graphene measured in graphite,  Science 330 (2010) 805–808.
  * Jellison Jr et al. [2007] G. E. Jellison Jr, J. D. Hunn, H. N. Lee,  Measurement of optical functions of highly oriented pyrolytic graphite in the visible,  Phys. Rev. B 76 (2007) 085125.
  * Elias et al. [2011] D. Elias, R. Gorbachev, A. Mayorov, S. Morozov, A. Zhukov, P. Blake, L. Ponomarenko, I. Grigorieva, K. Novoselov, F. Guinea, et al.,  Dirac cones reshaped by interaction effects in suspended graphene,  Nat. Phys. 7 (2011) 701–704.
  * Bostwick et al. [2010] A. Bostwick, F. Speck, T. Seyller, K. Horn, M. Polini, R. Asgari, A. H. MacDonald, E. Rotenberg,  Observation of plasmarons in quasi-freestanding doped graphene,  Science 328 (2010) 999–1002.
  * Bessler et al. [2019] R. Bessler, U. Duerig, E. Koren,  The dielectric constant of a bilayer graphene interface,  Nanoscale Adv. 1 (2019) 1702–1706.
  * Raabe and Sadus [2011] G. Raabe, R. J. Sadus,  Molecular dynamics simulation of the dielectric constant of water: The effect of bond flexibility,  J. Chem. Phys. 134 (2011) 234501.
  * Aragones et al. [2011] J. Aragones, L. MacDowell, C. Vega,  Dielectric constant of ices and water: a lesson about water interactions,  J. Phys. Chem. A 115 (2011) 5745–5758.
  * Sharma et al. [2007] M. Sharma, R. Resta, R. Car,  Dipolar correlations and the dielectric permittivity of water,  Phys. Rev. Lett. 98 (2007) 247401.
  * Ruiz Pestana et al. [2018] L. Ruiz Pestana, O. Marsalek, T. E. Markland, T. Head-Gordon,  The quest for accurate liquid water properties from first principles,  J. Phys. Chem. Lett. 9 (2018) 5009–5016.
  * Gigli et al. [2021] L. Gigli, M. Veit, M. Kotiuga, G. Pizzi, N. Marzari, M. Ceriotti,  Thermodynamics and dielectric response of $BaTiO\\_3$ by data-driven modeling,  arXiv preprint arXiv:2111.05129 (2021).
  * Elstner et al. [1998] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, G. Seifert,  Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties,  Phys. Rev. B 58 (1998) 7260.
  * Hourahine et al. [2020] B. Hourahine, B. Aradi, V. Blum, F. Bonafé, A. Buccheri, C. Camacho, C. Cevallos, M. Y. Deshaye, T. Dumitrică, A. Dominguez, S. Ehlert, M. Elstner, T. van der Heide, J. Hermann, S. Irle, J. J. Kranz, C. Köhler, T. Kowalczyk, T. Kubař, I. S. Lee, V. Lutsker, R. J. Maurer, S. K. Min, I. Mitchell, C. Negre, T. A. Niehaus, A. M. N. Niklasson, A. J. Page, A. Pecchia, G. Penazzi, M. P. Persson, J. Řezáč, C. G. Sánchez, M. Sternberg, M. Stöhr, F. Stuckenberg, A. Tkatchenko, V. W. z. Yu, T. Frauenheim,  DFTB+, a software package for efficient approximate density functional theory based atomistic simulations,  J. Chem. Phys. 152 (2020) 124101–1–124101–19.
  * Panosetti et al. [2021] C. Panosetti, S. B. Anniés, C. Grosu, S. Seidlmayer, C. Scheurer,  DFTB modeling of lithium-intercalated graphite with machine-learned repulsive potential,  J. Phys. Chem. A. 125 (2021) 691–699.
  * Anniés et al. [2021] S. Anniés, C. Panosetti, M. Voronenko, D. Mauth, C. Rahe, C. Scheurer,  Accessing Structural, Electronic, Transport and Mesoscale Properties of Li-GICs via a Complete DFTB Model with Machine-Learned Repulsion Potential,  Materials 14 (2021) 6633.
  * Panosetti et al. [2020] C. Panosetti, A. Engelmann, L. Nemec, K. Reuter, J. T. Margraf,  Learning to use the force: Fitting repulsive potentials in density-functional tight-binding with gaussian process regression,  J. Chem. Theory Comput. 16 (2020) 2181–2191.
  * Bahn and Jacobsen [2002] S. R. Bahn, K. W. Jacobsen,  An object-oriented scripting interface to a legacy electronic structure code,  Comput. Sci. Eng. 4 (2002) 56–66.
  * Shanno [1985] D. F. Shanno,  On Broyden-Fletcher-Goldfarb-Shanno method,  J. Optim. Theory Appl. 46 (1985) 87–94.
  * Berne et al. [1998] B. J. Berne, G. Ciccotti, D. F. Coker, Classical and quantum dynamics in condensed phase simulations: Proceedings of the International School of Physics, World Sci., 1998.
  * Henkelman and Jónsson [2000] G. Henkelman, H. Jónsson,  Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points,  J. Chem. Phys. 113 (2000) 9978.
  * Henkelman et al. [2000] G. Henkelman, B. P. Uberuaga, H. Jónsson,  A climbing image nudged elastic band method for finding saddle points and minimum energy paths,  J. Chem. Phys. 113 (2000) 9901.
  * Johnson [1988] D. D. Johnson,  Modified Broyden’s method for accelerating convergence in self-consistent calculations,  Phys. Rev. B 38 (1988) 12807.
  * Zhechkov et al. [2005] L. Zhechkov, T. Heine, S. Patchkovskii, G. Seifert, H. A. Duarte,  An efficient a posteriori treatment for dispersion interaction in density-functional-based tight binding,  J. Chem. Theory Comput. 1 (2005) 841.
  * Koskinen and Mäkinen [2009] P. Koskinen, V. Mäkinen,  Density-functional tight-binding for beginners,  Comput. Mater. Sci. 47 (2009) 237–253.
  * Valencia et al. [2006] F. Valencia, A. H. Romero, F. Ancilotto, P. L. Silvestrelli,  Lithium adsorption on graphite from density functional theory calculations,  J. Phys. Chem. B 110 (2006) 14832–14841.
  * Krishnan et al. [2013] S. Krishnan, G. Brenet, E. Machado-Charry, D. Caliste, L. Genovese, T. Deutsch, P. Pochet,  Revisiting the domain model for lithium intercalated graphite,  Appl. Phys. Lett. 103 (2013) 251904.
  * Song et al. [2001] M. K. Song, S. Do Hong, K. T. No,  The structure of lithium intercalated graphite using an effective atomic charge of lithium,  J. Electrochem. Soc. 148 (2001) A1159.
  * Rakotomahevitra et al. [1992] A. Rakotomahevitra, C. Demangeat, J. Parlebas, G. Moraitis, E. Razafindrakoto,  Electronic structure and properties of light atoms intercalated in graphite,  J. Phys. Condens. Matter 4 (1992) 4621.
  * Pande and Viswanathan [2018] V. Pande, V. Viswanathan,  Robust high-fidelity DFT study of the lithium-graphite phase diagram,  Phys. Rev. Mater. 2 (2018) 125401.