#  Pushing the Limits of Quantum Computing for Simulating PFAS Chemistry

Emil Dimitrov, Goar Sánchez-Sanz, James Nelson, Lee O’Riordan, Myles Doyle,
Seán Courtney, Venkatesh Kannan
Irish Centre for High-End Computing (ICHEC)
{emil.dimitrov<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
&Hassan Naseri, Alberto García García, James Tricker , Marisa Faraggi
Accenture
{hassan.naseri, a.b.garcia.garcia, james.tricker,
<EMAIL_ADDRESS>
&Joshua Goings, Luning Zhao
IonQ
{goings<EMAIL_ADDRESS>Corresponding author.

###### Abstract

Accurate and scalable methods for computational quantum chemistry can
accelerate research and development in many fields, ranging from drug
discovery to advanced material design. Solving the electronic Schrödinger
equation is the core problem of computational chemistry. However, the
combinatorial complexity of this problem makes it intractable to find exact
solutions, except for very small systems. The idea of quantum computing
originated from this computational challenge in simulating quantum-mechanics.
We propose an end-to-end quantum chemistry pipeline based on the variational
quantum eigensolver (VQE) algorithm and integrated with both HPC-based
simulators and a trapped-ion quantum computer. Our platform orchestrates
hundreds of simulation jobs on compute resources to efficiently complete a set
of _ab initio_ chemistry experiments with a wide range of parameterization.

Per- and poly-fluoroalkyl substances (PFAS) are a large family of human-made
chemicals that pose a major environmental and health issue globally. Our
simulations includes breaking a Carbon-Fluorine bond in trifluoroacetic acid
(TFA), a common PFAS chemical. This is a common pathway towards destruction
and removal of PFAS. Molecules are modeled on both a quantum simulator and a
trapped-ion quantum computer, specifically IonQ Aria. Using basic error
mitigation techniques, the 11-qubit TFA model (56 entangling gates) on IonQ
Aria yields near-quantitative results with milli-Hartree accuracy. Our novel
results show the current state and future projections for quantum computing in
solving the electronic structure problem, push the boundaries for the VQE
algorithm and quantum computers, and facilitates development of quantum
chemistry workflows.

###### Index Terms:

Variational Quantum Eigensolver (VQE) molecular simulation scalability bond
stretching PFAS

##  I Introduction

Per- and poly-fluoroalkyl substances (PFAS) are a large family of human-made
molecules that are extensively used for their water and oil repellent, heat
resistant, and non-stick properties in areas such as clothing, adhesives, fire
retardants, food packaging, cooking utensils, and paints. PFAS are highly
resistant to bio-degradation, and have high mobility in water and soils. PFAS
can contaminate drinking water [1] and bio-accumulate in wildlife. PFAS
exposure has detrimental health impacts on humans including cancer, thyroid
diseases, hormonal changes, and reduced vaccine response, among others [2, 3].
Therefore, remediation of PFAS from the ecosystem is crucial and has been
widely studied both from experimental and computational/theoretical points of
view [4].

Several techniques have been developed to physically remove PFAS from drinking
water, for example, absorption of PFAS in active carbon, water filtration, or
using ion exchange techniques[5, 6, 7]. However, once the molecules are
isolated, they must be still degraded and destroyed.

PFAS family presents an unusual thermal stability, mainly due to the large
electro-negativity of fluorine, which forms strong C-F bonds. These bonds
result in high surface activity and thermal resistance, as well as resistance
to many forms of physical, chemical, and biological degradation. From chemical
point of view, PFAS degradation presents a huge technical challenge, which
depends on the type of PFAS, the length of the chains, and additional
substituents within the compound. As an example, Figure 1 illustrates
perfluorocarbonylic acids, PFCAs (CnF2n+1-COO-), and perfluorosulfonic acids,
PFSAs (CnF2n+1-SO${}_{3}^{-}$) which are less stable than their telomeric
counterparts, TFCA (CnF2n+1-CH2-CH2-COO-) and TFSA
(CnF2n+1-CH2-CH2-SO${}_{3}^{-}$) [8].

Figure 1: Representative structures of PFAS molecules that are common in the
environment: PFCA (perfluorocarbonylic acid), PFSA (perfluorosulfonic acid),
TFCA (fluorotelomer carboxylic acid) and TFSA (trifluoromethanesulfonic acid)
with n=6; n indicating the variable length of CnF2n+1

In the last few years, there has been significant development in PFAS
destruction technologies. Examples include wet oxidation using heat activated
oxidants [9, 10] and ion exchange using plasma reactors [11, 12] among other
methods. Most of these techniques involve the use of strong oxidative species
(mostly a variety of radicals) and/or reductive species (e.g. hydrated
electron e${}^{-}_{aq}$ produced by using UV light) [13]. The main objective
of such techniques is to degrade PFAS molecules into smaller and treatable
units. From the structural perspective, two main strategies can be considered:
defluorination (dissociation of C-F bonds) and degradation of the backbone
chain into smaller parts. However, these methods are even more challenging
since C-F is among one of the strongest bonds in chemistry. Therefore, an
attractive and potentially useful challenge is to computationally investigate
PFAS degradation through defluorination and breaking C-F bonds.

For this computational chemistry challenge, there is a large number of
possible experiments involving the dissociation of a single fluorine atom
depending on the size of the PFAS, the number of C-F bonds and their
positions. Computational coupled-cluster-based methods such as CCSD or CCSD(T)
are preferable to achieve chemical accuracy, since these methods have been
established as reliable and accurate for the computation of chemical energies.
CCSD(T) is often referred to as the “gold-standard” of electronic structure.
However, these methods scale at $N^{6}$ and $N^{7}$ with the number of basis
functions $N$, leading to high computational cost. As the number of basis
functions grows with the number of atoms in a molecule, the computational
expense rapidly grows with atom count. Hence, applying these high-accuracy
methods to large molecules such as PFAS is impractical.

In this project we focus on the molecular simulation and analysis of the
electronic and chemical properties of PFAS-like molecular structures.
Currently, high-performance computing (HPC) models and tools are used for
simulation and analysis of molecular structures and their electronic/chemical
properties. However, current HPC simulations are computationally expensive and
not particularly scalable, which limits the accuracy, size, and number of
compounds subject to study. Approximations are necessary to reduce the
computational cost and also to map molecular and electronic properties, which
are inherently quantum, for execution on classical computation platforms.

With the advent of quantum computing hardware and algorithms, there has been
interest in applying quantum computational methods to the simulation and
analysis of the electronic structure systems. The scope of this project is to
develop a high-performance hybrid quantum computational workflow with the aim
of performing scalable and repeatable experiments using Variational Quantum
Eigensolver (VQE) -based methods[14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25]. We demonstrate the development of an end-to-end workflow to run hybrid
high-performance quantum computational chemistry experiments with a wide range
of choices in the parameter space (e.g. molecule types, active spaces, qubit
reduction methods, chemistry parameters, and VQE-related parameters). This
workflow can interface with quantum computing software simulators as well as
quantum computing hardware. The pipeline is used to compute the energies in a
bond-breaking experiment for PFAS-related molecules both in simulations and on
a real quantum hardware.

The rest of this manuscript is organized as follows. Section 2 describes the
implementation of the workflow in details. Section 3 contains the results,
which are divided into two subsections: the results obtained on a quantum
simulator, and the results from a real quantum computer. Finally, we summarize
the work in Section 4.

##  II Implementation

As outlined in the introduction, the dissociation energy of C-F bonds is
studied in this paper. Breaking a C-F bond is a fundamental step in the
degradation process of PFAS. To simulate bond dissociation, a target molecule
is selected and a dictionary of geometries is generated. Then, the position of
one of the atoms, usually the outermost (e.g. fluorine atom in TFA), is
shifted in steps to simulate bond breaking and evaluate the energy surface at
different distances. The electronic structure problem for these geometries are
first solved at the Hartree-Fock self-consistent field level. From these
results the one- and two-body electronic tensors are obtained to form the
fermionic Hamiltonian, which are then solved using VQE to obtain an estimate
of the ground state energy. The fermionic Hamiltonian must be mapped to qubit
operations to be executed on a quantum hardware. This transformation can be
efficiently performed classically using, e.g., the Jordan-Wigner, Bravyi-
Kitaev, and other transformations. Once the Hamiltonian is represented by a
set of qubit operations, they can be executed using VQE given a suitable trial
circuit on a real quantum computing hardware or a software simulator. In order
to orchestrate many such experiments and run them on different backend systems
(both quantum simulator and real quantum computer), the QPFAS workflow is
developed, which is described in the following.

###  II-A Workflow for scalable molecular simulation

QPFAS is a computational orchestration workflow that is built around the VQE
algorithm, and can be executed on gate model quantum computing hardware as
well as software simulators. QPFAS is integrated and tested with multiple
quantum simulators. The main simulators results presented in this paper are
from the QULACS [26] simulator. The results of the real quantum computer are
obtained using the IonQ computer based on trapped ion technology [27]. Figure
2 outlines the QPFAS workflow. The top row shows the parameters that are to be
specified at each step, while different modules are described in the second
row of the figure. The workflow employs the Python package Tequila [28] to
create the quantum circuits and implement the VQE algorithm.

Figure 2: Scheme of the QPFAS calculation workflow. The parameters to be used
in each stage are detailed in the upper part of the figure. The workflow steps
are represented in boxes in the lower part of the figure.

In the first step of the QPFAS workflow, a molecule object was instantiated to
store the molecule’s geometry, multiplicity, charge, basis set, and active
space. This step provides the information required to perform classical
chemistry calculations and fully specifies the chemistry problem. Next, the
molecule is transformed in Tequila format[28] specifying the qubit
transformation to be used. Both objects (model formats) are used in step three
to create the ansatz circuit and the qubit Hamiltonian. After creating these
two, it is possible to run VQE using a classic optimizer in the last step of
the workflow, where the number of samples and the backend (quantum computing
software simulator or hardware) need to be also specified as parameters.

The QPFAS workflow yields ground state energies of molecules. These VQE-
derived energies are benchmarked against those from classical methods such as
Coupled Cluster Single-Double (CCSD) or Full Configuration Interaction (FCI)
with Tequila. Here, it is important to mention that certain choices of active
space are not compatible with Tequila’s classical solver, so in these cases
the Python-based Simulations of Chemistry Framework (PySCF) [29] is used for
the classical calculations.

####  II-A1 VQE approaches in QPFAS

VQE [30] is a hybrid classical-quantum algorithm that is suited for noisy
intermediate-scale quantum (NISQ) computing devices. It is an approximate
method to solve optimization problems, more specifically the minimum
eigenvalue problem. Due to the large resource requirements of other quantum
algorithms, such as quantum phase estimation,[31] for computational chemistry
problems, VQE is likely to remain the de facto method in this field for
foreseeable future. In this work, in addition to the standard VQE, ADAPT-VQE
and VQE by tapering are also considered as they have certain advantages for
computational chemistry problems. The standard VQE and its variations are
illustrated in Figure 3, and desctibed in the following.

(a)

(b)

(c)

Figure 3:  VQE enabled in PQFAS workflow: (a) Standard VQE, (b) ADAPT-VQE
which grows one operator at a time in a manner set by the molecule being
simulated, (c) VQE with tapering to find symmetries in the qubit Hamiltonian.
Standard and ADAPT-VQE run a self-consistent loop placed inside the red box in
the figure.

##### Standard VQE algorithm

(Figure 3(a)) VQE is based on the _variational principle_ , which states that
the smallest valid energy solution for a state will always be an upper bound
to the energy of the ground state. In the VQE algorithm, the molecular wave
function is encoded as a parameterized circuit called the ansatz, $U(\theta)$.
A Hamiltonian, $H$, is transformed from fermionic operators to qubit
operators. To find the ground-state energy, the parameters are iteratively
updated to minimize the expectation value of the Hamiltonian (1), until
convergence is reached.

$E_{\text{ground-state}}=\min\langle 0|U(\theta)^{\dagger}HU(\theta)|0\rangle$
(1)

where $|0\rangle$ is the reference state the circuit acts upon usually the
Hartree-Fock ground state of the system. To generate the molecular
Hamiltonian, the geometry of the molecule and the basis set to be used must be
specified, and the necessary integrals to be computed. The geometry of the
molecules is stored in an xyz file that contains the specification and
coordinates of each atom. The second quantized Hamiltonian can then be written
as:

$H=h_{0}+\sum_{ij}h_{ij}a_{i}^{\dagger}a_{j}+\sum_{ijkl}h_{ijkl}a_{i}^{\dagger}a_{j}^{\dagger}a_{k}a_{l}$
(2)

where $h_{0}$ is a constant, $h_{ij}$ and $h_{ijkl}$ are the single and two-
body integrals respectively, and $a_{i}^{\dagger}$ is the creation operator
for orbital $i$ and $a$ is the annihilation operator for orbital $i$. To
implement the Hamiltonian on a quantum computing hardware or software
simulator, the fermionic Hamiltonian (2), should be mapped to a qubit
Hamiltonian. This is done with both Jordan-Wigner and Bravyi-Kitaev
transformations in the workflow.

##### ADAPT-VQE algorithm

(Figure 3(b)) While being based on the standard VQE algorithm, the ADAPT-VQE
uses an adaptive ansatz as introduced in [32]. A single operator is added to
the ansatz (which initially is the identity circuit) from a pool of candidate
operators. Operators are added until convergence is reached. Convergence is
deemed as the gradient of energy with respect to the pool of operators. Each
time an operator is added to the ansatz, the VQE algorithm is run to find the
optimal parameters.

##### VQE algorithm with tapering

(Figure 3(c)) This variation uses the tapering algorithm [33] to remove qubits
before running the VQE. The tapering algorithm is based on finding symmetries
in the qubit Hamiltonian, where symmetry is defined by terms that commute with
the Hamiltonian. For each independent symmetry found, one qubit can be removed
through a transformation of the Hamiltonian. While this method reduces the
number of qubits used, there are two downsides to be considered. In this VQE
variation the transformation of the Hamiltonian does not preserve the number
of particles and thus only the hardware ansatz can be used as it is the only
ansatz that does not conserve particle number. Additionally, when removing a
qubit, it is essential to know what vector to project it into ($0$ or $1$).
Thus, if $n$ qubits are removed, then the VQE must be run $2^{n}$ times,
resulting in the lowest energy value being the final result for the ground
state.

####  II-A2 Ansatz calculation methods in QPFAS

In order to run VQE, it is necessary to define an ansatz. There are generally
two paradigms for VQE ansatze design: a) chemistry-inspired ansatze which are
designed to capture the form of the molecular electronic wave function, and b)
hardware efficient ansatze (HEAs) which attempt to offer flexible and
expressive circuits that are shallower and more efficient to run on noisy
hardware. Early demonstrations of VQE on quantum hardware utilized the unitary
coupled cluster with singles and doubles (UCCSD) ansatz [30, 15, 16, 17, 18].
Despite its success in treating strongly correlated systems, UCCSD becomes
impractical on NISQ quantum computers due to its rapid increase in entangling
gates. Consequently, hardware efficient ansatze (HEAs) have gained attention
for requiring significantly shallower circuits. While HEAs have been used to
simulate various molecules on superconducting quantum computers [34, 20, 21,
22, 23], HEAs struggle to deal with noise and the so-called barren plateaus in
optimization [35]. Moreover, they too are subject to the measurement
bottleneck [36]. To maximize efficiency while maintaining chemical fidelity,
we use a modified ansatz derived from the UCCSD ansatz: the unitary pair CCD
(pUCCD) ansatz [37, 38].

In total, the QPFAS workflow allows for the use of five ansatz design methods.
Three of them are chemistry-inspired ansatze based on Unitary Coupled-Cluster
Single and Double excitations and two are hardware efficient ansatze. These
five ansatze well represent the current state of the art of VQE for chemistry
problems at the time this project. The available ansatze in QPFAS, including
pUCCD, are explained in the following.

  * •

UCCSD and k-UpCCGSD: These two ansatze are chemistry inspired. In UCCSD ansatz
only single and double excitations are considered, which preserve the Hartree-
Fock ground state. It is the original ansatz proposed for VQE in chemistry in
2014 [30]. For k-UpCCGSD ansatz all single excitations are taken into account
(even ones that destroy the Hartree Fock ground-state), but only double
excitations that move a pair of electrons from one spatial orbital to another
are allowed. This is a more efficient approximation to UCCSD, which has never
been tested on a real quantum computer – simulation only. Both these ansatze
are of the form, $U=e^{T-T^{\dagger}}$, where $T$ is an operator generating
excitations, and contain a depth parameter. In UCCSD this corresponds to the
number of trotter steps, while in k-UpCCGSD the depth ($k$) is the number of
times the set of excitations are repeated (with different angles each time).

  * •

pUCCD: This is also a chemistry inspired ansatz derived from the UCCSD ansatz,
but only contains double excitations [37, 38]. The pUCCD ansatz has several
appealing characteristics for use on near-term hardware [39]. First, it allows
for a hard-core boson representation and requires half the number of qubits
compared to other UCC ansatze. Second, at most three circuits are needed for
energy expectation value computation, eliminating the measurement bottleneck
for this ansatz (measurement is $\mathcal{O}(1)$). Finally, it has entangling
gate requirements that goes as $\mathcal{O}(N^{2})$ and entangling gate depth
that goes as $\mathcal{O}(N)$ [24].

  * •

“hardware” and “hardware conserving”: These two ansatze consider quantum
computing hardware. hardware ansatz [34] and hardware conserving ansatz [40]
both consist of repeating layers of rotation gates (acting on individual
qubits) followed by an entangling layer. In the case of hardware ansatz the
circuit acts on the null state (all qubits set to zero) and the entangling
layer is a set of CNOT gates. The hardware conserving ansatz acts on the
Hartree-Fock ground state and the entangling layer is designed to conserve the
number of particles.

####  II-A3 Orbital reduction methods in QPFAS

In quantum computer the number of qubits required to run a circuit is directly
related with the number of spin orbitals of the problem. Hence, using a subset
of the orbitals instead of the entire set requires a significantly fewer
qubits for implementation. The QPFAS workflow supports two automatic methods
to find such a subset of orbitals: frozen core approximation and Natural
Orbital Occupation Numbers (NOONs). It is also possible to provide a manual
input. These are referred to as orbital reduction methods. For background, the
orbitals are categorized as occupied, frozen, or virtual. Usually, the inner
occupied orbitals are called core orbitals. These can be assumed to be frozen,
i.e., not participant in the excitations. The unoccupied orbitals are called
virtual. The remaining occupied orbitals and the virtual orbitals are referred
to as the active set. Thus the Hamiltonian is expressed as

$H=\sum_{i\in F}h_{ii}+\sum_{ij\in A}\left(h_{ij}+\sum_{k\in
F}(h_{kijk}-h_{kikj}+h_{ikkj}-h_{ikjk})\right)c_{i}^{\dagger}c_{j}+\sum_{ijkl\in
A}h_{ijkl}c_{i}^{\dagger}c_{j}^{\dagger}c_{k}c_{l}$ (3)

where $F$ is the set of frozen orbitals and $A$ is the active set.

The frozen core approximation treats the lowest energy orbitals as being
occupied, assuming that the core electrons are not correlated. The number of
orbitals to freeze is determined based on the atomic species of the
molecule111ORCA Input Library, Frozen core calculations.
https://sites.google.com/site/orcainputlibrary/frozen-core-calculations. For
Natural Orbital Occupation Numbers (NOONs), the single particle density
matrix, $\rho(x,y)$ is defined by the set of basis functions
$\\{\phi_{i}(x)\\}$, given by

$\rho(x,y)=\sum_{ij}\phi_{i}(x)\gamma_{ij}\phi_{j}(y)^{*}.$ (4)

Since $\gamma_{ij}$ is a hermitian matrix we can diagonalize it, where the
resulting diagonal matrix contains the NOONs. For Hartree-Fock, the NOONs are
either zero or one. By computing the NOONs with a higher-order method (such as
MP2) we get a non-trivial result. If the NOONs for a given orbital is close to
zero (or two) within a specified tolerance, then the orbital is considered to
be unoccupied (doubly occupied).

##  III Experiments and results

Due to the limitations of quantum computing simulators and hardware, we study
the bond-breaking for smaller PFAS-like molecules suach as CH3F. This lays the
groundwork to study larger molecules with C-F bonds in the future. Carrying
out these experiments with QPFAS workflow allows integrating HPC models and
IonQ real quantum hardware. Another objective is to explore the current
computational limits for quantum chemistry calculations using an end-to-end
hybrid workflow. Hence, the following experiments seek to explore three main
characteristics:

  * •

Scalability of the problem in terms of the size of the molecules;

  * •

Functionality of the QPFAS workflow with the largest bond-breaking experiment
that can fit into our HPC or a real quantum computer;

  * •

VQE accuracy and its performance under a wide range of parameterizations.

The experiments executed on the HPC simulator serve to refine the parameters
for running a bond-breaking experiments in the real quantum computer. The
experiment and the results are divided in two subsection: first, calculations
in QULACS simulator on ICHEC HPC cluster, and second, calculations in the IonQ
quantum computer.

###  III-A Experiments using QULACS software simulator

The first part of the experimentation is focused on gaining knowledge on
scalability and parameter optimization running calculations on a group of
eleven molecules with increasing size. The second part is dedicated to explore
energy dissociation only on F2 and CH3F. Then the bond stretching calculations
of these two molecules are executed, as they are the biggest molecules we can
execute in the ICHEC’s HPC in a reasonable amount of time (less than 24
hours). F2 molecule is selected because it is one of the largest with an F-F
bond that we could run on the quantum emulator. CH3F is the largest molecule
that can be executed with a C-F bond in the quantum emulator.

####  III-A1 Scalability and parameterization study

VQE is a heuristic quantum algorithm that scales exponentially when the number
of variables of the problem increases. Hence, there is a practical hard limit
on the problem size when executed on a classical quantum simulator. To avoid
unbounded runtime of the problems on the simulator, a timeout limit of 24
hours is set in the executions. All molecules/runs that exceed this timeout
are excluded from the analysis of the results.

First study is to explore the active space selection methods as a function of
the size of the molecules, as presented in Figure 4. This shows the number of
qubits needed to represent the problem and the number of terms in the
Hamiltonian. The QULACS quantum circuit simulator runs on the Kay
cluster222ICHEC HPC cluster, Kay.
https://www.ichec.ie/about/infrastructure/kay typically on a single high-
memory node. Hartree-Fock, Frozen Core and NOONs are selected as active space
methods to to compare their error levels as a function of molecule size
(Figure 4 (a)). For the NOONs method, the upper threshold parameter of 0.002
is selected based on experiments yielding errors within the acceptable range.

Figure 4 shows how the error (compared against classically-computed CCSD as a
baseline), number of qubits and number of Hamiltonian terms vary with choice
of active space. In the top panel (a) compares the error levels to CCSD
results as a function of molecule size. HF and NOONs errors scale with number
of electrons while Frozen Core stay stable. The middle panel (b) shows the
number of qubits reduced by using the frozen core and NOONs approximations.
The bottom panel (c) show the impact on the number of terms in the
Hamiltonian, CH3F being the most impacted. The hatching refers to tapering on
top of the active space approximations.

Figure 4: Comparison of three active space methods, HF in red, Frozen Core in
orange and NOONs in green for (a) panel. Middle (b) and bottom (c) panels
correspond to number of qubits and number of Hamiltonian terms as function of
molecule size respectively. Blue bars refer to full set of orbitals included
in the calculation. The hatching refers to tapering on top of the active space
approximations.

For a group of selected molecules (only five to simplify the visualization),
Figure 5 illustrates the VQE error in upper panel (a) and the VQE runtime in
bottom panel (b) as a function of the number of active orbitals. For these
molecules we run experiments using three orbital reduction methods NOONs, Comb
and Hartree-Fock. In average, the VQE time scales up when the number of active
orbitals and the molecule size increase, as expected. The results for all
molecules tested is available on the QPFAS GitHub page333Qubit number scaling
https://github.com/ICHEC/QPFAS/blob/0.4.0-beta/py/examples/circuit_requirements/all_qubit_number_experiments.png.

Figure 5: Scaling up of the problem in terms of number of orbitals. Upper
panel (a) shows VQE error for O2, CO, CH4, F2 and CH3F as a function of the
number of active orbitals. VQE error is expressed in Hartrees (Ha) in left
vertical axis and in kcal/mol in right vertical axe. Three reduction orbital
methods where applied , NOONs, Comb and HF . Lower panel (b) is a plot of VQE
time (logarithm scale) to show the scaling up in terms of active orbitals for
the same group of molecules than upper panel

For the same group of selected molecules as above, Figure 6 shows how the
problem scales in terms of the number of parameters of the VQE ansatz for
increasing active orbitals and molecule size. The plot shows the ramp up in
calculation time by increasing the VQE ansatz parameters and the size of the
molecules, as expected. It also confirms that error levels can be
systematically lowered by increasing the complexity of the VQE ansatz, which
translates to increased computation times/resources. The results for all
molecules tested is available on the QPFAS GitHub page. 444Number of
parameters scaling
https://github.com/ICHEC/QPFAS/blob/0.4.0-beta/py/examples/circuit_requirements/all_parameter_number_experiments.png.

Figure 6: Scaling up of the problem in terms of  VQE ansatz parameters. Upper
panel (a) shows VQE error for O2, CO, CH4, F2 and CH3F as a function of VQE
ansatz parameters. VQE error is expressed in Hartree(Ha) in left vertical axe
and in kcal/mol in right vertical axe. Three reduction orbital methods where
applied NOONs, Comb and HF. Lower panel (b) shows VQE time (logarithm scale)
in terms of active orbitals for the same group of molecules than upper panel.
The error levels can be systematically lowered by increasing the complexity of
the VQE ansatz, which translates to increased computation times/resources.

##### Results for parameter optimisation

Parameter optimization is a set of tests using H2 and HF molecules to find
best parameters of the QPFAS workflow for larger simulations such as CH3F
molecule. The studied parameters include a list of optimizers (refer to the
figure in Github for details), the qubit transformations Jordan-Wigner and
Bravyi-Kitaev, a the range of five ansatze, and different sample sizes. We
concluded that hardware ansatze are less performant than chemistry inspired
ansatze. In addition, the results show that the use of one transformation or
the other does not privilege any optimizer. However, Powell is the optimizer
that presents the smallest error for kUpCCGSD, tapering and UCCSD in both
transformations. We also investigate the number of samples need to
successfully run the VQE algorithm. This study allow us to identify the number
of samples against final VQE energy for the H2 molecule. We run the workflow a
total of five times to measure the variance of the results. It can be inferred
from the results (refer to the figure in Github for details) that more than
10,000 samples must be used to achieve an error less than $0.001$ Ha. The
results also indicate that the BFGS optimizer performs better than COBYLA in
small sample sizes.

##### Results for bond stretching

Here we present the results on how well the VQE can describe the energy
profile of a molecule when bonds are stretched and broken. The previous
experiments helped to select a set of optimal parameters to run the bond-
breaking experiment. These experiment are done for F2 and CH3F molecules as
they can give us more information about the breakdown in PFAS substances
regarding C-F and F-F bonds. Standard VQE with UCCSD ansatz is used for F2,
and ADAPT-VQE with UCCSD is for CH3F. ADAPT-VQE is the most performant ansatz
at the time this computation is run. To map spin operators onto fermionic
operators we use Jordan-Wigner transformation for both molecules.

Figure 7 shows the dissociation curves of the F2 molecule as a function of
bond distance. The three curves of energy follow the same trend for lower bond
distances, however HF shows higher energy values when the bond starts to be
dissociated. FCI produces the same results as HF. The right vertical axe of
the figure shows the VQE error, which reaches its maximum value when the
molecule is in the dissociation stage. The error is quite small during all the
calculation steps.

Figure 7: Energy as a function of bond stretching for F2 molecule. Three
methods are calculated, blue dotted line corresponds to FCI, dashed black line
to HF and VQE standard UCCSD with Jordan-Wigner transformation in orange dot
points and full line. FCI and HF are superposed. Left vertical axe corresponds
to energy values expressed in Hartree(Ha), right vertical axe corresponds to
error of VQE calculation expressed in Ha units, the values of computed VQE
error are plotted with dashed blue line and dot points. The bond distance in
horizontal axe is in Angstrom(A).

Figure 8 shows the bond dissociation curve for CH3F by stretching the C-F
bond. The energy from HF is above the VQE result during the entire
dissociation process, which becomes much more significant when the C-F bond is
completely broken. However, FCI energy remains below the energy values
obtained with ADAPT-VQE. For CH3F we observe that the VQE error is maximum at
the ground state energy unlike what is observed for F2. Even if the error is
quite small during all the calculations in both molecules, we can observe that
the error is two orders of magnitude smaller for F2. The results show that VQE
can capture correlations and provide accuracy beyond HF, especially for a
bond-breaking experiment. Considering that FCI is not an scalable approach by
any means, this proves VQE as a valid alternative for post-HF quantum
chemistry calculation. This is specially important due to the suitability of
VQE for quantum computing technology, as explored in the next subsection.

Figure 8: CH3F using ADAPT-VQE, JW transformer, UCCSD and bond stretching .
Energy as a function of bond stretching for CH3F molecule. Three methods are
calculated, blue dotted line corresponds to FCI, dashed black line to HF and
ADAPT-VQE UCCSD with Jordan-Wigner transformation in orange dot points and
full line. Left vertical axe corresponds to energy values expressed in Hartree
(Ha), right vertical axe corresponds to error of VQE calculation expressed in
Ha units, the values of computed VQE error are plotted with dashed blue line
and dot points. The bond distance in horizontal axe is in Angstrom(A).

###  III-B Experiments using IonQ’s Aria trapped-ion quantum computer

A key objective of this work is to run experiments on currently available
quantum processing unit (QPU) devices, while looking ahead to advances in
quantum hardware that will enable the simulation of larger PFAS-like
molecules. The IonQ Aria quantum computer is integrated in the QPFAS workflow
as the quantum hardware where the end-to-end process can be executed. The
integration allows us to run shallower circuits with higher fidelity, as a
result of the accuracy and all-to-all connectivity of the trapped-ion
architecture. Hence, the results are less affected by the noise in quantum
computing devices. Given the capabilities of the current quantum hardware, we
tighten our focus to evaluating the scalability of the QPFAS workflow running
the C-F bond dissociation on a small PFAS-mimic molecule (CH3F). For this
experiment, we do not consider the convergence behaviour of the pUCCD ansatz
and instead obtain the classically-optimized parameters for the VQE algorithm.
This is a common approach on VQE studies both on trapped-ion[41, 24] and
superconducting devices [39]. VQE’s flexibility in circuit design allows fine
tuning based on quantum gate fidelity, number of qubits, and desired accuracy,
making it well-suited for the NISQ era. However, running shallower circuits in
VQE comes with its costs: predicted energy remains approximate, and many more
measurements are required.

At the time the simulations were performed (late 2022) IonQ had two publicly
available quantum computers to run the experiments proposed in this work; an
11-qubit trapped-ion QPU (Harmony) and a 25-qubit trapped-ion QPU named Aria.
As Aria’s gate fidelities are significantly higher than the previous
generation (Harmony) and it has sufficient qubit counts to simulate PFAS-like
molecules, we decided to perform the experiments on the Aria QPU in order to
achieve one of the largest quantum chemistry simulations on real quantum
hardware at the time.

As outlined in previous sections, the C-F bond is one of the strongest bonds
in chemistry. Thus, we evaluate PFAS-like molecules to break the C-F bond of
them. The CH3F molecule is chosen as the largest molecule for the experiments.
Hence, we carry out the next two experiments on the quantum computer Aria to
meet the objectives of the study:

  1. 1.

CH3F equilibrium geometry: we choose the equilibrium geometry of the CH3F
molecule to calculate its ground state. Since a sull VQE algorithm would
require a lot of computational time being executed in a real quantum hardware,
we perform most of the VQE iterations using IonQ’s cloud statevector
simulator[42] (to obtain the values of the optimized parameters of the
ansatz). Then, we execute the final energy evaluation using the IonQ Aria
quantum computer.

  2. 2.

CH3F stretched geometry: we set another experiment around the C-F bond
breaking of the CH3F molecule. For this experiment, a geometry of the molecule
is defined further away from the equilibrium geometry, trying to simulate how
the stretch would be. Similar to the above case, we execute almost all the
iterations of the VQE in the IonQ statevector simulator, then the final energy
evaluation is done using IonQ Aria.

These three experiments are performed under the same parameterization of the
QPFAS workflow, which is:

  * •

Basis set: STO-3G basis set is used to create the QPFAS molecule.

  * •

Active space: the frozen core approximation is used to create the QPFAS
molecule. In case of the CH3F molecule, it reduces the number of qubits by 2,
leaving the number of active orbitals on 11.

  * •

Transformation: Jordan-Wigner transformation is used to create the qubit
Hamiltonian.

  * •

Ansatz: the pUCCD ansatz is used as it is an efficient shallow ansatz, and has
shown high-performance in recent quantum experiments both at Google[39] and
IonQ[24]. In the pUCCD ansatz, the CH3F molecule has 11 qubits and 28
variational parameters after freezing the core electrons, with no more than
253 single qubit gates and 56 entangling gates.

####  III-B1 Results from IonQ quantum computer and error mitigation

TABLE I: CH3F energy of equilibrium geometry | threshold + symmetry post-selection |
---|---|---
Raw energy (Hartree) | 0.0% | 0.5% | 1.0% | 1.5% | 2.0% | exact
-135.74352 | -136.54827 | -136.76696 | -136.85525 | -137.14006 | -137.19404 | -137.19512
TABLE II: CH3F energy of stretched geometry | threshold + symmetry post-selection |
---|---|---
Raw energy (Hartree) | 0.0% | 0.5% | 1.0% | 1.5% | 2.0% | exact
-136.11887 | -136.68776 | -136.84109 | -136.85204 | -137.06390 | -137.06390 | -137.11248

Working in the noisy intermediate-scale quantum (NISQ) era, the results are
subject to errors. For VQE, where the expectation value of an electronic
Hamiltonian is to be estimated to high degrees of precision, noise can
significantly impact the computed results [30, 15, 16, 17, 18, 34, 20, 21, 22,
23, 24]. Although the purpose of this work is not to showcase error mitigation
protocols, it can be seen from the raw data that the stretched and equilibrium
geometries are in qualitative disagreement: the stretched geometry is _lower_
in energy than the equilibrium geometry, which is incorrect both theoretically
and physically. To mitigate this discrepancy and also to show a simple path to
recovering near-ideal results, we exploit two symmetries of the chemical
problem to mitigate the impact of noise on the results.

The first method we use is symmetry post-selection for the $\langle
Z\rangle$-basis measurements [43]. Given the Jordan-Wigner encoding and the
hard-core bosonic approximation inherent in pUCCD, the Hamming weight of all
symmetry-allowed bitstrings (output states) corresponds to the number of
electron pairs. Any output bitstring that deviates from the ideal Hamming
weight can be discarded, as it can only arise in error due to noise. Once the
erroneous bitstrings are discarded, the symmetry-corrected histogram may be
$L_{1}$-renormalized. Following this, we can see in Tables I and II that from
symmetry post-selection alone, a reduction in error on the order of 250-500
mHartree can be obtained. Despite this massive improvement in accuracy,
however, symmetry post-selection is not sufficient in this case to fix the
qualitative discrepancy between the energy ordering of the two CH3F
geometries.

The second error-mitigation technique is inspired by the fact that by setting
all variational parameters to zero, the Hartree-Fock energy—which is
efficiently known classically—can be computed for a given (noisy) circuit.
Moreover, in the $\langle Z\rangle$-basis, the ideal output corresponds to one
and only one bit-string, corresponding to the state preparation circuit.
Hence, by learning the ideal output for the quantum circuit, the empirical
impact of noise can be estimated as the proportion of the output that deviate
from the ideal [44]. In analogy with other signal processing techniques, such
as amplitude-limiting, clipping, and threshold-based peak detection, we
discard small states that cannot be reliably distinguished from noise. A
similar approach is described in Ref [45]. Therefore, applying different
truncation thresholds to the Hartree-Fock states (Fig. 10) we see that
discarding states with probability lower than 1.5% yields a good balance
between removing noise and (potentially) discarding signal. Since the circuits
and hardware are the same (differing only in the value of the variational
parameters) estimating the noise floor from the Hartree-Fock outputs can be
applied to the optimized circuit, where the output is not in general known _a
priori_. The impact of different noise-floor thresholds on the $\langle
Z\rangle$-basis can be observed in Fig. 9 as well as in Tables I and II. It
can be seen that the empirically-determined threshold of discarding states
below 1.5% is sufficient to restore the qualitative correctness of the
computed results on quantum hardware, as well as provide near-exact
quantitative accuracy. This method of truncation applies only to the $\langle
Z\rangle$-basis measurements, because the ideal output is a single bitstring.
The $\langle X\rangle$\- and $\langle Y\rangle$-basis measurements are
uniformly distributed in the Hartree-Fock case, preventing this type of noise
mitigation.

To further improve the results, more advanced noise mitigation techniques
could be used, as this is an active and ongoing area of research [43, 45, 46,
47, 20, 39, 48]. However, by using simple noise mitigation methods, we show
here that it is possible to deliver quantum value on today’s NISQ devices.

Figure 9: Results of thresholding-based noise mitigation and symmetry post-
selection on the $\langle Z\rangle$-basis measurements for (a) equilibrium and
(b) stretched CH3F geometries with optimized pUCCD parameters. Histogram
bitstrings with probability less than the threshold value are discarded, along
with bitstrings that are known to be wrong by symmetry (e.g. not particle
conserving). The histograms are then $L_{1}$-renormalized. From the known
Hartree-Fock results in Fig. 10, we can determine that the empirically optimal
noise-floor threshold is at probabilities $\leq 1.5\%$. Error bars correspond
to 2$\sigma$, where $\sigma$ is the standard error of the estimated
expectation values.  Figure 10: Results of thresholding-based noise mitigation
and symmetry post-selection on the $\langle Z\rangle$-basis measurements for
(a) equilibrium and (b) stretched CH3F geometries with Hartree-Fock
parameters. Histogram bitstrings with probability less than the threshold
value are discarded, along with bitstrings that are known to be wrong by
symmetry (e.g. not particle conserving). The histograms are then 1-norm
renormalized. Because the Hartree-Fock energy is efficiently and exactly
known, and moreover the circuit topology is identical to the optimized pUCCD
case, the Hartree-Fock hardware experiments allow us to estimate the noise
floor, which can then be applied to identical circuits where the energy is not
_a priori_ known. Error bars correspond to 2$\sigma$, where $\sigma$ is the
standard error of the estimated expectation values.

##  IV Conclusions

Remediation of global PFAS pollution is a an active research topic that aims
at finding new mechanisms to destroy these substances. A major pathway for
PFAS destruction is breaking carbon-fluorine bonds. Computational simulating
of this process can speed-up the research towards finding optimal processes
and catalysts. Solving the electronic Schrodinger problem is at the core of
computational chemistry, which has always been a technical challenge due to
its exponential complexity. The original idea of quantum computing came from
this need for simulating nature at quantum level. In this article, we proposed
the QPFAS workflow, a hybrid quantum-classical workflow for end-to-end
chemistry simulation based on VQE algorithm. The QPFAS workflow allows
flexible definition of many experiments, including bond-breaking calculations,
and orchestrates the execution of sub-tasks on both HPC-based simulators and
quantum computers.

We completed simulation experiments on a group of molecules to understand how
the VQE algorithm scales in terms of molecule size, number of active orbitals,
and optimization parameters. Two molecules that approach PFAS substances were
identified, F2 and CH3F, with which breaking-bond simulation experiments have
been carried out on a quantum computing simulator, reaching the limits of
classical computing at ICHEC’s HPC system. The results on quantum simulators
were consistent with the expected theoretical analysis, and helped to better
understand the scalability and accuracy of the VQE algorithm. They provided
important understanding on how to adjust the parameters and port the
calculations to a real quantum computer.

Moreover, we studied the bond breaking of CH3F molecule with IonQ’s trapped-
ion quantum computer. To address the limitations of noisy quantum computing
technology, we introduced two error mitigation techniques: (a) the symmetry
post-selection for the ⟨Z⟩-basis measurements, and (b) correcting the error
bias by calculating the classically-known Hartree-Fock energy in the noisy
device. This error mitigation approach returned near-quantitative results,
with the CH3F equilibrium energy at milli-Hartree precision with the expected
result. The accuracy of the results obtained for the CH3F molecule is
unprecedented on a quantum computer for any bond-breaking experiment of a
comparable size.

The creation of the QPFAS workflow allowed the hybrid study of the bond
dissociation problem with great flexibility in parameterization and
orchestration of large number of parallel tasks. This enabled us to optimize
the quantum computer resources by selecting the most efficient parameters,
leading to one of the largest-ever molecules studied on computer computer. The
workflow can be used to conduct a wide range of computational chemistry
experiments, and further develop the VQE algorithm for quantum computers.

## Acknowledgments

This project was supported in part by Enterprise Ireland Innovation
Partnership Programme (IP/2019/0835) and EuroCC in Ireland (grant agreement
no. 951732).

## References

  * Agency [2016] E. U. S. P. Agency, “Lifetime health advisories and health effects support documents for perfluorooctanoic acid and perfluorooctane sulfonate.” _Fed. Reg._ , pp. 33 250–33 251, 2016. [Online]. Available: https://www.federalregister.gov/documents/2016/05/25/2016-12361/lifetime-health-advisories-and-health-effects-support-documents-for-perfluorooctanoic-acid-and
  * Grandjean et al. [2017] P. Grandjean, C. Heilmann, P. Weihe _et al._ , “Estimated exposures to perfluorinated compounds in infancy predict attenuated vaccine antibody concentrations at age 5-years,” _Journal of Immunotoxicology_ , vol. 14, no. 1, pp. 188–195, 2017, pMID: 28805477. [Online]. Available: https://doi.org/10.1080/1547691X.2017.1360968
  * Looker et al. [2013] C. Looker, M. I. Luster, A. M. Calafat _et al._ , “Influenza Vaccine Response in Adults Exposed to Perfluorooctanoate and Perfluorooctanesulfonate,” _Toxicological Sciences_ , vol. 138, no. 1, pp. 76–88, 11 2013. [Online]. Available: https://doi.org/10.1093/toxsci/kft269
  * Dombrowski et al. [2018] P. M. Dombrowski, P. Kakarla, W. Caldicott _et al._ , “Technology review and evaluation of different chemical oxidation conditions on treatability of pfas,” _Remediation Journal_ , vol. 28, no. 2, pp. 135–150, 2018. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/rem.21555
  * Rahman et al. [2014] M. F. Rahman, S. Peldszus, and W. B. Anderson, “Behaviour and fate of perfluoroalkyl and polyfluoroalkyl substances (pfass) in drinking water treatment: A review,” _Water Research_ , vol. 50, pp. 318–340, 2014. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0043135413008518
  * Kumarasamy et al. [2020] E. Kumarasamy, I. M. Manning, L. B. Collins _et al._ , “Ionic fluorogels for remediation of per- and polyfluorinated alkyl substances from water,” _ACS Central Science_ , vol. 6, no. 4, pp. 487–492, 2020. [Online]. Available: https://doi.org/10.1021/acscentsci.9b01224
  * Xiao et al. [2017] X. Xiao, B. A. Ulrich, B. Chen, and C. P. Higgins, “Sorption of poly- and perfluoroalkyl substances (pfass) relevant to aqueous film-forming foam (afff)-impacted groundwater by biochars and activated carbon,” _Environmental Science & Technology_, vol. 51, no. 11, pp. 6342–6351, 2017, pMID: 28582977. [Online]. Available: https://doi.org/10.1021/acs.est.7b00970
  * Bentel et al. [2019] M. J. Bentel, Y. Yu, L. Xu _et al._ , “Defluorination of per- and polyfluoroalkyl substances (pfass) with hydrated electrons: Structural dependence and implications to pfas remediation and management,” _Environmental Science & Technology_, vol. 53, no. 7, pp. 3718–3728, 2019, pMID: 30874441. [Online]. Available: https://doi.org/10.1021/acs.est.8b06648
  * Bruton and Sedlak [2017] T. A. Bruton and D. L. Sedlak, “Treatment of aqueous film-forming foam by heat-activated persulfate under conditions representative of in situ chemical oxidation,” _Environmental Science & Technology_, vol. 51, no. 23, pp. 13 878–13 885, 2017, pMID: 29164864.
  * Park et al. [2016] S. Park, L. S. Lee, V. F. Medina _et al._ , “Heat-activated persulfate oxidation of pfoa, 6:2 fluorotelomer sulfonate, and pfos under conditions suitable for in-situ groundwater remediation,” _Chemosphere_ , vol. 145, pp. 376–383, 2016. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0045653515304288
  * Saleem et al. [2020] M. Saleem, O. Biondo, G. Sretenovia _et al._ , “Comparative performance assessment of plasma reactors for the treatment of pfoa; reactor design, kinetics, mineralization and energy yield,” _Chemical Engineering Journal_ , vol. 382, p. 123031, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S138589471932443X
  * Singh et al. [2020] R. K. Singh, N. Multari, C. Nau-Hix _et al._ , “Removal of poly- and per-fluorinated compounds from ion exchange regenerant still bottom samples in a plasma reactor,” _Environmental Science & Technology_, vol. 54, no. 21, pp. 13 973–13 980, 2020, pMID: 33103418.
  * Liu et al. [2021] Z. Liu, M. J. Bentel, Y. Yu _et al._ , “Near-quantitative defluorination of perfluorinated and fluorotelomer carboxylates and sulfonates with integrated oxidation and reduction,” _Environmental Science & Technology_, vol. 55, no. 10, pp. 7052–7062, 2021, pMID: 33950686.
  * Peruzzo et al. [2014a] A. Peruzzo, J. McClean, P. Shadbolt _et al._ , “A variational eigenvalue solver on a photonic quantum processor,” _Nat. Commun._ , vol. 5, p. 4213, Jul. 2014.
  * O’Malley et al. [2016] P. J. J. O’Malley, R. Babbush, I. D. Kivlichan _et al._ , “Scalable quantum simulation of molecular energies,” _Phys. Rev. X_ , vol. 6, no. 3, p. 031007, Jul. 2016.
  * Colless et al. [2018] J. I. Colless, V. V. Ramasesh, D. Dahlen _et al._ , “Computation of molecular spectra on a quantum processor with an Error-Resilient algorithm,” _Phys. Rev. X_ , vol. 8, no. 1, p. 011021, Feb. 2018.
  * McCaskey et al. [2019] A. J. McCaskey, Z. P. Parks, J. Jakowski _et al._ , “Quantum chemistry as a benchmark for near-term quantum computers,” _npj Quantum Information_ , vol. 5, no. 1, pp. 1–8, Nov. 2019.
  * Nam et al. [2020] Y. Nam, J.-S. Chen, N. C. Pisenti _et al._ , “Ground-state energy estimation of the water molecule on a trapped-ion quantum computer,” _npj Quantum Information_ , vol. 6, no. 1, pp. 1–6, Apr. 2020.
  * Kandala et al. [2017a] A. Kandala, A. Mezzacapo, K. Temme _et al._ , “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” _Nature_ , vol. 549, no. 7671, pp. 242–246, Sep. 2017.
  * Kandala et al. [2019] A. Kandala, K. Temme, A. D. Córcoles _et al._ , “Error mitigation extends the computational reach of a noisy quantum processor,” _Nature_ , vol. 567, no. 7749, pp. 491–495, Mar. 2019.
  * Gao et al. [2019] Q. Gao, H. Nakamura, T. P. Gujarati _et al._ , “Computational investigations of the lithium superoxide dimer rearrangement on noisy quantum devices,” Jun. 2019.
  * Rice et al. [2021] J. E. Rice, T. P. Gujarati, M. Motta _et al._ , “Quantum computation of dominant products in lithium-sulfur batteries,” _J. Chem. Phys._ , vol. 154, no. 13, p. 134115, Apr. 2021.
  * Gao et al. [2021] Q. Gao, G. O. Jones, M. Motta _et al._ , “Applications of quantum computing for investigations of electronic transitions in phenylsulfonyl-carbazole TADF emitters,” _npj Computational Materials_ , vol. 7, no. 1, pp. 1–9, May 2021.
  * Zhao et al. [2023] L. Zhao, J. Goings, K. Shin _et al._ , “Orbital-optimized pair-correlated electron simulations on trapped-ion quantum computers,” _npj Quantum Information_ , vol. 9, no. 1, pp. 1–9, Jun. 2023.
  * Goings et al. [2023] J. Goings, L. Zhao, J. Jakowski _et al._ , “Molecular symmetry in VQE: A dual approach for Trapped-Ion simulations of benzene,” Aug. 2023.
  * [26] “Qulacs simulator,” http://docs.qulacs.org/en/latest/, accessed: 2023-10-1.
  * [27] “Ionq,” https://ionq.com/, accessed: 2023-10-31.
  * Kottmann et al. [2021] J. S. Kottmann, S. Alperin-Lea, T. Tamayo-Mendoza _et al._ , “Tequila: A platform for rapid development of quantum algorithms,” _Quantum Science and Technology_ , vol. 6, no. 2, p. 024009, 2021.
  * Sun et al. [2018] Q. Sun, T. C. Berkelbach, N. S. Blunt _et al._ , “Pyscf: the python-based simulations of chemistry framework,” _Wiley Interdisciplinary Reviews: Computational Molecular Science_ , vol. 8, no. 1, p. e1340, 2018.
  * Peruzzo et al. [2014b] A. Peruzzo, J. McClean, P. Shadbolt _et al._ , “A variational eigenvalue solver on a photonic quantum processor,” _Nature communications_ , vol. 5, no. 1, pp. 1–7, 2014.
  * Goings et al. [2022] J. J. Goings, A. White, J. Lee _et al._ , “Reliably assessing the electronic structure of cytochrome P450 on today’s classical computers and tomorrow’s quantum computers,” _Proc. Natl. Acad. Sci. U. S. A._ , vol. 119, no. 38, p. e2203533119, Sep. 2022.
  * 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_ , vol. 10, no. 1, pp. 1–9, 2019.
  * Bravyi et al. [2017] S. Bravyi, J. M. Gambetta, A. Mezzacapo, and K. Temme, “Tapering off qubits to simulate fermionic hamiltonians,” _arXiv preprint arXiv:1701.08213_ , 2017\.
  * Kandala et al. [2017b] A. Kandala, A. Mezzacapo, K. Temme _et al._ , “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” _Nature_ , vol. 549, no. 7671, pp. 242–246, 2017.
  * McClean et al. [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy _et al._ , “Barren plateaus in quantum neural network training landscapes,” _Nat. Commun._ , vol. 9, no. 1, p. 4812, Nov. 2018.
  * Gonthier et al. [2020] J. F. Gonthier, M. D. Radin, C. Buda _et al._ , “Measurements as a roadblock to near-term practical quantum advantage in chemistry: resource analysis,” Dec. 2020.
  * Elfving et al. [2020a] V. E. Elfving, B. W. Broer, M. Webber _et al._ , “How will quantum computers provide an industrially relevant computational advantage in quantum chemistry?” Sep. 2020.
  * Elfving et al. [2020b] V. E. Elfving, M. Millaruelo, J. A. Gámez, and C. Gogolin, “Simulating quantum chemistry in the seniority-zero space on qubit-based quantum computers,” Jan. 2020.
  * O’Brien et al. [2022] T. E. O’Brien, G. Anselmetti, F. Gkritsis _et al._ , “Purification-based quantum error mitigation of pair-correlated electron simulations,” Oct. 2022\.
  * Barkoutsos et al. [2018] P. K. Barkoutsos, J. F. Gonthier, I. Sokolov _et al._ , “Quantum algorithms for electronic structure calculations: Particle-hole hamiltonian and optimized wave-function expansions,” _Physical Review A_ , vol. 98, no. 2, p. 022322, 2018.
  * Khan et al. [2022] I. T. Khan, M. Tudorovskaya, J. J. M. Kirsopp _et al._ , “Chemically aware unitary coupled cluster with ab initio calculations on system model h1: A refrigerant chemicals application,” Oct. 2022.
  * [42] “IonQ quantum cloud,” https://ionq.com/quantum-cloud, accessed: 2023-5-22.
  * Bonet-Monroig et al. [2018] X. Bonet-Monroig, R. Sagastizabal, M. Singh, and T. E. O’Brien, “Low-cost error mitigation by symmetry verification,” _Phys. Rev. A_ , vol. 98, no. 6, p. 062339, Dec. 2018.
  * Urbanek et al. [2021] M. Urbanek, B. Nachman, V. R. Pascuzzi _et al._ , “Mitigating depolarizing noise on quantum computers with Noise-Estimation circuits,” _Phys. Rev. Lett._ , vol. 127, no. 27, p. 270502, Dec. 2021.
  * Montanaro and Stanisic [2021] A. Montanaro and S. Stanisic, “Error mitigation by training with fermionic linear optics,” Feb. 2021.
  * Bennewitz et al. [2022] E. R. Bennewitz, F. Hopfmueller, B. Kulchytskyy _et al._ , “Neural error mitigation of Near-Term quantum simulations,” _Nature Machine Intelligence_ , vol. 4, no. 7, pp. 618–624, Jul. 2022.
  * Quek et al. [2022] Y. Quek, D. S. França, S. Khatri _et al._ , “Exponentially tighter bounds on limitations of quantum error mitigation,” Oct. 2022.
  * Weber et al. [2021] T. Weber, M. Riebisch, K. Borras _et al._ , “Modelling for quantum error mitigation,” Apr. 2021.


# Evidence of a dibaryon spectrum in coherent $\pi^{0}\pi^{0}d$
photoproduction at forward deuteron angles

T.C. Jude<EMAIL_ADDRESS>S. Alef111No longer employed in academia R.
Beck  A. Braghieri  P.L. Cole  D. Elsner  R. Di Salvo  A. Fantini  O.
Freyermuth  F. Frommberger  F. Ghio  S. Goertz  A. Gridnev  K. Kohl  N.
Kozlenko  A. Lapik  P. Levi Sandri  V. Lisin  G. Mandaglio  D.
Moricciani222Deceased V. Nedorezov333Deceased D. Novinskiy  P. Pedroni  A.
Polonskiy  B.-E. Reitz 444No longer employed in academia M. Romaniuk  G.
Scheluchin555No longer employed in academia H. Schmieden  A. Stuglev  V.
Sumachev666Deceased V. Tarakanov  Rheinische Friedrich-Wilhelms-Universität
Bonn, Physikalisches Institut, Nußallee 12, 53115 Bonn, Germany  Rheinische
Friedrich-Wilhelms-Universität Bonn, Helmholtz-Institut für Strahlen- und
Kernphysik, Nußallee 14-16, 53115 Bonn, Germany  Petersburg Nuclear Physics
Institute NRC “Kurchatov Institute”, Gatchina, Leningrad District, 188300,
Russia  INFN sezione di Pavia, Via Agostino Bassi, 6 - 27100 Pavia, Italy
Lamar University, Department of Physics, Beaumont, Texas, 77710, USA  INFN
Roma “Tor Vergata”, Via della Ricerca Scientifica 1, 00133, Rome, Italy
Università di Roma “Tor Vergata”, Dipartimento di Fisica, Via della Ricerca
Scientifica 1, 00133, Rome, Italy  INFN sezione di Roma La Sapienza, P.le Aldo
Moro 2, 00185, Rome, Italy  Istituto Superiore di Sanità, Viale Regina Elena
299, 00161, Rome, Italy  Russian Academy of Sciences Institute for Nuclear
Research, Prospekt 60-letiya Oktyabrya 7a, 117312, Moscow, Russia  INFN -
Laboratori Nazionali di Frascati, Via E. Fermi 54, 00044, Frascati, Italy
INFN sezione Catania, 95129, Catania, Italy  Università degli Studi di
Messina, Dipartimento MIFT, Via F. S. D’Alcontres 31, 98166, Messina, Italy
Institute for Nuclear Research of NASU, 03028, Kyiv, Ukraine

###### Abstract

The coherent reaction, $\gamma d\rightarrow\pi^{0}\pi^{0}d$ was studied with
the BGOOD experiment at ELSA from threshold to a centre-of-mass energy of 2850
MeV. A full kinematic reconstruction was made, with final state deuterons
identified in the forward spectrometer and $\pi^{0}$ decays in the central BGO
Rugby Ball. The strength of the differential cross section exceeds what can be
described by models of coherent photoproduction and instead supports the three
isoscalar dibaryon candidates reported by the ELPH collaboration at 2.38, 2.47
and 2.63 GeV/c2. A low mass enhancement in the $\pi^{0}\pi^{0}$ invariant mass
is also observed at the $d^{*}(2380)$ centre-of-mass energy which is
consistent with the ABC effect. At higher centre-of-mass energies, a narrow
peak in the $\pi^{0}d$ invariant mass at 2114 MeV/c2 with a width of 20 MeV/c2
supports a sequential two-dibaryon decay mechanism.

###### keywords:

BGOOD , dibaryon , coherent photoproduction

###### PACS:

13.60.Le,25.20.-x

††journal: Physics Letters B

##  1 Introduction

Our understanding of hadron structure and relevant degrees of freedom is a
crucial test for QCD in the low energy non-perturbative regime. Until
recently, the deuteron was considered the only bound dibaryon system, however
as early as 1964, Dyson and Xuong [1] predicted a sextet of dibaryon states
from SU(6) symmetry for baryons, denoted $\mathcal{D}_{IJ}$ for isospin $I$
and spin $J$. With the addition of the ground state deuteron,
$\mathcal{D}_{01}$, and the virtual unbound state, $\mathcal{D}_{10}$,
observed as structure in $NN$ scattering, four non-strange states,
$\mathcal{D}_{12}$ and $\mathcal{D}_{21}$ ($N\Delta$), and $\mathcal{D}_{03}$
and $\mathcal{D}_{30}$ ($\Delta\Delta$) were predicted, the masses of which
were determined from the deuteron mass and nucleon scattering data. These
calculations led to a plethora of searches for dibaryon states throughout the
1960s and 70s, with mixed interpretations (for a recent review see ref. [2]
and earlier reviews in refs. [3, 4]). In addition to genuine dibaryon states,
One Pion Exchange models (OPE) between nucleons, for example, were also able
to describe peaks observed in $\pi\pi d$ and $\pi d$ systems [5, 6]. Limited
experimental statistical precision prevented firm conclusions on this however.

Structures in the mass range of the suspected $d^{*}(2380)$ hexaquark, first
identified in the fusion reaction $pn\rightarrow d\pi^{0}\pi^{0}$ [7, 8] have
sparked renewed interest in dibaryon searches in the non-strange sector (see
ref. [9] for a review). A strong indication of the $d^{*}(2380)$, with
$IJ^{P}=03^{+}$, has been observed in a multitude of final states and
observables [10, 11, 12, 13, 14, 15, 16], and it has been considered the
$\mathcal{D}_{03}$ in a $\Delta\Delta$ configuration. This agrees with the
predictions of Goldman [17], a three-body $\pi NN$ calculation of Gal and
Garcilazo [18], and is further supported by partial wave analyses [19, 20, 21,
22, 23]. Recent models benefiting from high precision experimental data and
theoretical developments have also described the $d^{*}(2380)$ structure via
OPE mechanisms and triangle singularities [24, 25], however it has been
debated whether such models describe the data consistently well over the broad
range of kinematics and reaction channels that the $d^{*}(2380)$ may have been
observed in [26].

Numerous calculations have determined a sizeable colour confined hexaquark
component of the $d^{*}(2380)$ (see for example, ref. [27]) and if it exists,
its structure and properties may have important astrophysical implications
[28, 29]. Electromagnetic production of a dibaryon directly from the deuteron
ground state would impose constraints on its size and structure via transition
form factors. The ELPH collaboration measured the $\gamma
d\rightarrow\pi^{0}\pi^{0}d$ cross section [30] which was in agreement with
the coherent photoproduction model of Fix, Arenhövel and Egorov [31, 32]. The
lowest energy data point at approximately $E_{\gamma}=570$ MeV however
suggested a $d^{*}(2380)$ contribution, which was also supported by
preliminary data from the A2 collaboration [33]. An additional ELPH dataset at
higher energies [34] exhibited peaks at 2.47 GeV/c2 and 2.63 GeV/c2 which
could be described by the quasi-free excitation of one nucleon followed by a
coalescence of the nucleons to the deuteron. The angular distributions however
required two isoscalar dibaryons in the reaction mechanism. There was also an
indication of the 2.14 GeV/$c^{2}$ isovector $\mathcal{D}_{12}$ in the
$\pi^{0}d$ invariant mass spectrum.

The reaction $\gamma d\rightarrow\pi^{0}\pi^{0}d$ is an ideal channel to
search for dibaryons. The isoscalar final state is only sensitive to
intermediate isoscalar dibaryons, compared to $\gamma
d\rightarrow\pi^{+}\pi^{-}d$ which also has isovector coupling and background
contributions from the large $\gamma N\pi^{\pm}$ coupling in the Kroll-
Ruderman term (see for example ref. [31]). Isovector dibaryons can still be
identified in the $\pi^{0}d$ invariant mass via sequential dibaryon decays.
The suppressed cross section for conventional coherent processes, which
reduces quickly with momentum transfer to the deuteron, may also help to
identify structure originating from dibaryon formation.

Experimentally, a kinematic complete identification of the reaction with a
clean separation of deuterons and protons is essential to separate the small
number of events from coherent reactions compared to the large amount of non-
coherent quasi-free events from $\gamma p(n)\rightarrow\pi^{0}\pi^{0}p$. The
BGOOD experiment [35] at ELSA [36, 37], described in sec. 2 is ideal, as the
Forward Spectrometer cleanly separates charged particles from their mass
reconstruction and the BGO Rugby Ball identifies neutral meson decays.

This letter presents differential cross section data for $\gamma
p\rightarrow\pi^{0}\pi^{0}d$ at forward deuteron angles, which cannot be
described via a conventional coherent reaction as this is strongly suppressed
due to the large momentum transfer to the deuteron. The data is intended to
shed more light on the contrasting interpretations of earlier observed
structures and thus on the possible existence of genuine dibaryon resonances.

##  2 Experimental setup and analysis procedure

BGOOD [35] is comprised of two main parts: a central region, ideal for neutral
meson identification, and a forward spectrometer for charged particle
identification and momentum reconstruction. The BGO Rugby Ball is the main
detector over the central region, covering laboratory polar angles 25 to 155∘.
The detector is comprised of 480 BGO crystals for the reconstruction of photon
momenta via electromagnetic showers in the crystals. The separate time readout
per crystal enables a clean separation and identification of neutral meson
decays. Between the BGO Rugby Ball and the target are the Plastic
Scintillating Barrel for charged particle identification via $\Delta E-E$
techniques and the MWPC for charged particle tracking and vertex
reconstruction.

The forward spectrometer covers a laboratory polar angle 1-12∘. The tracking
detectors, MOMO and SciFi are used to track charged particles from the target.
Downstream of these is the Open Dipole Magnet operating at an integrated field
strength of 0.216 Tm. A series of eight double sided Drift Chambers track
charged particle trajectories after the curvature in the magnetic fields and
are used to determine particle momenta with a resolution of approximately 6
%777The resolution improves to 3 % if the Open Dipole Magnet is operating at
the maximum field strength of 0.432 Tm.. Three Time of Flight (ToF) Walls
downstream of the drift chambers determine particle $\beta$ and are used in
combination with the measured momentum for particle identification via mass
determination. Track reconstruction in the Forward Spectrometer is described
in ref. [35].

The small intermediate region between the central region and the Forward
Spectrometer is covered by SciRi, which consists of concentric rings of
plastic scintillators for charged particle detection.

The deuterium target data presented was taken over 26 days using an 11 cm long
target and an ELSA electron beam energy of 2.9 GeV. The electron beam was
incident upon a thin diamond radiator888A diamond radiator was used to produce
coherent, linearly polarised photon beam with a maximum polarisation at a beam
energy of 1.4 GeV, however the polarisation was not required for the presented
analysis. to produce an energy tagged bremsstrahlung photon beam which was
subsequently collimated. The photon beam energy, $E_{\gamma}$, was determined
per event by momentum analysing the post bremsstrahlung electrons in the
Photon Tagger. The integrated photon flux from $E_{\gamma}=450$ to 1200 MeV
was $7.5\times 10^{12}$. The hardware trigger used for the presented analysis
required an energy tagged incident photon and an energy deposition in the BGO
Rugby Ball of approximately 150 MeV (see ref. [39] for details).

Candidate events were selected where exactly four photons were identified in
the BGO Rugby Ball (via a veto with the Plastic Scintillating Barrel) and one
charged particle in the Forward Spectrometer, corresponding to the two
$\pi^{0}\rightarrow\gamma\gamma$ decays and a forward going deuteron. Events
were rejected if any additional charged particle was identified.

The invariant mass of each two photon system was required to be within 40 MeV
of the $\pi^{0}$ mass (corresponding to approximately 2.5$\sigma$). For a
given event, all combinations where two $\pi^{0}$ from the four photons could
be reconstructed were retained for further analysis.

Two mass determinations of the forward going particle were made per event. The
first was the mass reconstruction in the forward spectrometer via momentum and
$\beta$ measurements, referred to as the forward spectrometer mass. The second
was the missing mass recoiling from the $\pi^{0}\pi^{0}$ system, referred to
as the $\pi^{0}\pi^{0}$ missing mass. For the dominant quasi-free background
reaction, $\gamma p(n)\rightarrow\pi^{0}\pi^{0}p$, the forward spectrometer
mass is peaked at the proton mass, whereas the $\pi^{0}\pi^{0}$ missing mass
is higher than the deuteron mass and increases with $E_{\gamma}$ due to the
assignment of the target mass to that of the deuteron.

The blue line in Fig.1(a) shows the forward spectrometer mass after two
$\pi^{0}$ have been identified. There is a small peak at the $\pi^{+}$ and a
larger peak at the proton mass from quasi-free production. The peak
corresponding to deuterons is barely discernible in the high energy tail of
the proton candidates. Additional selection criteria were used to enhance the
signal from deuterons in order to visualise the mass distribution. The angle
between the missing momentum from the $\pi^{0}\pi^{0}$ system and the forward
going charged particle was required to be smaller than 7.5∘ and the
$\pi^{0}\pi^{0}$ missing mass was required to be lower than 1900 MeV/c2 to
remove background from quasi-free production off the proton. The red line in
Fig.1(a) (scaled to equal the number of protons in the blue spectrum) shows
the forward spectrometer mass after these selection criteria, where a peak at
the deuteron mass is clear. Fig.1(b) shows the same data (not scaled and with
finer binning), where a fit has been made using simulated $\gamma
d\rightarrow\pi^{0}\pi^{0}d$ events and real data of quasi-free production off
the proton, with good agreement to the data.

Figure 1: Forward spectrometer mass, with the requirement that 2$\pi^{0}$ are
identified in the BGO Rugby Ball and the forward going particle momentum is
between 800 to 950 MeV/c. (a) All events from the dataset are shown in blue.
The red line corresponds to events passing additional selection criteria
described in the text and scaled by a factor of 2000. (b) A fit to the events
(black circles) passing the additional selection criteria, including simulated
$\pi^{0}\pi^{0}d$ events (shaded green area) and quasi-free proton background
(blue line). The total fit is the red line.

The mass distribution in Fig.1(b) was used to determine a sensible deuteron
mass selection cut for the forward spectrometer mass, after which the
selection criterion on the $2\pi^{0}$ missing mass to enhance the deuteron
signal was removed. The forward spectrometer mass was subsequently required to
be between 1550 to 2500 MeV/c2, corresponding to $1.3\sigma$ below to
$2.6\sigma$ above the deuteron mass (approximating the distribution as
Gaussian). The asymmetric selection limited the contribution of background
from quasi-free protons.

The $\pi^{0}\pi^{0}$ missing mass after these selection criteria is shown in
Fig. 2. A peak at the deuteron mass of 1876 MeV/$c^{2}$ is evident, with
background from quasi-free production off the proton at higher mass. This
background contribution is described by changing the forward spectrometer mass
selection to between 600 to 1300 MeV/$c^{2}$ over the proton mass999An
alternative method to describe the quasi-free background was studied using
liquid hydrogen data and selecting forward particles with a forward
spectrometer mass over the deuteron mass region. The lower statistics and the
absence of the target Fermi motion however demonstrated that selecting the
quasi-free events from the liquid deuterium data was preferred.. A fit was
made to extract the yield of $\gamma d\rightarrow\pi^{0}\pi^{0}d$ events using
this quasi-free background and simulated $\pi^{0}\pi^{0}d$ events with a phase
space distribution. The two distributions, shown in fig. 2 as the shaded green
and blue regions are scaled to minimise the $\chi^{2}$ for each $E_{\gamma}$
interval.

Figure 2: $\pi^{0}\pi^{0}$ missing mass after the selection criteria described
in the text. Every other photon beam energy interval, $E_{\gamma}$ is shown up
to 1024 MeV, and labelled inset in units of MeV together with $W$ and the
reduced $\chi^{2}$ of the fit. The data are the black circles, and the fit to
the data is the thick red line. Fitted contributions from simulated
$\pi^{0}\pi^{0}d$ and background from quasi-free events off the proton are
shown as the shaded green and blue spectra respectively.

Fig. 3 shows the detection efficiency for the reconstruction of $\gamma
d\rightarrow\pi^{0}\pi^{0}d$ events assuming a phase space distribution and if
a sequential decay of two dibaryons is assumed (discussed later in sec. 3).
The differential cross section in sec. 3 used the weighted detection
efficiency of these two distributions, indicated by the purple stars.

Figure 3: $\gamma d\rightarrow\pi^{0}\pi^{0}d$ detection efficiency for
$\cos\theta_{\mathrm{CM}}^{d}$ $>0.8$ assuming a phase space distribution and
an intermediate isovector dibaryon (black circles and red squares
respectively) . The efficiency of detecting background channels despite
selection criteria removing most events are shown as the green and blue
triangles for $\gamma d\rightarrow\pi^{0}\pi^{0}\pi^{0}d$ and $\gamma
d\rightarrow\eta\pi^{0}d$ respectively.

Simulated coherent production of $3\pi^{0}d$ and $\eta\pi^{0}d$ events were
used to investigate possible background from these channels where the cross
sections are poorly known. The false efficiency of the extent that these
channels pass the selection criteria is also shown in fig. 3. With
approximately 25 % relative to the $\pi^{0}\pi^{0}d$ efficiency, it is clear
that they cannot be neglected. The $\pi^{0}\pi^{0}$ missing mass spectra were
therefore also fitted including the additional simulated background channels.
The fit could accommodate a significant fraction of $3\pi^{0}d$ background, an
example of which is shown in fig. 4(b). This additional background however was
predominantly under the quasi-free proton events and did not significantly
affect the extracted yield of the $\pi^{0}\pi^{0}d$ channel or the reduced
$\chi^{2}$ of the fit over most of the beam energy range. The systematic
uncertainty in the yield extraction arising from this background contribution
increased approximately exponentially with beam energy, and was estimated as 5
%, 12 % and 26 % for $E_{\gamma}$ = 900, 1000 and 1100 MeV respectively. The
systematic uncertainty arising from reactions off the target cell windows was
determined by an equivalent analysis of data using a hydrogen target. Shown in
Fig. 4(c), it is clear there is negligible signal with only background from
reactions off the proton. The uncertainty was therefore considered to be 1 %
or less.

Figure 4: An example of systematic uncertainties evident in $\pi^{0}\pi^{0}$
missing mass at $E_{\gamma}=885$ MeV. (a) A “standard” fit using the simulated
signal and quasi-free proton background, as is shown in fig. 2. (b) The
inclusion of additional background from simulated $\gamma d\rightarrow
3\pi^{0}d$ events (filled magenta). (c) Using a liquid hydrogen target.

Systematic uncertainties from the modelling of the detector efficiencies (3.6
%), forward track finding (1.0 %), timing cuts (2.0 %) and beam spot alignment
(4.0 %) were measured previously [39]. Systematic uncertainties from the
$\pi^{0}\rightarrow\gamma\gamma$ identification (3.5 % for both $\pi^{0}$) and
the angle selection cut between the missing and measured deuteron momentum
(1.0 %) were determined by an equivalent analysis of the reaction $\gamma
p\rightarrow\pi^{0}\pi^{0}p$ with a hydrogen target where the beam and trigger
conditions were identical to the deuterium target dataset. The measured
differential cross section (not shown) gave a good agreement to previous
results from the A2 Collaboration [40].

##  3 Results and interpretation

###  3.1 The $\gamma d\rightarrow\pi^{0}\pi^{0}d$ differential cross section
versus $W$

The differential cross section for $\cos\theta_{\mathrm{CM}}^{d}$ $>0.8$
(where $\theta^{d}_{\mathrm{CM}}$ is the centre-of-mass polar angle of the
deuteron) is shown in fig. 5. The data peaks at $W\sim 2650$ MeV with a cross
section of 4 nb/sr. This is approximately an order of magnitude higher than
the model prediction of Fix, Arenhövel and Egerov [31, 32] which assumed
coherent production off the deuteron, where at forward angles the cross
section falls very quickly due to the increasing momentum transfer. For
$\cos\theta_{\mathrm{CM}}^{d}$ $>0.8$ at $W=2300$ and 2800 MeV, the three-
momentum transfer to the deuteron is 0.4 and 1.0 GeV/c respectively, which is
much higher than the Fermi momentum of the constituent nucleons (typically 80
MeV/c) and therefore what can be transferred to the deuteron for it to remain
intact.

Figure 5: $\gamma d\rightarrow\pi^{0}\pi^{0}d$ differential cross section for
$\cos\theta_{\mathrm{CM}}^{d}$$>0.8$. Systematic errors are the grey bars on
the abscissa. Superimposed is the model prediction from Fix, Arenhövel and
Egorov [31, 32] scaled by a factor of five (green line), and the toy pickup
model set at an arbitrary scale (red line, see text for details).

Fig. 5 also shows the results from a toy pickup model depicted in fig. 6(a),
where $\Delta^{0}\pi^{+}$ are produced in a quasi-free reaction off the proton
and the $\pi^{+}$ subsequently rescatters off the spectator neutron, producing
a $\Delta^{+}$. This model uses the same diagram as a conventional
$\Delta\Delta$ excitation by $t$-channel meson exchange that was previously
suggested as a mechanism for the ABC effect [41, 42, 43, 44]: a low mass
enhancement of the $\pi^{0}\pi^{0}$ invariant mass first reported in ref.
[45], and later attributed to the $d^{*}(2380)$. In the model, the proton and
neutron from the decays of the $\Delta$s coalesce to form the deuteron in the
final state if their relative momentum is smaller than a momentum sampled from
the internal Fermi momentum distribution of the deuteron. On-shell momentum
and energy conservation is assumed, with the Breit-Wigner mass and width of
the $\Delta$ and the spin dependent $(1+3\cos^{2}\theta)$ decay angular
distribution. A uniform initial state photon energy distribution was used, and
the output was weighted according to the pion exchange propagator,
$1/(m_{\pi}^{2}+q^{2})^{2}$, where $q$ is the pion momentum in the centre-of-
mass frame of the $\pi^{+}$ $n_{\mathrm{spectator}}$ system. An additional
$q^{2}$ weighting was assumed for the dominant magnetic coupling at the
$\gamma p\Delta^{0}\pi^{+}$ vertex. The model is simplistic and used only to
judge if this (or a similar) mechanism could be responsible for the observed
spectrum. The toy model produces a broad distribution which peaks roughly at
twice the $\Delta$ mass. The distribution is at an arbitrary scale and would
need to be weighted by the $\Delta\pi$ photoproduction cross section for an
accurate determination, however it is clear that the width and shape cannot
describe the $\gamma d\rightarrow\pi^{0}\pi^{0}d$ excitation spectrum. Other
combinations of higher lying resonances were also investigated, where the mass
distributions became broader and moved to higher energies, but also could not
describe the data.

Figure 6: Possible mechanisms contributing to the $\gamma
d\rightarrow\pi^{0}\pi^{0}d$ reaction. (a) The Toy pickup model described in
the text. (b) A sequential dibaryon decay mechanism.

Fig. 7 shows the same differential cross section data but fitted with the
$d^{*}(2380)$ and the two additional isoscalar dibaryons reported by the ELPH
collaboration [34]. A Breit-Wigner function is assumed for each dibaryon
candidate, with masses and widths fixed from ref. [34] but not their relative
amplitudes. A term proportional to the square of the centre-of-mass momentum
is also included. This is a small contribution, rising slowly with $W$ and
approximates the phase space available to the deuteron and two $\pi^{0}$ in
the final state. The fit used the statistical and fitting systematic
uncertainties summed in quadrature and achieved a reduced $\chi^{2}$ of 0.96.
The fit therefore is consistent the three dibaryon scenario however with
limited statistical precision and resolution in $W$. The inclusion of the
$d^{*}(2380)$ is also tentative, as a reasonable reduced $\chi^{2}$ of 1.35 is
achieved if it is omitted and any small enhancement at low $W$ could also be
described by variations of the toy pickup model described above. This
excitation spectrum may also yet be described by OPE models driven by triangle
singularities [38], however such calculations are beyond the scope of this
paper.

The blue squares and fit showing a proposed sequential dibaryon decay are
discussed in sec. 3.3.

Figure 7: $\gamma d\rightarrow\pi^{0}\pi^{0}d$ differential cross section for
$\cos\theta_{\mathrm{CM}}^{d}$$>0.8$ (the same as in fig. 5). A fit including
three Breit-Wigner functions (BW) are shown as the red lines, with the fixed
masses and widths labelled inset. The additional small centre-of-mass momentum
term described in the text is not shown. The blue square data points are the
differential cross section for the first of the sequential dibaryon candidate
determined from the $\pi^{0}d$ invariant mass distributions shown in fig. 9
with a Breit-Wigner function fitted and the mass and width labelled inset.
Only the statistical uncertainties are shown.

###  3.2 $\pi^{0}\pi^{0}$ and $\pi^{0}d$ invariant mass distributions over the
$d^{*}(2380)$ mass range

Fig. 8 shows the invariant mass of the $\pi^{0}d$ and $\pi^{0}\pi^{0}$ systems
for $W$ from 2270 to 2441 MeV (corresponding to five tagger channels). This is
centred over the $d^{*}(2380)$ mass, covering approximately 90 % of the Breit-
Wigner function in fig. 7. Background from quasi-free production off the
proton was subtracted by scaling the background according to the fit to the
$\pi^{0}\pi^{0}$ missing mass spectra in fig. 2. Only events within $3\sigma$
of the deuteron mass in the missing mass spectra were included to reduce the
amount of quasi-free background required to be subtracted. An additional
systematic uncertainty of 10 % for the subtraction of this background was
estimated by varying the selection cut, giving a total of 13.5 % when combined
in quadrature with the other systematic uncertainties.

Figure 8: Differential cross section versus the invariant mass of the (a)
$\pi^{0}d$ and (b) $\pi^{0}\pi^{0}$ system for $\cos\theta_{\mathrm{CM}}^{d}$
$>0.8$ and $W$ from 2270 to 2441 MeV. The measured data are the black data
points and the systematic uncertainties are the grey bars on the absisca. The
green line is the phase space distribution with an integral equal to the
measured data. The $\pi^{0}\pi^{0}$ invariant mass includes the distribution
from the ABC effect [46] (magenta line), with a scale fixed by the second data
point at 0.29 GeV/c2.

The phase space distribution gives a reasonable description of the $\pi^{0}d$
invariant mass in fig. 8(a). The $\pi^{0}\pi^{0}$ invariant mass in fig. 8(b)
exhibits a low mass enhancement and a dip at approximately 0.34 GeV/c2 which
appears similar to the ABC effect in refs. [7, 8] which was attributed to the
$d^{*}(2380)$. The pink line (scaled to match the data at 0.29 GeV/c2) shows
the expected distribution in double pionic fusion to 4He from ref. [46], where
an intermediate $d^{*}(2380)$ is formed. Qualitatively there is a good
agreement and is preferred over the phase space spectrum, despite limited
statistical precision. A fit was made (not shown) to the data including this
“ABC spectrum” and a Gaussian distribution centred around 0.4 GeV/c2 to
describe the higher invariant mass range. A differential cross section for
$\gamma d\rightarrow d^{*}(2380)\rightarrow\pi^{0}\pi^{0}d$ was subsequently
determined as $(17\pm 7_{\mathrm{stat}}\pm 4_{\mathrm{sys}})$ nb/sr for
$\cos\theta_{\mathrm{CM}}^{d}$$>0.8$. The systematic uncertainty was estimated
by varying fit parameters to the higher mass part of the spectrum and also by
using a modified ABC spectrum from the reaction $pn\rightarrow d\pi\pi$ [47].
The angular distribution of the $d^{*}(2380)$ is already well determined in
fusion reactions (see for example, ref. [7]), and so the total cross section
can be extrapolated to $(8.7\pm 3.5_{\mathrm{stat}}\pm 2.1_{\mathrm{sys}})$
nb. With improved statistics, this fitting method to the ABC effect can enable
a particularly accurate $d^{*}(2380)$ photoproduction cross section
measurement.

Alternative descriptions of the ABC effect via OPE models driven by triangle
diagrams may emerge [38], and it is anticipated this data and planned future
datasets to improve statistical precision can help discriminate between these
different interpretations.

###  3.3 Proposed sequential dibaryon decay for $W>2500$ MeV

For $W$ higher than 2500 MeV, a double peaking structure is observed in the
$\pi^{0}d$ invariant mass (fig. 9 (left panels)),which is similar to what was
observed by the ELPH Collaboration [34], where it was interpreted as an
isovector dibaryon with a mass of $2140\pm 11$ MeV/c2 and a width of $91\pm
11$ MeV/c2 from the decay of an isoscalar dibaryon. This is depicted in Fig.
6(b), where the $N(1680)N$ and $\Delta N$ configurations are proposed as
dibaryons in the sequential decay. The reaction mechanism was input to the
BGOOD simulation where the mass and width of the isoscalar dibaryon was varied
until a comparison to the real data over multiple $W$ intervals achieved a
minimal $\chi^{2}$. Shown in fig. 9 (left panels) as the blue line with an
additional phase space contribution in green, a mass of 2117 MeV/c2 and a
width of 20 MeV/c2 proved optimal. The higher energy broader peak is the
reflection of the uncorrelated $\pi^{0}d$ combination. A systematic
uncertainty of 10 MeV/c2 is attributed to the mass due to the energy
calibration the BGO Rugby Ball. The mass is comparable to the mass of 2140
MeV/c2 determined by the ELPH collaboration when accounting for the
statistical and systematic uncertainties in both measurements. The measured 20
MeV/c2 width of the dibaryon candidate is approximately the same as the
experimental $\pi^{0}d$ mass resolution and can therefore be considered an
upper limit. This is much narrower than the width of 91 MeV/c2 reported by the
ELPH collaboration, which cannot be accounted for in this data and cannot give
a satisfactory fit. The fact that the peak is below the $N\Delta$ threshold
and is much narrower than the $\Delta$ may suggest this is a genuine dibaryon
state. It is worth noting a similar peak with a mass of 2130 MeV/c2 and width
50 MeV/c2 was observed in the $\pi^{+}d$ system in previous bubble chamber
data but with limited statistics [48], where it was also suggested to be a
$\pi d$ resonance. Alternative mechanisms however can not be excluded at
present, and may yet explain the small deviation in mass and significantly
different width of the observed structure.

Figure 9: Differential cross section versus the invariant mass of the (left
panels) $\pi^{0}d$ and (right panels) $\pi^{0}\pi^{0}$ system for
$\cos\theta_{\mathrm{CM}}^{d}$ $>0.8$ and $W$ intervals labelled inset
(corresponding to four tagger channels each). The measured data are the black
data points and the systematic uncertainties are the grey bars on the absisca.
The fitted distribution to the $\pi^{0}d$ invariant mass (red line) is
comprised of phase space (green line) and proposed sequential dibaryon decay
(blue line) contributions.

The dibaryon sequential decay and phase space contributions are superimposed
onto the $\pi^{0}\pi^{0}$ invariant mass distributions in fig. 9 (right
panels). Unsurprisingly, the two distributions are similar as the simulation
used relative $s$-wave distributions for the dibaryon decays. Whilst this
distribution describes the data reasonably well, in all three $W$ intervals
there appears some structure between 0.4 to 0.5 GeV/c2 which may hint at
dynamics beyond relative $s$-wave decays, albeit limited by statistical
precision. This structure in the $\pi^{0}\pi^{0}$ invariant mass distribution
can be more clearly seen in fig. 10, where $W$ extends over the full range of
the sequential dibaryon decay candidate. A dip at 0.45 GeV/c2 is observed,
where the strength reduces by approximately 50 %. The modelled distribution
gives a poor description with a reduced $\chi^{2}$ of 2.55, however an
improvement can be made by including dibaryon decays with higher angular
momenta. Shown in fig. 10, changing the angular momentum of both the first and
second dibaryon decays to relative $p$-wave for example gives an improved fit
with a reduced $\chi^{2}$ of 1.78. This is only a demonstration of the impact
the angular distributions can have and a detailed analysis requires higher
statistics.

Figure 10: $\pi^{0}\pi^{0}$ invariant mass distribution for 2523 $<W<2738$ MeV
(over the proposed sequential dibaryon decay) . The shaded red distribution is
comprised of the phase space (thick green line) and two sequential dibaryon
decays both in $s$-wave (shaded blue area). The thick red line is comprised of
the phase space and the two sequential dibaryon decays both in $p$-wave (thick
blue line) .

The $\pi^{0}d$ invariant mass for each $W$ interval was fitted to extract the
differential cross section of the proposed sequential decay in fig. 6(b),
shown as the blue squares in fig. 7. A Breit-Wigner function was fitted,
extracting a mass of ($2618\pm 14$) MeV/c2 and width of ($148\pm 29$) MeV/c2
with a reduced $\chi^{2}$ of 1.18. Both the mass and width agrees with the
proposed highest mass dibaryon from the ELPH collaboration [34] (also fitted
in fig. 7). It is therefore suggested that two decay modes of an $N\Delta$
isoscalar dibaryon are observed, either directly to $\pi^{0}\pi^{0}d$ or to
$\pi^{0}\mathcal{D}_{12}$, where the $\mathcal{D}_{12}$ is an $N\Delta$
configuration. Proposed in the diagram in fig. 6(b), the mass and width are
both compatible with a bound $N(1680)5/2^{+}N$ dibaryon system. Numerous other
four star resonances101010As defined by the Particle Data Group [49]. have a
similar mass to the $N(1680)$, however couple less strongly to $\pi N$ and
with smaller photoproduction cross sections. The positive parity also requires
a decay with odd relative angular momentum to the $N\Delta$ $\pi^{0}$ system,
which is both supported by the $\pi^{0}\pi^{0}$ invariant mass distribution,
where two relative $p$-wave decays give an improved description, and by the
change in spin required of the constituents in the two dibaryon candidates.

##  4 Conclusions

The differential cross section for the coherent reaction, $\gamma
d\rightarrow\pi^{0}\pi^{0}d$ has been measured at BGOOD for
$\cos\theta_{\mathrm{CM}}^{d}$$>0.8$ via the identification of
$\pi^{0}\rightarrow\gamma\gamma$ in the BGO Rugby Ball and the deuteron in the
Forward Spectrometer. In this kinematic regime, the momentum transferred to
the deuteron is significantly higher than the internal Fermi momentum and
prevents the strength and shape of the excitation spectra being described by
the conventional coherent photoproduction models of Fix, Arenhövel and Egorov
[31, 32], or by a simplistic toy model of quasi-free production where final
state nucleons coalesce to form the deuteron via $\Delta\pi$ rescattering.
Instead, the data supports the findings of the ELPH collaboration, with
proposed isoscalar dibaryon states at 2380, 2470 and 2630 MeV/c2. The lowest
of these states agrees in mass and width to the $d^{*}(2380)$ dibaryon. Finer
$W$ resolution and higher statistical precision is needed for confirmation of
these states however.

A low mass enhancement in the $\pi^{0}\pi^{0}$ invariant mass distribution at
low $W$ agrees with the ABC effect and adds credence to the $d^{*}(2380)$
contributing to the excitation spectrum. With knowledge of the $d^{*}(2380)$
angular distribution, fitting to this spectrum enables an extraction of the
$\gamma d\rightarrow d^{*}(2380)\rightarrow\pi^{0}\pi^{0}d$ cross section of
$(8.7\pm 3.5_{\mathrm{stat}}\pm 2.1_{\mathrm{sys}})$ nb.

A peak in the $\pi^{0}d$ invariant mass at 2114 MeV/c2 supports the case of
the isovector dibaryon reported by the ELPH collaboration. A striking feature
of this new data is the narrow width of less than 20 MeV/c2, conflicting with
the 91 MeV/c2 width reported at ELPH. The differential cross section with
respect to $W$ is described well by the proposed 2630 MeV/$c^{2}$ dibaryon and
supports the sequential decay of an isoscalar to an isovector dibaryon decay
mechanism. The $\pi^{0}\pi^{0}$ invariant mass also exhibits structure which
suggests that the $\pi^{0}$ are beyond relative $s$-wave, maybe indicating a
high-spin component of the isoscalar state, for example, a quasi-bound
N(1680)5/2+N system.

## Acknowledgements

We would like to thank the staff and shift-students of the ELSA accelerator
for providing an excellent beam. We thank Alexander Fix, Mikhail Bashkanov,
Daniel Watts and Eulogio Oset for insightful comments.

This work is supported by SFB/TR-16, DFG project numbers 388979758 and
405882627, the RSF grant number 19-42-04132, and the Third Scientific
Committee of the INFN. This publication is part of a project that has received
funding from the European Union’s Horizon 2020 research and innovation
programme under grant agreement STRONG – 2020 - No. 824093. P. L. Cole
gratefully acknowledges the support from the U.S. National Science Foundation
(NSF-PHY-1307340, NSF-PHY-1615146, and NSF-PHY-2012826) and the Fulbright U.S.
Scholar Program (2014/2015).

## References

  * [1] F. J. Dyson and N.-H. Xuong.  Phys. Rev. Lett., 13(26):815, 1964.
  * [2] H. Clement.  Prog. Part. Nucl. Phys., 93:195, 2017.
  * [3] M. P. Locher, M. E. Sainio and A. Svarc,  Adv. Nucl. Phys., 17:47, 1986.
  * [4] I. I. Strakovsky,  Fiz. Elem. Chast. At. Yadra, 22:615, 1991.
  * [5] Tsu Yao.  Phys. Rev. B, 134:454, 1964.
  * [6] I. Bar-Nir.  Nucl. Phys. B, 54:17, 1973.
  * [7] P. Adlarson et al.  Phys. Rev. Lett., 106:242302, 2011.
  * [8] M. Bashkanov.  Phys. Rev. Lett., 102:052301, 2009.
  * [9] H. Clement and T. Skorodko.  Chin. Phys. C, 45:022001, 2021.
  * [10] P. Adlarson et al.  Phys. Lett. B, 721:229, 2013.
  * [11] P. Adlarson et al.  Phys. Rev. Lett., 112:202301, 2014.
  * [12] P. Adlarson et al.  Phys. Rev. C, 90:035204, 2014.
  * [13] M. Bashkanov, S. Kay, D. P. Watts et al.  Phys. Lett. B, 789:7, 2019.
  * [14] M. Bashkanov, D. P. Watts, S. J. D. Kay et al.  Phys. Rev. Lett., 124:132001, 2020.
  * [15] P. Adlarson et al.  Phys Rev. C, 88:055208, 2013.
  * [16] P.Adlarson et al.  Phys. Lett. B, 743:325, 2015.
  * [17] T. Goldman, K. Maltman, G. J. Stephenson, K. E. Schmidt, and Fan Wang.  Phys Rev. C, 39(5):1889, 1989.
  * [18] A. Gal and H. Garcilazo.  Nuc. Phys. A, 928:73, 2014.
  * [19] A.V. Kravtsov, M.G. Ryskin, and I.I. Strakovsky.  J. Phys. G: Nucl. Part. Phys., 9:L187, 1983.
  * [20] N. Hoshizaki.  Phys. Rev. C, 45(R1424), 1992.
  * [21] N. Hoshizaki.  Prog. Theor. Phys., 89:563, 1993.
  * [22] R.A. Arndt, J.S. Hyslop, and L.D. Roper.  Phys. Rev. D, 35:128, 1987.
  * [23] C. .H. Oh, R. A. Arndt, I. I. Strakovsky, and R. L. Workman.  Phys. Rev. C 56:635, 1997.
  * [24] R. Molina, N. Ikeno, and E. Oset.  arXiv:2102.05575v3 [nucl-th], 2021.
  * [25] N. Ikeno, R. Molina, and E. Oset.  Phys Rev. C, 104:014614, 2021.
  * [26] M. Bashkanov, H. Clement and T. Skorodko.  arXiv:2106.00494v2 [nucl-th], 2021.
  * [27] P. Shen Y. Dong, F. Huang and Z. Zhang.  Phys Rev. C, 94:014003, 2016.
  * [28] I. Vidaña, M. Bashkanov, D. P. Watts, and A. Pastore.  Phys. Lett. B, 781:112, 2018.
  * [29] M. Bashkanov and D. P. Watts.  J. Phys. G: Nucl. Part. Phys, 47:03LT01, 2020.
  * [30] T. Ishikawa et al.  Phys. Lett. B, 772:398, 2017.
  * [31] A. Fix and H. Arenhövel.  Eur. Phys. J. A, 25:115, 2005.
  * [32] M. Egorov and A. Fix.  Nucl. Phys. A, 933:104, 2015.
  * [33] M. S. Guenther.  PoS(Hadron2017), 210:051, 2018.
  * [34] T. Ishikawa et al.  Phys. Lett. B, 789:413, 2019.
  * [35] S Alef et al.  Eur. Phys. J. A, 56:104, 2020.
  * [36] W. Hillert.  Eur. Phys. J. A, 28:139, 2006.
  * [37] W. Hillert et al.  EPJ Web Conf., 134:05002, 2017.
  * [38] E. Oset.  private communication.  2021\.
  * [39] S. Alef et al.  Eur. Phys. J. A, 57:80, 2021.
  * [40] V. L Kashevarov et al.  Phys. Rev. C, 85:064610, 2012.
  * [41] T. Risser and M. D. Shuster.  Phys. Lett. B, 43:68, 1973.
  * [42] T. Risser I. Bar-Nir and M. D. Shuster.  Nucl. Phys. B, 87(109), 1975.
  * [43] E. Oset L. Alvarez-Ruso and E. Hernandez.  Nucl. Phys. A, 633:519, 1998.
  * [44] L. Alvarez-Ruso.  Phys. Lett. B, 452:207, 1999.
  * [45] N. E. Booth, A. Abashian, and K. M. Crowe.  Phys. Rev. Lett., 7(1):35, 1961.
  * [46] P. Adlarson et al.  Phys. Rev. C, 86:032201, 2012.
  * [47] M. Bashkanov.  Private communication, 2021.
  * [48] D. C. Brunt, M. J. Clayton, and B. A. Westwood.  Phys. Lett. B, 26:317, 1968.
  * [49] P.A. Zyla et al. (Particle Data Group).  Prog. Theor. Exp. Phys., page 083C01, 2020.
  * [50] R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla.  Phys. Rev. C, 51(1):40, 1995.