# Combination of the Perturbation Theory with Configuration Interaction Method

M. G. Kozlov1,2, I. I. Tupitsyn1,3, A. I. Bondarev1,4, and D. V. Mironova2
1Petersburg Nuclear Physics Institute of NRC “Kurchatov Institute”, 188300
Gatchina, Russia  2St. Petersburg Electrotechnical University “LETI”, Prof.
Popov Str. 5, 197376 St. Petersburg, Russia  3Department of Physics, St.
Petersburg State University, Ulianovskaya 1, Petrodvorets, 198504 St.
Petersburg, Russia  4Center for Advanced Studies, Peter the Great St.
Petersburg Polytechnic University, Polytekhnicheskaja 29, 195251 St.
Petersburg Russia

###### Abstract

Present atomic theory provides accurate and reliable results for atoms with a
small number of valence electrons. However, most current methods of
calculations fail when the number of valence electrons exceeds four or five.
This means that we can not make reliable predictions for more than a half of
the periodic table. Here we suggest a modification of the CI+MBPT
(configuration interaction plus many-body perturbation theory) method, which
may be applicable to atoms and ions with filling $d$ and $f$ shells.

##  I Introduction

At present there are several methods of the relativistic correlation
calculations of atoms, such as multiconfiguration Dirac-Fock Desclaux (1975);
Jönsson et al. (2013); Froese Fischer et al. (2016), configuration interaction
(CI) Gu (2008); Fritzsche et al. (2002); Jiang et al. (2016); Tupitsyn et al.
(2003, 2005), many-body perturbation theory (MBPT) Dzuba et al. (1987);
Blundell et al. (1987); Sapirstein (1998), CI+MBPT Dzuba et al. (1996); Kahl
and Berengut (2019); Kozlov et al. (2015), coupled cluster Eliav and Kaldor
(2010); Saue et al. (2020); Blundell et al. (1989); Safronova and Johnson
(2008); Oleynichenko et al. (2020), and others. Calculations are usually done
in the no-virtual-pair approximation using Dirac-Coulomb, or Dirac-Coulomb-
Breit approximations Johnson (2007). QED corrections may be included using
radiative potential method developed by Flambaum and Ginges (2005), Ginges and
Berengut (2016), and QEDMOD potential Shabaev et al. (2013); Tupitsyn et al.
(2016).

The coupled cluster method is one of the most popular and effective methods
for calculation of atoms with a small number of open shell electrons (or
holes). Calculations of the spectra of atoms and ions with many valence
electrons (e. g. transition metals, lanthanides, and actinides) are very
difficult and usually not very accurate. The reason for that is a combination
of strong correlations and a very large configuration space. To account for
strong correlations one needs non-perturbative methods, such as CI. On the
other hand, a large configuration space makes such calculations very
expensive. As a compromise one can try to combine CI with perturbation theory
(PT). We will first assume that all closed atomic shells are considered
frozen. Then we are treating only valence correlations and consider a
combination of the valence CI with valence perturbation theory (VPT). Later we
will see that this approach can also be used to treat core-valence
correlations.

Recently there were several attempts Dzuba et al. (2017); Geddes et al.
(2018); Dzuba et al. (2019); Imanbaeva and Kozlov (2019) to develop an
effective and fast CI+VPT method to speed up calculations for such systems,
where straightforward CI calculations are impossible. Application of these
methods for systems with a large number of valence electrons was demonstrated
in Refs. Cheung et al. (2020); Li and Dzuba (2020). A general idea of all
these calculation schemes is to make CI in a smaller subspace $P$ and
calculate corrections from a complementary subspace $Q$ using VPT. In Refs.
Dzuba et al. (2017); Geddes et al. (2018); Dzuba et al. (2019) it is suggested
to neglect non-diagonal blocks of the CI matrix in the subspace $Q$, which is
equivalent to using VPT. All these methods require summation over all
determinants of the complementary subspace $Q$. Though calculating this sum is
much easier than calculating and diagonalizing the whole CI matrix, it is
still too expensive for the number of valence electrons approaching, or
exceeding ten.

In the paper Dzuba et al. (2019), the sum over determinants was partly
substituted by the sum over configurations that led to a significant increase
in calculation speed. Here we want to make another step in this direction. To
this end, we will partly substitute VPT with many-body perturbation theory
(MBPT). The method we propose here is similar to the old CI+MBPT method Dzuba
et al. (1996) but uses different splitting of the problem into the CI and MBPT
parts. In particular, we suggest to account for double excitations (D) from
the subspace $P$ by means of the MBPT and treat single excitations (S) within
VPT, or, if possible, include them directly in CI. We think that this variant
is not only more efficient for treating valence correlations, but may also be
used for the core-valence correlations.

##  II Formalism

###  II.1 Valence correlations

Consider many-electron atom, or ion with $N$ valence electrons, where $N\gg
1$. Let us first assume that other electrons always occupy closed core shells,
which is known as a frozen core approximation. Our aim is to solve the
$N$-electron Schrödinger equation and find the spectrum of this system.

We start with splitting $N$-electron configuration space in two orthogonal
subspaces $P$ and $Q$. The subspace $P$, which we call valence, includes the
most important shells. It may be not obvious from the start, which orbitals
are ‘important’. We definitely must include into subspace $P$ all orbitals
with occupation numbers of the order of unity in the physical states, we are
interested in. Complementary subspace $Q$ includes S, D, and so on excitations
from the valence shells to the virtual ones, thus, $Q=Q_{S}+Q_{D}+\dots$. We
start by solving the matrix equation in the subspace $P$,

$\displaystyle\hat{P}H\hat{P}\Psi_{a}$ $\displaystyle=E_{a}\hat{P}\Psi_{a}\,,$
(1)

where $H$ is the Hamiltonian for valence electrons and $\hat{P}$ is the
projector on the subspace $P$. We can find a correction from the complementary
subspace $Q$ using the second-order perturbation theory:

$\displaystyle\delta E_{a}$ $\displaystyle=\sum_{n\in
Q}\frac{\langle\Psi_{a}|\hat{P}H\hat{Q}|n\rangle\langle
n|\hat{Q}H\hat{P}|\Psi_{a}\rangle}{E_{a}-E_{n}},$ (2)

where $|n\rangle$ are $N$-electron determinants in the complementary subspace
$Q$ and $E_{n}=\langle n|\hat{Q}H\hat{Q}|n\rangle$.

The wavefunction $\Psi_{a}$ is a linear combination of the determinants:

$\displaystyle\Psi_{a}$ $\displaystyle=\sum_{m\in
P}C^{a}_{m}|m\rangle=\sum_{p,m_{p}}C^{a}_{p,m_{p}}|m_{p}\rangle\,.$ (3)

Here and below indexes $p$ and $q$ run over configurations in the subspaces
$P$ and $Q$ respectively and indexes $m_{p}$ and $n_{q}$ numerate determinants
within one configuration. Now Eq. (2) takes the form:

$\delta
E_{a}=\sum_{p,m_{p}}\sum_{p^{\prime},m_{p^{\prime}}}C^{a}_{p,m_{p}}C^{a}_{p^{\prime},m_{p^{\prime}}}\\\
\times\sum_{q,n_{q}}\frac{\langle m_{p}|H|n_{q}\rangle\langle
n_{q}|H|m_{p^{\prime}}\rangle}{E_{a}-E_{n_{q}}},$ (4)

where the sum over the subspace $Q$ is also split in two.

For an atom with $N\sim 10$, the dimension of space $Q$ is very large, which
makes evaluation of expression (4) very lengthy. Therefore, our aim is to
substitute double sum over $q$ and $n_{q}$ by a single sum over $q$. To this
end we do the following approximation: we substitute the energy $E_{n_{q}}$ in
the denominator by the configuration average:

$\displaystyle\bar{E}_{q}$
$\displaystyle=\frac{1}{N_{q}}\sum_{n_{q}=1}^{N_{q}}E_{n_{q}}\,,$ (5)

where $N_{q}$ is the number of determinants in configuration $q$. Using this
approximation we rewrite (4) in a form:

$\delta
E_{a}=\sum_{p,m_{p}}\sum_{p^{\prime},m_{p^{\prime}}}C^{a}_{p,m_{p}}C^{a}_{p^{\prime},m_{p^{\prime}}}\\\
\times\sum_{q}\frac{\langle m_{p}|H\left({\sum_{n_{q}}|n_{q}\rangle\langle
n_{q}|}\right)H|m_{p^{\prime}}\rangle}{E_{a}-\bar{E}_{q}}\,.$ (6)

Below we will show that in some very important cases one can get rid of the
internal sum over $n_{q}$.

Hamiltonian $H$ includes one-particle and two-particle parts. The former
consists of the kinetic term and the core potential, while the latter
corresponds to the Coulomb (or Coulomb-Breit) interaction between valence
electrons. Thus, in the sum over $q$ remain only configurations, which differ
by no more than two electrons from configurations $p$ and $p^{\prime}$. This
means that within this approximation the subspace $Q$ is actually truncated to
$Q_{S}+Q_{D}$. All non-zero contributions correspond to the diagrams, shown in
Fig. 1.

Figure 1: Set of connected second-order diagrams. Black dots correspond to the
core potential and wavy lines to the Coulomb interaction. Double and single
lines denote electrons in valence and virtual orbitals respectively. Non-
symmetric diagrams $(b)$ and $(e)$ have mirror twins.

According to our definition of the spaces $P$ and $Q$, the latter must include
at least one electron in the virtual shell. Diagrams $(a)$, $(b)$, and $(e)$
include only one intermediate line, so they describe single excitations from
the subspace $P$. Diagrams $(c)$ and $(d)$ include two intermediate lines, but
only diagram $(d)$ describes double (D) excitations, as both intermediate
lines correspond to the virtual shells.

Figure 1 shows that all many-electron matrix elements in Eq. (6) are reduced
to the effective one-electron, two-electron, and three-electron contributions.
Effective one-electron contributions are described by diagram $(a)$; diagrams
$(b)$, $(c)$, and $(d)$ correspond to the two-electron contributions; finally,
diagram $(e)$ describes effective three-electron contributions, see Figs. 2
and 3.

Figure 2: Many-electron second-order expression in Eq. (7) (left) is reduced
to the two-particle expression (middle), which, in turn, is reduced to the
effective two-particle interaction (right).  Figure 3: The case when many-
electron second-order expression in Eq. (6) (left) is reduced to the three-
particle expression (middle), which, in turn, is reduced to the effective
three-particle interaction (right). The initial configuration on the left
differs from the intermediate configuration by the upper two electrons. The
final configuration differs from the intermediate one by the second and third
electron from the top.

For combinatorial reasons the number of configurations with two excited
electrons is much bigger, than the number of those with only one such
electron. Therefore the vast majority of terms in Eq. (6) correspond to the
two-electron excitations from configurations $p$ and $p^{\prime}$. For these
terms in the Hamiltonian $H$ only the two-electron interaction $V$ can
contribute, so we can neglect the one-electron part and make substitution
$H\to V$. As we saw above, all such terms are described by the single diagram
$(d)$ from Fig. 1.

Let us consider the sum over doubly excited configurations. It can be written
as:

$\delta
E_{a}^{D}=\sum_{p,m_{p}}\sum_{p^{\prime},m_{p^{\prime}}}C^{a}_{p,m_{p}}C^{a}_{p^{\prime},m_{p^{\prime}}}\\\
\times\sum_{q\in Q_{D}}\frac{\langle
m_{p}|V\left({\sum_{n_{q}}|n_{q}\rangle\langle
n_{q}|}\right)V|m_{p^{\prime}}\rangle}{E_{a}-\bar{E}_{q}}\,,$ (7)

Non-zero contributions come from determinants $|n_{q}\rangle$, which differ
from both determinants $|m_{p}\rangle$ and $|m_{p^{\prime}}\rangle$ by two
electrons. It is clear that it must be the same two electrons, see Fig. 2. In
this case the second-order expression from (7) is reduced to the effective
two-particle interaction Dzuba et al. (1996):

$\displaystyle\delta E_{a}^{D}$
$\displaystyle=\sum_{p,m_{p}}\sum_{p^{\prime},m_{p^{\prime}}}C^{a}_{p,m_{p}}C^{a}_{p^{\prime},m_{p^{\prime}}}\langle
m_{p}|V_{\mathrm{eff}}|m_{p^{\prime}}\rangle\,.$ (8)

This effective interaction can be expressed in terms of the effective radial
integrals, which are similar to the Coulomb radial integrals. The latter
appear when we expand Coulomb interaction in spherical multipoles,

$\displaystyle V=\sum_{k,\varkappa}V^{k}_{\varkappa}.$ (9)

The matrix element of each multipole component $V^{k}_{\varkappa}$ has the
form:

$\langle
c,d|V_{\varkappa}^{k}|a,b\rangle=(-1)^{m_{c}+m_{b}+1}\delta_{p}\sqrt{(2j_{a}+1)(2j_{b}+1)(2j_{c}+1)(2j_{d}+1)}\\\
\left(\\!\begin{array}[]{ccc}j_{c}&j_{a}&k\\\ -m_{c}&m_{a}&\varkappa\\\
\end{array}\\!\right)\left(\\!\begin{array}[]{ccc}j_{b}&j_{d}&k\\\
-m_{b}&m_{d}&\varkappa\\\
\end{array}\\!\right)\left(\\!\begin{array}[]{ccc}j_{c}&j_{a}&k\\\
\frac{1}{2}&-\frac{1}{2}&0\\\
\end{array}\\!\right)\left(\\!\begin{array}[]{ccc}j_{b}&j_{d}&k\\\
\frac{1}{2}&-\frac{1}{2}&0\\\ \end{array}\\!\right)R^{k}_{a,b,c,d}\,,$ (10)

where round brackets denote 3j-symbols, $R^{k}_{a,b,c,d}$ denotes the Coulomb
radial integral and $\delta_{p}$ ensures parity selection rule:
$\delta_{p}=\xi(l_{a}+l_{c}+k)\xi(l_{b}+l_{d}+k)$ and $\xi(n)=1,0$ for
$n=\mbox{even, odd}$. A similar multipole expansion holds for the effective
interaction $V_{\mathrm{eff}}$, the effective radial integral being Dzuba et
al. (1996):

$R^{k,\mathrm{eff}}_{a,b,c,d}=\sum_{k_{1},k_{2}}\sum_{m,n}(-1)^{\chi}(2j_{m}+1)(2j_{n}+1)(2k+1)\left\\{\\!\begin{array}[]{ccc}j_{c}&j_{a}&k\\\
k_{1}&k_{2}&j_{m}\\\
\end{array}\\!\right\\}\left\\{\\!\begin{array}[]{ccc}j_{b}&j_{d}&k\\\
k_{2}&k_{1}&j_{n}\\\
\end{array}\\!\right\\}\left(\\!\begin{array}[]{ccc}j_{m}&j_{a}&k_{1}\\\
\frac{1}{2}&-\frac{1}{2}&0\\\
\end{array}\\!\right)\left(\\!\begin{array}[]{ccc}j_{b}&j_{n}&k_{1}\\\
\frac{1}{2}&-\frac{1}{2}&0\\\ \end{array}\\!\right)\\\
\left(\\!\begin{array}[]{ccc}j_{c}&j_{m}&k_{2}\\\
\frac{1}{2}&-\frac{1}{2}&0\\\
\end{array}\\!\right)\left(\\!\begin{array}[]{ccc}j_{n}&j_{d}&k_{2}\\\
\frac{1}{2}&-\frac{1}{2}&0\\\
\end{array}\\!\right)\left(\\!\begin{array}[]{ccc}j_{c}&j_{a}&k\\\
\frac{1}{2}&-\frac{1}{2}&0\\\
\end{array}\\!\right)^{-1}\\!\\!\left(\begin{array}[]{ccc}j_{b}&j_{d}&k\\\
\frac{1}{2}&-\frac{1}{2}&0\\\
\end{array}\right)^{-1}\frac{R^{k_{1}}_{a,b,m,n}R^{k_{2}}_{c,d,m,n}}{\Delta_{E}}\,,$
(11)

where curly brackets denote 6j-coefficients, the phase
$\chi={j_{a}+j_{b}+j_{c}+j_{d}+j_{m}+j_{n}+k_{1}+k_{2}+k+1}$, and $\Delta_{E}$
is energy denominator, which we will discuss later. For the effective
interaction there is no link between parity and multipolarity $k$, so for
$V_{\mathrm{eff}}$ we do not have factor $\delta_{p}$ as in Eq. (10). The sum
in (11) runs over multipolarities $k_{1}$ and $k_{2}$, which satisfies the
triangle rule $|k_{1}-k_{2}|\leq k\leq k_{1}+k_{2}$.

All single excitations are described by the remaining diagrams from Fig. 1.
The diagram $(a)$ has a form of the effective one-electron radial integral,
while diagrams $(b)$ and $(c)$ are reduced to the two-electron effective
radial integrals. In principle, these effective radial integrals can be
calculated and stored. However, the diagram $(e)$ corresponds to the effective
three-particle interaction. It is difficult to include such interactions into
CI matrix for several reasons:

  * •

When $N>3$ the number of such effective three-particle integrals is huge.

  * •

It is difficult to store them and find them.

  * •

The number of the non-zero matrix elements in the matrix drastically
increases. The matrix becomes less sparse and its diagonalization is much more
difficult and time-consuming.

Because of all that it is inefficient to use the MBPT approach for three-
particle diagrams and it is much easier to treat them within the determinant-
based PT. However, it is difficult then to separate them from other
contributions, which correspond to single excitations. Thus, it is better not
to use MBPT for single excitations at all. We suggest to use instead any form
of the determinant-based VPT described in Refs. Dzuba et al. (2017); Geddes et
al. (2018); Imanbaeva and Kozlov (2019); Dzuba et al. (2019). This means that
we do VPT in the subspace $Q_{S}$. Note that the dimension of this subspace is
incomparably smaller than the dimension of the $Q_{D}$ subspace. In some cases
it may be so small, that we can include $Q_{S}$ in the subspace $P$, where we
do CI.

###  II.2 Core-valence correlations

Figure 4: Set of one-electron second-order diagrams accounting for the
excitations from the core. Diagrams $(e)$ and $(f)$ have mirror twins.
Diagrams $(c)$ and $(d)$ describe double excitations from the core.

It is easy to use the scheme described above for the core-valence correlations
as well. Now $P$ subspace corresponds to the frozen-core approximation and the
subspaces $Q_{S}$ and $Q_{D}$ include single and double excitations from the
core respectively. This means that these subspaces include many-electron
states with one and two holes in the core. As before, the second-order MBPT
corrections are described by one-electron, two-electron, and three-electron
diagrams. All one-electron diagrams are given in Fig. 4. Excitations from the
core correspond to the hole lines with arrows looking to the left. It is easy
to see that only diagrams $(c)$ and $(d)$ describe double excitations.
Therefore, we need to calculate them and store as one-electron effective
radial integrals, see Fig. 5 (note that there are no one-electron
contributions for the valence excitations). Expressions for these diagrams
were given in Ref. Dzuba et al. (1996).

There is only one two-electron diagram, which corresponds to the double
excitations from the core. This diagram must be calculated and added to the
similar diagram for valence excitations, which was discussed in the previous
section, see Fig. 6. Finally, in analogy with the valence correlations, the
three-particle diagrams correspond to the single excitations from the core.

Figure 5: Diagrams, which correspond to the double excitations from closed
shells. These diagrams are described by the effective one-electron radial
integrals, designated by a black circle.  Figure 6: Diagrams contributing to
the effective two-electron radial integrals. First diagram accounts for the
double excitations to the virtual shells and second diagram accounts for the
double excitations from closed shells.

We conclude that in order to account for both valence and core-valence
correlations we need to calculate one-electron and two-electron effective
radial integrals, which corresponds to the diagrams from Figs. 5 and 6. At the
same time, we need to include all single excitations from the core shells and
all single excitations to the virtual shells either in the subspace $P$ or in
the subspace $Q_{S}$. After that, we make CI calculation with effective radial
integrals possibly followed by the VPT calculation in the $Q_{S}$ subspace.

###  II.3 Sketch of the possible calculation scheme

Let us describe a most general computational scheme.

  * •

Basis set orbitals are divided into four groups: inner core, outer core,
valence, and virtual orbitals. The inner core is kept frozen on all stages of
calculation.

  * •

Effective radial integrals are calculated for the valence orbitals, which
account for the double excitations from the outer core and the double
excitations from the valence orbitals to the virtual ones.

  * •

Full CI calculation is done for the valence electrons. The effective radial
integrals are added to the conventional radial integrals when the Hamiltonian
matrix is formed.

  * •

Determinant-based PT is used in the complementary subspace $Q_{S}$, which
includes single excitations from the outer core and single excitations to the
virtual states.

Depending on the number of the valence electrons and the size of the core this
scheme can be simplified. If there are only two valence electrons, one can
include all virtual basis states into valence space. Single excitation from
the core can be also added to the valence space. Double excitations from the
core are accounted for through the effective radial integrals, while single
excitations are included explicitly in the CI matrix. Formally this means that
we substitute $P$, $Q$ decomposition by the $P^{\prime}$, $Q_{D}$
decomposition:

$\displaystyle P+Q=P+Q_{S}+Q_{D}=P^{\prime}+Q_{D}\,,$ (12) $\displaystyle
P^{\prime}\equiv P+Q_{S}\,.$ (13)

In the new valence space $P^{\prime}$, we solve matrix equation with the
energy-dependent effective Hamiltonian Dzuba et al. (1996):

$\displaystyle H_{\mathrm{eff}}(E)=H+V_{\mathrm{eff}}(E)\,,$ (14)
$\displaystyle\hat{P^{\prime}}H_{\mathrm{eff}}(E_{a})\hat{P^{\prime}}\Psi_{a}=E_{a}\hat{P^{\prime}}\Psi_{a}\,,$
(15)

where $\hat{P^{\prime}}$ is the projector on the subspace $P^{\prime}$. When
the size of the matrix $H_{\mathrm{eff}}$ becomes too large, one can neglect
the non-diagonal part of the matrix in the $Q_{S}$ space, as in the emu CI
method Geddes et al. (2018).

##  III Energy denominators

Let us discuss the energy denominator $\Delta_{E}$ in Eq. (11). For simplicity
we will consider the Rayleigh-Schrödinger perturbation theory, where the
denominator in Eq. (6) would be $\bar{E}_{p}-\bar{E}_{q}$. Here $\bar{E}_{p}$
and $\bar{E}_{q}$ are average energies (5) for configurations $p$ and $q$.
Note that in order to return to the Brillouin-Wigner perturbation theory we
will need to add $E_{a}-\bar{E}_{p}$, which can be approximately done using
the method suggested in Dzuba et al. (1996).

In the conventional MBPT the denominator $\bar{E}_{p}-\bar{E}_{q}$ is reduced
to the difference of the Hartree-Fock energies of the orbitals
$\varepsilon_{i}$ which are different in these two configurations. That would
give the following energy denominator in Eq. (11):

$\displaystyle\Delta_{E}\equiv\Delta_{E}(ab\to
mn)=\varepsilon_{a}+\varepsilon_{b}-\varepsilon_{m}-\varepsilon_{n}\,,$ (16)

where we assume that configuration $q$ differs from $p$ by excitation of two
electrons from shells $a$ and $b$ to virtual shells $m$ and $n$ respectively.
This expression neglects the interaction of the electrons with each other and
depends on the choice of the Hartree-Fock potential. In order to improve this
approximation, we will consider general expression for the average energy of
the relativistic electronic configuration.

###  III.1 Average energy of the relativistic configuration

The average energy of the relativistic configuration $\bar{E}_{p}$ Mann
(1973); Grant (1970):

$\bar{E}_{p}=\sum_{a\in p}q_{a}\,I_{a}+\tfrac{1}{2}\,\sum_{a\in
p}q_{a}\,(q_{a}-1)\,U_{aa}\\\ +\sum_{a<b;\,a,b\in p}q_{a}\,q_{b}\,U_{ab},$
(17)

where $q_{a}$ and $q_{b}$ are occupation numbers for the shells $a$ and $b$ in
configuration $p$ and matrix elements of the potential $U$ are given by:

$\displaystyle U_{ab}=\left\\{\begin{array}[]{ll}\displaystyle
F^{0}(a,a)+\sum_{k>0}2\,f^{k}_{a,a}\,F^{k}(a,a)\,,&a=b\,,\\\\[14.22636pt]
\displaystyle F^{0}(a,b)+\sum_{k}g^{k}_{a,b}\,G^{k}(a,b)\,,&a\neq
b\,.\end{array}\right.$ (20)

In these equations $I_{a}$ is the one-electron radial integral, while
$F^{k}(a,b)$ and $G^{k}(a,b)$ are standard Coulomb and exchange two-electron
radial integrals Grant (1970). The angular factors $f^{k}_{a,a}$ and
$g^{k}_{a,b}$ are also defined in agreement with Ref. Grant (1970):

$\displaystyle\begin{array}[]{lll}f^{k}_{a,a}&=\displaystyle-\,\frac{1}{2}\,\frac{2j_{a}+1}{2j_{a}}\,\left(\begin{array}[]{llll}j_{a}&j_{a}&k\\\
\frac{1}{2}&-\frac{1}{2}&0\end{array}\right)^{2}\,,\\\\[14.22636pt]
\displaystyle
g^{k}_{a,b}&=\displaystyle-\left(\begin{array}[]{llll}j_{a}&j_{b}&k\\\
\frac{1}{2}&-\frac{1}{2}&0\end{array}\right)^{2}\,,\end{array}$ (27)

where $j_{a}$ and $j_{b}$ are the one-electron total angular momenta.

Let us use Eq. (17) to calculate the energy difference between configurations
$p$ and $q$ which differ by the excitation of two electrons from shells $a,b$
to shells $m,n$. In other words we need to calculate how the energy changes
when occupation numbers change in the following way: $\delta q_{a}=\delta
q_{b}=-1$ and $\delta q_{m}=\delta q_{n}=1$. To this end we, can use Taylor
expansion of Eq. (17) near the initial configuration $p$:

$\displaystyle\bar{E}_{q}=\bar{E}_{p}+\sum_{a}\frac{\partial\bar{E}_{p}}{\partial
q_{a}}\,\delta
q_{a}+\frac{1}{2}\sum_{a,b}\frac{\partial^{2}\bar{E}_{p}}{\partial
q_{a}\partial q_{b}}\,\delta q_{a}\,\delta q_{b}\,,$ (28)

where derivatives are given by:

$\displaystyle\frac{\partial\bar{E}_{p}}{\partial q_{a}}$
$\displaystyle=I_{a}+(q_{a}-\tfrac{1}{2})U_{aa}+\sum_{b\neq a}q_{b}U_{ab}$
$\displaystyle=I_{a}-\tfrac{1}{2}U_{aa}+\sum_{b}q_{b}U_{ab}\,,$ (29)
$\displaystyle\frac{\partial^{2}\bar{E}_{p}}{\partial q_{a}\partial q_{b}}$
$\displaystyle=U_{ab}\,.$ (30)

Note that all higher derivatives vanish, so expression (28) is exact. With its
help we get:

$\Delta_{E}(ab\to mn)=I_{a}+I_{b}-I_{m}-I_{n}\\\ +\sum_{c\in
p}q_{c}\,(U_{ac}+U_{bc}-U_{mc}-U_{nc})\\\ -U_{aa}-U_{bb}-U_{ab}-U_{mn}\\\
+U_{am}+U_{bn}+U_{an}+U_{bm}\,.$ (31)

This expression can be also used for the special cases $a=b,\,\delta q_{a}=-2$
and/or $m=n,\,\delta q_{m}=2$.

Equation (31) includes the sum over the occupied shells of the initial
configuration $p$. Let us introduce one-electron energies in respect to this
configuration as:

$\displaystyle\varepsilon_{a}$ $\displaystyle=I_{a}+\sum_{c\in
p}q_{c}\,U_{ac}-(1-\delta_{q_{a},0})\,U_{aa}\,.$ (32)

Then Eq. (31) is simplified to

$\Delta_{E}(ab\to
mn)=\varepsilon_{a}+\varepsilon_{b}-\varepsilon_{m}-\varepsilon_{n}\\\
-U_{ab}-U_{mn}+U_{am}+U_{bn}+U_{an}+U_{bm}\,.$ (33)

The first line here reproduces the conventional MBPT denominator (16), while
the second line gives corrections caused by the interactions of the electrons
with each other. It is important that in this form we do not have explicit
sums over all electrons, which significantly simplifies calculations.

In the relativistic calculations the non-relativistic configurations are
typically not used. However, sometimes one may need to find the average energy
of the non-relativistic configuration. In the Appendix A we derive the
necessary expressions for this case.

##  IV Numerical tests

We made four test calculations for very different systems. In the first two
calculations for He I and B I, there was no core and we tested our method for
the valence correlations. Then we applied our method for the highly charged
ion Fe XVII, where there is a very strong central field, correlation
corrections are rather small, and perturbation theory must be quite accurate.
In this system we had core $1s^{2}$, so we calculated core-valence correlation
corrections as well as valence ones. Finally, we made calculations for Sc I,
where valence $3d$ electrons have a large overlap with the core shell $3p^{6}$
and core-valence correlation corrections are as important as valence ones.

Table 1: Ground state binding energy of He I (in a.u.). CI calculations are made for three spaces: $P$, $P+Q_{S}$, and $P+Q$. $\Delta_{P+Q}$ is the difference from the CI result in the $P+Q$ space. Three variants of PT calculations are made based on the CI calculation in $P+Q_{S}$ space: (a) determinant-based PT; (b) effective Hamiltonian with Hartree-Fock denominators (16); (c) effective Hamiltonian with corrected denominators (33). Experimental binding energy is given for comparison in the last column Kramida et al. (2021). | $P$ | $P+Q_{S}$ | $P+Q$ | PT | NIST
---|---|---|---|---|---
|  |  |  | (a) | (b) | (c) | Ref. Kramida et al. (2021)
$E(1s^{2})$ | $2.8626$ | $2.8700$ | $2.9010$ | $2.9021$ | $2.9064$ | $2.9031$ | $2.9034$
$\Delta_{P+Q}$ | $0.0384$ | $0.0310$ | $0.0000$ | $-0.0011$ | $-0.0054$ | $-0.0021$ | $-0.0024$
Table 2: Ground state binding energy of B I (in a.u.). CI calculations are made for valence spaces $P$ and $\tilde{P}$, which included 3 and 4 lower shells respectively. Experimental binding energy is given for comparison in the last column Kramida et al. (2021). | $P$ | $P+Q_{S}$ | $\tilde{P}$ | $\tilde{P}+\tilde{Q}_{S}$ | NIST
---|---|---|---|---|---
|  | $H$ | $H_{\mathrm{eff}}$ |  | $H$ | $H_{\mathrm{eff}}$ | Ref. Kramida et al. (2021)
$E({}^{2}P_{1/2})$ | $24.5683$ | $24.5976$ | $24.6595$ | $24.5721$ | $24.5999$ | $24.6581$ | $24.6581$
$\Delta_{\mathrm{NIST}}$ | $0.0898$ | $0.0605$ | $-0.0014$ | $0.0860$ | $0.0582$ | $0.0000$ | $0.0000$
Table 3: Low-lying energy levels of Fe XVII in respect to the ground state (in cm-1). The subspace $Q_{S}$ includes single excitations to virtual shells $n=5-17$. The subspace $Q_{S}^{\prime}$ in addition includes single excitations from the $1s$ shell. Effective Hamiltonians account for the respective double excitations. For each calculation we also give relative accuracy in percent. Config. | Level | NIST | CI($P$) | CI${}_{\mathrm{emu}}(P+Q_{S})$ | CI${}_{\mathrm{emu}}(P+Q^{\prime}_{S})$
---|---|---|---|---|---
|  | Ref. Kramida et al. (2021) | $H$ | $H_{\mathrm{eff}}$ | $H$ | $H_{\mathrm{eff}}$ | $H^{\prime}_{\mathrm{eff}}$
$2p^{6}$ | ${}^{1}S_{0}$ | $0$ | $0$ |  | $0$ |  | $0$ |  | $0$ |  | $0$ |
$2p^{5}3p$ | ${}^{3}S_{1}$ | $6093450$ | $6076370$ | $-0.28\%$ | $6083540$ | $-0.16\%$ | $6088405$ | $-0.08\%$ | $6095600$ | $0.04\%$ | $6095086$ | $0.03\%$
$2p^{5}3p$ | ${}^{3}D_{2}$ | $6121690$ | $6105049$ | $-0.27\%$ | $6111933$ | $-0.16\%$ | $6117307$ | $-0.07\%$ | $6124215$ | $0.04\%$ | $6123709$ | $0.03\%$
$2p^{5}3p$ | ${}^{3}D_{3}$ | $6134730$ | $6118010$ | $-0.27\%$ | $6125056$ | $-0.16\%$ | $6130067$ | $-0.08\%$ | $6137137$ | $0.04\%$ | $6136602$ | $0.03\%$
$2p^{5}3p$ | ${}^{1}P_{1}$ | $6143850$ | $6127278$ | $-0.27\%$ | $6134193$ | $-0.16\%$ | $6139345$ | $-0.07\%$ | $6146283$ | $0.04\%$ | $6145772$ | $0.03\%$
$2p^{5}3s$ | $2^{o}$ | $5849490$ | $5830778$ | $-0.32\%$ | $5838679$ | $-0.18\%$ | $5842900$ | $-0.11\%$ | $5850823$ | $0.02\%$ | $5850330$ | $0.01\%$
$2p^{5}3s$ | $1^{o}$ | $5864770$ | $5846269$ | $-0.32\%$ | $5854109$ | $-0.18\%$ | $5858397$ | $-0.11\%$ | $5866260$ | $0.03\%$ | $5865678$ | $0.02\%$
$2p^{5}3s$ | $1^{o}$ | $5960870$ | $5942198$ | $-0.31\%$ | $5950103$ | $-0.18\%$ | $5954316$ | $-0.11\%$ | $5962244$ | $0.02\%$ | $5961601$ | $0.01\%$
$2p^{5}3d$ | ${}^{3}P_{1}^{o}$ | $6471800$ | $6455306$ | $-0.25\%$ | $6462010$ | $-0.15\%$ | $6463149$ | $-0.13\%$ | $6469882$ | $-0.03\%$ | $6468962$ | $-0.04\%$
$2p^{5}3d$ | ${}^{3}P_{2}^{o}$ | $6486400$ | $6470075$ | $-0.25\%$ | $6476738$ | $-0.15\%$ | $6477839$ | $-0.13\%$ | $6484531$ | $-0.03\%$ | $6483612$ | $-0.04\%$
$2p^{5}3d$ | ${}^{3}F_{4}^{o}$ | $6486830$ | $6471630$ | $-0.23\%$ | $6478532$ | $-0.13\%$ | $6478129$ | $-0.13\%$ | $6485057$ | $-0.03\%$ | $6484147$ | $-0.04\%$
$2p^{5}3d$ | ${}^{3}F_{3}^{o}$ | $6493030$ | $6477585$ | $-0.24\%$ | $6484338$ | $-0.13\%$ | $6484319$ | $-0.13\%$ | $6491101$ | $-0.03\%$ | $6490177$ | $-0.04\%$
$2p^{5}3d$ | ${}^{1}D_{2}^{o}$ | $6506700$ | $6491383$ | $-0.24\%$ | $6498026$ | $-0.13\%$ | $6498360$ | $-0.13\%$ | $6505032$ | $-0.03\%$ | $6504101$ | $-0.04\%$
Table 4: Low-lying energy levels of Sc I (in cm-1). For each calculation we also give the differences with NIST Kramida et al. (2021) and the average absolute difference $|\Delta|_{\mathrm{av}}=\frac{1}{k}\sum_{i=1}^{k}|\Delta_{i}|$. For the CI calculations in the $P+Q_{S}$ space we use the emu CI approach Geddes et al. (2018) where we neglect non-diagonal matrix elements in the $Q_{S}$ subspace. On the diagonal we use averaging over relativistic configurations, see Eq. (17). Config. | Level | NIST | CI($P$) | CI${}_{\mathrm{emu}}(P+Q_{S})$
---|---|---|---|---
|  | Ref. Kramida et al. (2021) | $H$ | $H$ | $H_{\mathrm{eff}}$
|  | $E$ | $E$ | $\Delta$ | $E$ | $\Delta$ | $E$ | $\Delta$
$3d4s^{2}$ | ${}^{2}D_{3/2}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
| ${}^{2}D_{5/2}$ | $168$ | $147$ | $-21$ | $157$ | $-11$ | $155$ | $-13$
$3d^{2}4s$ | ${}^{4}F_{3/2}$ | $11520$ | $14945$ | $3425$ | $7361$ | $-4159$ | $11786$ | $266$
| ${}^{4}F_{5/2}$ | $11558$ | $14968$ | $3410$ | $7422$ | $-4136$ | $11847$ | $290$
| ${}^{4}F_{7/2}$ | $11610$ | $15001$ | $3391$ | $7489$ | $-4121$ | $11914$ | $304$
| ${}^{4}F_{9/2}$ | $11677$ | $15047$ | $3370$ | $7541$ | $-4136$ | $11963$ | $285$
$3d^{2}4s$ | ${}^{2}F_{5/2}$ | $14926$ | $17368$ | $2442$ | $11331$ | $-3595$ | $15661$ | $735$
| ${}^{2}F_{7/2}$ | $15042$ | $17455$ | $2413$ | $11453$ | $-3589$ | $15781$ | $739$
$3d^{2}4s$ | ${}^{2}D_{5/2}$ | $17013$ | $19972$ | $2960$ | $14574$ | $-2439$ | $17475$ | $462$
| ${}^{2}D_{3/2}$ | $17025$ | $19980$ | $2955$ | $14601$ | $-2424$ | $17500$ | $475$
$3d^{2}4s$ | ${}^{4}P_{1/2}$ | $17226$ | $20329$ | $3103$ | $14606$ | $-2620$ | $17472$ | $246$
| ${}^{4}P_{3/2}$ | $17255$ | $20339$ | $3084$ | $14679$ | $-2576$ | $17552$ | $297$
| ${}^{4}P_{5/2}$ | $17307$ | $\quad 20380$ | $3073$ | $\quad 14739$ | $-2568$ | $\quad 17606$ | $299$
$3d4s4p$ | ${}^{4}F_{3/2}^{o}$ | $15673$ | $13921$ | $-1751$ | $16019$ | $346$ | $15872$ | $200$
| ${}^{4}F_{5/2}^{o}$ | $15757$ | $14002$ | $-1754$ | $16099$ | $342$ | $15953$ | $197$
| ${}^{4}F_{7/2}^{o}$ | $15882$ | $14139$ | $-1743$ | $16211$ | $330$ | $16064$ | $183$
| ${}^{4}F_{9/2}^{o}$ | $16027$ | $14290$ | $-1737$ | $16340$ | $314$ | $16194$ | $168$
$3d4s4p$ | ${}^{4}D_{1/2}^{o}$ | $16010$ | $14265$ | $-1745$ | $16318$ | $308$ | $16448$ | $438$
| ${}^{4}D_{3/2}^{o}$ | $16022$ | $14311$ | $-1711$ | $16351$ | $329$ | $16517$ | $495$
| ${}^{4}D_{5/2}^{o}$ | $16141$ | $14375$ | $-1766$ | $16403$ | $262$ | $16559$ | $418$
| ${}^{4}D_{7/2}^{o}$ | $16211$ | $14458$ | $-1753$ | $16503$ | $292$ | $16621$ | $410$
$3d4s4p$ | ${}^{2}D_{3/2}^{o}$ | $16023$ | $14172$ | $-1851$ | $16516$ | $493$ | $16442$ | $419$
| ${}^{2}D_{5/2}^{o}$ | $16097$ | $14189$ | $-1907$ | $16525$ | $428$ | $16449$ | $352$
$3d4s4p$ | ${}^{4}P_{1/2}^{o}$ | $18504$ | $16854$ | $-1650$ | $18528$ | $24$ | $18529$ | $25$
| ${}^{4}P_{3/2}^{o}$ | $18516$ | $16930$ | $-1586$ | $18538$ | $23$ | $18543$ | $27$
| ${}^{4}P_{5/2}^{o}$ | $18571$ | $17007$ | $-1565$ | $18577$ | $6$ | $18572$ | $1$
$|\Delta|_{\mathrm{av}}$ |  | 2247 | 1595 | 310

###  IV.1 Ground state of He I

Helium is the simplest system where correlation effects can be tested. We
calculate the ground state energy, where correlation corrections are the
largest. We choose the space $P$ to include shells $n=1\dots 3$. The space $Q$
includes virtual shells $s,p,d$ with $4\leq n\leq 20$. For this model problem,
we can easily do CI in the whole space $P+Q$ thus producing the “exact”
solution and compare these results with different variants of the perturbation
theory discussed above. Results are listed in Table 1.

One can see that the valence CI provides accuracy on the order of 1%. The
accuracy does not improve when we account for the single excitations to the
virtual shells. However, when we include double excitations the agreement with
the “exact” answer is significantly better. The determinant-based PT gives the
best result. The results obtained with the effective Hamiltonian are less
accurate, but corrections to the denominators reduce the discrepancy. Even the
uncorrected variant of the MBPT is closer to the “exact” answer by an order of
magnitude compared to the valence CI.

###  IV.2 Ground state of B I

B is a five electron system. The full CI calculation here is already very
expensive. The determinant-based PT is also rather lengthy, so we made
calculations only with the effective Hamiltonian and compared our results with
the experiment Kramida et al. (2021). The effective radial integrals were
calculated using the Hartree-Fock denominators. We tested two variants of the
valence space: the first one, $P$, included shells $n=1\dots 3$ and the second
one, $\tilde{P}$, included also the shell $n=4$. Corresponding $Q$ and
$\tilde{Q}$ spaces included $s,p,d,f,g$ shells up to $n=20$. Results of these
calculations for the ground state ${}^{2}P_{1/2}$ are given in Table 2. We see
that the accuracy of the CI calculation does not change much when we include
an extra shell in the subspace $P$. The accuracy of the CI calculation in the
subspaces $P+Q_{S}$ and $\tilde{P}+\tilde{Q}_{S}$ is only slightly better than
similar calculation in the subspaces $P$ and $\tilde{P}$. Only including
double excitations by means of the MBPT improves the agreement with the
experiment by more than an order of magnitude.

###  IV.3 Spectrum of Fe XVII

Ten-electron ion Fe XVII plays an important role in astrophysics and plasma
physics, see Ref. Kühn et al. (2020a) and references therein. The spectrum of
this ion was calculated within several different approaches Kühn et al.
(2020b) with relative accuracy of about 0.03%. Here we repeat these
calculations using the new method. We use basis set $[17spdfg]$. Virtual
orbitals starting from $4s$ and up are formed from B-splines using the method
from Ref. Kozlov and Tupitsyn (2019). Valence subspace $P$ includes shells
$2s,2p,3s,3p,3d,4s,4p,4d$, and $4f$, while the $1s$ shell is frozen. Single
excitations to all higher orbitals are included in the subspace $Q_{S}$ and
the subspace $Q_{S}^{\prime}$ in addition includes single excitations from the
$1s$ shell. We make two CI calculations in the spaces $P$ and $P+Q_{S}$
respectively. Then we repeat these calculations using the effective
Hamiltonian, which accounts for the excitations to the subspace $Q_{D}$.
Finally, we make CI calculation in the $P+Q_{S}^{\prime}$ for the effective
Hamiltonian $H_{\mathrm{eff}}^{\prime}$ which accounts for the double
excitations from $1s$ shell as well as for the double excitations to the
virtual shells with $n\geq 5$. Results of all these calculations are given in
Table 3.

One can see that already the CI calculation in the subspace $P$ is quite
accurate here, the relative errors being about 0.3%. This is not surprising
for such a strong central field. When we increase the size of the
configuration space by adding single excitations to the virtual shells
$n=5\dots 17$ the errors substantially decrease but remain of the same order
of magnitude. The same happens when we do CI for the effective Hamiltonian in
the subspace $P$. Only when we include both single and double excitations to
the virtual shells by doing CI for the effective Hamiltonian in the subspace
$P+Q_{S}$ we increase the accuracy by an order of magnitude, the errors being
0.04% or less. Adding S and D excitations from the $1s$ shell leads to
corrections to the transition energies within $0.01\%$. Our final accuracy is
similar to the accuracy obtained in Ref. Kühn et al. (2020b), where CI space
included all double and some triple excitations to all virtual shells (the
basis set there was different, but of the same length). In our present
calculation, the size of the space $P+Q_{S}$ is about 1.4 million
determinants, and the size of the space $P+Q_{S}^{\prime}$ is close to 2
million determinants, which is significantly less than the CI space of Ref.
Kühn et al. (2020a).

###  IV.4 Spectrum of Sc I

The ground state configuration for Sc I is $\mathrm{[Ar]}3d^{1}4s^{2}$ and
lowest excited states belong to the configurations $3d^{2}4s$ and $3d4s4p$.
The $3d$ shell has a large overlap with the core shells $3s$ and $3p$. Because
of that frozen core approximation can not reproduce even the lowest part of
the spectrum. Including $3s$ and $3p$ shells into the valence space makes its
size extremely large. Therefore, this is a good system to apply our method.

We use a short basis set $[9spdfgh]$, which is constructed as described in
Ref. Kozlov and Tupitsyn (2019). In the valence space $P$, the shells $n\leq
3$ are closed and the virtual shells $n\geq 8$ and all $h$ orbitals are empty.
The space $Q_{S}$ includes single excitations from the upper core shells $n=3$
and single excitations to the virtual shells. We keep core shells up to $n\leq
2$ frozen on all stages. Results of the calculation of the spectrum are
presented in Table 4, where excitation energies from the ground state in cm -1
are shown for approximately 10 lower levels of each parity. The sizes of the
valence space $P$ and $P+Q_{S}$ are about $6\times 10^{4}$ and $1\times
10^{6}$ determinants respectively. We list the results of three calculations:
the full CI in the valence space $P$ and emu CI Geddes et al. (2018) in the
space $P+Q_{S}$ for the bare and the effective Hamiltonians. The effective
radial integrals were calculated with the Hartree-Fock denominators. For each
of these calculations we also give differences from the experimental values
Kramida et al. (2021) and the averaged absolute difference.

One can see that all the levels in the CI calculation are shifted from their
experimental energies: the levels of the configuration $3d^{2}4s$ lie higher
by 3 thousand inverse centimeters, while the levels of the configuration
$3d4s4p$ lie lower by 2 thousand inverse centimeters. The picture changes
drastically when we add single excitations and solve the problem in the space
$P+Q_{S}$. Now the levels of the configuration $3d^{2}4s$ lie lower by 3
thousand inverse centimeters, while the levels of the configuration $3d4s4p$
are almost in place. Finally, when we use the effective Hamiltonian, which
accounts for the double excitations, the levels get closer to their places
with the average deviation about 300 cm-1, or 7 times smaller, than for the CI
calculation.

In this test calculation, we used a rather short basis set and were probably
rather far from saturation. Therefore we can not reliably estimate the
ultimate accuracy of the method for scandium. Looking at the results we see
that the size of the PT corrections is very large and there is also large
cancellation between contributions of the single and double excitations.
Therefore it is unlikely that converged results would be significantly better
than what we got here. On the other hand, we see systematic improvement in our
final results compared to the pure valence calculation. It is also worth
mentioning that if one would try to include all double excitations in CI
calculation the size of the configuration space would be much above $1\times
10^{8}$ even for the basis set as short as this one.

##  V Conclusions

We suggest a new version of the CI+MBPT method Dzuba et al. (1996) with the
different division of the many-electron space into parts where non-
perturbative and perturbative methods are used. This new division may be more
practical for the atoms with many valence electrons, where the size of the
valence space may be too big for solving the matrix eigenvalue problem. This
method can be used in the all-electron calculations for light atoms as well as
for the calculations with the frozen core. In the latter case, the single and
double excitations from (some of) the core shells can be treated
perturbatively. We ran four rather different tests which showed systematic
one-order-of-magnitude improvement of the results when we added MBPT
corrections to the CI calculations.

## Acknowledgements

We thank Marianna Safronova, Charles Cheung, and Sergey Porsev for their
constant interest in this work and very useful discussions. This work was
supported by the Russian Science Foundation (Grant No. 19-12-00157). I.I.T.
acknowledges the support from the Resource Center “Computer Center of SPbU”,
St. Petersburg, Russia.

##  Appendix A Average energy of the non-relativistic configuration

In the average over non-relativistic configuration (LS-average) Tupitsyn et
al. (2018); Lindgren and Rosén (1975), the occupation numbers for the
relativistic orbitals $q_{a}$ may be non-integer, while occupation numbers for
non-relativistic orbitals $q_{A}$ are still integer (we use capital letters
$A,B,M,N$ to designate non-relativistic orbitals). Below we show that properly
defining one-electron integrals $I_{A}$ and two-electron matrix elements
$U_{AB}$ we obtain expressions similar to Eqs. (32, 33).

The average energy of the non-relativistic configuration $R$ can be written
as:

$\bar{E}_{R}=\sum_{a}\tilde{q}_{a}\,I_{a}+\frac{1}{2}\sum_{a}\tilde{q}_{a}\,(\tilde{q}_{a}-w_{a})\,F^{0}(a,a)+\sum_{a<b}\tilde{q}_{a}\,\tilde{q}_{b}\,w_{ab}\,F^{0}(a,b)\\\
+\sum_{a,k>0}\tilde{q}_{a}(\tilde{q}_{a}-w_{a})\,f^{k}_{aa}\,F^{k}(a,a)+\sum_{a<b,k}\tilde{q}_{a}\,\tilde{q}_{b}\,w_{ab}\,g^{k}_{ab}\,G^{k}(a,b)\,,$
(34)

where

$\displaystyle\tilde{q}_{a}=\displaystyle\frac{2j_{a}+1}{4l_{a}+2}\,q_{A}\,,\qquad
w_{a}=\displaystyle\frac{q_{A}-\tilde{q}_{a}+2j_{a}}{4l_{a}+1}\,,\qquad
w_{ab}=\left\\{\begin{array}[]{cll}\frac{4l_{a}+2}{4l_{a}+1}\,\frac{q_{A}-1}{q_{A}}\,,&\,A=B\,,\\\
1\,,&\,A\neq B\,.\end{array}\right.$ (37)

Using expressions:

$\displaystyle\tilde{q}_{a}(\tilde{q}_{a}-w_{a})=\frac{2j_{a}+1}{4l_{a}+2}\,\frac{2j_{a}}{4l_{a}+1}\,q_{A}(q_{A}-1)\,,$
$\displaystyle
w_{ab}\,\tilde{q}_{a}\tilde{q}_{b}=\frac{2j_{a}+1}{4l_{a}+2}\,\frac{2j_{a^{\prime}}+1}{4l_{a}+1}\,q_{A}(q_{A}-1)\,,\quad
A=B,\,\,j_{a}\neq j_{a^{\prime}}\,,$ (38)

we rewrite the equation (34) in the form

$\bar{E}_{R}=\sum_{A}q_{A}\,\sum_{j_{a}}\frac{2j_{a}+1}{4l_{a}+2}I_{a}+\frac{1}{2}\sum_{A}q_{A}\,(q_{A}-1)\,\sum_{a,a^{\prime}\in
A}\frac{(2j_{a}+1)(2j_{a^{\prime}}+1-\delta_{a,a^{\prime}})}{(4l_{a}+2)(4l_{a}+1)}\,F^{0}(a,a^{\prime})+\\\
+\frac{1}{2}\sum_{A\neq B}q_{A}\,q_{B}\sum_{a\in A,b\in
B}\frac{(2j_{a}+1)(2j_{b}+1)}{(4l_{a}+2)(4l_{b}+2)}\,\left[F^{0}(a,b)+\sum_{k}g^{k}_{ab}\,G^{k}(a,b)\right]\\\
+\frac{1}{2}\sum_{A}\,q_{A}(q_{A}-1)\sum_{a,a^{\prime}\in
A}\sum_{k>0}\,\frac{2j_{a}+1}{4l_{a}+2}\,\frac{2j_{a^{\prime}}+1}{4l_{a}+1}\,g^{k}_{aa^{\prime}}\,G^{k}(a,a^{\prime})\,.$
(39)

In the last sum, the term $k=0$ is absent since $j_{a}\neq j_{a^{\prime}}$ and
$k\geq|j_{a}-j_{a^{\prime}}|$. Now we can introduce non-relativistic analogues
of the integrals $I_{a}$ and matrix elements $U_{ab}$ and rewrite Eq. (34)
alike Eq. (17):

$\displaystyle\bar{E}_{R}$
$\displaystyle=\sum_{A}q_{A}\,I_{A}+\frac{1}{2}\,\sum_{A}q_{A}\,(q_{A}-1)\,U_{AA}+\frac{1}{2}\,\sum_{A\neq
B}\,q_{A}\,q_{B}\,U_{AB}\,,$ (40) $\displaystyle I_{A}$
$\displaystyle=\sum_{a\in A}\frac{2j_{a}+1}{4l_{a}+2}\,I_{a}\,,$ (41)
$\displaystyle U_{AA}$ $\displaystyle=\sum_{a,a^{\prime}\in
A}\frac{(2j_{a}+1)(2j_{a^{\prime}}+1)}{(4l_{a}+2)(4l_{a}+1)}\,\left[F^{0}(a,a^{\prime})+\sum_{k>0}g^{k}_{aa^{\prime}}\,G^{k}(a,a^{\prime})\right]-\sum_{a\in
A}\frac{(2j_{a}+1)}{(4l_{a}+2)(4l_{a}+1)}\,F^{0}(a,a)\,,$ (42) $\displaystyle
U_{AB}$ $\displaystyle=\sum_{a\in A,b\in
B}\frac{(2j_{a}+1)(2j_{b}+1)}{(4l_{a}+2)(4l_{b}+2)}\,\left[F^{0}(a,b)+\sum_{k}g^{k}_{ab}\,G^{k}(a,b)\right]\,.$
(43)

Using Eq. (40) we get the following derivatives by analogy with Eqs. (29, 30):

$\displaystyle\frac{\partial\bar{E}_{R}}{\partial
q_{A}}=I_{A}+\sum_{B}q_{B}\,U_{AB}-\frac{1}{2}\,U_{AA}\,,$
$\displaystyle\frac{\partial^{2}\bar{E}_{R}}{\partial q_{A}\partial q_{B}}$
$\displaystyle=\displaystyle U_{AB}\,.$ (44)

The difference in energy between two configurations is

$\displaystyle\Delta\bar{E}$ $\displaystyle=\sum_{A}I_{A}\,\delta
q_{A}+\sum_{A}\left(q_{A}-\frac{1}{2}\right)\,U_{AA}\,\delta q_{A}+\sum_{B\neq
A}q_{B}\,U_{AB}\,\delta q_{B}+\frac{1}{2}\sum_{A,B}U_{AB}\,\delta
q_{A}\,\delta q_{B}\,.$ (45)

This equation allow us to find the energy of a double excitation $\delta
q_{A}=-1,\,\delta q_{B}=-1,\,\delta q_{N}=1,\,\delta q_{M}=1$:

$\Delta_{\bar{E}}(AB\to
NM)=I_{A}+I_{B}-I_{M}-I_{N}+\sum_{C}q_{C}\,(U_{AC}+U_{BC}-U_{MC}-U_{NC})\\\
-U_{AA}-U_{BB}-U_{AB}-U_{NM}+U_{AN}+U_{BN}+U_{AM}+U_{BM}\,.$ (46)

If we introduce an averaged one-electron energy by analogy with (32) we can
rewrite (46) as

$\displaystyle\varepsilon_{A}$
$\displaystyle=I_{A}+\sum_{B}q_{B}\,U_{AB}-(1-\delta_{q_{A},0})\,U_{AA}\,,$
(47) $\displaystyle\Delta_{\bar{E}}(AB\to NM)$
$\displaystyle=\varepsilon_{A}+\varepsilon_{B}-\varepsilon_{M}-\varepsilon_{N}-U_{AB}-U_{NM}+U_{AN}+U_{BN}+U_{AM}+U_{BM}\,.$
(48)

We obtained corrections to the standard MBPT energy denominator (16) using two
approximations. Averaging over relativistic configurations gives expression
(33) and averaging over non-relativistic configurations leads to expression
(48). These expressions differ only by the definitions of the one-electron
energies and two-electron matrix elements.

## References

  * Desclaux (1975) J. P. Desclaux, Computer Physics Communications 9, 31 (1975).
  * Jönsson et al. (2013) P. Jönsson, G. Gaigalas, J. Bieroń, C. F. Fischer, and I. P. Grant, 184, 2197 (2013).
  * Froese Fischer et al. (2016) C. Froese Fischer, M. Godefroid, T. Brage, P. Jönsson, and G. Gaigalas, Journal of Physics B Atomic Molecular Physics 49, 182004 (2016).
  * Gu (2008) M. F. Gu, Canadian Journal of Physics 86, 675 (2008).
  * Fritzsche et al. (2002) S. Fritzsche, C. F. Fischer, and G. Gaigalas, Computer Physics Communications 148, 103 (2002).
  * Jiang et al. (2016) J. Jiang, J. Mitroy, Y. Cheng, and M. W. J. Bromley, Phys. Rev. A 94, 062514 (2016).
  * Tupitsyn et al. (2003) I. I. Tupitsyn, V. M. Shabaev, J. R. Crespo López-Urrutia, I. Draganić, R. S. Orts, and J. Ullrich, Phys. Rev. A 68, 022511 (2003).
  * Tupitsyn et al. (2005) I. I. Tupitsyn, A. V. Volotka, D. A. Glazov, V. M. Shabaev, G. Plunien, J. R. C. López-Urrutia, A. Lapierre, and J. Ullrich, Phys. Rev. A 72, 062503 (2005).
  * Dzuba et al. (1987) V. A. Dzuba, V. V. Flambaum, P. G. Silvestrov, and O. P. Sushkov, J. Phys. B 20, 1399 (1987).
  * Blundell et al. (1987) S. A. Blundell, D. S. Guo, W. R. Johnson, and J. Sapirstein, Atomic Data and Nuclear Data Tables 37, 103 (1987).
  * Sapirstein (1998) J. Sapirstein, Reviews of Modern Physics 70, 55 (1998).
  * Dzuba et al. (1996) V. A. Dzuba, V. V. Flambaum, and M. G. Kozlov, Phys. Rev. A 54, 3948 (1996).
  * Kahl and Berengut (2019) E. Kahl and J. Berengut, Computer Physics Communications 238, 232 (2019), eprint arXiv:1805.11265.
  * Kozlov et al. (2015) M. G. Kozlov, S. G. Porsev, M. S. Safronova, and I. I. Tupitsyn, Comput. Phys. Commun. 195, 199 (2015), ISSN 0010-4655.
  * Eliav and Kaldor (2010) E. Eliav and U. Kaldor, _Relativistic Four-Component Multireference Coupled Cluster Methods: Towards A Covariant Approach_ (Springer Netherlands, Dordrecht, 2010), pp. 113–144, ISBN 978-90-481-2885-3.
  * Saue et al. (2020) T. Saue, R. Bast, A. S. P. Gomes, H. J. A. Jensen, L. Visscher, I. A. Aucar, R. Di Remigio, K. G. Dyall, E. Eliav, E. Fasshauer, et al., J. Chem. Phys. 152, 204104 (2020), eprint arXiv:2002.06121.
  * Blundell et al. (1989) S. A. Blundell, W. R. Johnson, Z. W. Liu, and J. Sapirstein, Phys. Rev. A 40, 2233 (1989).
  * Safronova and Johnson (2008) M. S. Safronova and W. R. Johnson, Advances in Atomic Molecular and Optical Physics 55, 191 (2008).
  * Oleynichenko et al. (2020) A. V. Oleynichenko, A. Zaitsevskii, L. V. Skripnikov, and E. Eliav, Symmetry 12 (2020).
  * Johnson (2007) W. R. Johnson, _Atomic Structure Theory. Lectures on Atomic Physics_ (Springer Berlin Heidelberg, 2007).
  * Flambaum and Ginges (2005) V. V. Flambaum and J. S. M. Ginges, Phys. Rev. A 72, 052115 (2005).
  * Ginges and Berengut (2016) J. S. M. Ginges and J. C. Berengut, Phys. Rev. A 93, 052509 (2016).
  * Shabaev et al. (2013) V. M. Shabaev, I. I. Tupitsyn, and V. A. Yerokhin, Phys. Rev. A 88, 012513 (2013).
  * Tupitsyn et al. (2016) I. I. Tupitsyn, M. G. Kozlov, M. S. Safronova, V. M. Shabaev, and V. A. Dzuba, Phys. Rev. Lett. 117, 253001 (2016).
  * Dzuba et al. (2017) V. A. Dzuba, J. C. Berengut, C. Harabati, and V. V. Flambaum, Phys. Rev. A 95, 012503 (2017).
  * Geddes et al. (2018) A. J. Geddes, D. A. Czapski, E. V. Kahl, and J. C. Berengut, Physical Review A 98, 042508 (2018), eprint arXiv:1805.06615.
  * Dzuba et al. (2019) V. A. Dzuba, V. V. Flambaum, and M. G. Kozlov, Phys. Rev. A 99, 032501 (2019), eprint ArXiv:1812.09826.
  * Imanbaeva and Kozlov (2019) R. T. Imanbaeva and M. G. Kozlov, Annalen der Physik 531, 1800253 (2019).
  * Cheung et al. (2020) C. Cheung, M. S. Safronova, S. G. Porsev, M. G. Kozlov, I. I. Tupitsyn, and A. I. Bondarev, Physical Review Letters 124, 163001 (2020), eprint arXiv:1912.08714.
  * Li and Dzuba (2020) J. Li and V. Dzuba, Journal of Quantitative Spectroscopy and Radiative Transfer 247, 106943 (2020), eprint arXiv:2002.08033.
  * Mann (1973) J. B. Mann, Atomic Data and Nuclear Data Tables 12, 1 (1973).
  * Grant (1970) I. P. Grant, Advances in Physics 19, 747 (1970).
  * Kramida et al. (2021) A. Kramida, Y. Ralchenko, J. Reader, and NIST ASD Team, _Nist atomic spectra database (version 5.9)_ (2021), URL https://physics.nist.gov/asd.
  * Kühn et al. (2020a) S. Kühn, C. Shah, J. R. C. López-Urrutia, K. Fujii, R. Steinbrügge, J. Stierhof, M. Togawa, Z. Harman, N. S. Oreshkina, C. Cheung, et al., Physical Review Letters 124, 225001 (2020a), eprint arXiv:1911.09707.
  * Kühn et al. (2020b) S. Kühn et al., _High resolution photoexcitation measurements exacerbate the long-standing Fe XVII oscillator strength problem: Supplemental material_ (2020b), URL 10.1103/physrevlett.124.225001.
  * Kozlov and Tupitsyn (2019) M. Kozlov and I. Tupitsyn, Atoms 7, 92 (2019).
  * Tupitsyn et al. (2018) I. I. Tupitsyn, N. A. Zubova, V. M. Shabaev, G. Plunien, and T. Stöhlker, Physical Review A 98, 022517 (2018), eprint arXiv:1804.09986.
  * Lindgren and Rosén (1975) I. Lindgren and A. Rosén, in _Case Studies in Atomic Physics_ , edited by E. McDaniel and M. McDowell (Elsevier, 1975), pp. 93 – 196, ISBN 978-0-7204-0331-2.


# Computational Approaches for Modeling Power Consumption on an
Underwater Flapping Fin Propulsion System

Brian Zhou1 , Jason Geder2, Alisha Sharma2,3, Julian Lee4, Marius Pruessner2,
Ravi Ramamurti2, Kamal Viswanath2 This work was performed at the U.S. Naval
Research Laboratory through a base program funded effort as part of the
Science and Engineering Apprenticeship (SEAP) Program.

###### Abstract

The last few decades have led to the rise of research focused on propulsion
and control systems for bio-inspired unmanned underwater vehicles (UUVs),
which provide more maneuverable alternatives to traditional UUVs in underwater
missions. Propulsive efficiency is of utmost importance for flapping-fin UUVs
in order to extend their range and endurance for essential operations. To
optimize for different gait performance metrics, we develop a non-dimensional
figure of merit (FOM), derived from measures of propulsive efficiency, that is
able to evaluate different fin designs and kinematics, and allow for
comparison with other bio-inspired platforms. We create and train
computational models using experimental data, and use these models to predict
thrust and power under different fin operating states, providing efficiency
profiles. We then use the developed FOM to analyze optimal gaits and compare
the performance between different fin materials. These comparisons provide a
better understanding of how fin materials affect our thrust generation and
propulsive efficiency, allowing us to inform control systems and weight for
efficiency on an inverse gait-selector model.

## Introduction

Unmanned underwater vehicles (UUVs) have a variety of industrial and research
applications including exploration, mapping, and minesweeping. Historically,
these operations have primarily been conducted with propeller-driven UUVs.
Because propeller-driven UUVs lack maneuverability and subsequently resistance
to turbulence, their operational domain is limited to relatively quiescent and
deeper waters. Marine animals offer promising solutions to expanding the
envelope of UUV operations because they swim with high propulsive efficiency
and have high maneuverability in water (Masud, La Mantia, and Dabnichki 2022;
Eloy 2012; Taylor, Nudds, and Thomas 2003; Rohr and Fish 2004; Rohr et al.
1998; Triantafyllou, Triantafyllou, and Yue 2000; Nedelcu et al. 2018), which
motivates replication of their fins and other appendages in robotic designs
(Techet 2008). As such, recent decades have seen the rise of research in bio-
inspired propulsion systems to fill the operational gap in littoral waters and
create systems with greater agility and maneuverability compared to
traditional propeller-based systems. Fin designs inspired by a variety of
animals have been studied, including dolphins (Rohr et al. 1998; Rohr and Fish
2004), penguins (Masud, La Mantia, and Dabnichki 2022), and snakes. Among
these, fish-inspired fins have driven the majority of research due to the
agility these species exhibit, outperforming capabilities of traditional
propulsion-based systems (Tangorra et al. 2007; Lauder and Madden 2006;
Mignano et al. 2019; Nedelcu et al. 2018).

Designing, optimizing, and replicating marine flapping fin motion with bio-
inspired fins requires testing of different parameters. Previous studies have
examined and tested the effects that different parameters such as material
properties, kinematics or fin gaits, and fin shape have on a flapping fin’s
thrust output to better replicate and understand fish hydrodynamic performance
(Sampath et al. 2020; Geder et al. 2013; Mignano et al. 2019; Yun, Kim, and
Kim 2015; Nguyen, Lee, and Ahn 2016; Geder et al. 2017). We identify and aim
to resolve the following gaps in literature:

  1. 1.

Prior studies have placed emphasis on studying the effects of designs and
gaits have on the thrust and lift forces (Sampath et al. 2020; Geder et al.
2013; Mignano et al. 2019; Yun, Kim, and Kim 2015; Nguyen, Lee, and Ahn 2016;
Geder et al. 2017). Less research has been conducted on the effects designs
and flapping gaits have on power consumption and UUV efficiency.

  2. 2.

Previous research on power has utilized CFD simulations to study hydrodynamic
power (Palmisano et al. 2013), but this is ineffective for control system
integration that requires taking into account power loss of integrated
actuators.

Analyzing practical power draw and building models will allow us to create
better design and gait recommendations while considering propulsive efficiency
and integrating power onto a model-driven control system. This allows a
control system to have different gait settings, such as a toggle-able setting
optimizing for force output or gait efficiency, drastically increasing the
possible mission duration for size-constrained littoral UUVs.

We explore the prediction of actuator power consumption and integration of
these predictions into a control system for a set of fins on a UUV. Propeller-
based systems have often been preferred due to both their repair-ability and
efficiency (Palmisano et al. 2013). While flapping fin systems may not exceed
propellers in propulsive efficiency, it is essential for a prototype UUV to
optimize for efficiency in thrust generation to ensure feasibility in a range
of naval missions. Optimizing for efficiency is complex; previous research
demonstrates that multiple variables affect thrust generation and efficiency
in tandem fin configurations including flow speed and fin phase (Mignano et
al. 2019), as well as other stroke kinematics. Due to the popularity of bio-
inspired fish fins, there have been many propulsors already created (Tangorra
et al. 2007), but few have developed the models necessary to use the fins on
platforms or offer a dimensionless metric that can be used to compare the
efficiency of different designs between different UUV control systems. We
generate a dimensionless figure of merit (FOM) that can be used for creating a
comparative to other designs and control systems and optimizing power
efficiency.

To generate a FOM, we develop forward models with mathematical and machine
learning approaches that use static and dynamic information about the fin
design and kinematics on the control system to output a low-power, time-
efficient, and accurate prediction of power consumption. Combined with another
integrated model for determining thrust (Lee et al. 2022), the model can take
in a desired thrust or location to find the most power-efficient fin
kinematics to accomplish its goal.

We use the validated thrust and power model to output kinematics
interpolations across the entire possible gait space, allowing us to analyze
FOM trends and develop a forward FOM model that can predict the efficiency of
a certain gait. We generate a FOM that accomplishes the following objectives:

  1. 1.

Adaptability: the FOM forward model can function even if physical properties
(actuator, fin design) are switched by simply retraining the model

  2. 2.

Practicality: the FOM can accurately model and account for an actuator’s power
loss

  3. 3.

Flexibility: the FOM can optimize for different outputs or forces by simply
feeding in a different model

## Materials and Methods

The research objectives guiding this study are to identify the most important
parameters that affect flapping fin power consumption and efficiency to
develop a forward model of propulsive efficiency for practical comparison and
analysis.

### Parameters and Outputs

There is one fin design parameter in this study, material flexibility, and
there are three fin materials in total (Table 1). There are also four
kinematics parameters that define a singular gait, including stroke phase
offset, stroke amplitude, pitch amplitude, and flapping frequency (Table 2).

We collect experimental data on forces, current draw and voltage for the
different fin gaits, which we use to develop our forward models for force and
power. This information is laid out in Table 3, and force vectors are defined
in Figure 1.

Table 1: Fin Material Properties Material |  Young’s Modulus
---|---
Rigid Nylon |  1 GPa
PDMS 1:10 |  850 kPa
PDMS 1:20 |  310 kPa
Table 2: Kinematic Gait Parameters Parameter |  Symbol |  Description
---|---|---
Static Kinematics
Frequency (Hz) |  $f$ |  Number of flap cycles per second
Stroke Pitch Offset ($\degree$) |  $\delta$ |  Phase offset of the pitch cycle to the stroke cycle, calculated as $\frac{1}{16}$th of one cycle
Stroke Amplitude ($\degree$) |  $\Phi$ |  Maximum stroke angle over one flap cycle
Pitch Amplitude ($\degree$) |  $\Theta$ |  Maximum pitch angle over one flap cycle
Dynamic Kinematics
Stroke Angle ($\degree$) |  $\phi$ |  Time history of stroke angle
Pitch Angle ($\degree$) |  $\theta$ |  Time history of pitch angle

Figure 1: CAD design of single flapping fin propulsor with coordinate reference frames. Table 3: Control System Measurements Parameter |  Symbol |  Description
---|---|---
XYZ Forces
Thrust (N) |  $T$ |  Force generated along stroke axis
Lift (N) |  $L$ |  Force generated perpendicular to both
Side Force (N) |  $S$ |  Force generated along pitch axis
Power Consumption
Stroke Current (A) |  $I_{\phi}$ |  Time history of current draw for the stroke actuator
Pitch Current (A) |  $I_{\theta}$ |  Time history of current draw for the pitch actuator
Voltage ($V$) |  $V$ |  Voltage of both actuators

### Data Collection

The parameter space for the kinematics variables is large, but we constrain
the tests to parameter values that are physically achievable, given an
operating frequency. At higher frequencies, the fins are physically unable to
reach certain stroke and pitch amplitudes. Equation 1 defines our range of
achievable strokes and pitches with respect to frequency:

$\text{Attainable gaits:}\begin{cases}\text{$0<\Phi<97-f*30$}\\\
\text{$0<\Theta<75-f*26$}\\\ \end{cases}$ (1)

With this equation, we collected data that covers the entire scope of the
achievable kinematics range at each frequency. In total, 864 unique gaits
listed in Table 4 were tested for each design. For each fin gait, 10 flap
cycles were run, and 5 of the middle cycles were used for analysis to account
for discrepancies when the actuator started and ended the cycle motions. Data
was collected in a zero velocity flow condition, which previous research has
demonstrated is important for low-speed maneuvering to station-keep and offset
buoyancy (Geder et al. 2021).

Table 4: Gait Combinations Parameter |  Symbol
---|---
Stroke Amplitude ($\degree$) |  $0,15,25,32.5,40,55$
Pitch Amplitude ($\degree$) |  $0,15,25,32,38,55$
Frequency (Hz) |  $0.75,1.00,1.25,1.50,1.75,2.00$
Stroke Pitch Offset ($\degree$) |  $-22.5,0,22.5,45$
Voltage ($V$) |  Constant at $4.98V$

This process was replicated for all three fin material designs. Taking the
experimental data as defined in Table 3, we compute the total power
consumption of both the stroke and pitch actuators with Equation 2.

$P=I_{\phi}*V+I_{\theta}*V$ (2)

### Experimental Setup

Data was collected on the experimental setup shown in Figures 2 and 3. As
prior research has demonstrated that tandem fin configurations have minimal
effects on thrust output and power consumption, we only run single-fin tests
to gather data.

Figure 2: CAD design of tandem fins mounted to instrumentation and control
platform.

Figure 3: Tandem fin platform in experimental test environment.

The control platform was mounted in a 2.41 x 0.76 x 0.76m (length, width,
height) glass tank. A microcontroller controlled the fin actuators to collect
data on the programmed gait combinations as laid out in Tables 2 and 4.
Potentiometers (TT Electronics P260) measure the stroke and pitch angles over
time, while load cells (Interface 3A60A) measure the generated forces (Geder
et al. 2021).

## Model

### Objectives

Developing a forward model that inputs gait information and outputs a thrust
or power prediction allows for a wider range that we can interpolate results
from our figure of merit. Additionally, a fast enough model can allow for
future integration onto the control system to take in power or the FOM as a
metric in an inverse model. These objectives guided the selection and
development of our model:

  1. 1.

Completion of a baseline model that can accurately take in kinematics data
(frequency, stroke amplitude, pitch amplitude, and offset) to output predicted
power

  2. 2.

Capability to take in static information such as material and flexibility to
use different models to maximize accuracy and usefulness of integrated model

  3. 3.

Run-time speed of 100 forward passes per second at minimum

First, the model takes in design-related parameters such as material, design,
rigidity, and tandem fin spacing to procure the most optimal model. Since the
fins installed are static and will not be replaced during missions, this
allows the model to only load what is relevant to the mission without wasting
excess computational power. Second, the loaded model will take in more dynamic
and changing kinematics-related information such as frequency, stroke angle,
pitch angle, angle offset, and tandem phasing.

### Model Results

We examined five approaches to model power; two were polynomial models and
three utilized ML approaches.

Our baseline model, the linear model, performed the fastest with an average
error across all 3 data sets of 0.3891 W. As the true value and predicted
values appeared to have an exponential relationship in the linear model, we
tested a quartic polynomial which produced better results (Figure 4) with an
average error of 0.1815 W. Out of the three data sets, the PDMS 1:10 fins fit
the best to the quartic regression.

Figure 4: Quartic polynomial model performance for synthetic data, comparing
predicted power values to the true power value. Colors indicate the material,
with blue (rigid) green (1:10) red (1:20)

Utilizing a Convolutional Neural Network (CNN) as a baseline ML approach, we
see an increase in accuracy with an average error of only 0.0907 W while still
being able to operate at a speed fast enough for model integration.

While the first three models were implemented and other approaches such as a
Multi Layer Perceptron Regression were considered, we ultimately chose to use
a time-series model for modeling thrust and power to understand patterns
present throughout a gait. As such, we used a Long Short-Term Memory model
(LSTM) similar to the thrust LSTM from previous literature (Lee et al. 2022)
that is able to substantially decrease the average error between interpolated
data and create a time series of power consumption during one flap cycle. Most
importantly for our purposes, the LSTM can accurately interpolate between gaps
of data, which is useful for our FOM analysis.

When trained on all experimental gaits, the thrust LSTM reached an average
error of 0.0083 N and the power LSTM reached an average error of 0.0408 W. The
most visible shortfall of the power LSTM model is its inability to grasp the
time history of power consumption at certain gaits. At both extremes of gaits
with lower amplitudes or higher amplitudes, the power history of each cycle
can be vastly different, leaving the model unable to fully understand the time
history. For the purpose of interpolating the power consumption across a gait,
this is not a setback.

Out of all five power models, the LSTM has the best performance but is also
the most time-consuming. Table 5 highlights the results and relative time
performances for a forward model.

Table 5: Model Performances Model |  Averaged Error (W) |  Runtime
---|---|---
Linear Polynomial |  0.3891 |  Low
Quartic Polynomial |  0.1815 |  Low
MLP |  0.1229 |  Medium
CNN |  0.0907 |  Medium
LSTM |  0.0072 |  High

### Generating Models

With the LSTM having the most accurate interpolations, we train each LSTM
model to 1000 epochs for modeling both power and thrust predictions. For the
thrust model, we utilize a previously developed model on the rigid fin data
set (Lee et al. 2022) and tune it to produce optimal results on the PDMS 1:10
and 1:20 data sets. The statistics for the LSTM models produced are found in
Table 6.

Table 6: Average Error for LSTM Models Interpolations |  Rigid |  PDMS 1:10 |  PDMS 1:20
---|---|---|---
Power (W) |  0.0236 |  0.072 |  0.0268
Thrust (N) |  0.0002 |  0.0186 |  0.0041

## Figure of Merit

### Objectives

We create a dimensionless FOM that takes in a gait or flap cycle along with
relevant design information and outputs a metric that rates the efficiency of
force generation.

### Development

A baseline FOM simply computing the force in question over the power would
provide a relative metric that can be used to compare gaits assuming the same
flow speed, displayed in Equation 3. For integration purposes, this would be
suitable and could return as the pure value or a percentage compared to the
highest recorded gait efficiency for a loaded design.

$\eta=\frac{F_{avg}}{P_{avg}}$ (3)

While this FOM has units, we can convert it to a dimensionless parameter by
multiplying by a velocity. Since all tests were conducted in a constant,
static flow, this FOM velocity term can be set to 1 m/s, as shown in Equation
4, effectively cancelling the velocity term for relative efficiency
comparisons between our tests.

$\eta=\frac{F_{avg}*1\text{ m/s}}{P_{avg}}$ (4)

To make our FOM dimensionless for all future flow conditions, we can iterate
off a previous FOM to develop a universal FOM that implements the flow speed
by calculating the tip velocity, as shown in Equations 5 and 6.

$\eta=\frac{F_{avg}v}{P_{avg}}$ (5) $v=\frac{F_{avg}}{P_{avg}}$ (6)

This FOM is easily able to adapt to different objectives. Because we collect
data for all force axes and the individual stroke and pitch actuators, we are
able to create new deviations of Equation 5 to compare an individual fin
efficiency or how efficient a gait is for a specific force axis (Equation 7).
The FOM can also be used to isolate how efficiency changes when changing
individual actuators (Equation 8), and optimize for different objectives in
future integration. Examples are listed below:
$T$ = Thrust, $L$ = Lift, $S$ = Side Force

$\eta_{Ft}=\frac{T_{avg}v}{P_{avg}};\quad\eta_{Fl}=\frac{L_{avg}v}{P_{avg}};\quad\eta_{Fs}=\frac{S_{avg}v}{P_{avg}}$
(7)

$P_{s}$ = Power consumed by stroke actuator
$P_{p}$ = Power consumed by pitch actuator

$\eta_{As}=\frac{F_{avg}v}{P_{s}};\quad\eta_{Ap}=\frac{F_{avg}v}{P_{p}}$ (8)

## Results

We used our generated LSTM models to generate predictions at points
interpolated between training data. In total, we generated interpolations
within the constraints given by our collected data. We interpolated data for
every stroke and pitch combination from 0 to 55 degrees with 1 degree
increments, frequency from 0.75Hz to 2Hz with 0.125Hz increments, and SPO from
-22.5 to 45 degrees with 5.625 degree increments. In total, 435,600
interpolations were calculated for each data set. A sample of the rigid fin
data space is shown in Figure 5. The interpolations filled gaps of data.

Figure 5: 3D visualization of the rigid fin data set. Figure 6: Figure of
merit, thrust, and power contours. The columns are different data sets (rigid,
PDMS 1:10, and PDMS 1:20) while the rows graph different gait results (FOM
value, Thrust, and Power). The PDMS 1:10 fin generates the highest possible
thrust and is higher overall in more gaits; the rigid fin is significantly
worse at thrust generation than either of the PDMS fins. The PDMS 1:10 and
1:20 fins are comparable in power consumption but differ in trends at higher
stroke and pitch combinations; the rigid fin consumes significantly more
power. The PDMS 1:10 fin has the largest FOM values, with the PDMS 1:20 fin
following. The rigid fin is significantly worse in all 3 metrics. Figure 7:
Figure of Merit results for a constant Stroke-Pitch-Offset of 0 degrees and
varying frequency. The rows graph all interpolated frequency values with
0.25Hz increments while the columns are different data sets (rigid, PDMS 1:10,
and PDMS 1:20). We verify that PDMS 1:10 is the best performing across all
frequencies with higher FOM values and averages compared to the PDMS 1:20 or
rigid fins. Additionally, increasing the frequency will improve the FOM
regardless of fin design.

To gain a better understanding of how our Figure of Merit correlates to the
thrust force and power consumption, we created a grid of contours for one
gait, shown in Figure 6. Here, the stroke and pitch range for a frequency of 2
Hz and Stroke-Pitch Offset of 0∘ is depicted. The grid contains columns with
the three material data sets and each row graphs part of the figure of merit
equation: FOM, thrust, and power in that order.

Beginning with the thrust, a few trends are immediately apparent. First, the
fin design that can generate the largest force is the PDMS 1:10 design,
generating a maximum force of around 1.6 N at a 40∘ stroke amplitude and
30-40∘ pitch amplitude. Both other designs trail behind, with the PDMS 1:20
fin being able to generate a thrust of 1.5 N and the rigid fin only managing
up to 1.1 N. The PDMS 1:10 fin also has more gaits at higher levels. An
observation of the contour reveals that there are more combinations of stroke
and pitch that produce higher thrusts when compared to the PDMS 1:20 fin. An
analysis of all gaits including other SPO and frequency confirms these
findings. The PDMS 1:10 fin’s maximum thrust is 2.1 N, while the PDMS 1:20 fin
produces a maximum force of 1.6 N and the rigid fin produces a maximum force
of 1.2 N. The PDMS 1:10 fin has the highest average thrust generation,
followed by the PDMS 1:20 fin.

Power consumption goes in the reverse order. The design that requires the
highest wattage is the rigid fin design, requiring 7.6 W for any gait with a
stroke amplitude of more than 40∘. The PDMS 1:10 and PDMS 1:20 fins are
similar, with a maximum power consumption of 7.1 W in this contour. However,
at cases above 40∘ stroke and 30∘ pitch, the PDMS 1:10 fin observes lower
wattage consumed at the same gait combination. The PDMS 1:20 fin appears to
depend less on the pitch amplitude, with power more dependent on the stroke
amplitude. These observations are verified when looking at all gaits. While
the PDMS fins have a maximum wattage of around 7.5 W, it occurs at much fewer
gaits than with the rigid data set.

Another interesting trend is visible when comparing the FOM and thrust charts,
which appear almost identical with only a few differences in their trends. The
explanation becomes evident when looking at the contours for power, which have
a near-linear trend across stroke amplitude. While pitch amplitude does affect
both the PDMS 1:10 and 1:20 designs, the stroke amplitude has the most
recognizable and significant effect. Future work will include generating
additional figures to verify that this trend exists across all frequencies and
stroke-pitch offsets.

This analysis allows us to conclude that the PDMS 1:10 fin design is the most
efficient out of all 3 designs, with the highest thrust generation and
efficiency. Following in second is the PDMS 1:20 fin, which has the second
largest thrust generation and efficiency. These trends are confirmed across
the ranges of stroke and pitch (Figure 6) as well as frequency (Figure 7) and
SPO (Figure 8). This suggests that the most efficient design that is able to
generate the largest thrust may lie between the two, and is something of
interest for future exploration.

From Figure 7, we observe that across the entire range of frequencies, the
PDMS 1:10 outperforms both the rigid and PDMS 1:20 fin designs. Additionally,
we observe that regardless of fin design, increasing the frequency will
improve the FOM metric.

From Figure 8, we observe that across the entire range of stroke-pitch offset,
the PDMS 1:10 outperforms both the rigid and PDMS 1:20 fin designs.
Additionally, we observe that regardless of fin design, a more negative offset
will slightly improve the FOM metric, although the difference is very
marginal. At high SPO, the PDMS 1:20 fin design appears to diverge from the
expected trends at 22.5 ∘ or higher and invites future exploration. In the
future, revisiting the model’s training data for high SPO will likely resolve
the issue.

Figure 8: Figure of Merit results for a constant frequency of 1.5Hz and
varying Stroke-Pitch-Offset. The rows graph all interpolated SPO values with
11.25∘ increments while the columns are different data sets (rigid, PDMS 1:10,
and PDMS 1:20). We verify that the PDMS 1:10 fin is the best performing across
all SPO values with higher FOM values and averages compared to the PDMS 1:20
or rigid fins. With the exception of 22.5-45 degrees for PDMS 1:20, a more
negative SPO will slightly improve the FOM regardless of fin design.

We conclude the following from the data interpolated:

  1. 1.

The best performing fin design is PDMS 1:10, followed by PDMS 1:20. Both
consume similar amounts of power, but PDMS 1:10 fins produce a higher thrust.
The rigid fin consumes more power and produces less thrust.

  2. 2.

The most optimal fin design likely lies in between the PDMS 1:10 and 1:20
fins.

  3. 3.

The most efficient gait will occur at a high stroke (40-55), centered pitch
(20-35), high frequency (2Hz) and low SPO (-22.5∘).

## Conclusion and Future Work

We use a forward-passing LSTM model to generate kinematic interpolations and
for fin gaits with the goal of integration onto a control system to optimize
for the efficiency of gaits and provide a better understanding of material
designs and their relation to efficiency.

We evaluate four different models for thrust and power to interpolate between
experimental data with high accuracy: the linear model, quartic polynomial
model, Convolutional Neural Network, and Long-Short-Term Memory model. All
four evaluated models accomplish all three criterion we laid out. They are
able to:

  * •

Complete a baseline model that inputs gait parameters to output either thrust
or power with minimal error

  * •

Retrain on different designs and maintain a similar or better accuracy

  * •

Run at a speed and size suitable for integration onto a control system ($>$100
computations per second)

Out of all four models, we find that the LSTM model is able to produce the
most accurate results on both the full data set and when we remove specific
gaits to create interpolations.

Using the generated interpolations, we develop a dimensionless Figure of Merit
that is able to compare our fin efficiency to other flapping systems and
evaluate the efficiency of gaits onboard our control system.

With the FOM, we conclude that both PDMS materials are more efficient than the
rigid fin, with the PDMS 1:10 fin generating the maximum thrust with the
lowest power consumption. There were observable trends consistent across all
materials, with a higher frequency, stroke amplitude, and negative Stroke-
Pitch-Offset all contributing to a greater efficiency. The most efficient gait
was concluded to be for the PDMS 1:10 data set with a -22.5∘ Stroke-Pitch-
Offset and frequency of 2 Hz. This understanding will allow us to both design
fins that generate a higher thrust and maintain the highest power efficiency,
and tune an inverse search model to search gaits it knows to be power-
efficient or optimal.

In the future, we aim to integrate the inverse search model with the FOM as a
weighting onto the control system to generate the most efficient or powerful
gaits. Before this integration is complete, we will combine both the thrust
and power LSTM models to decrease the model’s compute time on the control
system. Additionally, we will further verify our FOM by introducing physical
bounds on the thrust/power values a model can predict and conducting
additional tests with non-zero flow speeds. Introducing disturbances and
turbulence into the water ahead of time will aid in evaluation of model
efficiency prediction accuracy for dynamic environments similar to those the
UUV will swim in.

## References

  * Eloy (2012) Eloy, C. 2012.  Optimal Strouhal number for swimming animals.  _Journal of Fluids and Structures_ , 30: 205–218.
  * Geder et al. (2017) Geder, J.; Ramamurti, R.; Pruessner, M.; and Viswanath, K. 2017.  Underwater Thrust Performance of Tandem Flapping Fins: Effects of Stroke Phasing and Fin Spacing.  In _OCEANS 2017 - Anchorage_ , 1–7.
  * Geder et al. (2013) Geder, J. D.; Ramamurti, R.; Pruessner, M.; and Palmisano, J. 2013.  Maneuvering performance of a four-fin bio-inspired UUV.  In _2013 OCEANS - San Diego_ , 1–7.  ISSN: 0197-7385.
  * Geder et al. (2021) Geder, J. D.; Ramamurti, R.; Sampath, K.; Pruessner, M.; and Viswanath, K. 2021.  Fluid-Structure Modeling and the Effects of Passively Deforming Fins in Flapping Propulsion Systems.  In _OCEANS 2021: San Diego – Porto_ , 1–9.  ISSN: 0197-7385.
  * Lauder and Madden (2006) Lauder, G. V.; and Madden, P. G. A. 2006.  Learning from fish: Kinematics and experimental hydrodynamics for roboticists.  _International Journal of Automation and Computing_ , 3(4): 325–335.
  * Lee et al. (2022) Lee, J.; Viswanath, K.; Geder, J.; Sharma, A.; Pruessner, M.; and Zhou, B. 2022.  Data-Driven Machine Learning Models for a Multi-Objective Flapping Fin Unmanned Underwater Vehicle Control System.
  * Masud, La Mantia, and Dabnichki (2022) Masud, M. H.; La Mantia, M.; and Dabnichki, P. 2022.  Estimate of Strouhal and Reynolds numbers for swimming penguins.  _Journal of Avian Biology_ , 2022(2): e02886.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/jav.02886.
  * Mignano et al. (2019) Mignano, A. P.; Kadapa, S.; Tangorra, J. L.; and Lauder, G. V. 2019.  Passing the Wake: Using Multiple Fins to Shape Forces for Swimming.  _Biomimetics_ , 4(1): 23.  Number: 1 Publisher: Multidisciplinary Digital Publishing Institute.
  * Nedelcu et al. (2018) Nedelcu, A.-T.; Faităr, C.; Stan, L.-C.; and Buzbuchi, N. 2018.  Underwater Vehicle CFD Analyses and Reusable Energy Inspired by Biomimetic Approach.  preprint, ENGINEERING.
  * Nguyen, Lee, and Ahn (2016) Nguyen, P. L.; Lee, B. R.; and Ahn, K. K. 2016.  Thrust and swimming speed analysis of fish robot with non-uniform flexible tail.  _Journal of Bionic Engineering_ , 13(1): 73–83.
  * Palmisano et al. (2013) Palmisano, J. S.; Geder, J. D.; Pruessner, M. D.; and Ramamurti, R. 2013.  Power and Thrust Comparison of Bio-mimetic Pectoral Fins with Traditional Propeller-based Thrusters.  In _18th International Symposium on Unmanned Untethered Submersible Technology_ , 8.
  * Rohr and Fish (2004) Rohr, J. J.; and Fish, F. E. 2004.  Strouhal numbers and optimization of swimming by odontocete cetaceans.  _Journal of Experimental Biology_ , 207(10): 1633–1642.
  * Rohr et al. (1998) Rohr, J. J.; Hendricks, E. W.; Quiqley, L.; Fish, F. E.; and Gilpatrick, J. W. 1998.  Observations of Dolphin Swimming Speed and Strouhal Number.:.  Technical report, Defense Technical Information Center, Fort Belvoir, VA.
  * Sampath et al. (2020) Sampath, K.; Geder, J. D.; Ramamurti, R.; Pruessner, M. D.; and Koehler, R. 2020.  Hydrodynamics of tandem flapping pectoral fins with varying stroke phase offsets.  _Physical Review Fluids_ , 5(9): 094101.  Publisher: American Physical Society.
  * Tangorra et al. (2007) Tangorra, J. L.; Davidson, S. N.; Hunter, I. W.; Madden, P. G. A.; Lauder, G. V.; Dong, H.; Bozkurttas, M.; and Mittal, R. 2007.  The Development of a Biologically Inspired Propulsor for Unmanned Underwater Vehicles.  _IEEE Journal of Oceanic Engineering_ , 32(3): 533–550.
  * Taylor, Nudds, and Thomas (2003) Taylor, G. K.; Nudds, R. L.; and Thomas, A. L. R. 2003.  Flying and swimming animals cruise at a Strouhal number tuned for high power efficiency.  _Nature_ , 425(6959): 707–711.
  * Techet (2008) Techet, A. H. 2008.  Propulsive performance of biologically inspired flapping foils at high Reynolds numbers.  _Journal of Experimental Biology_ , 211(2): 274–279.
  * Triantafyllou, Triantafyllou, and Yue (2000) Triantafyllou, M. S.; Triantafyllou, G. S.; and Yue, D. K. P. 2000.  Hydrodynamics of Fishlike Swimming.  _Annual Review of Fluid Mechanics_ , 32(1): 33–53.
  * Yun, Kim, and Kim (2015) Yun, D.; Kim, K.-S.; and Kim, S. 2015.  Thrust characteristic of a caudal fin with spanwise variable phase.  _Ocean Engineering_ , 104: 344–348.