# On-the-fly Tailoring towards a Rational Ansatz Design for Digital Quantum
Simulations

Dibyendu Mondal, Sonaldeep Halder, Dipanjali Halder, Rahul Maitra
<EMAIL_ADDRESS>Department of Chemistry,
Indian Institute of Technology Bombay,
Powai, Mumbai 400076, India

###### Abstract

Recent advancements in quantum information and quantum technology has
stimulated a good deal of interest in the development of quantum algorithms
for energetics and properties of many-fermionic systems. While the variational
quantum eigensolver is the most optimal algorithm in the Noisy Intermediate
Scale Quantum era, it is imperative to develop low depth quantum circuits that
are physically realizable in quantum devices. Within the unitary coupled
cluster framework, we develop COMPASS, a disentangled ansatz construction
protocol that can dynamically tailor an optimal ansatz using the one and two-
body cluster operators and a selection of rank-two scatterers. The
construction of the ansatz may potentially be performed in parallel quantum
architecture through energy sorting and operator commutativity prescreening.
With significant reduction in the circuit depth towards the simulation of
molecular strong correlation, COMPASS is shown to be highly accurate and
resilient to the noisy circumstances of the near-term quantum hardware.

##  I Introduction

The determination of atomic and molecular energetics and properties is one of
the most anticipated applications of quantum computing in near term quantum
computers. With the advent of coherent multiple qubit quantum processors, the
problem with the exponential growth of the Hilbert space that one often
encounters in many-body physics and chemistry can be handled in a tractable
manner. While the limited coherence time and poor gate fidelity are major
impediments towards the implementation of quantum many-body methods in Noisy
Intermediate Scale Quantum (NISQ) devices, nonetheless, one cannot deny the
importance of designing methods and algorithms that are realizable on quantum
computers to span many-body Hilbert spaceMcArdle  _et al._ (2020); Ortiz  _et
al._ (2001); Tilly  _et al._ (2022).

Historically, the quantum phase estimation algorithm (QPEA), proposed by
Abrams and Lloyd Abrams and Lloyd (1997, 1999), was the first quantum
algorithm towards the simulation of state energies of a many-fermionic system.
QPEA relies on the construction of an input reference state that is projected
onto the target eigenstate via a unitary operation Halder  _et al._ (2023);
Wang  _et al._ (2008); Aspuru-Guzik  _et al._ (2005). Although theoretically
appealing, the long coherence time required for the such unitary evolution due
to extremely deep and complex quantum circuits warrants the availability of
sufficiently large number of fault-tolerant qubits which is still beyond
experimental realization. The variational quantum eigensolverPeruzzo  _et al._
(2014) provides a lucrative alternative platform for the digital simulation of
quantum many-body systems in NISQ devices.

VQE has been realized experimentally in various quantum hardware architectures
like the photonic quantum processors Peruzzo  _et al._ (2014), superconducting
quantum processorsColless  _et al._ (2018) and trapped ion architecturesHempel
_et al._ (2018). Significant theoretical developments of VQE tailored for
various hardware platforms immediately followed its experimental
demonstration. VQE presumes the prior knowledge of the many-body Hamiltonian
matrix elements from the classically computed values and aims to finding the
"best" variational approximation to the ground state wavefunction of the
Hamiltonian starting from a trial wavefunction Ansatz. The trial wavefunction
is constructed in terms of a series of tunable parameters which provide its
variational flexibility. The state of the system is prepared through a
parametrized quantum circuit, followed by computing the expectation value of
the Hamiltonian through repeated measurements. Both these steps are carried
out in quantum hardwares while the update of the parameters take place via
some classical optimization in classical hardware. The attractive feature of
VQE is that the state preparation can be done via shorter quantum circuits
compared to QPEA, making it a more desirable candidate for NISQ devices.
Needless to say, the VQE results and the associated energy landscape often
depend critically on the choice of the parametrized ansatz Metcalf  _et al._
(2020); Sugisaki  _et al._ (2022); Tilly  _et al._ (2022).

One of the most significant developments towards the practical realization of
the unitary evolution in a quantum computer directly originated from the
unitary variant of coupled cluster theory (UCC) Shen  _et al._ (2017); Romero
_et al._ (2018); Evangelista, Chan, and Scuseria (2019); Sokolov  _et al._
(2020); Anand  _et al._ (2022); Peruzzo  _et al._ (2014) with a given
truncation in the rank of the cluster expansion. While the UCC ansatz with
single and double excitation operators (UCCSD) provides a compact wavefunction
ansatz, they work best in the cases where the molecular non-dynamic
correlation is somewhat weak. Several efforts have been made to make the UCC
work for areas of strong correlation, either by increasing the excitation rank
at the expense of higher implementation overhead or by including generalized
hamiltonian matrix elements in the wavefunction parametrization. The latter
class of methods stems from Nakatsuji’s density theoremNakatsuji (1976) and
Nooijen’s conjectureNooijen (2000) about the possibility of expressing the
eigenstates of a Hamiltonian containing one and two body terms in terms of
generalized excitation operators. While theoretically exact, a direct VQE
implementation with all generalizedLee  _et al._ (2018) excitation operators
incurs high implementation cost which is far beyond NISQ realization. The
excitation list can further be pruned by keeping only generalized paired
terms, giving rise to k-UpCCGSD ansatzLee  _et al._ (2018) for which the
circuit depth grows only linearly with the system size. In a different
approach, the present authors proposed a partially disentangled form of the
unitary in which the high rank excitations are implicitly folded in through
nested commutators of two-body scattering operators with effective hole-
particle excitation rank of one and the standard cluster operatorsHalder,
Prasannaa, and Maitra (2022).

The goal for an expressive ansatz is to span the $N-$electron Hilbert space
through a set of connected operators, ideally of rank $N$ at max. While
keeping the maximum rank of the operators (that enter the parametrization) two
for better NISQ realization, one must be able to express any $N-$body
connected hole-particle excitation through low power cumulative actions of
different non-commuting operatorsMaitra, Akinaga, and Nakajima (2017);
Tribedi, Chakraborty, and Maitra (2020). The action of any generalized
operator with effective excitation rank zero on functions of the so-called
primary excitation subspace (spanned by the set of $n-$tuply excited
determinants; vide infra) either leads to nilpotent solution, or results in
redundant functions within the primary excitation subspace. The two-body
operator with effective excitation rank one (to be referred as scattering
operator or scatterers, henceforth) can of course recursively lead to
$(n+1),(n+2),(n+3)...)-$tuply excited functions to span the secondary
excitation subspace. However, for such non-trivial action to exist, the
specific scatterer (that leads to a determinant of the secondary excitation
subspace by its action on a given primary excited determinant) must not
commute with the cluster operator which generates the primary excitation
function. As such, in this work, we pre-screen the scatterers on-the-fly in
terms of their (non-)commutativity with the cluster operators to come up with
a dynamically optimal solution protocol: the COMmutativity Pre-screened
Automated Selection of Scatterers (COMPASS). COMPASS can be implemented in a
parallel quantum architecture to choose the energetically most significant
cluster operators to span the primary excitation manifold and select the
scatterers from an operator bath based on the commutativity criteria to span
the secondary excitation manifold. The entire ansatz may be constructed in a
parallel quantum environment and can be implemented with extremely shallow
quantum circuit.

Following the discussions on the genesis of COMPASS, we would discuss its
performance for the cases of molecular strong correlation in which
conventional UCCSD fails to achieve desirable accuracy. In particular, we
would study the accuracy of COMPASS in handling molecular strong correlation,
vis-a-vis the number of parameters needed. Furthermore, we would analyse the
accuracy of COMPASS in the case of noisy simulations to demonstrate its
expected performance in faulty quantum devices and would argue its suitability
as a leading candidate for NISQ realization.

##  II Theory

###  II.1 Choice of the Operators Class and Motivation towards the Genesis of
COMPASS:

In order to design a compact parametrized ansatz, the choice of the operators
play the pivotal role. As discussed in the introduction, we choose to work
with a set of one and two-body cluster operators along with a set of
scatterers. The scatterers are two-body operators with effective hole-particle
excitation rank one and have one quasi-hole or quasi-particle "destruction
operator". Each element of the scatterers is equipped with an in-built
projector such that it has non-trivial action on only a selected set of
primary excitation subspace determinants:

$\displaystyle
S_{h}=\frac{1}{2}\sum_{amij}s^{am}_{ij}\\{a^{\dagger}m^{\dagger}ji\\}|\cdot\cdot
i\cdot\cdot j\cdot\cdot\cancel{m}\cdot\cdot\cancel{a}\cdot\cdot\rangle$
$\displaystyle\langle\cdot\cdot\cancel{a}\cdot\cdot\cancel{m}\cdot\cdot
j\cdot\cdot i\cdot\cdot|$ (1) $\displaystyle
S_{p}=\frac{1}{2}\sum_{abie}s^{ab}_{ie}\\{a^{\dagger}b^{\dagger}ei\\}|\cdot\cdot
i\cdot\cdot e\cdot\cdot\cancel{a}\cdot\cdot\cancel{b}\cdot\cdot\rangle$
$\displaystyle\langle\cdot\cdot\cancel{b}\cdot\cdot\cancel{a}\cdot\cdot
e\cdot\cdot i\cdot\cdot|$ (2)

Here, a, b, c, …, etc. denote the set of unoccupied particle orbitals and i,
j, k, …, etc. are the set of occupied hole orbitals with respect to the HF
vacuum. Note that, ‘m’ is a hole state and ‘e’ is a particle state with
respect to HF vacuum and they together form a contractible set of orbitals
(CSO). The orbitals constituting CSO appear as the quasi-hole and quasi-
particle destruction operators in $S$. Thus, the operator $S$ that contains
the label $m$ (or $e$) may be denoted as $S^{m}$ (or $S_{e}$). Similarly, the
cluster operator $T$ with orbital label $m$ (or $e$) may be denoted as $T_{m}$
(or $S^{e}$). Note that for the cluster operators, both $m$ ($e$) appears as
quasi-hole (quasi-particle) creation operators whereas they appear as
destruction operators in $S$. As such, the contractions between $S$ and $T$
take place through them to simulate a higher excitation rank operator. Each
contraction between the scatterers and the cluster operators increases the
effective hole-particle excitation rank by one, for example,
${\contraction{}{S}{}{T_{2}}ST_{2}}\rightarrow T_{3}$,
$\contraction{}{S}{}{S}\contraction[2ex]{}{S\;}{T_{2}}{}S\;ST_{2}\rightarrow
T_{4}$… . The anti hermitian counterpart of the $S$ and $T$ operators
($\sigma=S-S^{\dagger};\tau=T-T^{\dagger}$) may be used to construct an
unitary evolution operator which can be implemented in quantum architecture.
With the knowledge of the (non-)commutativity among various
$\sigma^{m}=S^{m}-{S^{m}}^{\dagger}$ and $\tau_{m}=T_{m}-{T_{m}}^{\dagger}$
(and among $\sigma_{e}=S_{e}-{S_{e}}^{\dagger}$ and
$\tau^{e}=T^{e}-{T^{e}}^{\dagger}$), the effects of connected high rank
excitations can be built through nested commutators like [$\sigma,\tau_{2}$]
$\rightarrow\tau_{3}$; [[$\sigma,\tau_{2}$],$\tau_{2}$]$\rightarrow\tau_{4}$…
and so on. As such a partially disentangled (factorized) ansatz can be shown
to include the infinite commutators even when a finite Trotter order Das and
Chakrabarti (2005) is used to approximate $e^{\tau}$ and/or $e^{\sigma}$.

$e^{\sigma}.e^{\tau}=e^{\sigma+\tau+[\sigma,\tau]+[[\sigma,\tau],\tau]+...}$
(3)

Note that the action of the unitary $e^{\sigma}$ on the
$e^{\tau}|\Phi_{HF}\rangle$ is partially nilpotent due to the inbuilt
projector in the definition of $\sigma$. Expanding the entangled states
$e^{\tau}|\Phi_{HF}\rangle$ in terms of constituent zero, one, two,… body
excited determinants (which belong to the primary excitation manifold), one
may write

$e^{\sigma}\Big{(}e^{\tau}\ket{\Phi_{HF}}\Big{)}\rightarrow
e^{\sigma}\Big{(}|\Phi_{HF}\rangle+\sum c_{I}|\Phi_{I}\rangle\Big{)}$ (4)

Here $\Phi_{I}$’s are the various determinants belonging to the primary
excitation manifold generated by $e^{\tau}\ket{\Phi_{HF}}$. The index $I$
would generically denote the composite hole-particle indices associated with
the excited determinants or with the cluster operators. Note that each
$\sigma$ operator is characterized by a quasi-hole/particle destruction
operator that belongs to CSO. The action of $\sigma$ is nontrivial only on
certain set of primary subspace determinants. For example, a $\sigma$ operator
with quasi-hole destruction orbital ’$m$’ acts only on those excited primary
subspace determinants where the occupancy of ’$m$’, $n_{m}$ is zero. For each
$\sigma^{m}$, the occupancy of ’$m$’ may be used to divide the primary
subspace determinants on which $\sigma^{m}$ acts into two sets.

$\sum_{I}^{N}c_{I}|\Phi_{I}\rangle=\sum_{J=1}^{N_{1}}c_{J}|\Phi_{J}(n_{m}=0)\rangle+\sum_{K=1}^{N_{2}}c_{K}|\Phi_{K}(n_{m}=1)\rangle$
(5)

with $N=N_{1}+N_{2}$ represents the total number of determinants generated by
$e^{\tau}|\Phi_{HF}\rangle$ that span the primary subspace. Note that for a
hole state ’$m$’, $e^{\sigma^{m}}(\sum_{J}c_{J}|\Phi_{J}(n_{m}=0)\rangle)$ has
non-vanishing action while
$e^{\sigma^{m}}(\sum_{K}c_{K}|\Phi_{K}(n_{m}=1)\rangle)$ is nilpotent. These
determinants with $n_{m}=0$ are principally generated by $\tau_{m}$ and these
cluster operators do not commute with $\sigma^{m}$. Unfortunately, an
unrestricted construction of the entangled state
$\sum_{I}c_{I}|\Phi_{I}(n_{m}=0\oplus n_{m}=1)\rangle)$ utilizes high number
of parameters, resulting in a deep circuit. An exact similar analysis can be
done when the destruction operator is a quasi-particle state that belongs to
CSO.

Instead, one may selectively generate only those primary subspace determinants
by the cluster operators on which the scatterers have a non-trivial action.
The COMPASS dynamically chooses the "best" set of cluster amplitudes in a
factorized manner. Each of the primary subspace determinants, generated by the
individual cluster amplitudes, are scanned for whether a non-commuting
scatterer (with which the corresponding cluster operator share CSO) have
significant effect or not. Where the commutativity criteria is met, the
scatterers are immediately allowed to act upon the entangled state. In the
next section, we would present the genesis of COMPASS and discuss upon how it
directs to the construction of the most optimal ansatz.

###  II.2 Development of COMPASS: an Automated Toolkit for Dynamic Ansatz
Design

Figure 1: Schematic representation of COMPASS protocol. In step 1(a-c), each
cluster operator $\tau$ is screened (potentially in a massively parallel
quantum architecture) through a threshold criteria and kept in an operator
block as per their descending contribution to correlation energy. Each such
selected $\tau$ is checked if that contains any orbital(s) that belongs to CSO
(step 2(a)). For each operator block with selected $\tau$ operator,
appropriate non-commuting $\sigma$ operator(s) are fetched from the operator
bath. The blocks are expanded by including the energetically dominant $\sigma$
operators (denoted in this figure as $\sigma^{\prime}$), one at a time,
selected via the optimization of two-parameter energy functionals (step 2(b)).
This step can again be performed in a parallel quantum architecture. The
operator bath, in principle, may supply infinite number of $\sigma$ operators
and is never drained out. The final ansatz is constructed by concatenating the
various operator blocks, keeping the ordering unchanged (step 3(a)).

COMPASS relies on the choice of the "best" set of the cluster operators and
selection of the appropriate scatterers that have significant contribution by
their action on certain entangled states. The whole ansatz is constructed
dynamically, possibly in a parallel quantum architecture and the resulting
ansatz features as a disentangled product of various $e^{\sigma_{\mu}}$ and
$e^{\tau_{I}}$ in an interwoven manner. As mentioned before, the indices
${I,J,K,L,...}$ are used to denote the composite excitation labels associated
with $\tau$ whereas ${\mu,\nu,\lambda}$ would denote the composite orbital
labels for $\sigma$. COMPASS consists of three components: a parallel
selection of the "best" amplitudes and energy ordering, selection of the
appropriate scatterers from an operator bath via CSO scanning, and the
preparation of the final ansatz. A schematic figure of COMPASS is presented in
Fig 1 and the details of the same is discussed below.

####  II.2.1 Choice of the "best" cluster operators and the ordering of their
appearance:

The choice of the "most important" cluster operators is a key step to design a
shallow depth quantum circuit in the NISQ devices Fan  _et al._ (2021).
Towards this, we have considered only the double excitation operators
(generically to be denoted as $\tau_{d}$) while the single excitation cluster
operators (generically to be denoted as $\tau_{s}$) will be handled
separately. With the availability of multiple quantum devices, one may
parallelly optimize the various one parameter energy functional:

$\displaystyle E_{I}$ $\displaystyle=$
$\displaystyle\bra{\Psi(\theta_{I})}\hat{H}\ket{\Psi(\theta_{I})},$
$\displaystyle=$
$\displaystyle\bra{\Psi_{HF}}e^{-\tau_{I}}He^{\tau_{I}}\ket{\Psi_{HF}}$
$\displaystyle\forall$ $\displaystyle I\in(1,n_{o}^{2}n_{v}^{2})$

Each of these can be done through a one parameter circuit, with $\theta_{I}$
being the sole variational parameter. Here
$\ket{\Psi(\theta_{I})}=e^{\tau_{I}(\theta_{I})}\ket{\Psi_{HF}}$. If
sufficiently large number of quantum computers are available, each $E_{I}$ can
be evaluated on different devices simultaneously. Otherwise, on a single
quantum device one may repeat this multiple times and take an average of
$E_{I}$. Only those $\tau_{I}$ operators will be kept in the cluster operator
pool for which $\Delta E_{I}=|E_{I}-E_{HF}|>\epsilon_{1}$, where
$\epsilon_{1}$ is a predefined threshold and the rest of the cluster operators
are discarded.

With the $N_{T}$ number of cluster operators that pass through the energy
screening, we align them in a descending order of their contribution to
correlation energy. Energetically most contributing cluster operator is
allowed to act on the HF state directly, followed by the next and so on. This
implies that for $N_{T}$ cluster operators with $...\Delta E_{J}>...>\Delta
E_{K}>...>\Delta E_{L}...$, $e^{\tau_{J}}$ is placed in operator block
$\beta$, $e^{\tau_{K}}$ in operator block $\gamma$, … and so on in a
disentangled manner. Note that $...\delta>...>\gamma>...>\beta>...>1$

$\displaystyle|\Psi\rangle$ $\displaystyle=$
$\displaystyle...\Big{[}e^{\tau_{L}}\Big{]}_{\delta}...\Big{[}e^{\tau_{K}}\Big{]}_{\gamma}...\Big{[}e^{\tau_{J}}\Big{]}_{\beta}...|\Psi_{HF}\rangle$
$\displaystyle|\Psi\rangle$ $\displaystyle=$
$\displaystyle\prod_{\alpha=1}^{N_{T}}\Big{[}e^{\tau_{I}}\Big{]}_{\alpha}|\Psi_{HF}\rangle$
(7)

where the quantities inside parenthesis denote the order of the operator
blocks in which they act on $|\Psi_{HF}\rangle$. This ordering will be
maintained throughout. In general, if a given $\tau_{I}$ is placed in operator
block $[..]_{\alpha}$, we would denote it as $[\tau_{I}]_{\alpha}$. We mention
that each cluster operator is allowed to appear only once irrespective of the
block it is placed. In the next subsection, we will dynamically expand the
operator blocks by placing the scatterers appropriately based on the
commonality of CSO labels shared between the cluster operators and the
scatterers.

####  II.2.2 Selection of the Scatterers from Operator Bath

For each operator $[\tau_{I}]_{\alpha}$ selected through the energy screening,
one first checks if $\tau_{I}$ contains any orbital that belongs to CSO. For
the cluster operators which do not contain any index of CSO, the corresponding
operator block is not expanded any further. For all other cluster operators
with one or more orbitals that belongs to CSO, a pool of scatterers is created
from the scatterer bath through commutativity screening. This implies that the
cluster operator and the scatterers should contain the same set of
contractible orbital(s). Thus for a given $[{\tau_{{m}_{I}}}]_{\alpha}$, the
scatterer pool contains all the operators of the structure like $\sigma^{m}$.
With $N_{\alpha}$ number of such scatterer selected in the pool (corresponding
to the operator block $\alpha$ containing $\tau_{I}$), one optimizes the
following two parameter energy functional:

$\displaystyle E_{I\mu}$ $\displaystyle=$
$\displaystyle\bra{\Psi(\theta_{I},\theta_{\mu})}\hat{H}\ket{\Psi(\theta_{I},\theta_{\mu})},$
$\displaystyle=$
$\displaystyle\bra{\Psi_{HF}}\overline{e^{-\tau_{I}}e^{-\sigma_{\mu}}}H\overline{e^{\sigma_{\mu}}e^{\tau_{I}}}\ket{\Psi_{HF}}$
$\displaystyle\forall$ $\displaystyle\mu\in(1,N_{\alpha})$

Here the overline suggests the connected action of $\tau_{I}$ and the various
$\sigma_{\mu}$’s that share at least one common orbital index that belongs to
CSO. However, only those scatterers are screened for which
$|E_{I\mu}-E_{I}|>\epsilon_{2}$, where $\epsilon_{2}$ is a predefined
threshold whose value is usually taken to be order of magnitude less than
$\epsilon_{1}$. With the energy condition met, the operator block
$[...]_{\alpha}$ is expanded by including appropriate $\sigma_{\mu}$’s through
disentangled (factorized) unitary.

$\displaystyle U_{\alpha}$ $\displaystyle=$
$\displaystyle\Big{[}\prod_{\mu}\overline{e^{\sigma_{\mu}}e^{\tau_{I}}}\Big{]}_{\alpha}$
(9)

There may, of course, be specific cases where the $\tau$ operator in a given
operator block does not contain any orbital of CSO and in that case, as
mentioned before, no $\sigma$ gets attached to it. However, it still gets
placed according to its energy ordering. Note that, for general cases, with
each $\sigma_{\mu}$ being selected one-by-one for a given $\tau_{I}$, the
optimization in Eq. II.2.2 can be done with shallow two-parameter circuit with
$\theta_{I},\theta_{\mu}$ as the two variational parameters. One may further
note that the selection of scatterers for various operator blocks are
independent of each other and thus the evaluation of $E_{I\mu}$ in Eq. II.2.2
for various pairs of $\tau_{I}$ and $\sigma_{\mu}$ can be done with multi-
level parallelization. Furthermore, if sufficiently large number of quantum
computers are available, one may evaluate the energy functional $E_{I\mu}$ in
different devices simultaneously or otherwise, one may perform it multiple
times in a single quantum device and take an average value. We further note
that while a given $\tau_{I}$ appears only once in the ansatz, a given
$\sigma_{\mu}$ may be attached to various $\tau_{I}$’s in their respective
operator blocks and thus the $\sigma$ operators are never drained out from the
operator bath.

####  II.2.3 Weaving the Final Ansatz

The final ansatz is constructed by concatenating the various operator blocks.
As mentioned above, depending on the commonality of the contractible orbitals
between the $\tau$ and $\sigma$ operators, several $\tau_{I}$ operator may
attach same $\sigma$ operator to appear in their respective operator blocks.
Finally, all the one-body cluster operators, $\tau_{s}$, are placed at the end
in lexical ordering.

$U=\prod_{s}e^{\tau_{s}}\prod_{\alpha}U_{\alpha}$ (10)

Once the full ansatz is constructed for a given molecule with a fixed nuclear
arrangement, the parameters involved in the ansatz, Eq. 10 are optimized using
the standard VQE hybrid quantum-classical framework till all the parameters
are converged.

One may note that the proposed ansatz may span the entire $N-$electron Hilbert
space. However, in COMPASS each $\sigma$ operator is chosen and clubbed with
one or more non-commuting $\tau$ operator(s) in a way that lowers the
correlation energy beyond a certain predefined threshold. The selected
$\sigma$ operators are allowed to act immediately upon certain entangled
states that are generated by the action of the associated non-commuting $\tau$
on the Hartree Fock determinant. Note that only those $\sigma$ operators are
kept in a given operator block (having a pre-selected cluster operator $\tau$)
that can generate at least a triply excited secondary excitation subspace
function.

##  III Results:

COMPASS has been implemented with an interface to qiskit-nature Abraham et.
al. (2021) which imports the one and two-electron integrals from PySCF Sun
_et al._ (2018). All the calculations performed in this study employed STO-3G
basis set Hehre, Stewart, and Pople (1969) with a direct spinorbital to qubit
mapping. The Jordan-Wigner transformation Seeley, Richard, and Love (2012) was
employed to encode second quantized operators to qubit operators. For all our
calculations, we chose the L-BFGS-B Morales (2002); Byrd  _et al._ (1995)
optimizer for the classical optimization. Also, we had initialized the qubits
to the Hartree-Fock reference state and each parameters of the ansatz was
initialized to the optimized values obtained from the minimization of energy
functional of Eq.II.2.2 (and Eq. II.2.1).

While the $\sigma$ operators generically used so far towards the development
have all possible spin- and spatially-unrestricted terms in it, in the actual
implementation, the list is further significantly pruned by keeping only a few
specific kinds of the same. Two different cases are considered which are
different from each other in the choice of the $\sigma$ operator taken in the
operator bath. In particular, we have worked with the (a) opposite spin (OP)
sector of $\sigma$ and (b) the partially paired (PP) sector of $\sigma$.

The OP sector incorporates the specific low-spin channel of $\sigma$. This
implies that the spins in the excitation vertex and the scattering vertex of
are different.

$\displaystyle\sigma_{OP}\in\\{{(\sigma_{h})}_{i_{\alpha}j_{\beta}}^{a_{\alpha}u_{\beta}},{(\sigma_{h})}_{i_{\beta}j_{\alpha}}^{a_{\beta}u_{\alpha}};{(\sigma_{p})}_{i_{\alpha}v_{\beta}}^{a_{\alpha}b_{\beta}},{(\sigma_{p})}_{i_{\beta}v_{\alpha}}^{a_{\beta}b_{\alpha}}\\}$
(11)

$u$ and $v$ are the spatial active hole and virtual orbitals that form the
CSO, and $\alpha$ and $\beta$ denote the spin-up and spin-down electrons,
respectively.

The PP sector of $\sigma$ includes (i) only those quasi-hole creation
operators which originate from the same spatial orbitals (for $\sigma_{h}$)
and (ii) only those quasi-particle creation operators which share the same
spatial orbitals (for $\sigma_{p}$).

$\displaystyle\sigma_{PP}\in\\{{(\sigma_{h})}_{i_{\alpha}i_{\beta}}^{a_{\alpha}u_{\beta}},{(\sigma_{h})}_{i_{\beta}i_{\alpha}}^{a_{\beta}u_{\alpha}};{(\sigma_{p})}_{i_{\alpha}v_{\beta}}^{a_{\alpha}a_{\beta}},{(\sigma_{p})}_{i_{\beta}v_{\alpha}}^{a_{\beta}a_{\alpha}}\\}$
(12)

Irrespective of the choice of the sector of $\sigma$, the orbitals
constituting CSO are restricted only to those that span the chemically active
region.

Figure 2: Accuracy as a function of the bond length parameter for various
versions of COMPASS with respect to FCI: (a) $BH$, (b) $H_{2}O$, and (c)
linear $BeH_{2}$. The region shaded in pale blue indicates chemical accuracy
in kcal/mol. (d), (e) and (f) estimate parameter counts for $BH$, $H_{2}O$ and
linear $BeH_{2}$, respectively, along the potential energy profile. Figure 3:
NPE (white symbols) and average error (black symbols) (in kcal/mol) across the
potential energy profile for (a) $BH$, (b) $H_{2}O$, and (c) linear $BeH_{2}$
models. For $H_{2}O$, the NPE shows two orders of magnitude improvement over
the UCCSD while for the linear $BeH_{2}$, all variants of COMPASS have several
orders of magnitude improvement in NPE, estimating a uniform description of
many-body correlation effects over the entire energy profile. The color-coded
horizontal axis denotes the order of NPE and average errors.

As discussed previously, COMPASS uses two threshold parameters ($\epsilon_{1}$
and $\epsilon_{2}$) which ensure to include the most dominant cluster
amplitudes and the scatterers in the final ansatz. As such, COMPASS
calculations with thresholds of $\epsilon_{1}$ and $\epsilon_{2}$ would be
denoted as COMPASS($-log(\epsilon_{1}),-log(\epsilon_{2})$).

###  III.1 Accuracy vs Parameter Count over the Potential Energy Profile of
Strongly Correlated Molecules:

Towards the study of the accuracy of COMPASS, three difficult test cases have
been identified: the potential energy surface for the the stretching of $BH$
single bond, the symmetric simultaneous single bond stretching of $H_{2}O$ and
linear $BeH_{2}$. In all cases, the results are compared against UCCSD. In all
our calculations, all the amplitudes corresponding to the spin-complemented
operators are treated as independent.

Stretching of the single bond in $BH$ is one of the most difficult test cases
for assessing the accuracy of any quantum many-body theory. The system shows
the signature of strong molecular correlation due to the interplay of the
ground and excited roots when the $B--H$ bond is stretched. We model its
dissociation profile with COMPASS and UCCSD. BH consists of 6 electrons in 12
spinorbitals with a total Hilbert space dimension of 4096. The UCCSD energy
profile, when plotted as a function of $B--H$ distance shows energy error
$\sim 1$ kcal/mol with respect to the classically exact Full Configuration
Interaction (FCI) method throughout various molecular arrangements. COMPASS
with both the variants, OP and PP, outperform UCCSD throughout the energy
surface and particularly when the molecule is sufficiently stretched, it shows
improvement by an order of magnitude. Interestingly enough, COMPASS(5,7) with
OP variant takes about 67 parameters at bond length of 2.5 and 2.75Å, which is
the highest (compared to other nuclear arrangements) in terms of the number of
parameters required for COMPASS, but still this is an order of magnitude less
than the 117 parameters (without combining the spin-complementary excitations)
taken up by UCCSD. The accuracy throughout the potential energy profile is
also illustratively evident from the non-parallelity error (NPE) and the
average error for various schemes and is shown in Fig. 3(a).

Figure 4: Accuracy as a function of bond length parameters for $H_{2}O$
($H-O-H=104.4776^{\circ}$) under gaussian noise model characterized by
standard deviation $SD$ : (a) $SD=10^{-2}$, (b) $SD=10^{-3}$, (c)
$SD=10^{-4}$. For both (b) and (c) the average energy predicted by COMPASS are
well within chemical accuracy. For all the variants, COMPASS predicts energy
which are at least an order of magnitude accurate than UCCSD across the energy
profile; demonstrating its resilience to noisy environment. The vertical lines
across each point suggest the standard deviation of the observed energy values
with 100 independent samples.

$H_{2}O$ in STO-3G basis, with one core spatial orbital frozen, renders to be
a system with 8 electrons in 12 spinorbitals. In Fig. 2b, we have plotted the
energy error with respect to FCI method as a function of the bond length
parameter keeping the $H-O-H$ angle fixed at $104.4776^{\circ}$. While the
UCCSD ansatz requires 92 parameters (without combining the spin- complementary
excitations), it fails to achieve energy error within chemical accuracy when
the bonds are sufficiently stretched. Dramatically enough, all the variants of
COMPASS with different threshold parameters significantly outperforms UCCSD
and resulting COMPASS energy values are just a few micro-Hartree away from FCI
throughout the energy surface. Even more impressive is the significantly fewer
number of parameters selected by most of the COMPASS variants to achieve this
level of accuracy throughout the potential energy profile: COMPASS(5,6) with
the PP variant barely needs about 66 parameters, resulting in a much compact
form of the wavefunction ansatz. Beyond 1.4 Å, OP variants require higher
number of parameters than the PP counterpart; however, this also results in at
least an order of magnitude improvement in accuracy. The overall superiority
of COMPASS in comparison to UCCSD is again measured in terms of NPE and the
average error and is deomnstrated in Fig. 3(b).

The efficacy of COMPASS is much more pronounced in strongly correlated systems
like a linear $H-Be-H$ model when we simultaneously stretch the $Be-H$ bonds
in a symmetric manner. Due to the strong electronic correlation on symmetric
stretching of the bonds, the UCCSD behaves somewhat poorly and beyond
$R_{Be-H}=1.75$Å, its energy error with respect to FCI goes beyond 1 kcal/mol.
The error further increases on further stretching the bonds, marking a clear
signature of strong correlation. COMPASS with all the variants show remarkable
improvement over UCCSD results with at least one order of magnitude in the
region of short $R_{Be-H}$ length and by several orders of magnitude when the
bonds are further stretched. This comes at a remarkable reduction in the
number of parameter count over UCCSD: in the worst case scenario, COMPASS(5,7)
with OP variant utilises about 54 parameters around $R_{Be-H}=2.00-2.25$Å,
while UCCSD requires 92 parameters. In fact, the OP variant shows several
order of magnitude improvement in accuracy over the PP variant, particularly
in regions beyond 2 Å. However, for two different sets of $\epsilon_{1}$ and
$\epsilon_{2}$, both COMPASS-OP and COMPASS-PP require similar number of
parameters throughout the potential energy profile. The extremely high
accuracy over the potential energy profile is again measured in terms of the
NPE and average error (Fig. 3(c)) and COMPASS variants show clear two orders
of magnitude improvement of these metrics over UCCSD. The accuracy and
parameter count in linear $BeH_{2}$ model clearly demonstrate the performance
of COMPASS in compactifying the ansatz for the simulation of strong electronic
correlation, making it more suitable for NISQ realization.

###  III.2 Simulation with Gaussian Noise Model:

The simulation so far has presumed a noiseless implementation in an ideal
quantum environment. However, the NISQ devices are not fault-tolerant and as
such in practical scenario, one needs to account for the uncertainties due to
imperfect implementation. Towards this, for each set of optimal parameters
generated in an ideal quantum setting, we randomly produce 100 samples of
noisy parameters derived via a Gaussian distribution model Meitei  _et al._
(2021); Stair and Evangelista (2021). For each parameter, the distribution is
centered at the optimal parameter value and the noisy parameters
($\theta_{m}$) are generated as:

$\theta_{m}=exp\Big{(}-\frac{(\theta-\theta_{n})^{2}}{2(SD)^{2}}\Big{)}$ (13)

Here $\theta_{n}$ is the optimal value of the parameters and $SD$ is the
standard deviation. Three different cases are considered with $SD=10^{-2}$,
$SD=10^{-3}$ and $SD=10^{-4}$. With each of the 100 samples taken into account
for various $SD$, we estimate the average energy.

To estimate the performance of COMPASS with imprecise implementation of the
parameters, we choose to work with $H_{2}O$ with the same set of geometrical
parameters as in the case discussed earlier. We took a conservative COMPASS
parameters $\epsilon_{1}=10^{-5},\epsilon_{2}=10^{-6}$. In all various noise
strengths, the different variants of COMPASS seem to be more resilient than
UCCSD. With low noise strength ($SD=10^{-4},SD=10^{-3}$), in the regions where
the bonds are sufficiently stretched and the strong correlation dominates,
COMPASS variants are clearly a few orders of magnitude more accurate than the
UCCSD counterpart (Fig. 4b, c). Only in the case of reasonably strong noise
with $SD=10^{-2}$, COMPASS fails to achieve chemical accuracy though it is
about an order of magnitude more accurate than UCCSD. One may note that in
cases of such noisy simulations, the optimal parameters are of the same order
of magnitude as the noise itself, and this may not potentially cause a
hindrance to its realistic implementation.

##  IV Discussion and Future Outlook:

Being a dynamic structured ansatz, it is difficult to comment on the resource
requirements and the computational cost of COMPASS. The measure for the
resource requirement for a NISQ device is often taken to be the circuit depth
and the number of measurements. While the number of measurement is an
important parameter that determines the overall time to the solution, for a
NISQ realization, the circuit depth may be taken as the most important metric.
One may note that the latter determines the feasibility of implementation of
the algorithm in devices with short coherence time.

The parent VQE algorithm suffers from the limitation of large number of
measurements required to reach the solution. The COMPASS algorithm is likely
to suffer from the same drawback, if not slightly more. This mainly is due to
the (a) the ansatz construction (step 1,2, Fig. 1) and (b) optimization of the
final ansatz. One may note that the ansatz construction only requires several
one and two-parameter energy functional optimization. Each of such
optimization require very few number of iterations and hence they involve few
measurements. Most importantly, each of these energy optimizations are
independent and thus can be performed in parallel quantum architecture. Thus,
with the development of hardware capabilities, the requirement of increased
number of shot count during the ansatz construction is likely to pose limited
challenges. However, the amplitude optimization towards the final state
preparation may take somewhat more number of steps (than standard UCCSD) due
to the generalized operators present in the ansatz, resulting in a slower
convergence. This is principally attributed to the classical optimizer and one
may adopt various strategies to accelerate. The development of the fastest
convergence algorithm would be studied in near future as it requires several
aspects to take care of.

We must point out that COMPASS is conceptually entirely different than the
sequential growth ansatz like ADAPT-VQE Grimsley  _et al._ (2019) where the
operators are chosen based on the energy gradients with respect to individual
parameters from an operator pool. COMPASS, on the other hand, is a dynamically
expandable ansatz construction protocol where no such energy gradient needs to
be computed to select the operators; rather, one selects the "best" set of
parameters through energy optimization of one and two-parameter energy
functionals in parallel quantum architecture. The tailoring protocol
automatically ensures to include the most significant cluster amplitudes
according to energetically optimal ordering. Appropriate scatterers are
immediately allowed to act upon once an entangled state is generated through
the action of cluster operators on the Hartree Fock determinant. This has
shown to create a significantly compact and accurate wavefunction ansatz while
excitations of arbitrarily high rank can be captured through implicit
commutativity between various class of operators. The COMPASS protocol tailors
the ansatz depending on the system and its associated electronic complexity
under consideration, while keeping the parameter count to the minimum. The
compact representation of the wavefunction makes it a desirable candidate
towards its realization on a NISQ platform.

One major advantage of COMPASS is its ability to build the wavefunction ansatz
in a parallel architecture, making the ansatz construction significantly less
impacted by the noisy environment. The accuracy is largely controlled by the
two tuneable parameters. Even with somewhat conservative choices of these
parameters, the numerical results with various complex molecular systems amply
demonstrate its consistent superiority over UCCSD in terms of accuracy, gate
count and the ease of implementation. With high degree of quantum parallelism
towards the preparation of the ansatz and its extremely low execution gate
depth towards the state preparation through the parametrized circuit, COMPASS
would be an extremely desirable candidate for digital molecular simulations in
near-term quantum devices.

##  V Acknowledgement:

DM thanks Prime Minister’s Research Fellowship (PMRF), Government of India for
his research fellowship. SH thanks Council of Scientific and Industrial
Research (CSIR), Government of India, and DH thanks Industrial Research and
Consultancy Center (IRCC), IIT Bombay for their research fellowships.

## References

  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, “Quantum computational chemistry,” Rev. Mod. Phys. 92, 015003 (2020).
  * Ortiz  _et al._ (2001) G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, “Quantum algorithms for fermionic simulations,” Phys. Rev. A 64, 022319 (2001).
  * Tilly  _et al._ (2022) J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, and J. Tennyson, “The variational quantum eigensolver: A review of methods and best practices,” Physics Reports 986, 1–128 (2022).
  * Abrams and Lloyd (1997) D. S. Abrams and S. Lloyd, “Simulation of many-body fermi systems on a universal quantum computer,” Phys. Rev. Lett. 79, 2586–2589 (1997).
  * Abrams and Lloyd (1999) D. S. Abrams and S. Lloyd, “Quantum algorithm providing exponential speed increase for finding eigenvalues and eigenvectors,” Phys. Rev. Lett. 83, 5162–5165 (1999).
  * Halder  _et al._ (2023) D. Halder, V. S. Prasannaa, V. Agarawal, and R. Maitra, “Iterative quantum phase estimation with variationally prepared reference state,” International Journal of Quantum Chemistry 123, e27021 (2023), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.27021 .
  * Wang  _et al._ (2008) H. Wang, S. Kais, A. Aspuru-Guzik, and M. R. Hoffmann, “Quantum algorithm for obtaining the energy spectrum of molecular systems,” Phys. Chem. Chem. Phys. 10, 5388–5393 (2008).
  * Aspuru-Guzik  _et al._ (2005) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, “Simulated quantum computation of molecular energies,” Science 309, 1704–1707 (2005), https://www.science.org/doi/pdf/10.1126/science.1113479 .
  * 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, “A variational eigenvalue solver on a photonic quantum processor,” Nature Communications 5 (2014), 10.1038/ncomms5213.
  * 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, “Computation of molecular spectra on a quantum processor with an error-resilient algorithm,” Phys. Rev. X 8, 011021 (2018).
  * 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, “Quantum chemistry calculations on a trapped-ion quantum simulator,” Phys. Rev. X 8, 031022 (2018).
  * Metcalf  _et al._ (2020) M. Metcalf, N. P. Bauman, K. Kowalski, and W. A. de Jong, “Resource-efficient chemistry on quantum computers with the variational quantum eigensolver and the double unitary coupled-cluster approach,” Journal of Chemical Theory and Computation 16, 6165–6175 (2020), pMID: 32915568, https://doi.org/10.1021/acs.jctc.0c00421 .
  * Sugisaki  _et al._ (2022) K. Sugisaki, T. Kato, Y. Minato, K. Okuwaki, and Y. Mochizuki, “Variational quantum eigensolver simulations with the multireference unitary coupled cluster ansatz: a case study of the c2v quasi-reaction pathway of beryllium insertion into a h2 molecule,” Phys. Chem. Chem. Phys. 24, 8439–8452 (2022).
  * Shen  _et al._ (2017) Y. Shen, X. Zhang, S. Zhang, J.-N. Zhang, M.-H. Yung, and K. Kim, “Quantum implementation of the unitary coupled cluster for simulating molecular electronic structure,” Physical Review A 95 (2017), 10.1103/physreva.95.020501.
  * Romero  _et al._ (2018) J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, “Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz,” Quantum Science and Technology 4, 014008 (2018).
  * Evangelista, Chan, and Scuseria (2019) F. A. Evangelista, G. K.-L. Chan, and G. E. Scuseria, “Exact parameterization of fermionic wave functions via unitary coupled cluster theory,” The Journal of Chemical Physics 151, 244112 (2019), https://doi.org/10.1063/1.5133059 .
  * Sokolov  _et al._ (2020) I. O. Sokolov, P. K. Barkoutsos, P. J. Ollitrault, D. Greenberg, J. Rice, M. Pistoia, and I. Tavernelli, “Quantum orbital-optimized unitary coupled cluster methods in the strongly correlated regime: Can quantum algorithms outperform their classical equivalents?” The Journal of Chemical Physics 152, 124107 (2020), https://doi.org/10.1063/1.5141835 .
  * Anand  _et al._ (2022) A. Anand, P. Schleich, S. Alperin-Lea, P. W. K. Jensen, S. Sim, M. Díaz-Tinoco, J. S. Kottmann, M. Degroote, A. F. Izmaylov, and A. Aspuru-Guzik, “A quantum computing view on unitary coupled cluster theory,” Chem. Soc. Rev. 51, 1659–1684 (2022).
  * Nakatsuji (1976) H. Nakatsuji, “Equation for the direct determination of the density matrix,” Phys. Rev. A 14, 41–50 (1976).
  * Nooijen (2000) M. Nooijen, “Can the eigenstates of a many-body hamiltonian be represented exactly using a general two-body cluster expansion?” Phys. Rev. Lett. 84, 2108–2111 (2000).
  * Lee  _et al._ (2018) J. Lee, W. J. Huggins, M. Head-Gordon, and K. B. Whaley, “Generalized unitary coupled cluster wave functions for quantum computation,” Journal of Chemical Theory and Computation 15, 311–324 (2018).
  * Halder, Prasannaa, and Maitra (2022) D. Halder, V. S. Prasannaa, and R. Maitra, “Dual exponential coupled cluster theory: Unitary adaptation, implementation in the variational quantum eigensolver framework and pilot applications,” The Journal of Chemical Physics 157, 174117 (2022), https://doi.org/10.1063/5.0114688 .
  * Maitra, Akinaga, and Nakajima (2017) R. Maitra, Y. Akinaga, and T. Nakajima, “A coupled cluster theory with iterative inclusion of triple excitations and associated equation of motion formulation for excitation energy and ionization potential,” The Journal of Chemical Physics 147, 074103 (2017), https://doi.org/10.1063/1.4985916 .
  * Tribedi, Chakraborty, and Maitra (2020) S. Tribedi, A. Chakraborty, and R. Maitra, “Formulation of a dressed coupled-cluster method with implicit triple excitations and benchmark application to hydrogen-bonded systems,” Journal of Chemical Theory and Computation 16, 6317–6328 (2020), pMID: 32794747, https://doi.org/10.1021/acs.jctc.0c00736 .
  * Das and Chakrabarti (2005) A. Das and B. Chakrabarti, “Quantum annealing and related optimization methods,” Quantum Annealing and Related Optimization Methods, Edited by A. Das and B.K. Chakrabarti. 2005 XIV, 378 p. 124 illus. Also available online. ISBN 3-540-27987-3. Berlin: Springer, 2005. (2005).
  * Fan  _et al._ (2021) Y. Fan, C. Cao, X. Xu, Z. Li, D. Lv, and M.-H. Yung, “Circuit-depth reduction of unitary-coupled-cluster ansatz by energy sorting,” (2021).
  * Abraham et. al. (2021) H. Abraham et. al., “Qiskit: An open-source framework for quantum computing,” (2021).
  * Sun  _et al._ (2018) Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, “Pyscf: the python-based simulations of chemistry framework,” WIREs Computational Molecular Science 8, e1340 (2018), https://wires.onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1340 .
  * Hehre, Stewart, and Pople (1969) W. J. Hehre, R. F. Stewart, and J. A. Pople, “Self-consistent molecular-orbital methods. i. use of gaussian expansions of slater-type atomic orbitals,” The Journal of Chemical Physics 51, 2657–2664 (1969), https://doi.org/10.1063/1.1672392 .
  * Seeley, Richard, and Love (2012) J. T. Seeley, M. J. Richard, and P. J. Love, “The bravyi-kitaev transformation for quantum computation of electronic structure,” The Journal of Chemical Physics 137, 224109 (2012), https://doi.org/10.1063/1.4768229 .
  * Morales (2002) J. L. Morales, “A numerical study of limited memory bfgs methods,” Applied Mathematics Letters 15, 481–487 (2002).
  * Byrd  _et al._ (1995) R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, “A limited memory algorithm for bound constrained optimization,” SIAM Journal on Scientific Computing 16, 1190–1208 (1995), https://doi.org/10.1137/0916069 .
  * Meitei  _et al._ (2021) O. R. Meitei, B. T. Gard, G. S. Barron, D. P. Pappas, S. E. Economou, E. Barnes, and N. J. Mayhall, “Gate-free state preparation for fast variational quantum eigensolver simulations,” npj Quantum Information 7, 155 (2021).
  * Stair and Evangelista (2021) N. H. Stair and F. A. Evangelista, “Simulating many-body systems with a projective quantum eigensolver,” PRX Quantum 2, 030301 (2021).
  * Grimsley  _et al._ (2019) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, “An adaptive variational algorithm for exact molecular simulations on a quantum computer,” Nature Communications 10 (2019), 10.1038/s41467-019-10988-2.


# Radiative $E1$ transitions between ${}^{3}P_{1}$ and ${}^{3}S_{1}$
quarkonium states

Vaishali Guleria1, Eshete Gebrehana2, Shashank Bhatnagar1

###### Abstract

In this work we study the E1 decay processes, ${}^{3}P_{1}$ $\rightarrow$
${}^{3}S_{1}\gamma$, and ${}^{3}S_{1}$ $\rightarrow$ ${}^{3}P_{1}\gamma$ in
the framework of Bethe-Salpeter equation and calculate their decay widths. We
have used algebraic forms of Salpeter wave functions obtained through analytic
solutions of mass spectral equations for ground and excited states of
${}^{3}S_{1}$, and ${}^{3}P_{1}$ equal mass quarkonia in approximate harmonic
oscillator basis to do analytic calculations of their decay widths. These
decay widths have been compared with data and other models.

1\. Department of Physics, University Institute of Sciences, Chandigarh
University, Mohali-140413, India
2\. Department of Physics, Woldia University, Ethiopia

Key words: Bethe-Salpeter equation, radiative decays, decay widths, charmonia,

##  1\. Introduction

Radiative transitions offer valuable insights into the dynamics and properties
of quarkonia. Understanding the electromagnetic transitions in quarkonia is
crucial to the field of theoretical and experimental particle physics,
shedding light on the fundamental forces that govern the behavior of these
exotic particles. Further, among quarkonium decays, radiative decays are
particularly valuable as tests of various models, since the photon is directly
observed and the nature of the electromagnetic transition is well understood.
$E1$ transitions between $Q\bar{Q}$ states are the simplest radiative
transitions of the form $J^{++}\rightarrow 1^{--}\gamma$ (with $J=0,1,2$). And
since the E1 transitions are characterized by $|\Delta L|=1$, thus in these
transitions there is a change in parity between the initial and final hadronic
states. Electric dipole transitions are much stronger than magnetic dipole
transitions, and involve transitions between excited states.

The present work involves radiative E1 transitions involving ${}^{3}S_{1}$ and
${}^{3}P_{1}$ charmonium states. We wish to mention that though a lot of
investigation has been carried out on S-wave vector mesons such as, $0^{-+}$,
and $1^{--}$, however comparatively lesser investigation has been carried out
on P-wave mesons such as $1^{++},1^{+-}$ and $0^{++}$. Now, regarding the
$1^{++}$ mesons, they were first seen in pp collisions by R704 Collaboration.
However, not many decays of these mesons are experimentally observed as can be
checked from PDG tables [1, 2, 3]. In 2013, the Belle Collaboration[4]
reported measurements of $B\rightarrow\chi_{c1}\gamma K$ decays, and
$B\rightarrow\chi_{c2}\gamma K$ decays, where the $\chi_{c1}$ and $\chi_{c2}$
decay to $J/\Psi\gamma$. These results were obtained from a data sample of
$776\times 106$ $B\bar{B}$ events collected with the Belle detector at the
KEKB asymmetric-energy $e^{-}e^{+}$ collider operating at the $Y(4S)$
resonance.

Further, an indirect way of producing P-wave states is through $e^{-}e^{+}$
annihilation, which first produces ${}^{3}S_{1}$ charmonium states such as
$J/\Psi$, and $\Psi(2S)$, followed by their radiative M1 and E1 decays which
produce charmonium states, ${}^{1}S_{0}$ (such as $\eta_{c}$), and
${}^{3}P_{1}$ states (such as $\chi_{c1}$). We list here some of the indirect
ways of producing $1^{++}$ state: In [5], from $e^{-}e^{+}$ collision data
acquired with the CLEO detector at CESR, they observed the non-$D\bar{D}$
decay $\Psi(3770)\rightarrow\gamma\chi_{c1}$ with a statistical significance
of $6.6\sigma$, using the two-photon cascades to $J/\Psi$ and
$J/\Psi\rightarrow l+l$. We further wish to mention that in 2015, the BESIII
collaboration[6] independently confirmed the existence of the charmonium-like
state $X(3823)$ in the $\chi_{c1}\gamma$ system with a statistical
significance of $6.2\sigma$. This observation was made in the $e^{-}e^{+}$
collision process, specifically
$e^{-}e^{+}\rightarrow\pi^{+}\pi^{-}\chi_{c1}\gamma$ [6]. The measured mass of
$X(3823)$ in this study is $3823.7\pm 1.8\pm 0.7$ MeV[6], which is consistent
with the mass measurement previously reported by the Belle Collaboration, thus
providing further validation for the existence of the $X(3823)$ particle. The
characteristics of $X(3823)$ align with $1^{3}D_{2}$ state, which is now
officially recognised as $\Psi_{2}(3823)$[1]. The primary decay modes of
$\Psi_{2}(3823)$ are expected to involve radiative and hadronic transitions
into other charmonium states. According to phenomenological studies, it is
anticipated that the partial widths for
$\Gamma_{\psi_{2}\rightarrow\chi_{c1}\gamma}$ are in the range of
approximately 200 to 300 keV, while
$\Gamma_{\psi_{2}\rightarrow\chi_{c2}\gamma}$ is estimated to be around 60
keV. Additionally, $\Gamma_{\Psi_{2}\rightarrow\ J/{\psi}\pi\pi}$ is predicted
to have a partial width of approximately 160 keV.

However, it’s worth noting that experimental results from BESIII have provided
upper limits, indicating that the ratios of $\Gamma_{\Psi_{2}(3823)\rightarrow
J/{\Psi}(\pi\pi)}$ to $\Gamma_{\Psi_{2}(3823)\rightarrow\chi_{c1}(\gamma)}$
are less than 0.06 and $\Gamma_{\Psi_{2}(3823)\rightarrow
J/{\Psi}(\pi^{0}\pi^{0})}$ to
$\Gamma_{\Psi_{2}(3823)\rightarrow\chi_{c1}\gamma}$ are less than 0.11. These
experimental limits are in stark contrast to the theoretical predictions. This
suggests that the decay channel $\psi_{2}\rightarrow\chi_{c1}\gamma$ may be
the dominant mode for $\psi_{2}$ decay. In experimental observations, LHCb has
established an upper limit of $\Gamma_{\psi_{2}}<5.2MeV$, while a more recent
measurement by BESIII has refined this constraint to $\Gamma_{\psi_{2}}<2.9$
MeV with $90\%$ confidence.

Also, recently, [7] searched for processes,
$e^{-}e^{+}\rightarrow\phi\chi_{c0,c1,c2}$ using a data sample collected with
the BESIII detector operating at the BEPCII storage ring at a centre-of-mass
energy, $\sqrt{s}=4.6$ GeV. They further observed the processes,
$e^{-}e^{+}\rightarrow\phi\chi_{c1}$, and
$e^{-}e^{+}\rightarrow\phi\chi_{c2}$, with a statistical significance of more
than $10\sigma$.

Now, we have earlier studied the mass spectrum, and the leptonic decays of
ground and excited states of heavy-light axial vector mesons such as $1^{++}$
and $1^{+-}$ in [8], along with the radiative E1 decay processes[9] involving
$1^{+-}$ mesons, such as, $1^{+-}\rightarrow 0^{-+}\gamma$, and
$0^{-+}\rightarrow 1^{+-}\gamma$. Thus in the present work, we were interested
in calculating the radiative decay processes involving $1^{++}$ mesons such
as, $1^{++}\rightarrow 1^{--}\gamma$, and $1^{--}\rightarrow 1^{++}\gamma$.
Here we wish to mention that radiative decays ${}^{3}P_{1}$ $\leftrightarrow$
${}^{3}S_{1}$ receive contributions from both the $E1$ and $M2$ transitions,
where $M2$ transitions have obvious corrections to the radiative transitions
of $\chi_{c1}\rightarrow J/\Psi\gamma$, by interfering with the E1
transitions. The corrections from the $M2$ transition might be $(10-20)\%$
[10, 11] of the total partial decay widths.

There have been a number of theoretical studies of these E1 transitions based
on a range of different models[13, 14, 15, 16, 17]. Radiative transitions
between $1^{--}$ and $1^{++}$ which proceed through the emission of a photon,
are characterized by $|\Delta L|=1$, Thus in these transitions there is a
change in parity between the initial and final hadronic states, though the
total parity remains conserved. We determine the radiative decay widths of
heavy-light quarkonia for the processes mentioned above using a $4\times 4$
Bethe-Salpeter equation (BSE). This BSE approach is fully relativistic, taking
into account the relativistic effects of quark spins and providing a
consistent treatment of the internal motion of constituent quarks within the
hadron due to its covariant structure. We have wave functions that fulfill the
three-dimensional (3D) Bethe-Salpeter equation (BSE). This equation is derived
from the four-dimensional (4D) BSE by employing the covariant instantaneous
ansatz (CIA),which is a Lorentz-invariant extension of the instantaneous
approximation. Consequently, our wave functions incorporate relativistic
effects.

This paper is organized as follows: In Sec. II, we give the general
formulation of the process ${H\rightarrow H^{\prime}+\gamma}$ in the framework
of a $4\times 4$ Bethe-Salpeter equation under the covariant instantaneous
ansatz. In Sec. III, we calculate the single photon decay widths for the
process $A^{+}\rightarrow V\gamma$ . In Sec. IV, we deal with the process
$V\rightarrow A^{+}\gamma$. In Sec. V, we give numerical results and
discussions.

##  2\. Radiative decays of quarkonia through $V\rightarrow A^{+}\gamma$

Radiative M1 and E1 decays proceed through the process, $H\rightarrow
H^{\prime}\gamma$, where $H$ and $H^{\prime}$ are initial and final quarkonia,
that can be considered as two Lorentz frames. At leading order, two triangle
quark-loop diagram shown in Fig.1 contribute to this process. The second
diagram is obtained from the first by reversing the direction of internal
fermion lines. We calculate this process in the $4\times 4$ representation of
BSE. We had earlier calculated some of the M1 and E1 decays of quarkonium in
[18, 9]. Thus we will retain the same notation as in [18, 9].

Figure 1: Radiative decays of equal mass quarkonia

Now, the EM transition amplitude of the process, $V\rightarrow A^{+}\gamma$
(where $V=1^{--}$, and $A^{++}=1^{++}$) is

$M_{fi}=-2e\int\frac{d^{4}q}{(2\pi)^{4}}Tr[\overline{\Psi}_{A}(P^{\prime},q^{\prime}){\not}\epsilon^{\lambda^{\prime}}\Psi_{V}(P,q)S_{F}^{-1}(-p_{2})].$
(1)

where the total $M_{fi}$ is two times the contribution from Diagram 1, since
both diagrams contribute equally. Here,
$\overline{\Psi}_{A}(P^{\prime},q^{\prime})=\gamma_{4}\Psi_{A}^{{\dagger}}(P^{\prime},q^{\prime})\gamma_{4}$
is adjoint 4D BS wave function for final axial meson, while $\Psi_{V}(P,q)$ is
BS wave function of initial vector meson. Here we have taken, $P,q$ and
$\epsilon^{\lambda}$ as the total momentum, internal momentum, and
polarization vector of initial hadron, $H$, while $P^{\prime},q^{\prime}$, and
$\epsilon^{\lambda^{\prime\prime}}$ as the corresponding variables of the
final hadron, $H^{\prime}$. Further, $k$, and $\epsilon^{\lambda^{\prime}}$
are momentum and polarization vectors of emitted photon. Thus if $p_{1,2}$,
and $p^{\prime}_{1,2}$ are the momenta of the two quarks in initial and final
hadron respectively, then, the momentum relations for initial hadron are:
$P=p_{1}+p_{2}$, and $p_{1,2}=\frac{1}{2}P\pm q$, while for final hadron,
$P^{\prime}=p^{\prime}_{1}+p^{\prime}_{2}$, and
$p^{\prime}_{1,2}=\frac{1}{2}P^{\prime}\pm q^{\prime}$, while the relation
between their internal momentum is: $q^{\prime}=q-\frac{1}{2}k$ (see [9] for
details).

To simplify calculations, we prefer to work in the rest frame of the initial
hadron, where we decompose the internal momentum
$q^{\prime}=(\hat{q}^{\prime},iM\sigma^{\prime})$ of the final hadron into two
components with $\hat{q}^{\prime}=q^{\prime}-\sigma^{\prime}P$, and
$\sigma^{\prime}=\frac{q^{\prime}.P}{P^{2}}$, with $\hat{q}^{\prime}$ being
transverse to the initial hadron momentum $P$, and $\sigma^{\prime}$ being
longitudinal to $P$. Thus, $P.\hat{q}^{\prime}=0$. The relationship between
the transverse components, and longitudinal components of internal momenta of
the two hadrons is given as in Eq.(13) and (15) of [9]. From the Feynman
diagrams we see that conservation of momentum demands that, $P=P^{\prime}+k$.
Now, for the first diagram, we have the kinematical relations,
$p_{1}=p^{\prime}_{1}+k$, and $-p_{2}=-p^{\prime}_{2}$, where $k=P-P^{\prime}$
is the momentum of the emitted photon. From kinematical relations, it can be
shown that the energy of final meson can be expressed as [18, 9],
$E^{\prime}=\frac{M^{2}+M^{\prime 2}}{2M}$. It is to be noted that 4D BS wave
functions of initial vector meson involved in the process (in variable,
$\hat{q}$) is $\Psi_{V}(P,q)=S_{F}(p_{1})\Gamma_{V}(\hat{q})S_{F}(-p_{2})$,
while the wave function of the final axial vector meson is,
$\Psi_{A}(P^{\prime},q^{\prime})=S_{F}(p^{\prime}_{1})\Gamma_{A}(\hat{q}^{\prime})S_{F}(-p^{\prime}_{2})$,
where, $\hat{q}^{\prime}=q^{\prime}-\frac{q^{\prime}.P}{P^{2}}P$ is transverse
to initial hadron momentum, $P$, and $\Gamma_{A}(\hat{q}^{\prime})$ is the
hadron-quark vertex function for the final meson. Now, we reduce the above
equation to the effective 3D form by integrating over the longitudinal
component, $Md\sigma$. This can be expressed as,

$M_{fi}=-2e\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\int\frac{iMd\sigma}{(2\pi)}Tr[\bar{\Gamma}_{A}(\hat{q}^{\prime})S_{F}(p^{\prime}_{1}){\not}\epsilon^{\prime\prime}S_{F}(p_{1})\Gamma_{V}(\hat{q})S_{F}(-p_{2})].$
(2)

To calculate further, we split the quark and anti-quark propagators in terms
of projection operators, and carrying out in integrations in complex $\sigma$
plane over the poles of the quark propagators, we obtain the expression for
effective 3D form of transition amplitude, $M^{1}_{fi}$ under Covariant
Instantaneous Ansatz for Diagram 1, as[18, 9],

$M^{1}_{fi}=-ie\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\frac{1}{M^{2}}Tr\bigg{[}\alpha_{1}{\not}P\overline{\psi}_{A}^{++}(\hat{q}^{\prime}){\not}\epsilon^{\prime\prime}\psi_{V}^{++}(\hat{q})+\alpha_{2}{\not}P\overline{\psi}_{A}^{++}(\hat{q}^{\prime}){\not}\epsilon^{\prime\prime}\psi_{V}^{--}(\hat{q})\\\
+\alpha_{3}{\not}P\overline{\psi}_{A}^{--}(\hat{q}^{\prime}){\not}\epsilon^{\prime\prime}\psi_{V}^{++}(\hat{q})+\alpha_{4}{\not}P\overline{\psi}_{A}^{--}(\hat{q}^{\prime}){\not}\epsilon^{\prime\prime}\psi_{V}^{--}(\hat{q})\bigg{]}$
(3)

Here, $\alpha_{1},...,\alpha_{4}$ [18, 9] are the results of contour integrals
over $Md\sigma$. $M_{fi}$ above is a generalized expression for transition
amplitude expressed as a linear superposition of terms involving all possible
combinations of $++$, and $--$ components of Salpeter wave functions of final
and initial hadrons. Now on lines on [12], to simplify calculations, we ignore
the terms involving $\psi^{+-},\psi^{-+}$ and $\psi^{--}$ that contribute less
than $1\%$, and thus consider only the terms involving $\psi^{++}$ that give
dominant contribution, and write $M^{1}_{fi}$ as,

$M^{1}_{fi}=-ie\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\frac{1}{M^{2}}Tr\bigg{[}\alpha_{1}{\not}P\overline{\psi}_{A}^{++}(\hat{q}^{\prime}){\not}\epsilon^{\prime\prime}\psi_{V}^{++}(\hat{q})\bigg{]}.$
(4)

where the $++$ components of the B.S. wave function for axial vector meson,
and vector meson [24, 26] respectively, can be obtained from the equations,

$\psi_{A}^{\pm\pm}(\hat{q}^{\prime})=\Lambda^{\pm}_{1}(\hat{q}^{\prime})\frac{{\not}P}{M}\psi_{A}(\hat{q}^{\prime})\frac{{\not}P}{M}\Lambda^{\pm}_{2}(\hat{q}^{\prime}),$
(5)

and

$\psi_{V}^{\pm\pm}(\hat{q})=\Lambda^{\pm}_{1}(\hat{q})\frac{{\not}P}{M}\psi_{V}(\hat{q})\frac{{\not}P}{M}\Lambda^{\pm}_{2}(\hat{q}).$
(6)

Now, to calculate the process, we need the 3D BS wave functions,
$\psi_{V}(\hat{q})$, and $\psi_{A}(\hat{q})$ for vector and axial vector
mesons respectively. We then start with the general 4D form of BS wave
functions [19, 20], expressed as a superposition of various Dirac structures,
each multiplied by scalar amplitudes, $f_{i}(P,q)$, which are all independent.
We then make use of their 3D decomposition under Covariant Instantaneous
Ansatz [22, 21, 24, 25, 22, 23]. We further make use of the 3D Salpeter
equations that are obtained by 3D reduction of the 4D Bethe-Salpeter equation
under Covariant Instantaneous Ansatz following a sequence of steps outlined in
[25, 22] that are expressed as,

$\displaystyle(M-\omega_{1}-\omega_{2})\psi^{++}(\hat{q})=\Lambda_{1}^{+}(\hat{q})\Gamma(\hat{q})\Lambda_{2}^{+}(\hat{q}),$
$\displaystyle(M+\omega_{1}+\omega_{2})\psi^{--}(\hat{q})=-\Lambda_{1}^{-}(\hat{q})\Gamma(\hat{q})\Lambda_{2}^{-}(\hat{q}),$
$\displaystyle\psi^{+-}(\hat{q})=0,$ $\displaystyle\psi^{-+}(\hat{q})=0.$ (7)

These 3D Salpeter equations in Eq.(7) depend on the Lorentz-invariant variable
$\hat{q}^{2}=q^{2}-(q.P)^{2}/P^{2}$ which is a scalar, whose validity extends
over the entire time-like region of the 4D space, while also keeping contact
with the surface, $P.q=0$ (hadron rest frame). Putting these 3D wave functions
into the last two Salpeter equations, $\Psi^{+-}(\hat{q})=0$, and
$\Psi^{-+}(\hat{q})=0$ provide constraint relations between the amplitudes of
Dirac structures for any meson. Then, putting this wave function into the
first two Salpeter equations leads to two coupled integral equations in two
independent Salpeter amplitudes. With the structure of the interaction kernel
incorporated, these coupled equations are then decoupled and reduced to two
identical algebraic equations, which are in turn reduced to mass spectral
equations. The 3D wave functions for $1^{--}$ and $1^{++}$ (for equal mass
quarkonia) are expressed as [22, 25]:

$\displaystyle\psi_{V}(\hat{q})=N_{V}\bigg{[}iM{\not}\epsilon+{\not}\epsilon{\not}P+\frac{i}{2m}{\not}P{\not}\epsilon{\not}\hat{q}+\frac{i}{2m}\hat{q}.\epsilon{\not}P\bigg{]}\phi_{V}(\hat{q});$
$\displaystyle\psi_{A}(\hat{q}^{\prime})=N_{A}\gamma_{5}\bigg{[}iM^{\prime}{\not}\epsilon^{\prime}+{\not}\epsilon^{\prime}{\not}P^{\prime}+2i\frac{{\not}\epsilon^{\prime}{\not}P^{\prime}{\not}\hat{q}^{\prime}}{M}\bigg{]}\phi_{A^{+}}(\hat{q}^{\prime}),$
(8)

where as mentioned above, $\epsilon^{\lambda}$ and
$\epsilon^{\lambda^{\prime\prime}}$ are the polarization vectors of vector and
axial vector mesons respectively. In [24, 21, 25, 9], we had derived the
analytic forms of 3D radial wave functions for $S,P,V$ and $A$ mesons from
their respective mass spectral equations in an approximate harmonic oscillator
basis, that were used to calculate various M1 and E1 transitions in [18, 9].
These wave functions for $V$ and $A$ quarkonia are:

$\begin{split}\phi_{V}(1S,\hat{q})&=\frac{1}{\pi^{3/4}}\frac{1}{\beta_{V}^{3/2}}e^{-\frac{\hat{q}^{2}}{2\beta_{V}^{2}}}\\\
\phi_{V}(2S,\hat{q})&=\sqrt{\frac{3}{2}}\frac{1}{\pi^{3/4}}\frac{1}{\beta_{V}^{3/2}}\bigg{(}1-\frac{2\hat{q}^{2}}{3\beta_{V}^{2}}\bigg{)}e^{-\frac{\hat{q}^{2}}{2\beta_{V}^{2}}}\\\
\phi_{V}(1D,\hat{q})&=\sqrt{\frac{4}{15}}\frac{1}{\pi^{3/4}}\frac{1}{\beta_{V}^{7/2}}\hat{q}^{2}e^{-\frac{\hat{q}^{2}}{2\beta_{V}^{2}}}\\\
\phi_{V}(3S,\hat{q})&=\sqrt{\frac{15}{8}}\frac{1}{\pi^{3/4}}\frac{1}{\beta_{V}^{3/2}}\bigg{(}1-\frac{4\hat{q}^{2}}{3\beta_{V}^{2}}+\frac{4\hat{q}^{4}}{15\beta_{V}^{4}}\bigg{)}e^{-\frac{\hat{q}^{2}}{2\beta_{V}^{2}}}\\\
\phi_{A}(1p,\hat{q})&=\sqrt{\frac{2}{3}}\frac{1}{\pi^{3/4}}\frac{1}{\beta_{A}^{5/2}}\hat{q}e^{-\frac{\hat{q}^{2}}{2\beta_{A}^{2}}}\\\
\phi_{A}(2p,\hat{q})&=\sqrt{\frac{5}{3}}\frac{1}{\pi^{5/2}}\hat{q}(1-\frac{2\hat{q}^{2}}{5\beta_{A}^{2}})e^{-\frac{\hat{q}^{2}}{2\beta_{A}^{2}}}\\\
\end{split}$ (9)

where $\beta_{A,V}$ are their inverse range parameters, with
$\beta_{A}=(\frac{2}{3}M\omega_{q\bar{q}}^{2})^{\frac{1}{4}}$, and
$\beta_{V}=(m\omega_{q\bar{q}}^{2})^{\frac{1}{4}}$[25, 8].

Now, the $++$ components of vector ($1^{--}$), and axial ($1^{++}$) vector
meson wave functions are:

$\displaystyle\Psi_{V}^{++}(\hat{q})=N_{V}\phi_{V}(\hat{q})[i{\not}\epsilon
a_{1}+{\not}\epsilon{\not}Pa_{2}+i{\not}\epsilon{\not}P{\not}\hat{q}a_{3}+{\not}\epsilon{\not}\hat{q}a_{4}+i{\not}Pa_{5}+{\not}P{\not}\hat{q}a_{6}+i{\not}\hat{q}{\not}P{\not}\epsilon
a_{7}+i{\not}\hat{q}a_{8}]$
$\displaystyle\overline{\Psi}_{A}^{++}(\hat{q}^{\prime})=N_{A}\phi_{V}(\hat{q}^{\prime})\gamma_{5}[i{\not}\epsilon^{\prime}b_{1}+{\not}\epsilon^{\prime}{\not}P^{\prime}b_{2}+i{\not}P^{\prime}b_{5}+i{\not}P^{\prime}{\not}\hat{q}^{\prime}{\not}\epsilon^{\prime}b_{3}].$
(10)

where the coefficients associated with various Dirac structures for $V$ meson
are:

$\displaystyle
a_{1}=\frac{M}{4}-\frac{2mM}{\omega}-\frac{M\hat{q}^{2}}{8\omega
m}+\frac{m^{2}M}{4\omega^{2}}+\frac{M\hat{q}^{2}}{4\omega^{2}}-\frac{\hat{q}^{2}}{4\omega}$
$\displaystyle
a_{2}=-\frac{3m}{8\omega}+\frac{1}{4}+\frac{m^{2}}{4\omega^{2}}+\frac{\hat{q}^{2}}{2\omega^{2}}$
$\displaystyle
a_{3}=\frac{1}{4\omega}-\frac{\hat{q}^{2}}{8\omega^{2}m}+\frac{1}{8m}+\frac{m}{8\omega^{2}}$
$\displaystyle a_{4}=0$ $\displaystyle
a_{5}=(\epsilon.P)[\frac{\hat{q}^{2}}{4m\omega
M}+\frac{m^{2}}{4\omega^{2}M}-\frac{m}{2\omega
M}+\frac{\hat{q}^{2}}{2\omega^{2}M}]-(\hat{q}.\epsilon)[\frac{1}{8m}-\frac{m}{8\omega^{2}}-\frac{3\hat{q}^{2}}{8\omega^{2}m}]$
$\displaystyle a_{6}=(\epsilon.P)[\frac{m}{4\omega^{2}M}+\frac{1}{2\omega
M}+(\hat{q}.\epsilon)[-\frac{9}{8\omega^{2}}]$ $\displaystyle
a_{7}=-\frac{1}{4\omega}+\frac{m}{4\omega^{2}}$ $\displaystyle
a_{8}=(\hat{q}.\epsilon)[-\frac{M}{2\omega^{2}}+\frac{1}{4\omega}]$ (11)

Similarly we write down the coefficients of $b_{1},...,b_{7}$ associated with
$A$ meson are:

$\displaystyle
b_{1}=\frac{M^{\prime}}{4}+\frac{mM^{\prime}}{2\omega^{\prime}}+\frac{m^{2}M^{\prime}}{4\omega^{\prime
2}}-\frac{M^{\prime}\hat{q}^{\prime 2}}{4\omega^{\prime
2}}-\frac{\hat{q}^{\prime 2}}{m}$ $\displaystyle
b_{2}=\frac{1}{4}-\frac{m^{2}}{4\omega^{\prime 2}}+\frac{\hat{q}^{\prime
2}}{4\omega^{\prime 2}};~{}b_{3}=\frac{1}{2M^{\prime}}+\frac{\hat{q}^{\prime
2}}{2M^{\prime}\omega^{\prime 2}}$ $\displaystyle
b_{4}=\frac{3m}{2\omega^{\prime}};~{}b_{5}=-\frac{\hat{q^{\prime}}.\epsilon^{\prime}}{\omega^{\prime}}$
$\displaystyle
b_{6}=-(\hat{q}^{\prime}.\epsilon^{\prime})\frac{\hat{q}^{\prime
2}}{M^{\prime}\omega^{\prime
2}};~{}b_{7}=(\hat{q}^{\prime}.\epsilon^{\prime})\frac{m}{M^{\prime}\omega^{\prime
2}}.$ (12)

Now, our mass spectrum, and the 3D wave functions $\phi(\hat{q})$ in Eq.(9)
were calculated from the mass spectral equations, that are in turn derived
from the 3D Salpeter equations in Eq.(7) in the rest frame of the hadron
(please see [24, 25].

The Bethe-Salpeter normalizers, $N_{V}$, and $N_{A}$ are obtained through
current conservation condition,

$2iP_{\mu}=\int\frac{d^{4}q}{(2\pi)^{4}}Tr\bigg{[}\overline{\psi}(P,q)[\frac{\partial}{\partial
P_{\mu}}S_{F}^{-1}(p_{1})]\psi(P,q)S_{F}^{-1}(-p_{2})\bigg{]}+(1\leftrightarrow
2),$ (13)

This leads their algebraic forms,

$\displaystyle
N_{A}^{-2}=\frac{4}{M}\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\phi^{2}_{A}(\hat{q})[2mM-M^{2}+2\hat{q}^{2}];$
$\displaystyle
N_{V}^{-2}=\frac{4}{M}\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\phi^{2}_{V}(\hat{q})[-\frac{5M^{2}\hat{q}^{2}}{6m^{2}}+\frac{20M\hat{q}^{2}}{3m}].$
(14)

Here, it is to be mentioned that, the transverse component of internal
momentum of the final meson can be expressed as,
$\hat{q}^{\prime}=\hat{q}+\hat{m}_{2}\hat{P}^{\prime}$[18], where
$\hat{m}_{2}$ act as momentum partitioning parameters. Now squaring both sides
of this relation, and making use of the fact that $\hat{P}^{\prime}$, and
$\hat{q}$ are both transverse to the initial hadron momentum, along with the
momentum of emitted meson,
$|\hat{P}^{\prime}|=|\overrightarrow{P}|=\frac{M^{2}-M^{\prime 2}}{2M}$, we
can express the relationship between $\hat{q}^{\prime 2}$, and $\hat{q}^{2}$
as [18]:

$\hat{q}^{\prime 2}=\hat{q}^{2}+2\hat{m}_{2}\frac{(M^{2}-M^{\prime
2})}{2M}|\hat{q}|+\hat{m}_{2}^{2}\frac{(M^{2}-M^{\prime 2})^{2}}{4M^{2}}$ (15)

where,$|\hat{q}|$ is the length of the 3-D vector, $\hat{q}$, defined as
$|\hat{q}|=\sqrt{\hat{q}^{2}}=\sqrt{q^{2}-(q.P)^{2}/P^{2}}$, and is a Lorentz-
invariant [26, 22]) variable.

The transition amplitude, $M_{fi}$ is expressed as,

$M_{fi}=-2ie_{Q}\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\frac{1}{M^{2}}Tr[\alpha_{1}{\not}P\bar{\Psi}^{++}_{A}(\hat{q}^{\prime}){\not}\epsilon^{\prime\prime}\Psi^{++}_{V}(\hat{q})]$
(16)

The trace over the gamma matrices can be expressed as,

$[Tr]=\alpha_{1}\bigg{[}X^{\prime}_{1}\epsilon_{\mu\nu\alpha\beta}\epsilon^{\prime\prime}_{\mu}\epsilon^{\prime}_{\nu}\epsilon_{\alpha}P_{\beta}+X^{\prime}_{2}\epsilon_{\mu\nu\alpha\beta}k_{\mu}\epsilon_{\nu}\epsilon^{\prime}_{\alpha}\epsilon^{\prime\prime}_{\beta}+X^{\prime}_{3}\epsilon_{\mu\nu\alpha\beta}P_{\mu}k_{\nu}\epsilon^{\prime}_{\alpha}\epsilon_{\beta}+X^{\prime}_{4}\epsilon_{\mu\nu\alpha\beta}k_{\mu}P_{\nu}\epsilon_{\alpha}\epsilon^{\prime\prime}_{\beta}+X^{\prime}_{5}\epsilon_{\mu\nu\alpha\beta}P_{\mu}k_{\nu}\epsilon^{\prime}_{\alpha}\epsilon^{\prime\prime}_{\beta}\bigg{]},$
(17)

where,

$\displaystyle
X^{\prime}_{1}=b_{1}a_{1}-M^{2}b_{2}a_{2}-(\frac{M^{2}}{2}-\frac{1}{2}P^{\prime}.P)b_{3}a_{1}-M^{2}\hat{q}^{2}b_{3}a_{3}+M^{2}\hat{q}^{2}b_{3}a_{7}-(\frac{M^{2}}{2}-\frac{1}{2}P^{\prime}.P)b_{4}a_{2}-\hat{q}^{2}b_{4}a_{4}$
$\displaystyle
X^{\prime}_{2}=M^{2}b_{2}a_{2}-(M^{2}-\frac{1}{2}P^{\prime}.P)b_{3}a_{1}+M^{2}\hat{q}^{2}b_{3}a_{3}-M^{2}\hat{q}^{2}b_{3}a_{7}+\frac{M^{2}}{2}b_{4}a_{2}$
$\displaystyle X^{\prime}_{3}=[-1-\frac{P^{\prime}.P}{M^{2}}]b_{3}a_{1}$
$\displaystyle
X^{\prime}_{4}=b_{6}a_{1}+(M^{2}-\frac{P^{\prime}.P}{M^{2}})b_{7}a_{2}-\hat{q}^{2}b_{7}a_{4}$
$\displaystyle
X^{\prime}_{5}=(-M^{2}-\frac{1}{2}P^{\prime}.P)b_{3}a_{5}-\hat{q}^{2}b_{3}a_{8}.$
(18)

Thus, we can write $M_{fi}$ as,

$\displaystyle
M_{fi}=-2ie_{Q}\frac{N_{V}N_{A}}{M^{2}}[G_{1}~{}\epsilon_{\mu\nu\alpha\beta}\epsilon^{\prime\prime}_{\mu}\epsilon^{\prime}_{\nu}\epsilon_{\alpha}P_{\beta}+G_{2}~{}\epsilon_{\mu\nu\alpha\beta}k_{\mu}\epsilon_{\nu}\epsilon^{\prime}_{\alpha}\epsilon^{\prime\prime}_{\beta}+(P.\epsilon^{\prime\prime})G_{3}~{}\epsilon_{\mu\nu\alpha\beta}P_{\mu}k_{\nu}\epsilon^{\prime}_{\alpha}\epsilon_{\beta}+$
$\displaystyle(I.\epsilon^{\prime})G_{4}~{}\epsilon_{\mu\nu\alpha\beta}k_{\mu}P_{\nu}\epsilon_{\alpha}\epsilon^{\prime\prime}_{\beta}+(I.\epsilon)G_{5}~{}\epsilon_{\mu\nu\alpha\beta}P_{\mu}k_{\nu}\epsilon^{\prime}_{\alpha}\epsilon^{\prime\prime}_{\beta}],$
(19)

where,

$\displaystyle
G_{1}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\alpha_{1}X^{\prime}_{1}(\hat{q})\phi_{A}(\hat{q})\phi_{V}(\hat{q})$
$\displaystyle
G_{2}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\alpha_{1}X^{\prime}_{2}(\hat{q})\phi_{A}(\hat{q})\phi_{V}(\hat{q})$
$\displaystyle G_{3}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\alpha
1X^{\prime}_{3}\phi_{A}(\hat{q})\phi_{V}(\hat{q})$ $\displaystyle
G_{4}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\alpha_{1}X^{\prime}_{4}\phi_{A}(\hat{q})\phi_{V}(\hat{q})$
$\displaystyle
G_{5}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\alpha_{1}X^{\prime}_{5}\phi_{A}(\hat{q})\phi_{V}(\hat{q}).$
(20)

Now, since $\hat{q}$ is an effective 3D vector, we can write,
$\hat{q}.\epsilon^{\prime}=|\hat{q}|(I.\epsilon^{\prime})$, where we take
$I_{\mu}$ as a unit vector along the direction of $\hat{q}$, and is expressed
as, $I=\frac{\hat{q}}{|\hat{q}|}$, where $|\hat{q}|$ is a Lorentz-invariant
variable[22] defined above. Taking the initial vector meson to be in its rest
frame, where $P=(\overrightarrow{0},iM$), and since the polarization vector of
final axial meson
$\epsilon^{\prime\prime}=(\overrightarrow{\epsilon}^{\prime\prime},i0)$, due
to which the dot product, $P.\epsilon^{\prime\prime}=0$. Thus the term with
$G_{3}$ in $M_{fi}$ will not contribute. Further, $G_{4}$ and $G_{5}$
correspond to magnetic dipole $M2$ transitions, and make sub-dominant
contributions to $M_{fi}$. These sub- dominant contributions are expected to
contribute $10-20\%$ [10, 11] of the decay widths, and are thus ignored in
present calculation. The dominant contributions to decay widths, that
correspond to $E1$ transitions, should come from $G_{1}$, and $G_{2}$ terms,
in terms of which we write $M_{fi}$ as,

$M_{fi}=-2ie_{Q}\frac{N_{V}N_{A}}{M^{2}}[G_{1}~{}\epsilon_{\mu\nu\alpha\beta}~{}\epsilon^{\lambda^{\prime\prime}}_{\mu}\epsilon^{\lambda^{\prime}}_{\nu}\epsilon^{\lambda}_{\alpha}P_{\beta}+G_{2}~{}\epsilon_{\mu\nu\alpha\beta}~{}\epsilon^{\lambda^{\prime\prime}}_{\mu}\epsilon^{\lambda^{\prime}}_{\nu}\epsilon^{\lambda}_{\alpha}k_{\beta}].$
(21)

Thus, we can write $M_{fi}$ as,

$\displaystyle
M_{fi}=-2ie_{Q}\frac{N_{V}N_{A}}{M^{2}}\epsilon^{\lambda^{\prime\prime}}_{\mu}\epsilon^{\lambda^{\prime}}_{\nu}\epsilon^{\lambda}_{\alpha}M_{\mu\nu\alpha},$
$\displaystyle
M_{\mu\nu\alpha}=G_{1}~{}\epsilon_{\mu\nu\alpha\beta}P_{\beta}+G_{2}~{}\epsilon_{\mu\nu\alpha\beta}k_{\beta}.$
(22)

Now condition for gauge invariance, $k_{\mu}M_{\mu\nu\alpha}=0$, leads to,
$G_{1}=0$. Thus, we obtain,

$M_{fi}=-2ie_{Q}\frac{N_{V}N_{A}}{M^{2}}G_{2}~{}\epsilon_{\mu\nu\alpha\beta}~{}\epsilon^{\lambda^{\prime\prime}}_{\mu}\epsilon^{\lambda^{\prime}}_{\nu}\epsilon^{\lambda}_{\alpha}k_{\beta}.$
(23)

Now, to calculate the decay widths, we need to calculate the spin averaged
amplitude modulus square, $|\overline{M}_{fi}|^{2}$, where
$|\overline{M}_{fi}|^{2}=\frac{1}{2j+1}\sum_{\lambda,\lambda^{\prime}}|{M}_{fi}|^{2}$,
where we average over the initial polarization states $\lambda$ of V-meson,
and sum over the final polarization states, $\lambda^{\prime\prime}$ of
A-meson, and $\lambda^{\prime}$ of photon. We make use of the normalizations,
$\Sigma_{\lambda}\epsilon_{\mu}^{\lambda}\epsilon_{\nu}^{\lambda}=\frac{1}{3}(\delta_{\mu\nu}+\frac{P_{\mu}P_{\nu}}{M^{2}})$
for vector meson,
$\Sigma_{\lambda^{\prime\prime}}\epsilon_{\mu}^{\lambda^{\prime\prime}}\epsilon_{\nu}^{\lambda^{\prime\prime}}=\frac{1}{3}(\delta_{\mu\nu}+\frac{P^{\prime\prime}_{\mu}P^{\prime\prime}_{\nu}}{M^{\prime\prime
2}})$ for final axial meson, and
$\Sigma_{\lambda^{\prime}}\epsilon_{\mu}^{\lambda^{\prime}}\epsilon_{\nu}^{\lambda^{\prime}}=\delta_{\mu\nu}$,
for the emitted photon, with $M_{fi}$ taken from the previous equation.

We obtain,

$|\overline{M}_{fi}|^{2}=\frac{4}{9}\frac{e_{Q}^{2}N_{A}^{2}N_{V}^{2}}{M^{4}}|G_{2}|^{2}\bigg{[}\frac{(P.k)^{2}}{M^{2}}+\frac{(P^{\prime}.k)^{2}}{M^{\prime
2}}-\frac{(P^{\prime}.k)(P^{\prime}.P)(P.k)}{M^{2}M^{\prime 2}}\bigg{]}.$ (24)

Taking charge of charmed quark, $e_{Q}=\frac{2}{3}e$, and the dot products of
momenta calculated in the rest frame of initial vector meson:
$P^{\prime}.P=-\frac{(M^{2}+M^{\prime 2})}{2}$, $P.k=-\frac{M^{2}-M^{\prime
2}}{2}$, $P^{\prime}.k=-\frac{M^{2}-M^{\prime 2}}{2}$ we can write,

$|\overline{M}_{fi}|^{2}=\frac{8\pi\alpha_{em}}{81}N_{A}^{2}N_{V}^{2}|G_{2}|^{2}\frac{(M^{2}-M^{\prime
2})^{2}(M^{2}+M^{\prime 2})}{M^{6}M^{\prime 2}}.$ (25)

The decay width of the process ($V\rightarrow A\gamma$) in the rest frame of
the initial vector meson is expressed as

$\Gamma_{V\rightarrow A\gamma}=\frac{|\overline{M}_{fi}|^{2}}{8\pi
M^{2}}|\overrightarrow{P^{\prime}}|,$ (26)

where we make use of the fact that modulus of the momentum of the emitted
axial vector meson can be expressed in terms of masses of particles as,
$|\overrightarrow{P^{\prime}}|=|\overrightarrow{k}|=\omega_{k}=\frac{1}{2M}(M^{2}-M^{\prime
2})$, where, $\omega_{k}$ is the kinematically allowed energy of the emitted
photon.

The vector and axial vector meson masses calculated in our framework[25, 8]
are given in Table 1.

|  BSE-CIA[25, 8] |  Expt.[2]
---|---|---
$M_{J/\Psi}(1S)$ |  3.0970 |  3.0969$\pm$0.000011
$M_{\Psi}(2S)$ |  3.6676 |  3.6861$\pm$0.00034
$M_{\Psi}(1D)$ |  3.7716 |  3.7773$\pm$0.00033
$M_{\chi_{c1}}(1P)$ |  3.4783 |  3.510$\pm$ 0.0007
$M_{\chi_{c1}}(2P)$ |  3.9560 |  3.871.69$\pm$ 0.00017
Table 1: Mass spectra of ground and excited states of vector ($1^{--}$) and
axial vector ($1^{++}$) quarkonia (in GeV) in BSE-CIA used in the transitions
studied in this work, along with data

The E1 decay widths for $V\rightarrow A\gamma$ involving vector ($1^{++}$),
and axial vector ($1^{++}$) $c\bar{c}$ mesons are calculated in our model
along with results of other models are given Table 2.

|  BSE-CIA |  Expt.[1] |  CQM[11] |  RQM[13]
---|---|---|---|---
$\Gamma_{\psi(2S)\rightarrow\chi_{c1}(1P)\gamma}$ |  23.54 |  28.665$\pm$0.78 |  20 |  22.9
$\Gamma_{\psi(3S)\rightarrow\chi_{c1}(1P)\gamma}$ |  22.72 |  $<$ 272$\pm$ 34 |  |
$\Gamma_{\psi(1D)\rightarrow\chi_{c1}(1P)\gamma}$ |  48.89 |  67.77$\pm$2.49 |  70 |
Table 2: Radiative decay widths of heavy-light $1^{--}$ charmonia (in keV) for
E1 transitions in BSE, along with experimental data and results of other
models.

##  3\. Radiative transition, $A\rightarrow V\gamma$

We now consider the radiative decay process, $A\rightarrow V\gamma$. At
leading order, this process again involves a quark-triangle loop, but with two
hadron-quark vertices, the initial hadron, $A=1^{++}$, and the final hadron,
$V=1^{--}$. There are two diagrams that contribute to this process in Fig.2,
both of which contribute equally. In each diagram, there are two Lorentz
frames involved, the rest frame of the initial meson, $H$, and the rest frame
of final meson, $H^{\prime}$.

As in previous section, we again write the total momentum and the internal
momentum of initial hadron, $H(=A)$ as $P$, and $q$, while the corresponding
momentum variables for final hadron, $H^{\prime}(=V)$ as, $P^{\prime}$, and
$q^{\prime}$. And let $k$, and $\epsilon^{\lambda^{\prime}}$ be momentum and
polarization vectors of emitted photon, while $\epsilon^{\lambda}$ be the
polarization vector of final vector meson. Thus if $p_{1,2}$, and
$p^{\prime}_{1,2}$ are the momenta of the two quarks in initial and final
hadron respectively, then, the kinematical relations between the momenta of
initial and final hadrons can be expressed as, $q^{\prime}=q-\frac{1}{2}k$
(see [9] for details). To simplify calculations, we again work in rest frame
of initial hadron. The amplitude for the decay process, $A\rightarrow
V\gamma$, is two times the amplitude for Diagram 1. Following Section 2, we
write the invariant amplitude for this process as,

$M_{fi}=-2ie_{Q}\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\frac{1}{M^{2}}Tr\bigg{[}\alpha_{1}{\not}P\overline{\psi}_{V}^{++}(\hat{q}^{\prime}){\not}\epsilon^{\prime\prime}\psi_{A}^{++}(\hat{q})\bigg{]},$
(27)

where the $++$ components of vector and axial vector meson wave functions are
given as in Eq.(10) of the previous section, except that we now take
$(\epsilon,P,\hat{q})$ as the polarization vector and momentum variables of
initial axial meson, while,
($\epsilon^{\prime\prime},P^{\prime},\hat{q}^{\prime}$) as the corresponding
variables of the final vector meson.

Now trace over the Dirac gamma matrices can be expressed as,

$TR=\alpha_{1}[X_{1}\epsilon_{\mu\nu\alpha\beta}P_{\mu}\epsilon_{\nu}\epsilon^{\prime}_{\alpha}\epsilon^{\prime\prime}_{\beta}+X_{2}\epsilon_{\mu\nu\alpha\beta}k_{\mu}\epsilon_{\nu}\epsilon^{\prime}_{\alpha}\epsilon^{\prime\prime}_{\beta}+X_{3}\epsilon_{\mu\nu\alpha\beta}P_{\mu}k_{\nu}\epsilon^{\prime\prime}_{\alpha}\epsilon_{\beta}+X_{4}\epsilon_{\mu\nu\alpha\beta}P_{\mu}k_{\nu}\epsilon^{\prime}_{\alpha}\epsilon^{\prime\prime}_{\beta}+X_{5}\epsilon_{\mu\nu\alpha\beta}P_{\mu}k_{\nu}\epsilon_{\alpha}\epsilon^{\prime\prime}_{\beta}]$
(28)

where,

$\displaystyle
X_{1}=-a_{1}b_{1}+M^{2}a_{2}b_{2}+(\hat{m}_{2}M^{2}-\frac{P^{\prime}.P}{P^{2}}\hat{m}_{2}M^{2})a_{3}b_{1}-M^{2}\hat{q}^{2}a_{3}b_{3}+\hat{m}_{2}(1-\frac{P^{\prime}.P}{P^{2}})a_{4}b_{2}-\hat{q}^{2}a_{4}b_{4}+$
$\displaystyle(\hat{m}_{2}M^{2}(1-\frac{P^{\prime}.P}{P2})-\hat{q}^{2}M^{2})a_{7}b_{1},$
$\displaystyle
X_{2}=-M^{2}a_{2}b_{2}-\frac{P^{\prime}.P}{P^{2}}\hat{m}_{2}M^{2}a_{3}b_{1}+M^{2}\hat{q}^{2}a_{3}b_{3}-M^{2}\hat{m}_{2}a_{4}b_{2}+M^{2}\hat{m}_{2}(-2+\frac{P^{\prime}.P}{P^{2}})a_{7}b_{1}+M^{2}\hat{q}^{2}a_{7}b_{3},$
$\displaystyle X_{3}=2\hat{m}_{2}[2-\frac{P^{\prime}.P}{P^{2}}]a_{3}b_{1},$
$\displaystyle
X_{4}=\hat{q}^{2}a_{3}b_{5}+M^{2}\hat{m}_{2}\frac{P^{\prime}.P}{P^{2}}a_{3}b_{6}-\hat{q}^{2}a_{7}b_{5}-M^{2}\hat{m}_{2}(1-\frac{P^{\prime}.P}{P^{2}})a_{7}b_{6},$
$\displaystyle
X_{5}=-2M^{2}\hat{m}_{2}a_{6}b_{2}+\hat{q}^{2}a_{6}b_{4}+\hat{m}_{2}a_{8}b_{1},$
(29)

where the coefficients, $a_{1},...,a_{8}$, and $b_{1},...,b_{7}$ for
$A\rightarrow V\gamma$ are of the same form as in Eqs.(11)-(12) for
$V\rightarrow A\gamma$ in previous section, except that, we have to
interchange $M\rightleftharpoons
M^{\prime},\omega\rightleftharpoons\omega^{\prime},\hat{q}\rightleftharpoons\hat{q}^{\prime}$.

Thus, we can write $M_{fi}$ as,

$\displaystyle
M_{fi}=-2ie_{Q}\frac{N_{A}N_{V}}{M^{2}}[F_{1}\epsilon_{\mu\nu\alpha\beta}P_{\mu}\epsilon_{\nu}\epsilon^{\prime}_{\alpha}\epsilon^{\prime\prime}_{\beta}+F_{2}\epsilon_{\mu\nu\alpha\beta}k_{\mu}\epsilon_{\nu}\epsilon^{\prime}_{\alpha}\epsilon^{\prime\prime}_{\beta}+(P.\epsilon^{\prime})F_{3}\epsilon_{\mu\nu\alpha\beta}P_{\mu}k_{\nu}\epsilon^{\prime\prime}_{\alpha}\epsilon_{\beta}+(I.\epsilon)F_{4}\epsilon_{\mu\nu\alpha\beta}P_{\mu}k_{\nu}\epsilon^{\prime}_{\alpha}\epsilon^{\prime\prime}_{\beta}+$
$\displaystyle(I.\epsilon^{\prime})F_{5}\epsilon_{\mu\nu\alpha\beta}P_{\mu}k_{\nu}\epsilon_{\alpha}\epsilon^{\prime\prime}_{\beta}];$
$\displaystyle
F_{1}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\alpha_{1}X_{1}(\hat{q})\phi_{A}(\hat{q})\phi_{V}(\hat{q})$
$\displaystyle
F_{2}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\alpha_{1}X_{2}(\hat{q})\phi_{A}(\hat{q})\phi_{V}(\hat{q})$
$\displaystyle
F_{3}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\alpha_{1}X_{3}\phi_{A}(\hat{q})\phi_{V}(\hat{q})$
$\displaystyle
F_{4}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\alpha_{1}X_{4}\phi_{A}(\hat{q})\phi_{V}(\hat{q})$
$\displaystyle
F_{5}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\alpha_{1}X_{5}\phi_{A}(\hat{q})\phi_{V}(\hat{q}).$
(30)

As mentioned earlier, since $\hat{q}$ is an effective 3D variable, we have
written
$\hat{q}.\epsilon^{\prime}=\hat{q}_{\nu}\epsilon^{\prime}_{\nu}=|\hat{q}|(I.\epsilon^{\prime})$,
where, $I_{\mu}$ is a unit vector along the direction of $\hat{q}_{\nu}$, and
is expressed as, $I=\frac{\hat{q}}{|\hat{q}|}$, where $|\hat{q}|$ is a
Lorentz-invariant variable [18, 26, 22]. Similarly we write
$\hat{q}.\epsilon^{\prime\prime}=|\hat{q}|(I.\epsilon^{\prime\prime})$, and
$\hat{q}.\epsilon=|\hat{q}|(I.\epsilon)$.

In the above expression for $M_{fi}$ the term associated with the form factor,
$F_{3}$ vanishes on account of $P.\epsilon^{\prime}=0$, since
$P=(\overrightarrow{0},iM)$ is momentum of initial axial vector meson at rest,
while $\epsilon^{\prime}=(\overrightarrow{\epsilon^{\prime}},i0)$ is the
polarization vector of the emitted photon, that is transversely polarized.
Further, we wish to mention that the amplitudes, $F_{4}$ and $F_{5}$
contribute to M2 transitions, and make sub-dominant contributions to decay
width, and are ignored. And the dominant contributions arise from E1
transitions that come from the amplitudes, $F_{1}$ and $F_{2}$, in terms of
which we can write the amplitude, $M_{fi}$ as,

$M_{fi}=-2ie_{Q}\frac{N_{A}N_{V}}{M^{2}}[F_{1}\epsilon_{\mu\nu\alpha\beta}\epsilon^{\prime\prime}_{\mu}\epsilon^{\prime}_{\nu}\epsilon_{\alpha}P_{\beta}+F_{2}\epsilon_{\mu\nu\alpha\beta}\epsilon^{\prime\prime}_{\mu}\epsilon^{\prime}_{\nu}\epsilon_{\alpha}k_{\beta}]$
(31)

Thus, we can write, $M_{fi}$ as,

$\displaystyle
M_{fi}=\epsilon^{\prime\prime}_{\mu}\epsilon^{\prime}_{\nu}\epsilon_{\alpha}M_{\mu\nu\alpha};$
$\displaystyle
M_{\mu\nu\alpha}=F_{1}\epsilon_{\mu\nu\alpha\beta}P_{\beta}+F_{2}\epsilon_{\mu\nu\alpha\beta}k_{\beta}.$
(32)

Now, electromagnetic gauge invariance demands,

$k_{\mu}M_{\mu\nu\alpha}=0,$ (33)

which leads to the form factor, $F_{1}=0$. Thus we can express the amplitude,
$M_{fi}$, that corresponds to E1 transition and makes dominant contribution
as,

$\displaystyle
M_{fi}=-2ie_{Q}\frac{N_{A}N_{V}}{M^{2}}[F_{2}\epsilon_{\mu\nu\alpha\beta}\epsilon^{\prime\prime}_{\mu}\epsilon^{\prime}_{\nu}\epsilon_{\alpha}k_{\beta}]$
$\displaystyle
F_{2}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\alpha_{1}X_{2}(\hat{q})\phi_{A}(\hat{q})\phi_{V}(\hat{q}).$
(34)

Here, $X_{2}$ appearing in the expression for $F_{2}$ is given in Eqs.(21).
Now, to calculate the decay widths, we need to calculate the spin averaged
amplitude modulus square, $|\overline{M}_{fi}|^{2}$, where
$|\overline{M}_{fi}|^{2}=\frac{1}{2j+1}\sum_{\lambda,\lambda^{\prime}}|{M}_{fi}|^{2}$,
where we average over the initial polarization states $\lambda$ of A-meson,
and sum over the final polarization $\lambda^{\prime}$ of photon. We make use
of the normalizations,
$\Sigma_{\lambda}\epsilon_{\mu}^{\lambda}\epsilon_{\nu}^{\lambda}=\frac{1}{3}(\delta_{\mu\nu}+\frac{P_{\mu}P_{\nu}}{M^{2}})$
for vector meson, and
$\Sigma_{\lambda^{\prime}}\epsilon_{\mu}^{\lambda^{\prime}}\epsilon_{\nu}^{\lambda^{\prime}}=\delta_{\mu\nu}$,
for the emitted photon, with $M_{fi}$ taken from the previous equation, we
get,
$\sum_{\lambda^{\prime}}\sum_{\lambda}|\epsilon^{\lambda^{\prime}}.\epsilon^{\lambda}|^{2}=1$.

We obtain,

$|\overline{M}_{fi}|^{2}=\frac{4e_{Q}^{2}N_{A}^{2}N_{V}^{2}}{9M^{4}}|F_{2}|^{2}\bigg{[}\frac{(P^{\prime}.k)^{2}}{M^{\prime
2}}+\frac{(P.k)^{2}}{M^{2}}-\frac{(P^{\prime}.k)(P^{\prime}.P)(P.k)}{M^{2}M^{\prime
2}}\bigg{]},$ (35)

where, the dot products of momenta calculated in the rest frame of initial
vector meson are: $P^{\prime}.P=-\frac{(M^{2}+M^{\prime 2})}{2}$,
$P.k=-\frac{M^{2}-M^{\prime 2}}{2}$, $P^{\prime}.k=-\frac{M^{2}-M^{\prime
2}}{2}$. Thus we notice that the general expression for $|\bar{M}_{fi}|^{2}|$
for the process $1^{++}\rightarrow 1^{--}\gamma$ is the same as the
corresponding expression for the process $1^{--}\rightarrow 1^{++}\gamma$ in
Eq.(23), and (24), except for the difference in the form factors $F_{2}$, and
$G_{2}$ respectively.

The decay width of the process ($A^{+}\rightarrow V\gamma$) in the rest frame
of the initial meson is expressed as in Eq.(25), and we make use of the fact
that modulus of the momentum of the emitted vector meson can be expressed in
terms of masses of particles as,
$|\overrightarrow{P^{\prime}}|=|\overrightarrow{k}|=\omega_{k}=\frac{1}{2M}(M^{2}-M^{\prime
2})$ (where $M/M^{\prime}$ is mass of initial axial/final vector meson, where,
$\omega_{k}$ is the kinematically allowed energy of the emitted photon.

|  BSE-CIA |  Expt.[1] |  QM[11] |  LFM[10] |  BSE[12]
---|---|---|---|---|---
$\Gamma_{\chi_{c1}(1P)\rightarrow\gamma J/\psi(1S)}$ |  291.11 |  288$\pm$11.52 |  275 |  324$\pm$20 |  306
$\Gamma_{\chi_{c1}(2P)\rightarrow\gamma\psi(2S)}$ |  69.703 |  53.55$\pm$23.8 |  |  |  146.0
$\Gamma_{\chi_{c1}(2P)\rightarrow\gamma J\psi(1S)}$ |  18.959 |  9.52$\pm$0.84 |  |  |  33.0
Table 3: Radiative decay widths of $1^{++}$ charmonia (in keV) in present
model along with experimental data and results of other models.

##  4\. Conclusion

In present work, we have studied the radiative E1 transitions
${}^{3}P_{1}\rightarrow^{3}S_{1}\gamma$, and
${}^{3}S_{1}\rightarrow^{3}P_{1}\gamma$ in the framework of Bethe-Salpeter
equation and calculated their decay widths. The 3D Salpeter wave functions
$\phi_{V}(\hat{q})$, and $\phi_{A}(\hat{q})$ used for calculation of the
transition amplitudes, were derived as analytic solutions of mass spectral
equations for vector and axial quarkonia[25, 8] in an approximate harmonic
oscillator basis. Our results of decay widths are compared with data and other
models. In our framework, ratio of decay widths,
$\frac{\Gamma_{\chi_{c1}(2P)\rightarrow\gamma\psi(2S)}}{\Gamma_{\chi_{c1}(2P)\rightarrow\gamma
J\psi(1S)}}=3.676$ (c.f. Expt.= 5.625[1]). The corresponding ratio of these
decay widths in [12] is 4.424.

## References

  * [1] R. L. Workman et al. (Particle Data Group), PTEP2022, 083C01 (2022).
  * [2] P.A.Zyla et al., (Particle Data Group), Prog. Theo. Expt. Phys. 2020, 083C01 (2020).
  * [3] K.A.Olive et al., (Particle Data Group), Chin. Phys. C38, 090001 (2014).
  * [4] V. Bhardwaj et al. (Belle), Phys. Rev. Lett. 111, 032001 2013.
  * [5] T.E.Coan et al., hep-ex/0509030.
  * [6] M. Ablikim et al. (BESIII), Phys. Rev. Lett. 115, 011803(2015); arXiv:1503.08203 [hep-ex].
  * [7] M.Ablikim et al., Phys. Rev. D97, 032008 (2018).
  * [8] V.Guleria, S.Bhatnagar, Intl. J. Theor. Phys. 60, 3143 (2021).
  * [9] V.Guleria, E.Gebrehana, S.Bhatnagar, Phys. Rev. D104, 094045 (2021).
  * [10] Y.L.Shi, arxiv:1611.03712[hep-ph].
  * [11] Wei-Jun Deng, Li-Ye Xiao, Long-Cheng Gui, and Xian-Hui Zhong, arxiv: 1510.08269[hep-ph].
  * [12] T.H. Wang, G.L. Wang, Phys. Lett. B697, 233 (2011).
  * [13] D. Ebert, R. N. Faustov and V. O. Galkin, Phys. Rev. D67, 014027 (2003).
  * [14] S. Godfrey and K. Moats, Phys. Rev. D 92, 054034(2015).
  * [15] J. Segovia, P. G. Ortega, D. R. Entem and F. Fernandez, Phys. Rev. D 93, 074027 (2016).
  * [16] J.K. He, C. J. Fan, Phys. Rev. D103, 114006 (2021).
  * [17] J.K. He, Y. D. Yang, Nucl. Phys. B943, 114627 (2019).
  * [18] S.Bhatnagar, E.Gebrehana, Phys. Rev. D102, 094024 (2020).
  * [19] C. H. L. Smith, Ann. Phys. (N.Y.) 53, 521 (1969).
  * [20] R. Alkofer, L.V.Smekel, Phys. Rep. 353, 281 (2002).
  * [21] S.Bhatnagar, L. Alemu, Phys. Rev. D97, 034021 (2018).
  * [22] S.Bhatnagar, V.Guleria, Nucl. Phys. A1041, 122783 (2024).
  * [23] M.Narang, S.Bhatnagar, Few-Body Syst. (2023).
  * [24] H.Negash, S.Bhatnagar, Intl. J. Mod. Phys. E25, 1650059 (2016).
  * [25] E.Gebrehana, S.Bhatnagar, H.Negash, Phys. Rev. D100, 054034 (2019).
  * [26] C.H.Chang, J.K.Chen, G.L.Wang, Commun. Theo. Phys.(Beijing) 46, 467 (2006).
  * [27] Y.L.Shi, Eur. Phys. J. C77, 253 (2017).

50