# Kondo Effect in Defect-bound Quantum Dots Coupled to NbSe2

T. R. Devidas  Tom Dvir  The Racah Institute of Physics, The Center for
Nanoscience and Nanotechnology, The Hebrew University, Jerusalem 91904, Israel
Enrico Rossi  Department of Physics, William & Mary, Williamsburg, VA 23187,
USA  Hadar Steinberg<EMAIL_ADDRESS>The Racah Institute of Physics,
The Center for Nanoscience and Nanotechnology, The Hebrew University,
Jerusalem 91904, Israel

###### Abstract

> We report the fabrication of a van der Waals tunneling device hosting a
> defect-bound quantum dot coupled to NbSe2. We find that upon application of
> magnetic field, the device exhibits a zero-bias conductance peak. The peak,
> which splits at higher fields, is associated with a Kondo effect. At the
> same time, the junction retains conventional quasiparticle tunneling
> features at finite bias. Such coexistence of a superconducting gap and a
> Kondo effect are unusual, and are explained by noting the two-gap nature of
> the superconducting state of NbSe2, where a magnetic field suppresses the
> low energy gap associated with the Se band. Our data shows that van der
> Waals architectures, and defect-bound dots in them, can serve as a novel and
> effective platform for investigating the interplay of Kondo screening and
> superconducting pairing in unconventional superconductors.

##  I Introduction

The Kondo effect [1] is responsible for the low temperature resistivity upturn
in metals with dilute magnetic impurities [2] , and is a paradigmatic problem
in Condensed Matter Physics. In the Kondo effect, the coupling of a single
spin to a metallic environment causes the formation of a magnetic screening
cloud. It embodies the complexity arising when a single particle interacts
with a many-body environment, and was studied extensively using electronic
transport through GaAs quantum dots (QDs) [3] , and later carbon nanotubes [4]
and single molecules [5] . In these systems, the effect manifests itself as an
enhancement of the electronic conductance due to correlations between the
source and drain mediated by an unpaired spin localized within the QD.

A natural extension of the Kondo problem involves coupling a QD to one or two
superconducting leads. The first experimental study where a QD in the Kondo
regime was coupled to superconducting (S) source and drain electrodes was
reported by Buitelaar et al. [6] . In this ‘S-QD-S’ geometry, the electron
number in the impurity level was tuned to an odd value via electronic gating.
The SU(2) Kondo signature – zero-energy state - was observed whenever the
Kondo temperature $T_{K}$ became the dominant energy scale in the system. Such
an observation falls in line with the theory of many-body singlet state for
the Kondo ground state, which can create a Kondo signature whenever the
binding energy of the Kondo singlet exceeds the superconducting energy gap
$\Delta$ [7] .

Later studies have shown that even when the zero-bias Kondo feature is
suppressed, an applied source-drain bias $V_{SD}=\Delta$ reveals the formation
of a Kondo resonance within the QD - the strength of which appears to follow a
scaling law defined by the competition between $\Delta$ and $k_{B}T_{K}$ [8,
7] , where $k_{B}$ is the Boltzmann constant. This competition underlies a
rich set of physical phenomena involving the Andreev reflection at the QD-S
interface and the Kondo singlet formation at the QD-N (normal) interface [7,
9, 10] . Theoretical studies have suggested that adopting a hybrid N-QD-S
geometry would allow a more rigorous approach in studying these dynamics [7,
11, 12, 13] .

Experimental realizations of N-QD-S devices typically rely on clean
semiconductor nanowires (GaAs, InAs) [14] , coupled to superconductor
materials that could be evaporated as electrodes (Pb, Nb, Al, MoRe). The
advent of van der Waals (vdW) materials now enables the study of a whole
library of materials with diverse ground states [15] . Specifically, the
Transition Metal Dichalcogenides (TMD) family contains candidate 2D
superconductors (NbSe2, NbS2, PdTe2) that exhibit superconductivity from the
bulk regime to the monolayer limit [16] . Insulators such as MoS2, WSe2, and
hexagonal Boron Nitride (hBN) have been successfully implemented as barriers
in tunnel junctions [17, 18, 19] and have been shown to host atomic defects
[20, 21, 22] which could be used as QDs [23, 24, 25] .

In this work we demonstrate that a Kondo effect can also be realized by
coupling a QD to a layered superconductor. To place a QD in close proximity to
a layered superconductor, we make use of naturally-occurring defects in a TMD
semiconductor tunnel barrier placed on top of the superconductor NbSe2. As we
have shown in previous studies, such defect-bound dots may couple strongly to
the SC, giving rise to Andreev bound state (ABS) sub-gap conductance features
[21] . Conversely, when the QD is weakly coupled, it can serve as a sensitive
spectral [22] or compressibility [25] probe. In the present case, we find a
tell-tale zero-bias conductance peak which sets in at finite in-plane and out-
of-plane magnetic field. The peak splits at higher magnetic fields, consistent
with a Kondo feature. In NbSe2, tunneling measurements resolve two gaps [26,
17] . Here we find that the Kondo feature is correlated with a suppression in
the spectral signature of the lower of these two energy gaps, associated with
the Se band. The correlation of the peak value with the suppression of the
second gap tunneling feature suggests that the Kondo effect is mediated by
carriers belonging to the Se-derived band, which turns normal at a low
magnetic field, while the larger Nb-derived gap remains stable. Our results
suggest that defect-bound QDs can be used to study the Kondo effect in
unconventional, van der Waals, layered superconductors.

##  II Experimental Details

The schematic of the MoS2-graphene-NbSe2 vdW heterostructure is shown in
Figure 1(a) and the optical image of the measured device reported in this work
is shown in Figure 1(b). MoS2, graphene and NbSe2 are exfoliated independently
on Si/SiO2 (285 nm oxide) substrates. The desired flakes of 3-4 layers MoS2
(tunnel barrier), monolayer graphene and bulk NbSe2 are chosen by optical
contrast. The MoS2 barrier flake is first picked up using the polycarbonate
(PC) technique. Graphene is subsequently picked up using the van der Waals
interaction between MoS2 and graphene. The picked up heterostructure is then
transferred on to a 25 nm thick bulk NbSe2 flake. The process is carried out
in a glovebox under an Argon atmosphere. Standard e-beam lithographic
techniques are used to pattern tunnel electrodes and ohmic contacts on MoS2
and NbSe2 respectively. Ti/Au electrodes are evaporated using an e-beam
evaporator with an additional Argon ion milling step prior to the ohmic
contact deposition step, so as to obtain better contact on the NbSe2 flake. No
contacts are made to the graphene flake.

##  III Results

Figure 1: (a) Schematic of the device measured in the current work; (b)
Optical image of the measured device. NbSe2, graphene and MoS2 flakes are
outlined and labelled in different colours for easier identification. S and D
indicate the source and drain electrodes respectively. (c) Left $y$-axis -
Differential conductance ($dI/dV$) spectrum of the NbSe2-graphene-MoS2 tunnel
junction depicted in (b), Right $y$-axis - Derivative of the differential
conductance ($d^{2}I/dV^{2}$) spectrum measured at 28 mK. The black dashed
lines track the two NbSe2 gaps present in the spectrum at 0.6 mV and 1.25 mV;
(d) and (e) show the evolution of the tunneling spectrum as a function of
parallel ($B_{||}$) and perpendicular magnetic field ($B_{\perp}$)
respectively, for low magnetic field values.

The tunneling differential conductance ($dI/dV$) of the MoS2-graphene-NbSe2
stack at 28 mK and B = 0 T, is shown on the left $y$-axis of Figure 1(c). The
spectrum exhibits well-defined quasiparticle peaks and a hard gap - manifest
in a ratio of over 100 between the conductance outside the gap ($G_{N}$) and
zero-bias conductance ($G_{0}$). Such a hard gap, seen in our earlier studies
on NbSe2 tunneling devices [17, 18] , attests to the quality of the tunnel
junction, which suppresses two-particle tunneling contributions. The $2^{nd}$
derivative of the tunneling current ($d^{2}I/dV^{2}$), plotted in the right
$y$-axis of Figure 1(c) shows a clear separation between two distinct features
marked by black dashed lines. These two features correspond to the higher
energy gap, $\Delta_{1}=1.25$ mV, and the lower energy gap $\Delta_{2}=0.6$
mV. These were identified in the past as related to the two superconducting
bands in NbSe2 – the 2D-like niobium $d$-orbitals bands around the $\Gamma,K$
points, associated with $\Delta_{1}$, and selenium-derived $p$-orbitals which
form a small 3D Fermi surface around the $\Gamma$-point, associated with
$\Delta_{2}$ [17, 18, 27, 26] .

Interestingly, the tunneling spectrum does not reveal any signature which can
be associated with the graphene layer. This can be understood considering that
only for a small range of twist angles graphene’s Fermi pockets can have some
overlap with the Fermi pockets of a TMD like NbSe2 [28] . As a result, the
graphene spectrum may be decoupled from the NbSe2.

Figure 1(d) shows the evolution of the tunneling spectra as a function of in-
plane magnetic field ($B_{\parallel}$) applied parallel to the sample plane.
The observed spectra evolve in a way which is not characteristic of NbSe2
devices measured in the past [17] : First, we find that the superconducting
gap is filled even by the application of a very low field. Second, a zero-bias
conductance peak (ZBCP) develops at $B_{\parallel}=0.275$ T. In what follows,
we discuss possible mechanisms for these features.

Zero-bias spectral features are rather common in proximity superconducting
devices, and may have a number of possible origins. They could be associated
with Andreev bound states (ABS) residing on the surface of a $d$-wave
superconductor [29] , or with constructive superposition of bound states
formed by reflectionless tunneling at diffusive N-S junctions [30, 31] .

Zero-bias states may appear when a QD, proximity-coupled to a superconductor,
forms a local Andreev bound state which undergoes a singlet-doublet transition
[32, 33] . Conversely, when the QD is weakly coupled to both SC
($\Gamma_{SC}$) and normal ($\Gamma_{N}$) leads, i.e. $\Gamma_{SC,N}<<\Delta$,
QD-SC transport is dominated by single electron resonant tunneling. In this
regime, the spectrum exhibits sharp conductance peaks, but these can appear
both above the gap and below it, when the density of states (DOS) is not zero.
In a recent publication [22] we have reported the use of such a weakly-coupled
QD as a sensitive spectrometer.

The zero-bias conductance peak also appears as a response to out-of-plane
magnetic field $B_{\perp}$ as seen in Figure 1(e). To further study the nature
of this peak, we track it’s evolution as magnetic fields are applied in
directions parallel and perpendicular to the sample plane. Figure 2(a)
represents a color scale map of the tunneling spectra vs. $B_{\parallel}$
extending to 9 T. Following the evolution of the zero-energy state that
switches on at 0.275 T (Figure 1(d)), the peak persists at the same energy
until a magnetic field induced splitting of the state into two distinct
features is observed at 1 T. The splitting is symmetrical in $V_{SD}$ with
respect to zero.

Although such linear dispersions are regularly observed in ABS, here we argue
that these features are unlikely to be related with ABS, on two grounds.
First, ABS features in NbSe2 usually appear at finite bias, and may converge
to zero at finite field. They are observed starting at zero magnetic field and
do not require a finite field to be visible. Second, ABS features in NbSe2
tunnel devices are observed only in very thin flakes. In any flake thicker
than a few layers, the in-plane magnetic field introduces a sub-gap tunneling
signal which obscures them [21] .

We are thus led to propose that the observed feature - the zero-bias
conductance peak at a finite field, is associated with a Kondo origin.

The split states are separated by 2$E_{Z}$ [3, 34] , $E_{Z}=\pm g\mu_{B}B$
being the Zeeman energy, with a Landé $g$-factor (1.67$\pm$0.04). A similar
trend is observed in the evolution of the spectrum in $B_{\perp}$ with a field
induced splitting appearing at 0.6 T, shown in Figure 2(b). The value of Landé
$g$-factor obtained from out-of-plane data is 1.94$\pm$0.04. We note that
similar values have been seen in the past for atomic defect QDs [21, 22] .

Figure 2: Color map of the $dI/dV$ spectrum in (a) in-plane magnetic field
($B_{||}$) and (b) out-of-plane magnetic field ($B_{\perp}$). The black arrows
in (a) and (b) indicate the Zeeman split features.

The Kondo effect sets in below a Kondo temperature $T_{K}$ which can be
obtained through temperature dependence or through the profile of the zero-
bias conductance peak. Temperature dependence measurements of the zero-bias
conductance peak ($B_{\perp}$ = 25 mT) are carried out from 28 mK to 419 mK,
beyond which the peak is no longer detectable above the background. A
selection of these plots is shown in Figure 3(a). The maximal zero-bias
conductance at 28 mK is 11.1 $\mu S$ = 0.14 $e^{2}/h$. Figure 3(b) shows a
semi-log plot of zero-bias conductance peak heights as a function of
temperature, where peak height values are obtained after suitable background
subtraction. The data is fit to an empirical equation (1) derived from the
numerical renormalization group (NRG) theory for the Kondo ground state [6, 3]
.

$\centering
G(T)=\frac{G_{max}}{\bigg{[}1+\bigg{(}2^{\nicefrac{{1}}{{s}}}-1\bigg{)}\bigg{(}\frac{T}{T_{K}}\bigg{)}^{2}\bigg{]}^{s}}\@add@centering$
(1)

Where $G_{max}$ is the maximum conductance at the lowest temperature measured,
$T_{K}$ is the Kondo temperature and $s$ is a dimensionless value related to
the spin-state of the electron in the QD. The best fit to the data (solid red
line) yields $G_{max}$ = 0.0695$\pm$0.003 $e^{2}/h$; $T_{K}$ = 153$\pm$6 mK
(13.2 $\mu$eV) and $s$ = 1.4$\pm$0.3.

Clearly, the exponent $s$ deviates from the value $s$ = 0.22, expected for a
spin 1/2 system at the Kondo regime [3, 6] . Values removed from 0.22 could
indicate that the QD is in the mixed-valence regime, which is seen when it’s
energy is close to resonance, or, alternatively, that the spin is not 1/2.

Figure 3: Evolution of the sub-gap conductance ($dI/dV$) spectra as a function
of temperature in the QD-NbSe2 junction (a) Scans at select temperatures; (b)
Kondo peak heights as a function of temperature, after suitable background
subtraction, fitted (red solid line) to the NRG derived empirical equation [6,
3] ; (c) The zero-bias conductance peak at 28 mK, fit to the Fano function
[35] to extract the value of HWHM which defines $T_{K}$ [36] ; (d) The zero-
bias conductance peak at 28 mK, fit to the equation derived by Kretinin et al.
[37] to different values of $T_{K}$ and $s$.

It is also possible to evaluate $T_{K}$ from the conductance feature line-
shape. In a N-QD-S system, the half-width at half maximum (HWHM) of the peak
at the lowest temperature is directly related to $T_{K}$ [36] . Since the
Kondo state involves an interaction between a discrete state (QD) and
continuum (N lead), we use the Fano function [35] to extract the HWHM (Figure
3(c)). The function is defined as

$\centering
G(V_{SD})=A\frac{(\varepsilon+q)^{2}}{(1+\varepsilon^{2})}+B;\varepsilon=\frac{(V_{SD}-\varepsilon_{0})}{\Gamma_{d}}\@add@centering$
(2)

Where $\varepsilon_{0}$ is the resonance energy (0 eV in the present data), A
and B are constants, $\Gamma_{d}$ is the HWHM, and $q$ is a phenomenological
dimensionless shape parameter. The HWHM obtained is $\Gamma_{d}$ = 25.8$\pm$6
$\mu$eV which corresponds to $T_{K}$ = 300$\pm$70 mK. The other parameters
obtained from the fit are A = -5.5$\pm$0.7 $\mu$S, B = 0.11$\pm$0.06 $\mu$S, q
= 0.10$\pm$0.05.

The ambiguity in the extracted Kondo temperature has been observed earlier on
N-QD-N devices too [37, 38] . Kretinin et al. [37] improved on the earlier NRG
formula by suggesting that analysing the zero energy state at the lowest
temperature as a function of a normalised energy scale $\nu\equiv
eV_{SD}/(k_{B}T_{K})$ is a more reliable method to extract the correct Kondo
temperature of the system. It was based on the argument that temperature
dependence data might have additional features arising from non-Kondo origin
incorporated in the zero-energy state as a result of which the extracted
$T_{K}$ would not define the Kondo state alone. The argument might hold for
the anomaly in the value of $s$ parameter we obtained using the empirical
formula fit.

Using the equation derived by Kretinin et al.:

$\centering
G(T=0,\nu)=\frac{G_{max}}{\bigg{[}1+\bigg{(}2^{\nicefrac{{1}}{{s}}}-1\bigg{)}\frac{\nu^{2}}{\pi}\bigg{]}^{s}}\@add@centering$
(3)

and assuming our base temperature is close to the condition T = 0, we plot our
28 mK data (Figure 3(d)) along with two curves generated using Equation (3).
We notice that the curve generated by the $T_{K}$ and $s$ values obtained from
the temperature dependence fit (blue dash dot) doesn’t replicate the behaviour
of the measured data. However, for the values $T_{K}$ = 270 mK (23.26 $\mu$eV)
and $s$ = 0.22 (red dash dot), the curve closely follows the measured data.

##  IV Discussion

The results presented so far suggest that the zero-bias conductance peak we
observe is indeed associated with a Kondo feature. Yet the larger gap
$\Delta_{1}$ = 1.25 mV is far greater than $k_{B}T_{K}$, regardless of the
method we choose to evaluate $T_{K}$. Although theoretical models do suggest
that such a co-existence is possible, with a zero-bias feature retained even
when $k_{B}T_{K}\ll\Delta$ [7, 11, 39, 40] , most experiments carried out so
far on N-QD-S and S-QD-S systems by various groups [6, 8, 9, 10, 41, 42, 36,
43] do not observe Kondo features at this limit.

Here we suggest that the co-existence between the Kondo effect and
superconductivity is a consequence of the 2-gap nature of NbSe2. Since we
observe the Kondo zero-energy state even at fields far from critical fields of
bulk 2H-NbSe2, $H_{c2}^{\perp}\approx$ 4 T, $H_{c2}^{\parallel}\approx$ 17 T
[44] , the normal electrons must be associated with a superconducting band
that turns normal even at such low fields. As discussed previously, bulk NbSe2
has two superconducting gaps - one associated with niobium based bands around
the $\Gamma$ and $K$ points and a second proximitized inner gap from the small
selenium based bands around the $\Gamma$ point [17, 27] whose signatures are
observed at 1.25 mV and 0.6 mV respectively.

To elucidate the role of the Se-derived band in mediating the Kondo effect, we
examine the correlation between the appearance of a zero-energy state and the
second band in out-of-plane magnetic field (Figure 4). The out-of-plane
magnetic field is chosen since the inner gap exhibits stronger response to its
onset [18] . Figure 4(a) shows the zero-bias conductance as a function of
$B_{\perp}$ as it is swept in the positive direction (indicated by the red
arrow). A complete dataset tracing the zero-bias conductance $G_{0}$ from
negative to positive fields is shown in Supplementary Figure 1. The
conductance reaches the lowest value of 0.16 $\mu$S at $B_{\perp}$ = -0.2 mT,
where the tunneling spectrum reaches a sub-gap conductance typical of a hard-
gap NbSe2 junction [17] . At this limit, there is no indication of a zero-bias
conductance peak.

As $B_{\perp}$ is increased, $G_{0}$ increases in discrete steps of $\approx$
1 $\mu$S, observed at 1.0 mT, 2.2 mT and 3.1 mT. As $B_{\perp}$ is further
increased, the $G_{0}$ increase becomes continuous. Tunneling spectra are
measured for fields that mark a discrete increase in $G_{0}$ (Figure 4(b)) and
also at two higher fields. The second derivative of tunneling current as a
function of magnetic field is plotted in Figure 4(c). The data in Figure
4(b)-(c) are offset along the $y$-axis for clarity. We observe that the first
signature of a distinct zero-bias peak appears at a field of 3.1 mT and
persists at higher fields. Simultaneously, the feature at 0.6 mV corresponding
to the Se derived superconducting bands loses prominence until it is difficult
to resolve from its background.

Figure 4: (a) Differential conductance at zero-bias ($G_{0}$) as a function of
$B_{\perp}$ at mT resolution. The red arrow indicates the sweep direction of
the $B_{\perp}$; (b) sub-gap tunneling spectra at various field values
indicated by the coloured dotted lines in (a); (c) second derivative of the
tunneling current at field values indicated by dotted lines in (a); (d) zero-
bias conductance peak height (left $y$-axis) and height of the feature at 0.6
mV in second derivative of tunneling current (right $y$-axis) as a function of
$B_{\perp}$

Figure 4(d) shows the evolution of the zero-bias conductance peak height (left
$y$-axis) and the height of the feature in the second derivative of tunneling
current (right $y$-axis) corresponding to the inner Se superconducting band at
0.6 mV. An inverse correlation is evident from the data: The peak at 0.6 mV is
the strongest at -0.2 mT where the conductance at zero-bias reaches the lowest
value. With reduction in height of the feature in magnetic field, the zero-
bias conductance feature grows in strength. The reduced height of the zero-
bias conductance peak at higher fields can be attributed to the higher
background conductance with magnetic field. Supplementary Figure 2 shows the
individual sub-gap tunneling spectra shown in Figure 4(b) with their
background conductance. The discrete jumps ($\approx 1\mu$S) in the zero-bias
conductance with magnetic field can be associated with the discrete entry of
vortices in the tunnel junction area of NbSe2 [18] . This correlation thus
supports the interpretation that the Se band turning normal could assist in
the spin-flip co-tunneling event, as depicted in the schematic Figure 5(a).
Interestingly, we have shown in a previous work that in the weak coupling
limit, a defect-bound QD favors coupling to the Se-band over coupling to the
Nb-bands [22] .

Figure 5: (a) Schematic of spin-flip co-tunneling event – (above) prohibited
at zero applied magnetic field due to absence of states within the
superconducting gap, (below) allowed via filled second band in NbSe2 in an
applied magnetic field; (b) Sub-gap density of states of NbSe2 as a function
of $B_{||}$, obtained via a large N expansion with a broadening $\delta$ =
0.04 meV. The inset shows the experimental differential tunneling conductance
data for the same values of $B_{||}$.

To better understand the origin of the observed zero-bias conductance peak we
study a Kondo model in which we treat the weakly coupled QD as an impurity,
with magnetic moment ${\bf S}$, placed at position ${\bf R}$. The Hamiltonian
for the electrons in NbSe2 can be written as $H=H_{0}+H_{J}$, where $H_{0}$ is
the Hamiltonian when no QD is present, and $H_{J}$ describes the coupling of
the QD to the electrons in NbSe2. In the Kondo limit
$H_{J}=J\sum_{\sigma\sigma^{\prime}}{c}^{\dagger}_{{\bf
R}\sigma}{\bm{\tau}}_{\sigma\sigma^{\prime}}{c}_{{\bf
R}\sigma^{\prime}}\cdot{\bf S}$, where $J>0$ is the antiferromagnetic coupling
between the QD and NbSe2 electronic states, $\tau_{i}$ are the Pauli matrices
in spin space, ${c}^{\dagger}_{{\bf R}\sigma}$ and ${c}_{{\bf R}\sigma}$ are
the creation and annihilation operators, respectively, for an electron at
position ${\bf R}$ with spin $\sigma$. To treat the interaction between ${\bf
S}$ and the electrons we use the large-N expansion [45, 46] . The large-N
expansion has been shown to give accurate results, in agreement with other
methods such as the numerical renormalization group [47, 48] , when the
impurity is isolated and the Kondo problem has only one channel. The
electrons’ density of states $\rho(\varepsilon)$ is the key property that
determines the features of the Kondo effect. Here, instead of assuming a
constant $\rho(\varepsilon)$, or a prescribed form [49, 50, 51] , we use the
density of states extracted for different values of the magnetic field from
spectroscopic measurements on few-layer superconducting NbSe2 [17] . One
important aspect of using the experimentally obtained $\rho(\varepsilon)$ is
that it allows us to correctly capture the effect of $B$ on it, and therefore
the unusual evolution of the Kondo peak with $B$. In the large-N expansion
${\bf S}$ is expressed in terms of auxiliary creation (annihilation) fermionic
operators ${f}^{\dagger}_{\sigma}$ (${f}_{\sigma}$) satisfying the constraint
$n_{f}=\sum_{\sigma}{f}^{\dagger}_{\sigma}{f}_{\sigma}=1$, with
$\sigma=1,\ldots,N$, so that $H_{J}=J\sum_{{\bf k},{\bf
k}^{\prime},\sigma,\sigma^{\prime}}{c}^{\dagger}_{{\bf k}\sigma}{c}_{{\bf
k}^{\prime}\sigma^{\prime}}{f}^{\dagger}_{\sigma^{\prime}}{f}_{\sigma}$, with
${c}^{\dagger}_{{\bf k}\sigma}$ (${c}_{{\bf k}\sigma}$) the creation
(annihilation) operators for an electron with momentum ${\bf k}$ and spin
$\sigma$. In the remainder we assume $|{\bf S}|=1/2$ and therefore set $N=2$.
We decouple the quartic interaction term $H_{J}$ via the mean-field
$m\sim\sum_{{\bf k},\sigma}\langle{\hat{f}}^{\dagger}_{\sigma}{\hat{c}}_{{\bf
k}\sigma}\rangle$.

The constrain $|{\bf S}|=1/2$ implies $n_{f}=1$ and this is enforced via a
Lagrange multiplier, $\mu_{f}$, which plays the role of the chemical potential
of the $f$-electrons. By minimizing the effective action corresponding to
$H=H_{0}+H_{J}$, within the saddle-point approximation, we obtain [46]

$\displaystyle\sum_{\sigma=\pm 1}\int_{-D}^{D}d\varepsilon
n_{F}(\varepsilon)\frac{\rho(\varepsilon)(\varepsilon-\mu_{f}+h\sigma)}{(\varepsilon-\mu_{f}+h\sigma)^{2}+[\pi\rho(\varepsilon)m^{2}]^{2}}+\frac{2}{J}=0;$
(4) $\displaystyle\sum_{\sigma=\pm 1}\int_{-D}^{D}d\varepsilon
n_{F}(\varepsilon)\frac{\rho(\varepsilon)m^{2}}{(\varepsilon-\mu_{f}+h\sigma)^{2}+[\pi\rho(\varepsilon)m^{2}]^{2}}-1=0$
(5)

where $n_{F}(\varepsilon)$ is the Fermi function, $D$ is the electrons’
bandwidth, and $h\approx(1/2)g_{f}\mu_{B}B$, describes the Zeeman effect, with
$g_{f}$ the effective g-factor for the QD extracted from the experimentally
observed Zeeman splitting. Equations (4), (5) can be inverted to find $m^{2}$
and $\mu_{f}$ for a given value of $J$ and $B$.

From Eqs. (4), (5) we obtain the values of $m$ and $\mu_{f}$ in the limit of
$T\to 0$ (we set $T=6\times 10^{-4}$ K). For $\varepsilon\gg\Delta_{1}$,
$\rho(\varepsilon)$, to very good approximation, can be taken to be constant.
Let $\rho_{0}$ be the value of $\rho(\varepsilon)$ for
$\varepsilon\gg\Delta_{1}$. The value of $\rho_{0}$ that enters Eqs. (4), (5)
depends on the effective spatial extension of the magnetic impurity, i.e., in
our case, of the defect-bound QD. $\rho_{0}$ and $J$ are very difficult to
estimate, both theoretically and experimentally. On the other hand, $m$ and
$\mu_{f}$, depend mostly on the product $\rho_{0}J$, and not on the separate
values of $\rho_{0}$ and $J$. We studied the dependence of the Kondo peak on
the product of $\rho_{0}J$ and found the general result that the peak grows
with magnetic field for $\rho_{0}J\lesssim 0.2$ and decreases for
$\rho_{0}J\gtrsim 0.2$. We find that the scaling of the Kondo peak with $B$
observed experimentally is best approximated when $\rho_{0}J\approx 0.14$. In
the remainder we therefore set $\rho_{0}J=0.14$.

We obtain the electrons’ Green’s function at position ${\bf R}$ and energy
$\omega$, $G({\bf R},{\bf R},\omega)$, renormalized by the coupling to the
effective magnetic impurity, via the equations

$\displaystyle G({\bf R},{\bf R},\omega)=G_{0}({\bf R},{\bf
R},\omega)+m^{2}G_{0}({\bf R},{\bf R},\omega)F(\omega)G_{0}({\bf R},{\bf
R},\omega),$ (6) $\displaystyle
F(\omega)=[(\omega+i\delta-\mu_{f}-m^{2}G_{0}({\bf R},{\bf
R},\omega))\tau_{0}-(1/2)g\mu_{B}B\tau_{z}]^{-1},$ (7)

where $G_{0}({\bf R},{\bf R},\omega)=-i(\pi/2)\rho(\omega)\tau_{0}$, is the
electrons’ bare local Green’s function, $F(\omega)$ is the Green’s function
for the effective fermionic degrees of freedom $f$, and $\delta$ is a small
broadening that we introduce to take into account the effect on the spectra of
thermal fluctuations and non-magnetic disorder. The local DOS at the position
of the QD, at energy $\omega$, $\rho_{\rm tot}(\omega)$ is then obtained as
$\rho_{\rm tot}(\omega)=-{\imaginary}[{\rm Tr}(G)]/\pi$.

Figure 5 (b) shows $\rho_{\rm tot}(\omega)$ for different values of $B$
obtained via the large-N expansion assuming $\delta=0.04~{}\text{meV}$. We see
that for $B=0$, given the vanishing of the DOS as $\varepsilon\to 0$ for
superconducting NbSe2, no peak for $\varepsilon\approx 0$ is present in
$\rho_{\rm tot}(\omega)$, indicating the absence of Kondo screening. NRG
results [52] show that when $N(\varepsilon)\sim|\varepsilon|^{a}$ with a $>$
1/2, and perfect particle-hole symmetry is present, no Kondo effect can be
realized. In superconducting NbSe2 the presence of a magnetic field $B$
induces a finite density of states at low energies. This makes possible the
establishment of a Kondo cloud, and, when $B$ is not too large, dominates over
the suppression of $T_{K}$ due to the polarization of ${\bf S}$ induced by the
magnetic field. For $B\gtrsim 0.25$ T a peak develops: the softening, at low
energies, of the superconducting gap allows the establishment of Kondo
screening. At larger $B$ the Kondo peak splits, as expected [53] .

##  V Summary

Our work shows that defect-bound QDs can be used as platforms to probe the
interplay between the energy scales associated with the Kondo effect and with
superconductivity. This opens a number of interesting future possibilities.
First, as discussed above, NbSe2 is a 2-gap superconductor - offering a richer
phase space involving both order parameters. Interestingly, NbSe2 is expected
to develop a triplet order parameter at high $B_{\parallel}$ [54] . This could
give rise to a field-driven transition in the type of Kondo effect. Finally,
it is interesting to consider coupling QDs to unconventional SCs such as
FeTe0.55Se0.45 [55] , which can be produced using exfoliation, or even to the
superconducting states appearing in twisted graphene systems [56] .

##  VI Methods

The tunneling measurements are performed using standard low frequency lockin
techniques in a BlueFors dilution cryostat with a base temperature of 20 mK.
An AC excitation of 30$\mu V$ is applied across the device between the
tunneling electrodes and the ohmic drain electrode using a Zurich instruments
MFLI digital lock-in amplifier (LIA). The superconductor is tuned in and out
of the superconducting gap by applying a DC bias voltage ($V_{SD}$), also
obtained from the LIA. The tunneling current reaching the drain is fed into a
FEMTO current amplifier. The output of the current amplifier is then fed into
a Keithley 2000 digital multimeter to measure $I_{DC}$ and into the LIA input
to measure the differential conductance $dI/dV$ as a function of $V_{SD}$. The
magnetic field measurements are performed with the 2 axis 9T-3T (Z-Y)
superconducting magnet attached to the dilution. The magnet is powered by the
commercial power supply from American Magnetic Inc. (AMI), with a resolution
of 1 mT (supplied). For resolutions smaller than that, the magnet is powered
by a Keithley 2400 SMU in the current source mode. The applied field is
evaluated from the current to field ratio provided by AMI.

##  VII Acknowledgements

Devices for this study were fabricated at the Center for Nanoscience and
Nanotechnology, The Hebrew University. Funding for this work was provided by a
European Research Council Starting Grant (No. 637298, TUNNEL), Israel Science
Foundation grant 861/19, and BSF grant 2016320. E.R. thanks the Aspen Center
for Physics, which is supported by National Science Foundation grant
PHY-1607611, where part of this work was conducted, and Kevin Ingersent for
very helpful discussions.

## References

  * Kondo [1964] J. Kondo, Prog. Theor. Phys. 32, 37 (1964).
  * de Haas  _et al._ [1934] W. J. de Haas, J. de Boer, and G. J. van dën Berg, Physica 1, 1115 (1934).
  * Goldhaber-Gordon  _et al._ [1998] D. Goldhaber-Gordon, J. Göres, M. A. Kastner, H. Shtrikman, D. Mahalu, and U. Meirav, Phys. Rev. Lett. 81, 5225 (1998).
  * Nygård  _et al._ [2000] J. Nygård, D. H. Cobden, and P. E. Lindelof, Nature 408, 342 (2000).
  * Park  _et al._ [2002] J. Park, A. N. Pasupathy, J. I. Goldsmith, C. Chang, Y. Yalsh, J. R. Petta, M. Rinkoski, J. P. Sethna, H. D. Abruña, P. L. McEuen, and D. C. Ralph, Nature 417, 722 (2002).
  * Buitelaar  _et al._ [2002] M. R. Buitelaar, T. Nussbaumer, and C. Schönenberger, Phys. Rev. Lett. 89, 256801 (2002).
  * Clerk  _et al._ [2000] A. A. Clerk, V. Ambegaokar, and S. Hershfield, Phys. Rev. B - Condens. Matter Mater. Phys. 61, 3555 (2000).
  * Buizert  _et al._ [2007] C. Buizert, A. Oiwa, K. Shibata, K. Hirakawa, and S. Tarucha, Phys. Rev. Lett. 99, 136806 (2007).
  * Yeyati  _et al._ [2003] A. L. Yeyati, A. Martín-Rodero, and E. Vecino, Phys. Rev. Lett. 91, 266802 (2003).
  * Sand-Jespersen  _et al._ [2007] T. Sand-Jespersen, J. Paaske, B. M. Andersen, K. Grove-Rasmussen, H. I. Jørgensen, M. Aagesen, C. B. Sørensen, P. E. Lindelof, K. Flensberg, and J. Nygård, Phys. Rev. Lett. 99, 126603 (2007).
  * Cuevas  _et al._ [2001] J. C. Cuevas, A. Levy Yeyati, and A. Martín-Rodero, Phys. Rev. B 63, 094515 (2001).
  * Li  _et al._ [2016] L. Li, Z. Cao, T. F. Fang, H. G. Luo, and W. Q. Chen, Phys. Rev. B 94, 165144 (2016).
  * Žitko  _et al._ [2015] R. Žitko, J. S. Lim, R. López, and R. Aguado, Phys. Rev. B 91, 045441 (2015).
  * Jespersen  _et al._ [2006] T. S. Jespersen, M. Aagesen, C. Sørensen, P. E. Lindelof, and J. Nygård, Phys. Rev. B 74, 233304 (2006).
  * Ajayan  _et al._ [2016] P. Ajayan, P. Kim, and K. Banerjee, Phys. Today 69, 38 (2016).
  * Ugeda  _et al._ [2016] M. M. Ugeda, A. J. Bradley, Y. Zhang, S. Onishi, Y. Chen, W. Ruan, C. Ojeda-Aristizabal, H. Ryu, M. T. Edmonds, H. Z. Tsai, A. Riss, S. K. Mo, D. Lee, A. Zettl, Z. Hussain, Z. X. Shen, and M. F. Crommie, Nat. Phys. 12, 92 (2016).
  * Dvir  _et al._ [2018a] T. Dvir, F. Massee, L. Attias, M. Khodas, M. Aprili, C. H. L. Quay, and H. Steinberg, Nat. Commun. 9, 598 (2018a).
  * Dvir  _et al._ [2018b] T. Dvir, M. Aprili, C. H. L. Quay, and H. Steinberg, Nano Lett. 18, 7845 (2018b).
  * Khestanova  _et al._ [2018] E. Khestanova, J. Birkbeck, M. Zhu, Y. Cao, G. L. Yu, D. Ghazaryan, J. Yin, H. Berger, L. Forró, T. Taniguchi, K. Watanabe, R. V. Gorbachev, A. Mishchenko, A. K. Geim, and I. V. Grigorieva, Nano Letters 18, 2623 (2018).
  * Zhou  _et al._ [2013] W. Zhou, X. Zou, S. Najmaei, Z. Liu, Y. Shi, J. Kong, J. Lou, P. M. Ajayan, B. I. Yakobson, and J.-C. Idrobo, Nano Lett. 13, 2615 (2013).
  * Dvir  _et al._ [2019] T. Dvir, M. Aprili, C. H. L. Quay, and H. Steinberg, Phys. Rev. Lett. 123, 217003 (2019).
  * Devidas  _et al._ [2021] T. R. Devidas, I. Keren, and H. Steinberg, Nano Letters 21, 6931 (2021).
  * Chandni  _et al._ [2015] U. Chandni, K. Watanabe, T. Taniguchi, and J. P. Eisenstein, Nano Lett. 15, 7329 (2015).
  * Greenaway  _et al._ [2018] M. T. Greenaway, E. E. Vdovin, D. Ghazaryan, A. Misra, A. Mishchenko, Y. Cao, Z. Wang, J. R. Wallbank, M. Holwill, Y. N. Khanin, S. V. Morozov, K. Watanabe, T. Taniguchi, O. Makarovsky, T. M. Fromhold, A. Patanè, A. K. Geim, V. I. Fal’ko, K. S. Novoselov, and L. Eaves, Commun. Phys. 1, 94 (2018).
  * Keren  _et al._ [2020] I. Keren, T. Dvir, A. Zalic, A. Iluz, D. LeBoeuf, K. Watanabe, T. Taniguchi, and H. Steinberg, Nat. Commun. 11, 3408 (2020).
  * Noat  _et al._ [2015] Y. Noat, J. A. Silva-Guillén, T. Cren, V. Cherkez, C. Brun, S. Pons, F. Debontridder, D. Roditchev, W. Sacks, L. Cario, P. Ordejón, A. García, and E. Canadell, Phys. Rev. B 92, 134510 (2015).
  * Kiss  _et al._ [2007] T. Kiss, T. Yokoya, A. Chainani, S. Shin, T. Hanaguri, M. Nohara, and H. Takagi, Nat. Phys. 3, 720 (2007).
  * Gani  _et al._ [2019] Y. S. Gani, H. Steinberg, and E. Rossi, Phys. Rev. B 99, 235404 (2019).
  * Millo and Koren [2018] O. Millo and G. Koren, Phil. Trans. R. Soc. A 376, 20140143 (2018).
  * Van Wees  _et al._ [1992] B. J. Van Wees, P. De Vries, P. Magnée, and T. M. Klapwijk, Phys. Rev. Lett. 69, 510 (1992).
  * Marmorkos  _et al._ [1993] I. K. Marmorkos, C. W. J. Beenakker, and R. A. Jalabert, Phys. Rev. B 48, 2811 (1993).
  * Lee  _et al._ [2014] E. J. H. Lee, X. Jiang, M. Houzet, R. Aguado, C. M. Lieber, and S. De Franceschi, Nat. Nanotechnol. 9, 79 (2014).
  * Scherübl  _et al._ [2020] Z. Scherübl, G. Fülöp, C. P. Moca, J. Gramich, A. Baumgartner, P. Makk, T. Elalaily, C. Schönenberger, J. Nygård, G. Zaránd, and S. Csonka, Nat. Commun. 11, 1834 (2020).
  * Pustilnik and Glazman [2004] M. Pustilnik and L. Glazman, J. Phys.: Condens. Matter 16, 513 (2004).
  * Fano [1961] U. Fano, Phys. Rev. 124, 1866 (1961).
  * Deacon  _et al._ [2010a] R. S. Deacon, Y. Tanaka, A. Oiwa, R. Sakano, K. Yoshida, K. Shibata, K. Hirakawa, and S. Tarucha, Phys. Rev. B 81, 121308 (2010a).
  * Kretinin  _et al._ [2012] A. V. Kretinin, H. Shtrikman, and D. Mahalu, Phys. Rev. B 85, 201301 (2012).
  * Van der Wiel  _et al._ [2000] W. G. Van der Wiel, S. De Franceschi, T. Fujisawa, J. M. Elzerman, S. Tarucha, and L. P. Kouwenhoven, Science 289, 2105 (2000).
  * Tanaka  _et al._ [2007] Y. Tanaka, N. Kawakami, and A. Oguri, J. Phys. Soc. Japan 76, 074701 (2007).
  * Domański and Donabidowicz [2008] T. Domański and A. Donabidowicz, Phys. Rev. B 78, 073105 (2008).
  * Gräber  _et al._ [2004] M. R. Gräber, T. Nussbaumer, W. Belzig, and C. Schönenberger, Nanotechnology 15, S479 (2004).
  * Deacon  _et al._ [2010b] R. S. Deacon, Y. Tanaka, A. Oiwa, R. Sakano, K. Yoshida, K. Shibata, K. Hirakawa, and S. Tarucha, Phys. Rev. Lett. 104, 076805 (2010b).
  * Kanai  _et al._ [2010] Y. Kanai, R. S. Deacon, A. Oiwa, K. Yoshida, K. Shibata, K. Hirakawa, and S. Tarucha, Phys. Rev. B 82, 054512 (2010).
  * Xi  _et al._ [2016] X. Xi, Z. Wang, W. Zhao, J. H. Park, K. T. Law, H. Berger, L. Forró, J. Shan, and K. F. Mak, Nature Physics 12, 139 (2016).
  * Read and Newns [1983] N. Read and D. M. Newns, Journal of Physics C: Solid State Physics 16, 3273 (1983).
  * Bickers [1987] N. E. Bickers, Rev. Mod. Phys. 59, 845 (1987).
  * Bulla  _et al._ [2008] R. Bulla, T. A. Costi, and T. Pruschke, Reviews of Modern Physics 80, 395 (2008).
  * Gonzalez-Buxton and Ingersent [1998] C. Gonzalez-Buxton and K. Ingersent, Physical Review B 57, 14254 (1998).
  * Rossi and Morr [2006] E. Rossi and D. K. Morr, Physical Review Letters 97, 236602 (2006).
  * Mastrogiuseppe  _et al._ [2014] D. Mastrogiuseppe, A. Wong, K. Ingersent, S. E. Ulloa, and N. Sandler, Phys. Rev. B 90, 035426 (2014).
  * Principi  _et al._ [2015] A. Principi, G. Vignale, and E. Rossi, Physical Review B 92, 041107 (2015), arXiv:1410.8532 .
  * Ingersent [1996] K. Ingersent, Phys. Rev. B 54, 11936 (1996).
  * Hewson [1993] A. C. Hewson,  _The Kondo Problem to Heavy Fermions_ (Cambridge University Press, 1993).
  * Möckli and Khodas [2020] D. Möckli and M. Khodas, Phys. Rev. B 101, 014510 (2020).
  * Zalic  _et al._ [2019] A. Zalic, S. Simon, S. Remennik, A. Vakahi, G. D. Gu, and H. Steinberg, Phys. Rev. B 100, 064517 (2019).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018).

## Supplementary Information

### Supplementary Figure 1

Figure S1: (a) Differential conductance at zero bias ($G_{0}$) as a function
of $B_{\perp}$ at mT resolution over an extended range. Red arrow indicates
the scan direction. $G_{0}$ is the lowest at -0.2 mT, indicated by a black
arrow; Sub-gap differential conductance spectra dI/dV for magnetic field
values (b) lower than -0.2 mT and (c) higher than -0.2 mT are shown. The
colored dash-lines indicate the magnetic field at which the spectrum was
measured.

### Supplementary Figure 2

Figure S2: (a)-(f) Sub-gap differential conductance (dI/dV) spectra of QD-
NbSe2 junction along with the respective background conductance.

The individual sub-gap differential conductance spectra of the QD-NbSe2
junction, for the magnetic fields shown in Figure 4 in main text, in the range
$-0.3\text{ mV}\leq\Delta\leq 0.3\text{ mV}$ are shown in Figure S2. The red
dot-dash lines are the background conductance estimated using the function
$a_{0}+\alpha\absolutevalue{x}+\beta x^{2}$; where $a_{0}$ is the conductance
at $V_{SD}=0\text{ V}$, $\alpha$ and $\beta$ are coefficients. The values for
the $a_{0}$, $\alpha$ and $\beta$ are tabulated in Table 1.

Magnetic Field | $a_{0}$ | $\alpha$ | $\beta$
---|---|---|---
(mT) | ($\mu S$) |  |
-0.2 | 0.16 | 0 | 0
1.0 | 1.05 | 0 | 16
2.2 | 1.05 | 0.6 | 21
3.1 | 1.05 | 2.1 | 16
5.7 | 1.45 | 5 | 18
9.0 | 1.45 | 8 | 18.5
Table 1: Coefficients used to obtain the background conductance for $dI/dV$
spectra in Figure S2(a)-(f)


# High-Performance Computing in Battery Development: From Pore Scale to
Continuum

Benjamin Kellers1,2,111Corresponding author, Email<EMAIL_ADDRESS>Martin P. Lautenschlaeger1,2 Julius Weinmiller1,2 Lukas Krumbein1,2 Simon
Hein1,2 Timo Danner1,2 Arnulf Latz1,2,3 1 German Aerospace Center (DLR),
Institute of Engineering Thermodymanics, 89081 Ulm, Germany.
2 Helmholtz Institute Ulm for Electrochemical Energy Storage (HIU), 89081 Ulm,
Germany.
3 Ulm University, Institute of Electrochemistry, 89081 Ulm, Germany.

###### Abstract

An application for high-performance computing (HPC) is shown that is relevant
in the field of battery development. Simulations of electrolyte wetting and
flow are conducted using pore network models (PNM) and the lattice Boltzmann
method (LBM), while electrochemical simulations are conducted using the tool
BEST. All aforementioned software packages show an appropriate scaling
behavior. A workflow for optimizing battery performance by improving the
filling of battery components is presented. A special focus is given to the
unwanted side effect of gas entrapment encountered during filling. It is also
known to adversely affect the electrochemical performance of batteries and can
be partially prevented by appropriate microstructure design such as electrode
perforation.

##  1 Introduction

Lithium-ion batteries play a major role in electronics with applications
ranging from mobile phones to electric vehicles. Although the technology could
be considered mature, there is still room for improvement in several aspects
of, e.g., the production. One of them is the filling with liquid electrolyte.
This process is time-critical, cost-intensive and known to have a significant
influence on the electrochemical performance of batteries.

The underlying reason is that battery components typically consist of porous
media to maximize reactive surface areas and the energy density in a given
volume. However, the porosity including pore morphology and topology as well
as differences in wettability of electrolyte and ambient air induce issues
during the filling. They can result in partial gas entrapment, i.e. incomplete
electrolyte saturation. Thus, only parts of the surface of the active material
in the electrode are submersed with electrolyte, while the rest is considered
electrochemically inactive. Moreover, the entrapped air is an insulator for
ionic and electric conductivity. Technically speaking, the effective
tortuosity, i.e. the length of the shortest path between two opposite sides of
the structure, is increased. This finally leads to a decrease of the battery
performance.

Capturing these phenomena through experiments is challenging and sometimes
even impossible. Thus, theoretical modeling and simulations are instead
applied here. Electrolyte wetting is simulated using pore network models (PNM)
and the lattice Boltzmann method (LBM). The software tool BEST is used to
predict the electrochemical performance. All tools are connected in a
workflow. A critical aspect of this study is the multi-criteria optimization.
Structural parameters and material properties which improve the filling might
still worsen the battery performance. Thus, creating a feedback loop in the
process optimization is key for the proposed workflow.

The present study is structured as follows: In Sec. 2 an overview of the
workflow is given and the computational methods and the software packages are
introduced. Simulation results are presented and discussed in Sec. 3.
Information on HPC scalability of the software tools are given in Sec. 4.
Finally, conclusions are drawn in Sec. 5, where also a short outlook is given.

##  2 Workflow and software packages

The workflow is schematically shown in Fig. 1. It is described in the
following, going from top to bottom and left to right.

Figure 1: Flowchart of the workflow. The data flow and the relationships
between the different steps are illustrated. The boxes are colored with
respect to the software used. Gray labels depict steps that are not directly
related to any specific software tool.

Since pore-scale simulations are applied, 3D images of porous battery
electrodes or separators are required. In this work, lithium-ion battery
cathodes are used as an example. The corresponding structural data is either
virtually generated using physically-based algorithms [1] or comes as
grayscale images from tomography techniques (FIB-SEM or micro-CT) of real
components which requires an additional binarization step in the workflow. The
latter, for brevity, is however not described here.

In the workflow, as mentioned before, three different software tools and also
simulation techniques are used and connected. They are depicted by the three
vertical strands in Fig. 1. From the left to the right, these correspond to
PNM using OpenPNM [2], LBM using Palabos [3], and electrochemical simulations
using BEST [4], i.e. a finite volume implementation of Li-ion battery models
developed by our group [13].

For the PNMs, based on the binarized image data the pore network information
is extracted. Therefore, the open-source Python library PoreSpy [5] is used.
It reconstructs the geometrical and topological properties of the
structurally-resolved pore space using simplified geometries, which typically
means spherical pores connected by cylindrical throats (cf. Fig. 2). The
underlying procedure is based on a corrected watershed algorithm to segment
the image data. In the current work, the so-called SNOW-algorithm [6] is used
to also avoid oversegmentation; an issue often observed when segmenting energy
materials. Fig. 2 shows examples of such extracted networks. On the left-hand
side, a simplified 2D case consisting of only two pores is schematically
depicted. On the right-hand side, a full network for a 3D battery electrode
structure consisting of several thousand pores is given. After the network
extraction, the network properties, such as for example the median pore
diameter $d_{\text{50}}$ can be computed. It can also be used to determine
other physical properties of interest. The method is computationally very
efficient, since only analytical expressions are solved for pores and throats.
In the current work, a percolation into a pore network is simulated to
determine the pressure-saturation behavior during the filling process.

Figure 2: Network extraction from structurally-resolved image data. Left: 2D
pore space segmented into two pores (yellow and green). The disk sizes (blue)
correspond to the inscribed diameters of the pores and are connected by a
throat (opaque rectangle). Right: Extracted full network embedded into the
corresponding 3D image of a battery cathode. The pore color corresponds to the
pore volume. The throat diameters are depicted for topological reasons only
and do not relate to the physical size.

For the LBM simulations an additional preprocessing step of the binarized
image data is needed in which solid interfaces are differentiated from the
solid bulk. This helps to speed up the simulations. The preprocessing is done
using the Python libraries NumPy and SciPy. The filling process itself is
modeled as a two-phase flow using the Shan-Chen method [7]. Palabos [3] is
used for the simulations from which also pressure-saturation curves are
determined. They are used for validating the results determined using PNM.
Moreover, the local distribution of entrapped residual gas is finally
analyzed. For further details, the reader is referred to Ref. [8].

With respect to Fig. 1, in the workflow, LBM and PNM are connected through an
information exchange step. It is necessary to convert the LBM-specific unit
system to physical units. Therefore, geometrical information about the pore
space is needed, i.e. the median pore diameter $d_{\text{50}}$ is used in the
present study. In exchange, one pressure-saturation data set determined using
LBM is used for calibration of the PNM. Details regarding this calibration and
necessary corrections in the PNM will be discussed in a future publication of
our group. The calibrated PNM can then be used to cover huge parameter spaces
in percolation simulations to gain understanding of the influence of different
input parameters on the pressure-saturation behavior. Since the PNM
computations are extremely lightweight, this speeds up the optimization of the
filling process by several orders of magnitude.

The binarized image data of the porous electrodes also serves as input for the
electrochemical simulations. The corresponding domain consists of the
electrode, a separator, and a counter-electrode. On this geometrical basis a
set of coupled partial differential and algebraic equations consisting of
mass, charge and energy balances is solved. The parametrization of the
conservation equations is done through dedicated experiments or using
literature data.

Finally, simulation results regarding pressure-saturation data and
electrochemical performance are considered in a feedback loop, that aims to
improve both aspects while minimizing negative effects on each of them.

##  3 Results and discussion

In the following, a selection of representative results are discussed only.
For the sake of completeness the reader is referred to related publications
from our group [8, 9, 10].

Fig. 3 illustrates the final state of an LBM filling simulation and the
pressure-saturation curves corresponding to the simulation setup. It can be
seen from Fig. 3 that a relatively large amount of entrapped gas remains
inside the porous structure. This is also obvious from the LBM pressure-
saturation curves, which show that the final degree of saturation for
different contact angles between electrolyte and active material ranges
between 75 % and 95 %. The PNM does not yet model gas entrapment, such that
all PNM simulations reach 100 % saturation in this study.

Figure 3: Simulation results for electrolyte filling using LBM and PNM. Left:
Entrapped gas (gray) in the pore space at the end of the filling process.
Electrolyte and active material are fully transparent. Right: Pressure-
saturation curves determined using LBM (dots) and PNM (lines) simulations for
a selection of contact angles $\theta$.

Fig. 4 shows results of the electrochemical simulations. Two types of
electrodes are studied: (1) unstructured electrodes and (2) structured
electrodes with improved wettability. In experiments the structure is
realized, by selectively removing electrode material in a cylindrical hole
pattern with a short-pulsed laser. This procedure is mainly applied to thick
electrodes and known to improve rate capability, cycle stability, tendency to
lithium plating [10] as well as the wettability and electrolyte saturation
[11].

Fig. 4a) shows the state of charge (SOC) distribution within unit cells of the
unstructured and structured electrodes at different times. The structuring by
perforation improves the wettability and thereby provides additional transport
channels for lithium ions with low tortuosity. Both effects are able to
improve the utilization of the active material illustrated by the increase in
red areas shown in Fig. 4a). Fig. 4b) shows the accessible capacity of an
electrode with $6\,\text{mAh/cm}^{2}$ theoretical capacity: The perforation
improves the electrochemical performance for high operating currents $j$.

Figure 4: Electrochemical aspects. a) SOC distribution of lithium in
unstructured electrode (top) and structured electrode (bottom). The
structuring improves the homogeneity of the SOC distribution. b)
Electrochemical performance for operation currents of $j\geq
6\text{mA/cm}^{2}$.

##  4 Necessity and benefits of using HPC infrastructures

Multiple Python libraries and some custom in-house code are involved in the
described workflow. More specifically, the scientific stack of Python
consisting of NumPy and SciPy as well as Scikit-Image [12] were used. These
partly rely on precompiled C-code which often comes with multi-threading
capabilities. However, they do not utilize MPI optimization out-of-the-box.
Thus, all steps in the aforementioned workflow, for which Python was used,
were run on a single CPU only, though making use of multiple cores or threads.
The same applies for the network extraction using PoreSpy as well as the PNM
simulations using OpenPNM, as both are based on the scientific stack of
Python, too.

Therefore, performance characteristics of precompiled Python C-code are shown
in Fig. 5 using the example of a matrix multiplication with NumPy. However,
also with respect to PoreSpy, full use of the HPC environment was made
regarding the usage of RAM size and memory bandwidth. Its network extraction
algorithm can be very RAM-demanding, since multiple copies of the original
image are required for the necessary computations. Palabos, on the other hand,
hugely benefits from an HPC environment, since it is fully MPI optimized. The
models used require only local information which is extremely favorable for
multi-core parallelization. Fig. 5 shows the huge difference the MPI
optimization makes when comparing to regular multi-threaded scaling.

Figure 5: Performance Test. Left: Scaling of fully resolved and homogenized
LBM electrolyte filling simulations using Palabos. Right: Scaling of a matrix
multiplication in Python using NumPy in a multi-threaded environment to
illustrate the scaling of non-MPI computations in the workflow.

##  5 Conclusions and Outlook

This study shows the potential power of combining different simulation tools
into a workflow by information exchange and feedback loops. Often in process
optimization, parameters have a negative effect on each other such that trade-
offs are necessary to improve the process as a whole. This is showcased here
by the interdependence of electrolyte filling and electrochemical performance
of batteries. Another aspect highlighted in this study, is the fact that
software which is not HPC-optimized per se, can still benefit from HPC
environments. Additionally, coupling of low and high fidelity models result in
the necessary efficiency to use such a workflow on large parameter spaces.

As part of our future work, this study will be extended in two directions.
First, by performing electrochemical simulations on partially saturated
electrodes determined using LBM simulations. Second, by simulating the effect
of structured or laser-perforated electrodes on the electrolyte filling.

## Acknowledgements

We gratefully acknowledge financial support from the European Union’s Horizon
2020 Research and Innovation Programme within the project "DEFACTO" [grant
number 875247] for parts of the LBM and PNM activities of our research as well
as from the ‘Bundesministerium für Bildung und Forschung’ within the project
HiStructures under the reference No. 03XP0243E. This work also contributes to
the research performed at CELEST (Center for Electrochemical Energy Storage
Ulm-Karlsruhe). The simulations were carried out on the Hawk at the High
Performance Computing Center Stuttgart (HLRS) [grant LaBoRESys], and on JUSTUS
2 at Ulm University [grant INST 40/467-1 FUGG].

## References

  * [1] D. Westhoff, et al. Generation of virtual lithium-ion battery electrode microstructures based on spatial stochastic modeling. Comp. Mat. Sc. 151 (2018) 53.
  * [2] J. Gostick, et al. OpenPNM: A Pore Network Modeling Package. Comp. Sc. & Eng. 18 (2016) 60.
  * [3] J. Latt, et al. Palabos: Parallel Lattice Boltzmann Solver. Comp. & Math. w. Appl. 81 (2021) 334.
  * [4] ITWM. BEST - Battery and Electrochemistry Simulation Tool. 2014 URL: https://www.itwm.fraunhofer.de/best
  * [5] J. T. Gostick, et al. PoreSpy: A Python Toolkit for Quantitative Analysis of Porous Media Images. JOSS 4 (2019) 1296.
  * [6] J. T. Gostick. Versatile and efficient pore network extraction method using marker-based watershed segmentation. Phys. Rev. E 96 (2017) 023307.
  * [7] X. Shan and H. Chen. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev. E 47 (1993) 1815–1819.
  * [8] M. P. Lautenschlaeger, et al. Understanding Electrolyte Filling of Lithium-Ion Battery Electrodes on the Pore Scale Using the Lattice Boltzmann Method. Batter. Supercaps 5 (2022) e202200090.
  * [9] M. P. Lautenschlaeger, et al. Homogenized lattice Boltzmann model for simulating multi-phase flows in heterogeneous porous media. Adv. Water Resour. 170 (2022) 104320.
  * [10] V. De Lauri et al. Beneficial Effects of Three-Dimensional Structured Electrodes for the Fast Charging of Lithium-Ion Batteries. ACS Appl. Energy Mater., 4 (2021) 13847.
  * [11] M. J. Kleefoot, et al. Enhancement of the wettability of graphite-based lithium-ion battery anodes by selective laser surface modification using low energy nanosecond pulses. Int. J. Adv. Manuf. Technol. 118 (2022) 1987–1997.
  * [12] S. van der Walt, et al. scikit-image: Image processing in Python. PeerJ 2 (2014) e453.
  * [13] A. Latz, et al. Multiscale modeling of lithium ion batteries: thermal aspects. BJNANO 6 (2015) 987.