# Joint operation of a fast-charging EV hub with a stand-alone independent
battery storage system under fairness considerations

Diwas Paudel, Luke Wolf, and Tapas K. Das
Department of Industrial and Management Systems Engineering,
University of South Florida,
Tampa, FL
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

The need for larger-scale fast-charging electric vehicle (EV) hubs is on the
rise due to the growth in EV adoption. Another area of power infrastructure
growth is the proliferation of independently operated stand-alone battery
storage systems (BSS), which is fueled by improvements and cost reductions in
battery technology. Many possible uses of the stand-alone BSS are being
explored including participation in the energy and ancillary markets, load
balancing for renewable generations, and supporting large-scale load-consuming
entities like hospitals. In this paper, we study a novel usage of the stand-
alone BSS whereby in addition to participating in the electricity reserve
market, it allows an EV hub to use a part of its storage capacity, when
profitable. The hub uses the BSS storage capacity for arbitrage consequently
reducing its operating cost. We formulate this joint operation as a bi-
objective optimization model. We then reformulate it into a second-order cone
Nash bargaining problem, the solution of which guarantees fairness to both the
hub and the BSS. A sample numerical case study is formulated using actual
prices of electricity and simulated data for the reserve market and EV
charging demand. The Nash bargaining solution shows that both participants can
benefit from the joint operation.

_Keywords_ EV charging $\cdot$ battery storage system $\cdot$ bi-objective
$\cdot$ Nash bargain $\cdot$ Second-order cone problem (SOCP)

##  1 Introduction

Rapid growth of electric vehicles (EVs) has prompted a critical need for an
extensive infrastructure of EV charging facilities Conzade et al. (2021). A
substantial portion of these facilities is anticipated to be hubs housing
several (often many) fast charging stations. The growth in EV charging
infrastructure will play a pivotal role in minimizing range anxiety among EV
owners further promoting EV adoption. The existing Tesla Supercharging
stations are examples of fast-charging hubs Tesla . The charging hubs can
source the needed power in multiple ways. A primary source of power is the
day-ahead (DA) electricity market where hub operators can commit to buying a
certain amount of the needed power for every hour of the following day. The
real-time (RT) electricity market serves as one of the secondary sources to
supply any additional power requirement beyond the day-ahead commitment. The
other secondary source of power can be a battery storage system (BSS), which
can be built-in or external to the hub. In addition to serving as a secondary
source of power, BSS can allow the hub operators to use the storage capacity
for arbitrage (i.e., to store power at cheaper prices and discharge when
prices are high).

With the advancement of battery technology, large-scale stand-alone battery
storage systems are becoming commonplace. Such BSS can partake in supporting
one or more of the entities including the electric grid (by supplying power
for the reserve/ancillary markets for load balancing and frequency
regulation), renewable electricity-generating farms, and a variety of large
load-consuming entities like hospitals, factories, etc. As the cost of
batteries continues to fall, large stand-alone BSS is becoming a more
mainstream element of our newly built power infrastructure NREL (2023). In the
rest of the paper, we use the word hub to refer to the fast-charging EV hub
and the word BSS to refer to the stand-alone battery storage system.

EV charging hubs with numerous stations will soon join the ranks of
significant load-consuming centers and thus will benefit from working with
independently operated stand-alone BSS. In this paper, we consider a joint
operation between a fast-charging EV hub and a stand-alone BSS. The hub uses
BSS capacity to arbitrage by storing power from the DA and RT markets and
discharging power to either charge EVs and/or sell power back to the RT
market. The BSS is considered to participate in the reserve market, by
submitting reserve up and reserve down bids, and also to allow the hub
operator to use a part of its storage capacity. It is considered that the
objective of the hub operator is to reduce its cost of power in meeting the EV
charging demand for the day by optimally sourcing DA and RT power and deciding
on the arbitrage strategy. The BSS, on the other hand, aims to maximize its
gross profit by optimally using the storage capacity among its two
participating activities. To our knowledge, a study examining joint operation
between a fast-charging hub and a stand-alone multipurpose BSS, as explained
above, has not been presented in the literature. We undertake such a study by
developing a bi-objective optimization model. We then reformulate the model to
obtain a Nash bargaining solution that is fair to both participants. The model
is implemented on a sample case study problem built with real-world electric
power market data and simulated data for EV charging demand and reserve market
operations. The results demonstrate the value of cooperation and thus will
help to establish new ventures in cooperative operation between hubs and BSS.

The rest of this paper is organized as follows. A brief review of the related
literature is provided in Section 2. We describe the problem in Sections 3 and
provide the mathematical formulation in Section 4. The details of our solution
approach are provided in Section 5 followed by a numerical case study in
Section 6. Finally, the concluding remarks are provided in Section 7.

##  2 Related literature

In this section, we present a discussion of the relevant literature that
explores the following aspects of our problem: power management of EV charging
hubs and the operation of the stand-alone BSS. Based on the review, we
highlight the gap in the literature that our paper aims to address.

The power management of EV charging hubs is a widely studied topic in the open
literature. These studies primarily consider using electricity from DA and RT
markets, while some also consider the use of in-house BSS, to meet the EV
charging demands. The study presented in Zheng et al. (2020) approaches hub
power management via a profit-maximizing day-ahead power procurement strategy
accounting for RT market volatility and no BSS availability. The strategy is
derived using a scenario-based stochastic optimization approach. A
distributionally robust DA commitment strategy for hub power management is
presented in Paudel et al. (2023), which in addition to the DA and RT power
considers the in-house BSS availability. A system state-dependent real-time
power management strategy is obtained in Paudel and Das (2023) via a deep-
reinforcement learning model that considers both DA and RT markets as well as
in-house BSS. There exists a significant collection of literature that
considers EV charging hubs which are also parking garages where EVs remain
parked for a significant length of time while also receiving charge. These
papers consider effectively scheduling the charging time of each EV depending
on their arrival and departure schedules and market price variations to reduce
the cost of charging while also supporting network load balancing Subramanian
and Das (2019). Some of these papers also consider arbitrage using available
power in the parked EVs by participating in a variety of power markets
including the ancillary market DeForest et al. (2018) Sarker et al. (2015). To
our knowledge, there are no hub power management studies that consider working
jointly with independently operated stand-alone battery storage systems, which
are becoming commonplace in modern-day power infrastructure. In the following
paragraph, we review some of the papers that explore some of the potential
services that stand-alone BSS can provide to the power grid.

The literature on the operation of stand-alone BSS primarily focuses on their
participation in the reserve/frequency markets. A robust formulation for the
operation of a BSS participating in spinning and reserve markets is presented
in Kazemi et al. (2017). The BSS generates revenue through accepted bids,
received as capacity payments, and also through the deployment of accepted
bids at the prevailing real-time market prices. A similar study considering
the reserve market participation by stand-alone BSS is also provided in
Padmanabhan et al. (2019) and Xu et al. (2017). A bidding strategy for battery
storage systems in the secondary control reserve market is examined in Merten
et al. (2020). The study indicates that the feasibility of revenue generation
through reserve market participation by the BSS in stand-alone mode is
limited. However, it becomes viable when operating jointly with a virtual
power plant that includes wind, photovoltaic, and thermal power sources. An
economic study conducted in Germany explored bidding on both day-ahead and
automatic frequency restoration reserve markets by a stand-alone BSS Nitsch et
al. (2021). The study shows an improved economic potential for stand-alone BSS
in 2030 than it was in 2019. In a similar economic study of BSS Hu et al.
(2022), the authors indicate that present BSS cost and lifetime make it not
suitable for the energy market participation but has significant economic
potential in the reserve markets.

Based on the reviewed papers, contemporary Battery Storage Systems (BSS)
demonstrate the ability to offer cost-effective services to electricity
markets, especially in the reserve market. However, there’s a noticeable gap
in studies assessing the synergy between a stand-alone BSS participating in
the reserve market and fast-charging EV hubs. Therefore, our study aims to
address this gap by investigating the joint operation of BSS with a fast-
charging hub. Furthermore, our research introduces a novel pathway for
expanding the range of services a BSS can provide through its collaboration
with a fast-charging EV hub.

##  3 Problem Description

Figure 1: A schematic for joint operation of fast charging EV hub and battery
storage system

We consider a hub with a large number of identical DC-fast charging stations
where EVs arrive only for charging and leave promptly after charging. The hub
draws power from the day-ahead and real-time markets of the electricity grid
to meet EV charging needs. The primary function of the stand-alone independent
BSS is considered to be supporting the reserve market by offering the network
operator the use of its storage capacity via reserve up and down services
bids. The reserve up capacity represents the amount of stored power that the
BSS is willing to supply in real-time in case of an unexpected shortfall in
the network. The reserve down capacity, on the other hand, represents the
amount of storage capacity that the BSS is willing to offer for storing any
excess power in the network. The acceptance of up and/or down bid offers does
not imply the actual deployment, i.e., supply or storage of power to/from the
grid by the BSS. The deployment occurs sparingly, the probability of which can
be accessed from the historical market data.

The focus of this paper is to study the joint operation of a hub with a BSS
that is connected to the same bus of power network for their mutual benefit.
In the joint operation mode, the hub uses part of the BSS storage capacity to
support its charging needs through arbitrage, i.e. by storing power at a low
cost and discharging it when cost-economic, thereby reducing the cost of its
operation The BSS benefits from the additional revenue generated by allowing
the hubs to utilize part of its capacity when not deployed for the reserve
purposes. The total storage capacity is comprised of multiple compartments
each of which can be operated independently by the BSS operator. The key
problem addressed in this paper is to determine the optimal joint operational
strategies for the hub and the BSS that yield a fair distribution of the added
benefits generated from cooperation.

##  4 Mathematical model

In this section, we first present two separate models for the independent
operation of the hub and the battery storage system. The model P1 is focused
on minimizing the hub operator’s cost of meeting the charging demand by
procuring an optimal amount of power from the day-ahead market for each time
interval depending on the expected price volatilities of the real-time market.
The model P2 maximizes the profit of the BSS operator by optimally selecting
up and down bid quantities in the reserve market for each time interval.
Thereafter, we present a joint operation model P3 where in addition to
participating in the reserve market, for each time period, the BSS operator
also leases to the hub operator a part of the storage capacity. The hub
operator uses the leased capacity for arbitrage and makes a payment for each
unit of the stored/discharged quantities. The bi-objective model P3 allows
both the hub and the BSS operators to improve their objectives through
cooperation.

###  4.1 P1: Fast-charging EV hub power management model

In this model, it is assumed that the hub operator acts independently of BSS
and aims to minimize its cost of charging by optimally choosing to use power
from the DA and RT markets based on their expected prices. The model
determines the optimal hourly quantities for the hub operator to commit in the
DA market.

$da_{t}^{com}=da_{t}^{ev}+da_{t}^{rt},\quad\forall t\in\mathcal{T}.$ (1a)
$da_{t}^{com}\leq U_{t}^{com},\quad\forall t\in\mathcal{T}.$ (1b)
$ev_{t}^{l}=da_{t}^{ev}+rt_{t}^{ev},\quad\forall t\in\mathcal{T}.$ (1c) $\min
f_{p_{1}}^{a}=\min\sum_{t\in\mathcal{T}}\Bigg{[}\lambda_{t}^{da}da_{t}^{ev}+(\lambda_{t}^{da}-\lambda_{t}^{rt})da_{t}^{rt}+\lambda_{t}^{rt}rt_{t}^{ev}\Bigg{]}$
(1d)

Among the constraints of P1, (1a) ensures that the DA commitment
$(da_{t}^{com})$ by the hub operator is fully consumed by charging the EVs in
the hub $(da_{t}^{ev})$ and selling in real-time market $(da_{t}^{rt})$. The
upper bound for the day-ahead commitment quantities is provided by (1b).
Similarly, (1c) ensures that the EV charging demand in the hub $(ev_{t}^{l})$
is fulfilled by using the power from the DA commitment and RT market
$(rt_{t}^{ev})$. Note that we consider aggregated EV charging demand for each
hub. The goal of the hub operator is to satisfy the total charging demand in
the hub in the most cost-effective way as represented by the objective
function in (1d), where $\lambda_{t}^{da}$ and $\lambda_{t}^{rt}$ are the
prices of DA and RT power, respectively. The first term in (1d) represents the
total cost of using DA power to charge the EVs. The cost of selling DA power
in the RT market is captured by the second term. The last term represents the
cost of charging the EVs with RT power. Problem P1 is represented in a compact
form as follows.

$\displaystyle\begin{split}\min_{\bm{x}}f_{p_{1}}^{a}(\bm{x}),\\\
\bm{x}\in\bm{X_{f}},\end{split}$ (2)

where, $\bm{x}$ and $\bm{X_{f}}$ represent the decision variables and the
decision space, respectively, of P1. This compact notation is used in the
description of the model solution approach in Section 5.

###  4.2 P2: BSS model for reserve market participation

This model assists an independent BSS operator, participating only in the
reserve market, in determining the hourly bids for reserve-up and reserve-down
capacities to maximize profit. The reserve-up capacity denotes the power
quantity that the BSS is ready to deliver in real-time to address an
unforeseen surge in network demand. Conversely, the reserve-down capacity
signifies the power amount that the BSS is prepared to absorb (store) in
instances of surplus supply within the network. It is considered that the BSS
comprises multiple compartments and within each compartment, there is an
independent battery unit. These battery units are capable of independent
operation, allowing for greater flexibility and control. Each compartment can
either charge, discharge, or idle at any time. The total reserve-up and down
bids made by the BSS operator are aggregations of compartment bids, which are
bounded by the maximum discharge/charge rates of the respective compartments.
The up and down bids are assumed to be accepted and deployed with given
probabilities, as is considered in Sarker et al. (2015) DeForest et al.
(2018).

$x_{t}^{k}+y_{t}^{k}\leq 1,\quad\forall k\in\mathcal{K},\forall
t\in\mathcal{T}.$ (3a) $\widehat{p}_{t}^{up,k}<=M^{d,k}y_{t}^{k},\quad\forall
k\in\mathcal{K},t\in\mathcal{T}.$ (3b)
$\widehat{p}_{t}^{dn,k}<=M^{c,k}x_{t}^{k},\quad\forall
k\in\mathcal{K},t\in\mathcal{T}.$ (3c)
$\sum_{k\in\mathcal{K}}p_{t}^{up,k}=\pi_{t}^{d,up}\sum_{k\in\mathcal{K}}\widehat{p}_{t}^{up,k},\quad\forall
t\in\mathcal{T}.$ (3d)
$\sum_{k\in\mathcal{K}}p_{t}^{dn,k}=\pi_{t}^{d,dn}\sum_{k\in\mathcal{K}}\widehat{p}_{t}^{dn,k},\quad\forall
t\in\mathcal{T}.$ (3e) $p_{t}^{dn,k}+rt_{t}^{b,k}\leq
M^{c,k}x_{t}^{k},\quad\forall k\in\mathcal{K},\forall t\in\mathcal{T}.$ (3f)
$p_{t}^{up,k}\leq M^{d,k}y_{t}^{k},\quad\forall k\in\mathcal{K},\forall
t\in\mathcal{T}.$ (3g) $\phi_{t}^{b,k}\leq M^{k},\quad\forall
k\in\mathcal{K},\forall t\in\mathcal{T}.$ (3h)
$\phi_{t}^{b,k}=\phi_{t-1}^{b,k}+p_{t}^{dn,k}+rt_{t}^{b,k}-p_{t}^{up,k},\quad\forall
k\in\mathcal{K},\forall t\in\mathcal{T}.$ (3i) $\max
f_{p_{2}}^{b}=\max\bigg{[}r^{cap}+r^{dep}-c^{\phi}-c^{deg}\bigg{]}$ (3j)

Let, $x_{t}^{k}$ and $y_{t}^{k}$ be the binary variables, where in each time
period $t\in\mathcal{T}$, a battery compartment $k\in\mathcal{K}$ can operate
in one of the three modes: charging $(x_{t}^{k}=1$ & $y_{t}^{k}=0)$,
discharging $(x_{t}^{k}=0$ & $y_{t}^{k}=1)$, and idling $(x_{t}^{k}=0$ &
$y_{t}^{k}=0)$ as shown in (3a). Let, $\widehat{p}_{t}^{up,k}$ and
$\widehat{p}_{t}^{dn,k}$ be the reserve-up and reserve-down bid quantities
with corresponding bid prices $\lambda_{t}^{up}$ and $\lambda_{t}^{dn}$,
respectively. The upper limits for the quantity bids are the maximum allowable
charging $M^{c,k}$ and discharging $M^{d,k}$ rates, respectively (see (3b) and
(3c)). It is assumed that the up and down bids are accepted, as offered, with
probabilities $(\pi_{t}^{a,up})$ and $(\pi_{t}^{a,dn})$, respectively.

The BSS receives a capacity payment for the offered quantity on standby and
receives an additional payment if called upon to deploy. The probabilities for
up and down deployment are $(\pi_{t}^{d,up})$ and $(\pi_{t}^{d,dn})$,
respectively. The estimates of the actual dispatch of the up $(p_{t}^{up,k})$
and down $(p_{t}^{dn,k})$ capacities are given in (3d) and (3e). Each
compartment can be charged from the actual deployment of the reserve-down
capacity and/or with the purchase of power from the RT market $(rt_{t}^{b,k})$
as given in (3f). The BSS is considered to discharge only to deploy the
reserve-up capacities, as in (3g). The maximum amount of power that can be
stored in a compartment is controlled by (3h). The balance of the stored power
is provided in (3i), which ensures that the power in compartment $k$ in time
$t$, $(\phi_{t}^{b,k})$, is equal to the power added to
$(p_{t}^{dn,k}+rt_{t}^{b,k})$ and substracted from $(p_{t}^{up,k})$ its power
level in the previous time period $(\phi_{t-1}^{b,k})$.

As stated earlier, the BSS participates in the reserve market with the aim to
maximize its profit as provided in the objective function (3j). The first term
in $f_{p_{2}}^{b}$ represents the revenue from providing the reserve
capacities $(r^{cap})$ in the day-ahead market, of which the components are
shown in (4).

$r^{cap}=\sum_{t\in\mathcal{T}}\Bigg{[}\lambda_{t}^{up}\pi_{t}^{a,up}\sum_{k\in\mathcal{K}}\widehat{p}_{t}^{up,k}+\lambda_{t}^{dn}\pi_{t}^{a,dn}\sum_{k\in\mathcal{K}}\widehat{p}_{t}^{dn,k}\Bigg{]},$
(4)

where the first and second terms represent the revenue from accepted reserve
up and down bids, respectively, with corresponding prices $\lambda_{t}^{up}$
and $\lambda_{t}^{dn}$. The second term in $f_{p_{2}}^{b}$ represents the
revenue from the actual deployment of the reserve capacities $(r^{dep})$, the
components of which are shown in (5).

$r^{depl}=\sum_{t\in\mathcal{T}}\Bigg{[}\lambda_{t}^{rt}\pi_{t}^{d,up}\sum_{k\in\mathcal{K}}p_{t}^{up,k}+\lambda_{t}^{rt}\pi_{t}^{d,dn}\sum_{k\in\mathcal{K}}p_{t}^{dn,k}\Bigg{]},$
(5)

where $\lambda_{t}^{rt}$ is the real-time price of electricity at time $t$.
The cost of purchasing RT power to charge the BSS ($c^{\phi}$) and the battery
degradation cost ($c^{deg}$) in $f_{p_{2}}^{b}$ are shown in (6) and (7),
respectively. The method of calculating the degradation cost is adopted from
Ortega-Vazquez (2014).

$c^{\phi}=\sum_{t\in\mathcal{T}}\lambda_{t}^{rt}\sum_{k\in\mathcal{K}}rt_{t}^{b,k}.$
(6)

$c^{deg}=\sum_{t\in\mathcal{T}}\sum_{k\in\mathcal{K}}C^{K}\Big{|}\frac{m^{k}}{100}\Big{|}\Bigg{[}\frac{\sum_{t\in\mathcal{T}}\big{(}p_{t}^{up,k}+p_{t}^{dn,k}\big{)}}{BC^{k}}\Bigg{]},$
(7)

where $m^{k}$ is the slope of the linear approximation of the battery life as
a function of the cycles, $BC^{k}$ is the capacity, and $C^{k}$ is the cost of
compartment $k\in\mathcal{K}$.

We represent P2 in a compact form as (8) with $\bm{y}$ being the its decision
variables and $\bm{Y_{f}}$ being its decision space.

$\displaystyle\begin{split}\max_{\bm{y}}f_{p_{2}}^{b}(\bm{y})\\\
\bm{y}\in\bm{Y_{f}}\end{split}$ (8)

###  4.3 P3: Model for joint operation of fast-charging EV hub and BSS

We formulate the joint operation of the charging hub and BSS as a bi-objective
cooperative optimization model. The joint operation results from the hub’s use
of a portion of the BSS capacity. The hub operator minimizes the cost of
meeting the charging demand by optimally choosing the hourly DA commitment
levels and arbitrage strategy for the leased BSS capacity. Whereas, the BSS
operator maximizes profit by optimally biding in the DA reserve market and
leasing a portion of BSS capacity to the hub operator.

$x_{t}^{k}+y_{t}^{k}\leq 1,\quad\forall k\in\mathcal{K},\forall
t\in\mathcal{T}.$ (9a) $\widehat{p}_{t}^{up,k}<=M^{d,k}y_{t}^{k},\quad\forall
t\in\mathcal{T}.$ (9b) $\widehat{p}_{t}^{dn,k}<=M^{c,k}x_{t}^{k},\quad\forall
t\in\mathcal{T}.$ (9c)
$\sum_{k\in\mathcal{K}}p_{t}^{up,k}=\pi_{t}^{d,up}\sum_{k\in\mathcal{K}}\widehat{p}_{t}^{up,k},\quad\forall
t\in\mathcal{T}.$ (9d)
$\sum_{k\in\mathcal{K}}p_{t}^{dn,k}=\pi_{t}^{d,dn}\sum_{k\in\mathcal{K}}\widehat{p}_{t}^{dn,k},\quad\forall
t\in\mathcal{T}.$ (9e) $\phi_{t}^{a,k}+\phi_{t}^{b,k}\leq M^{k},\quad\forall
k\in\mathcal{K},\forall t\in\mathcal{T}.$ (9f)
$\phi_{t}^{a,k}=\phi_{t-1}^{a,k}+da_{t}^{a,k}+rt_{t}^{a,k}-a_{t}^{ev,k}-a_{t}^{rt,k},\quad\forall
k\in\mathcal{K},\forall t\in\mathcal{T}.$ (9g)
$\phi_{t}^{b,k}=\phi_{t-1}^{b,k}+p_{t}^{dn,k}+rt_{t}^{b,k}-p_{t}^{up,k},\quad\forall
k\in\mathcal{K},\forall t\in\mathcal{T}.$ (9h)
$da_{t}^{a,k}+rt_{t}^{a,k}+p_{t}^{dn,k}+rt_{t}^{b,k}\leq
M^{c,k}x_{t}^{k},\quad\forall k\in\mathcal{K},\forall t\in\mathcal{T}.$ (9i)
$a_{t}^{ev,k}+a_{t}^{rt,k}+p_{t}^{up,k}\leq M^{d,k}y_{t}^{k},\quad\forall
k\in\mathcal{K},\forall t\in\mathcal{T}.$ (9j)
$da_{t}^{com}=\sum_{k\in\mathcal{K}}da_{t}^{a,k}+da_{t}^{ev}+da_{t}^{rt},\quad\forall
t\in\mathcal{T}.$ (9k)
$ev_{t}^{l}=da_{t}^{ev}+\sum_{k\in\mathcal{K}}a_{t}^{ev,k}+rt_{t}^{ev},\forall
t\in\mathcal{T}.$ (9l) $\begin{split}\min
f_{p_{3}}^{a}=&\min\sum_{t\in\mathcal{T}}\Bigg{[}\Bigg{\\{}\lambda_{t}^{da}da_{t}^{ev}+(\lambda_{t}^{da}-\lambda_{t}^{rt})da_{t}^{rt}+\lambda_{t}^{rt}rt_{t}^{ev}\Bigg{\\}}+\\\
&\lambda_{t}^{da}\sum_{k\in\mathcal{K}}da_{t}^{a,k}+\lambda_{t}^{rt}\sum_{k\in\mathcal{K}}rt_{t}^{a,k}-\lambda_{t}^{rt}\sum_{k\in\mathcal{K}}a_{t}^{rt,k}+c^{deg}\bigg{(}(1+\alpha)\sum_{k\in\mathcal{K}}\Big{(}a_{t}^{ev,k}+a_{t}^{rt,k}\Big{)}\bigg{)}\Bigg{]}.\end{split}$
(9m) $\begin{split}\max
f_{p_{3}}^{b}=&\max\sum_{t\in\mathcal{T}}\sum_{k\in\mathcal{K}}\Bigg{[}\Bigg{\\{}\pi_{t}^{a,up}\lambda_{t}^{u}\widehat{p}_{t}^{up,k}+\pi_{t}^{a,dn}\lambda_{t}^{d}\widehat{p}_{t}^{dn,k}+\\\
&\pi_{t}^{d,up}\lambda_{t}^{rt}p_{t}^{up,k}+\pi_{t}^{d,dn}\lambda_{t}^{rt}p_{t}^{dn,k}-\lambda_{t}^{rt}rt_{t}^{b,k}-c^{deg}p_{t}^{up,k}\Bigg{\\}}+c^{deg}\alpha\Big{(}a_{t}^{ev,k}+a_{t}^{rt,k}\Big{)}\Bigg{]}.\end{split}$
(9n)

The constraints in (9a), (9b), (9c), (9d), and (9e) are same as in P2. The
total power stored in compartment $k$ by the hub $(\phi_{t}^{a,k})$ and the
BSS $(\phi_{t}^{b,k})$ must be less than or equal to the total storage
capacity of the compartment, see (9f). The balance of power stored by the hub
in compartment $k$ is maintained using (9g). This ensures that the power
stored by the hub ($\phi_{t}^{a,k}$) at time $t$ is equal to the power in the
previous time ($\phi_{t-1}^{a,k}$) plus the power added from DA commitment
($da_{t}^{a,k}$) and/or the RT market ($rt_{t}^{a,k}$) minus the power
discharged to meet EV charging demand ($a_{t}^{ev,k}$) and/or the power sold
back to the grid ($a_{t}^{rt,k}$). Similarly, the balance of power stored by
the BSS in a compartment is maintained by (9h). The upper limits of the power
added to and discharged from a compartment by the hub and BSS are maintained
by (9i) and (9j), respectively. The balance of hourly usage of DA commitment
and fulfillment of EV charging demand are ensured in (9k) and (9l),
respectively. The objective functions for the hub and the BSS operators in
joint operation mode are provided by (9m) and (9n), respectively. The
objective function for the hub in (9m), in addition to the terms of the
objective function of P1 shown within braces, has four other terms: cost of
storing DA power, cost of storing RT power, cost of selling BSS power to RT
grid, and the cost of leasing battery compartments. The objective function for
BSS in (9n) contains all the terms, within braces, from the objective function
of P2. The only additional term accounts for the profit from leasing BSS
capacity to the hub. Equation (10) below represents the compact form of P3,
where $\bm{z}$ denotes the vector of decision variables and $\bm{Z_{f}}$
denotes the decision space.

$\displaystyle\begin{split}\min_{\bm{z}}f_{p_{3}}^{a},\max_{\bm{z}}f_{p_{3}}^{b}\\\
\bm{z}\in\bm{Z_{f}}.\end{split}$ (10)

##  5 Solution approach

The formulation P1 is a linear program and P2 is a mixed integer linear (MILP)
program, both of which can be solved using any commercial solver. Let $d_{1}$
be the solution of problem P1, i.e.,
$d_{1}=\\{\min_{\bm{x}}f_{p_{1}}^{a}:\bm{x}\in\bm{X_{f}}\\}$, and $d_{2}$ be
the solution of problem P2, i.e.,
$d_{2}=\\{\max_{\bm{y}}f_{p_{2}}^{b}:\bm{y}\in\bm{Y_{f}}\\}$. The problem P3
is a bi-objective mixed integer linear program and can be solved using
algorithms such as $\epsilon$-constraint method Mavrotas (2009), the outcome
of which would be a Pareto frontier where each point on the frontier is a
vector comprising objective function values for both the players i.e., the hub
and the BSS. The Pareto frontier with multiple solutions does not provide a
rule to select a solution that is convenient to both the players. A Nash
bargaining solution (NBS) of P3 can identify a single solution from the
frontier that is fair to both participants Nash Jr (1950) Charkhgard et al.
(2022). This solution satisfies a set of axioms to obtain a fair bargain for
all the players. The axioms, as proposed by Nash Nash Jr (1950) include
individual rationality, Pareto optimality, symmetry, linear invariance, and
independence of irrelevant alternatives. This is the advantage of NBS compared
to any other method (e.g., Khanal and Charkhgard ) that provides a single
Pareto optimal solution. The Nash bargaining problem for P3 can be formulated
as follows, where $d_{1}$ and $d_{2}$ are the disagreement points that
represent the payoff of each player under no cooperation.

$\displaystyle\begin{split}\max\quad&(d_{1}-f_{p_{3}}^{a})(f_{p_{3}}^{b}-d_{2})\\\
&\bm{z}\in\bm{Z_{f}},\\\ &f_{p_{3}}^{a}\leq d_{1},\\\ &f_{p_{3}}^{b}\geq
d_{2}.\end{split}$ (11)

The terms $(d_{1}-f_{p_{3}}^{a})$ and $(f_{p_{3}}^{b}-d_{2})$, in the
objective function above, capture the added benefits from cooperation for the
cost-minimizing hub operator and the profit-maximizing BSS operator with
respect to their disagreement points. The above formulation (11) is a maximum
multiplicative program and is not solvable as it is. Hence we present a
second-order cone program (SOCP) reformulation as follows.

$\begin{split}\max\gamma\\\ 0\leq\gamma\leq\sqrt{\tau_{1}^{0}\tau_{2}^{0}}\\\
0\leq\tau_{1}^{0}=d_{1}-f_{p_{3}}^{a}\\\
0\leq\tau_{2}^{0}=f_{p_{3}}^{b}-d_{2}\\\ \bm{z}\in\bm{Z_{f}}\end{split}$ (12)

The above formulation is a SOCP since any constraint of the form $\\{u,v,w\geq
0:u\leq\sqrt{vw}\\}$ is equivalent to $\Big{\\{}u,v,w\geq
0:\sqrt{u^{2}+(\frac{(v-w)^{2}}{2}}\leq\frac{v+w}{2}\Big{\\}}$. The
formulation in (12) can be solved to optimality using most commercial solvers.
Interested readers are referred to Ben-Tal and Nemirovski (2001) for more
details on SOCP reformulation.

##  6 Numerical Case Study

Figure 2: Hub and BSS location near a major street intersection in Tampa,
Florida, USA, as considered in the numerical case study

We build a representative case study to demonstrate the effectiveness of our
proposed model. The hub and the BSS are assumed to be connected to the same
bus of the power network and hence are subjected to the same electricity
prices. The BSS specification considered is similar to one of the battery
storage systems operated by an electric power provider in the Tampa Bay area,
Florida, USA. It is comprised of six compartments, each with a maximum storage
capacity of 4000 kW. Each compartment maintains a minimum power of 500 kW and
has a maximum charge/discharge rate of 3000 kW/hr. The charging hub is
considered to be located near a busy roadway intersection in the Tampa Bay
area (see figure 2) and is assumed to house 150 fast charging stations capable
of charging the EVs at the maximum rate of 100kW/hr. It is assumed that, to
limit arbitrage by the hub, maximum DA commitment is capped at two times the
hourly expected EV charging demand.

The number of EVs seeking to charge in the hub at any hour is determined as
follows. The hourly average traffic flow data through the intersection
($N_{t}$), as depicted in Figure 2, is obtained from the Florida Department of
Transportation FDOT (2023). We adopted traffic data for the month of September
2023. It is assumed that the $\alpha\%$ of the traffic are EVs and $\beta\%$
of those EVs use public charging facilities. Then the average number of the
EVs passing through the intersection in each hour $t$ that might seek to
charge at the hub is $\alpha\beta N_{t}$. These average values are used as the
rate parameters for Poisson distributions to generate the hourly numbers of
EVs that are potential candidates to charge in the hub $(\hat{N_{t}})$. The
actual hourly number of EVs that charge in the hub $(n_{t})$ is found by using
a binomial distribution with parameters $(\hat{N_{t}},p_{t})$, where $p_{t}$
is the probability of an EV with potential to charge actually receiving
charge. The values of $p_{t}$ are adopted from the EV charging behavior study
by Idaho National Lab Deng et al. (2018). We consider $\alpha$ = 0.25 Conzade
et al. (2021) and $\beta$ = 0.42 Narayan et al. (2022). EVs arriving at the
hub are assumed to have three different battery sizes, 50 kW, 75 kW, and 100
kW with probabilities of 0.3, 0.4, and 0.3, respectively. As per the study
conducted by Idaho National Laboratory in the U.S. Idaho National Laboratory
(2015), the amount of charge that an EV seeks to receive varies between 5% to
95% of its battery size. The hourly values of the total demand for charge in
the hub for all 30 days are shown in Figure 3(a). The blue lines in the figure
represent the $10^{th}$, $50^{th}$, and $90^{th}$ percentiles. Since the Tampa
Bay area electricity supply is regulated and does not have DA and RT markets,
the hourly DA and RT prices are adopted from the electricity market (ERCOT) in
the State of Texas, U.S. ERCOT . The DA and RT price profiles are given in
Figures 3(b) and 3(c), respectively, with the $10^{th}$, $50^{th}$, and
$90^{th}$ percentiles marked in blue.

(a)

(b)

(c)

Figure 3: (a) Total EV charging demands, (b) Day-ahead prices, and (c) Real-
time prices

The battery storage system operator estimates the probability of acceptance of
the reserve up and down bids, as well as the probability of deployment of the
accepted bids via simulation of the market clearing mechanism. The simulation
uses published historical data from the ERCOT market for the price/quantity
bids for all market participants, total accepted quantities in the DA reserve
market, and the total deployed up and down quantities. The simulation yields
the market clearing prices and the accepted bids. The number of accepted bids
divided by the total number of bids yields the probability of acceptance of a
reserve bid. Similarly, the total quantity of reserve power deployed divided
by the total quantity of accepted bids yields the probability of deployment.
The simulation results are summarized in Figure 4. Figures 4(a) and 4(b)
depict the probability of acceptance of reserve up and reserve down bids,
respectively. The price of accepted up and down bids are shown in figure 4(c)
and 4(d), respectively. The deployment probability for the accepted reserve up
and down bids are shown in figures 4(e) and 4(f), respectively.

(a)  Up bid acceptance probability

(b) Down bid acceptance probability

(c) Accepted price of up bids

(d) Accepted price of down bids

(e) Deployment probability of up bids

(f) Deployment probability of down bids

Figure 4: Results of reserve market clearing simulation using ERCOT data

The mathematical models are solved using Gurobi 9.5.2 in Python 3.8 using a
computer with an Intel<EMAIL_ADDRESS>processor supported by 32 GB RAM. All
the models are solved to an optimality gap of 0.05% or less. In what follows,
we describe the scenarios for which the optimization models (P1, P2, and P3,
described earlier) are solved and then present the corresponding results.

###  6.1 Results

We first solve our models P1, P2, and P3 with all the input parameters set at
their median values as shown in Figures 3 and 4. We solve P3 using a Nash
bargaining approach (as presented in (12)). We also solve P3 using a total
cost minimization (TCM) approach (13) to demonstrate the benefit of the NBS in
terms of objective function value fairness for both the hub and the BSS in the
joint operation mode. In the TCM approach, we minimize the combined cost of
the hub and the BSS (negative of its profit) given as

$\displaystyle\begin{split}\min_{\bm{z}}f_{p3}^{a}-f_{p3}^{b}\\\
\bm{z}\in\bm{Z_{f}},\end{split}$ (13)

which is a MILP and can be solved using any commercial solver. The results are
shown in Table 1. When the hub and the BSS act independently (as in P1 and P2,
respectively), both the hub cost and BSS profit are inferior to those obtained
from the joint operation (TCM and NBS). In the TCM solution, the hub being a
larger entity influences the solution and receives the most benefit from
cooperation with 11.02% cost reduction while the BSS increases its profit only
by 0.43% compared to the solutions from independent operations. Whereas, the
NBS offers a fair distribution of the benefits of cooperation with a 7.32%
cost reduction for the hub and a 5.23% profit increase for the BSS.

| Independent Operation | Joint Operation per TCM (P3) | Joint Operation per NBS (P3)
---|---|---|---
Hub cost ($) | 30,262.87 (from P1) | 26,929.22 (-11.02%) | 28,048.70 (-7.32%)
BSS profit ($) | 7,073.65 (from P2) | 7,103.95 (0.43%) | 7,443.91 (5.23 %)
Table 1: Hub costs and BSS profits

Since DA commitment is a primary source to procure power by the hub to control
cost in meeting the charging demand, we examined the DA commitment behavior
and DA power usage in independent as well as joint operation modes. Figure 5
shows the hourly DA commitments and DA usage when all the model parameters
(prices, charging demand, and reserve market probabilities) are at their
median values. The outer boundary line of each graph is obtained by joining
the hourly DA commitment values. The shaded areas under the line represent the
amounts of DA power used for charging EVs, storing in the BSS, and selling
back to the RT market. We note that since the median RT prices are generally
higher than DA prices, the hub’s choices for DA commitment in all operation
modes are similar, except for hour 18. It can be observed that the hub
generally commits to more than the EV charging requirement and sells the
excess mostly to the RT market and only sparingly stores DA power in the BSS.
Figure 6 also shows that when all parameters are at their median values, most
of the EV charging power is sourced from the DA market except for a few hours
when charging is supplemented by power from RT and BSS. Figure 7 depicts the
DA commitment pattern for various RT prices when DA price and other parameters
are their median values. As expected, when the RT price is lower than DA (RT
at the 10th percentile and DA at the median), the hub chooses to minimize its
DA commitment. In both of the other cases, RT prices are higher than DA and
hence the DA commitments are similar.

Figure 5: DA commitment and utilization at the median values of model
parameters Figure 6: EV charging demand and sources for meeting the demand at
median values of model parameters Figure 7: DA commitment with different RT
prices with DA price at the median

Since the BSS aims to maximize its profit by cooperating with the hub while
also participating in the reserve market, we examined BSS’ reserve market
bidding behavior in independent and joint operation modes (see Figure 8). It
can be seen that when BSS cooperates with the hub (P3), it offers less bid to
the reserve market compared to its independent operation (P2), which allows
for the use of BSS storage capacity by the hub. It can also be noted that in
the TCM approach, the hub, being the larger entity, skews the solution in its
favor by using a larger share of the storage capacity, thus reducing the bid
offers by the BSS. In the more equitable approach (NBS), the BSS is allowed to
use more of its capacity, when appropriate, to increase its profit, and hence
we see more bid offers compared to TCM. We also note that BSS makes more down
bids in the early hours of the day to store power, and submits more up bids in
later hours to discharge when prices are higher. Figure 9 depicts the BSS
capacity usage under independent and joint operation modes. Similar to that
noted from Figure 8, the NBS solution in Figure 9 allows for more fair usage
of storage capacity by the BSS as opposed to being dominated by the hub, as in
the TCM solution.

Figure 8: Reserve bid patterns by BSS Figure 9: BSS power

###  6.2 Sensitivity Analysis

The hub and the BSS while operating jointly guided by the NBS must consider a
multitude of parameters (including electricity and reserve market prices, EV
charging demand, and down/up acceptance and deployment probabilities) in
making respective cost-minimizing and profit-maximizing decisions. Even if we
consider only three values for each of the parameters (low (10th percentile),
median, and high (90th percentile)), there are $3^{9}$ possible combinations
to examine. Hence, we group the parameters into two subsets, each having a
more direct influence on a decision-maker.

EV charging demand
---
|  | Low | Median | High
|  | RT price | RT price | RT price
|  | Low | Median | High | Low | Median | High | Low | Median | High
| Low | 4.17 | 42.54 | 170.87 | 2.50 | 25.32 | 110.45 | 1.75 | 17.93 | 76.90
| Median | 5.16 | 11.92 | 52.96 | 4.35 | 7.32 | 31.97 | 3.76 | 5.14 | 22.33
DA price | High | 12.02 | 11.33 | 23.04 | 11.67 | 7.58 | 14.01 | 10.63 | 5.76 | 9.78
Table 2: Sensitivity of DA and RT prices and EV charging demand on hub’s cost
reduction (in percentage) from NBS compared to its individual solution (P1)

We examine the sensitivity of DA and RT prices and EV charging demand on the
hub’s cost by considering three possible values (low, median, and high) of the
parameters, which resulted in $3^{3}$ different combinations. For each of
these combinations, all other model parameters that are concerned with the
reserve market and hence primarily affecting the BSS’ profit are kept at their
median values. The sensitivity is measured by the percentage reduction in the
hub’s NBS cost compared to its cost when operating independently (P1). This is
calculated as {(NBS cost - P1 cost)/(P1 cost)}$\times$ 100 and are presented
in Table 2. Since the NBS differs from the P1 solution due to its ability to
arbitrage using BSS’ storage capacity, all the numbers in the table reflect
the incremental cost reduction benefits by the increase in revenue from
arbitrage. Clearly, for some DA/RT prices and EV charging demand combinations,
the arbitrage potentials for the hub are either very high or very low. The
results show that the very high potentials are generally for the combinations
when the DA prices are low and RT prices are high. The arbitrage potential is
low when both the DA and RT prices are low. It may also be noted that the
hub’s arbitrage potential decreases, irrespective of DA and RT price
combinations, with an increase in EV charging demand.

Figure 10: BSS usage by the hub under different DA prices, RT prices, and EV
charging demand

Among the combinations studied in Table 2, the highest cost reduction
potential (170.87 %) occurs when the DA price is low, the RT price is high,
and the EV charging demand is low. To further examine this case, we looked at
the BSS usage by the hub across all hours of the day. This is depicted in the
leftmost part of Figure 10, where the line graphs show RT prices to be
significantly higher than the DA prices, and the bar graphs show multiple
charge/discharge activities by the hub. We note that the significant price
differential between DA and RT prices encourages the hub operator to increase
its level of arbitrage via high usage of BSS storage capacity and thus reduce
its cost. We also examine two other cases: low DA and RT prices and high EV
charging demand (cost reduction of 1.75%); and high DA prices, low RT prices,
and low EV charging demand (cost reduction of 12.02%). The BSS usage by the
hub for these two cases is depicted in the middle and the rightmost segment of
Figure 10, respectively. For the case with a price reduction of 1.75%, the low
differential between DA and RT prices prompts a limited hub engagement with
the BSS (i.e. lower level of charge/discharge activity) yielding a lower
potential for cost reduction. Whereas in the case of a 12.02 % cost reduction,
we note that the RT prices are significantly lower than the DA prices and also
do not vary widely across the hours. This yields a moderate amount of
arbitrage by the hub using the BSS capacity.

Parameters | F-statistic | F-critical | P-value | Decision
---|---|---|---|---
Up bid market clearing price $(\lambda^{up})$ | 11.92 | 4.279 | 0.002 | Significant
Down bid market clearing price $(\lambda^{dn})$ | 1.82 | 4.279 | 0.18 | Not Significant
Probability of acceptance of up bid $(\pi^{a,up})$ | 8.82 | 4.279 | 0.006 | Significant
Probability of acceptance of down bid $(\pi^{a,dn})$ | 0.96 | 4.279 | 0.33 | Not Significant
Probability of deployment of up bid $(\pi^{d,up})$ | 83.79 | 4.279 | 3.94 $\times$ $10^{-9}$ | Significant
Probability of deployment of down bid $(\pi^{d,dn})$ | 32.10 | 4.279 | 9.07 $\times$ $10^{-6}$ | Significant
| Up bid market clearing price $(\lambda^{up})$ $\times$
---
probability of deployment of up bid $(\pi^{d,up})$
6.210 | 4.279 | 0.020 | Significant
| Probability of deployment of up bid $(\pi^{d,up})$ $\times$
---
probability of deployment of down bid $(\pi^{d,dn})$
6.36 | 4.279 | 0.019 | Significant
Table 3: Analysis of variance (ANOVA) to examine the sensitivity of reserve
market parameters on BSS’s profit increase (in percentage) from NBS compared
to its individual solution (P2)

Hereafter, we study the sensitivity of six reserve market parameters on the
BSS’s ability to increase profit by allowing the hub to use its storage
capacity. This profit increase potential (in percentage) is calculated as
{(NBS profit - P2 profit)/ P2 profit} $\times$ 100 %. The parameters studied
are up and down bid market clearing prices, probabilities of acceptance of up
and down bids, and probabilities of deployment of up and down bids. We used a
$2^{6-1}$ fractional factorial design-based analysis of variance (ANOVA)
study, see chapter 8 of Montgomery (2017) for the details of $2^{k}$ factorial
design. The 10th and 90th percentile values are considered as the two levels
of the parameters and the highest interaction term is used as the fractional
factorial design generator. In the thirty-two experiments for the $2^{6-1}$
fractional factorial design, the percentage profit varied between 1.52% and
45%. The ANOVA results are presented in Table 3. Of the six reserve market
parameters, the down-bid market clearing price and its probability of
acceptance are found not significant with a type-I error level of 5%. We also
found two of the two-parameter interactions to be significant, although with a
small margin. Of the four statistically significant parameters, the
probabilities of deployment of up and down bids are the most significant ones
with very high F-statistic values of 83.79 and 32.10, respectively. This is
expected as the actual up-and-down deployments are the major contributor to
BSS revenue from the reserve market. Examination of the BSS behavior from NBS
reveals that for high values of up and down deployment probabilities, the BSS
uses a higher proportion of its storage capacity for the reserve market,
leaving a lower capacity for the hub.

##  7 Concluding remarks

In this paper, we develop a bi-objective joint operation model between a cost-
minimizing fast-charging EV hub and a profit-maximizing grid-connected stand-
alone battery storage system. The BSS primarily uses its capacity to
participate in the reserve market. It also cooperates with the hub by allowing
the hub to use its storage capacity for arbitrage. The bi-objective problem is
reframed as a Nash bargaining problem, which is then reformulated into a
second-order cone program. We develop a sample case study problem adopting DA
and RT electricity prices from the ERCOT market and generating simulated
values for the reserve market prices and their probabilities as well as the EV
charging demand. We demonstrate the benefits/fairness of our proposed Nash
bargaining solution with the solution from a total cost minimization approach
as well as the independent solutions without cooperation. A detailed study of
parameter sensitivity reveals operational guidance for both the hub and the
BSS under various parameter value combinations.

Our work demonstrates that the EV charging hubs can significantly benefit
financially by utilizing the storage capacity of grid-connected stand-alone
BSS. We also show that the joint operation under NBS yields added benefits for
the BSS. The results from the numerical case study, for a hub with 150 DC
fast-charging stations and a BSS with 24 MW storage capacity, show that under
specific operating conditions, the hub and the BSS can financially gain up to
170% and 45%, respectively. These findings create new opportunities for the
growth of stand-alone BSS in support of large-scale fast-charging hubs,
ultimately increasing the rate of EV adoption.

Our study has the following limitations. Our model does not explicitly
consider uncertainties in the energy market and the reserve market parameters.
We investigated the performance of our model for various parameter value
combination scenarios. A distributionally robust approach to the Nash
bargaining problem could be a way to further accommodate these uncertainties
Liu et al. (2018) Peng et al. (2021). In our problem, we have considered
participation in the reserve up and down markets as the primary activities of
the battery storage system. Whereas in most electricity markets, BSS also
participates in other reserve market activities such as responsive reserve,
non-spinning responsive reserve, and contingency reserve. Consideration of
these additional activities will add more variables and constraints without
changing the fundamental nature of the model as we have developed. For
simplicity, we have considered only the hourly variation in the RT market
prices in formulating our model, which can be extended to accommodate more
frequent changes in the RT market.

## References

  * Conzade et al. [2021] Julian Conzade, Andreas Cornet, Patrick Hertzke, Russell Hensley, Ruth Heuss, Timo Möller, Patrick Schaufuss, Stephanie Schenk, Andreas Tschiesner, and Karsten von Laufenberg.  Why the automotive future is electric.  https://www.mckinsey.com/industries/automotive-and-assembly/our-insights/why-the-automotive-future-is-electric, 2021.  (accessed: 26 July 2022).
  * [2] Tesla.  Supercharger.  https://www.tesla.com/en_eu/supercharger.  (accessed : 2 July 2022).
  * NREL [2023] NREL.  Utility-scale battery storage.  https://atb.nrel.gov/electricity/2023/utility-scale_battery_storage, 2023.  (accessed: 01 December 2023).
  * Zheng et al. [2020] Yanchong Zheng, Hang Yu, Ziyun Shao, and Linni Jian.  Day-ahead bidding strategy for electric vehicle aggregator enabling multiple agent modes in uncertain electricity markets.  _Applied Energy_ , 280:115977, 2020.
  * Paudel et al. [2023] Diwas Paudel, Nicolas Bustos, and Tapas K Das.  A distributionally robust approach for day-ahead power procurement by ev charging hubs.  In _2023 IEEE Power & Energy Society General Meeting (PESGM)_, pages 1–5. IEEE, 2023.
  * Paudel and Das [2023] Diwas Paudel and Tapas K Das.  A deep reinforcement learning approach for power management of battery-assisted fast-charging ev hubs participating in day-ahead and real-time electricity markets.  _Energy_ , 283:129097, 2023.
  * Subramanian and Das [2019] Vignesh Subramanian and Tapas K Das.  A two-layer model for dynamic pricing of electricity and optimal charging of electric vehicles under price spikes.  _Energy_ , 167:1266–1277, 2019.
  * DeForest et al. [2018] Nicholas DeForest, Jason S MacDonald, and Douglas R Black.  Day ahead optimization of an electric vehicle fleet providing ancillary services in the los angeles air force base vehicle-to-grid demonstration.  _Applied energy_ , 210:987–1001, 2018.
  * Sarker et al. [2015] Mushfiqur R Sarker, Yury Dvorkin, and Miguel A Ortega-Vazquez.  Optimal participation of an electric vehicle aggregator in day-ahead energy and reserve markets.  _IEEE transactions on power systems_ , 31(5):3506–3515, 2015.
  * Kazemi et al. [2017] Mostafa Kazemi, Hamidreza Zareipour, Nima Amjady, William D Rosehart, and Mehdi Ehsan.  Operation scheduling of battery storage systems in joint energy and ancillary services markets.  _IEEE Transactions on Sustainable Energy_ , 8(4):1726–1735, 2017.
  * Padmanabhan et al. [2019] Nitin Padmanabhan, Mohamed Ahmed, and Kankar Bhattacharya.  Battery energy storage systems in energy and reserve markets.  _IEEE Transactions on Power Systems_ , 35(1):215–226, 2019.
  * Xu et al. [2017] Bolun Xu, Yishen Wang, Yury Dvorkin, Ricardo Fernández-Blanco, Cesar A Silva-Monroy, Jean-Paul Watson, and Daniel S Kirschen.  Scalable planning for energy storage in energy and reserve markets.  _IEEE Transactions on Power systems_ , 32(6):4515–4527, 2017.
  * Merten et al. [2020] Michael Merten, Christopher Olk, Ilka Schoeneberger, and Dirk Uwe Sauer.  Bidding strategy for battery storage systems in the secondary control reserve market.  _Applied Energy_ , 268:114951, 2020.
  * Nitsch et al. [2021] Felix Nitsch, Marc Deissenroth-Uhrig, Christoph Schimeczek, and Valentin Bertsch.  Economic evaluation of battery storage systems bidding on day-ahead and automatic frequency restoration reserves markets.  _Applied Energy_ , 298:117267, 2021.
  * Hu et al. [2022] Yu Hu, Miguel Armada, and María Jesús Sánchez.  Potential utilization of battery energy storage systems (bess) in the major european electricity markets.  _Applied Energy_ , 322:119512, 2022.
  * Ortega-Vazquez [2014] Miguel A Ortega-Vazquez.  Optimal scheduling of electric vehicle charging and vehicle-to-grid services at household level including battery degradation and price uncertainty.  _IET Generation, Transmission & Distribution_, 8(6):1007–1016, 2014.
  * Mavrotas [2009] George Mavrotas.  Effective implementation of the $\varepsilon$-constraint method in multi-objective mathematical programming problems.  _Applied mathematics and computation_ , 213(2):455–465, 2009.
  * Nash Jr [1950] John F Nash Jr.  The bargaining problem.  _Econometrica: Journal of the econometric society_ , pages 155–162, 1950.
  * Charkhgard et al. [2022] Hadi Charkhgard, Kimia Keshanian, Rasul Esmaeilbeigi, and Parisa Charkhgard.  The magic of nash social welfare in optimization: Do not sum, just multiply!  _The ANZIAM Journal_ , 64(2):119–134, 2022.
  * [20] Ashim Khanal and Hadi Charkhgard.  A criterion space search feasibility pump heuristic for solving maximum multiplicative programs.  _Available at SSRN 4658986_.
  * Ben-Tal and Nemirovski [2001] Aharon Ben-Tal and Arkadi Nemirovski.  On polyhedral approximations of the second-order cone.  _Mathematics of Operations Research_ , 26(2):193–205, 2001.
  * FDOT [2023] FDOT.  Traffic information.  https://www.fdot.gov/statistics/trafficdata, 2023.  (accessed : 03 March 2023).
  * Deng et al. [2018] Qian Deng, Sujit Tripathy, Daniel Tylavsky, Travis Stowers, and Jeff Loehr.  Demand modeling of a dc fast charging station.  In _2018 North American Power Symposium (NAPS)_ , pages 1–6. IEEE, 2018.
  * Narayan et al. [2022] Ajay Narayan, Aakash Krishna, Prasant Misra, Arunchandar Vasan, and Venkatesh Sarangan.  A dynamic pricing system for electric vehicle charging management using reinforcement learning.  _IEEE Intelligent Transportation Systems Magazine_ , 14(6):122–134, 2022.
  * Idaho National Laboratory [2015] Idaho National Laboratory.  Plug-in electric vehicle and infrastructure analysis.  https://inldigitallibrary.inl.gov/sites/sti/sti/6799570.pdf, 2015.
  * [26] ERCOT.  Market information.  https://www.ercot.com/mktinfo.  (accessed : 28 September 2023).
  * Montgomery [2017] Douglas C Montgomery.  _Design and analysis of experiments_.  John wiley & sons, 2017.
  * Liu et al. [2018] Yongchao Liu, Huifu Xu, Shu-Jung Sunny Yang, and Jin Zhang.  Distributionally robust equilibrium for continuous games: Nash and stackelberg models.  _European Journal of Operational Research_ , 265(2):631–643, 2018.
  * Peng et al. [2021] Shen Peng, Abdel Lisser, Vikas Vikram Singh, Nalin Gupta, and Eshan Balachandar.  Games with distributionally robust joint chance constraints.  _Optimization Letters_ , pages 1–23, 2021.


# LiNiO2 as a high-entropy charge- and bond-disproportionated glass

Kateryna Foyevtsova∗ Department of Physics & Astronomy, University of British
Columbia, Vancouver, BC V6T 1Z1, Canada  Stewart Blusson Quantum Matter
Institute, University of British Columbia, Vancouver, BC V6T 1Z4, Canada  Ilya
Elfimov  Department of Physics & Astronomy, University of British Columbia,
Vancouver, BC V6T 1Z1, Canada  Stewart Blusson Quantum Matter Institute,
University of British Columbia, Vancouver, BC V6T 1Z4, Canada  Joerg Rottler
Department of Physics & Astronomy, University of British Columbia, Vancouver,
BC V6T 1Z1, Canada  Stewart Blusson Quantum Matter Institute, University of
British Columbia, Vancouver, BC V6T 1Z4, Canada  George A. Sawatzky
Department of Physics & Astronomy, University of British Columbia, Vancouver,
BC V6T 1Z1, Canada  Stewart Blusson Quantum Matter Institute, University of
British Columbia, Vancouver, BC V6T 1Z4, Canada

###### Abstract

Understanding microscopic properties of LiNiO2, a Li-ion battery cathode
material with extraordinarily high reversible capacity, has remained a
challenge for decades. Based on extensive electronic structure calculations,
which reveal a large number of nearly degenerate phases involving local Jahn-
Teller effect as well as bond and oxygen-based charge disproportionation, we
propose that LiNiO2 exists in a high-entropy charge-glass like state at and
below ambient temperatures. Recognizing the glassy nature of LiNiO2 does not
only explain its key experimental features, but also opens a new path in
designing entropy-stabilized battery cathodes with superb capacities.

###### pacs:

61.43.Fs,64.70.P-,65.40.Gr,71.20.-b,75.47.Lx,71.20.-b,71.30.+h,71.45.Lr

##  I Introduction

LiNiO2 is well known as a promising cathode material for rechargeable Li-ion
batteriesGoodenough (1955, 1958). One of the most fascinating but also
puzzling properties of this oxide is its record-breaking stability upon
repeated charge and discharge cycles. It is indeed quite surprising that the
LiNiO2 system can be cycled more than a few times changing the Li
concentration by up to 80% on each cycleLiu et al. (2015); Arai et al. (1995).
This indicates that there must be something very special about the crystal
structure of LiNiO2 that keeps the original structure in place, and indeed
LiNiO2 has demonstrated a range of structural features that are hard to
understand in a conventional way.

LiNiO2 consists of layers of edge-sharing NiO6 octahedra on a triangular
lattice alternating with Li layers. The formally trivalent Ni ions are in a
low-spin configuration $t_{2g}^{6}e_{g}^{1}$ which makes them Jahn-Teller (JT)
activeGoodenough (1958); Hirakawa et al. (1985); Hirota et al. (1991); Yamaura
et al. (1996). Despite this, however, the NiO6 octahedra do not undergo a
cooperative JT distortion even at the lowest temperatures, as follows from the
numerous Rietveld refinements of the LiNiO2 powder diffraction dataGoodenough
(1958); Arai et al. (1995); Chung et al. (2005). This is in a surprising
contrast with the behavior of its sister compound NaNiO2 that shows a ferro-
orbital ordering of occupied $d_{z^{2}}$ orbitals and an associated
cooperative elongation of the NiO6 octahedra below 460 K as illustrated in
Fig. 1(a)Borgers and Enz (1966); Holzapfel et al. (2004). On the other hand,
structural analyses based on the neutron pair distribution function
(nPDF)Chung et al. (2005), extended x-ray absorption fine structureRougier et
al. (1995), and electron spin resonanceBarra et al. (1999); Reynaud et al.
(2001) measurements, strongly suggest local, but apparently disordered, Ni-O
bond disproportionation in LiNiO2. Disorder effects in this system, which also
include a spin-glass-like behavior at low temperaturesHirota et al. (1991),
are often linked to its unavoidable deviation from stoichiometryArai et al.
(1995); Petit et al. (2006). Yet, the exact nature of the disordered state in
LiNiO2 has so far remained unclear.

Apart from the above, LiNiO2 also possesses a rather peculiar electronic
property, which, although being often dismissed in literature, might be the
key to understanding its structural properties. As was first pointed out by
Kuiper et al.Kuiper et al. (1989), this system is in the negative charge-
transfer regime which brings the NiO6 octahedra closer to the
$(t^{6}_{2g}e^{2}_{g})\underline{L}$ (rather than the formal
$t^{6}_{2g}e^{1}_{g}$) average configuration, where $\underline{L}$ denotes a
ligand hole residing on a molecular orbital formed by the six
oxygen-$2p_{\sigma}$ orbitals in an octahedron. In this regime, the Jahn-
Teller effect might be in competition with another degeneracy lifting
mechanism, observed in rare-earth nickel perovskitesJohnston et al. (2014);
Green et al. (2016); Foyevtsova et al. (2015), which involves formation of two
types of size-disproportionated NiO6 octahedra with configurations
$t^{6}_{2g}e^{2}_{g}$ and $(t^{6}_{2g}e^{2}_{g})\underline{L}^{2}$ with
symmetries and spins corresponding to those of Ni2+ and Ni4+. Not
surprisingly, therefore, a strongly competing bond-disproportionated phase of
LiNiO2, with ordered stripes of collapsed and expanded NiO6 octahedra [Fig.
1(e)], was indeed recently predicted within density functional theory
(DFT)Chen et al. (2011).

Our present theoretical study offers a unifying explanation of the puzzling
structural and electronic properties of LiNiO2 outlined above as well as its
superior performance in Li-ion batteries by recognizing the fundamental role
of entropy in this system. Based on extensive electronic structure
calculations, we demonstrate that LiNiO2 has a large number of nearly
degenerate states involving NiO6 octahedra with $t^{6}_{2g}e^{2}_{g}$,
$(t^{6}_{2g}e^{2}_{g})\underline{L}$ and
$(t^{6}_{2g}e^{2}_{g})\underline{L}^{2}$ configurations, leading to basically
a large number of combinations of these into glassy-like structures which
increases the entropy resulting in a high-entropy material. It is shown that,
as a result of the charge-transfer energy being negative in LiNiO2, this state
involves charge disproportionation on oxygen rather than on nickel atoms.
Although entropy is known to play important roles in polymer science and
protein folding as well as glasses, multi-component oxides, and metal
alloysRost et al. (2015); Yao et al. (2018), these are mostly structurally
driven while our proposal involves an electronic-structure driven mechanism
for the pure material.

Figure 1: Possible orderings of short and long Ni-O bonds in the triangular
lattice of NiO6 octaherda. Nickel (oxygen) atoms are represented by small
(large) circles. The short Ni-O bonds are shown as black connecting lines. In
(a), a three-dimensional view of the lattice is shown on the left-hand side
for a reference. Panels (a) and (b) represent phases where only Jahn-Teller
(JT) distorted NiO6 octahedra are present, whereas panels (c), (d), and (e)
represent phases where, respectively, 50%, 67%, and 100% of the NiO6 octahedra
are size-disproportionated (SD).

##  II Methods

We performed density functional theory (DFT) calculations with the pseudo-
potential code VASPKresse and Furthmüller (1996); Paier et al. (2005) and the
all-electron code WIEN2kBlaha et al. (2001). All of our structural relaxations
were done in VASP, using a high energy cut-off of 550 Ry and fine $k$-grids
with densities close to that of a $12\times 12\times 12$ $k$-grid in a
$5.5\AA\times 5.5\AA\times 5.5\AA$ unit cell. A combination of the local
density approximation (LDA)Perdew and Wang (1992) and the LDA+U method as
designed by Anisimov et al. Anisimov et al. (1991, 1993); Liechtenstein et al.
(1995), with the Ni-$3d$ electrons’ on-site interactions $U=6$ eV and
$J_{\text{H}}=1$ eV, was used to treat exchange and correlation effects.

The neutron pair distribution functions (nPDF), $G(r)$, of the relaxed LiNiO2
phases were simulated with the PDFfit2 programFarrow et al. (2007) using the
following expressionChung et al. (2005):

$G(r)=\frac{1}{rN}\sum_{i,j}\left[\frac{b_{i}b_{j}}{\langle
b\rangle^{2}}\delta(r-r_{ij})\right]-4\pi\rho_{0},$

where $r$ is distance, $b_{i}$ and $r_{ij}$ are the scattering length of the
$i$th atom and the distance between the $i$th and $j$th atom, respectively,
$N$ is the number of atoms, and $\rho_{0}$ is the average atomic density.
Thermal broadening was accounted for through introducing the Debye-Waller
factors. We used theoretical structural parameters, scaled with an empirical
factor of 1.024, and theoretical thermal displacement parameters. In the
PDFfit2 program, we chose 0.08 for the $Q_{\text{damp}}$ parameter and 0.5 for
the linear atomic correlation factor $\delta_{1}$ for all the considered model
LiNiO2 phases except the rhombohedral one where we used $\delta_{1}=0.2$. The
neutron time-of-flight (TOF) powder diffraction patterns were simulated using
the GSAS packageLarson and Dreele (2000); Toby (2001). For both the nPDF and
the TOF powder diffraction simulations, we used the following theoretically
calculated atomic isotropic displacement parameters $U\equiv
U_{11}=U_{22}=U_{33}$, obtained for the temperature of 10 K employing VASP and
the Phonopy packageTogo and Tanaka (2015): $U$(Li) = 0.008 $\AA^{2}$, $U$(Ni)
= 0.0011 $\AA^{2}$, $U$(O) = 0.0025 $\AA^{2}$. The thermodynamic properties
were calculated using the phonon spectra obtained with VASP and Phonopy via
the dynamical matrix method.

##  III Results

Our following discussion will be based on a comparison of a number of
different bond-disproportionated phases of LiNiO2 shown in Fig. 1, which all
turn out to be easily achievable in a DFT structural relaxation. They include
two JT phases [(a) and (b)], the stripy phase with size-disproportionated NiO6
octahedra from Ref. Chen et al., 2011 (e), and two mixed phases where JT-
distorted and size-disproportionated octahedra are simultaneously present [(c)
and (d)]. The high-symmetry rhombohedral $R\bar{3}m$ phase of LiNiO2, which
has been traditionally assigned to LiNiO2 based on powder diffraction
refinement and has no Ni-O bond disproportionation, will also be considered
for completeness. For each phase, we use DFT and the local density
approximation $+U$ (LDA+U) method to perform full lattice relaxation and to
calculate the ground state electronic structure. Ferromagnetic alignment of Ni
spins is adopted throughout all the calculations, as it was established that
the total energy differences between different spin configurations are much
smaller than those between the different structural phases of LiNiO2, as
detailed in the Appendix. This finding of ours agrees with the conclusions of
Mostovoy and KhomskiiMostovoy and Khomskii (2002) who demonstrated decoupling
of spin and orbital degrees of freedom in LiNiO2.

Phase |  $E_{\text{tot}}$ (meV/f. u.) | Gap (eV) | $N_{\text{short}}/N_{\text{long}}$ |  $d_{\text{Ni-O}}$ (Å)
---|---|---|---|---
Rhomb. | 141.3 (147.5) | 0 | - | 1.98
ferro JT | 18.1 (0) | 0.43 | 2 | 1.90 2.12
zigzag JT |  0 (19.1) | 0.56 | 2 | 1.90 2.11
50% SD | 49.6 (60.1) | 0.35 | 1.4 | 1.90 2.06
67% SD | 32.1 (43.5) | 0.51 | 1.25 | 1.90 2.05
100% SD | 26.3 (63.3) | 0.58 | 1 | 1.89 2.06
Table 1: Calculated properties of the LiNiO2 model phases. $E_{\text{tot}}$
are relative total energies, $N_{\text{short}}/N_{\text{long}}$ are the ratios
of the numbers of short and long Ni-O bonds, and $d_{\text{Ni-O}}$ are the
average nearest-neighbor Ni-O bond lengths. A scaling factor of 1.024 has been
applied to the original $d_{\text{Ni-O}}$ values obtained within LDA+U. For
comparison, $E_{\text{tot}}$ values relevant for NaNiO2 are provided in
parentheses.  Figure 2:  Neutron pair distribution functions (nPDF) $G(r)$ of
LiNiO2 and NaNiO2. In (a)-(e), the experimental (black symbols)Chung et al.
(2005) and various theoretically simulated (red lines) nPDFs at the
temperature of 10 K are compared for LiNiO2. Panel (f) shows nPDFs of NaNiO2
simulated using its experimental crystal structure (black line)Holzapfel et
al. (2004) and using a relaxed structure from LDA+U calculations (red line). A
scaling factor of 1.024 has been applied to $r$ in $G(r)$ calculated using
results of LDA+U structural relaxations. Figure 3: Neutron time-of-flight
powder diffractograms of LiNiO2. Experiment (black symbols)Chung et al. (2005)
is compared to a simulation (red lines) where (a) the rhombohedral and (b) the
67% SD model phase of LiNiO2 at 10 K is used. A scaling factor of 1.024 has
been applied to the $d$-spacing in the simulated diffractograms.

We note first that for NaNiO2 (the sister compound) DFT correctly predicts the
JT phase with a ferro-orbital ordering of occupied $d_{z^{2}}$ orbitals [Fig.
1(a)] to be the lowest energy stateMeskine and Satpathy (2005); Chen et al.
(2011). For LiNiO2, on the other hand, we find that the lowest energy state is
the JT phase with a zigzag ordering of occupied $d_{z^{2}}$ orbitals [Fig.
1(b)]Chen et al. (2011). Its energy is by 18.1 meV per formula unit (f. u.)
lower than that of the ferro-orbitally ordered phase and by 141.3 meV/f. u.
lower than that of the rhombohedral phase (see Table 1). The fully and
partially size-disproportionated (SD) phases shown in Figs. 1(c)-(e) have
energies that are not too far from the energy of the ferro-orbital JT phase.
For example, the three-fold rotationally symmetric phase shown in (d), with
two thirds (or 67%) of Ni sites being SD, is higher in energy by only 14
meV/f. u. As one can also see in Table 1, all the considered bond-
disproportionated phases have a small charge gap in LDA+U, which qualitatively
agrees with experimentHirota et al. (1991). Interestingly, the same SD phases
can be obtained for NaNiO2, but, as also shown in Table 1, they are
energetically more strongly removed from the lowest energy state than in the
case of LiNiO2.

Before discussing the important implications following from these total energy
calculations, let us point out that LiNiO2 being in a mixed phase, such as the
67% SD one, could explain amazingly well the findings from both, the pair
distribution function and the powder diffraction analyses, which previously
have been regarded as mutually contradicting. In Figs. 2(a)-(e), the
experimental nPDF of LiNiO2Chung et al. (2005) is compared with simulated
nPDFs of the five model phases from Fig. 1. Both the experiment and the
simulations are performed at the temperature of 10 K. In order to correct for
the underestimation of interatomic distances in LDA+U, an empirical scaling
factor of 1.024 was applied to $r$ in the simulated nPDF $G(r)$. This
procedure gives very good agreement between theory and experiment for NaNiO2,
whose experimental crystal structure is unambiguously known [Fig.
2(f)]Holzapfel et al. (2004). For LiNiO2, it is the 67% SD phase’s nPDF that
has equally good agreement with experiment [Fig. 2(d)], while the zigzag JT
(lowest energy state) phase’s nPDF shows considerable deviations [Fig. 2(b)].
In the bond-disproportionated phases [Figs. 2(b)-(e)], the two nPDF peaks
around 2 Å are particularly sensitive to the presence of SD NiO6 sites and
their evolution from (b) to (e) is easy to understand. Indeed, the relative
heights of these peaks reflect the ratio between the numbers of short and long
nearest-neighbor Ni-O bonds, $N_{\text{short}}/N_{\text{long}}$. As the
fraction of SD sites grows, $N_{\text{short}}/N_{\text{long}}$ gradually
decreases from 2 in the fully JT phases to 1 in the fully SD phase (Table 1)
and so does the intensity of the first nPDF peak relative to that of the
second. Furthermore, only by having SD sites in the system can one accurately
reproduce the measured position of the second peak or, equivalently, the
average length of the long Ni-O bonds (Table 1). We would like to emphasize
that these results have been obtained with keeping the nPDF simulations as ab
initio as possible, which is a step forward compared with the original
analysis of Chung et al. in Ref. Chung et al., 2005. We also note that, in
general, the fraction of SD octahedra may probably vary from one sample to
another depending on their thermal histories and growth conditions, which
should be a subject of further investigations.

As for the powder diffraction experiments, we find that the 67% SD phase gives
a diffraction pattern that is strikingly similar to that of the rhombohedral
phase (Fig. 3). Moreover, our general observation is that all bond-
disproportionated phases with a preserved $C_{3}$ symmetry of the lattice tend
to have very similar diffraction patternspow . In other words, the powder
diffraction pattern of LiNiO2 appears to be barely sensitive to Ni-O bond
disproportionation as long as the $C_{3}$ symmetry is preserved, which in fact
can explain the long-standing confusion around the LiNiO2 crystal structure.

##  IV Discussion

###  IV.1 Role of entropy

Figure 4: (a) The five building block configurations of the NiO6 octahedra in
LiNiO2: the three JT distorted configurations and the two SD configurations.
In color shown are the density iso-surfaces of holes occupying $E_{g}$ oxygen
molecular orbitals. (b) Orbitally projected densities of states (PDOS) in the
zigzag JT phase of LiNiO2. The right panel features projections onto oxygen
molecular orbitals of various symmetries. Both plots are normalized to a
single NiO6 octahedron with only O-$p_{\sigma}$ orbitals participating.

At this point, the question arises as to why LiNiO2 would adopt a metastable
phase like the 67% SD one rather than the lowest energy JT zigzag phase, which
is what the above comparison with experiment strongly suggests, and this is
where the entropy factor comes into play. In fact, given the accuracy
limitations of the LDA+U method, the only robust conclusion that comes from
our total energy calculations is that there exist phases with varying amounts
of JT distorted and SD NiO6 octahedra whose energies lie within only 20 to 30
meV/f. u. above the ground-state energy. While the ferro JT, 67% SD, 50% SD,
and 100% SD phases are just some examples of such phases, further examples can
be easily come up with and tested in a calculation by exploring various ways
of positioning the five building block configurations of a NiO6 octahedron on
a triangular lattice: the three JT distorted and the two SD configurations
[Fig. 4(a)]. With so many metastable phases being so close in energy, a state
where the system is able to fluctuate between them can have a high enough
entropy to be stabilized at ambient temperatures by the entropy term in the
free energy. In this fluctuating state, each NiO6 site can assume any of the
building block configurations, and in an extreme limit would do it in a random
way. In reality, of course, there will be constraints due to local
correlations between different types of sites mostly determined by Pauling’s
principle of electroneutralityPauling (1948). Thus, an expanded NiO6
octahedral site with electronic configuration $d^{8}$ would “attract”
collapsed $d^{8}\underline{L}^{2}$ sites (each having $2e$ less charge due to
the 2 ligand holes than the reference $d^{8}$ site) in its immediate vicinity,
but due to the triangular geometry those will be somehow diluted by
$d^{8}\underline{L}$ and/or $d^{8}$ sites. One can envision this mechanism
driving a propagation of local $d^{8}$, $d^{8}\underline{L}$, and
$d^{8}\underline{L}^{2}$ configurations from a nucleation center ensured by
the entropy to be as random as possible at a given temperature. We believe
that the final entropy-stabilized state will involve local distributions of
$d^{8}$, $d^{8}\underline{L}$, and $d^{8}\underline{L}^{2}$ sites
corresponding to all and any of those depicted in Fig. 1. We can in fact make
the discussion more quantitative by computing and comparing the
configurational entropy of this entropy-stabilized state and that of the
lowest-energy zigzag Jahn-Teller phase. As detailed in the Appendix, the
impact of the vibrational (or phonon) free energy $F^{\text{phon}}$ on the
stability of various LiNiO2 phases is small compared to that of the
configurational entropy and therefore will not be considered here. The
configurational entropy of the zigzag Jahn-Teller phase is due to spin
fluctuations between the $S_{z}=\pm\frac{1}{2}$ states of the Ni ions’ spins:
$S^{\text{conf}}=Nk_{\text{B}}\ln 2=0.693Nk_{\text{B}}$, where $N$ is the
number of formula units in the crystal. For the entropy-stabilized state, let
us assume for simplicity that the $d^{8}$, $d^{8}\underline{L}$, and
$d^{8}\underline{L}^{2}$ sites have equal concentrations of $1/3$, and that
they are randomly distributed. In addition, the JT $d^{8}\underline{L}$ sites
have degeneracy of $g_{1}=6$ corresponding to three possible orientations of
the long Ni-O bond [see Fig. 4(a)] and two $S_{z}$ components of spin
$S=\frac{1}{2}$, while the large octahedron sites $d^{8}$ have degeneracy of
$g_{2}=3$ corresponding to three $S_{z}$ components of spin $S=1$. Neglecting
local correlations between the different types of sites, which otherwise would
somewhat reduce our estimate of $S^{\text{conf}}$, we get that
$S^{\text{conf}}=k_{\text{B}}\ln(3^{N}(g_{1}^{N/3}g_{2}^{N/3}))=2.062Nk_{\text{B}}$.
This upper bound is in fact very large and exceeds the value of $\ln 5=1.61$
in five-component mixtures at equal concentrations, which are conventionally
considered “high-entropy” materialsMurty et al. (2014). We note that, aside
from the chemical composition aspects used in the design of high-entropy
alloysRost et al. (2015); Sarkar et al. (2018), our analysis puts a new twist
on the determination of the entropy to include Jahn-Teller as well as
electronic and local disproportionation.

It would be very interesting to look for experimental signatures of the glass
transition in LiNiO2, at which the electronic and orbital fluctuations
discussed above should freeze. For example, it was found that the heat
capacity of LiNiO2 is anomalously enhanced below 300 KKawaji et al. (2002).
Clearly, further heat capacity measurements exploring much higher temperatures
are strongly desired in order to observe the full glass transition. Another
experimental observation which is strongly supportive of the glassy state in
LiNiO2 is the spin-glass like transition observed at around 20 K. Indeed,
having the Ni sites electronically in a random mixture of
$d^{7}\underline{L}$, $d^{8}\underline{L}^{2}$, and $d^{8}$ configurations
means having a random mixture of $S=\frac{1}{2}$, $S=0$, and $S=1$ spin
magnetic moments, respectively, with similarly randomized exchange
interactions between them, which is an important prerequisite for the
formation of a spin glass.

From the application point of view, the high configurational entropy of the
glassy state in LiNiO2 is very beneficial in terms of battery operation as the
system has many ways to adjust itself locally to a missing Li ion, which is
apparently what makes LiNiO2 such a good cathode. In this regard, further
increase of entropy upon chemical substitution of Ni with Co and Mn may be the
main reason for the increased stability of LiNi1-x-yCoxMnyO2 cathodes. Another
example of this concept can be found in a recent study on multi-component
monoxides which demonstrates that five-component systems perform much better
in batteries compared with four- or less component systems and attributes this
effect to entropy stabilizationSarkar et al. (2018).

###  IV.2 Oxygen holes

Finally, in order to better understand the microscopic nature of the
disordered glassy state in LiNiO2, let us examine its electronic structure. We
first note that LDA+U correctly describes LiNiO2 as a negative charge-transfer
systemKuiper et al. (1989), as far as the Ni oxidation state is concerned.
Figure 4(b) compares the Ni-$3d$ and O-$2p$ orbitally projected densities of
states (DOS) in the zigzag JT phase. One can see a considerable amount of the
O-$2p$ character in the states above the Fermi level, indicating that the true
Ni oxidation state is far less than 3+. We find close to 8 electrons in the
Ni-$3d$ shell inside the muffin-tin sphere. What obscures the picture, though,
is the strong hybridization between the Ni-$e_{g}$ orbitals and the
oxygen-$2p$ orbitals of a respective symmetry. As one can see in the right
panel of Fig. 4(b), it results in holes occupying $E_{g}$-symmetric molecular-
like orbitals formed by the oxygen-$2p_{\sigma}$ orbitals in an octahedral
cage. Although LDA+U is a mean-field method not capable of properly describing
a many-body wave-function, we still can associate the states observed in our
calculations on the expanded, JT distorted, and collapsed octahedra with,
respectively, the $t_{2g}^{6}e_{g}^{2}$, $(t_{2g}^{6}e_{g}^{2})\underline{L}$,
and $(t_{2g}^{6}e_{g}^{2})\underline{L}^{2}$ configurations. Hole density iso-
surfaces for each configuration are schematically shown in Fig. 4(a). Note
that a single hole occupying the $x^{2}-y^{2}$-symmetric $E_{g}$ molecular
orbital results in a shortening of only four Ni-O bonds, while two holes
occupying both the $x^{2}-y^{2}$\- and the $z^{2}$-symmetric $E_{g}$ orbitals
result in a collapse of all Ni-O bonds. We further find 1.68, 0.90, and 0.08
$\mu_{\text{B}}$ for the magnetic moments of the Ni ions inside the expanded,
JT distorted, and collapsed octahedra, respectively. Again, this is a result
of LDA+U being a mean-field method; in a true many-body calculation, one would
find for $(t_{2g}^{6}e_{g}^{2})\underline{L}^{2}$, for example, that the $S=1$
moment on the nickel forms a singlet state with the $S=1$ moment on the oxygen
molecular orbital and the total moment on the octahedron is zero.

Figure 5: Oxygen charge disproportionation in LiNiO2. (a) A schematic
illustration of the density iso-surfaces of holes occupying individual oxygen
molecular orbitals in the 67% SD phase. (b) Oxygen charge disproportionation
in the 67% SD phase versus a homogeneous hole distribution in the zigzag JT
phase.

Yet another peculiar property of the SD phases in LiNiO2 is that there occur
both bond disproportionation and oxygen charge disproportionation. This is in
contrast with the rare-earth nickel perovskites, where only bond
disproportionation occurs, which is a result of a different lattice geometry
and a different average number of holes per oxygen octahedron (three in LiNiO2
versus two in the rare-earth nickel perovskites). In Fig. 5(a), we
schematically show density iso-surfaces of holes occupying individual
molecular orbitals in the 67% SD phase of LiNiO2. One can clearly distinguish
oxygen sites where neighboring molecular orbitals overlap from those where
they do not. As a result, the two types of oxygen sites have different hole
concentrations, which is further illustrated in Fig. 5(b).

The extended nature of the oxygen molecular orbitals and the fact that they
can overlap result in a strong nearest-neighbor interaction between holes,
which is generally known to facilitate charge ordering effects. Although
formation of a glassy state may in general depend on various factors, we can
speculate that in such an environment localized charge impurities can serve as
nucleation centers for charge disproportionation. In LiNiO2, charge impurities
are readily available in a form of doped electrons entering the Ni layers as a
result of the unavoidable deviation from stoichiometryArai et al. (1995). The
doped electrons would attract holes according to Pauling’s principlePauling
(1948) and thus locally modify the energetics in favor of charge- (and size-)
disproportionated phases. Further investigations are, however, required in
order to validate this scenario.

We believe that spectroscopic probes, such as x-ray absorption spectroscopy
(XAS), could be crucial in verifying the mixed size-disproportionated state in
LiNiO2. Especially interesting would be to compare spectroscopic data from
LiNiO2 and its sister compound NaNiO2 where the ground state is contrastingly
uniform. Unfortunately, unambiguous interpretation of the existing Ni $2p$ XAS
dataAbbate et al. (1991); van Veenendaal and Sawatzky (1994); Montoro et al.
(1999) is not easy due to questionable sample quality and also the fact that
their proper theoretical description requires to go beyond the standard single
cluster calculationsGreen et al. (2016). We hope therefore that our
theoretical proposal of the glass-like electronic ground state of LiNiO2 will
stimulate experimental efforts in preparing and systematically measuring high-
quality LiNiO2 and NaNiO2 samples complemented with a proper theoretical
analysis of spectroscopic data in the spirit of Ref. Green et al., 2016. We
also note that in the case of the negative charge-transfer gap rare-earth
nickelates the detailed analysis of their XAS data identifying the origin of
the two-peak structure was only made after the single crystalline epitaxial
thin film studies with resonant inelastic x-ray scattering by Bisogni et
al.Bisogni et al. (2016). This kind of study on epitaxial thin films of LiNiO2
and NaNiO2 would be extremely helpful.

##  V Conclusion

In summary, we have used electronic structure methods to demonstrate that the
Li-ion battery material LiNiO2 is in a high-entropy charge-glass like state
characterized by a disordered mixture of JT distorted and SD NiO6 octaherda.
This state is associated with charge disproportionation on the oxygen rather
than on the nickel atoms due to LiNiO2 being in the negative charge-transfer
regime. Supported by rigorous calculations, our proposal explains the
previously highly debated nPDF and powder diffraction measurements on LiNiO2
as well as its enhanced heat capacity and the spin-glass like behavior. Most
importantly, however, we can conclude that it is the glassy nature that
renders LiNiO2 its extraordinarily high reversible capacity. We believe that
recognizing the role of entropy in stabilizing cathodes during the charge and
discharge cycles can greatly advance battery research.

## Appendix

###  V.1 Dependence of energy on a spin configuration

Figure 6 demonstrates that the energy differences between different structural
phases of LiNiO2 are very weakly dependent on the Ni spin configuration. In
this plot, we have considered the zigzag JT, ferro-orbital JT, and fully size-
disproportionated structural phases as well as one ferromagnetic and three
anti-ferromagnetic Ni spin configurations labeled as FM, AFM1, AFM2, and AFM3.

###  V.2 Calculation of thermodynamic properties

In order to estimate how important is the impact of the vibrational free
energy on the stability of various LiNiO2 phases, we have computed the
vibrational free energies $F^{\text{phon}}=E^{\text{phon}}-TS^{\text{phon}}$
of the zigzag JT and the 67% SD phases of LiNiO2 using the phonon spectra
obtained with VASP and Phonopy via the dynamical matrix method. The computed
vibrational entropies $S^{\text{phon}}$, the harmonic phonon energies
$E^{\text{phon}}$, as well as the total entropies and total free energies
$F=E^{T=0}+F^{\text{phon}}-TS^{\text{conf}}$ are shown in Fig. 7 as a function
of temperature $T$. $E^{T=0}$ is the $T=0$ energy, which is 0 and +32.1
meV/f.u. in the zigzag JT and the 67% SD phases, respectively. It is assumed
that the configurational entropy of the 67% SD phase is that of the entropy-
stabilized state with equal concentrations of $d^{8}$, $d^{8}\underline{L}$,
and $d^{8}\underline{L}^{2}$ sites in order to illustrate its significant
contribution compared with the phonon contribution.

Figure 6: Relative total energies as a function of a structural phase and Ni
spin configuration in LiNiO2. FM stands for the ferromagnetic spin
configuration. The anti-ferromagnetic configuration AFM1 corresponds to
alternating layers of up and down spins. The anti-ferromagnetic configurations
AFM2 and AFM3 correspond to two different arrangements of up and down spins
within a layer. Note that AFM2 and AFM3 are not equivalent due to the broken
$C_{3}$ symmetry in the considered structural phases. Figure 7:  Temperature
dependence of the thermodynamic properties of LiNiO2 in the zigzag Jahn-Teller
(zz JT) and the 67% size-disproportiotaled (67% SD) phases: (a)
configurational, vibrational, and total entropies $S$; (b) harmonic phonon
energy $E^{\text{phon}}$; (c) total free energy $F$.

## Acknowledgments

The authors are grateful to Daniel Khomskii for insightful discussions. K. F.
thanks Jae-Ho Chung and Thomas E. Proffen for sharing their neutron data and
Leopoldo Suescun for helping with GSAS simulations.

## References

  * Goodenough (1955) J. B. Goodenough, Phys. Rev. 100, 564 (1955).
  * Goodenough (1958) J. B. Goodenough, J. Phys. Chem. Solids 6, 287 (1958).
  * Liu et al. (2015) W. Liu, P. Oh, X. Liu, M.-J. Lee, W. Cho, S. Chae, Y. Kim, and J. Cho, Angewandte Chemie International Edition 54, 4440 (2015), ISSN 1521-3773, URL http://dx.doi.org/10.1002/anie.201409262.
  * Arai et al. (1995) H. Arai, S. Okada, H. Ohtsuka, M. Ichimura, and J. Yamaki, Solid State Ionics 80, 261 (1995), ISSN 0167-2738.
  * Hirakawa et al. (1985) K. Hirakawa, H. Kadowaki, and K. Ubukoshi, Journal of the Physical Society of Japan 54, 3526 (1985), eprint http://dx.doi.org/10.1143/JPSJ.54.3526, URL http://dx.doi.org/10.1143/JPSJ.54.3526.
  * Hirota et al. (1991) K. Hirota, Y. Nakazawa, and M. Ishikawa, Journal of Physics: Condensed Matter 3, 4721 (1991), URL http://stacks.iop.org/0953-8984/3/i=25/a=017.
  * Yamaura et al. (1996) K. Yamaura, M. Takano, A. Hirano, and R. Kanno, Journal of Solid State Chemistry 127, 109 (1996), ISSN 0022-4596, URL http://www.sciencedirect.com/science/article/pii/S0022459696903631.
  * Chung et al. (2005) J.-H. Chung, T. Proffen, S. Shamoto, A. M. Ghorayeb, L. Croguennec, W. Tian, B. C. Sales, R. Jin, D. Mandrus, and T. Egami, Phys. Rev. B 71, 064410 (2005), URL http://link.aps.org/doi/10.1103/PhysRevB.71.064410.
  * Borgers and Enz (1966) P. Borgers and U. Enz, Solid State Communications 4, 153 (1966), ISSN 0038-1098, URL http://www.sciencedirect.com/science/article/pii/0038109866900019.
  * Holzapfel et al. (2004) M. Holzapfel, S. de Brion, C. Darie, P. Bordet, E. Chappel, G. Chouteau, P. Strobel, A. Sulpice, and M. D. Núñez Regueiro, Phys. Rev. B 70, 132410 (2004), URL http://link.aps.org/doi/10.1103/PhysRevB.70.132410.
  * Rougier et al. (1995) A. Rougier, C. Delmas, and A. Chadwick, Solid State Communications 94, 123 (1995), ISSN 0038-1098, URL http://www.sciencedirect.com/science/article/pii/0038109895000208.
  * Barra et al. (1999) A.-L. Barra, G. Chouteau, A. Stepanov, A. Rougier, and C. Delmas, The European Physical Journal B - Condensed Matter and Complex Systems 7, 551 (1999), ISSN 1434-6036, URL https://doi.org/10.1007/s100510050648.
  * Reynaud et al. (2001) F. Reynaud, D. Mertz, F. Celestini, J.-M. Debierre, A. M. Ghorayeb, P. Simon, A. Stepanov, J. Voiron, and C. Delmas, Phys. Rev. Lett. 86, 3638 (2001), URL https://link.aps.org/doi/10.1103/PhysRevLett.86.3638.
  * Petit et al. (2006) L. Petit, G. M. Stocks, T. Egami, Z. Szotek, and W. M. Temmerman, Phys. Rev. Lett. 97, 146405 (2006), URL https://link.aps.org/doi/10.1103/PhysRevLett.97.146405.
  * Kuiper et al. (1989) P. Kuiper, G. Kruizinga, J. Ghijsen, G. A. Sawatzky, and H. Verweij, Phys. Rev. Lett. 62, 221 (1989), URL http://link.aps.org/doi/10.1103/PhysRevLett.62.221.
  * Johnston et al. (2014) S. Johnston, A. Mukherjee, I. Elfimov, M. Berciu, and G. A. Sawatzky, Phys. Rev. Lett. 112, 106404 (2014), URL http://link.aps.org/doi/10.1103/PhysRevLett.112.106404.
  * Green et al. (2016) R. J. Green, M. W. Haverkort, and G. A. Sawatzky, Phys. Rev. B 94, 195127 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.195127.
  * Foyevtsova et al. (2015) K. Foyevtsova, A. Khazraie, I. Elfimov, and G. A. Sawatzky, Phys. Rev. B 91, 121114 (2015), URL http://link.aps.org/doi/10.1103/PhysRevB.91.121114.
  * Chen et al. (2011) H. Chen, C. L. Freeman, and J. H. Harding, Phys. Rev. B 84, 085108 (2011), URL http://link.aps.org/doi/10.1103/PhysRevB.84.085108.
  * Rost et al. (2015) C. M. Rost, E. Sachet, T. Borman, A. Moballegh, E. C. Dickey, D. Hou, J. L. Jones, S. Curtarolo, and J.-P. Maria, Nature Communications 6, 8485 (2015).
  * Yao et al. (2018) Y. Yao, Z. Huang, P. Xie, S. D. Lacey, R. J. Jacob, H. Xie, F. Chen, A. Nie, T. Pu, M. Rehwoldt, et al., Science 359, 1489 (2018), ISSN 0036-8075, eprint http://science.sciencemag.org/content/359/6383/1489.full.pdf, URL http://science.sciencemag.org/content/359/6383/1489.
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Computational Materials Science 6, 15 (1996), ISSN 0927-0256, URL http://www.sciencedirect.com/science/article/pii/0927025696000080.
  * Paier et al. (2005) J. Paier, R. Hirschl, M. Marsman, and G. Kresse, The Journal of Chemical Physics 122, 234102 (2005), URL http://scitation.aip.org/content/aip/journal/jcp/122/23/10.1063/1.1926272.
  * Blaha et al. (2001) P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, and J. Luitz, _WIEN2K, An Augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties_ (Karlheinz Schwarz, Techn. Universität Wien, Austria, 2001).
  * Perdew and Wang (1992) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992), URL http://link.aps.org/doi/10.1103/PhysRevB.45.13244.
  * Anisimov et al. (1991) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B 44, 943 (1991), URL http://link.aps.org/doi/10.1103/PhysRevB.44.943.
  * Anisimov et al. (1993) V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk, and G. A. Sawatzky, Phys. Rev. B 48, 16929 (1993), URL http://link.aps.org/doi/10.1103/PhysRevB.48.16929.
  * Liechtenstein et al. (1995) A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467 (1995), URL http://link.aps.org/doi/10.1103/PhysRevB.52.R5467.
  * Farrow et al. (2007) C. L. Farrow, P. Juhas, J. W. Liu, D. Bryndin, E. S. Božin, J. Bloch, T. Proffen, and S. J. L. Billinge, Journal of Physics: Condensed Matter 19, 335219 (2007), URL http://stacks.iop.org/0953-8984/19/i=33/a=335219.
  * Larson and Dreele (2000) A. Larson and R. V. Dreele, Los Alamos National Laboratory Report LAUR 86, 748 (2000).
  * Toby (2001) B. H. Toby, J. Appl. Cryst. 34, 210 (2001).
  * Togo and Tanaka (2015) A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * Mostovoy and Khomskii (2002) M. V. Mostovoy and D. I. Khomskii, Phys. Rev. Lett. 89, 227203 (2002), URL http://link.aps.org/doi/10.1103/PhysRevLett.89.227203.
  * Meskine and Satpathy (2005) H. Meskine and S. Satpathy, Journal of Applied Physics 97, 10A314 (2005), eprint https://doi.org/10.1063/1.1854414, URL https://doi.org/10.1063/1.1854414.
  * (35) For instance, such is a JT phase where the direction of the long Ni-O bond in an octahedron alternates between successive layers.
  * Pauling (1948) L. Pauling, J. Chem. Soc. pp. 1461–1467 (1948), URL http://dx.doi.org/10.1039/JR9480001461.
  * Murty et al. (2014) B. S. Murty, J.-W. Yeh, and S. Ranganathan, _High-entropy alloys_ (Butterworth-Heinemann, 2014).
  * Sarkar et al. (2018) A. Sarkar, L. Velasco, D. Wang, Q. Wang, G. Talasila, L. de Biasi, C. Kübel, T. Brezesinski, S. S. H. Bhattacharya, H. H. Hahn, et al., Nature Communications 9, 3400 (2018).
  * Kawaji et al. (2002) H. Kawaji, T. Oka, T. Tojo, T. Atake, A. Hirano, and R. Kanno, Solid State Ionics 152-153, 195 (2002), ISSN 0167-2738, URL http://www.sciencedirect.com/science/article/pii/S0167273802003004.
  * Abbate et al. (1991) M. Abbate, F. M. F. de Groot, J. C. Fuggle, A. Fujimori, Y. Tokura, Y. Fujishima, O. Strebel, M. Domke, G. Kaindl, J. van Elp, et al., Phys. Rev. B 44, 5419 (1991), URL https://link.aps.org/doi/10.1103/PhysRevB.44.5419.
  * van Veenendaal and Sawatzky (1994) M. A. van Veenendaal and G. A. Sawatzky, Phys. Rev. B 50, 11326 (1994), URL https://link.aps.org/doi/10.1103/PhysRevB.50.11326.
  * Montoro et al. (1999) L. Montoro, M. Abbate, E. Almeida, and J. Rosolen, Chemical Physics Letters 309, 14 (1999), ISSN 0009-2614, URL http://www.sciencedirect.com/science/article/pii/S0009261499006508.
  * Bisogni et al. (2016) V. Bisogni, S. Catalano, R. J. Green, M. Gibert, S. Raoul, Y. Huang, V. N. Strocov, P. Zubko, S. Balandeh, J.-M. Triscone, et al., Nature Communications 7, 13017 (2016).