# A fast solver for the pseudo-two-dimensional model of lithium-ion batteries

Rachel Han<EMAIL_ADDRESS>Colin Macdonald<EMAIL_ADDRESS>Brian Wetton
<EMAIL_ADDRESS>
###### Abstract

The pseudo-two-dimensional (P2D) model is a complex mathematical model that
can capture the electrochemical processes in Li-ion batteries. However, the
model also brings a heavy computational burden. Many simplifications to the
model have been introduced in the literature to reduce the complexity. We
present a method for fast computation of the P2D model which can be used when
simplifications are not accurate enough. By rearranging the calculations, we
reduce the complexity of the linear algebra problem. We also employ automatic
differentiation, using an open source package JAX for robustness, while also
allowing easy implementation of changes to coefficient expressions. The method
alleviates the computational bottleneck in P2D models without compromising
accuracy.

###### keywords:

Li-ion battery , P2D model , finite differences , automatic-differentiation

[description=Electrolyte salt concentration [mol/m3]]$u(x,t)$$u(x,t)$
[description=Solid-phase surface concentration
[mol/m3]]$c^{*}(x,t)$$c^{*}(x,t)$ [description=Electrolyte diffusivity
[m2/s]]$D_{i}$$D_{i}$ [description=Solid-phase diffusivity
[m2/s]]$D^{s}_{i}$$D^{s}_{i}$ [description=Reaction rate
[m2.5/(mol0.5s)]]$k_{i}$$k_{i}$ [description=Thickness [m]]$l_{i}$$l_{i}$
[description=Particle radius [m]]$R_{p,i}$$R_{p,i}$ [description=Density
[kg/m3]]rho$\rho_{i}$ [description=Specific heat [J/(kg$\cdot$K)
]]$C_{p}$$C_{p,i}$ [description=Thermal conductivity
[W/(m$\cdot$K)]]lambda$\lambda_{i}$ [description=Solid-phase conductivity
[S/m]]sigma$\sigma_{i}$ [description=Porosity [-]]epsilon$\epsilon_{i}$
[description=Particle surface area to volume [m2/m3]]a$a_{i}$
[description=Reaction constant activation energy [J/mol]]Ek$E^{k_{i}}_{a}$
[description=Solid-phase diffusion activation energy
[J/mol]]Ek$E^{D^{s}_{i}}_{a}$ [description=Bruggeman’s coefficient
[-]]bruggbrugg [description=Faraday’s constant [C/mol]]F$F$
[description=Universal gas constant [J/(mol$\cdot$K)]]R$R$
[description=Porosity [-]]epsilon$\epsilon_{i}$ [description=Transference
number [-]]t$t_{+}$ [description=Filler fraction [-]]epsilonf$\epsilon_{f,i}$

[1]organization=University of British Columbia, addressline=6356 Agricultural
Road, postcode=V6T 1Z2, city=Vancouver, country=Canada

##  1 Introduction

Lithium-ion (Li-ion) batteries are essential in modern energy storage and
widely used in devices from portable electronics to large electric vehicles.
Li-ion batteries have high energy efficiency and high power density compared
to their predecessors such as lead-acid and zinc-carbon batteries [1]. They
garnered great research interest since their commercial inception in the 1980s
[2].

However, the battery needs to be monitored for control and optimization
purposes, especially for complex systems. This is due to the nonlinear
internal behaviour in the battery that is difficult to capture with external
observations. The battery management system (BMS) plays a key part in this
purpose. Its functions include predicting the state of charge (SOC) (the level
of charge given the battery capacity), the state of health (SOH) (a
qualitative measure of how much the battery has degraded from its original
specifications) and the remaining useful life (RUL) [3]. These characteristics
cannot be directly measured and can only be inferred by models.

The battery management system battery models can largely be classified into
two categories: equivalent circuit models and electrochemical models. In more
recent literature, data-driven methods such as deep learning have been
developed that try to predict internal battery behaviour, as done in [4], [5]
and [6]. We focus on a particular electrochemical model called the pseudo-two-
dimensional (P2D) model, also known as the Doyle-Fuller-Newman model [7]. It
is a full physics based model described by partial differential algebraic
equations. This complex model allows for a more accurate description of the
battery but also brings a heavy computational burden on a practical BMS. In
the literature, there are many electrochemical models with reduced complexity
such as the single particle model [8], [9] and approximations of solid-phase
diffusion [10], [11]. The simplified models are used in applications such as
control, on-line monitoring, optimization, parameter estimation and age
prediction [12]. In this work, we present a method for fast computation of the
P2D model which can be used when the simplifications are not sufficiently
accurate.

The nonlinear partial differential algebraic equations (PDAE) that describe
the P2D model are often implemented by Finite Difference Methods (FDM), Finite
Volume Methods (FVM) and Finite Elements Methods (FEM) [13]. Commercial
software such as COMSOL Inc. Multiphysics software [14] has FDM and FEM
implementations for the P2D model. LIONSIMBA [15] uses FVM to implement the
P2D model in MATLAB with reductions made in the solid-phase diffusion. It is
typical in literature ([16], [15]) to use black box differential algebraic
solvers, such as the DASSL solver [17] and SUNDIALS [18] for time stepping.

In this paper, we develop a fast solver for the thermal P2D model taken from
[15] using FDM, backward Euler and Newton’s method in Python. Additionally, we
also employ automatic differentiation, using an open-source package JAX [19].
This framework allows terms in the model to be changed easily, such as the
diffusive and conductivity expressions. The method alleviates the
computational bottleneck inherent in P2D models, without compromising
accuracy.

##  2 Thermal P2D model equations

The pseudo-two-dimensional (P2D) model is an electrochemical model that
captures the kinetics, transport processes and thermodynamics of a lithium-ion
battery [7]. The name pseudo-2D comes from the interesting simplification of
the battery geometry. In the $x$ dimension, the processes in electrolyte phase
across the battery are described. The second dimension, $r$, comes from the
following assumption: in the porous electrodes, the host intercalation
materials such as lithium cobalt oxide and graphite are represented by
spherical particles at each channel location in $x$, where $r$ is in the
direction normal to the surface of each particle. Moreover, $x$ is of
macroscopic scale whereas $r$ is microscopic. P2D is not a strict two-
dimensional model but rather a pseudo-2D model because the two dimensions with
different scales are coupled. Figure 2 illustrates the simplified geometry of
the P2D model.

Figure 1: A representation of chemical properties of lithium-ion battery with
lithium cobalt oxide and graphite. The purple spheres represent Lithium (Li)
that are intercalated in cathode and anode active materials. Figure 2: The P2D
model abstraction of Figure 1

The P2D model is described by the following quantities: lithium-ion
concentrations, potentials, ionic flux across the interface of the spherical
particle, current density and temperature. These quantities describe the
transport of lithium ion species, conservation of charge, kinetics of charge
transfer and thermal behaviour of the battery. The lithium-ion concentrations
are denoted by $u(x,t)$ in the electrolyte (liquid phase) and intercalated
lithium concentration $c(x,r,t)$ in the solid particle (solid phase). The
electric potentials in the liquid and solid phase are $\Phi(x,t)$ and
$\Psi(x,t)$ respectively. The interfacial ionic flux is denoted by $j(x,t)$,
$I$ is the current density and $T(x,t)$ is the temperature. We take the model
form and coefficients from [15]. The governing equations form a PDAE system,
where the system is described by partial differential equations and algebraic
constraints. We present the form of the equations in detail below, so that the
reader can see how we exploit the structure of the nonlinear PDAE to obtain
the fast solver.

Throughout the paper, we will use the notation $i\in\\{a,p,s,n,z\\}$ in
subscript to denote the battery section: the positive current collector,
cathode, separator, anode and the negative current collector respectively. The
ends of the battery are at $x=0$ and $x=L$, where $L$ is the length of the
battery. Figure 1 summarizes the domain. The subscript “eff” is used to denote
the effective coefficients based on Bruggeman theory [20], which accounts the
for varying conductivity or diffusivity due to different materials in
composite materials. The model parameters, coefficients and symbols are taken
from [15] and also listed in A and B.

###  2.1 Electrodes and Separator

The intercalation of Li+ in the solid particles drive the charge and discharge
in Li-ion batteries. The intercalation is driven by diffusion of the lithium
concentration,

$\frac{\partial c}{\partial t}=\frac{1}{r^{2}}\dfrac{\partial}{\partial
r}\left(r^{2}D^{s}\frac{\partial c(r,t)}{\partial r}\right),$ (2.1a)
$\displaystyle\frac{\partial c(r,t)}{\partial
r}\bigg{\rvert}_{r=0}=0,\qquad\frac{\partial c(r,t)}{\partial
r}\bigg{\rvert}_{r=R_{p}}=-\frac{j(x,t)}{D^{s}_{\text{eff},i}(T(x,t))}.$
(2.1b)

Here, $D^{s}$ is a constant solid-phase diffusion coefficient and
$D^{s}_{\text{eff},i}$ for $i\in\\{p,n\\}$ is the effective solid-phase
diffusion coefficient which depends on temperature in $x$. The kinetics of
charge transfer at the electrodes are described by the Butler-Volmer equation,

$j(x,t)=2k_{\text{eff}}\sqrt{u(x,t)(c^{\text{max}}_{s,i}-c^{*}(x,t))c^{*}(x,t)}\sinh\left[\dfrac{0.5F}{RT}\eta_{i}(x,t)\right]$
(2.2)

for $i\in\\{p,n\\}$, where $\eta_{i}(x,t)$ denotes the overpotential. The
equation (2.2) is an algebraic constraint in the PDAE, along with the equation
for $\eta_{i}(x,t)$ in (2.10). In the porous electrodes $i\in\\{p,n\\}$, the
accumulation of Li+ in the electrolyte phase is described in terms of
diffusion of the ions and $j$:

$\epsilon_{i}\frac{\partial{u(x,t)}}{\partial t}=\frac{\partial}{\partial
x}[\mathbf{D_{eff,i}}\frac{\partial u(x,t)}{\partial
x}]+a_{i}(1-t_{+})j(x,t),$ (2.3a) with equal flux conditions enforced at the
boundaries, $\displaystyle\frac{\partial u(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}_{0},\hat{x}_{n}}=0,$ (2.3b)
$\displaystyle-\mathbf{D_{eff,p}}\frac{\partial u(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}^{-}_{p}}=-\mathbf{D_{eff,s}}\frac{\partial
u(x,t)}{\partial x}\bigg{\rvert}_{x=\hat{x}^{+}_{p}},$ (2.3c)
$\displaystyle-\mathbf{D_{eff,s}}\frac{\partial u(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}^{-}_{s}}=-\mathbf{D_{eff,n}}\frac{\partial
u(x,t)}{\partial x}\bigg{\rvert}_{x=\hat{x}^{+}_{s}}.$ (2.3d)

The coefficient $\bm{D}_{\text{eff},i}$ is the effective diffusive coefficient
of Li+ in the electrolyte. The interface positions $\hat{x}_{i}$ are shown in
Figure 1.

In the separator where there are no solid particles, we have $j=0$ which
yields

$\epsilon_{i}\frac{\partial u(x,t)}{\partial t}=\frac{\partial}{\partial
x}\left[D_{\text{eff},s}\frac{\partial u(x,t)}{\partial x}\right],$ (2.4a)
with boundary conditions $\displaystyle-\bm{D}_{\text{eff},p}\frac{\partial
u(x,t)}{\partial
x}\rvert_{x=\hat{x}^{-}_{p}}=-\bm{D}_{\text{eff},s}\frac{\partial
u(x,t)}{\partial x}\rvert_{x=\hat{x}^{+}_{p}},$ (2.4b)
$\displaystyle-\bm{D}_{\text{eff},s}\frac{\partial u(x,t)}{\partial
x}\rvert_{x=\hat{x}^{-}_{s}}=-\bm{D}_{\text{eff},n}\frac{\partial
u(x,t)}{\partial x}\rvert_{x=\hat{x}^{+}_{s}}.$ (2.4c)

In the electrodes $i\in\\{p,n\\}$, the transport of Li+ is described by
$\Phi$, $T$, $u$ and $j$:

$a_{i}Fj(x,t)=-\frac{\partial}{\partial
x}\left[\mathbf{\kappa_{eff,i}}\frac{\partial\Phi(x,t)}{\partial
x}\right]+\frac{\partial}{\partial x}\left[\mathbf{\kappa_{eff,i}}\Upsilon
T(x,t)\frac{\partial\ln u(x,t)}{\partial x}\right],$ (2.5a) with boundary
conditions $\displaystyle\frac{\partial\Phi(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}_{0}}=0,$ (2.5b)
$\displaystyle-\mathbf{\kappa_{eff,p}}\frac{\partial\Phi(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}_{p}^{-}}=-\mathbf{\kappa_{eff,s}}\frac{\partial\Phi(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}_{p}^{+}},$ (2.5c)
$\displaystyle-\mathbf{\kappa_{eff,s}}\frac{\partial\Phi(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}_{s}^{-}}=-\mathbf{\kappa_{eff,n}}\frac{\partial\Phi(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}_{s}^{+}},$ (2.5d)
$\displaystyle\Phi(x,t)\bigg{\rvert}_{x=\hat{x}_{n}}=0,$ (2.5e) where
$\bm{\kappa}_{\text{eff},i}$ is the effective electrolyte phase conductivity.

Again, in the separator, the ionic flux $j$ is $0$. Therefore, the ionic
transport in the separator is described by:

$0=-\frac{\partial}{\partial
x}\left[\bm{\kappa}_{\text{eff},s}\frac{\partial\Phi(x,t)}{\partial
x}\right]+\frac{\partial}{\partial
x}\left[\bm{\kappa}_{\text{eff},s}T(x,t)\Upsilon\frac{\partial\ln
u(x,t)}{\partial x}\right]$ (2.6a) with boundary conditions
$\displaystyle-\bm{\kappa}_{\text{eff},p}\frac{\partial\Phi(x,t)}{\partial
x}\rvert_{x=\hat{x}^{-}_{p}}=-\bm{\kappa}_{\text{eff},s}\frac{\partial\Phi(x,t)}{\partial
x}\rvert_{x=\hat{x}^{+}_{p}},$ (2.6b)
$\displaystyle-\bm{\kappa}_{\text{eff},s}\frac{\partial\Phi(x,t)}{\partial
x}\rvert_{x=\hat{x}^{-}_{s}}=-\bm{\kappa}_{\text{eff},n}\frac{\partial\Phi(x,t)}{\partial
x}\rvert_{x=\hat{x}^{+}_{s}}.$ (2.6c)

The movement of electrons in the electrodes is governed by Ohm’s law [21],
which involves the solid phase potential $\Psi$,

$\displaystyle\frac{\partial}{\partial
x}\left[\mathbf{\sigma_{eff,i}}\frac{\partial\Psi(x,t)}{\partial
x}\right]=a_{i}Fj(x,t),$ (2.7a)
$\displaystyle\mathbf{\sigma_{eff,i}}\frac{\partial\Psi(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}_{0},\hat{x}_{n}}=-I_{\text{app}}(t),\qquad\mathbf{\sigma_{eff,i}}\frac{\partial\Psi(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}_{p},\hat{x}_{s}}=0.$ (2.7b)

For each electrode $i\in\\{p,n\\}$, $\sigma_{\text{eff},i}$ is the effective
solid-phase conductivity and $I_{\text{app}}$ is the applied current.

The heat created in the battery is described by:

$\rho_{i}C_{p,i}\frac{\partial T(x,t)}{\partial t}=\frac{\partial}{\partial
x}\left[\lambda_{i}\frac{\partial T(x,t)}{\partial
x}\right]+\bm{Q}_{\text{ohm}}+\bm{Q}_{\text{rxn}}+\bm{Q}_{\text{rev}},$ (2.8a)
which includes different heat sources. As mentioned in [15], the ohmic
generation rate $\bm{Q}_{\text{ohm}}$ takes account of heat generated from the
movement of Li+, the reaction generation rate $\bm{Q}_{\text{rxn}}$ describes
the heat resulting from the ionic flux and overpotentials, and finally the
reversible generation rate $\bm{Q}_{\text{rev}}$ accounts for the change of
entropy in the electrodes. Details of these expressions can be found in the A.
The thermal flux continuity interface conditions for temperature are enforced,
$\displaystyle-\lambda_{z}\frac{\partial T(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}_{0}^{-}}=-\lambda_{p}\frac{\partial
T(x,t)}{\partial x}\bigg{\rvert}_{x=\hat{x}_{0}^{+}},$ (2.8b)
$\displaystyle-\lambda_{n}\frac{\partial T(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}_{n}^{-}}=\lambda_{z}\frac{\partial T(x,t)}{\partial
x}\bigg{\rvert}_{x=\hat{x}_{n}^{+}}$ (2.8c)

Similarly, the temperature effects in the separator is described but without
$\bm{Q}_{\text{rxn}}$ and $\bm{Q}_{\text{rev}}$:

$\rho_{i}C_{p,i}\frac{\partial T(x,t)}{\partial t}=\frac{\partial}{\partial
x}\left[\lambda_{i}\frac{\partial T(x,t)}{\partial x}\right]+\bm{Q}_{ohm},$
(2.9a) $\displaystyle-\lambda_{p}\frac{\partial T(x,t)}{\partial
x}\rvert_{x=\hat{x}^{-}_{p}}=-\lambda_{s}\frac{\partial T(x,t)}{\partial
x}\rvert_{x=\hat{x}^{+}_{p}},$ (2.9b) $\displaystyle-\lambda_{s}\frac{\partial
T(x,t)}{\partial x}\rvert_{x=\hat{x}^{-}_{s}}=-\lambda_{n}\frac{\partial
T(x,t)}{\partial x}\rvert_{x=\hat{x}^{+}_{s}}.$ (2.9c)

Lastly, we give an equation for the surface overpotential $\eta(x,t)$ of the
electrode $i\in\\{p,n\\}$:

$\eta_{i}(x,t)=\Psi(x,t)-\Phi(x,t)-\mathbf{U_{i}}(T(x,t),c(x,r,t))$ (2.10)

where $\mathbf{U_{i}}$ is the open circuit potential which is fitted
experimentally (see A).

###  2.2 Current Collectors

The current collectors consist of conductive metal materials for transfer of
electrons. We model the thermal effects in these regions:

$\displaystyle\rho_{i}C_{p,i}\frac{\partial T(x,t)}{\partial
t}=\frac{\partial}{\partial x}\left(\lambda_{i}\frac{\partial T(x,t)}{\partial
x}\right)+\frac{I^{2}_{app}(t)}{\sigma_{\text{eff},i}},$ (2.11a) with boundary
conditions $\displaystyle-\lambda_{a}\frac{\partial T(x,t)}{\partial
x}\rvert_{x=0}$ $\displaystyle=h(T_{ref}-T(x,t)),$ (2.11b)
$\displaystyle-\lambda_{z}\frac{\partial T(x,t)}{\partial x}\rvert_{x=L}$
$\displaystyle=h(T(x,t)-T_{ref}),$ (2.11c)

with heat exchange coefficient $h$. The boundary conditions represent Newton’s
law of cooling at the two ends of the battery, with ambient temperature
$T_{ref}$.

##  3 Implementation

The system of equations in Section 2 is discretised in space by cell centred
finite differences, and in time using backward Euler method for the time-
dependent variables, $u(x,t)$, $c(x,r,t)$ and $T(x,t)$. The complete
discretizations are given in Tables 1 and 2. We describe the structure of the
discrete nonlinear system below.

###  3.1 Solving the nonlinear system

We first describe how to solve the resulting nonlinear system from the PDAE
without considering the specific structure of the problem. We will call this
basic implementation as naïve method. After discretising in space, we obtain a
set of differential algebraic equations (DAE). In particular, we obtain a
semi-explicit index-1 DAE, where one can distinguish the differential and
algebraic variables. Variables without time derivatives are called algebraic
as in standard PDAE literature [22]. We denote the discrete versions of the
electrochemical quantities $c(x,r,t),j(x,t),u(x,t),\Phi(x,t),\Psi(x,t),T(x,t)$
and $\eta(x,t)$ as

$\displaystyle\bm{c}$
$\displaystyle=(c_{p},c_{n}),\quad\bm{j}=(j_{p},j_{n}),\quad\bm{u}=(u_{p},u_{s},u_{n}),$
$\displaystyle\bm{\Phi}$
$\displaystyle=(\Phi_{p},\Psi,\Phi_{n}),\quad\bm{\Psi}=(\Psi_{p},\Psi,_{n}),$
$\displaystyle\bm{T}$
$\displaystyle=(T_{a},T_{p},T_{s},T_{n},T_{z}),\quad\bm{\eta}=(\eta_{p},\eta_{n}).$

We can group the differential variables in $x$ as $\hat{U}=[\bm{u},\bm{T}]$
and algebraic variables as $\tilde{U}=[\bm{\Phi},\bm{\Psi},\bm{j},\bm{\eta}]$.
The solid particle concentration $\bm{c}$ is also a differential variable in
$r$. After discretising in time using backward Euler, we form the following
set of nonlinear equations,

$\displaystyle\bm{0}$ $\displaystyle=\bm{c}^{k}-\Delta
t\bm{G}(\bm{c}^{k},\hat{U}^{k},\tilde{U}^{k})-\bm{c}^{k-1}$ (3.1a)
$\displaystyle\bm{0}$ $\displaystyle=\hat{U}^{k}-\Delta
t\bm{\hat{F}}(\bm{c}^{k},\hat{U}^{k},\tilde{U}^{k})-\hat{U}^{k-1}$ (3.1b)
$\displaystyle\bm{0}$
$\displaystyle=\bm{\tilde{F}}(\bm{c}^{k},\hat{U}^{k},\tilde{U}^{k}),$ (3.1c)

where $k$ denotes the time step. From (3.1), we form a larger multi-
dimensional vector valued function $\bm{F}(U)$, where $U$ is the aggregate
vector containing all the unknowns,

$\begin{split}U&=[\begin{matrix}c_{p}&c_{n}&u_{p}&u_{s}&u_{n}&j_{p}&j_{n}&\eta_{p}&\eta_{n}\end{matrix}\\\
&\qquad\qquad\begin{matrix}\Phi_{p}&\Psi&\Phi_{n}&\Psi_{p}&\Psi_{n}&T_{a}&T_{p}&T_{s}&T_{n}&T_{z}\end{matrix}]^{T}.\end{split}$
(3.2)

Solving for (3.1) is equivalent to finding the root of the function

$\bm{F}(U^{k})=\bm{0}$ (3.3)

at every $k$th time step. We solve for the root of (3.3) iteratively using
Newton’s iteration,

$U^{p+1}=U^{p}-[\mathcal{J}^{-1}(U^{p})]\bm{F}(U^{p})$ (3.4)

where $p$ denotes the $p$th iterate of the Newton’s method and $\mathcal{J}$
is the Jacobian of $\bm{F}$ at $U^{p}$. Solving for the inverse of
$\mathcal{J}(U^{p})$ is computationally expensive, especially when we have a
large system. Therefore, we use a linear solver to solve for
$\delta^{p}=[\mathcal{J}^{-1}(U^{p})]\bm{F}(U^{p})$,

$\mathcal{J}(U)\delta^{p}=\bm{F}(U^{p}).$ (3.5)

We generate the iterates $U^{p+1}$ by

$U^{p+1}=U^{p}-\delta^{p}$ (3.6)

at every iteration until we converge to the root. In this way, we can
approximate the solution to PDAE by repeatedly solving large linear systems.

To give an idea of the size of the nonlinear system in this application, we
outline the number of unknowns. In each battery section $i\in\\{p,n\\}$, there
are $M_{i}$ points in $x$. For each point in $x$, there are $N_{i}$ points in
$r$. Then, the unknowns representing the solid-phase concentration $\bm{c}$
has $M_{p}(N_{p}+2)+M_{n}(N_{n}+2)$ points, where $N_{i}+2$ accounts for the
ghost points [23] needed to impose boundary conditions at $r=0$ and $r=R_{p}$.
Similarly, each of $u_{p},u_{s}$ and $u_{n}$ has $M_{p}+2$, $M_{s}+2$ and
$M_{n}+2$ unknowns respectively, and the same applies for the unknowns
$\bm{\Phi}$ and $\bm{T}$. For $\bm{j}$ and $\bm{\eta}$, there are $M_{p}$ and
$M_{n}$ points respectively in the positive and negative electrode, since
there are no boundary conditions for these variables. In total, there are

$M_{p}(N_{p}+2)+M_{n}(N_{n}+2)+3(M_{p}+M_{s}+M_{n}+6)+(M_{p}+M_{n}+4)+2(M_{p}+M_{n})$

unknowns. For example, when we set the number of grid points in each battery
section and for each solid particle to be $M=N=30$, we have $2456$ unknowns in
the resulting system. The number grows quadratically because of the dimension
$r$ of the particles. This emphasizes the need for fast and robust solver for
the full P2D model. The main contribution of this work is a fast solver
strategy to the linear systems using some of the special structure of the
discretised P2D model. For speed and robustness, we also need to consider how
to efficiently construct the Jacobians in Newton iterations.

###  3.2 Automatic-differentiation for Jacobians

The fastest method of computing the Jacobian is to manually differentiate the
function components with respect to every element in $U^{p}$ a priori. Then,
the cost of computation is just the function evaluation of the differentiated
components. This can be done for relatively simple models. However, for the
fully-coupled P2D model, manual differentiation is inconvenient due to complex
equations and prone to errors in the implementation. As alternatives, one can
use numerical, symbolic or automatic differentiation (AD). Numerical
differentiation, which involves finite differences, is easy to implement and
fast. With numerical differentiation, we must choose the step size $h$ which
is especially difficult for systems with multiple scales. We refer the readers
to [24] for the details on how round-off and truncation errors depend on the
step size $h$.

On the other hand, symbolic and automatic differentiation (AD) compute exact
Jacobians. Software for AD records every operation done on the inputs of a
function, then builds a graph that encodes all the operations [25]. Then the
rules of differentiation are applied on every node of the graph. Symbolic
differentiation produces a full expression for the derivative, then the
numerical values are substituted. This approach, however, has a problem of
expression swell causing slower performance, as pointed out by [25]. We use AD
for better performance, which uses intermediate values of computations to
calculate the derivatives on the fly. There are two ways to do this, the
forward and reverse mode AD. For functions
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ where $n>m$, the reverse mode is
preferred because the Jacobian can be constructed in $m$ passes but also uses
additional space. In contrast, the forward mode is preferred for when $m>n$
because it takes $n$ passes to build the Jacobian. When $m=n$, as is the case
for our system, the forward mode is usually preferred due to the large memory
usage with the reverse mode. However, in Section 4.2, we show how we can use
the backward mode AD to efficiently build the Jacobians. For further details
on how AD works, we refer the reader to [25]. Note that the use of
symbolic/automatic differentiation allows changes in terms in the equations
(effective diffusivity expressions, for example) with no user effort needed to
modify the Jacobian.

A particular AD software that is suitable for the Python environment is a
library named JAX [19], which we use extensively in our implementation.
Section 4.2 discusses how JAX gives an extra performance boost.

Table 1: Discretisation of PDAE from Section 2 Electrodes, $i\in\\{p,n\\}$
---
| Solid phase concentration
---
$\displaystyle
0=c^{k-1}_{n,m}-c^{k}_{n,m}+D^{s}\left[\frac{1}{r_{n}^{2}}\left(\frac{r_{n+1/2}^{2}\frac{c^{k}_{n+1,m}-c^{k}_{n,m}}{\Delta
r}-r_{n-1/2}^{2}\frac{c^{k}_{n,m}-c^{k}_{n-1,m}}{\Delta r}}{\Delta
r}\right)\right],$
$\displaystyle\qquad\frac{c^{k}_{1,m}-c^{k}_{0,m}}{\Delta
r}=0,\qquad\frac{c^{k}_{N+1,m}-c^{k}_{N,m}}{\Delta
r}+\frac{j^{k}_{m}}{D^{s}_{\text{eff}}}=0.$
| Electrolyte concentration
---
$\displaystyle 0=u^{k-1}_{m}-u^{k}_{m}-\frac{\Delta
t}{\epsilon_{i}}\left(\frac{{D_{\text{eff},m}}^{r}(u^{k}_{m+1}-u^{k}_{m})}{\Delta
x^{2}}-\frac{{D_{\text{eff},m}}^{l}(u^{k}_{m}-u^{k}_{m-1})}{\Delta
x^{2}}+a_{i}(1-t_{+})j^{k}_{m}\right),$
$\displaystyle\frac{u^{k}_{p_{1}}-u^{k}_{p_{0}}}{\Delta
x_{p}}=0,\qquad\frac{u^{k}_{n_{M_{n}+1}}-u^{k}_{n_{M_{n}}}}{\Delta x_{n}}=0$,
$\displaystyle\frac{u^{k}_{p_{M_{p}}}+u^{k}_{p_{M_{p}+1}}}{2}=\frac{u^{k}_{s_{0}}+u^{k}_{s_{1}}}{2},\qquad\frac{u^{k}_{s_{M_{s}}}+u^{k}_{s_{M_{s+1}}}}{2}=\frac{u^{k}_{n_{0}}+u^{k}_{n_{1}}}{2}$.
| Electrolyte potential
---
$\displaystyle
a_{i}Fj^{k}_{m}+(\kappa^{r}_{\text{eff}}\frac{\Phi^{k}_{m+1}-\Phi^{k}_{m}}{\Delta
x^{2}}-\kappa^{l}_{\text{eff}}\frac{\Phi^{k}_{m}-\Phi^{k}_{m-1}}{\Delta
x^{2}})-\gamma(\kappa^{r}_{\text{eff}}T^{k}_{m+1}\frac{\ln{u^{k}_{m+1}-\ln{u^{k}_{m}}}}{\Delta
x^{2}})$
$\displaystyle\frac{\Phi_{{p_{1}}}-\Phi_{{p_{0}}}}{\Delta
x}=0,\qquad\frac{\Phi_{n_{N}}+\Phi_{n_{N+1}}}{2}=0,$
$\displaystyle\frac{\Phi^{k}_{p_{M_{p}}}+\Phi^{k}_{p_{M_{p}+1}}}{2}=\frac{\Phi^{k}_{s_{0}}+\Phi^{k}_{s_{1}}}{2},\qquad\frac{\Phi^{k}_{s_{M_{s}}}+\Phi^{k}_{s_{M_{s+1}}}}{2}=\frac{\Phi^{k}_{n_{0}}+\Phi^{k}_{n_{1}}}{2}$.
| Solid phase potential
---
$\displaystyle\sigma_{\text{eff},i}\frac{\Psi^{k}_{m+1}-2\Psi^{k}_{m}+\Psi^{k}_{m-1}}{\Delta
x_{i}^{2}}=0,$
$\displaystyle\sigma_{\text{eff},p}\frac{\Psi^{k}_{p_{1}}-\Psi^{k}_{p_{0}}}{\Delta
x_{p}}=-I_{app},\qquad\sigma_{\text{eff},p}\frac{\Psi^{k}_{p_{M_{p}+1}}-\Psi^{k}_{p_{M_{p}}}}{\Delta
x_{p}}=0,$
$\displaystyle\sigma_{\text{eff},n}\frac{\Psi^{k}_{n_{1}}-\Psi^{k}_{n_{0}}}{\Delta
x_{n}}=0,\qquad\sigma_{\text{eff},n}\frac{\Psi^{k}_{n_{M_{n}+1}}-\Psi^{k}_{n_{M_{n}}}}{\Delta
x_{n}}=0.$
Table 2: Discretisation of PDAE from Section 2 Electrodes, $i\in\\{p,n\\}$
---
| Temperature
---
$\displaystyle\rho_{i}C_{p,i}\frac{T^{k}_{m}-T^{k-1}_{m}}{\Delta
t}-\left(\lambda_{i}\left[\frac{T^{k}_{m+1}-2T^{k}_{m}+T^{k}_{m-1}}{\Delta
x^{2}}\right]+Q_{ohm}+Q_{rxn}+Q_{rev}\right)=0,$
$\displaystyle-\lambda_{a}\frac{T^{k}_{a_{M_{a}+1}}-T^{k}_{a_{M_{a}}}}{\Delta
x_{a}}=-\lambda_{p}\frac{T^{k}_{p_{1}}-T^{k}_{p_{0}}}{\Delta
x_{p}},\qquad-\lambda_{n}\frac{T^{k}_{n_{M_{n}+1}}-T^{k}_{n_{M_{n}}}}{\Delta
x_{n}}=-\lambda_{z}\frac{T^{k}_{z_{1}}-T^{k}_{z_{0}}}{\Delta x_{z}}.$
$\displaystyle\frac{T^{k}_{p_{M_{p}}}+T^{k}_{p_{M_{p}+1}}}{2}=\frac{T^{k}_{s_{0}}+T^{k}_{s_{1}}}{2},\qquad\frac{T^{k}_{s_{M_{s}}}+T^{k}_{s_{M_{s+1}}}}{2}=\frac{T^{k}_{n_{0}}+T^{k}_{n_{1}}}{2}$.
| Ionic flux
---
$\displaystyle
j^{k}_{m}-2_{k_{\text{eff}},i}\sqrt{u^{k}_{m}(c^{max}-c^{*k}_{m})c^{*k}_{m}}\sinh{\frac{0.5F}{RT^{k}_{m}}\eta^{k}_{m}=0}$
| Overpotential
---
$\displaystyle\eta^{k}_{m}-\Psi^{k}_{m}+\Phi^{k}_{m}-U_{i}=0$
Current collectors, $i\in\\{a,z\\}$
| Temperature
---
$\displaystyle\rho_{i}C_{p,i}\frac{T^{k}_{m+1}-T^{k}_{m}}{\Delta
t}-\left(\lambda_{i}\left[\frac{T^{k}_{m+1}-2T^{k}_{m}+T^{k}_{m-1}}{\Delta
x^{2}}\right]+\frac{I^{2}_{app}}{\sigma_{\text{eff},i}}\right)=0,$
$\displaystyle-\lambda_{a}\frac{T^{k}_{a_{1}}-T^{k}_{a_{0}}}{\Delta
x_{a}}=h(T_{ref}-T^{k}_{a_{1/2}}),\qquad-\lambda_{z}\frac{T^{k}_{z_{M_{z}+1}}-T^{k}_{z_{M_{z}}}}{\Delta
x_{z}}=h(T^{k}_{z_{M_{z}+1/2}}-T_{ref}).$
Separator, $i=s$
| Electrolyte concentration
---
$\displaystyle 0=u^{k-1}_{m}-u^{k}_{m}-\frac{\Delta
t}{\epsilon_{i}}\left(\frac{{D_{\text{eff},m}}^{r}(u^{k}_{m+1}-u^{k}_{m})}{\Delta
x^{2}}-\frac{{D_{\text{eff},m}}^{l}(u^{k}_{m}-u^{k}_{m-1})}{\Delta
x^{2}}\right),$
$\displaystyle-D^{r}_{\text{eff},p_{M_{p}}}\frac{u_{M_{p}+1}-u_{M_{p}}}{\Delta
x_{p}}=-D^{r}_{\text{eff},s_{0}}\frac{u_{s_{1}}-u_{s_{0}}}{\Delta
x_{s}},\qquad-D^{r}_{\text{eff},s_{M_{s}}}\frac{u_{M_{s}+1}-u_{M_{s}}}{\Delta
x}=-D^{r}_{\text{eff},n_{0}}\frac{u_{n_{1}}-u_{n_{0}}}{\Delta x_{n}}.$
| Electrolyte potential
---
$\displaystyle-(\kappa^{r}_{\text{eff}}\frac{\Phi^{k}_{m+1}-\Phi^{k}_{m}}{\Delta
x^{2}}-\kappa^{l}_{\text{eff}}\frac{\Phi^{k}_{m}-\Phi^{k}_{m-1}}{\Delta
x^{2}})+\gamma(\kappa^{r}_{\text{eff}}T^{k}_{m+1}\frac{\ln{u^{k}_{m+1}-\ln{u^{k}_{m}}}}{\Delta
x^{2}})=0,$
$\displaystyle-\kappa^{r}_{\text{eff},p_{M_{p}}}\frac{\Phi_{M_{p}+1}-\Phi_{M_{p}}}{\Delta
x_{p}}=-\kappa^{r}_{\text{eff},s_{0}}\frac{\Phi_{s_{1}}-\Phi_{s_{0}}}{\Delta
x_{s}},\qquad-\kappa^{r}_{\text{eff},s_{M_{s}}}\frac{\Phi_{M_{s}+1}-\Phi_{M_{s}}}{\Delta
x_{s}}=-D^{r}_{\text{eff},n_{0}}\frac{u_{n_{1}}-u_{n_{0}}}{\Delta x_{n}}.$
| Temperature
---
$\displaystyle\rho_{i}C_{p,i}\frac{T^{k+1}_{m}-T^{k}_{m}}{\Delta
t}-\left(\lambda_{i}\left[\frac{T^{k}_{m+1}-2T^{k}_{m}+T^{k}_{m-1}}{\Delta
x^{2}}\right]+Q_{ohm}\right)=0,$
$\displaystyle-\lambda_{p}\frac{T^{k}_{p_{M_{p}+1}}-T^{k}_{p_{M_{p}}}}{\Delta
x_{p}}=-\lambda_{s}\frac{T^{k}_{s_{1}}-T^{k}_{s_{0}}}{\Delta
x_{s}},\qquad-\lambda_{n}\frac{T^{k}_{s_{1}}-T^{k}_{s_{0}}}{\Delta
x_{s}}=-\lambda_{n}\frac{T^{k}_{n_{M_{n}+1}}-T^{k}_{n_{M_{n}}}}{\Delta
x_{n}}.$

##  4 Decoupled method

As mentioned in Section 3.1, the major bottleneck of the computation comes
from solving the linear system

$\mathcal{J}(U^{p})\delta=\bm{F}(U^{p})$ (4.1)

at every Newton iteration at every time step. $\mathcal{J}$ is of size
$\mathcal{O}(NM)\times\mathcal{O}(NM)$. In this section, we present a
decoupled method to solve the linear system in two steps which results in a
smaller Jacobian matrix of size $\mathcal{O}(M)\times\mathcal{O}(M)$. We do
this by decoupling the equation for $\bm{c}$ from $\bm{F}(U)$. The technique
results in an exact solve without any approximations. This way, we can reduce
the computation cost of solving the linear system at every Newton iteration.
From (2.1) we see that equation for lithium-ion concentration $c$ is linear
except at the boundary of the sphere $r=R_{p}$. Writing the discretised (2.1)
in a matrix form, we obtain

$A\bm{c_{m}}^{k}=\begin{bmatrix}0\\\ c_{1,m}^{k-1}\\\ \vdots\\\
c_{N,m}^{k-1}\\\
0\end{bmatrix}-\frac{j_{m}^{k}}{D^{s}_{eff}(T^{k}_{m})}\begin{bmatrix}0\\\
0\\\ \vdots\\\ 0\\\ 1\end{bmatrix},$ (4.2)

where $\bm{c}^{k}_{m}$ is the particle equation at point $x_{m}$; therefore,
for each $m=1,2,\dots,M$, we would have to solve (4.2). Without the
decoupling, which will be described shortly, (4.2) is solved as part of the
larger system, as it is coupled to every other variable through $j$. $A$ is
the tridiagonal coefficient matrix for and $\bm{c}^{k}_{m}$. Let
$a_{n}=\frac{-D^{s}\Delta t}{r^{2}_{n}\Delta r^{2}}r^{2}_{n-1/2}$,
$b_{n}=1+\frac{D^{s}\Delta t}{r^{2}_{n}\Delta r^{2}}$,
$c_{n}=\frac{-D^{s}\Delta t}{r^{2}_{n}\Delta r^{2}}r^{2}_{n+1/2}$. Then,

$A=\begin{pmatrix}-1&1&0&0&0&\dots&0\\\ 0&a_{1}&b_{1}&c_{1}&0&\dots&0\\\
0&0&a_{2}&b_{2}&c_{2}&\dots&0\\\
\vdots&\vdots&\ddots&\ddots&\ddots&\ldots&\vdots\\\
0&0&\ldots&\ldots&\ldots&\frac{-1}{\Delta r}&\frac{1}{\Delta r}\end{pmatrix}.$
(4.3)

We now show how the system can be decoupled. Rearranging (4.2), we can split
$\bm{c}^{k}_{m}$ into two parts

$\displaystyle\bm{c}^{k}_{I_{m}}=A^{-1}\begin{bmatrix}0&c^{k-1}_{0,m}&\ldots&c_{N,m}^{k-1}&0\end{bmatrix}^{T},$
$\displaystyle\bm{c}^{k}_{II_{m}}=\frac{j_{m}^{k}}{D_{s}(T^{k}_{m})}A^{-1}\begin{bmatrix}0&0&\ldots&0&1\end{bmatrix}^{T}.$

The constant vector $A^{-1}\begin{bmatrix}0&0&\ldots&0&1\end{bmatrix}^{T}$ in
$\bm{c}^{k}_{II_{m}}$, which we denote $\gamma$, is a vector that can be
precomputed once. By formulating as such, we have partially decoupled the
particle equations from the whole system. This is because we can compute
$\bm{c}^{k}_{I_{m}}$ by solving

$A\bm{c}^{k}_{I_{m}}=\begin{bmatrix}0\\\ c_{1,m}^{k-1}\\\ \vdots\\\
c_{N,m}^{k-1}\\\ 0\end{bmatrix}$ (4.4)

cheaply via tridiagonal LU factorisation. Thus, we avoid the full
$\mathcal{O}(NM)\times\mathcal{O}(NM)$ Jacobian of $F(U)$. The vector
$\bm{c}^{k}_{II_{m}}$ can be obtained by scalar-vector multiplication once we
have a solution for $\bm{j}$ and $\bm{T}$. Equations that are functions of
$\bm{c}$, (2.2), (2.8) and (2.10), only require the interfacial concentration
$c^{*}(x,t)$. In our original system, $c^{*}$ is approximated as

$\displaystyle c^{*}(x_{m},t)\approx c^{*^{k}}_{m}$
$\displaystyle=\frac{c^{k}_{N,m}+c^{k}_{N+1,m}}{2}.$

This can be rewritten as

$\displaystyle c^{*^{k}}_{m}$
$\displaystyle=\underbrace{\frac{c^{k}_{1_{N,m}}+c^{k}_{1_{N+1,m}}}{2}}_{\textstyle=\alpha^{k}_{m}}-\frac{j^{k}_{m}}{D^{s}_{\text{eff}}(T^{k}_{m})}\underbrace{\frac{\gamma_{N+1}+\gamma_{N}}{2}}_{\textstyle=\beta}.$

Treating $\alpha^{k}_{m}$ and $\beta$ as constants, $c^{*}(x,t)$ is a function
of $j$ and $T$,
$c^{*k}_{m}=\alpha^{k}_{m}-\frac{j^{k}_{m}}{D^{s}(T^{k}_{m})}\beta$. In this
way, none of the channel equations depend explicitly on variable $\bm{c}$.
Therefore, $\bm{c}$ does not need to be computed in the Newton iterations. We
can proceed with the Newton iteration with
$U=(\bm{u},\bm{j},\bm{\eta},\bm{T}\bm{\phi},\bm{\psi})$ where we solve a
linear system of a $\mathcal{O}(M)\times\mathcal{O}(M)$ matrix. Table 3
compares the size of the Jacobian for each method. We stress here that the
decoupled method solves the same system as the naïve method. We obtain the
same solution to machine precision. If there is a nonlinear dependence in the
particle diffusion, the decoupling will still apply, although the tridiagonal
solves (4.4) for each particle will need to be updated at each Newton step.

In Section 4.1, we further improve the efficiency by reordering the entries of
Jacobian and computing only the nonzero entries of the reduced Jacobian.

$N_{p}=N_{n}$ | $M_{p}=M_{n}$ | $M_{s}=M_{a}=M_{z}$ |  $N_{\text{jac}}$: naive |  $N_{\text{jac}}$: decoupled
---|---|---|---|---
10 | 10 | 5 | 411 | 171
20 | 20 | 5 | 1171 | 291
30 | 30 | 5 | 2331 | 411
40 | 40 | 5 | 3891 | 531
50 | 50 | 5 | 5851 | 651
Table 3: The table shows the size of the Jacobian for each method at different
resolutions. The same number of grid points in $r$ and $x$ for positive and
negative electrode was used. The number of grid points was set fixed to $5$ in
the separator, positive and negative current collector. The size of the
Jacobian grows by $\mathcal{O}(NM)$ for the naive method, whereas only
$\mathcal{O}(M)$ for decoupled.

###  4.1 Reordering the Jacobian

The vector $U$ from our decoupled system can be reordered such that the
Jacobian is a banded diagonal matrix. By reordering, we can use an appropriate
banded linear solver to enhance the performance. We rearrange $U$ in an
alternating fashion,

$\begin{split}U&=[\begin{matrix}\bm{u}_{p}&\bm{u}_{s}&\bm{u}_{n}&\bm{j}_{p}&\bm{j}_{n}&\bm{\eta}_{p}&\bm{\eta}_{n}\end{matrix}\\\
&\qquad\qquad\begin{matrix}\bm{\Psi}_{p}&\bm{\Psi}_{n}&\bm{\Phi}_{p}&\bm{\Phi}_{s}&\bm{\Phi}_{n}&\bm{T}_{a}&\bm{T}_{p}&\bm{T}_{s}&\bm{T}_{n}&\bm{T}_{z}\end{matrix}]^{T}\end{split}$
(4.5) $\downarrow$
$\begin{split}U_{r}&=[\begin{matrix}\bm{T}_{a}&{u}_{p_{0}}&{\Psi}_{p_{0}}&{\Phi}_{p_{0}}&{T}_{p_{0}}\end{matrix}\\\
&\qquad\qquad\begin{matrix}u_{p_{1}}&j_{p_{0}}&\eta_{p_{0}}&\Phi_{s_{p_{1}}}&\Phi_{e_{p_{1}}}&T_{p_{1}}&u_{p_{2}}\cdots&T_{p_{Mp}}\end{matrix}\\\
&\qquad\qquad\qquad\begin{matrix}{u}_{p_{Mp+1}}{\Psi}_{p_{Mp+1}}&{\Phi}_{p_{Mp+1}}&{T}_{p_{Mp+1}}\end{matrix}\\\
&\qquad\qquad\qquad\qquad\begin{matrix}{u}_{s_{0}}&{\Phi}_{s_{0}}&{T}_{s_{0}}&u_{s_{1}}&\Phi_{e_{s_{1}}}&T_{s_{1}}&\cdots&{u}_{s_{Ms+1}}&{\Phi}_{s_{Ms+1}}&{T}_{s_{Ms+1}}\end{matrix}\\\
&\qquad\qquad\qquad\qquad\qquad\begin{matrix}\cdots&u_{n_{1}}&j_{n_{0}}&\eta_{n_{0}}&\Phi_{s_{n_{1}}}&\Phi_{e_{n_{1}}}&T_{n_{1}}&\cdots&\bm{T}_{z}\end{matrix}]^{T}\end{split}$
(4.6)

The variables are ordered such that $u_{i},j_{i},\eta_{i},\Psi_{i},\Phi_{i}$
and $T_{i}$ at each $i-$th grid point are next to each other, rather than the
original set up described in (3.2). Figure 3 illustrates how the sparsity
pattern of the Jacobian changes.

Figure 3: The sparsity pattern of Jacobian from the unordered system on the
left. On the right is the sparsity pattern after the Jacobian is reordered.

#### Banded matrix

The reordered Jacobian, which we denote by $\mathcal{J}_{r}$, can be stored in
a diagonal ordered form. There are $23$ nonzero lower and $23$ upper diagonals
in $\mathcal{J}_{r}$. Then, $\mathcal{J}_{r}$ can be represented as a matrix
of size $(23+23+1)\times N_{tot}$ where $N_{tot}$ is the size of $U$. We
denote this banded matrix by $\mathcal{J}_{b}$. The entries are stored in the
following way,

$\displaystyle\mathcal{J}_{b}[u+i-j,j]=\mathcal{J}_{r}[i,j],\quad
i=0,1,\cdots,(l+u),\quad j=0,1,\cdots,N_{tot}.$

The number of rows in $\mathcal{J}_{b}$ is fixed for a given model, as the
number of rows correspond to the number of partial derivatives in the system.
$\mathcal{J}_{b}$ is passed to the banded linear solver, which uses $LU$
decomposition to solve the linear system.

###  4.2 Constructing the Jacobian

We describe some of the performance gains that can be achieved by using JAX in
Python. JAX is especially useful for high-performance computations because it
uses a domain-specific compiler XLA (Accelerated Linear Algebra) that can
optimize computations at a lower level with improved memory usage and speed.
Compiling code just-in-time (jit) via XLA in Python accelerates the
performance because the source code is directly compiled to machine code as
the program runs. This is often faster than interpreting, which is typical for
languages like Python, making the performance comparable to that of compiled
languages like C. The readers can refer to [19] and [26] for more information.

In our problem, a fast compiler that optimizes differentiation computation is
very useful because we need to compute the system Jacobian at $U$ at every
Newton step. In addition to differentiation, the compiler can speed up any
array manipulations. Any functions that are compiled using XLA incur initial
overhead when first called, and is cached to be used later. This is useful for
functions that are called many times, as is the case in our problem with
multiple Newton updates and time stepping. For example, the first time we
build the Jacobian, there is an overhead of roughly 15 seconds (independent of
the size of the vector). The proceeding computations are fast.

However, there are limitations of using JAX’s built-in forward-mode Jacobian
function. While JAX can automatically output the Jacobian when $U$ is fed into
its built-in function, the output is dense. Currently, the sparse output of
Jacobians is not supported. This creates an extra overhead of converting the
matrix into a sparse form to use sparse linear solvers. This drawback also
causes a computational bottleneck. The evaluation of Jacobian with JAX’s
built-in function is quite costly for large vectors, even when compiled with
XLA. This is likely coming from the redundant computation effort computing
partial derivatives that evaluate to zeros. By the method discussed in
Sections 4.1, we bypass this problem by directly constructing the banded
reordered Jacobian, $\mathcal{J}_{b}$, without using JAX’s built-in Jacobian
function.

For each discretised PDE equation, we take the partial derivatives with
respect to their arguments via backward mode AD on a scalar level
($\mathbb{R}^{\text{Ntot}}\rightarrow\mathbb{R}$). This avoids computing
redundant zeros. Many of the discretised equations are repeated over the grid
points, but they can be vectorized. We predetermine the sparsity of
$\mathcal{J}_{b}$ prior to the simulation, which is unique to the reordering
scheme and the model. At each Newton update, we populate the matrix with
partial derivatives with the predetermined sparsity. The array update process
is sped up by compilation via XLA. This method provides an efficient way
constructing the Jacobian that takes advantage of the tools provided by JAX.

##  5 Experimental setup

All simulations were done on Mac CPU, with 3.5 GHz Dual-Core Intel Core i7
processor and 8GB RAM, in Python. Our code is freely available on Github at
https://github.com/hanrach/p2d_fast_solver. We ran the simulations until the
full battery discharge at $1$C, which is observed at $t=3620$ seconds. At the
current collectors and the separator, coarse grids ($M_{a},M_{s},M_{z}=5$)
were used for speed, because of the simpler governing equations in those
sections. We varied the grid sizes in the electrodes, with same number of grid
points at each positive and negative electrode in $x$ and $r$
($M_{p}=N_{p}=M_{n}=N_{n}$). The time step in the backward Euler method was
fixed to $\Delta t=10$.

###  5.1 Voltage curves under different C rates

We experimentally validate our solver by comparing the voltage curves against
Figure 9 in [15]. We run the simulation until the full discharge at different
C rates with the heat exchange coefficient $h$ fixed to $1$W/(m2K). We plot
how the average temperature and the battery voltage changes over time in
Figure 4. The battery voltage is defined as the difference in solid phase
potential at the two ends of the electrodes,
$V=\Psi_{p_{0}}-\Psi_{n_{M_{n}}}$. As expected, at higher C rates the
temperature rises rapidly and the voltage drops quickly to the cut-off voltage
of $2.5$V. These figures are qualitatively similar to Figure 9 in [15]. The
script `examples/different_C_rates.py` was used to generate Figure 4.

Figure 4: Temperature and voltage changes during the discharge cycle under
0.5C, 1C and 2C.

###  5.2 Comparsion between the naïve method and the fast solver

We compare the performance between the naïve method, which does not include
the decoupling and uses the built-in Jacobian function in JAX, and the fast
method where we decouple $c$ as well as construct the Jacobian in a diagonal
ordered form directly. When solving the linear system
$J\delta^{p}=\mathbf{U}^{p}$ for the naïve method, we convert the dense
Jacobian into a sparse matrix and apply a generic sparse linear solver. Table
4 summarizes the methods.

| Decoupled | Jacobian Function | Linear Solver
---|---|---|---
Naïve | No | JAX | Generic sparse
Fast | Yes | Custom | Tridiagonal and banded
Table 4: The fast method improves upon the naïve method by decoupling the
solid particle equations from the electrolyte equations, and employing the
custom Jacobian function and linear solvers that exploit the structure of the
matrices.

Figure 5 compares the cumulative time spent on solving the linear systems
within the Newton iterations. The overhead of conversion for the naïve method
is not costly compared to the actual time to solve the linear systems. This
factor is not included in the measurements in Figure 5. At a moderately fine
grid $M=30$, the naïve method spends a total of about $3.5$ seconds. The fast
solver is roughly $12$ times faster than the naïve method, spending about
$0.3$ seconds.

Figure 5: The total time of linear solves after the solver is run up to
$t=3520$. The $x$-axis corresponds to the size of the solution (particle and
channel). $M=N$ tested are $5,10,20$ and $30$. The time to assemble the
coefficients of the linear systems (4.2) are not included, since they can be
precomputed.

Figure 6 compares the total time of simulation between the naive solver and
the fast solver for varied $M=N$. We especially observe that the time spent on
evaluating the Jacobian is drastically reduced. At $M=30$, the fast solver
outperforms the naive method by $\approx 5000$ times.

Figure 6: The total simulation time until full discharge at $1$C. $M=N$ tested
are $5,10,20$ and $30$.

The code used to generate the performance measurements can be found in the
script `performance/full_simulation_comparison.py`. The data in Table 6 was
also generated by the same code. The script
`performance/reordered_performance.py` was used to obtain data for Table 7 and
Table 5.

###  5.3 Comparison between the fast solver and other solvers

We also compare the solver against LIONSIMBA [15], DUALFOIL [7] and COMSOL
[14]. We take the P2D model simulation data of LIONSIMBA, DUALFOIL and COMSOL
from [15], which simulates a discharge cycle of 1C under isothermal
conditions, done on a Windows 7 PC with 8GB of RAM. We also include the
runtimes of LIONSIMBA obtained on our machine (3.5 GHz Dual-Core Intel Core i7
with 8GB of RAM). For these simulations, the grid size was kept the same
across all battery sections. Table 5 shows that the simulation times from our
fast solver stay roughly constant for the resolutions tested, whereas this is
not true for other solvers.

$M$ | 10 | 20 | 30 | 40 | 50
---|---|---|---|---|---
Fast solver† | 4.09 s | 4.24 s | 3.98 s | 4.15 s | 4.49 s
LIONSIMBA† | 7.46 s | 9.85 s | 15.48 s | 26.80 s | 54.60 s
LIONSIMBA* | 28 s | 57 s | 105 s | 134 s | 223 s
DUALFOIL* | 28 s | 69 s | 97 s | 137 s | 185 s
COMSOL* | 96 s | 114 s | 143 s | 189 s | 244 s
Table 5: Comparison of simulation time between solvers. †The data was obtained
from 3.5 GHz Dual-Core Intel Core i7 with 8GB of RAM. *The data was taken from
[15] collected on a different machine.

##  6 Results

The fast method improves the efficiency of the naive solver in two ways.
Firstly, by decoupling the system, we reduce the time to solve the linear
systems, as shown in Figure 5. The figure compares performance between the
general sparse linear solver and the banded solver (scipy.linalg.solve_banded)
along with tridiagonal solves with $LU$ factorisation. By reformulating the
problem and taking advantage of the linear structure in (4.2), we improve the
linear algebraic efficiency.

Secondly, we work around the major bottleneck of the computation by
implementing our own Jacobian building function. Under the AD framework with
JAX, a large portion of the simulation time is spent on evaluating the
Jacobian and $\bm{F}(U)$. This dominates the total simulation time shown in
Figure 6. JAX’s Jacobian function takes up to six hours to evaluate the
Jacobians with the full channel and particle system. Decoupling alone without
using our custom Jacobian building function (the decoupled method)
significantly cuts the time to evaluate the Jacobian. Evaluating the Jacobian
of size $\mathcal{O}(M)\times\mathcal{O}(M)$ is much cheaper with JAX’s
Jacobian function, as opposed to $\mathcal{O}(NM)\times\mathcal{O}(NM)$. For
all grid points tested, the decoupled method is faster than the naive method,
drastically outperforming as $M$ increases. The decoupled formulation achieves
performance that is $\approx 870$ times faster than the naïve method at
$M=30$.

The fast method further improves the efficiency. JAX is only used to compute
the partial derivatives via backward mode differentiation. Compared to the
decoupled method, the custom function reduces the time to evaluate the
Jacobian up to $\approx 7$ times at $M=30$. This is evidently due to the fact
that the banded matrix in diagonal ordered form (of size $47\times N_{tot}$)
has a fixed number of rows and only the number of columns increases with
$M_{a},M_{p},M_{s},M_{n}$ and $M_{z}$. Table 6 compares the time to auto
differentiate and evaluate $\mathcal{J}(U)$ and $\bm{F}(U)$ and the total
simulation time of the three methods discussed.

| Evaluation (s) | Total (s)
---|---|---
Naive |  $2.22\text{\times}{10}^{4}$ s |  $2.31\text{\times}{10}^{4}$ s
Decoupled only (decoupled method) | 15.72 s | 21.63 s
Decoupled & Banded (fast method) | 0.3943 s | 3.898 s
Table 6: Evaluation time is the time spent on evaluating $\mathcal{J}(U)$ and
$\bm{F}(U)$. The simulations were compared at $M=30$.

Table 7 summarises the performance of the fast solver. There is an initial
overhead ($\approx 15$ seconds) that comes with using XLA compiler, which can
be precomputed before the actual simulation. Therefore, we do not factor it
into the total time. We observe that the total time spent on the simulation
stays roughly constant between three to five seconds. The grids at all battery
domains were set equal, with $N_{p}=N_{n}=10,20,30,40,50$ and
$M_{a}=M_{p}=M_{s}=M_{n}=M_{z}=10,20,30,40,50$. For these grid resolutions, we
observe that the evaluation and solve times are not the dominating factors but
the miscellaneous calculation overhead. This explains why the simulation time
seem constant with respect to $N$ and $M$. Choosing $M=N\gtrapprox 100$, a
clearer upward trend emerges where the solve time dominates the computation
(results not included). For practical applications, such fine resolution may
not be needed and $M=N\leq 50$ would suffice.

$M$ | 10 | 20 | 30 | 40 | 50
---|---|---|---|---|---
Solve | 0.213 | 0.289 | 0.346 | 0.4174 | 0.532
Evaluation | 0.376 | 0.412 | 0.386 | 0.426 | 0.483
Total | 4.092 | 4.243 | 3.987 | 4.1536 | 4.494
Table 7: Averaged times of the fast solver over 20 runs. Solve time indicates
the total time spent on linear solves. Evaluation time is the time spent on
evaluating $\mathcal{J}(U)$ and $\bm{F}(U)$.

##  7 Conclusion

We have implemented a fast and robust implementation of the full P2D model
without reducing the model. We give a thorough description of the
discretisation of the PDAE involving cell-centred finite differences and
backward Euler time stepping. The resulting system is a set of nonlinear
equations, which can be solved by Newton’s method.

In order to robustly compute the Jacobian at every Newton’s iteration,
automatic differentiation is used via the Python package JAX. Because the full
P2D model has a linear structure in the particle equations (2.1), we are able
to decouple the particle and channel equations at each time step. After
decoupling, we solve one large tridiagonal linear system per time step and
smaller sparse systems (channel equations) within the Newton iterations. The
full system is thus solved more efficiently.

The decoupled method enables us to run the simulation 1C discharge cycle in a
reasonable time under the AD framework, whereas such framework is infeasible
when we use the naive method. This was made possible by rearranging the
variables in the channel equations to create a matrix that is banded, which
can be solved efficiently. The computation was made more efficient by directly
building the banded Jacobian using the XLA compiler.

An additional way we can increase the efficiency is implementing a adaptive
time stepping scheme, based on error on $\bm{c}$ only. We have implemented
this for the scaled P2D model based on [27], which is a work in progress by
Hennessey et al. The adaptive time stepping takes larger time steps in the
beginning and gradually adjusts to smaller steps as the problem gets harder
near the end of the cycle. For the scaled P2D model, it improved the
efficiency by roughly a factor of two compared to the fixed time stepping
scheme. This is future work to be explored.

We have demonstrated that by exploiting the structure of the problem, and
utilizing state-of-art software, the complex P2D model can be solved
accurately in under a few seconds. The fast solver can be potentially used as
part of larger computations, such as optimizing multiple-battery pack
behaviours. The solver would be especially useful where the P2D model has to
be solved repeatedly over many cycles for estimation, optimization or control
purposes.

##  Appendix A Symbols and Additional equations

[title=List of Symbols,type=symbols,style=long]

Open Circuit Potential
---
$\bm{U}_{p}=\bm{U}_{\text{p,ref}}+(T(x,t)-T_{\text{ref}})\frac{\partial\bm{U}_{p}}{\partial\bm{T}}\rvert_{T_{\text{ref}}}$
$\bm{U}_{n}=\bm{U}_{\text{n,ref}}+(T(x,t)-T_{\text{ref}})\frac{\partial\bm{U}_{n}}{\partial\bm{T}}\rvert_{T_{\text{ref}}}$
Entropy Change
$\frac{\partial\bm{U}_{p}}{\partial\bm{T}}\rvert_{T_{\text{ref}}}=-0.001\left(\frac{0.199521039-0.928373822\theta_{p}+1.364550689000003\theta_{p}^{2}-0.6115448939999998\theta_{p}^{3}}{1-5.661479886999997\theta_{p}+11.47636191\theta_{p}^{2}-9.82431213599998\theta_{p}^{3}+\
3.046755063\theta_{p}^{4}}\right)$
$\frac{\partial\bm{U}_{n}}{\partial\bm{T}}\rvert_{T_{\text{ref}}}=0.001\left(\frac{\begin{multlined}0.005269056+3.299265709\theta_{n}-91.79325798\theta_{n}^{2}+1004.911008\theta_{n}^{3}-5812.278127\theta_{n}^{4}\\\
+19329.7549\theta_{n}^{5}-37147.8947\theta_{n}^{6}+38379.18127\theta_{n}^{7}-16515.05308\theta_{n}^{8}\end{multlined}0.005269056+3.299265709\theta_{n}-91.79325798\theta_{n}^{2}+1004.911008\theta_{n}^{3}-5812.278127\theta_{n}^{4}\\\
+19329.7549\theta_{n}^{5}-37147.8947\theta_{n}^{6}+38379.18127\theta_{n}^{7}-16515.05308\theta_{n}^{8}}{\begin{multlined}1-48.09287227\theta_{n}+1017.234804\theta_{n}^{2}-10481.80419\theta_{n}^{3}+59431.3\theta_{n}^{4}\\\
-195881.6488\theta_{n}^{5}+374577.3152\theta_{n}^{6}-385821.1607\theta_{n}^{7}+165705.8597\theta_{n}^{8}\end{multlined}1-48.09287227\theta_{n}+1017.234804\theta_{n}^{2}-10481.80419\theta_{n}^{3}+59431.3\theta_{n}^{4}\\\
-195881.6488\theta_{n}^{5}+374577.3152\theta_{n}^{6}-385821.1607\theta_{n}^{7}+165705.8597\theta_{n}^{8}}\right)$
$\theta_{p}=\frac{c^{*}_{s,p}(x,t)}{c^{\text{max}}_{s,p}}$
$\theta_{n}=\frac{c^{*}_{s,n}(x,t)}{c^{\text{max}}_{n,p}}$
Open circuit potential
$U_{\text{p,ref}}=\frac{-4.656+88.669\theta_{p}^{2}-401.119\theta_{p}^{4}+342.909\theta_{p}^{6}-462.471\theta_{p}^{8}+433.434\theta_{p}^{1}0}{-1+18.933\theta_{p}^{2}-79.532\theta_{p}^{4}+37.311\theta_{p}^{6}-73.083\theta_{p}^{8}+95.96\theta_{p}^{1}0}$
$U_{\text{n,ref}}=0.7222+0.1387\theta_{n}+0.029\theta_{n}^{0.5}-\frac{0.0172}{\theta_{n}}+\frac{0.0019}{\theta^{1.5}}+0.2808e^{0.9-15\theta_{n}}-0.7984e^{0.4465\theta_{n}-0.4108}$
Heat source terms (electrodes)
$\bm{Q}_{\text{ohm}}=\sigma_{\text{eff},i}\left(\frac{\partial\Psi(x,t)}{\partial
x}\right)^{2}+\kappa_{\text{eff},i}\left(\frac{\partial\Phi_{e}(x,t)}{\partial
x}\right)^{2}+\frac{2\kappa_{\text{eff},i}RT(x,t)}{F}(1-t_{+})\frac{\partial\ln
u(x,t)}{\partial x}\frac{\Phi_{e}(x,t)}{\partial x}$
$\bm{Q}_{\text{rxn}}=Fa_{i}j(x,t)\eta_{i}(x,t)$
$\bm{Q}_{\text{rev}}=Fa_{i}j(x,t)T(x,t)\frac{\partial U_{i}}{\partial
T}\rvert_{T_{\text{ref}}}$
Heat source terms (separator)
$\bm{Q}_{\text{ohm}}=\kappa_{\text{eff},i}\left(\frac{\partial\Phi_{e}(x,t)}{\partial
x}\right)^{2}+\frac{2\kappa_{\text{eff},i}RT(x,t)}{F}(1-t_{+})\frac{\partial\ln
u(x,t)}{\partial x}\frac{\Phi_{e}(x,t)}{\partial x}$
Coefficients
$\bm{D}_{\text{eff},i}=\epsilon_{i}^{\text{brugg}_{i}}\times 10^{-4}\times
10^{-4.43-\frac{54}{T(x,t)-229-5\times 10^{-3}u(x,t)}-0.22\times
10^{-3}u(x,t)}$
$\begin{aligned} \kappa_{\text{eff},i}=\epsilon_{i}^{\text{brugg}_{i}}\times
10^{-4}\times u(x,t)\left(-10.5+0.668\times 10^{-3}u(x,t)+0.494\times
10^{-6}u(x,t)^{2}+\right.\\\ \left.T(x,t)(0.074-1.78\times
10^{-5}u(x,t)-8.86\times 10^{-10}u(x,t)^{2})+\right.\\\
\left.T(x,t)^{2}(-6.96\times 10^{-5}+2.8\times
10^{-8}u(x,t))\right)^{2}\end{aligned}$
$k_{\text{eff}}=k_{i}e^{\frac{E^{k_{i}}_{a}}{R}\left(\frac{1}{T(x,t)}-\frac{1}{T_{\text{ref}}}\right)}$
$D^{s}_{\text{eff},i}=D^{s}_{i}e^{\frac{E^{D^{s}_{i}}_{a}}{R}\left(\frac{1}{T(x,t)}-\frac{1}{T_{\text{ref}}}\right)}$
$\sigma_{\text{eff},i}=\sigma_{i}(1-\epsilon_{i}-\epsilon_{f,i})$
$\Upsilon=\frac{2(1-t_{+})R}{F}$

##  Appendix B Parameters used in simulation

|  | Al CC | Cathode | Separator | Anode | Carbon CC
---|---|---|---|---|---|---
$c_{e}^{\text{init}}$ |  [mol/m3] | - | 1000 | 1000 | 1000 | -
$c_{s}^{\text{avg, init}}$ |  [mol/m3] | - | 25751 | - | 26128 | -
$c_{s}^{\text{max}}$ |  [mol/m3] | - | 51554 | - | 30555 | -
$D_{i}$ |  [m${}^{2}/s$] | - | $7.5\times 10^{-10}$ | $7.5\times 10^{-10}$ | $7.5\times 10^{-10}$ | -
$D^{s}_{i}$ |  [m${}^{2}/s$] | - | $10^{-14}$ | - | $3.9\times 10^{-14}$ | -
$k_{i}$ |  [m${}^{2.5}/$ (mol0.5s)] | - | $2.334\times 10^{-11}$ | - | $5.301\times 10^{-11}$ | -
$l_{i}$ | [m] | $10^{-5}$ | $8\times 10^{-5}$ | $2.5\times 10^{-5}$ | $8.8\times 10^{-5}$ | $10^{-5}$
$R_{p,i}$ | [m] | - | $2\times 10^{-6}$ | - | $2\times 10^{-6}$ | -
$\rho_{i}$ |  [kg/m3] | 2700 | 2500 | 1100 | 2500 | 8940
$C_{p,i}$ | [J/(kg K)] | 897 | 700 | 700 | 700 | 385
$\lambda_{i}$ | [W/(m K)] | 237 | 2.1 | 0.16 | 1.7 | 401
$\sigma_{i}$ | [S/m] | $3.55\times 10^{7}$ | 100 | - | 100 | $5.96\times 10^{7}$
$\epsilon_{i}$ | - | - | 0.385 | 0.724 | 0.485 | -
$a_{i}$ |  [m2/m3] | - | 885000 | - | 723600 | -
$E^{D^{s}_{i}}_{a}$ | [J/mol] | - | 5000 | - | 5000 | -
$E^{k_{i}}_{a}$ | [J/mol] | - | 5000 | - | 5000 | -
brugg | - | - | 4 | 4 | 4 | -
$F$ | 96485 [C/mol] | - | - | - | - | -
$R$ | 8.314472 [J/(mol K)] | - | - | - | - | -
$t_{+}$ | 0.354 |  |  |  |  |
$\epsilon_{f,i}$ | - | - | 0.025 | - | 0.0326 | -

## References

  * [1] N. Nitta, F. Wu, J. T. Lee, G. Yushin, Li-ion battery materials: present and future, Materials Today 18 (5) (2015) 252 – 264\.  doi:https://doi.org/10.1016/j.mattod.2014.10.040.
URL http://www.sciencedirect.com/science/article/pii/S1369702114004118

  * [2] J. Heelan, E. Gratz, Z. Zheng, Q. Wang, M. Chen, D. Apelian, Y. Wang, Current and prospective li-ion battery recycling and recovery processes, The Journal of The Minerals, Metals & Materials Society (TMS) 68 (10) (2016) 2632–2638.
  * [3] B. Pattipati, K. Pattipati, J. P. Christopherson, S. M. Namburu, D. V. Prokhorov, Liu Qiao, Automotive battery management systems, in: 2008 IEEE AUTOTESTCON, 2008, pp. 581–586.
  * [4] K. Liu, K. Li, J. Deng, A novel hybrid data-driven method for li-ion battery internal temperature estimation, in: 2016 UKACC 11th International Conference on Control (CONTROL), IEEE, 2016, pp. 1–6.
  * [5] K. A. Severson, P. M. Attia, N. Jin, N. Perkins, B. Jiang, Z. Yang, M. H. Chen, M. Aykol, P. K. Herring, D. Fraggedakis, M. Z. Bazant, S. J. Harris, W. C. Chueh, R. D. Braatz, Data-driven prediction of battery cycle life before capacity degradation, Nature Energy 4 (5) (2019) 383–391.  doi:10.1038/s41560-019-0356-8.
URL https://doi.org/10.1038/s41560-019-0356-8

  * [6] M.-F. Ng, J. Zhao, Q. Yan, G. J. Conduit, Z. W. Seh, Predicting the state of charge and health of batteries using data-driven machine learning, Nature Machine Intelligence 2 (3) (2020) 161–170.  doi:10.1038/s42256-020-0156-7.
URL https://doi.org/10.1038/s42256-020-0156-7

  * [7] M. Doyle, Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell, Journal of The Electrochemical Society 140 (6) (1993) 1526.  doi:10.1149/1.2221597.
  * [8] M. Guo, G. Sikha, R. E. White, Single-particle model for a lithium-ion cell: Thermal behavior, Journal of The Electrochemical Society 158 (2) (2010) A122.
  * [9] E. Prada, D. Di Domenico, Y. Creff, J. Bernard, V. Sauvant-Moynot, F. Huet, Simplified electrochemical and thermal model of LiFePO4-graphite Li-ion batteries for fast charge applications, Journal of The Electrochemical Society 159 (9) (2012) A1508.
  * [10] V. Ramadesigan, V. Boovaragavan, J. C. Pirkle Jr, V. R. Subramanian, Efficient reformulation of solid-phase diffusion in physics-based lithium-ion battery models, Journal of The Electrochemical Society 157 (7) (2010) A854.
  * [11] X. Han, M. Ouyang, L. Lu, J. Li, Simplification of physics-based electrochemical model for lithium ion battery on electric vehicle. Part I: Diffusion simplification and single particle model, Journal of Power Sources 278 (2015) 802–813.
  * [12] A. Jokar, B. Rajabloo, M. Désilets, M. Lacroix, Review of simplified pseudo-two-dimensional models of lithium-ion batteries, Journal of Power Sources 327 (2016) 44–55.
  * [13] V. Ramadesigan, P. W. C. Northrop, S. De, S. Santhanagopalan, R. D. Braatz, V. R. Subramanian, Modeling and simulation of lithium-ion batteries from a systems engineering perspective, Journal of The Electrochemical Society 159 (3) (2012) R31–R45.  doi:10.1149/2.018203jes.
URL https://doi.org/10.1149%2F2.018203jes

  * [14] L. Cai, R. E. White, Mathematical modeling of a lithium ion battery with thermal effects in COMSOL Inc. Multiphysics (MP) software, Journal of Power Sources 196 (14) (2011) 5985 – 5989.  doi:https://doi.org/10.1016/j.jpowsour.2011.03.017.
URL http://www.sciencedirect.com/science/article/pii/S0378775311005994

  * [15] M. Torchio, L. Magni, R. B. Gopaluni, R. D. Braatz, D. M. Raimondo, LIONSIMBA: a Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control, Journal of The Electrochemical Society 163 (7) (2016) A1192.
  * [16] G. G. Botte, V. R. Subramanian, R. E. White, Mathematical modeling of secondary lithium batteries, Electrochimica Acta 45 (15-16) (2000) 2595–2609.
  * [17] L. R. Petzold, Description of DASSL: a differential/algebraic system solver, Tech. rep., Sandia National Labs., Livermore, CA (USA) (1982).
  * [18] A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, C. S. Woodward, Sundials: Suite of nonlinear and differential/algebraic equation solvers, ACM Transactions on Mathematical Software (TOMS) 31 (3) (2005) 363–396.
  * [19] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, S. Wanderman-Milne, JAX: composable transformations of Python+NumPy programs (2018).
URL http://github.com/google/jax

  * [20] V. D. Bruggeman, Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen. I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen, Annalen der physik 416 (7) (1935) 636–664.
  * [21] J. Newman, K. Thomas-Alyea, Electrochemical Systems, The ECS Series of Texts and Monographs, Wiley, 2004.
URL https://books.google.ca/books?id=vArZu0HM-xYC

  * [22] U. Ascher, L. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, Other Titles in Applied Mathematics, Society for Industrial and Applied Mathematics, 1998.
URL https://books.google.ca/books?id=VL51G5JYYAYC

  * [23] J. W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, Vol. 22 of Texts in Applied Mathematics, Springer-Verlag, New York, 1995\.
  * [24] R. L. Burden, J. D. Faires, Numerical analysis (2011).
  * [25] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, J. M. Siskind, Automatic differentiation in machine learning: a survey, The Journal of Machine Learning Research 18 (1) (2017) 5595–5637.
  * [26] X. team within Google, XLA: Optimizing compiler for machine learning.
URL https://www.tensorflow.org/xla

  * [27] M. Hennesy, I. Moyles, T. Myers, B. Wetton, The P2D model for Li-ion battery electrodes: When is it really needed?


††thanks: Corresponding author

# Parity-deformed $su(2)$ and $so(3)$ Algebras: a Basis for Quantum Optics and
Quantum Communications Applications

W. S. Chung<EMAIL_ADDRESS>Department of Physics and Research Institute of
Natural Science, College of Natural Science, Gyeongsang National University,
Jinju 660-701, Korea  H. Hassanabadi<EMAIL_ADDRESS>Faculty of
Physics, Shahrood University of Technology, Shahrood, Iran  L. M. Nieto
<EMAIL_ADDRESS>https://ftao.uva.es/lmnieto/ Departamento de
Física Teórica, Atómica y Óptica, and Laboratory for Disruptive
Interdisciplinary Science (LaDIS), Universidad de Valladolid, 47011
Valladolid, Spain  S. Zarrinkamar<EMAIL_ADDRESS>Departamento de
Física Teórica, Atómica y Óptica, and Laboratory for Disruptive
Interdisciplinary Science (LaDIS), Universidad de Valladolid, 47011
Valladolid, Spain  Department of Basic Sciences, Garmsar Branch, Islamic Azad
University, Garmsar, Iran

###### Abstract

Having in mind the significance of parity (reflection) in various areas of
physics, the single-mode and two-mode Wigner algebras are considered adding to
them a reflection operator. The associated deformed $su(2)$ algebra,
$su_{\nu}(2)$, and the deformed $so(3)$ algebra, $so_{\nu}(3)$, are
constructed for the widely used Jordan-Schwinger and Holstein-Primakoff
realizations, commenting on various aspects and ingredients of the formalism
for both single-mode and two-mode cases. Finally, in this framework the
parity-deformed $so_{\nu}(3)$ representation is analyzed, due to its potential
application in the study of qubit and qutrit systems.

Keywords: Deformed $su(2)$ algebra, deformed $so(3)$ algebra, parity, Jordan-
Schwinger realization, Holstein-Primakoff realization, quantum communications.

##  1 Introduction

After a century of the advent of quantum mechanics, this theory still has
outstanding challenges and open questions to debate in many basic concepts and
formulas of the theory. As a result, a variety of modifications to the related
fundamental relationships have been proposed, primarily in an attempt to
investigate not only problems in physics, but also in other domains of science
and technology beyond our current experimental techniques and facilities. On
the other hand, it is well known that theoretical physics has always been
inspiring and that the history of physics is full of examples of theoretical
predictions prior to their experimental verification. This makes theoretical
and mathematical physics continue to be an attractive field of study for the
scientific community, competing with the incredible technological achievements
that are being achieved in recent times.

Given the current great technological importance of quantum theory in modern
technologies, in what is called second quantum revolution, which covers
applications in various fields, such as quantum information, quantum
communication, quantum teleportation and quantum technologies in general, we
are interested in analyzing possible generalizations of some fundamental basic
results that can be contrasted using current or future experimental techniques
in development. In particular, remembering here the importance that the
analysis of the harmonic oscillator problem has had in both optics and quantum
technologies, we are going to theoretically consider a generalization of
Wigner’s algebra modified by parity inclusion. To be precise, what is
considered here is a kind of parity-warped operator algebra, sometimes called
Wigner, Dunkl, or Wigner-Dunkl algebra. The pillars of the formulation of the
problem and the consideration of the reflection operator go back to the works
of Holstein and Primakoff, Wigner, Baxter and later Dunkl [1, 2, 3, 4, 5].
Generalization of algebra is usually done by the Jordan-Schwinger and
Holstein-Primakoff realizations of the associated creation and annihilation
operators and their extensions [6, 7, 8, 9, 10]. In physics, however, the
applications are quite interesting and include a variety of concepts such as
coherent states, compressed states, cat states, fidelity, entanglement
transfer, etc. Let us now review, in chronological order, some of the most
notable works that are closely related to the objectives of our work.

In pioneering work, with a goal quite different from that pursued here,
Plyushchay investigated the deformed Heisenberg algebra with reflection in
[11], demonstrating that the normalized deformed algebra has the form of a
gluon algebra. Using a reflection-deformed Heisenberg algebra, a unified
description of anyons, bosons and fermions was analyzed in the instructive
survey [12]. The hidden supersymmetry and related tri-supersymmetric structure
of some physical examples were investigated in full detail by Jakubský et al.
[13], who also studied the role of parity. In interesting novel work, Dehghani
and others studied the Wigner-Heisenberg algebra in relation to so-called cat
states [14]. The same group later analyzed a parity-deformed version of the
Jaynes-Cummings model, studying in detail the entangled states in a lossy
cavity interacting with an external field [15]. The hidden superconformal
symmetry of the harmonic oscillator problem was obtained by a non-local Foldy-
Wouthuysen transformation on [16]. The consequences of parity deformation on
entanglement transfer to distant atomic qubits were considered in [17], where
the qubits were assumed to be embedded in two lossy cavities connected by a
leaky fiber.

In [18] a form of generalized Dunkl oscillator has been formulated in both
relativistic and non-relativistic cases. Quantum error tracking with parity
measurements was discussed in [19]. The analysis of dynamical transformations
between quantum reference frames, including the parity properties, was
formulated by Ballesteros et al. [20]. The optical properties of a pumped-
dissipative quantum dot cavity with a parity-deformed Heisenberg algebra were
investigated in [21] and it was shown that the deformation produces a
collective effect on the photonic emission of the system. Fakhri et al.
studied a two-photon parity-deformed Jaynes-Cummings model and discussed in
detail the effect of deformation on important physical quantities, including
fidelity [22]. Fakhri and Sayyah-Fard [23] derived and discussed in depth the
$x$ representations of the coherent and cat states of the para-bosonic
oscillator algebra. Dong and his collaborators [24] studied in depth the
supersymmetric structure of the Wigner-Dunkl problem. The effect of parity on
quantum optimization and computation was well covered in [25], [26] and [27].
The problem of two qubits moving in a parity-deformed field was considered in
[28], where it was shown that parity warping can give us greater protection
against entanglement. The gauge invariance of a charged Wigner-Dunkl system
was analyzed in [29] and it was shown that the gauge invariance holds in case
of minimal coupling to a vector potential. In reference [30] a generalization
of the Wigner-Dunkl problem was considered and the related coherent states
were studied, as well as the corresponding Mandel parameter and the
bunching/anti-bunching effect. Very instructive and quite up-to-date
applications, as well as the importance of parity in the study of fidelity and
quantum error correction, are investigated in references [31, 32].

In [33] the Jaynes-Cummings model with intensity-dependent coupling with a
Holstein-Primakoff $SU(1,1)$ state was considered, and it was shown that the
revivals of the radiation squeezing exhibit a periodic character. The
generalization of the Jaynes-Cummings model was studied using the deformed
harmonic oscillator problem and the Holstein-Primakoff realization of
$su(1,1)$ by Chaichian and coworkers [34]. The Holstein–Primakoff
transformation is also present in many other fields, including the laser-
cooled trapped ion interaction [35], the $N$-body Lipkin model [36], and
equilibrium electron-magnon systems [37]. The transfer of quantum information
between a nanomechanical resonator and a superconducting transmission line
resonator via a charged qubit in the case of two non-degenerate modes was
considered using the Jordan-Schwinger realization in [38]. The generalization
of the Jaynes-Cummings and Dicke models was formulated using a generalized
Jordan–Shchwinger realization and therefore the associated integrability and
superintegrability were discussed in [39]. The realization of three-
dimensional Jordan-Schwinger problems was studied in [40]. Finally let us
mention that the $SO(3)$ group has been used to investigate modern
applications of quantum theory, including the maximally entangled states of
qubits [41, 42] and qutrits [43]. The perfect transfer of path-entangled and
single-photon states in the so-called ${J_{x}}$ photonic lattices was
discussed in the novel paper by Perej-Leija and co-workers [44].

In the present work, due to the wide application of parity-dependent systems
and the already mentioned applications in quantum communication and computing,
we aim to investigate the effect of parity in a special form of generalized
Wigner algebra.

This article is organized as follows. In Section 2, the single-mode Wigner
algebra and its realizations are reviewed, while the two-mode Wigner algebra
is discussed in Section 3. In both cases, what will be called $\nu$-deformed
Wigner algebra will be defined and studied. Two realizations of the algebra
$su_{\nu}(2)$, the Jordan-Schwinger realization and the Holstein-Primakoff
realization, are investigated in Section 4, while a representation of the
$so_{\nu}(3)$ algebra is studied in Section 5. Finally, in Section 6 we
summarize the main conclusions of the present work.

##  2 Single-mode Wigner Algebra and its realizations

The single mode Wigner algebra [1] is given by the commutation relations

$[a,a^{\dagger}]=1+2\nu R,\qquad[N,a^{\dagger}]=a^{\dagger},\qquad[N,a]=-a,$
(1)

where $a$ and $a^{\dagger}$ are annihilation and creation operators,
respectively, and $2\nu$ denotes a real positive constant called Wigner
parameter. The reflection operator $R$ is defined by

$R=(-1)^{N},$ (2)

and obeys the following anticommutation relations

$\\{R,a\\}=\\{a^{\dagger},R\\}=0,$ (3)

as well as

$R^{\dagger}=R^{-1}=R,\qquad R^{2}=I.$ (4)

Let us introduce the eigenvectors of the number operator $N$ as
$N|n\rangle=n|n\rangle$, $n=0,1,2,\dots$ [18], where $N$ is a Hermitian
operator and the ground state $|0\rangle$ satisfies $a|0\rangle=0.$ From (1),
we have the following representation

$\displaystyle a|n\rangle=\sqrt{[n]_{\nu}}\ |n-1\rangle,\qquad
a^{\dagger}|n\rangle=\sqrt{[n+1]_{\nu}}\ |n+1\rangle,$ (5)

where the $\nu$-deformed numbers are defined as

$[n]_{\nu}=n+\nu(1-(-1)^{n}).$ (6)

In fact, the first $\nu$-deformed numbers are [18]

$[0]_{\nu}=0,\qquad[1]_{\nu}=1+2\nu,\qquad[2]_{\nu}=2,\qquad[3]_{\nu}=3+2\nu,\qquad[4]_{\nu}=4.$
(7)

In general, we have

$[2k]_{\nu}=2k,\qquad[2k+1]_{\nu}=2k+1+2\nu,\qquad k=0,1,2,\dots$ (8)

Some properties for the $\nu$-deformed numbers are the following

$\displaystyle[n]_{\nu}+[n+1]_{\nu}$ $\displaystyle=2n+1+2\nu,$
$\displaystyle[n+2]_{\nu}-[n]_{\nu}$ $\displaystyle=2,$
$\displaystyle[n]_{\nu}+[n+2k+1]_{\nu}$ $\displaystyle=2n+2k+1+2\nu,\quad
k\in\mathbb{Z},$ $\displaystyle[n+2k]_{\nu}-[n]_{\nu}$ $\displaystyle=2,\quad
k\in\mathbb{Z},$ $\displaystyle[m]_{\nu}[n+1]_{\nu}-[n]_{\nu}[m+1]_{\nu}$
$\displaystyle=m-n-\nu(2n+1)(-1)^{m}$
$\displaystyle\quad+\nu(2m+1)(-1)^{n}-2\nu^{2}((-1)^{m}-(-1)^{n}).$ (9)

If we spilt the set of integers as
$\mathbb{Z}=\mathbb{Z}_{0}\oplus\mathbb{Z}_{1}$, with
$\mathbb{Z}_{0}=\\{0,2,4,\dots\\}$ and $\mathbb{Z}_{1}=\\{1,3,5,\dots\\}$,
Equation (9) can be rewritten as

$[m]_{\nu}[n+1]_{\nu}-[n]_{\nu}[m+1]_{\nu}=\left\\{\begin{array}[]{ll}(m-n)(1+2\nu),&n\in\mathbb{Z}_{0},m\in\mathbb{Z}_{0},\\\
m-n+2\nu(m+n+1)+4\nu^{2},&n\in\mathbb{Z}_{0},m\in\mathbb{Z}_{1},\\\
m-n-2\nu(m+n+1)-4\nu^{2},&n\in\mathbb{Z}_{1},m\in\mathbb{Z}_{0},\\\
(m-n)(1-2\nu),&n\in\mathbb{Z}_{1},m\in\mathbb{Z}_{1}.\end{array}\right.$ (10)

Now we will look for the coordinate realization for the single-mode Wigner
algebra (1). The realization depends on the choice of the basis functions. For
the ordinary boson algebra, we have the realization [18]

$a=\partial_{x},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ a^{\dagger}=x,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ N=a^{\dagger}a=x\partial_{x}.$ (11)

For the monomial basis given by $f_{n}(x)=x^{n}$ ($n=0,1,2,\dots$), this
realization is

$af_{n}(x)=nf_{n-1}(x),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ a^{\dagger}f_{n}(x)=f_{n+1}(x),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ Nf_{n}(x)=nf_{n}(x).$ (12)

Since $af_{0}=0$, this representation is bounded from below and $f_{0}$ is
called the ground state. Considering now the so-called quasi-polynomial basis
$\phi_{n}(x)$ given by

$\phi_{n}(x)=\frac{\Gamma(x+1)}{\Gamma(x+1-n)},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ n=0,1,2,\dots,$ (13)

a new realization for the ordinary boson algebra (1) appears as

$a=e^{\partial_{x}}-1,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ a^{\dagger}=xe^{-\partial_{x}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
N=x\left(1-e^{-\partial_{x}}\right)=a^{\dagger}a,$ (14)

whose corresponding representation is

$a\phi_{n}(x)=n\phi_{n-1}(x),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ a^{\dagger}\phi_{n}(x)=\phi_{n+1}(x),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ N\phi_{n}(x)=n\phi_{n}(x).$ (15)

###  2.1 Monomial basis Wigner algebra realization

For the monomial basis $f_{n}(x)=x^{n}$, the realization of the Wigner algebra
(1) is

$a=\partial_{x}+\frac{1}{x}\left(\nu-\nu(-1)^{x\partial_{x}}\right),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
a^{\dagger}=x,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
N=x\partial_{x},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
R=(-1)^{x\partial_{x}}.$ (16)

The representation of this algebra in this basis is

$af_{n}(x)=(n+\nu-\nu(-1)^{n})f_{n-1}(x),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
a^{\dagger}f_{n}(x)=f_{n+1}(x),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ Nf_{n}(x)=nf_{n}(x).$ (17)

###  2.2  Quasi-polynomial basis Wigner algebra realization

In this subsection we derive the quasi-polynomial basis for the Wigner
algebra, which is related to $Z_{2}$-grading, which implies that
$Rf_{n}(x)=(-1)^{n}f_{n}(x)$. Thus, we change the quasi polynomial for an
ordinary boson algebra into the $Z_{2}$-graded quasi polynomial for the Wigner
algebra. The $Z_{2}$-graded quasi polynomial basis is defined as

$\phi^{\nu}_{n}(x):=\prod_{k=0}^{n-1}(x-k-\nu(-1)^{k}),\quad n\geq 1,$ (18)

with $\phi^{\nu}_{0}(x):=1$. Therefore, the realization is

$a=T_{\nu}e^{\partial_{x}}-1,\qquad
a^{\dagger}=(x-\nu)T_{\nu}e^{-\partial_{x}},$ (19)

where $\partial_{x}:=\frac{d}{dx}$ is the ordinary derivative,
$e^{\pm\partial_{x}}f(x):=f(x\pm 1)$ and

$T_{\nu}F(\nu):=(-1)^{\nu\frac{d}{d\nu}}F(\nu)=F(-\nu).$

One can easily check that the above realization obeys

$\displaystyle a\phi^{\nu}_{n}(x)$ $\displaystyle=[n]_{\nu}\
\phi^{\nu}_{n-1}(x),$ $\displaystyle a^{\dagger}\phi^{\nu}_{n}(x)$
$\displaystyle=\phi^{\nu}_{n+1}(x),$ $\displaystyle\phi^{\nu}_{n}(x)$
$\displaystyle=(a^{\dagger})^{n}1,$ (20)
$\displaystyle[a,a^{\dagger}]\phi^{\nu}_{n}(x)$ $\displaystyle=(1+2\nu
R)\,\phi^{\nu}_{n}(x),$

where

$N=a^{\dagger}a-\nu+\nu R.$ (21)

Therefore, the set of polynomials $\\{{\phi}^{\nu}_{n}(x)|n=0,1,\dots\\}$
provides a basis for the realization of the Wigner algebra.

##  3 Two-mode Wigner algebra

We can extend the single-mode Wigner oscillator algebra into the two-mode
Wigner oscillator algebra as follows [4]

$\displaystyle[a_{i},a^{\dagger}_{j}]$ $\displaystyle=\delta_{ij}(1+2\nu
R_{i}),\qquad[a_{i},a_{j}]=0,\qquad[a^{\dagger}_{i},a^{\dagger}_{j}]=0,$ (22)
$\displaystyle[N_{i},a^{\dagger}_{j}]$
$\displaystyle=\delta_{ij}a^{\dagger}_{i},\qquad\qquad\quad\ \
[N_{i},a_{j}]=-\delta_{ij}a_{i},$ (23)

where $i=1,2$, and

$R_{i}=(-1)^{N_{i}},\qquad\\{R_{i},a_{i}\\}=\\{a^{\dagger}_{i},R_{i}\\}=0,\qquad[R_{i},a_{j}]=[R_{i},a^{\dagger}_{j}]=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ i\neq
j.$ (24)

The number operators and the associated ladder operators satisfy (6)

$a^{\dagger}_{i}a_{i}=[N_{i}]_{\nu}.$ (25)

The eigenvectors of the number operators $N_{i}$ are defined as

$N_{i}|n_{1},n_{2}\rangle=n_{i}|n_{1},n_{2}\rangle,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ n_{i}=0,1,2,\dots,$ (26)

where $N_{i}$ are Hermitian, and the ground state $|0,0\rangle$ obeys
$a_{i}|0,0\rangle=0$, $i=1,2$. From (23), we have

$\displaystyle a_{1}|n_{1},n_{2}\rangle$ $\displaystyle=\sqrt{[n_{1}]_{\nu}}\
|n_{1}-1,n_{2}\rangle,\quad$ $\displaystyle
a^{\dagger}_{1}|n_{1},n_{2}\rangle=\sqrt{[n_{1}+1]_{\nu}}\
|n_{1}+1,n_{2}\rangle,$ (27) $\displaystyle a_{2}|n_{1},n_{2}\rangle$
$\displaystyle=\sqrt{[n_{2}]_{\nu}}\ |n_{1},n_{2}-1\rangle,\quad$
$\displaystyle a^{\dagger}_{2}|n_{1},n_{2}\rangle=\sqrt{[n_{2}+1]_{\nu}}\
|n_{1},n_{2}+1\rangle.$ (28)

##  4 Realizations of $su(2)$ and $su_{\nu}(2)$ Wigner algebras

We are going to consider two different realizations of the Wigner algebras
$su(2)$ and $su_{\nu}(2)$: the Jordan-Schwinger realization and the Holstein-
Primakoff realization.

###  4.1 Jordan-Schwinger realization based on the two-mode Wigner algebra

From (23), we have that the Jordan-Schwinger realization of the ordinary
$su(2)$ algebra is given by [40]

$J_{+}=a^{\dagger}_{1}a_{2},\qquad J_{-}=a_{1}a^{\dagger}_{2},\qquad
J_{0}=\frac{1}{2}(N_{1}-N_{2}),$ (29)

with

$[J_{+},J_{-}]=2J_{0},\qquad[J_{0},J_{\pm}]=\pm J_{\pm}.$ (30)

If we now consider the two-mode Wigner oscillator algebra (22)–(23), we will
obviously need more generators to construct a Jordan-Schwinger realization due
to the existence of reflection operators, and we can therefore write

$[J_{+},J_{-}]=2J_{0}-\nu(2N_{2}+1)R_{1}+\nu(2N_{1}+1)R_{2}+2\nu^{2}(R_{2}-R_{1}).$
(31)

Therefore, if we now introduce the following three generators

$P=N_{1}R_{2}-N_{2}R_{1},\qquad K=\frac{1}{2}(R_{2}-R_{1}),\qquad
Q=\frac{1}{2}(R_{2}+R_{1}),$ (32)

then, the algebra $su(2)$ plus the reflection symmetry (i.e., the algebra
$su_{\nu}(2)$) can be described with the above operators as follows

$\displaystyle[J_{0},J_{\pm}]$ $\displaystyle=\pm
J_{\pm},\qquad\qquad\qquad\qquad\quad[J_{+},J_{-}]=2J_{0}+2\nu
P+2\nu(2\nu+1)K,$ $\displaystyle[K,Q]$
$\displaystyle=[K,P]=[K,J_{0}]=0,\qquad\,\\{K,J_{\pm}\\}=0,$
$\displaystyle[Q,P]$ $\displaystyle=[Q,J_{0}]=0,\qquad\qquad\qquad\
\\{Q,J_{\pm}\\}=0,$ (33) $\displaystyle[P,J_{0}]$
$\displaystyle=0,\qquad\qquad\qquad\qquad\qquad\ \ \\{P,J_{\pm}\\}=\pm
2QJ_{\pm}\leavevmode\nobreak\ .$

The representation of the algebra $su_{\nu}(2)$ is then given as

$\displaystyle J_{0}|n_{1},n_{2}\rangle$ $\displaystyle=\frac{n_{1}-n_{2}}{2}\
|n_{1},n_{2}\rangle,$ $\displaystyle P|n_{1},n_{2}\rangle$
$\displaystyle=(n_{1}(-1)^{n_{2}}-n_{2}(-1)^{n_{1}})|n_{1},n_{2}\rangle,$
$\displaystyle Q|n_{1},n_{2}\rangle$
$\displaystyle=\frac{1}{2}\left((-1)^{n_{1}}+(-1)^{n_{2}}\right)\
|n_{1},n_{2}\rangle,$ $\displaystyle K|n_{1},n_{2}\rangle$
$\displaystyle=\frac{1}{2}\left((-1)^{n_{2}}-(-1)^{n_{1}}\right)\
|n_{1},n_{2}\rangle,$ (34) $\displaystyle J_{+}|n_{1},n_{2}\rangle$
$\displaystyle=\sqrt{[n_{1}+1]_{\nu}[n_{2}]_{\nu}}\ |n_{1}+1,n_{2}-1\rangle,$
$\displaystyle J_{-}|n_{1},n_{2}\rangle$
$\displaystyle=\sqrt{[n_{1}]_{\nu}[n_{2}+1]_{\nu}}\ |n_{1}-1,n_{2}+1\rangle,$

with $n_{i}=0,1,2,\dots$ Introducing

$n_{1}=j+m,\quad\qquad n_{2}=j-m,$ (35)

and

$|\widetilde{j,m}\rangle=|j+m,j-m\rangle=|n_{1},n_{2}\rangle,$ (36)

we have the following representation

$\displaystyle J_{0}|\widetilde{j,m}\rangle$
$\displaystyle=m\,|\widetilde{j,m}\rangle,$ $\displaystyle
P|\widetilde{j,m}\rangle$ $\displaystyle=[(j+m)(-1)^{j-m}-(j-m)(-1)^{j+m}]\
|\widetilde{j,m}\rangle,$ $\displaystyle Q|\widetilde{j,m}\rangle$
$\displaystyle=\frac{1}{2}[(-1)^{j-m}+(-1)^{j+m}]\ |\widetilde{j,m}\rangle,$
$\displaystyle K|\widetilde{j,m}\rangle$
$\displaystyle=\frac{1}{2}[(-1)^{j-m}-(-1)^{j+m}]\ |\widetilde{j,m}\rangle,$
(37) $\displaystyle J_{+}|\widetilde{j,m}\rangle$
$\displaystyle=\sqrt{[j+m+1]_{\nu}[j-m]_{\nu}}\ |\widetilde{j,m+1}\rangle,$
$\displaystyle J_{-}|\widetilde{j,m}\rangle$
$\displaystyle=\sqrt{[j+m]_{\nu}[j-m+1]_{\nu}}\ |\widetilde{j,m-1}\rangle.$

This representation has dimension $(2j+1)$, which implies that the allowed
value of $m$ are

$m=-j,-j+1,\dots,j-1,j,\qquad\text{with}\qquad
J_{+}|\widetilde{j,j}\rangle=J_{-}|\widetilde{j,-j}\rangle=0.$ (38)

Therefore, the value of $2j$ is an integer.

####  4.1.1  Representation for odd $2j$

In this case $2m$ is also odd. Thus, we have

$Q=0,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
K=R_{J},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
P=2jR_{J},$ (39)

where we set $R_{J}=(-1)^{j-J_{0}}.$ Thus, the $su_{\nu}(2)$ algebra has three
generators:

$\displaystyle[J_{0},J_{\pm}]$ $\displaystyle=\pm
J_{\pm},\qquad\quad[J_{+},J_{-}]=2J_{0}+2\nu(2\nu+j+1)R_{J},$ (40)
$\displaystyle R_{J}^{2}$ $\displaystyle=I,\qquad\qquad\ \
[R_{J},J_{0}]=0,\qquad\quad\ \\{R_{J},J_{\pm}\\}=0.$ (41)

####  4.1.2 Representation for even $2j$

In this case $2m$ is also even and we may write

$K=0,\qquad\qquad Q=R_{J},\qquad\qquad P=2J_{0}R_{J}.$ (42)

Thus, the $su_{\nu}(2)$ algebra has three generators satisfying

$[J_{0},J_{\pm}]=\pm J_{\pm},\ \ [J_{+},J_{-}]=2J_{0}\left(1+2\nu
R_{J}\right),\ \ R_{J}^{2}=I,\ \ [R_{J},J_{0}]=0,\ \ \\{R_{J},J_{\pm}\\}=0.$
(43)

####  4.1.3 Examples

  * •

$j=1/2$: In this case, we have the following matrix representation

$J_{0}=\begin{pmatrix}1/2&0\cr 0&-1/2\end{pmatrix},\qquad
J_{+}=\begin{pmatrix}0&[1]_{\nu}\cr 0&0\end{pmatrix},\qquad
J_{-}=\begin{pmatrix}0&0\cr[1]_{\nu}&0\end{pmatrix}.$ (44)

If we assign to $J_{0},J_{\pm}$ the $\nu$-deformed Pauli matrices
$\frac{1}{2}\sigma_{z}^{\nu},\sigma_{\pm}^{\nu}$, then we can write

$[\sigma_{+}^{\nu},\sigma_{-}^{\nu}]=(1+3\nu+4\nu^{2})\sigma_{z}^{\nu},$ (45)

where we used $R_{J}=2J_{0}.$

  * •

$j=1$: The matrix representation in this case is

$\hskip-14.22636ptJ_{0}=\begin{pmatrix}1&0&0\cr 0&0&0\cr 0&0&-1\end{pmatrix},\
J_{+}=\begin{pmatrix}0&\sqrt{[2]_{\nu}!}&0\cr 0&0&\sqrt{[2]_{\nu}!}\cr
0&0&0\end{pmatrix},\ J_{-}=\begin{pmatrix}0&0&0\cr\sqrt{[2]_{\nu}!}&0&0\cr
0&\sqrt{[2]_{\nu}!}&0\end{pmatrix}.$ (46)

If we asign $J_{0},J_{\pm}$ to the the $\nu$-deformed angular momentum
generators $L_{z}^{\nu},L_{\pm}^{\nu}$, respectively, we have the following
relationships

$[L_{+}^{\nu},L_{-}^{\nu}]=2L_{z}^{\nu}(1+2\nu L_{z}^{\nu}),$ (47)

where we used $R_{J}=L_{z}^{\nu}.$ The algebra (47) is called quadratic
algebra.

  * •

$j=3/2$: Here, the representation takes the form

$\displaystyle\hskip-8.5359ptJ_{+}=\begin{pmatrix}0&\sqrt{[3]_{\nu}[1]_{\nu}}&0&0\cr
0&0&\sqrt{[2]_{\nu}[2]_{\nu}}&0\cr 0&0&0&\sqrt{[3]_{\nu}[1]_{\nu}}\cr
0&0&0&0\cr\end{pmatrix}\\!,\hskip-8.5359pt$ $\displaystyle
J_{0}=\begin{pmatrix}3/2&0&0&0\cr 0&1/2&0&0\cr 0&0&-1/2&0\cr
0&0&0&-3/2\cr\end{pmatrix}\\!,\ \ \ \ $
$\displaystyle\hskip-8.5359ptJ_{-}=\begin{pmatrix}0&0&0&0\cr\sqrt{[3]_{\nu}[1]_{\nu}}&0&0&0\cr
0&\sqrt{[2]_{\nu}[2]_{\nu}}&0&0\cr
0&0&\sqrt{[3]_{\nu}[1]_{\nu}}&0\cr\end{pmatrix}\\!,\hskip-8.5359pt$
$\displaystyle R_{J}=\begin{pmatrix}1&0&0&0\cr 0&-1&0&0\cr 0&0&1&0\cr
0&0&0&-1\cr\end{pmatrix}\\!.$ (48)

  * •

$j=2$: In this case, the matrices appear in the form

$\displaystyle J_{+}$
$\displaystyle=\begin{pmatrix}0&\sqrt{[4]_{\nu}[1]_{\nu}}&0&0&0\cr
0&0&\sqrt{[3]_{\nu}[2]_{\nu}}&0&0\cr 0&0&0&\sqrt{[3]_{\nu}[2]_{\nu}}&0\cr
0&0&0&0&\sqrt{[4]_{\nu}[1]_{\nu}}\cr 0&0&0&0&0\cr\end{pmatrix},$
$\displaystyle J_{-}$
$\displaystyle=\begin{pmatrix}0&0&0&0&0\cr\sqrt{[4]_{\nu}[1]_{\nu}}&0&0&0&0\cr
0&\sqrt{[3]_{\nu}[2]_{\nu}}&0&0&0\cr 0&0&\sqrt{[3]_{\nu}[2]_{\nu}}&0&0\cr
0&0&0&\sqrt{[4]_{\nu}[1]_{\nu}}&0\cr\end{pmatrix},$ (49) $\displaystyle J_{0}$
$\displaystyle=\begin{pmatrix}2&0&0&0&0\cr 0&1&0&0&0\cr 0&0&0&0&0\cr
0&0&0&-1&0\cr 0&0&0&0&-2\cr\end{pmatrix},\qquad\qquad
R_{J}=\begin{pmatrix}1&0&0&0&0\cr 0&-1&0&0&0\cr 0&0&1&0&0\cr 0&0&0&-1&0\cr
0&0&0&0&1\cr\end{pmatrix}.$

###  4.2 The Holstein-Primakoff realization of $su_{\nu}(2)$ algebra

For the ordinary $su(2)$ algebra, the Holstein-Primakoff single-mode
realization has the form [1]

$J_{+}=\sqrt{2j-N}a,\qquad J_{-}=a^{\dagger}\sqrt{2j-N},\qquad J_{0}=j-N.$
(50)

The latter is related to the $(2j+1)$-dimensional Fock space spanned by
$|n\rangle$, which obeys $N|n\rangle=n|n\rangle,n=0,1,\dots,2j$. A comparison
with equation (29) gives

$\frac{N_{1}-N_{2}}{2}\Longleftrightarrow j-N,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ R_{1}=R_{2}=R.$ (51)

Thus, the algebra $su_{\nu}(2)$ reads

$[J_{+},J_{-}]=2J_{0}(1+2\nu R),\qquad[J_{0},J_{\pm}]=\pm J_{\pm}.$ (52)

Now, we can find the Holstein-Primakoff realization of algebra $su_{\nu}(2)$,
which is given by

$\displaystyle J_{0}$ $\displaystyle=j-N,\qquad\qquad\qquad\qquad R=(-1)^{N},$
$\displaystyle J_{+}$
$\displaystyle=\sqrt{\frac{(2j-N)(N+1)+\nu(1+2j+R(-1+2j-2N))}{N+1+\nu(1+R)}}\
a,$ (53) $\displaystyle J_{-}$ $\displaystyle=a^{\dagger}\
\sqrt{\frac{(2j-N)(N+1)+\nu(1+2j+R(-1+2j-2N))}{N+1+\nu(1+R)}}.$

##  5 The $so_{\nu}(3)$ algebra

Finally, in this section we are going to give some new results about the
deformed algebra $so_{\nu}(3)$, which is obtained from the algebra
$su_{\nu}(2)$ considering

$L_{z}=J_{0},\qquad L_{x}=\frac{1}{2}(J_{+}+J_{-}),\qquad
L_{y}=\frac{i}{2}(J_{-}-J_{+}).$ (54)

As a result, the algebra $so_{\nu}(3)$ is given by

$\displaystyle[L_{z},L_{x}]$
$\displaystyle=iL_{y},\qquad[L_{z},L_{y}]=-iL_{x},$
$\displaystyle[L_{x},L_{y}]$ $\displaystyle=2L_{z}+2\nu P+2\nu(2\nu+1)K,$
$\displaystyle[K,Q]$
$\displaystyle=[K,P]=[K,L_{z}]=0,\qquad\\{K,L_{x}\\}=\\{K,L_{y}\\}=0,$ (55)
$\displaystyle[Q,P]$
$\displaystyle=[Q,L_{z}]=0,\qquad\\{Q,L_{x}\\}=\\{Q,L_{y}\\}=0,$
$\displaystyle[P,L_{z}]$
$\displaystyle=0,\qquad\\{P,L_{x}\\}=2iQL_{y},\qquad\\{P,L_{y}\\}=-2iQL_{x}.$

###  5.1  Representation for odd $2j$

In this case, denoting $R_{L}=(-1)^{j-L_{z}}.$, the $so_{\nu}(3)$ algebra has
three generators that satisfy:

$[L_{z},L_{x}]=iL_{y},\qquad[L_{z},L_{y}]=-iL_{x},\qquad[L_{x},L_{y}]=2L_{z}+2\nu(2\nu+j+1)R_{L}.$
(56)

###  5.2  Representation for even $2j$

Now the algebra $so_{\nu}(3)$ also has three generators:

$[L_{z},L_{x}]=iL_{y},\qquad[L_{z},L_{y}]=-iL_{x}\qquad[L_{x},L_{y}]=2L_{z}(1+2\nu
R_{L}).$ (57)

##  6 Conclusions

Starting from very recent applications of some concepts from various fields of
theoretical physics and quantum technologies, we have developed a new parity-
deformed type of Wigner algebras, these being the essential reasons for
introducing these new mathematical structures:

  * •

The broad application of parity-deformed systems in quantum information and
communication, including the study of coherent and compressed states, cat
states, data transfer, and error correction.

  * •

The application of the Holstein-Primakoff and Jordan-Schwinger realizations in
various fields, including quantum communication.

  * •

The usefulness of $SO(3)$ group representations in the study of multi-qubit
systems.

  * •

The experimental limitations and the importance of theoretical calculations in
the modeling and preparation of experiments in the areas already mentioned.

  * •

The application of the groups considered in this work in the study of qutrits.

  * •

The possible application and generalization of traditional light-matter
Hamiltonians, including the conventional Rabi model and its recent
generalizations.

As a result of this work, and based on the single mode and two mode Wigner
algebra, the deformed algebra $su(2)$ ($su_{\nu}(2)$) and the deformed algebra
$so(3)$ ($so_{\nu}(3)$) were constructed. First, the single-mode Wigner
algebra and two types of its realizations were analyzed: one based on the
monomial basis and another based on the graded quasi-polynomial basis $Z_{2}$.
The single-mode Wigner algebra was then extended to the two-mode case, which
was used to obtain the algebra $su_{\nu}(2)$. The following results were
found: (a) for $j=1/2$ we have the same algebra as $su(2)$, (b) the case $j=1$
produces the quadratic algebra and (c) the cases $j=3/2,2,5/2,\dots$ result in
the algebra $su_{\nu}(2)$, which is neither the algebra $su(2)$ nor the
quadratic algebra. The Holstein-Primakoff realization of the algebra
$su_{\nu}(2)$ was also discussed, and the algebra $so_{\nu}(3)$ was derived
from the algebra $su_{\nu}(2)$.

Having prepared the required deformed basis, our next goal is to generalize
the present work and investigate its consequences on the temporal evolution of
certain related quantum systems, including the analysis of the fidelity of
entanglement transfer as well as the entanglement transfer in the so-called
${J_{x}}$ photonic networks. We hope that the present formulation helps
generalize the conventional Rabi model, a topic we are currently working on.

## Aknowledgments

The research of L.M.N. and S.Z. was supported by Q-CAYLE project, funded by
the European Union-Next Generation UE/MICIU/Plan de Recuperacion,
Transformacion y Resiliencia/Junta de Castilla y Leon (PRTRC17.11), and also
by RED2022-134301-T and PID2020-113406GB-I00, both financed by
MICIU/AEI/10.13039/501100011033. The authors are grateful to Dr. I.A.
Bocanegra-Garay and Dr. D.G. Bussandri for useful comments on the manuscript.

## References

  * [1] T. Holstein and H. Primakoff, Field Dependence of the Intrinsic Domain Magnetization of a Ferromagnet, Phys. Rev. 58 (1940) 1058.
  * [2] E. P. Wigner, Do the equations of motion determine the quantum mechanical commutation relations?, Phys. Rev. 77 (1950) 711.
  * [3] L. Yang, A note on the quantum rule of the harmonic oscillator, Phys. Rev. 84 (1951) 788.
  * [4] M. T. Raiford, Statistical dynamics of quantum oscillators and parametric amplification in a single mode, Phys. Rev. A 24 (1970) 1541.
  * [5] C. F. Dunkl, Integral kernels with reflection group invariance, Can. J. Math. 43 (1991)1213. C. F. Dunkl and Y. Xu (2014) Orthogonal Polynomials of Several Variables (Encyclopedia of Mathematics and its Applications) 2nd edn (Cambridge: Cambridge University Press)
  * [6] M. Arik and D. D. Coon, Hilbert spaces of analytic functions and generalized coherent states, J. Math. Phys. 17 (1976) 524.
  * [7] L.C. Biedenharn, The quantum group $SU_{q}(2)$ and a $q$-analogue of the boson operators, J. Phys. A 22 (1989) L873.
  * [8] J. Macfarlane, On $q$ Analogs of the Quantum Harmonic Oscillator and the Quantum Group $SU_{q}(2)$, J. Phys. A 22 (1989) 4581.
  * [9] T.D. Palev and N. I. Stoilova, Wigner Quantum Systems (Lie superalgebraic approach), Rep. Math. Phys. 49 (2002) 395. R. C. King and T. D. Palev, N.I. Stoilova and J. Van der Jeugt, The noncommutative and discrete space structure of the $sl(1|3)$ Wigner quantum oscillator, J. Phys. A 36 (2003) 4337. R.C. King et al., A noncommutative $n$-particle 3D Wigner quantum oscillator, J. Phys. A 36 (2003) 1199.
  * [10] P. A. Horvathy, M. S. Plyushchay and M. Valenzuela, Bosonized supersymmetry of anyons and supersymmetric exotic particle on the non-commutative plane, Nucl. Phys. B 768 (2007) 247.
  * [11] M. S. Plyushchay, Deformed Heisenberg algebra with reflection, Nucl. Phys. B 491 (1997) 619.
  * [12] P. A. Horvathy, M. S. Plyushchay and M. Valenzuela, Bosons, fermions and anyons in the plane, and supersymmetry, Ann. Phys. 325 (2010) 1931.
  * [13] V. Jakubský, L. M. Nieto and M. S. Plyushchay, The origin of the hidden supersymmetry, Phys. Lett. B 692 (2010) 51.
  * [14] A. Dehghani et al., Cat-states in the framework of Wigner–Heisenberg algebra, Ann. Phys. 362 (2015) 659.
  * [15] A. Dehghani et al., Parity Deformed Jaynes-Cummings Model: “Robust Maximally Entangled States”, Sci. Rep. 6 (2016) 38069.
  * [16] L. Inzunza and M. S. Plyushchay, Hidden superconformal symmetry: Where does it come from?, Phys. Rev. D 97 (2018) 045002.
  * [17] A. Dehghani et al., Entanglement transfer in a noisy cavity network with parity-deformed fields, J. Opt. Soc. Am. B 36 (2019) 1858.
  * [18] S. Sargolzaeipour, H. Hassanabadi and W. S. Chung, Effect of the Wigner Dunkl algebra on the Dirac equation and Dirac harmonic oscillator, Mod. Phys. Lett. A 25 (2018) 1850146. W. S. Chung and H. Hassanabadi, One-dimensional quantum mechanics with Dunkl derivative, Mod. Phys. Lett. A 34 (2019) 1950190.
  * [19] R. Mohseninia et al., Always-On Quantum Error Tracking with Continuous Parity Measurements, Quantum 4 (2020) 358.
  * [20] A. Ballesteros, F. Giacomini and G. Gubitosi, The group structure of dynamical transformations between quantum reference frames, Quantum 5 (2021) 470.
  * [21] L. M. Hurtado-Molina et al., Formation of spectral triplets induced by parity deformation in a quantum dot–cavity system, Physica B 604 (2021) 412698.
  * [22] H. Fakhri, S. Mirzaei and M. Sayyah-Fard, Two-photon Jaynes–Cummings model: a two-level atom interacting with the para-Bose field, Quant. Inf. Proc. 20 (2021) 398.
  * [23] H. Fakhri and M. Sayyah-fard, Para-Bose oscillator algebras of odd orders: $x$-representations and Wigner functions for coherent and cat states and their photon-added and photon-subtracted counterparts, Eur. Phys. J. Plus 137 (2022) 1199.
  * [24] S. H. Dong et al., Supersymmetric Wigner–Dunkl quantum mechanics, Results Phys. 39 (2022) 105664.
  * [25] M. Drieb-Schön et al., Parity Quantum Optimization: Encoding Constraints, Quantum 7 (2023) 951.
  * [26] M. Fellner et al., Parity Quantum Optimization: Benchmarks, Quantum 7 (2023) 952.
  * [27] R. Raussendorf et al., The role of cohomology in quantum computation with magic states, Quantum 7 (2023) 979.
  * [28] B. Mojaveri and J. Taghipour, Entanglement protection of two qubits moving in an environment with parity-deformed fields, Eur. Phys. J. Plus 138 (2023) 263
  * [29] G. Junker et al., On the gauge invariance of Wigner–Dunkl quantum mechanics in the presence of a constant magnetic field, Ann. Phys. 454 (2023) 169336.
  * [30] P. Sedaghatnia et al., Investigation of the generalised Wigner–Dunkl harmonic oscillator and its coherent states, Ann. Phys. 458 (2023) 169445.
  * [31] X. Y. Zhu et al., Fast high-fidelity parity measurement of spin qubits in quantum dots using modulated pulses, Phys. Rev. A 107 (2023) 033708.
  * [32] Y. Lu et al., High-fidelity parametric beam splitting with a parity-protected converter, Nat. Commun. 14 (2023) 5767.
  * [33] V. Buzek, Jaynes-Cummings model with intensity-dependent coupling interacting with Holstein-Primakoff $SU(1,1)$ coherent state, Phys. Rev. 39 (1989) 3196,
  * [34] M. Chaichian, D. Ellinas and P. Kulish, Quantum Algebra as the Dynamical Symmetry of the Deformed Jaynes-Cummings Model, Phys. Rev. Lett. 65 (1990) 980.
  * [35] V. I. Koroli and V. V. Zalamai, Dynamics of a laser-cooled and trapped radiator interacting with the Holstein–Primakoff $SU(1,1)$ coherent state, J. Phys. B: At. Mol. Opt. Phys. 42 (2009) 035505.
  * [36] J. G. Hirsch et al., Virtues and limitations of the truncated Holstein–Primakoff description of quantum rotors, Phys. Scr. 87 (2013) 038106.
  * [37] U. Bajpai, A. Suresh and B. K. Nikolic, Quantum many-body states and Green’s functions of equilibrium electron-magnon systems: Localized spin operators versus their mapping to Holstein-Primakoff bosons, Phys. Rev. B 104 (2021) 184425.
  * [38] C. P. Sun et al., Quantum information transfer for two nondegenerate modes, Phys. Rev. A 73 (2006) 022318. C. P. Sun et al., Quantum transducers: Integrating transmission lines and nanomechanical resonators via charge qubits, Phys. Rev. A 73 (2006) 022318.
  * [39] T. Skrypnyk, Integrability and superintegrability of the generalized two-level many-mode Jaynes–Cummings and Dicke models, J. Math. Phys. 50 (2009) 103523.
  * [40] V. Sunil Kumar, A. B. Bindu and R. Jagannathan, Jordan–Schwinger-type realizations of three-dimensional polynomial algebras, Mod. Phys. Lett. A 17 (2002) 1559.
  * [41] F. De Zola, Topological phase for entangled two-qubit states and the representation of the $SO(3)$ group, J. Opt. B: Quantum Semiclass. Opt. 7 (2005) 372.
  * [42] A. R. P. Rau, Symmetries and Geometries of Qubits, and Their Uses, Symmetry 13 (2021) 1732.
  * [43] D. B. Uskov and P. M. Alsing, Vector properties of entanglement in a three-qubit system, Phys. Rev. A 102 (2020) 032401.
  * [44] A. Perej-Leija et al., Perfect transfer of path-entangled photons ${J_{x}}$ photonic lattices, Phys. Rev. A 87 (2013) 022303.