# Molecular excited state calculations with the QEB-ADAPT-VQE

Yordan S. Yordanov1,2 Crispin H. W. Barnes1 David R. M. Arvidsson-Shukur2,1 $\
{}^{1}$ Cavendish Laboratory, Department of Physics, University of Cambridge,
Cambridge CB3 0HE, United Kingdom  $\ {}^{2}$ Hitachi Cambridge Laboratory, J.
J. Thomson Avenue, CB3 0HE, Cambridge, United Kingdom

###### Abstract

Calculations of molecular spectral properties, like photodissociation rates
and absorption bands, rely on knowledge of the excited state energies of the
molecule of interest. Protocols based on the variational quantum eigensolver
(VQE) are promising candidates to calculate such energies on emerging noisy
intermediate scale quantum (NISQ) computers. The successful implementation of
these protocols on NISQ computers, relies on ansätze that can accurately
approximate the molecular states and that can be implemented by shallow
quantum circuits. In this paper, we introduce the excited qubit-excitation-
based adaptive (e-QEB-ADAPT)-VQE protocol to calculate molecular excited state
energies. The e-QEB-ADAPT-VQE constructs efficient problem-tailored ansätze by
iteratively appending evolutions of qubit excitation operators. The e-QEB-
ADAPT-VQE is an adaptation of the QEB-ADAPT-VQE protocol, which is designed to
be independent on the choice of an initial reference state. We perform
classical numerical simulations for LiH and BeH2 to benchmark the performance
of the e-QEB-ADAPT-VQE. We demonstrate that the e-QEB-ADAPT-VQE can construct
highly accurate ansätze that require at least an order of magnitude fewer
$CNOT$s than standard fixed UCC ansätze, such as the UCCSD and the GUCCSD.

##  I Introduction

Quantum molecular simulations with the variational quantum eigensolver (VQE)
McArdle  _et al._ (2020); McClean  _et al._ (2016); O’Malley  _et al._ (2016);
Wang  _et al._ (2019); Arute  _et al._ (2020) are a promising application for
emerging noisy intermediate-scale quantum (NISQ) Preskill (2018); Arute  _et
al._ (2019); Elfving  _et al._ (2020) computers. The VQE is a hybrid quantum-
classical algorithm that utilizes the Rayleigh Ritz variational principle to
determine the lowest eigenvalue of a Hamiltonian operator, by optimizing an
ansatz. In particular the VQE can be used to solve the electronic structure
problem Helgaker  _et al._ (2014), and find the ground state configuration and
energy of a molecule. By utilizing both a quantum and a classical computer,
the VQE is less quantum hardware demanding at the expense of requiring more
quantum measurements and classical post-processing, as compared to purely
quantum algorithms, like the Quantum Phase Estimation (QPE) algorithm Dorner
_et al._ (2009). A major challenge for the practical realization of a
molecular VQE simulation on NISQ computers is to construct a variationally
flexible ansatz that: (1) accurately approximates the eigenstates of $H$; (2)
is easy to optimize; and (3) can be implemented by a shallow circuit that uses
few $2$-qubit entangling gates, e.g. $CNOT$s, which are the current bottleneck
of NISQ computers.

The most widely used type of ansätze for molecular VQE simulations, are the
unitary coupled cluster (UCC) ansätze Hempel  _et al._ (2018); Romero  _et
al._ (2018); Lee  _et al._ (2018); Dallaire-Demers  _et al._ (2019); Sokolov
_et al._ (2020); Peruzzo  _et al._ (2014). They were motivated by the
classical coupled cluster theory Helgaker  _et al._ (2014), and correspond to
products of fermionic excitation evolutions. Due to their fermionic structure,
UCC ansätze preserve many of the physical symmetries of electronic
wavefunctions, which makes them accurate and easy to optimize at the same
time. State of the art ADAPT-VQE protocols Grimsley  _et al._ (2019); Tang
_et al._ (2019); Sim  _et al._ (2020); Daniel Claudino and Humble (2020),
iteratively construct problem-tailored ansätze, which consists of close to
optimal number of fermionic excitation evolutions. These ansätze are
implemented by shallow ansatz circuits, and have few variational parameters,
while at the same time can be highly accurate.

In a previous work Yordanov  _et al._ (2021), we introduced the qubit-
excitation based (QEB)-ADAPT-VQE. Unlike the original fermionic-ADAPT-VQE
Grimsley  _et al._ (2019), the QEB-ADAPT-VQE constructs a problem-tailored
ansatz by appending “qubit excitation evolutions” (unitary evolutions of
excitation operators that satisfy “qubit commutation relations” Yordanov  _et
al._ (2021); Wu and Lidar (2002); Yordanov  _et al._ (2020); Xia and Kais
(2020)). Qubit excitation evolutions are implemented by simpler circuits than
the standard fermionic excitation evolutions Yordanov  _et al._ (2020). Also,
as demonstrated in Refs. Tang  _et al._ (2019); Yordanov  _et al._ (2021); Xia
and Kais (2020), both types of excitation evolutions approximate electronic
wavefunctions comparably well. Thus, in Ref. Yordanov  _et al._ (2021) we
demonstrated that the QEB-ADAPT-VQE can simulate molecular ground states,
using shallower quantum circuits than the previous state of the art fermionic-
ADAPT-VQE Grimsley  _et al._ (2019) and qubit-ADAPT-VQE Tang  _et al._ (2019)
protocols.

In this work, we present a modified version of the QEB-ADAPT-VQE designed to
simulate molecular excited states, which we call the excited(e)-QEB-ADAPT-VQE.
The QEB-ADAPT-VQE, similarly to other VQE protocols, largely relies on an
initial reference state, $|\psi_{0}\rangle$, that has a large overlap with the
true ground state, $|E_{0}\rangle$. In particular, the QEB-ADAPT-VQE evaluates
the energy gradients of individual qubit excitation evolutions in order to
decide the best way to grow its ansatz at each iteration. These energy
gradients are evaluated for zero values of the variational parameters on the
presumption that $|\psi_{0}\rangle$ is close to $|E_{0}\rangle$. However,
choosing an initial reference state that is guaranteed to be close an unknown
excited state is challenging. We modify the ansatz growing strategy of the
QEB-ADAPT-VQE, so that it does not require $|\psi_{0}\rangle$ to be close to
the target (ground or excited) state. Additionally, we present a comparison of
qubit and fermionic excitation evolutions in their ability to construct
ansätze for excited states.

The paper is organized as follows: In Sec. II we include a theoretical
introduction. Section III describes the steps of the e-QEB-ADAPT-VQE protocol,
and in Sec. IV we benchmark its performance for the first excited states of
LiH and BeH2. Lastly, Sec. V presents a comparison between qubit and fermionic
excitation evolutions.

##  II Theory

###  II.1 The electronic structure problem and the variational quantum
eigensolver

Finding the ground-state $|E_{0}\rangle$ and corresponding energy $E_{0}$ of a
molecule is known as the “electronic structure problem” Helgaker  _et al._
(2014). This problem can be solved by solving the time-independent Schrödinger
equation $H|E_{0}\rangle=E_{0}|E_{0}\rangle$, where $H$ is the electronic
Hamiltonian of the molecule. Within the Born-Oppenheimer approximation, where
the nuclei of the molecule are assumed to be motionless, $H$ can be written in
second quantized form as

${H}=\sum_{i,k}^{N_{{}_{MO}}}h_{i,k}^{\ }{a}_{i}^{\dagger}{a}_{k}^{\
}+\sum_{i,j,k,l}^{N_{{}_{MO}}}h_{i,j,k,l}^{\
}{a}_{i}^{\dagger}{a}_{j}^{\dagger}{a}_{k}^{\ }{a}_{l}^{\ }.$ (1)

Here, $N_{{}_{MO}}$ is the number of considered spin-orbitals,
${a}^{\dagger}_{i}$ and ${a}_{i}$ are the fermionic ladder operators,
corresponding to the $i^{th}$ molecular spin-orbital, and the factors $h_{ij}$
and $h_{ijkl}$ are one- and two-electron integrals, written in a spin-orbital
basis Helgaker  _et al._ (2014). The Hamiltonian expression in Eq. 1 can be
mapped to quantum-gate operators using a qubit encoding method, e.g. the
Jordan-Wigner (JW) Wigner and Jordan (1928) or the Bravyi-Kitaev Bravyi and
Kitaev (2002) methods. Throughout this work, we assume the more
straightforward JW encoding, where the occupancy of the $i^{th}$ molecular
spin-orbital is represented by the state of the $i^{th}$ qubit.

The fermionic ladder operators ${a}^{\dagger}_{i}$ and ${a}_{i}$ satisfy anti-
commutation relations

$\\{a_{i},a^{\dagger}_{j}\\}=\delta_{i,j},\ \
\\{a_{i},a_{j}\\}=\\{a_{i}^{\dagger},a^{\dagger}_{j}\\}=0.$ (2)

Within the JW encoding, ${a}^{\dagger}_{i}$ and ${a}_{i}$ can be written in
terms of quantum gate operators as

$a_{i}^{\dagger}=Q_{i}^{\dagger}\prod_{r=0}^{i-1}Z_{r}=\frac{1}{2}(X_{i}-iY_{i})\prod_{r=0}^{i-1}Z_{r}\;\;\mathrm{and}$
(3)
$a_{i}=Q_{i}\prod_{r=0}^{i-1}Z_{r}=\frac{1}{2}(X_{i}+iY_{i})\prod_{r=0}^{i-1}Z_{r},$
(4)

where

$Q_{i}^{\dagger}\equiv\frac{1}{2}(X_{i}-iY_{i})\;\;\mathrm{and}\
Q_{i}\equiv\frac{1}{2}(X_{i}+iY_{i}).$ (5)

We refer to ${Q_{i}^{\dagger}}$ and ${Q_{i}}$ as qubit creation and
annihilation operators, respectively. The operators $Q^{\dagger}_{i}$ and
$Q_{i}$ satisfy the “qubit commutation relations” Wu and Lidar (2002);
Yordanov  _et al._ (2020, 2021)

$\displaystyle\\{Q_{i},Q_{i}^{\dagger}\\}=I,[Q_{i},Q^{\dagger}_{j}]=0\text{ if
}\ i\neq j,\;\;\mathrm{and}$
$\displaystyle\;\;[Q_{i},Q_{j}]=[Q_{i}^{\dagger},Q^{\dagger}_{j}]=0\text{ for
all }i,j.$ (6)

We will use ${Q_{i}^{\dagger}}$ and ${Q_{i}}$ in Sec. II.2 to define qubit
excitation evolutions.

Substituting Eqs. (3) and (4) into Eq. (1), $H$ can be written as

${H}=\sum_{r}h_{r}\prod_{s=0}^{N_{MO}-1}\sigma_{s}^{r},$ (7)

where $\sigma_{s}$ is a Pauli operator ($X_{s}$, $Y_{s}$, $Z_{s}$ or $I_{s}$)
acting on qubit $s$, and $h_{r}$ (not to be confused with $h_{ik}$ and
$h_{ijkl}$) is a real scalar coefficient. The number of terms in Eq. (7)
scales as $O(N_{{}_{MO}}^{4})$. The expectation value of $H$ can be evaluated
by individually measuring on a quantum computer the expectation values of all
Pauli strings in Eq. (7).

Once $H$ is mapped to a Pauli string representation, the VQE can be used to
minimize the expectation value
$E(\boldsymbol{\theta})=\langle\psi(\boldsymbol{\theta})|H|\psi(\boldsymbol{\theta})\rangle$,
where $|\psi(\boldsymbol{\theta})\rangle$ is a trial state. The VQE relies
upon the Rayleigh-Ritz variational principle

$\langle\psi(\boldsymbol{\theta})|H|\psi(\boldsymbol{\theta})\rangle\geq
E_{0},$ (8)

to find an lower-bounded estimate for $E_{0}$. The VQE is a hybrid-quantum-
classical algorithm that uses a quantum computer to prepare the trial state
$|\psi(\boldsymbol{\theta})\rangle$ and evaluate $E(\boldsymbol{\theta})$, and
a classical computer to process the measurement data and update
$\boldsymbol{\theta}$ at each iteration. The trial state
$|\psi(\boldsymbol{\theta})\rangle=U(\boldsymbol{\theta})|\psi_{0}\rangle$ is
generated by the ansatz $U(\boldsymbol{\theta})$ applied to an initial
reference state $|\psi_{0}\rangle$.

###  II.2 Qubit and fermionic excitation evolutions

The most widely used type of ansätze for molecular VQE simulations, are the
unitary coupled cluster (UCC) ansätze, motivated by the classical coupled
cluster method Helgaker  _et al._ (2014). UCC ansätze are constructed as
products of fermionic excitation evolutions. The QEB-ADAPT-VQE and the e-QEB-
ADAPT-VQE instead construct ansätze that consist of qubit excitation
evolutions. This subsection provides definitions of the two types of unitary
operations.

####  II.2.1 Fermionic excitation evolutions

Single and double fermionic excitation operators, are defined, respectively,
by the skew-Hermitian operators

$T_{ik}\equiv a^{\dagger}_{i}a_{k}-a^{\dagger}_{k}a_{i}\;\;\mathrm{and}$ (9)
$T_{ijkl}\equiv
a^{\dagger}_{i}a^{\dagger}_{j}a_{k}a_{l}-a^{\dagger}_{k}a^{\dagger}_{l}a_{i}a_{j}.$
(10)

Single and double fermionic excitation evolutions are thus given,
respectively, by the unitary operators

$A_{ik}(\theta)=e^{\theta
T_{ik}}=\exp\left[\theta(a^{\dagger}_{i}a_{k}-a^{\dagger}_{k}a_{i})\right]\;\;\mathrm{and}$
(11) $A_{ijkl}(\theta)=e^{\theta
T_{ijkl}}=\exp\left[\theta(a^{\dagger}_{i}a^{\dagger}_{j}a_{k}a_{l}-a^{\dagger}_{k}a^{\dagger}_{l}a_{i}a_{j})\right].$
(12)

Using equations (3) and (4), for $i<j<k<l$, $A_{ik}$ and $A_{ijkl}$ can be
expressed in terms of quantum gate operators as

$A_{ij}(\theta)=\exp\left[i\frac{\theta}{2}(X_{i}Y_{k}-Y_{i}X_{k})\prod_{r=i+1}^{k-1}Z_{r}\right]\;\;\mathrm{and}$
(13) $\displaystyle
A_{ijkl}(\theta)=\exp\Bigg{[}i\frac{\theta}{8}(X_{i}Y_{j}X_{k}X_{l}+Y_{i}X_{j}X_{k}X_{l}+Y_{i}Y_{j}Y_{k}X_{l}$
$\displaystyle+Y_{i}Y_{j}X_{k}Y_{l}-X_{i}X_{j}Y_{k}X_{l}-X_{i}X_{j}X_{k}Y_{l}$
$\displaystyle-
Y_{i}X_{j}Y_{k}Y_{l}-X_{i}Y_{j}Y_{k}Y_{l})\prod_{r=i+1}^{j-1}Z_{r}\prod_{r^{\prime}=k+1}^{l-1}Z_{r^{\prime}}\Bigg{]}.$
(14)

$CNOT$-efficient circuits to implemented single and double fermionic
excitation evolutions were derived in Ref. Yordanov  _et al._ (2020).

####  II.2.2 Qubit excitation evolutions

Qubit excitation evolutions are generated by “qubit excitation operators”.
Single and double qubit excitation operators are defined, respectively, by the
skew-Hermitian operators

$\tilde{T}_{ik}=Q^{\dagger}_{i}Q_{k}-Q^{\dagger}_{k}Q_{i}\;\;\mathrm{and}$
(15)
$\tilde{T}_{ijkl}=Q^{\dagger}_{i}Q^{\dagger}_{j}Q_{k}Q_{l}-Q^{\dagger}_{k}Q^{\dagger}_{l}Q_{i}Q_{j}.$
(16)

Hence, single and double qubit excitation evolutions are defined,
respectively, by the unitary operators

$\tilde{A}_{ik}(\theta)=e^{\theta\tilde{T}_{ik}}=\exp\big{[}\theta(Q^{\dagger}_{i}Q_{k}-Q^{\dagger}_{k}Q_{i})\big{]}\;\;\mathrm{and}$
(17)
$\tilde{A}_{ijkl}(\theta)=e^{\theta\tilde{T}_{ijkl}}=\exp\big{[}\theta(Q^{\dagger}_{i}Q^{\dagger}_{j}Q_{k}Q_{l}-Q^{\dagger}_{k}Q^{\dagger}_{l}Q_{i}Q_{j})\big{]}.$
(18)

Using Eq. (5), $\tilde{A}_{ik}$ and $\tilde{A}_{ijkl}$ can be re-expressed in
terms of quantum gate operators:

$\tilde{A}_{ik}(\theta)=\exp\Big{[}i\frac{\theta}{2}(X_{i}Y_{k}-Y_{i}X_{k})\Big{]}\;\;\mathrm{and}$
(19)
$\displaystyle\tilde{A}_{ijkl}(\theta)=\exp\Big{[}i\frac{\theta}{8}(X_{i}Y_{j}X_{k}X_{l}+Y_{i}X_{j}X_{k}X_{l}$
$\displaystyle+Y_{i}Y_{j}Y_{k}X_{l}+Y_{i}Y_{j}X_{k}Y_{l}-X_{i}X_{j}Y_{k}X_{l}$
$\displaystyle-
X_{i}X_{j}X_{k}Y_{l}-Y_{i}X_{j}Y_{k}Y_{l}-X_{i}Y_{j}Y_{k}Y_{l})\Big{]}.$ (20)

$CNOT$-efficient circuits to implement single and double qubit excitation
evolutions were derived in Ref. Yordanov  _et al._ (2020).

###  II.3 Finding excited state energies with the VQE

The two most common methods to calculate excited state energies with the VQE
are the quantum subspace expansion McClean  _et al._ (2020, 2017); Colless
_et al._ (2018), and the overlap-based method Higgott  _et al._ (2019); Jones
_et al._ (2019). In this work we use the latter overlap-based method.

As the name suggest the overlap-based method works by including the overlap
between previously found, eigenstates of $H$ and the state we are currently
searching for in the cost function of the VQE. For example, after we find the
ground state $|E_{0}\rangle$ we modify H as

$H\rightarrow H_{1}=H+\alpha_{0}|E_{0}\rangle\langle E_{0}|,$ (21)

where $\alpha_{0}$ is a real positive scalar coefficient. If
$\alpha_{0}+E_{0}>E_{1}$ 111To make sure that $\alpha_{0}+E_{0}>E_{1}$, prior
to knowing $E_{1}$, we can choose an arbitrary value for $\alpha_{0}$, which
is large compared to the energy scale of the problem., where $E_{1}$ is the
energy of the first excited state of $H$, the lowest energy eigenvalue of the
modified Hamiltonian $H_{1}$ will be shifted to $E_{1}$. Hence, running the
VQE for $H_{1}$ should output an estimate for $E_{1}$:

$E_{1}=\min_{\boldsymbol{\theta_{1}}}\langle\psi_{0}|U_{1}^{\dagger}(\boldsymbol{\theta_{1}})H_{1}U_{1}(\boldsymbol{\theta_{1}})|\psi_{0}\rangle=\\\
\min_{\boldsymbol{\theta_{1}}}\langle\psi_{0}|U_{1}^{\dagger}(\boldsymbol{\theta_{1}})\big{(}H+\alpha_{0}|E_{0}\rangle\langle
E_{0}|\big{)}U_{1}(\boldsymbol{\theta_{1}})|\psi_{0}\rangle$ (22)

After $E_{1}$ is found the same procedure can be repeated recursively multiple
times to find the next energy eigenvalues of $H$. The $k^{\text{th}}$ energy
eigenvalue of $H$, $E_{k}$, can be estimated as

$E_{1}=\min_{\boldsymbol{\theta_{k}}}\langle\psi_{0}|U_{k}^{\dagger}(\boldsymbol{\theta_{k}})H_{k}U_{k}(\boldsymbol{\theta_{k}})|\psi_{0}\rangle=\\\
\min_{\boldsymbol{\theta_{k}}}\langle\psi_{0}|U_{k}^{\dagger}(\boldsymbol{\theta_{k}})\big{(}H+\sum_{r=0}^{k-1}\alpha_{r}|E_{r}\rangle\langle
E_{r}|\big{)}U_{k}(\boldsymbol{\theta_{k}})|\psi_{0}\rangle$ (23)

The first term in Eq. (23), is calculated as described above, by measuring the
expectation values of the Pauli string terms in the expression for $H$ [Eq.
(7)]. The overlap terms in Eq. (23) can be calculated with the SWAP test
Cincio  _et al._ (2018).

##  III The excited-QEB-ADAPT-VQE

The e-QEB-ADAPT-VQE algorithm finds an estimate for the $k^{\text{th}}$
excited state of an electronic Hamiltonian, $H$, by constructing a problem-
tailored ansatz. The ansatz is constructed by iteratively appending qubit
excitation evolutions, based on the greedy strategy to obtain the lowest
estimate for the energy, $E(\boldsymbol{\theta})$, at each iteration. In this
subsection, we describe the preparation components, and the iterative ansatz-
constructing loop, of the e-QEB-ADAPT-VQE algorithm. For a description of the
original QEB-ADAPT-VQE, we refer the reader to Ref. Yordanov  _et al._ (2021).

First, we transform the $H$ to a quantum-gate-operator representation as
described in Sec. II.1. This involves the calculation of the one- and two-
electron integrals $h_{ik}$ and $h_{ijkl}$ [Eq. (1)], which can be done
efficiently (in time polynomial in $N_{MO}$) on a classical computer. If we
want to find the energy of the excited state $|E_{k}\rangle$, $E_{k}$, we add
the additional overlap terms to $H$, as described in Sec. II.3, to get
$H_{k}$.

Second, we define an ansatz element pool $\mathbb{P}(\tilde{A},N_{MO})$ of all
unique single and double qubit excitations, $\tilde{A}_{ik}(\theta)$ and
$\tilde{A}_{ijkl}(\theta)$, respectively, for $i,j,k,l\in\\{0,N_{MO}-1\\}$.
The size of this pool is
$|\mathbb{P}(\tilde{A},N_{MO})|=\binom{N_{MO}}{2}+3\binom{N_{MO}}{4}$.

Third, we choose an initial reference state $|\psi_{0}\rangle$. As we
mentioned in the introduction, the performance of the e-QEB-ADAPT-VQE is
intended to be independent of the choice of $|\psi_{0}\rangle$. Hence, we
arbitrarily choose the initial reference state to be the Hartree-Fock state,
but in principle any computational basis state with a Hamming weight equal to
the number of electrons, $N_{e}$, will do.

Now we are ready to begin constructing the ansatz. We set the iteration number
to $m=1$, the ansatz to the identity $U\rightarrow U^{(0)}=I$, and initiate
the iterative ansatz-constructing loop. We describe the five steps of the
$m^{th}$ iteration below. Afterwards we comment on these steps.

  1. 1.

Prepare trial state
$|\psi^{[m-1]}\rangle=U\left(\boldsymbol{\theta}^{[m-1]}\right)|\psi_{0}\rangle$,
with values for $\boldsymbol{\theta}^{[m-1]}$ as determined in the previous
iteration.

  2. 2.

For each (single or double) qubit excitation evolution
$\tilde{A_{p}}(\theta_{p})\in\mathbb{P}(\tilde{A},N_{MO})$:

    1. (a)

Run the VQE to find

$\min\limits_{\theta_{p}}E(\theta_{p})=\min\limits_{\theta_{p}}\langle\psi^{[m-1]}|\tilde{A_{p}}^{\dagger}(\theta_{p})H_{k}\tilde{A_{p}}(\theta_{p})|\psi^{[m-1]}\rangle.$

    2. (b)

Calculate the energy reduction magnitude
$E^{[m-1]}-\min\limits_{\theta_{p}}E(\theta_{p})$

  3. 3.

Identify the set of $n$ qubit excitation evolutions,
$\mathbb{\tilde{A}}^{[m]}(n)$, corresponding to the $n$ largest energy
reductions in the previous step. For
$\tilde{A_{p}}(\theta_{p})\in\mathbb{\tilde{A}}^{[m]}(n)$:

    1. (a)

Run the VQE to find
$\min\limits_{\boldsymbol{\theta}^{[m-1]},\theta_{p}}E\left(\boldsymbol{\theta}^{[m-1]},\theta_{p}\right)=$

$\min\limits_{\boldsymbol{\theta}^{[m-1]},\theta_{p}}\langle\psi_{0}|U^{\dagger}\big{(}\boldsymbol{\theta}^{[m-1]}\big{)}\tilde{A_{p}}^{\dagger}(\theta_{p})H_{k}\tilde{A_{p}}(\theta_{p})U\big{(}\boldsymbol{\theta}^{[m-1]}\big{)}|\psi_{0}\rangle.$

    2. (b)

Calculate the energy reduction
$\Delta{E}^{[m]}_{p}=E^{[m-1]}-\min\limits_{\boldsymbol{\theta}^{[m-1]},\theta_{p}}E\left(\boldsymbol{\theta}^{[m-1]},\theta_{p}\right)$.

    3. (c)

Save the (re)optimized values of
$\boldsymbol{\theta}^{[m-1]}\cup\\{\theta_{p}\\}$ as
$\boldsymbol{\theta}^{[m]}_{p}$.

  4. 4.

Identify the qubit excitation evolution
$\tilde{A}^{[m]}\left(\theta^{[m]}\right)\equiv\tilde{A}_{p^{\prime}}(\theta_{p^{\prime}})$
corresponding to the largest energy reduction
$\Delta{E}^{[m]}\equiv\Delta{E}_{p^{\prime}}^{[m]}$ in the previous step.

If $\Delta E^{[m]}<\epsilon$, where $\epsilon>0$ is an energy threshold:

    1. (a)

Exit

Else:

    1. (a)

Append $\tilde{A}^{[m]}\left(\theta^{[m]}\right)$ to the ansatz:

$U\big{(}\boldsymbol{\theta}^{[m-1]}\big{)}\rightarrow
U\big{(}\boldsymbol{\theta}^{[m]}\big{)}=\tilde{A}^{[m]}\big{(}\theta^{[m]}\big{)}U\big{(}\boldsymbol{\theta}^{[m-1]}\big{)}$
(24)

    2. (b)

Set $E^{[m]}=E^{[m-1]}-\Delta{E}^{[m]}_{p}$

    3. (c)

Set the values of the new set of variational parameters,
$\boldsymbol{\theta}^{[m]}=\boldsymbol{\theta}^{[m-1]}\cup\\{\theta_{p^{\prime}}\\}$,
to $\boldsymbol{\theta}^{[m]}_{p^{\prime}}$

  5. 5.

Enter the $m+1$ iteration by returning to step $1$

We now elaborate on the steps outlined above.

The $m^{\text{th}}$ iteration of the e-QEB-ADAPT-VQE starts by preparing the
trial state $|\psi^{[m-1]}\rangle$ obtained in the $(m-1)^{\text{th}}$
iteration.

At each iteration we aim to identify the qubit excitation evolution that, when
appended to the ansatz, would maximise the energy reduction
$E(\boldsymbol{\theta}^{[m-1]})-E(\boldsymbol{\theta}^{[m]})$. To identify
such a qubit excitation evolution, first we calculate (step 2) the individual
energy reduction contributions of each qubit excitation evolution
$\tilde{A}_{p}(\theta_{p})\in\mathbb{P}(\tilde{A},N_{MO})$. Each of these
energy reductions is calculated by a single-parameter VQE optimization
performed to minimize the energy expectation value
$\langle\psi^{[m-1]}|\tilde{A}_{p}^{\dagger}(\theta_{p})H\tilde{A}_{p}(\theta_{p})|\psi^{[m-1]}\rangle$.
In this way, we get an indication by how much the energy expectation value
will be reduced when each $\tilde{A}_{p}(\theta_{p})$ is appended to the
ansatz $U(\boldsymbol{\theta}^{[m-1]})$ and the energy is minimized along the
full set of parameters $\boldsymbol{\theta}^{[m-1]}\cup\\{\theta_{p}\\}$, at a
reduced cost of performing $O(N_{MO}^{4})$ single-parameter VQE optimizations
(instead of $O(N_{MO}^{4})$ $m$-parameter VQE optimizations).

In the original QEB-ADAPT-VQE such an indication is obtained by measuring the
energy gradients
$\left\\{\frac{\partial}{\partial\theta_{p}}\langle\psi^{[m-1]}|\tilde{A}_{p}^{\dagger}(\theta_{p})H\tilde{A}_{p}(\theta_{p})|\psi^{[m-1]}\rangle\right\\}$.
However, these energy gradients must be evaluated for some values of the
respective parameters $\left\\{\theta_{p}\right\\}$. If we are approximating
the ground state $|E_{0}\rangle$ and we have a high level of confidence that
$|\psi_{0}\rangle$ has a large overlap with $|E_{0}\rangle$, then the
gradients can be evaluated conveniently for $\left\\{\theta_{p}=0\right\\}$.
But if we are approximating the excited state $|E_{k}\rangle$, where we do not
know how close $|\psi_{0}\rangle$ and $|E_{k}\rangle$ are, measuring the
energy gradients for $\left\\{\theta_{p}=0\right\\}$, or any other arbitrary
values of $\\{\theta_{p}\\}$, can be a poor indicator, leading to an
inefficient and slow construction of the ansatz, and a possibility of getting
stuck in local energy minimum (see Sec. IV.1). Instead the e-QEB-ADAPT-VQE
measures the individual energy reduction contribution for each qubit
excitation evolution, which can be obtained by a single-parameter VQE
optimization, using a direct search minimizer Kokail  _et al._ (2019), e.g.
the Nelder Mead Nelder and Mead (1965), irrespectively of the overlap of
$|\psi_{0}\rangle$ and $|E_{k}\rangle$.

If we denote by $N_{H}$ the number of terms in the Pauli string representation
of $H$ [Eq. (7)], calculating a single energy-gradient requires measuring
$2N_{H}$ expectation values. On the other hand, running a single-parameter VQE
requires measuring $\gamma N_{H}$ expectation values, where $\gamma$ is the
number of function evaluations required for a single parameter minimization.
Hence, the additional cost of the technique pursued here in comparison to
using energy gradients as in the original QEB-ADAPT-VQE, is a factor of
$\gamma/2$ more quantum computer measurements in step 2222For the Nelder Mead
$\gamma$ is on the order of $10$.. In addition, once the ansatz reaches some
critical size, so that $|\psi^{[m-1]}\rangle$ can be assumed to have a
significant overlap with $|E_{k}\rangle$, we can switch to the cheaper QEB-
ADAPT-VQE.

The individual energy reductions calculated in step 2, indicate how much each
qubit excitation evolution can decrease $E^{(m-1)}$ when appended to
$U(\boldsymbol{\theta}^{(m-1)})$. However, the largest individual energy
reduction does not necessarily correspond to the largest energy reduction when
the ansatz is optimized over all variational parameters. In step 3, we
identify the set of $n$ qubit excitation evolutions with the individual energy
reductions: $\mathbb{\tilde{A}}^{[m]}(n)$. We assume that
$\mathbb{\tilde{A}}^{[m]}(n)$ likely contains the qubit excitation evolution
that reduces $E^{[m-1]}$ the most. For each of the $n$ qubit excitation
evolutions in $\mathbb{\tilde{A}}^{[m]}(n)$, we run the VQE, for all
variational parameters, with the ansatz from the previous iteration to
calculate how much it contributes to the energy reduction. If multiple quantum
devices are available, step 3 can be parallelized.

In step 4, we pick the qubit excitation evolution,
$\tilde{A}^{[m]}\left(\theta^{[m]}\right)$, corresponding to the largest
energy reduction, $\Delta{E}^{[m]}$. If $\Delta{E}^{[m]}$ is below some
threshold $\epsilon>0$, we exit the iterative loop. If instead the
$|\Delta{E}^{[m]}|>\epsilon$, we add
$\tilde{A}^{[m]}\left(\theta^{[m]}\right)$ to the ansatz and begin the next
iteration.

Figure 1: Energy convergence plots for the first excited states of LiH and
BeH2 in the STO-3G basis at equilibrium bond distances of
$r_{\text{Li-H}}=1.546{\textup{\AA}}$ and
$r_{\text{Be-H}}=1.316{\textup{\AA}}$, respectively. The plots are obtained
with the QEB-ADAPT-VQE and the e-QEB-ADAPT-VQE for $n=10$.

##  IV Benchmarking the e-QEB-ADAPT-VQE

In this section, we benchmark the performance of the e-QEB-ADAPT-VQE by
finding the excited state energies of LiH and BeH2, which have been simulated
on real quantum computers Peruzzo  _et al._ (2014); Hempel  _et al._ (2018).
The results presented here are based on classical numerical VQE simulations,
performed with a custom in-house code333The code and the numerical data are
available upon request from the authors. . The Hamiltonians for LiH and BeH2
are represented in the STO-3G orbital basis set Ditchfield  _et al._ (1971),
with no frozen orbitals assumed. Hence, the wavefunction of LiH is represented
by a $12$ qubit state, and that of BeH2 by a $14$ qubit state. In step 2 of
the e-QEB-ADAPT-VQE, the VQE optimization is performed with the direct-search
Nelder Mead method Nelder and Mead (1965), and in step 3 with the gradient
descent BFGS method Fletcher (2013) .

###  IV.1 The QEB-ADAPT-VQE vs the e-QEB-ADAPT-VQE

In this subsection we compare the convergence rates, to the first excited
state energies of LiH and BeH2, of the QEB-ADAPT-VQE and the e-QEB-ADAPT-VQE.
With this comparison we demonstrate that the energy-gradient-based ansatz-
growing strategy of the QEB-ADAPT-VQE is unsuitable for excited states.

Figure 1 presents energy convergence plots for the first excited states of LiH
and BeH2 in the STO-3G basis at equilibrium bond distances of
$r_{\text{Li-H}}=1.546{\textup{\AA}}$ and
$r_{\text{Be-H}}=1.316{\textup{\AA}}$, respectively. The plots are obtained
with each of the two protocols. In the case of LiH the QEB-ADAPT-VQE converges
slower, requiring more than twice as many ansatz-constructing iterations, than
the e-QEB-ADAPT-VQE. Since each ansatz-constructing iteration corresponds to a
qubit excitation evolution in the ansatz, the ansatz constructed by the QEB-
ADAPT-VQE is roughly twice as large as the one constructed by the e-QEB-ADAPT-
VQE. As suggested in Sec. III above, the underlying reason for this is the
inability of the energy-gradient-based ansatz-growing strategy of the QEB-
ADAPT-VQE to grow an efficient ansatz when the initial reference state,
$|\psi_{0}\rangle$, does not have a significant overlap with the target
excited state, $|E_{k}\rangle$. In the case of BeH2, the QEB-ADAPT-VQE
completely fails to converge to chemical accuracy, likely getting stuck in a
local energy minimum. Nevertheless, the e-QEB-ADAPT-VQE, being independent on
the overlap of $|\psi_{0}\rangle$ and $|E_{k}\rangle$, converges successfully
for both molecules.

###  IV.2 Energy dissociation curves

In this section, we benchmark the performance of the e-QEB-ADAPT-VQE, for
$\epsilon=10^{-6}$ Hartree and $\epsilon=10^{-8}$ Hartree, by obtaining energy
dissociation plots for the first excited states of LiH and BeH2. For a
comparison, we also include energy dissociation plots obtained with the VQE
using the standard UCCSD Hempel  _et al._ (2018); Romero  _et al._ (2018)
ansatz, which consists of single and double fermionic excitation evolutions
above the Hartree-Fock state, and the larger generalized UCCSD (GUCCSD) ansatz
Lee  _et al._ (2018), which consists of all unique single and double fermionic
excitation evolutions.

The energy dissociation curves for the two molecules are given in Figs. 2a and
2b. Figures 2c and 2d show the errors of each method with respect to the full
configuration interaction (FCI) energy Helgaker  _et al._ (2014, 2012). In the
case of LiH, all methods achieve chemical accuracy. Figures 2e and 2f show the
numbers of variational parameters, which are also equivalent to the number of
qubit/fermionic excitation evolutions, used in the ansatz of each method. The
ansätze constructed by the e-QEB-ADAPT-VQE, for LiH, are extremely compact,
consisting of at most $27$ qubit evolutions, whereas the UCCSD and the GUCCSD
consist of $200$ and $1521$ fermionic evolutions, respectively.

In the case of BeH2 the results are more complicated. First, the UCCSD-VQE
fails to achieve chemical accuracy for the majority of bond distances (Fig.
2d). This result is not surprising as the simple UCCSD ansatz is not suitable
to approximate strongly correlated states such as excited states. Second, in
Fig. 2d we see that the e-QEB-ADAPT-VQE fails to achieve chemical accuracy for
bond distances of $r_{Be-H}=1.5\textup{\AA}$ and $r_{Be-H}=1.75\textup{\AA}$.
We comment on the reason for this in Sec. IV.3. The only method that achieves
chemical accuracy at all bond distances is the GUCCSD-VQE, because of the
highly variationally flexible GUCCSD ansatz. Nevertheless, for the majority of
bond distances the e-QEB-ADAPT-VQE constructs ansätze that are more accurate
and consist of nearly $50$ times fewer ansatz elements than the GUCCSD (Figs.
2d and 2f).

These results together with the fact that qubit evolutions are implemented by
simpler circuits than fermionic evolutions, implies that the ansätze
constructed by the e-QEB-ADAPT-VQE are implemented by much shallower circuits
that have much fewer $CNOTs$ than the UCCSD and the GUCCSD ansätze. Table 1
summarizes the $CNOT$ counts for the UCCSD and the GUCCSD ansätze, and the
ansätze constructed by the e-QEB-ADAPT-VQE($\epsilon=10^{-8}$) for each
molecule. For the e-QEB-ADAPT-VQE($\epsilon=10^{-8}$), the maximum $CNOT$
counts over all bond distances for each molecule are shown. As we can see in
the table, the ansätze constructed by the e-QEB-ADAPT-VQE use at least $10$
times fewer $CNOT$s than the UCCSD, and more than $60$ times fewer $CNOT$s
than the GUCCSD.

| UCCSD | GUCCSD | e-QEB-ADAPT-VQE
---|---|---|---
LiH | 3496 | 29447 | 311
BeH2 | 8980 | 64064 | 896
Table 1: $CNOT$ counts for the UCCSD, the GUCCSD, and the ansätze constructed
by the e-QEB-ADAPT-VQE($\epsilon=10^{-8}$) for the first excited states of LiH
and BeH2 in the STO-3G basis.  Figure 2:  Energy dissociation curves for LiH
and $\text{BeH}_{2}$ in the STO-3G orbital basis. a,b Absolute values for the
estimated first excited state energies. c,d Error in the estimated energy
values with respect to the exact FCI energy. e,f Number of variational
parameters, also equal to the number of qubit/fermionic evolutions, of the
ansatz used by each method.

###  IV.3 Finding the wrong excited state energy

As we saw in Fig. 2d, the e-QEB-ADAPT-VQE fails to achieve chemical accuracy
for BeH2 at bond distances $r_{Be-H}=1.5\textup{\AA}$ and
$r_{Be-H}=1.75\textup{\AA}$. In order to investigate the reason for this see
Fig. 3, which depicts the FCI energies for the $10$ lowest energy states of
BeH2 (in the STO-3G basis) as function of bond distance.

The $9$ excited states shown in Fig. 3 are ordered in three degenerate energy
levels. The table below summarises the energy levels and their corresponding
excited states at bond distances $r_{Be-H}=1.5\textup{\AA}$ and
$r_{Be-H}=1.75\textup{\AA}$:

$r_{Be-H}=1.5\text{\AA}$
---
| exc. states | energy, Hartree
---|---
{1,2} | $E_{I}(1.5\text{\AA})\approx-15.3343$
{3,4,5,6,7,8} | $E_{II}(1.5\text{\AA})\approx-15.3331$
{9} | $E_{III}(1.5\text{\AA})\approx-15.3026$
$r_{Be-H}=1.75\text{\AA}$
| exc. states | energy, Hartree
---|---
{1,2} | $E_{I}(1.75\text{\AA})\approx-15.3131$
{3,4,5} | $E_{II}(1.75\text{\AA})\approx-15.3059$
{6,7,8,9} | $E_{III}(1.75\text{\AA})\approx-15.3056$
Table 2: Energy levels for the $9$ lowest excited states of BeH2 in the STO-3G
basis. The ground state (not shown) is counted as the $0^{\text{th}}$.

For $r_{Be-H}=1.5\text{\AA}$ we can see that $E_{I}(1.5\text{\AA})$ and
$E_{II}(1.5\text{\AA})$ are very close. Upon expectation of the energy
estimate of the e-QEB-ADAPT-VQE($\epsilon=10^{-8}$), for
$r_{Be-H}=1.5\text{\AA}$, it turns out that the algorithm actually finds
$E_{II}(1.5\text{\AA})$ with accuracy of $3\times 10^{-8}$ Hartree. Therefore,
the problem is not a lack of accuracy, but converging to the wrong excited
state.

Two possible reasons for this are: (1) the classical minimizers employed by
the e-QEB-ADAPT-VQE get stuck in a local minimum corresponding to
$E_{II}(1.5\text{\AA})$, and/or (2) the e-QEB-ADAPT-VQE fails to construct an
ansatz that can approximate the degenerate excited states corresponding to
energy level $E_{I}(1.5\text{\AA})$. Since the e-QEB-ADAPT-VQE uses a
combination of the BFGS and the Nelder-Mead optimization methods, where the
Nelder-Mead is a direct-search method, which is unlikely to get stuck in local
minimum, the first possible reason can be ruled out. Hence, we are left with
the second one.

Upon explicit inspection of the degenerate states corresponding to energy
level $E_{I}(1.5\text{\AA})$, it is found that they have strong static
correlations Helgaker  _et al._ (2012), where two dominant Slater determinants
contributed equally to the wavefunction. On the other hand for each of the
degenerate states corresponding to energy level $E_{II}(1.5\text{\AA})$ there
is one dominant Slater determinant. Hence, the first few ansatz-constructing-
iterations of the e-QEB-ADAPT-VQE grow the ansatz with qubit excitation
evolutions, which map between the initial reference state, $|\psi_{0}\rangle$,
and the dominant Slater determinant of one of the states corresponding to
$E_{II}(1.5\text{\AA})$, because these qubit excitation evolutions decreases
the estimate for the energy by the most. From that point on, the e-QEB-ADAPT-
VQE continues to construct an ansatz that approximates a state corresponding
to $E_{II}(1.5\text{\AA})$, instead a state corresponding to
$E_{I}(1.5\text{\AA})$.

In the case of $r_{Be-H}=1.75\text{\AA}$, the situation is similar. However,
in this case the three energy levels $E_{I}(1.75\text{\AA})$,
$E_{II}(1.75\text{\AA})$ and $E_{III}(1.75\text{\AA})$ are close, and the QEB-
ADAPT-VQE converges to an excited state corresponding to
$E_{III}(1.75\text{\AA})$. Again the states corresponding to
$E_{I}(1.75\text{\AA})$ and $E_{II}(1.75\text{\AA})$ have two (not equally)
dominant Slater determinants, whereas the states corresponding to
$E_{III}(1.75\text{\AA})$ have weaker static correlations, and just one
dominant Slater determinant.

Figure 3: FCI energies for the $10$ lowest energy states of BeH2 in the STO-3G
basis. The numbers of the excited states correspond to their order at the
ground state equilibrium bond distance of $r_{Be-H}=1.316\AA$.

Overall the following conclusion can be made: A generic feature of the e-QEB-
ADAPT-VQE is to be more “willing” to construct ansätze for less statically
correlated states (where one Slater determinant is dominant). This owes to its
greedy strategy to grow its ansatz with the qubit excitation evolution that
decrease the energy estimate by the most at each ansatz-growing iteration.
Hence, in the case when (1) the lowest energy eigenstates of an electronic
Hamiltonian are ordered in two or more closely spaced (possibly degenerate)
energy levels, and (2) the lowest energy level corresponds to eigenstates that
have stronger static correlations than the eigenstates above them, the e-QEB-
ADAPT-VQE might converge to one of the less statically correlated eigenstates,
instead to one with the lowest energy.

Figure 4: Energy convergence plots for the first excited states of LiH and
BeH2 in the STO-3G basis at equilibrium bond distances of
$r_{\text{Li-H}}=1.546{\textup{\AA}}$ and
$r_{\text{Be-H}}=1.316{\textup{\AA}}$, respectively. The plots are obtained
with the e-QEB-ADAPT-VQE for $n_{qe}=10$.

##  V Qubit vs fermionic excitation evolutions

The use of qubit excitation evolutions by the QEB-ADAPT-VQE was motivated by
the fact that they are implemented by simpler quantum circuits than fermionic
excitation evolutions Yordanov  _et al._ (2020, 2021), and the presumption
that a qubit-excitation-based ansatz can approximate an electronic
wavefunction as well as a fermionic-excitation-based ansatz. In Ref. Yordanov
_et al._ (2021) we compared the two types of ansätze in approximating the
ground states for small molecules, and found that the ansätze perform
similarly, in terms of accuracy. In this section, we generalise this
comparison for (more strongly correlated) excited states.

Figure 4 shows energy convergence plots for the first excited states of LiH
and BeH2, at bond distances of $r_{\text{Li-H}}=1.546{\textup{\AA}}$ and
$r_{\text{Be-H}}=1.316{\textup{\AA}}$, respectively. The blue plots are
obtained with the e-QEB-ADAPT-VQE for $n=10$ (the same plots as in Fig. IV.1),
and the red plots are obtained with the e-QEB-ADAPT-VQE ansatz-constructing
routine, for $n=10$, using a pool of single and double fermionic (instead of
qubit) excitation evolutions. We can observe close similarity between the two
types of ansätze. In the case of LiH (Fig. 4a) the qubit-excitation-based
ansatz is slightly more accurate per number of ansatz elements, while in the
case of BeH2 (Fig. 4b) the fermionic-excitation-based ansatz is a bit more
accurate instead. These differences are relatively small, so we derive the
same conclusion as in Refs. Yordanov  _et al._ (2021); Xia and Kais (2020),
that qubit and fermionic excitation evolutions can approximate electronic
wavefunctions, comparably well.

##  VI Conclusion

In this paper we proposed a modified version of the QEB-ADAPT-VQE algorithm,
the e-QEB-ADAPT-VQE, designed to simulate low-lying excited molecular states.

The QEB-ADAPT-VQE relies on an initial reference state, which has a large
overlap with the target ground state in order to construct an efficient
ansatz. However, choosing an initial reference state that has a significant
overlap with a target excited state is not as straightforward as for a target
ground state. Therefore, the e-QEB-ADAPT-VQE is designed to be less dependent
on the choice of an initial reference state. This is achieved by an ansatz
growing strategy that does not rely on energy-gradient evaluations, like the
QEB-ADAPT-VQE, but instead on individual energy-reduction evaluations. The
modified ansatz growing strategy comes at a cost of up to a constant factor of
more quantum computer measurements, as compared to the original QEB-ADAPT-VQE.

We benchmarked the performance of the e-QEB-ADAPT-VQE with classical numerical
simulations, by constructing ansätze for the first excited states of LiH and
BeH2. We found that the e-QEB-ADAPT-VQE can construct highly accurate ansätze
that are implemented by multiple times shallower circuits, which require
multiple times fewer $CNOT$s, than fixed UCC ansätze, like the UCCSD and the
GUCCSD. Therefore, the e-QEB-ADAPT-VQE is especially suitable for NISQ
computers, for which the number of $2$-qubits entangling gates, such as
$CNOT$s, that can be applied reliably is the current bottleneck Preskill
(2018).

We also found that the e-QEB-ADAPT-VQE might fail to find excited-state
energies in order of increasing energy, if the excited states are ordered in
closely spaced degenerate energy levels, and the lower lying excited states
are more statically correlated than the states above them. This problem
derives from the fact that the e-QEB-ADAPT-VQE pursues the greedy strategy to
achieve a lowest estimate for the energy at each ansatz-growing iteration.
Hence, the algorithm runs the risk of going along the wrong path and converge
to a state, which is not necessarily the one with the lowest energy. In fact,
this problem will be present in any iterative-VQE protocol that relies on such
a greedy strategy, including the QEB-ADAPT-VQE, the ADAPT-VQE Grimsley  _et
al._ (2019) and the qubit-ADAPT-VQE Tang  _et al._ (2019). Nevertheless, this
problem is unlikely to be significant, because in practice one would be
interested in finding the whole spectrum of low lying excited state energies,
so the order in which they are found is not necessarily important.

Lastly, we used the e-QEB-ADAPT-VQE ansatz-constructing routine to compare
qubit-excitation-based and fermionic-excitation-based ansätze in approximating
excited states. The two types of ansätze were found to perform similarly,
achieving a particular accuracy with approximately the same number of
qubit/fermionic excitation evolutions. These observations generalize previous
results Yordanov  _et al._ (2021); Xia and Kais (2020) about the equivalence
of qubit and fermionic excitation evolutions in approximating electronic
wavefunctions.

###### Acknowledgements.

The authors wish to thank Dr N. Mertig for useful discussions. Y.S.Y.
acknowledges financial support from the EPSRC and Hitachi via CASE
studentships RG97399.

## References

  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Rev. Mod. Phys. 92, 015003 (2020).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New Journal of Physics 18, 023023 (2016).
  * O’Malley  _et al._ (2016) P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis, Phys. Rev. X 6, 031007 (2016).
  * Wang  _et al._ (2019) D. Wang, O. Higgott, and S. Brierley, Phys. Rev. Lett. 122, 140504 (2019).
  * Arute  _et al._ (2020) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, B. B. Buckley, D. A. Buell, _et al._ , arXiv preprint arXiv:2004.04174 (2020), 10.1126/science.abb9811.
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Nature 574, 505 (2019).
  * Elfving  _et al._ (2020) V. E. Elfving, B. W. Broer, M. Webber, J. Gavartin, M. D. Halls, K. P. Lorton, and A. Bochevarov,  (2020), arXiv:2009.12472 .
  * Helgaker  _et al._ (2014) T. Helgaker, P. Jorgensen, and J. Olsen,  _Molecular electronic-structure theory_ (John Wiley & Sons, 2014).
  * Dorner  _et al._ (2009) U. Dorner, R. Demkowicz-Dobrzanski, B. J. Smith, J. S. Lundeen, W. Wasilewski, K. Banaszek, and I. A. Walmsley, Phys. Rev. Lett. 102, 040403 (2009).
  * Hempel  _et al._ (2018) C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz, H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush, A. Aspuru-Guzik, R. Blatt, and C. F. Roos, Phys. Rev. X 8, 031022 (2018).
  * Romero  _et al._ (2018) J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, Quantum Science and Technology 4, 014008 (2018).
  * Lee  _et al._ (2018) J. Lee, W. J. Huggins, M. Head-Gordon, and K. B. Whaley, Journal of chemical theory and computation 15, 311 (2018).
  * Dallaire-Demers  _et al._ (2019) P.-L. Dallaire-Demers, J. Romero, L. Veis, S. Sim, and A. Aspuru-Guzik, Quantum Science and Technology 4, 045005 (2019).
  * Sokolov  _et al._ (2020) I. O. Sokolov, P. K. Barkoutsos, P. J. Ollitrault, D. Greenberg, J. Rice, M. Pistoia, and I. Tavernelli, The Journal of Chemical Physics 152, 124107 (2020).
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, Nature communications 5, 4213 (2014).
  * Grimsley  _et al._ (2019) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nature communications 10, 1 (2019).
  * Tang  _et al._ (2019) H. L. Tang, E. Barnes, H. R. Grimsley, N. J. Mayhall, and S. E. Economou, arXiv preprint arXiv:1911.10205 (2019).
  * Sim  _et al._ (2020) S. Sim, J. Romero, J. F. Gonthier, and A. A. Kunitsa, arXiv preprint arXiv:2010.00629 (2020).
  * Daniel Claudino and Humble (2020) A. J. M. Daniel Claudino, Jerimiah Wright and T. S. Humble, arXiv preprint arXiv:2011.01279 (2020).
  * Yordanov  _et al._ (2021) Y. S. Yordanov, V. Armaos, C. H. Barnes, and D. R. Arvidsson-Shukur, Commun Phys 4, 228 (2021), https://doi.org/10.1038/s42005-021-00730-0.
  * Wu and Lidar (2002) L.-A. Wu and D. Lidar, Journal of Mathematical Physics 43, 4506 (2002).
  * Yordanov  _et al._ (2020) Y. S. Yordanov, D. R. M. Arvidsson-Shukur, and C. H. W. Barnes, Phys. Rev. A 102, 062612 (2020).
  * Xia and Kais (2020) R. Xia and S. Kais, Quantum Science and Technology (2020), https://doi.org/10.1088/2058-9565/abbc74.
  * Wigner and Jordan (1928) E. Wigner and P. Jordan, Z. Phys 47, 631 (1928).
  * Bravyi and Kitaev (2002) S. B. Bravyi and A. Y. Kitaev, Annals of Physics 298, 210 (2002).
  * McClean  _et al._ (2020) J. R. McClean, Z. Jiang, N. C. Rubin, R. Babbush, and H. Neven, Nature Communications 11, 1 (2020).
  * McClean  _et al._ (2017) J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong, Phys. Rev. A 95, 042308 (2017).
  * Colless  _et al._ (2018) J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. de Jong, and I. Siddiqi, Phys. Rev. X 8, 011021 (2018).
  * Higgott  _et al._ (2019) O. Higgott, D. Wang, and S. Brierley, Quantum 3, 156 (2019).
  * Jones  _et al._ (2019) T. Jones, S. Endo, S. McArdle, X. Yuan, and S. C. Benjamin, Phys. Rev. A 99, 062304 (2019).
  * Cincio  _et al._ (2018) L. Cincio, Y. Subaşı, A. T. Sornborger, and P. J. Coles, New Journal of Physics 20, 113022 (2018).
  * Kokail  _et al._ (2019) C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, _et al._ , Nature 569, 355 (2019).
  * Nelder and Mead (1965) J. A. Nelder and R. Mead, The Computer Journal 7, 308 (1965), https://academic.oup.com/comjnl/article-pdf/7/4/308/1013182/7-4-308.pdf .
  * Ditchfield  _et al._ (1971) R. Ditchfield, W. J. Hehre, and J. A. Pople, The Journal of Chemical Physics 54, 724 (1971).
  * Fletcher (2013) R. Fletcher,  _Practical methods of optimization_ (John Wiley & Sons, 2013).
  * Helgaker  _et al._ (2012) T. Helgaker, S. Coriani, P. Jorgensen, K. Kristensen, J. Olsen, and K. Ruud, Chemical reviews 112, 543 (2012).


# Calculation of binding energy and wave function for exotic hidden-charm
pentaquark

F. Chezani Sharahi1,2, M. Monemzadeh1
1Department of Physics, University of kashan, Iran
2 Department of Physics, University of Arak, Iran
<EMAIL_ADDRESS>author: Majid Monemzadeh<EMAIL_ADDRESS>
###### Abstract

In this study, pentaquark $P_{c}(4380)$ composed of a baryon $\Sigma_{c}$ and
a $\overline{D}^{*}$ meson is considered. Pentaquark is as a bound state of
two-body systems composed of a baryon and a meson. The calculated potential
will be expanded and replaced in the Schrödinger equation until tenth
sentences of expansion. Solving the Schrödinger equation with the expanded
potential of Pentaquark leads to an analytically complete approach. As a
consequence, the binding energy $E_{B}$ of pentaquark $P_{c}$ and wave
function are obtained. The results $E_{B}$ will presented in the form of
tables so that we can review the existence of pentaquark $P_{c}$. Then, the
wave function will be shown on diagrams. Finally, the calculated results are
compared with another obtained results and the mass of observing pentaquark
$P_{c}$ and the radius of pentaquark are estimated.

Keywords: Pentaquark , Exotic Hidden-charm, Binding energy,

##  1 Introduction

In the first researches, the existence of multi-quark states are illustrated
in its simplest possible form, in which baryons made of three fundamental
quarks and mesons from a quark and an antiquark[1]. Indeed, from the point of
view of the mathematics and physics, there was no QCD theorem opposing the
existence of exotic multi-quark. In gauge field theory, the QCD principle
allows the existence of multi - quarks and hybrids, which include quark and
gluonic degrees of freedom[2].

Searching to find the pentaquark and its probing has a long history. About ten
years ago, in the study of pentaquark, a big progress was occurred when LEPS
collaborations getting started with a claim to find out strong evidence of the
pentaquark with the mass of about $1.540GeV$[3]. Then, many theoretical and
experimental approaches were pursued and many ideas were proposed in this
field. For example, Zou and his colleagues suggested that the components of
pentaquark be included with nucleon. Since the heavy quarks play an important
role in stabilization of multi - quark systems and these play exactly the same
role as the hydrogen molecule in QED[4, 5, 6, 7], so there are theoretical
predictions about exotic hidden-charm pentaquark. In particular, the
possibility existence of hidden-charm molecular baryons composed of an anti-
charmed meson and a charmed baryon was systematically studied within the one
boson exchange model[3].

In 2003, LEPS collaboration reported the evidence of pentaquark $s$ state with
content of quark $uudd\bar{s}$ and very narrow width[8]. Unfortunately, this
exotic flavor was not confirmed in subsequent experiments[9, 10]. In fact, the
possible theoretical arguments are presented for the non-existence of a stable
pentaquark $s$ in the references[11, 12]. Also, this mode has not been found
with a light flavor yet. However, the baryons with light flavor may be able to
have a significant pentaquark components[2].

One decade ago, numerous researches were done in around of the world to find
exotic particles. The result of these efforts was observation the mesonic X ,
Y , Z particles in Belle, BESIII, BABAR and LHCb. Some of these were
considered as candidates for exotic states because they don’t fit in a regular
mesonic structure[1]. The common point of the exotic states is that all of
these contain the heavy quarks and antiquarks.

Due to the heavy quarks, the exotic states can be stable and light modes can
be combined with a regular mode[13]. This guess is consistent with the fact
that all of the exotic states have a hidden $c$ or $b$ which this is
experimentally observed. If this claim is valid that the heavy components
stabilize multi-quark systems by particle physics scientists, it will a
natural proposal from their experimental colleagues to be search for exotic
states consist of clear quarks $c$ and $b$, for example $b\bar{c}qq^{\prime}$
and etc[2].

Figure 1: (Color line)The s-wave charmed baryons with $J^{p}=1/2^{+}$ and
pseudo- scalar or vector s-wave charmed mesons which form the molecular
baryons with double-charm[6].

Here, one hidden-charm molecular baryon composed of an anti-charmed meson and
a charmed baryon has been studied. This can have one of the two flavors i.e.
symmetric $6_{F}$ or anti-symmetric flavor $\bar{3}_{F}$ as shown in Fig. 1.
So spin-parity of s-wave charmed baryon is $J^{p}=1/2^{+}$ or $J^{p}=3/2^{+}$
for $6_{F}$ and $J^{p}=1/2^{+}$ for $\bar{3}_{F}$. The pseudo- scalar or
vector anti-charmed meson is made of s-wave anti-charmed mesons. In Figure 1,
hidden-charm molecular states composed of the anti-charmed mesons and the
charmed baryons are placed inside the green range[6].

Recently, LHCb collaboration has observed two resonance structures
$P_{c}(4380)$ and $P_{c}(4450)$ with mass and width decay
$M_{P_{c}(4380)}=4380\pm 8\pm 29MeV$, $\Gamma_{P_{c}(4380)}=205\pm 18\pm
86MeV$ and $M_{P_{c}(4450)}=4449.8\pm 1.7\pm 2.5MeV$,
$\Gamma_{P_{c}(4450)}=39\pm 5\pm 19MeV$ in the invariant mass spectrum $J/\psi
p$ from $\Lambda_{b}\rightarrow J/\psi pK$. According to the final state
$J/\psi p$ is concluded that two observing states $P_{c}$, unable to be
isosinglet and these are consist of hidden-charm quantum numbers. Also, it is
suggested structures for each of the states $P_{c}$ which is considered for
$P_{c}(4380)$, $\Sigma_{c}(2455)\overline{D}^{*}$ and $P_{c}(4450)$,
$\Sigma_{c}^{*}(2520)\overline{D}^{*}$[2, 3].

Binding energy may be calculated analytically by solving Schrödinger equation
with the expanded potential of Pentaquark. This approach offers advantages
over numerical solution of the Schrödinger equation for pentaquark. First all
values for binding energy can be calculated and this is done more accurately.
Second wave function can be computed and presented graphically ($\phi(r)$
versus $r$), from which useful particle characteristics and data could be
extracted.

In this work, it is described by analytical solution of pentaquark
$\Sigma_{c}\overline{D}^{*}$ in four sections: where apart from the
introduction as section 1, the potential of pentaquark $P_{c}(4380)$ is
discussed in section 2. Then, the Schrödinger equation using the potential of
pentaquark is analytically solved in section 3. Finally, important conclusions
are discussed in section 4.

##  2  The potential of pentaquark $P_{c}$ and its expansion

Pentaquark is considered to be composed of one baryon and one meson. For
pentaquark, with such a structure, the potential is given as follows[2]:

$V_{\Sigma_{c}\overline{D}^{*}}(r)=\frac{1}{3}\ \frac{gg_{1}}{f_{\pi}}\
\nabla^{2}Y(\Lambda,m_{\pi},r)\mathcal{J}_{0}\mathcal{G}_{0},$ (1)

where, the coupling constant $g=0.59\pm 0.07\pm 0.01$ is extracted from the
width of $D^{*}$ [14, 15], $g_{1}=0.75,0.94,1.95$ [6, 16]. Also, it is the
mass of pion $m_{\pi}=135MeV$ and pion decay constant $f_{\pi}=132MeV$ [2].
The amount of phenomenological cutoff parameter is considered
$\Lambda=0.8GeV-2.5GeV$ [2, 12]. Finally, the $Y(\Lambda,m_{\pi},r)$ is [2]:

$\displaystyle Y(\Lambda,m_{\pi},r)=\frac{1}{4\pi r}\
(e^{-m_{\pi}r}-e^{-\Lambda r})-\frac{\Lambda^{2}-m_{\pi}^{2}}{8\pi\Lambda}\
e^{-\Lambda r}.$ (2)

Now, we are calculating and expanding $\nabla^{2}Y$:

$\displaystyle\nabla^{2}Y(\Lambda,m_{\pi},r)$ $\displaystyle=$
$\displaystyle\frac{m_{\pi}^{2}}{4\pi r}\ (e^{-m_{\pi}r}-e^{-\Lambda
r})-\frac{\Lambda^{3}-m_{\pi}^{2}\Lambda}{8\pi}\ e^{-\Lambda r}$ (3)
$\displaystyle=$
$\displaystyle\frac{1}{4\pi}(b_{0}+b_{1}r+b_{2}r^{2}+b_{3}r^{3}+b_{4}r^{4}+b_{5}r^{5}$
$\displaystyle+b_{6}r^{6}+b_{7}r^{7}+b_{8}r^{8}+b_{9}r^{9}+...).$

Where:

$\displaystyle b_{0}$ $\displaystyle=$
$\displaystyle(-m^{3}+\frac{\Lambda^{3}}{2}+\frac{m^{2}\Lambda}{2});\ \ \ \ \
\ \ \ \ \ b_{1}=(\frac{m^{4}}{2}-\frac{\Lambda^{4}}{2});$ $\displaystyle
b_{2}$ $\displaystyle=$
$\displaystyle(-\frac{m^{5}}{3!}+\frac{\Lambda^{5}}{2\times
2!}-\frac{m^{2}\Lambda^{3}}{6\times 2!});\ \
b_{3}=(\frac{m^{6}}{4!}-\frac{\Lambda^{6}}{2\times
3!}+\frac{m^{2}\Lambda^{4}}{4\times 3!});$ $\displaystyle b_{4}$
$\displaystyle=$ $\displaystyle(-\frac{m^{7}}{5!}+\frac{\Lambda^{7}}{2\times
4!}-\frac{3m^{2}\Lambda^{5}}{6\times 4!});\ \
b_{5}=(\frac{m^{8}}{6!}-\frac{\Lambda^{8}}{2\times
5!}+\frac{m^{2}\Lambda^{6}}{3\times 5!});$ $\displaystyle b_{6}$
$\displaystyle=$ $\displaystyle(-\frac{m^{9}}{7!}+\frac{\Lambda^{9}}{2\times
6!}-\frac{5m^{2}\Lambda^{7}}{14\times 6!});\ \
b_{7}=(\frac{m^{10}}{8!}-\frac{\Lambda^{10}}{2\times
7!}+\frac{3m^{2}\Lambda^{8}}{8\times 7!});$ $\displaystyle b_{8}$
$\displaystyle=$ $\displaystyle(-\frac{m^{11}}{9!}+\frac{\Lambda^{11}}{2\times
8!}-\frac{7m^{2}\Lambda^{9}}{18\times
8!});b_{9}=(\frac{m^{12}}{10!}-\frac{\Lambda^{12}}{2\times
9!}+\frac{2m^{2}\Lambda^{10}}{5\times 9!}).$ (4)

##  3 Solving the Schrödinger equation for pentaquark $P_{c}$

To investigate the existence of one bound state of pentaquark, we solve the
Schrödinger equation with calculated potential in previous section for
pentaquark $P_{c}$.

The radial Schrödinger equation for two-body systems is:

$(\frac{d^{2}}{dr^{2}}+\frac{2}{r}\frac{d}{dr}-\frac{l(l+1)}{r^{2}})R_{n,l}(r)+\frac{2\mu}{\hbar}(E-V(r))R_{n,l}(r)=0$
(5)

Taking $\hbar=1$ and changing the variable $\phi(r)=rR_{n,l}(r)$, Eq. (5)
becomes:

$\frac{d^{2}}{{dr}^{2}}\ \phi(r)+2\mu(E-V(r)-\frac{l(l+1)}{2\mu
r^{2}})\phi(r)=0\ \ \ .$ (6)

By placing the potential $V(r)=V_{\Sigma_{c}\overline{D}^{*}}(r)$ and the
expanded shape of $\nabla^{2}Y$ in Eq. (6), an equation is expressed as
follows:

$\displaystyle\frac{d^{2}}{dr^{2}}\phi(r)+2\mu(E-C_{0}-C_{1}r-C_{2}r^{2}-C_{3}r^{3}-C_{4}r^{4}-C_{5}r^{5}$
$\displaystyle-
C_{6}r^{6}-C_{7}r^{7}-C_{8}r^{8}-C_{9}r^{9}+...-\frac{l(l+1)}{2\mu
r^{2}})\phi(r)=0.$ (7)

Where:

$\displaystyle C_{n}=B_{1}b_{n};\ \ \ \ n=0,1,...,9;\ \ \ \
B_{1}=\frac{1}{12\pi}\frac{gg_{1}}{f_{\pi}^{2}}\mathcal{J}_{0}\mathcal{G}_{0};\
\ .$ (8)

In Eq. (8), For $\Sigma_{c}\overline{D}^{*}$ With $I=1/2,J=3/2$, is placing
the numerical value of $\mathcal{J}_{0}\mathcal{G}_{0}=1$ product from
literature[2].

By considering the following proposed reply (cf. [17, 18]) for the
differential Eq. (7) yields:

$\phi(r)=N(r)e^{M(r)}=r^{n}e^{M(r)}.$ (9)

Differentiating second degree of $\phi$ in Eq.(9) gives:

$\phi^{\prime\prime}(r)=(n(n-1)r^{-2}+2nM^{\prime}r^{-1}+M^{\prime\prime}+{M^{\prime}}^{2})r^{n}e^{M}.$
(10)

Here, we solving differential equation (7) by considering expansion $V(r)$ up
to the $10^{th}$ order for calculating binding energy $E_{B}$ of pentaquark.
Hence, an approximation was attempted up to the $10^{th}$ order, not only
expanded potential behavior up to the $10^{th}$ order is the same as potential
behavior in equation (1), but also resulting of binding energy have sufficient
precision compared to numerical literature[2, 19]. This is indicating adequacy
of approximation.

Now, we can be considered two position for $\phi(r)$. The first $\phi(0)=0$ is
for $n=1$ and the second $\phi(0)=cte$ is considered for $n=0$. Thus, we study
them in separated subsection.

###  3.1 The position $\phi(0)=0$

To considered $\phi(0)=0$ then $M(r)$ is Yields[19]:

$\displaystyle M(r)$ $\displaystyle=$ $\displaystyle
a_{1}r^{2}+a_{2}r^{3}+a_{3}r^{4}+a_{4}r^{5}+a_{5}r^{6}+a_{6}r^{7}$ (11)
$\displaystyle+a_{7}r^{8}+a_{8}r^{9}+a_{9}r^{10}.$

Replacing Eq. (11) into the Eq. (10) and comparing two Eqs. of (7), (10), the
following expression (Eq.(12)) for $r^{-2}$ is obtained in terms of $l(l+1)$
in Eq.(8), as well as, a system of 11 non-linear equations expressed later:

$n(n-1)=l(l+1)$ (12)

In the base state, two values for $n$ namely $n=0$ and $n=1$ is obtained that
according to Eq. (9) and the condition $\phi(0)=0$, the value $n=0$ could be
unacceptable and $N(r)=r$ will be. Therefore, $\phi(r)$ ,
$\phi^{\prime\prime}(r)$ are obtained as follows:

$\displaystyle\phi(r)=re^{M(r)}$
$\displaystyle\phi^{\prime\prime}(r)=(2M^{\prime}r^{-1}+M^{\prime\prime}+{M^{\prime}}^{2})re^{M}.$
(13)

After replacing Eq.(11) in Eq. (13) and equal to the different powers of r,
Separately, the following nonlinear equations are obtained:

$\displaystyle 6a_{1}=-2\mu(E-C_{0})$ $\displaystyle 12a_{2}=2\mu C_{1}$
$\displaystyle 20a_{3}+4a_{1}^{2}=2\mu C_{2}$ $\displaystyle
30a_{4}+12a_{1}a_{2}=2\mu C_{3}$ $\displaystyle
42a_{5}+9a_{2}^{2}+16a_{1}a_{3}=2\mu C_{4}$ $\displaystyle
56a_{6}+20a_{1}a_{4}+24a_{2}a_{3}=2\mu C_{5}$ $\displaystyle
72a_{7}+16a_{3}^{2}+24a_{1}a_{5}+30a_{2}a_{4}=2\mu C_{6}$ $\displaystyle
90a_{8}+28a_{1}a_{6}+36a_{2}a_{5}+40a_{3}a_{4}=2\mu C_{7}$ $\displaystyle
110a_{9}+25a_{4}^{2}+32a_{1}a_{7}+42a_{2}a_{6}+48a_{3}a_{5}=2\mu C_{8}$
$\displaystyle
132a_{10}+36a_{1}a_{8}+48a_{2}a_{7}+56a_{3}a_{6}+60a_{4}a_{5}=2\mu C_{9}.$
(14)

Here, by replacing the numerical values of constants for pentaquark
$P_{c}(4380)$, $\Sigma_{c}=2455MeV$ and $\overline{D}^{*}=2008.32MeV$, we
obtained the binding energy of pentaquark[2]. In table 1, 2 have been shown
different values of binding energy from Eq. (14) for $P_{c}(4380)$.

Table 1:  $(E_{B}(MeV),M_{P_{c}}(MeV))$ for $P_{c}(4380)$ and $g_{1}=0.94$ $\Lambda(MeV)$ | $g=0.51$ | $g=0.59$ | $g=0.67$
---|---|---|---
$800$ | $(-63,\ 4400.32)$ | $(-40,\ 4423.32)$ | $(-17,\ 4446.32)$
$850$ | $(-62,\ 4401.32)$ | $(-34,\ 4429.32)$ | $(-6.7,\ 4456.62)$
$900$ | $(-59,\ 4404.32)$ | $(-26,\ 4437.32)$ | -
$1000$ | $(-45,\ 4418.32)$ | $(-0.11,\ 4463.21)$ | -
$1100$ | $(-22,\ 4441.32)$ | - | -
Table 2:  $(E_{B}(MeV),M_{P_{c}}(MeV))$ for $P_{c}(4380)$ and $g_{1}=0.75$ $\Lambda(MeV)$ | $g=0.51$ | $g=0.59$ | $g=0.67$
---|---|---|---
$800$ | $(-90,\ 4373.32)$ | $(-73,\ 4390.32)$ | $(-56,\ 4407.32)$
$850$ | $(-95,\ 4368.32)$ | $(-75,\ 4390.32)$ | $(-53,\ 4407.32)$
$900$ | $(-99,\ 4364.32)$ | $(-74,\ 4389.32)$ | $(-49,\ 4414.32)$
$1000$ | $(-102,\ 4361.32)$ | $(-67,\ 4369.32)$ | $(-32,\ 4431.32)$
$1100$ | $(-98,\ 4365.32)$ | $(-51,\ 4412.32)$ | $(-3.14,\ 4460.18)$
$1200$ | $(-86,\ 4377.32)$ | - | -

In the tables above $M_{P_{c}}$ is calculated as follows:

$\displaystyle M_{P_{c}}=M_{\Sigma_{c}}+M_{\overline{D}^{*}}+E_{B}.$ (15)

According to the obtained values, it is observed that in $g_{1}=0.75$, the
results for pentaquark mass are much closer to the $M_{P_{c}(4380)}$. Also, to
check the results, one of the main differences between this paper and the
other references [2, 20] is that the acceptable results are obtained for the
$M_{P_{c}(4380)}$ only in the $800MeV\leq\Lambda\leq 1200MeV$ and these can’t
be found in $\Lambda\geq 1200MeV$.

###  3.2 The position $\phi(0)=cte$

Now, we considered $\phi(0)=cte$, thus $M(r)$ will be as follows[19]:

$\displaystyle M(r)=$ $\displaystyle
a_{1}r+a_{2}r^{2}+a_{3}r^{3}+a_{4}r^{4}+a_{5}r^{5}+a_{6}r^{6}$ (16)
$\displaystyle+a_{7}r^{7}+a_{8}r^{8}+a_{9}r^{9}.$

Here, replacing Eq. (15) into the Eq. (10) and comparing to Eq. (7), a similar
equation is obtained by Eq. (12). In this position, i.e. $\phi(0)=cte$, the
value $n=0$ could be acceptable and $N(r)=1$ will be. Therefore, $\phi(r)$ ,
$\phi^{\prime\prime}(r)$ are obtained as follows:

$\displaystyle\phi(r)=e^{M(r)}$
$\displaystyle\phi^{\prime\prime}(r)=(M^{\prime\prime}+{M^{\prime}}^{2})e^{M}.$
(17)

The same as before, by replacing Eq.(15) in Eq. (16), 10 nonlinear equations
are obtained, we have:

$\displaystyle 2a_{2}+a_{1}^{2}=2\mu(C_{0}-E)$ $\displaystyle
6a_{3}+4a_{1}a_{2}=2\mu C_{1}$ $\displaystyle
12a_{4}+4a_{2}^{2}+6a_{1}a_{3}=2\mu C_{2}$ $\displaystyle
20a_{5}+8a_{1}a_{4}+12a_{2}a_{3}=2\mu C_{3}$ $\displaystyle
30a_{6}+9a_{3}^{2}+10a_{1}a_{5}+16a_{2}a_{4}=2\mu C_{4}$ $\displaystyle
42a_{7}+12a_{1}a_{6}+20a_{2}a_{5}+24a_{3}a_{4}=2\mu C_{5}$ $\displaystyle
56a_{8}+16a_{4}^{2}+14a_{1}a_{7}+24a_{2}a_{6}+30a_{3}a_{5}=2\mu C_{6}$
$\displaystyle
72a_{9}+16a_{1}a_{8}+28a_{2}a_{7}+36a_{3}a_{6}+40a_{4}a_{5}=2\mu C_{7}$
$\displaystyle
25a_{5}^{2}+18a_{1}a_{9}+32a_{2}a_{8}+42a_{3}a_{7}+48a_{4}a_{6}=2\mu C_{8}$
$\displaystyle 36a_{2}a_{9}+48a_{3}a_{8}+56a_{4}a_{7}+60a_{5}a_{6}=2\mu
C_{9}.$ (18)

Also, in this position for $P_{c}(4380)$, we are obtained the binding energy
of pentaquark and the numerical coefficients of wave function. In table 3, has
been shown the amount of binding energy from Eq. (17) and $\Lambda=800MeV$ for
$P_{c}(4380)$.

Table 3:  $(E_{B}(MeV),M_{P_{c}}(MeV))$ for $P_{c}(4380)$ and $\Lambda=800MeV$ $g_{1}$ | $g=0.51$ | $g=0.59$ | $g=0.67$
---|---|---|---
$0.75$ | $(-15.99,\ 4447.33)$ | $(-15.80,\ 4447.52)$ | $(-15.62,\ 4447.7)$
$0.94$ | - | $(-15.47,\ 4447.85)$ | $(-15.27,\ 4448.05)$
$1.95$ | - | $(-14.01,\ 4449.32)$ | -

To confirm the existence of pentaquark states, binding energy must be negative
i.e. $E_{B}<0$. Also, the total mass of particles contributing of
pentaquark(i.e. the sum of baryon $\Sigma_{c}$ and meson $\overline{D}^{*}$
masses) in addition to binding energy is closer to the mass of pentaquark
$P_{c}(4380)$. Here, the obtained results indicate binding energy ranging
$-102\leq E_{B}\leq-0.11$ for pentaquark $P_{c}(4380)$ which are negative.
Also, that conform to calculated results in the literatures[2, 20, 21, 22,
23]. According to mentioned above, obtained results are acceptable to a great
extent and they could be considered a clear evidence for the existence of a
bound five-quark state.

Figure 2, demonstrates the wave function’s diagrams for pentaquark $P_{c}$ in
$\Lambda=800MeV$, at different values of $g$ and $g_{1}$. These charts tend to
zero at the given value. As shown in the graphs, wave functions become zero in
$x\sim 0.12-0.13$, indicating that the maximum pentaquark radius is ranges
from 23.67 to 25.64 fm.

Figure 2: $P_{c}(4380)$ at $r$ in $\Lambda=800MeV$ for different values of $g$
and $g_{1}$

##  4 CONCLUSIONS

In this article, pentaquark $P_{c}(4380)$ system consisting of baryon
$\Sigma_{c}$ and $\overline{D}^{*}$ meson has been considered. The obtained
potential for pentaquark in reference[2] was presented and expanded. Then,
expanded potential was replaced in the Schrödinger equation and that was
solved as a bound state of two-body systems. By solving this to analytically
approach and according to the values of constants and cutoff, 10 nonlinear
differential equations and binding energy $E_{B}$ of pentaquark $P_{c}$ and
wave function coefficients were obtained. Results $E_{B}$ and wave function
were presented as tables and diagrams in the previous section, which could be
confirmed the existence of a bound state of pentaquark $P_{c}(4380)$. Then, it
specified that the wave function plots tend to be zero at a given value.
Therefore, the maximum radius of pentaquark $P_{c}$ were found out which
ranged from $x=23.67fm$ to $x=25.64fm$. We observed that the calculated values
are matches good with others and the mass of $P_{c}(4380)$. Also, we can say
that advantages of this study are that results obtained analytically while
other references are numerically calculated. Hence, the results are more
comprehensive and complete than other articles.

##  5 Data Availability

The authors confirm that the data supporting the findings of this study are
available within the article and its supplementary materials.

##  6 Disclosure

The research was performed as part of the employment of the authors.

##  7 Conflicts of Interest

The authors also declare that there is no conflict of interests regarding the
publication of this paper.

## References

  * [1] M. Gell-Mann; Phys. Lett. 8, 214 (1964).
  * [2] R. Chen, L. Xiang, L. Xue-Qian and Z. Shi-Lin; Phys. Rev. lett. 115, No. 13, (2015) 132002.
  * [3] J. He; Physics Letters B 753 (2016): 547-551.
  * [4] M. Karliner and J. L. Rosner ; Phys. rev. lett, 115, No.12, 122001(2015)[arXiv:1506.06386].
  * [5] J. J. Wu, R. Molina, E. Oset, and B. S. Zou; Phys. Rev. Lett. 105, 232001 (2010).
  * [6] Z. C. Yang, Z. F. Sun, J. He, X. Liu, and S. L. Zhu; Chin. Phys. C 36, 6 (2012).
  * [7] T. Uchino, W. H. Liang, and E. Oset; arXiv:1504.05726.
  * [8] LEPS Collaboration (T. Nakano et al.); Phys. Rev. Lett. 91, (2003) 012002.
  * [9] K. H. Hicks; Prog. Part. Nucl. Phys. 55, 647 (2005).
  * [10] T. Liu, Y. Mao, and B. Q. Ma; Int. J. Mod. Phys. A 29,1430020 (2014).
  * [11] C. Gignoux, B. Silvestre-Brac, and J. M. Richard; Phys. Lett. B 193, 323 (1987).
  * [12] D. O. Riska and N. N. Scoccola; Phys. Lett. B 299, 338 (1993).
  * [13] X. Q. Li and X. Liu; Eur. Phys. J. C 74, 3198 (2014).
  * [14] C. Isola, M. Ladisa, G. Nardulli, and P. Santorelli; Phys. Rev. D 68, 114001 (2003).
  * [15] X. Liu, Y.-R. Liu, W.-Z. Deng, and S.-L. Zhu; Phys. Rev. D 77, 094015 (2008).
  * [16] W. L. Wang, F. Huang, Z. Y. Zhang and B. S. Zou; Phys. Rev. C, 84, No.1, 015203(2011).
  * [17] F. Chezani Sharahi, M. Monemzadeh, A. Abdoli Arani; Mod. Phys. Lett. A, No.1, 1950220(2019).
  * [18] F. Chezani Sharahi, M. Monemzadeh, A. Abdoli Arani; jrmbs 9, No. 4, 10.22055(2020).
  * [19] E. Kreyszig, Advanced engineering mathematics, John Wiley and Sons (2010).
  * [20] Y.Shimizu, S.Daiki and H. Masayasu; Phys. Rev. D 93, No. 11, 114003 (2016).
  * [21] M. Monemzadeh, N. Tazimi, S. Babaghodrat; Advances in High Energy Physics, (2016).
  * [22] M. Radin, S. Babaghodrat, M. Monemzadeh; Physical Review D 90, No. 4, 047701 (2014).
  * [23] M. Monemzadeh, M. Hadizadeh, N. Tazimi; International Journal of Theoretical Physics 50, No. 3, 737-743 (2011).