# Ca-dimers and solvent layering determine electrochemically active species in
Ca(BH4)2 in THF

Ana Sanz-Matias  Joint Center for Energy Storage Research, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA  The Molecular Foundry, Lawrence
Berkeley National Laboratory, Berkeley, CA 94720, USA  Fabrice Roncoroni
Joint Center for Energy Storage Research, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, USA  The Molecular Foundry, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA  Siddharth Sundararaman  Joint
Center for Energy Storage Research, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA  The Molecular Foundry, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, USA  David Prendergast  Joint Center for
Energy Storage Research, Lawrence Berkeley National Laboratory, Berkeley, CA
94720, USA  The Molecular Foundry, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA<EMAIL_ADDRESS>
###### Abstract

Divalent ions, such as Mg, Ca, and Zn, are being considered as competitive,
safe, and earth-abundant alternatives to Li-ion electrochemistry. However, the
challenge remains to match electrode and electrolyte materials that stably
cycle with these new formulations, based primarily on controlling interfacial
phenomena. We explore the formation of electroactive species in the
electrolyte Ca(BH4)2 in THF through molecular dynamics simulation. Free-energy
analysis indicates that this electrolyte has a majority population of neutral
Ca dimers and monomers, albeit with diverse molecular conformations as
revealed by unsupervised learning techniques, but with an order of magnitude
lower concentration of possibly electroactive charged species, such as the
monocation, CaBH${}_{4}^{+}$, which we show is produced via disproportionation
of neutral Ca(BH4)2 complexes. Dense layering of THF molecules within 1 nm of
the electrode surface (modeled here using graphite) hinders the approach of
reducible species to within 0.6 nm and instead enhances the local
concentration of species in a narrow intermediate-density layer from 0.7-0.9
nm. A dramatic increase in the monocation population in this intermediate
layer is induced at negative bias, supplied by local dimer disproportionation.
We see no evidence to support any functional role of fully-solvated Ca2+ in
the electrochemical activity of this electrolyte. The consequences for
performance and alternative formulations are discussed in light of this
molecular-scale insight.

## Introduction

To increase the rate of conversion to renewable energy sources,
electrification of various energy-intensive aspects of society is underway.
The concomitant demand for electrochemical energy storage solutions
increasingly highlights the limits of Li-ion technologies with respect to
performance, safety and sustainability. Multivalent ions such as Mg2+, Ca2+,
Zn2+, or even Al3+, offer more earth-abundant alternatives, some with higher
theoretical specific capacity and reduced safety concerns due to self-
passivation of metal anodes.1, 2, 3, 4, 5

However, realizing performant electrochemical cells using these ions is
hindered due to various issues driven by interfacial phenomena: low power
output and charging rates due to large overpotentials and associated
electrolyte decomposition and interphase growth.6, 7, 8, 9

At issue is our lack of understanding of the specific complex nature of
solvation in suitable electrolytes for multivalent ions and the identification
of which of these species are active at the electrode-electrolyte interface
and why.10 We may be biased by familiarity with aqueous solutions and the
ability of water to generate perfect electrolytes, with fully dissociated and
solvated ions, for many salts. However, organic solvents (required for a
sufficiently wide window of electrochemical stability in batteries), with
typically lower dielectric constants and larger molecular sizes, have
increased residence times for coordinating highly charged cations and have
relatively little interaction with anions, unlike ambipolar water molecules
which more easily solvate both charges.

Molecular dynamics provides a window to the inner workings of electrolytes,
revealing details of coordination of cations by solvent molecules and anions.
However, in the study of highly-charged species in poor dielectrics, we must
take care to avoid sampling only a limited set of coordination states due to
their long lifetimes and unavoidable limitations in computing time and
complexity. Free-energy sampling allows us the opportunity to pose fundamental
questions regarding the chemical composition of a poor electrolyte and the
mechanisms by which its solvated species interconvert.11, 12, 13, 14, 15, 16
Furthermore, mining the large quantities of compositional and conformational
data produced by these simulations presents its own challenge. Here we rely on
recently developed unsupervised learning approaches17 to provide a faster path
to extracting molecular-scale insight and guidance for future experiments to
validate our predictions.

We study Ca(BH4)2 in THF as a promising electrolyte candidate,18 spurred on by
recent studies of the bulk electrolyte 19 and the anode interface.6, 20 We can
also contrast its behavior with Mg(TFSI)2 in THF, studied recently using free
energy sampling.21 Ref. 19 proposed that neutral Ca dimers (Ca2(BH4)4)
facilitate the disproportionation of neutral monomers into (active)
monocations CaBH${}_{4}^{+}$ and anions Ca(BH4)${}_{3}^{-}$. In the same work,
molecular cluster calculations (embedded in a polarizable continuum model)
indicated that the dimer is the second most stable conformation after the
neutral monomer, but lacked specific solvent interactions beyond the first
coordination shell and any Debye screening from finite ion concentrations.
Interfacial characterization reveals the ready formation of solid-electrolyte
interphases incorporating oxidized boron and even embedded calcium hydride.6
And debate continues as to the presence or electrochemical relevance of the
fully-solvated Ca2+ dication. 18, 22, 23, 24, 25, 26

Figure 1: Analysis of the bulk electrolyte. (a) One time step sampled from our
molecular dynamics model of the bulk electrolyte showing two Ca2+ dications
and 4 BH${}_{4}^{-}$ anions in THF at room temperature (RT). (b) Populations
of electrolyte species derived from their respective free energies $\Delta G$
(kT), which were obtained by integration of (c) the free-energy surface
sampled using metadynamics with respect to Ca-Ca distance and Ca-B(H4)
coordination number. Minimum energy pathways for dimer disproportionation from
the neutral species, Ca(BH4)2 are indicated by dashed lines. (d) Population
analysis with respect to anion and THF coordination and conformation about Ca
ions obtained via unsupervised learning, indicating the diversity of species
umbrella-sampled from points on the free-energy surface corresponding to the
neutral monomer (CN(Ca-B) = 1.9 and d(Ca-Ca) = 11 Å) and the long (L) and
short (S) dimers (d(Ca-Ca) $<$ 7 Å). Atomic structures of dominant structures
(populations larger than 7.5%) are shown.

In this work, we reveal intricate details of the population of species in the
bulk of this electrolyte and striking differences within a nanometer of its
electrode interfaces, which are further exaggerated by potential differences.
We discuss the consequences of these predictions for functioning cells.

## Results

### In the bulk electrolyte

We employ free-energy analysis for a model of the bulk electrolyte at room
temperature (RT), comprising two Ca2+ ions (and four borohydride anions)
dissolved in THF (Fig. 1 a), using empirical force fields (see Methods and
Supporting Information). The free energy surface (Fig. 1 c) is sampled using
metadynamics27 with respect to two collective variables: the Ca-B coordination
number, which controls the charge of complexes, and the Ca-Ca distance, which
for a system with only two Ca ions, distinguishes between monomers and dimers.
Integration over the free energy surface indicates that the most favorable
species in solution are Ca-Ca dimers (Fig. 1 b). These are neutral complexes
with the formula Ca2(BH4)4 and come in two varieties: long dimers (36%)
bridged by a single BH${}_{4}^{-}$ and short dimers (21%) bridged by two
anions (see below for more detail). The next most common species is the
neutral monomer complex, Ca(BH4)2 (40%). It is notable that only a small
percentage (3%) of the solvated population is predicted to be charged, the
monocation, CaBH${}_{4}^{+}$, with its complementary anion,
Ca(BH4)${}_{3}^{-}$. The free energy for the formation of fully solvated
dications, Ca2+, in this model, is too high ($\sim$19 kT) to support a
significant population. These results provide some reordering with respect to
static estimates which predict the neutral monomer as most favorable, followed
by a single (short) dimer conformation, the monocation and the anion. 19 From
the relative populations of each species, it is clear that neglecting the long
dimer would lead to this conclusion. Furthermore, it is clear from RT sampling
that there are multiple possible conformations for the dimer at finite
temperature (as we explore below). However, both sets of calculations agree
that the fully solvated dication Ca2+ is the least favored in this set of
possibilities (1-1.2 eV higher in energy by static estimates 19, 0.48 eV from
free-energy sampling at RT).

Our 2D free energy surface (Fig. 1 c) reveals that direct interconversion of
these solvated objects, by removal or addition of borohydride anions, is
prevented by quite steep free energy barriers (24.4 kT or $\sim$0.6 eV). The
easier path to disproportionation (forming charged species from neutral
monomers) is through the formation of dimers, with exchange of borohydride
anions before dissociation into the complex monocation and anion (with
activation energies of 2.4 - 10 kT, Fig. 1 c and Fig. S1):

$\ce{2Ca(BH4)2->Ca2(BH4)4->CaBH4^{+}+Ca(BH4)3^{-}}$

At a minimum, this explains the prevalence of dimers in spectroscopic analysis
of the bulk electrolyte (using EXAFS and Raman)19 and the observation of
saturating ionic conductivity with increasing concentration as more neutral
dimers form (prior to precipitation). 19, 28

Strikingly, we find that bulk populations of each complex are conformationally
quite diverse. Our metadynamics simulations project the full free energy
landscape onto only a few collective variables, however, multiple molecular
conformations may satisfy these constraints (Ca-B coordination number and Ca-
Ca distance, in this case). Data-mining techniques (dimensionality reduction,
hierarchical clustering and permutation-invariant alignment 17 – details in
Methods and SI) applied to selective umbrella sampling of local minima in the
free energy landscape reveal a rich variety of solvated isomeric structures,
summarized in Fig. 1 d.

For example, the neutral monomer, Ca(BH4)2, when additionally coordinated by
three THF molecules ($\sim$ 60% of its population) adopts mostly bent
borohydride arrangements with a large dipole moment (8 D, see SI). In
addition, we found a significant portion ($\sim$ 30%) of monomers coordinated
by 4 THF molecules with an axial borohydride arrangement and low dipole moment
(2 D). These two structures had been proposed separately as the minimum energy
structure from quantum-chemical cluster calculations 19 and molecular dynamics
simulations, 28 respectively.

By contrast, the monocation occurs mostly with 5 coordinating THF molecules,
in agreement with previous reports.19

We find that dimers are always neutral (with four borohydride anions) and are
split in two main spatial configurations characterized as short (SD) and long
(LD) dimers with average Ca-Ca equilibrium distances of 4.48 and 5.25 Å,
respectively (Fig. S2). Furthermore, each was found to have sub-populations
with 4-7 solvent molecules, and, within those, several stereoisomers (Fig. 1
d). Key differences between the long and short dimers are the presence of a
single-anion or double-anion bridge and predominant 6 THF coordination or
mixed 5-6 THF coordination, respectively. Short dimers are in excellent
agreement with a previously proposed double-bridged dimer structure with a 4.4
ÅCa-Ca distance (from EXAFS fitting data19). The set of structures shown here
expands upon and underscores the configurational flexibility of calcium. 19
And, based on our understanding of the efficient disproportionation pathways
to form charged species via dimerization (discussed above), it makes sense
that the dominant dimer conformations form from combinations of bent and axial
borohydride arrangements of the neutral monomers (e.g., long dimer isomers 1,
3 and 7 with 6 THF molecules in Fig. 1 d are bent-bent, axial-bent and axial-
bent combinations, respectively).

Although all dimers here are neutral, each Ca ion within a given dimer may be
locally coordinated by 1 to 4 anions, with 1 (long) or 2 (short) shared
between them. Most commonly we observe [3,2] or [3,3] anion arrangements for
long or short dimers, respectively. Small populations of [1,4] dimers are
found and are key to some interfacial disproportionation processes discussed
below and in Fig. 5.

Figure 2: Effect of solvent layering at the interface. (a) The free-energy
profile (red) of a single THF molecule in THF at RT with respect to distance
from a graphite interface and the corresponding oxygen density profile
(black), with slight adjustment (dashed lines) at negative bias. The Dense
Layer (DL), Intermediate Density Layer (IDL) and bulk regions are color-coded
hereafter in red, blue and gray. The inset shows a snapshot of the interfacial
region sampled from molecular dynamics indicating THF molecules lying flat
against the graphite surface and somewhat constrained at the DL-IDL interface.
(b) The 3D free energy landscape derived from metadynamics with respect to Ca-
Ca distance, C-B coordination number, and distance from the graphite surface
for the neutral interface. Minimum free-energy paths for (c) neutral and (d)
charged species in the electrolyte arriving at the interface from the bulk
under zero (solid) and negative (dashed) bias conditions indicated by the
graphite surface charge density, $\sigma$. Sudden jumps in the free-energy
profiles reflect variations with respect to other collective variables in the
free-energy landscape. Under negative bias the long and short dimers are not
distinguished in our sampling. In general, negative bias stabilizes charged
species at the interface. Figure 3: Distribution of species at the neutral
interface. Molecular model of the simulation box, obtained from a snapshot of
the equilibrated MD trajectory (a), showing a dimer at the IDL (in blue), near
the DL (in red). Integrated free-energy and corresponding population per layer
(b), color-coded as red, blue and gray for the DL, IDL and bulk, respectively.
Unsupervised clustering analysis of umbrella sampling trajectories at the DL
and IDL (z = 5.75 and 8.25 Å) show the relative populations of dimer isomers
per layer (c,d). Representative structures have discrete orientations, which
depend on the layer (insets in e,f). Color lines show the orientation of the
Ca-Ca axis with respect to the surface normal, where 90° is parallel to the
surface, as shown in IDL dominating dimer LD 6 THF - 7 (also in a), while DL
dimer SD 5 THF - 4 is mostly perpendicular. Dipoles have also discrete
orientations with respect to the surface normal for given isomers (shown in
black, dotted lines) and are roughly aligned to the Ca-Ca axis in long, flat
dimers (IDL) or perpendicular, as in DL dimer SD 6 THF - 2.

### At the electrode-electrolyte interface

What happens to species in the electrolyte as they approach an interface, such
as the electrode surface? Firstly, as expected from simple statistical
mechanics modeling of molecules at hard interfaces, 29 the solvent, THF,
adopts a layered molecular structure 30, 21 near the surface with a dense
layer (DL) at 3-6 Å, followed by a low-density ‘gap’, and an intermediate
density layer (IDL) at $\sim$ 7-10 Å from the surface, with 2.5, 0.3 and 1.5
times the bulk THF density, respectively (Fig. 2 a). A third collective
coordinate, calcium distance from the interface, sampled with metadynamics
allows us to obtain the interfacial, 3D free-energy landscape (Fig. 2 b), of
which the most distant slice is the 2D bulk free-energy landscape in Fig. 1 c.
Dissolved species in the electrolyte near the electrode surface respect this
underlying solvent structure. Now, free-energy minima are distributed between
the bulk, IDL, and DL, and separated by barriers, as indicated by minimum
energy pathways of neutral and charged species approaching the interface in
the 3D landscape (Fig. 2 c and d).

A key observation is that the IDL defines an attractive basin for most
species, especially dimers and monocations, with minima in free energy that
are lower than in the bulk, implying that this is a narrow interfacial region
for enhanced concentration of solutes. Conversely, the DL defines a region
from which solutes may be excluded due to additional free energy costs,
without assistance of some applied bias (see below). This has some generally
important consequences for electrochemistry in this electrolyte at a non-
interacting electrode (such as graphite). Negligible populations of species
that might be specifically adsorbed or adjacent to the electrode would imply
the absence of inner-sphere electron transfer events during reduction.
Similarly, if the electroactive species are likely to be highly coordinated in
the IDL, then the necessary outer-sphere electron transfer events may very
well lead to reductive decomposition of the coordinating species (solvent or
anions) leading to low Coulombic efficiency and the formation of a solid-
electrolyte interphase (SEI). Based on this study, our claim is that this
inefficiency is as much a consequence of the strong solvent layering at the
interface as it is of the salt that defines the electrolyte.

Specifically, we find that, next to an unbiased non-interacting electrode
(Fig. 3 a), the IDL is dominated ($\sim$ 75%) by dimers (long dimers in
particular, 55%), with reduced populations of neutral monomers ($\sim$ 10%)
and a very slight increase to 8% in the monocation population (Fig. 3 b). Only
6% of the interfacial population is in the DL, according to thermodynamic
integration of the free energy. In the DL, the monocation population decreases
and dimers dominate ($\sim$ 85%) – especially short dimers, in contrast to
bulk and IDL populations.

Unsupervised clustering analysis of structures from umbrella-sampling
trajectories at the IDL and DL (z = 8.25 and 5.75 Å, respectively) reveals a
reduced number of favored dimer isomers compared to the bulk, with one given
isomer making $>$20 % of the population in each case (Fig. 3 c-d). In the IDL,
this is a bent-axial long dimer with 6 solvating THF molecules (indexed as
Isomer 7 or LD 6 THF - 7), already a favored species in the bulk. In the DL,
on the other hand, a short dimer also with 6 THF molecules (SD 6 THF - 2)
dominates. Furthermore, we find that the orientation of dimers at the
interface is discretized (Fig. 3 e,f). Favored IDL dimers have mostly flat
orientations of the Ca-Ca vector relative to the graphite surface, with two
dense layer THF molecules involved in solvation. Similarly, the dominant dimer
orientation in the DL is flat (i.e., in a plane parallel to the surface), with
both calcium ions embedded in that dense region, and with two THF molecules
from the IDL contributing to solvation. Additionally, perpendicular dimers
form at least 11% of the DL population. Isomer 4 of the short dimer with 5 THF
molecules (SD THF 5 - 4) is an asymmetric dimer, with calcium ions solvated by
2 and 3 THF molecules that sit in the DL and on the edge of the IDL,
respectively (Fig. S3). Note that the specific conformation of these
coordination complexes defines their effective dipole moment, which may, in
some cases, be orthogonal to the Ca-Ca vector of the dimer. An example in the
DL is the favored isomer SD 6 THF - 2, with its Ca-Ca vector parallel to the
surface but three BH${}_{4}^{-}$ sitting closer to the interface, near the
corresponding DL free energy minimum for isolated borohydride in THF next to
graphite (z= $\sim$4 Å, see Fig. S4). Dipole orientation will be discussed in
more detail below in the context of biased interfaces.

As observed in the bulk electrolyte, the absence of fully solvated dications
is more pronounced at the interface. This speaks to the rigidity of the THF
coordination sphere around this multivalent ion (similar to that around Mg as
seen in Ref. 21) whose size (an effective radius $\sim$ 5 Å) prevents the
close approach of this species to the electrode. The preferred, relatively
flat orientation of THF molecules in the dense layer is incompatible with the
preferred radial coordination of the cation and likely forces a reduction in
coordination number that raises the free energy of fully-solvated Ca2+
(reducing its relative population) relative to complexes coordinated by the
more compact borohydride anions.

### Biased Interfaces

So far, it seems that Ca(BH4)2 in THF is a poor electrolyte, with only a small
fraction of the salt concentration (3%) present as charged species in the
bulk, albeit with a noticeable enhancement to 10% in the intermediate density
layer (IDL). Otherwise this electrolyte is dominated by neutral species
(monomers and dimers). This is consistent with previous discussion of
undissociated neutral species as dominant in Ca-based electrolytes with boron-
containing anions10. However, the bulk free energy landscape (Fig. 1 c)
indicates that interconversion of species is possible (more details below) and
some equilibrium exists between charged and neutral species. We explore the
impact of biased/charged electrodes on the population of solutes by evenly
distributing opposing charges on either face of the two-layer graphite
electrode model, which, under periodic boundary conditions, polarizes the
electrolyte. We considered two specific charge states with (1) 0.065 and (2)
0.13 e/nm2. We estimated that these surface charge densities raise the Fermi
level of a graphene model by 0.56 and 0.74 eV, respectively (see Sec.
Supplementary Information 2). Screening at the interface by the solvent and
electric double layer will result in a much smaller potential difference
across the interface. A lower bound for this bias is 0.038 and 0.065 V,
respectively, provided by the Grahame equation, assuming a simulated bulk
concentration of 0.12 M.31 (Note that at higher bulk concentrations, with a
shorter Debye screening length, this estimated potential difference should be
even smaller – at 1.65 M 6 the calculation would indicate 0.01 V and 0.02 V,
respectively). The larger of these electrode surface charge densities is not
enough to draw monocations into the dense layer, which are held outside this
region by a barrier of approximately 0.25 eV.

Figure 4:  Populations at a biased interface. Snapshot obtained from the
equilibrated trajectory of charged state 2 ($\sigma$ = $\pm$ 0.13 e nm-2) with
highlighted layering (red/blue for DL/IDL) at the negative interface showing
the IDL solvated monocation (a). Free-energy and population distribution ((b)
at the negatively charged interface is dominated by the monocation, most
favored at the IDL. Unsupervised clustering analysis of Umbrella Sampling
trajectories of the monocation at the IDL, close to the edge of the DL (z =
5.8 Å) and close to the center of the DL (z = 4.8 Å) (c) classifies the
structures into two isomers of the five-fold and four-fold THF coordinated
monocation, which differ by slight changes in the local geometry of the
coordinating THF molecules. The corresponding Ca-B dipole orientation with
respect to the surface normal ($\vec{\mu}$), with 180∘ being the B pointing
away from the surface, indicate discrete orientation at the interface (d).
Representative structures of the main isomers at each layer in their most
likely orientation, with the surface on the left side (e).

From Fig. 2 (a) we see that the solvent layering remains practically unchanged
upon charging the electrode. Dimers are still favored in the intermediate
density layer (IDL) and further from the interface in charge state 1 (Fig.
S6). However, in charge state 2 some charged species become strongly
stabilized at the interface. Specifically, at this larger negative bias the
most favored species within the intermediate (IDL) and dense(DL) layers is now
the monocation, CaBH${}_{4}^{+}$, which almost entirely displaces the neutral
dimers and monomers (Fig. 4 a, b). Due to the higher free energy of species in
the DL, the interfacial population is largely limited to the IDL (99.4%),
according to thermodynamic integration. At this potential, the approach of the
monocation to the electrode, through the DL, would have to be an endergonic
process, requiring 5.7 kT of free energy, overcoming a barrier of 7.5 kT.
Although, this is some improvement over the case of the neutral electrode,
which presents a barrier of 9.9 kT to enter the DL, with a required input of
9.1 kT of free energy.

We can understand these free energy costs by following the evolution in
solvation and associated dipole orientation of the monocation. Unsupervised
learning analysis of various umbrella sampling trajectories: in the IDL, at
the edge of the DL and in the DL (z = 8.25, 5.8 and 4.8 Å, Fig. 4 c-e),
reveals that five-fold THF coordination dominates the IDL population, with the
BH${}_{4}^{-}$ anion pointing away from the surface – as one might expect
given the direction of the electric field at the negative electrode. However,
this dipole reorients at the edge of the DL, likely due to mixed DL/IDL
solvent coordination, and a small four-fold coordinated population appears.
Ultimately, in the center of the DL, four-fold coordination dominates, with
the dipole of the predominant isomer pointing away from the surface again.

This complicated and costly path for the monocation to reach the electrode
further emphasizes the important role of the THF solvent layering and specific
coordination in determining the electrochemical activity. By the same token,
the fully-solvated dication, Ca2+, with its somewhat rigid first solvation
shell, is still too unfavorable to define a noticeable population at the
negative electrode, despite its higher electrostatic charge. In fact, we find
that the dication is only slightly more stable than at the neutral interface
(Fig. 2 d).

Figure 5: Disproportionation pathways Scheme depicting the main two steps of
dimer disproportionation (a): reorganization (1-5), followed by separation
into ions (6) through distinctive paths A and B. Minimum energy pathways of
disproportionation at the neutral IDL and bulk (b) and at the charge state 2
IDL and bulk (c), obtained from the 3D free-energy landscapes. Since each
point can be labelled according to the sampled collective variables, we have
color-coded the distance to the surface in the IDL and bulk tones (blue,
gray).

### Generation of active species

Thus far, it seems that the monocation, CaBH${}_{4}^{+}$, is the strongest
candidate for the electroactive species in this electrolyte, given its high
population in the vicinity of the electrode upon charging to a sufficiently
high potential difference. We do not explore reduction potentials in this
work. However, as we have seen, the appearance of the monocation with close
proximity to the interface, in the IDL, indicates that it is a likely
candidate for an outer-sphere electron transfer process. With increasing
potential difference, to overcome the free energy costs associated with
penetration of the dense THF layer (DL) and reduced solvation, we would
predict an inner-sphere reduction process. However, the fact remains that this
poor electrolyte (only 3% of species in solution are charged, as mentioned
above) must supply the IDL with monocations in the first place, via some
disproportionation mechanism from neutral species (most likely dimers), and
replenish the same species while electroreduction and electrode deposition
consumes them. Any barriers in this supply chain should be evident in the
kinetics of the electrochemistry as an observed deposition overpotential or
associated rate limitations in charging cells with this electrolyte.

To shed light on the underlying processes, we approach the generation of
charged species (monocations) as a two-step process, involving dimer
reorganization and disproportionation. As we show below, the most stable or
prevalent dimer species are not readily disposed to disproportionate. Some
molecular rearrangement of borohydride anions and solvent molecules is first
required. The subsequent disproportionation follows two major pathways, which
can occur in the bulk electrolyte or at the interface, and which we have
investigated both at neutral and negative biases.

In the neutral cell, increased dimer concentrations in the intermediate
density layer (IDL) are readily explained through analysis of the free energy
landscape (Fig. 2). The minimum energy path (MEP) for disproportionation in
the presence of the interface indicates that neutral species in the bulk can
easily flow, without a significant free energy barrier, into the IDL.
Migration from the bulk to the IDL is essentially barrierless for all species
except Ca2+. As summarized in Fig. 5 a, to prepare for disproportionation of
the most favorable dimers, some reorganization into an intermediate (less
stable) dimer is required and ultimately Ca-Ca separation into ionic species
follows two main paths, A or B.

For Path A, the dominant long dimer configuration (1 in Fig. 5 a) can undergo
a low-barrier reorganization, through a short dimer (2), to another long dimer
with similar coordination of Ca ions with borohydrides (3). Then, the Ca-Ca
distance increases (4) up to 6.7 Å at the transition state so that the nascent
monocation Ca ion can increase its number of coordinating THF molecules, from
the original 3-4 to the preferred 5, upon full dissociation. Path B branches
from the short (2) or intermediate long (3) conformations, through further
borohydride and solvent reorganization, to a higher energy dimer (5), with a
[1,4] anion coordination, which then disproportionates to produce the
monocation (6).

Figures 5 b and c outline the MEP for disproportionation at the electrode
interface or in the bulk electrolyte under neutral or negative bias. Overall,
at the electrode interface, in the intermediate density layer (IDL),
disproportionation follows Path A, whereas Path B is preferred in the bulk,
likely due to configuration 5 being more favorable in the bulk than at the
interface (Fig. S7).

We find that interfacial (IDL) disproportionation at zero bias via Path A is
favored, since it has slightly lower barriers (Ea,3→4 = 8.4 kT and Ea,4→6 =
9.8 kT) and a lower free energy cost than bulk disproportionation via Path B
(Fig. 5 b). This is in agreement with the slight increase in monocation
population observed at the IDL in Fig. 3 (Slight variations in barriers
between bulk and interface models can be due to differences in the collective
variables and grid-spacing employed in our metadynamics simulations, SI Table
1).

At the negatively charged electrode, we have seen (Fig. 4) that the monocation
is the most favored species in the IDL, completely displacing the previously
dominant dimers with increasing negative charge on the electrode (charge state
2). Due to the size limits of our simulations, we may well expect that the
bulk thermodynamics are somewhat different from those under neutral
conditions, however, disproportionation still follows Path B in our
simulations (Fig. 5 c), albeit with an additional step involving the formation
of a long dimer conformation (3). Similarly, Path A is still preferred in the
IDL. Although disproportionation occurs in the bulk, barrier-less pathways to
the IDL suggest that dimers can approach the charged IDL and undergo
disproportionation quite favorably. Therefore, upon negative charging, the
concentration of monocations should increase in the IDL, based on a favorable
free energy and necessary supply from the local (IDL) dimer population via
interfacial disproportionation.

This free-energy analysis shows little activity in the dense layer (DL) due to
its low population of dissolved species. However, as noted above, the slight
stabilization of the monocation in the DL at negative bias results in a slight
increase within the overall interfacial (DL plus IDL) population of this
species – from $\sim$0.2% in the neutral DL to $\sim$0.6% in the DL at charge
state 2. Structural analysis indicates that these monocations have their
anionic end pointing away from the negative surface and have lower solvent
coordination. We have stated before that electrochemical activity is a strong
(exponentially decaying) function of distance of the reducible species from
the electrode. It is also very likely enhanced by reduced cation coordination.
The question remains whether these benefits outweigh the low population of
such species in the net electron transfer rate at similar potential
differences.

## Discussion

Based on our analysis of the Ca(BH4)2—THF electrolyte and its interfacial
speciation, we can propose the following phenomenology that may explain
existing observation and provide guidance and interpretation of future
characterization efforts.

First and foremost, the solvent, THF, and its strong solvation and interfacial
layering dictate much of what we have observed. The strong solvation of the
dication essentially prevents it from taking part in electrochemistry. And the
monocation requires a threshold negative bias to begin to dominate as an
interfacial species, albeit confined to the intermediate density layer (IDL),
7-10 Å from the electrode surface. The dense solvent layer (DL) is only
sparsely populated with solvated or partially solvated species.

That negative bias can begin to drive the interfacial population of charged
species significantly above that present in the bulk indicates that this poor
electrolyte is “activated” upon charging. Therefore, meaningful
characterization of the electrochemical activity of this electrolyte requires
operando measurements that are sensitive to within a nanometer of the
interface. The rich isomer subpopulations with distinct orientations in the
IDL make it an excellent playground for interfacially-sensitive, polarization-
dependent spectroscopies that can capture these differences. Furthermore, the
bias-dependence switch in local population from oriented dimers to monocations
should be observable with chemically-sensitive vibrational32, 33, 34 and
electronic probes.35, 36, 37

The required potential differences to enrich the IDL with charged species and,
presumably, the disruption and population of the DL with more charged species
at even higher potentials, may be consistent with known overpotentials for Ca
deposition. Plating and stripping of Ca using this electrolyte carries a
first, short duration overpotential of $\sim$250 mV, followed by a $\sim$100
mV overpotential in subsequent cycles.18 The remaining free-energy cost (5.7
kT) for monocations to access the DL in our higher charge state simulations
(estimated to be at 0.065 V), implies that the final potential difference to
enrich the DL with monocations for reduction is $\sim$ 0.2 V. It remains to be
determined how a more strongly interacting electrode surface (which may be
present after the first cycle) alters the THF solvent layering and, thereby,
the thermodynamics of electrolyte species within the first nanometer of the
evolving electrode surface. Strong solvent and anion coordination of
electrochemically active species, which has dominated our analysis, is very
likely the source of solid-electrolyte interphase (SEI) formation due to
electroreduction and decomposition of these ligands at the interface.

Although all small molecule polar solvents should present dense layers at
interfaces, the inaccessibility of the DL is also, in part, determined by the
choice of anion, whose electrode surface affinity, solvophobicity, size and
cation coordination strength may facilitate overcoming the molecular packing
in the DL and reduce associated overpotentials, ideally without the production
of a thick interphase due to its own electrochemical instability. An example
of a solvophobic anion that can bridge the THF dense layer is TFSI-, which,
according to previous free-energy studies, is more favorable in the DL and IDL
regions than in the bulk solvent ($\sim$2 kT)21. This contrasts with
BH${}_{4}^{-}$ (Fig. 2 b) which has reduced populations at the interface,
especially in the DL. We have seen that Mg, another divalent cation, which
when fully-solvated dication is also hindered from close approach to the
neutral electrode, but when combined with TFSI- can arrive stably at the
interface, effectively in the dense layer (barrier $<$4 kT, $\Delta$G $<$ 3
kT, z=6.3 Å). A non-sparsely populated dense layer with significant monocation
population seems likely in the MgTFSI2—THF electrolyte. Yet, we know, TFSI is
both inherently and electrochemically unstable, both for Mg38 and for Ca.39

The instability of TFSI may be driven by its tight coordination of cations
through its sulfone groups. Increased electroreductive stability in large
anions that simultaneously span the small solvent molecule dense layer may be
afforded by considering non-cation-coordinating species. For example,
closoboranes have been studied with Mg in tetraglyme.40 These bulkier anions
may also lead to better electrolytes overall (for example, preventing the
formation of dimers), readily producing charged electroactive species.
Although, strong solvent coordination, as we have seen for the fully-solvated
Ca2+, may still lead to significant overpotentials for electrodeposition and
associated solvent decomposition and SEI growth.

From the computational perspective, we have highlighted the value of free-
energy exploration and unsupervised learning to reveal the complexity of this
nominally simple electrolyte and tried to connect our simulations to observed
electrochemical behavior and characterization. However, as in all theoretical
models, our study has some inherent limitations. The complexity of the system
and the time-scales required to explore different coordination complexes
required the use of empirical force fields rather than ab initio methods. The
finite number of dissolved species in our simulations mimics only low
concentrations, with limited Debye screening, and our force-fields employ
static partial charges which must be scaled appropriately to attempt to
reproduce dielectric screening due to solvent or anion polarizability.

The notable outcome of this study is a reinforcement of the notion that
performant nonaqueous multivalent electrolytes (with high ionic conductivity
and low overpotential) have competing requirements for multivalent ion
coordination: to be strongly coordinating to keep the salt dissolved and the
ionic conductivity high while not so strongly coordinating that the ion cannot
break free from solvation during electrodeposition. The need for strong
solvent coordination is driven by competition with favorable ionic bonds
between counterions. If we want to maintain the advantages of earth-abundant,
multivalent ions, then one option would be to switch to larger anions without
specific coordinating moieties. However, this study highlights that there may
be other options to consider that sideline the isolated multivalent ion
altogether, which may even be irrelevant in terms of electrochemical activity.
Firstly, that coordination with counterions may actually help bring active
species closer to the electrode – that the ability of the species to respect
or disrupt the solvent inhomogeneous structure at the interface can dictate
which species approaches the electrode closest. Secondly, that incomplete
dissolution and the activity of oligomers (dimers in this case) could be key
to improved electrochemical activity albeit balanced by some power limitations
due to reduced bulk ionic conductivity and the need to regenerate
electroactive species through a disproportionation equilibrium. Clearly we
have more work to do, both in terms of understanding the reduction of these
clusters and their dissociation in the reduced state leading to Ca deposition;
the potential negative side-effects of reductive instability of the
coordinating species, already connected to SEI formation;6 and the ultimate
origins of measured overpotentials and currents in experiments. However, we
have highlighted the importance of free energy sampling to attack such complex
problems, even within relatively ideal conditions and contexts, and look
forward to seeing more studies of this kind in the future.

## Methods

Metadynamics sampling (MTD) with a classical force-field was used to obtain
free-energy landscapes as a function of $n$ collective variables. Equilibrium
populations at critical points of a given landscape were then collected with
Umbrella Sampling (US). Structural analysis of the US trajectories was then
performed with a python-based unsupervised learning protocol.

#### Free-energy sampling.

Metadynamics free-energy sampling27 was carried out using the COLVARS module41
implemented in LAMMPS.42 Systems (see Table S1 for a full list) were generated
using Packmol43. Concentration values were chosen to avoid forcing
aggregation. Dimerization free energy surfaces show that at Ca-Ca distances
larger than $\sim$ 10 Å, the free energy converges with respect to Ca-Ca
separation. That is the cutoff we consider between ”dissociated” and
”aggregated”, giving an effective radius of 5 Å for the first coordination
shell of the dication. Additionally, g(r) shows that the second solvation
shell (O-THF) settles at around 6 Å. Assuming an optimal close-packing of
ions, we obtain that SSIPs would be unavoidable at concentrations below 0.5 M
for a 6 Å radius and 0.7 M for a 5 Å radius. Hence, we selected concentrations
$\leq$ 0.03M, well below these limits. System equilibration consisted of
conjugate-gradient minimization to avoid steric clashes, followed by a short
NVT warm-up to room temperature (298 K) using a 1 fs timestep. Box-size
equilibration was achieved by continuing the trajectory under NPT conditions
at 1 atm with a 2 fs timestep. A final NVT step with the equilibrated lattice
parameters ( 20 ns) was sufficient to bring the systems to equilibrium. Force-
field parameters44, 45 were validated in our previous work on the same
system.17 The graphene was frozen in place by setting the forces to zero in
order to ensure neutral and charged simulations were comparable. These MD
setup was kept for the MTD and US simulations. Metadynamics calculations
parameters – the width of the grid along a collective variable (W), height of
the Gaussian ’hills’ used to bias the potential (H) and the frequency at which
they are added (Freq) – can be found in Table S1 and were chosen to ensure
convergence, namely, that the simulation reached diffusive regime in the given
collective variable space.21 Faster completion times were achieved by taking
advantage of multiple-walker metadynamics, which allows to parallelize
sampling among several trajectories (replicas) that update their biased
potential at a given frequency (RepFreq) with the total biased potential.
Despite this, the grid used in the three-dimensional MTD calculations was
necessarily coarser. The minima explored here tend to be separated by more
than 0.7 A (e.g., between the long and short dimer, or between solvent layers,
which is larger than the coarser grid resolution. In order to keep cell
neutrality and consistence with the neutral simulation, charged states 1 and 2
were generated by adding equal and opposite charges ($\pm$ 2 $\mu$C/cm2)
distributed evenly among two graphene layers, emulating a positively charged
and a negatively charged electrode.

Free-energy sampling was performed with combinations of the following
collective variables: the distance between the calcium and the center of mass
of the top graphene layer (dZ); the coordination number between the calcium
and the boron atom in a given BH4 (CN(Ca-B), r0 = 3.8 Å); and, to track
dimerization, the distance between two calcium atoms (dCa) or a Ca-Ca
coordination number (CN(Ca-Ca), r0 = 6.5 Å).

Note that the state free-energies in Fig. 1 b) were obtained by thermodynamic
integration of our 2D free energy surface (Fig. 1 c) over regions delimited by
given collective variable values (e.g., for the Ca(BH4)${}_{4}^{2-}$ state, a
Ca-Ca distance larger that 7 Å and a CN(Ca-B) larger than 3.5). Collective
variables apply to only one of the calcium atoms, and hence the free energy is
averaged over all (remaining) possible coordinations/distances for the other
atom. Hence, the room-temperature Boltzmann probability speaks of the
likelihood of finding a state formed by the constrained species (e.g.,
Ca(BH4)${}_{4}^{2-}$) in the environment of the remaining, unconstrained
species. Since our system contains two calciums and four borohydrides, the un-
constrained space is different depending on the value of the collective
variables. In the Ca(BH4)${}_{4}^{2-}$ basin, the un-constrained Ca can only
be fully solvated. On the other hand, in the Ca${2+}$ basin, the other calcium
can exist in five different ion coordination states (Ca2+, Ca(BH1)+, Ca(BH4)2,
Ca(BH4)${}_{3}^{-}$, and Ca(BH4)${}_{4}^{2-}$). This is the reason why no
symmetry is expected on the free energy surface along the coordination axis,
e.g., between CN=0 and CN=4; and CN=1 and CN=3.

#### Population sampling

The structures and equilibrium population of points of interest in the free
energy surface were collected using Umbrella Sampling. Initial structures were
obtained from metadynamics trajectory snapshots at the desired point in the CV
space. The collective variable coordinates were restrained with a harmonic
potential centered at the desired value, with force constant 1/w2 where w is
the width of the collective variable grid (see Table S1). Then, trajectories
of 150-200 ns were calculated under similar MD parameters as the final
equilibration step and MTD run.

#### Data analysis.

Umbrella sampling trajectories were analyzed with an unsupervised learning
methodology recently developed by us.17 In this protocol, between 5000 and
10000 local atomic arrangements were extracted from each US trajectory, which
when there aligned while taking into account possible permutations between
similar elements (e.g. THF - O).46, 47 Classification based on their
structural similarity was carried out using dimensionality reduction48, 49 and
clustering50 algorithms, in an ASE-compatible51 environment. The
$n$-dimensional free energy landscapes obtained from the MTD sampling were
explored with a Jupyter-adapted version of the MEPSAnd module52 in order to
find critical points and minimum energy pathways.

## Data availability

Additional computational details and free-energy information referenced in the
text can be found in the Supporting Information. The datasets generated and
analysed during the current study are available from the corresponding author
on reasonable request.

## Author Contributions

ASM performed the calculations, analyzed the data and wrote the original draft
with support from DP. ASM and FR developed and tested the clustering
algorithm. SS contributed to force field development. DP supervised and
managed the project. All authors contributed to editing the manuscript.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

This work was fully supported by the Joint Center for Energy Storage Research
(JCESR), an Energy Innovation Hub funded by the U.S. Department of Energy,
Office of Science, Basic Energy Sciences. The theoretical analysis in this
work was supported by a User Project at The Molecular Foundry and its
computing resources, managed by the High Performance Computing Services Group
at Lawrence Berkeley National Laboratory (LBNL), supported by the Director,
Office of Science, Office of Basic Energy Sciences, of the United States
Department of Energy under Contract DE-AC02-05CH11231.

## References

  * Dong  _et al._ [2020] H. Dong, O. Tutusaus, Y. Liang, Y. Zhang, Z. Lebens-Higgins, W. Yang, R. Mohtadi, and Y. Yao, High-power mg batteries enabled by heterogeneous enolization redox chemistry and weakly coordinating electrolytes, Nature Energy 5, 1043 (2020).
  * Mohtadi  _et al._ [2021] R. Mohtadi, O. Tutusaus, T. S. Arthur, Z. Zhao-Karger, and M. Fichtner, The metamorphosis of rechargeable magnesium batteries, Joule 5, 581 (2021).
  * Lv  _et al._ [2022] Y. Lv, Y. Xiao, L. Ma, C. Zhi, and S. Chen, Recent advances in electrolytes for ’beyond aqueous’ zinc-ion batteries, Advanced Materials 34, 2106409 (2022), https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.202106409 .
  * Gu  _et al._ [2017] S. Gu, H. Wang, C. Wu, Y. Bai, H. Li, and F. Wu, Confirming reversible al3+ storage mechanism through intercalation of al3+ into v2o5 nanowires in a rechargeable aluminum battery, Energy Storage Materials 6, 9 (2017).
  * Arroyo-de Dompablo  _et al._ [2020] M. E. Arroyo-de Dompablo, A. Ponrouch, P. Johansson, and M. R. Palacìn, Achievements, challenges, and prospects of calcium batteries, Chemical Reviews 120, 6331 (2020), doi: 10.1021/acs.chemrev.9b00339.
  * McClary  _et al._ [2022] S. A. McClary, D. M. Long, A. Sanz-Matias, P. G. Kotula, D. Prendergast, K. L. Jungjohann, and K. R. Zavadil, A heterogeneous oxide enables reversible calcium electrodeposition for a calcium battery, ACS Energy Letters 7, 2792 (2022), https://doi.org/10.1021/acsenergylett.2c01443 .
  * Liang  _et al._ [2020] Y. Liang, H. Dong, D. Aurbach, and Y. Yao, Current status and future directions of multivalent metal-ion batteries, Nature Energy 5, 646 (2020).
  * Deng  _et al._ [2022] X. Deng, L. Li, G. Zhang, X. Zhao, J. Hao, C. Han, and B. Li, Anode chemistry in calcium ion batteries: A review, Energy Storage Materials 53, 467 (2022).
  * Lee  _et al._ [2023] E.-S. Lee, S.-H. Huh, S.-H. Lee, and S.-H. Yu, On the road to stable electrochemical metal deposition in multivalent batteries, ACS Sustainable Chemistry & Engineering 11, 2014 (2023), https://doi.org/10.1021/acssuschemeng.2c05173 .
  * Forero-Saboya  _et al._ [2022] J. D. Forero-Saboya, D. S. Tchitchekova, P. Johansson, M. R. Palacìn, and A. Ponrouch, Interfaces and interphases in ca and mg batteries, Advanced Materials Interfaces 9, 2101578 (2022), https://onlinelibrary.wiley.com/doi/pdf/10.1002/admi.202101578 .
  * Roy  _et al._ [2016] S. Roy, M. D. Baer, C. J. Mundy, and G. K. Schenter, Reaction rate theory in coordination number space: An application to ion solvation, The Journal of Physical Chemistry C 120, 7597 (2016), doi: 10.1021/acs.jpcc.6b00443.
  * Baer and Mundy [2016] M. D. Baer and C. J. Mundy, Local aqueous solvation structure around ca2+ during ca2+-cl- pair formation, The Journal of Physical Chemistry B 120, 1885 (2016), doi: 10.1021/acs.jpcb.5b09579.
  * Joswiak  _et al._ [2018] M. N. Joswiak, M. F. Doherty, and B. Peters, Ion dissolution mechanism and kinetics at kink sites on nacl surfaces, Proceedings of the National Academy of Sciences 115, 656 (2018), https://www.pnas.org/doi/pdf/10.1073/pnas.1713452115 .
  * Baskin and Prendergast [2019] A. Baskin and D. Prendergast, ‘ion solvation spectra’: Free energy analysis of solvation structures of multivalent cations in aprotic solvents, The Journal of Physical Chemistry Letters 10, 4920 (2019), doi: 10.1021/acs.jpclett.9b01569.
  * Baskin and Prendergast [2020] A. Baskin and D. Prendergast, Ion solvation engineering: How to manipulate the multiplicity of the coordination environment of multivalent ions, The Journal of Physical Chemistry Letters 11, 9336 (2020), pMID: 33090799, https://doi.org/10.1021/acs.jpclett.0c02682 .
  * Silvestri  _et al._ [2022] A. Silvestri, P. Raiteri, and J. D. Gale, Obtaining consistent free energies for ion binding at surfaces from solution: Pathways versus alchemy for determining kink site stability, Journal of Chemical Theory and Computation 18, 5901 (2022), doi: 10.1021/acs.jctc.2c00787.
  * Roncoroni  _et al._ [2023] F. Roncoroni, A. S. Matias, S. Sundararaman, and D. Prendergast, Unsupervised learning of representative local atomic arrangements in molecular dynamics data,  (2023), arXiv:2302.01465 [cond-mat] .
  * Wang  _et al._ [2018] D. Wang, X. Gao, Y. Chen, L. Jin, C. Kuss, and P. G. Bruce, Plating and stripping calcium in an organic electrolyte, Nature Materials 17, 16 (2018).
  * Hahn  _et al._ [2020] N. T. Hahn, J. Self, T. J. Seguin, D. M. Driscoll, M. A. Rodriguez, M. Balasubramanian, K. A. Persson, and K. R. Zavadil, The critical role of configurational flexibility in facilitating reversible reactive metal deposition from borohydride solutions, J. Mater. Chem. A 8, 7235 (2020).
  * Bodin  _et al._ [2023] C. Bodin, J. Forero Saboya, P. Jankowski, K. Radan, D. Foix, C. Courrèges, I. Yousef, R. Dedryvère, C. Davoisne, M. Lozinšek, and A. Ponrouch, Boron-based functional additives enable solid electrolyte interphase engineering in calcium metal battery, Batteries & Supercaps 6, e202200433 (2023), https://chemistry-europe.onlinelibrary.wiley.com/doi/pdf/10.1002/batt.202200433 .
  * Baskin  _et al._ [2021] A. Baskin, J. W. Lawson, and D. Prendergast, Anion-assisted delivery of multivalent cations to inert electrodes, The Journal of Physical Chemistry Letters 12, 4347 (2021), pMID: 33929859, https://doi.org/10.1021/acs.jpclett.1c00943 .
  * Ta  _et al._ [2019] K. Ta, R. Zhang, M. Shin, R. T. Rooney, E. K. Neumann, and A. A. Gewirth, Understanding ca electrodeposition and speciation processes in nonaqueous electrolytes for next-generation ca-ion batteries, ACS Applied Materials & Interfaces 11, 21536 (2019), pMID: 31117456, https://doi.org/10.1021/acsami.9b04926 .
  * Melemed  _et al._ [2020] A. M. Melemed, A. Khurram, and B. M. Gallant, Current understanding of nonaqueous electrolytes for calcium-based batteries, Batteries & Supercaps 3, 570 (2020), https://chemistry-europe.onlinelibrary.wiley.com/doi/pdf/10.1002/batt.201900219 .
  * Jie  _et al._ [2020] Y. Jie, Y. Tan, L. Li, Y. Han, S. Xu, Z. Zhao, R. Cao, X. Ren, F. Huang, Z. Lei, G. Tao, G. Zhang, and S. Jiao, Electrolyte solvation manipulation enables unprecedented room-temperature calcium-metal batteries, Angewandte Chemie International Edition 59, 12689 (2020), https://onlinelibrary.wiley.com/doi/pdf/10.1002/anie.202002274 .
  * Liepinya and Smeu [2021] D. Liepinya and M. Smeu, A computational comparison of ether and ester electrolyte stability on a ca metal anode, Energy Material Advances 2021, 9769347 (2021).
  * Lu  _et al._ [2021] Y.-T. Lu, A. R. Neale, C.-C. Hu, and L. J. Hardwick, Divalent nonaqueous metal-air batteries, Frontiers in Energy Research 8, 357 (2021).
  * Laio and Parrinello [2002] A. Laio and M. Parrinello, Escaping free-energy minima, Proceedings of the National Academy of Sciences 99, 12562 (2002), https://www.pnas.org/doi/pdf/10.1073/pnas.202427399 .
  * Hahn  _et al._ [2021] N. T. Hahn, J. Self, K. S. Han, V. Murugesan, K. T. Mueller, K. A. Persson, and K. R. Zavadil, Quantifying species populations in multivalent borohydride electrolytes, The Journal of Physical Chemistry B 125, 3644 (2021), doi: 10.1021/acs.jpcb.1c00263.
  * Henderson  _et al._ [1976] D. Henderson, F. F. Abraham, and J. A. Barker, The ornstein-zernike equation for a fluid in contact with a surface, Molecular Physics 31, 1291 (1976), https://doi.org/10.1080/00268977600101021 .
  * Feng  _et al._ [2010] G. Feng, J. Huang, B. G. Sumpter, V. Meunier, and R. Qiao, Structure and dynamics of electrical double layers in organic electrolytes, Phys. Chem. Chem. Phys. 12, 5468 (2010).
  * Butt  _et al._ [2006] H.-J. Butt, K. Graf, and M. Kappl,  _Physics and Chemistry of Interfaces_ (Wiley-VCH, Germany, 2006) pp. 45,44,56,76–83.
  * Lu  _et al._ [2019] Y.-H. Lu, J. M. Larson, A. Baskin, X. Zhao, P. D. Ashby, D. Prendergast, H. A. Bechtel, R. Kostecki, and M. Salmeron, Infrared nanospectroscopy at the graphene-electrolyte interface, Nano Letters 19, 5388 (2019), doi: 10.1021/acs.nanolett.9b01897.
  * Yang  _et al._ [2022] S. Yang, X. Zhao, Y.-H. Lu, E. S. Barnard, P. Yang, A. Baskin, J. W. Lawson, D. Prendergast, and M. Salmeron, Nature of the electrical double layer on suspended graphene electrodes, Journal of the American Chemical Society 144, 13327 (2022), doi: 10.1021/jacs.2c03344.
  * He  _et al._ [2022] X. He, J. M. Larson, H. A. Bechtel, and R. Kostecki, In situ infrared nanospectroscopy of the local processes at the li/polymer electrolyte interface, Nature Communications 13, 1398 (2022).
  * Velasco-Velez  _et al._ [2014] J.-J. Velasco-Velez, C. H. Wu, T. A. Pascal, L. F. Wan, J. Guo, D. Prendergast, and M. Salmeron, The structure of interfacial water on gold electrodes studied by x-ray absorption spectroscopy, Science 346, 831 (2014), https://www.science.org/doi/pdf/10.1126/science.1259437 .
  * Wu  _et al._ [2018] C. H. Wu, T. A. Pascal, A. Baskin, H. Wang, H.-T. Fang, Y.-S. Liu, Y.-H. Lu, J. Guo, D. Prendergast, and M. B. Salmeron, Molecular-scale structure of electrode-electrolyte interfaces: The case of platinum in aqueous sulfuric acid, Journal of the American Chemical Society 140, 16237 (2018), doi: 10.1021/jacs.8b09743.
  * Prabhakaran  _et al._ [2023] V. Prabhakaran, G. Agarwal, J. D. Howard, S. Wi, V. Shutthanandan, D.-T. Nguyen, L. Soule, G. E. Johnson, Y.-S. Liu, F. Yang, X. Feng, J. Guo, K. Hankins, L. A. Curtiss, K. T. Mueller, R. S. Assary, and V. Murugesan, Coordination-dependent chemical reactivity of tfsi anions at a mg metal interface, ACS Applied Materials & Interfaces 15, 7518 (2023), doi: 10.1021/acsami.2c18477.
  * Yu  _et al._ [2017] Y. Yu, A. Baskin, C. Valero-Vidal, N. T. Hahn, Q. Liu, K. R. Zavadil, B. W. Eichhorn, D. Prendergast, and E. J. Crumlin, Instability at the electrode/electrolyte interface induced by hard cation chelation and nucleophilic attack, Chemistry of Materials 29, 8504 (2017), doi: 10.1021/acs.chemmater.7b03404.
  * Hahn  _et al._ [2022] N. T. Hahn, S. A. McClary, A. T. Landers, and K. R. Zavadil, Efficacy of stabilizing calcium battery electrolytes through salt-directed coordination change, The Journal of Physical Chemistry C 126, 10335 (2022), doi: 10.1021/acs.jpcc.2c02587.
  * Jay  _et al._ [2019] R. Jay, A. W. Tomich, J. Zhang, Y. Zhao, A. De Gorostiza, V. Lavallo, and J. Guo, Comparative study of mg(cb11h12)2 and mg(tfsi)2 at the magnesium/electrolyte interface, ACS Applied Materials & Interfaces 11, 11414 (2019), doi: 10.1021/acsami.9b00037.
  * Fiorin  _et al._ [2013] G. Fiorin, M. L. Klein, and J. Hénin, Using collective variables to drive molecular dynamics simulations, Molecular Physics 111, 3345 (2013), https://doi.org/10.1080/00268976.2013.813594 .
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Lammps - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Computer Physics Communications 271, 108171 (2022).
  * Martínez  _et al._ [2009] L. Martínez, R. Andrade, E. G. Birgin, and J. M. Martínez, Packmol: A package for building initial configurations for molecular dynamics simulations, Journal of Computational Chemistry 30, 2157 (2009), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.21224 .
  * Sambasivarao and Acevedo [2009] S. V. Sambasivarao and O. Acevedo, Development of opls-aa force field parameters for 68 unique ionic liquids, Journal of Chemical Theory and Computation 5, 1038 (2009), pMID: 26609613, https://doi.org/10.1021/ct900009a .
  * de Andrade  _et al._ [2002] J. de Andrade, E. S. Böes, and H. Stassen, Computational study of room temperature molten salts composed by 1-alkyl-3-methylimidazolium cationsforce-field proposal and validation, The Journal of Physical Chemistry B 106, 13344 (2002), https://doi.org/10.1021/jp0216629 .
  * Gunde  _et al._ [2022] M. Gunde, N. Salles, A. Hémeryck, and L. Martin-Samos, Iterative rotations and assignments (ira): A shape matching algorithm for atomic structures, Software Impacts 12, 100264 (2022).
  * Griffiths  _et al._ [2017] M. Griffiths, S. P. Niblett, and D. J. Wales, Optimal alignment of structures for finite and periodic systems, Journal of Chemical Theory and Computation 13, 4914 (2017), pMID: 28841314, https://doi.org/10.1021/acs.jctc.7b00543 .
  * McInnes  _et al._ [2018] L. McInnes, J. Healy, and J. Melville, Umap: Uniform manifold approximation and projection for dimension reduction (2018).
  * Becht  _et al._ [2019] E. Becht, L. McInnes, J. Healy, C.-A. Dutertre, I. W. H. Kwok, L. G. Ng, F. Ginhoux, and E. W. Newell, Dimensionality reduction for visualizing single-cell data using umap, Nature Biotechnology 37, 38 (2019).
  * McInnes  _et al._ [2017] L. McInnes, J. Healy, and S. Astels, hdbscan: Hierarchical density based clustering, Journal of Open Source Software 2, 205 (2017).
  * Larsen  _et al._ [2017] A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacobsen, The atomic simulation environment—a python library for working with atoms, Journal of Physics: Condensed Matter 29, 273002 (2017).
  * Marcos-Alcalde  _et al._ [2019] I. Marcos-Alcalde, E. López-Viñas, and P. Gómez-Puertas, MEPSAnd: minimum energy path surface analysis over n-dimensional surfaces, Bioinformatics 36, 956 (2019), https://academic.oup.com/bioinformatics/article-pdf/36/3/956/32369715/btz649.pdf .
  * Ufimtsev and MartÃ?nez [2008] I. S. Ufimtsev and T. J. MartÃ?nez, Quantum chemistry on graphical processing units. 1. strategies for two-electron integral evaluation, Journal of Chemical Theory and Computation 4, 222 (2008), doi: 10.1021/ct700268q.
  * Stephens  _et al._ [1994] P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields, The Journal of Physical Chemistry 98, 11623 (1994), https://doi.org/10.1021/j100096a001 .
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements h-pu, The Journal of Chemical Physics 132, 154104 (2010), https://doi.org/10.1063/1.3382344 .
  * Soler  _et al._ [2002] J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, The siesta method for ab initio order-n materials simulation, Journal of Physics: Condensed Matter 14, 2745 (2002).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Hamada  _et al._ [2013] I. Hamada, O. Sugino, N. Bonnet, and M. Otani, Improved modeling of electrified interfaces using the effective screening medium method, Phys. Rev. B 88, 155427 (2013).
  * Otani and Sugino [2006] M. Otani and O. Sugino, First-principles calculations of charged surfaces and interfaces: A plane-wave nonrepeated slab approach, Phys. Rev. B 73, 115407 (2006).

Supplemental Materials: Ca-dimers and solvent layering determine
electrochemically active species in Ca(BH4)2 in THF

Figure S1: Bulk free-energy surface of dimerization from ionic and neutral
species (a). Bulk disproportionation (b) versus direct dissociation (c)
pathways.  Figure S2:  Determination of short and long dimers equilibrium RT
distances. Free-energy of dimerization (distance Ca-Ca) as a function of the
Ca-B(H4) coordination number (left) and Ca-O(THF) coordination number (r0 =
3.6 Å, right). One-dimensional profiles obtained through thermodynamic
integration over each coordination number allow us to interpolate the
equilibrium distances of the short and long dimer.  Figure S3:  Solvent
coordination of perpendicular dimers. Orientation of the Ca-Ca axis with
respect to the surface normal of structures of the isomer SD 5 THF - 4 at the
DL (see Fig. 3 in the main text) depending on the THF-Ca coordination of the
Ca that resides in the DL. At strongly perpendicular orientations (larger that
150 ∘ with respect to the surface normal), two-fold solvent coordination is
favored.  Figure S4:  Free-energy of BH${}_{4}^{-}$ at the (biased) interface.
The free-energy profile of a single THF molecule in THF (top) at RT with
respect to distance from a graphite interface (red) and the corresponding
oxygen density profile (black), with slight adjustment (dashed lines) at
negative bias. The Dense Layer (DL), Intermediate Density Layer (IDL) and bulk
regions are color-coded hereafter in red, blue and gray. The bottom figure
shows the analogous free-energy profile of a single BH${}_{4}^{-}$ molecule in
THF at RT with respect to distance from a graphite interface (red) and the
corresponding boron density profile (black).

##  Supplementary Information 1 Dipole moments of Ca(BH4)2 isomers

In order to compute the dipole moment of specific Ca(BH4)2 isomer, an average
structure of the whole group was calculated (including THF solvating
molecules). Then, using Density Functional Theory as implemented in
Terachem,53 the atomic positions were optimized to obtain the dipole moment.
Typically, a decrease in the dipole moment of about 1 D is observed upon
minimization. The B3LYP54 exchange-correlation functional with a dispersion
correction (D3)55 and a 6-311G basis set were used.

##  Supplementary Information 2 Fermi energy changes for charged surfaces

Estimations of Fermi energy changes with electronic charge were carried out on
a model graphene supercell of dimensions 0.426 x 0.246 nm to mimic the
surface. Electronic structure calculations using Density functional theory
(DFT) as implemented in SIESTA56, with the PBE functional57 and a double-zeta
with polarization (DZP) basis set. In addition, orbitals-confinement was
reduced by using a PAO Energy shift of 0.001 Ry in order to capture the
electron density decay at the surface. A 10x10x1 k-point grid together with a
400 Ry charge density grid cutoff were used to ensure numerical convergence.
Additionally, the Effective Screening Medium (ESM)58, 59 model was used to
ensure a common reference vacuum level. In the ESM, the slab is placed between
two semi-inifinite mediums that can be vacuum or metal-like (to mimic an
electrode or electrolyte). We have compared two screening setups: one in which
the slab is sandwiched between metal mediums (known as bc2) or between a metal
and a vacuum media (bc3). Both cases yield the same result: a change in
workfunction of 0.56 and 0.74 eV for a surface charge of 0.065 and 0.13 e-
nm-2, respectively.

(a)

(b)

Figure S5:  Fermi energies of charged graphene layers. Aligned total density of states (DOS) of graphene with neutral and negative charge, with each Fermi level energy marked by a vertical line (a). The corresponding ESM potential as a function of the distance from the graphene (at z=0 Å) in the Metal $|$ slab $|$ Metal ESM mode (b).  Fig. | THF | Ca2+ | BH${}_{4}^{-}$ | G | Box pars. (Å) | CVs | qG (e-) | W | H | Freq | Reps. | RepFreq | Time (ns)
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Bulk | 200 | 2 | 4 | - | 30.3 x 30.3 x30.3 | dCa [0:16] Å | - | 0.05 | 0.02 | 200 | 20 | 50000 | 1660
1 |  |  |  |  |  | CN(Ca-B) [0-4] |  | 0.05 |  |  |  |  |
Neutral | 336 | 2 | 4 | 2 | 33.9 x 31.8 x 49.8 | dCa [0:16] Å | G1 = G2 = 0 | 0.2 | 0.02 | 200 | 70 | 50000 | 2082
2,3, |  |  |  |  |  | CN(Ca-B) [0-4] |  | 0.2 |  |  |  |  |
5 |  |  |  |  |  | dZ [0:25] |  | 0.5 |  |  |  |  |
2 | 336 | - | - | 2 | 33.9 x 31.8 x 49.8 | dZ [0:25] Å | G1 = G2 = 0 | 0.05 | 0.02 | 200 | 15 | 50000 | 1049
2 | 336 | - | - | 2 | 33.9 x 31.8 x 49.8 | dZ [0:25] | G1 = -1.4; G2 = 0 | 0.05 | 0.02 | 200 | 10 | 50000 | 620
CS1 | 514 | 2 | 4 | 2 | 33.9 x 31.8 x 71.8 | CN(Ca-Ca) [0:1] Å | G1 = -0.7; | 0.2 | 0.04 | 200 | 70 | 50000 | 1438
|  |  |  |  |  | CN(Ca-B) [0-4] | G2 = +0.7 | 0.2 |  |  |  |  |
|  |  |  |  |  | dZ [2:70] |  | 0.25 |  |  |  |  |
CS2 | 514 | 2 | 4 | 2 | 33.9 x 31.8 x 71.8 | CN(Ca-Ca) [0:1] Å | G1 = -1.4; | 0.2 | 0.04 | 200 | 70 | 50000 | 1464
2,4, |  |  |  |  |  | CN(Ca-B) [0-4] | G2 = +1.4 | 0.2 |  |  |  |  |
5 |  |  |  |  |  | dZ [2:70] |  | 0.25 |  |  |  |  |
Table S1:  Metadynamics simulation details. Figure number or label, box
composition (Ca2+, BH${}_{4}^{-}$ and graphene) and equilibrated lattice
parameters, collective variables and their ranges, surface charges, and
metadynamics parameters (hill width and height, frequency, number of
replicals, replica update frequency, and total simulation time).  Figure S6:
Species free-energy and populations in charged state 1. Populations of
electrolyte species derived from their respective free energies $\Delta$G
(kT), which were obtained by integration of the free-energy surface sampled
using metadynamics with respect to Ca-surface distance, Ca-Ca distance and
Ca-B(H4) coordination number with a surface charge of 0.065 e nm-2. At this
charge, the dominating species at the IDL are still dimers, but at the DL, the
neutral monomer Ca(BH4)2 becomes more favorable.  Figure S7:  Interfacial
free-energy of [1,4] dimers. Free energy pathways as a function of distance
from the neutral surface for dimers with [1,4] Ca-B(H4)- coordination numbers
at different Ca-Ca distances. The most favorable, with d(Ca-Ca) = 4.9-5.1 Å,
preferably sits far away from the layered interface, at $\sim$ 17.5 Å.