#  Variational Quantum Eigensolver with Reduced Circuit Complexity

Yu Zhang<EMAIL_ADDRESS>Theoretical Division, Los Alamos National Laboratory,
Los Alamos, NM, 87545, USA  Lukasz Cincio  Theoretical Division, Los Alamos
National Laboratory, Los Alamos, NM, 87545, USA  Christian F. A. Negre
Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, 87545,
USA  Piotr Czarnik  Theoretical Division, Los Alamos National Laboratory, Los
Alamos, NM, 87545, USA  Patrick Coles  Theoretical Division, Los Alamos
National Laboratory, Los Alamos, NM, 87545, USA  Petr M. Anisimov
Accelerators and Electrodynamics Group, Los Alamos National Laboratory, Los
Alamos, New Mexico 87545, USA  Susan M. Mniszewski  Computer, Computational
and Statistical Sciences Division, Los Alamos National Laboratory, Los Alamos,
New Mexico 87545, USA  Sergei Tretiak  Theoretical Division, Los Alamos
National Laboratory, Los Alamos, NM, 87545, USA  Center for Integrated
Nanotechnologies, Los Alamos National Laboratory, Los Alamos, New Mexico
87545, USA  Pavel A. Dub  Chemistry Division, Los Alamos National Laboratory,
Los Alamos, NM, 87545, USA

###### Abstract

The variational quantum eigensolver (VQE) is one of the most promising
algorithms to find eigenvalues and eigenvectors of a given Hamiltonian on
noisy intermediate-scale quantum (NISQ) devices. A particular application is
to obtain ground or excited states of molecules. The practical realization is
currently limited by the complexity of quantum circuits. Here we present a
novel approach to reduce quantum circuit complexity in VQE for electronic
structure calculations. Our algorithm, called ClusterVQE, splits the initial
qubit space into subspaces (qubit clusters) which are further distributed on
individual (shallower) quantum circuits. The clusters are obtained based on
quantum mutual information reflecting maximal entanglement between qubits,
whereas entanglement between different clusters is taken into account via a
new “dressed” Hamiltonian. ClusterVQE therefore allows exact simulation of the
problem by using fewer qubits and shallower circuit depth compared to standard
VQE at the cost of additional classical resources. In addition, a new gradient
measurement method without using an ancillary qubit is also developed in this
work. Proof-of-principle demonstrations are presented for several molecular
systems based on quantum simulators as well as an IBM quantum device with
accompanying error mitigation. The efficiency of the new algorithm is
comparable to or even improved over qubit-ADAPT-VQE and iterative Qubit
Coupled Cluster (iQCC), state-of-the-art circuit-efficient VQE methods to
obtain variational ground state energies of molecules on NISQ hardware. Above
all, the new ClusterVQE algorithm simultaneously reduces the number of qubits
and circuit depth, making it a potential leader for quantum chemistry
simulations on NISQ devices.

##  I Introduction

One of the major goals of computational chemistry is the development of
methods and algorithms for the calculation of the molecular electronic ground
and excited state energies and corresponding wave functions from first-
principles. Such eigenvalues and eigenvectors can be obtained from the
solution of the time-independent electronic Schrödinger equation. However,
since the invention of classical digital computers in the early 1940s, the
exact numerical solution of this central quantum mechanical equation remains
infeasible for systems having more than 12 electrons distributed on 184 spin-
orbitals Gan and Harrison (2005). The reason is that the solution space in
this approximation (i.e., the Fock space) grows factorially with the system
size (e.g., the number of electrons and basis functions). One of the most
promising and immediate applications of quantum computers is solving
classically intractable quantum chemistry problems Cao  _et al._ (2019);
McArdle  _et al._ (2020).

The quantum phase estimation (QPE) algorithm Kitaev (1995) represents the
natural translation of the full configuration interaction (FCI) procedure to
quantum computers Aspuru-Guzik  _et al._ (2005). However, QPE requires
millions of qubits and quantum gates even for relatively small systems, making
the algorithm unsuitable for practical applications on near-term noisy
intermediate-scale quantum (NISQ) devices Preskill (2018). Full quantum
eigensolver (FQE), in which the complexity of basic gate operators is
polylogarithmic in the number of spin-orbitals, represents another fully
quantum algorithm to simulate molecular systems Wei  _et al._ (2020). However,
FQE might not be suitable for advantageous quantum chemistry demonstrations
due to the limitations of NISQ devices. Indeed, leading digital quantum
computers based on superconducting qubit technology have limited coherence
times ($\sim$ 100 $\mu$s) and gate error rates ($\sim$ 2 $\times 10^{-2}$ for
a two-qubit gate), limiting the possible number of operations that can be
executed to evolve the quantum state. Under these conditions, leveraging
classical resources as much as possible through a hybrid quantum-classical
approach seems to be the most promising route toward achieving quantum
chemical advantage on NISQ devices Moll  _et al._ (2018); Cerezo  _et al._
(2020); Bharti  _et al._ (2021).

The variational quantum eigensolver (VQE) is one of the leading hybrid
quantum-classical algorithms developed specifically to simulate molecular
systems in their ground Peruzzo  _et al._ (2014); McClean  _et al._ (2016); Li
_et al._ (2019a); Fedorov  _et al._ (2021) and ultimately excited states Shen
_et al._ (2017); McClean  _et al._ (2017); Nakanishi  _et al._ (2019); Higgott
_et al._ (2019); Kawai and Nakagawa (2020); Greene-Diniz and Muñoz Ramo
(2021). In a typical VQE setup, a variational parameterized ansatz is used to
represent the trial wave function prepared with a quantum circuit and the
expectation value of the qubit Hamiltonian is measured. Then, the parameters
of the ansatz are iteratively optimized on a classical computer using the
Rayleigh-Ritz variational principle. Although VQE simulations of small
molecules have been performed on various quantum architectures such as photons
Peruzzo  _et al._ (2014), superconducting qubits Kandala  _et al._ (2019,
2017) and trapped ions Tang  _et al._ (2021a), major efforts are needed to
scale up this approach to larger molecular systems of chemical interest. Here,
albeit not as large as QPE circuits, VQE would also suffer from the large size
of the quantum circuits, which coupled to a classical optimization requiring
many variational parameters, can render calculations intractable.
Consequently, the construction of an efficacious ansatz and improved classical
optimizer is an active area of research Bharti  _et al._ (2021); Cerezo  _et
al._ (2020).

There are two main approaches for ansatz design for electronic structure
calculations. Hardware-efficient ansatzes Kandala  _et al._ (2019, 2017) are
constructed from repeated, dense blocks of a limited set of parameterized
gates that are easy to implement on a quantum device. While such an ansatz
requires fewer resources on the quantum processor, because of its totally
chemically agnostic nature it might face difficulties in the parameter
optimization due to the barren-plateau problem McClean  _et al._ (2018); Wang
_et al._ (2021); Cerezo  _et al._ (2021); Abbas  _et al._ (2020).
Consequently, simulation of large molecular systems with hardware-efficient
ansatz might become practically challenging and even impossible. Chemically-
inspired ansatzes, such as the Unitary Coupled-Cluster Singles and Doubles
(UCCSD) Bartlett  _et al._ (1989); Bartlett and Musiał (2007); Seeley  _et
al._ (2012a); Barkoutsos  _et al._ (2018); Romero  _et al._ (2018), are
designed by using the domain knowledge from classical quantum chemistry. The
standard version of UCCSD, however, has unfavorable scaling in the number of
gates required for larger molecular systems Kühn  _et al._ (2019).
Consequently, various approaches are being developed to improve on the UCCSD
method to obtain shorter circuits or obtain higher than UCCSD accuracy at
comparable circuit depth Fedorov  _et al._ (2021); Ryabinkin  _et al._ (2018);
Lee  _et al._ (2019); Matsuzawa and Kurashige (2020); Xia and Kais (2020);
Tkachenko  _et al._ (2021). Among them, “iterative VQE” algorithms Grimsley
_et al._ (2019); Tang  _et al._ (2021b); Ryabinkin  _et al._ (2020); Lang  _et
al._ (2021); Ryabinkin  _et al._ (2021); Yordanov  _et al._ (2020), which
instead of using a fixed ansatz, construct problem-tailored ansatzes on-the-
fly, have gained attention because of their controllable circuit-size.

Two notable iterative VQE algorithms to reduce quantum circuit depth are
(fermionic and qubit) ADAPT-VQE Grimsley  _et al._ (2019); Tang  _et al._
(2021b) and iterative Qubit Coupled Cluster (iQCC) Ryabinkin  _et al._ (2020);
Lang  _et al._ (2021); Ryabinkin  _et al._ (2021). Fermionic-ADAPT-VQE employs
a predetermined pool of spin-adapted fermionic (generalized) singles and
doubles excitation operators from which the ansatz is dynamically constructed
Grimsley  _et al._ (2019). The more measurement- and circuit-efficient (e.g.
less CNOT gates in particular) qubit-ADAPT-VQE uses selected Pauli words
(obtained from mapping fermionic operators into Pauli operators) to form the
“qubit” pool, but the algorithm comes with a larger number of ADAPT iterations
because of the larger number of variational parameters (more gradient
calculations required) Tang  _et al._ (2021b). IQCC uses arbitrarily shallow
quantum circuit depth at the expense of an iterative canonical Hamiltonian
dressing technique employing the qubit coupled cluster (QCC) ansatz. The
exponential growth of the “dressed” Hamiltonian with the size of the QCC
transformation in an anti-commuting set could be mitigated using the
involuntary linear combinations of Paulis technique Lang  _et al._ (2021). But
the algorithm is exponential with respect to the dressing steps.

Although circuit depth is one of the most critical metrics for NISQ devices.
The number of physical qubits on the actual quantum chip (circuit width) is
the necessary condition that limits the problem that can be solved. Several
techniques based on symmetries presented in the Hamiltonian have been used to
reduce the number of qubits Bravyi  _et al._ (2017); Setia  _et al._ (2020);
Yen  _et al._ (2019); Zhang  _et al._ (2021a). But these techniques can only
reduce a few qubits and have certain symmetry requirements. In this work, we
present the first-ever algorithm which reduces both quantum circuit depth and
width in VQE at the cost of additional classical resources. Our novel
algorithm, called ClusterVQE, uses the mutual information (measure of
correlation between spin-orbitals Rissler  _et al._ (2006); Huang and Kais
(2005)) to group the qubits into different clusters, which can be distributed
to different quantum circuits. The correlation between different clusters is
minimized by the mutual information optimization method and is taken into
account by building a “dressed” Hamiltonian. Even shallower circuit depths can
be achieved compared to the state-of-the-art qubit-ADAPT-VQE. Consequently,
the new algorithm is more robust to noise making it particularly attractive
for NISQ devices. Because the correlations between different clusters are
minimized, the dressed Hamiltonian protocol only introduces minor additional
classical resources requirements compared to the iQCC method. Most
importantly, the new ClusterVQE algorithm removes the entanglement between
different clusters via the dressed Hamiltonian. Consequently, ClusterVQE can
solve the original problem in a much smaller qubit space (clusters), reducing
the number of qubits significantly for simulating large molecules on NISQ
devices.

Figure 1:  Finding the lowest eigenvalue of a given qubit Hamiltonian
$\hat{H}$ with VQE and ClusterVQE: VQE uses a quantum circuit defined on $m$N
qubits ($m=2$ for the example used), whereas ClusterVQE uses N qubits and
shorter circuits at the expense of a dressed Hamiltonian $\hat{H}_{d}$.

##  II Methodology

###  II.1 Mutual Information based Clustering of Qubits

The first essential step of ClusterVQE is to split the qubits into different
clusters and minimize the correlations between different clusters. With this
aim, an entanglement map based on the mutual information (MI) that reflects
the electronic correlations among all pairs of qubits is calculated Tkachenko
_et al._ (2021). Previous results obtained for the orbital ordering problem in
the Density Matrix Renormalization Group (DMRG) method in classical quantum
chemistry calculations Rissler  _et al._ (2006) showed that the MI is a
reliable parameter to quantify the correlation between two quantum particles.
Our previous work also demonstrated that an approximate MI is sufficient for
the qubit re-arrangement Tkachenko  _et al._ (2021) for the reduction of
circuit depth. The MI has also been used to construct a compact entangler pool
for ADAPT-VQE Zhang  _et al._ (2021b). According to quantum information
theory, the MI between qubits $i$ and $j$ is defined as follows Rissler  _et
al._ (2006):

$I_{ij}=\frac{1}{2}(S_{i}+S_{j}-S_{ij})(1-\delta_{ij}),$ (1)

where $S_{i}$ and $S_{ij}$ are the single- and two-qubit Von Neumann
entropies, respectively, and $\delta_{ij}$ is the Kronecker delta. $S_{i}$ and
$S_{ij}$ are obtained from the corresponding one- and two-body density matrix
Tkachenko  _et al._ (2021).

Once the MI is obtained, it can be used to cluster the qubits into different
groups by minimizing the MI between different clusters which is defined as the
summation over MIs between inter-cluster qubits. The qubit clustering can be
achieved by using a classical graph partitioning method, such as the graph
partition function of the Metis library https://metis.readthedocs.io/en/latest
(Appendix. I.1). Alternatively, this can be mapped into a quadratic
unconstrained binary optimization (QUBO) problem and solved on a quantum
annealer by either the MI-based clustering method (details are given Appendix
I.2) or the quantum community detection method (QCD) Mniszewski  _et al._
(2021) (details are given in Appendix I.3). Problems with up to 64 variables
or nodes can be solved directly on the D-Wave 2000Q quantum annealer. Larger
sizes require a quantum-classical approach using _qbsolv_ Booth  _et al._
(2017) \+ D-Wave.

###  II.2 The ClusterVQE algorithm

Within the VQE algorithm, a parameterized ansatz
$\ket{\Phi}=\hat{U}(\theta)\ket{\Phi_{r}}$ is suggested and encoded on the
quantum circuit, where $\ket{\Phi_{r}}$ is the reference state. In this work,
the mean-field Hartree-Fock (HF) state is used as a reference. The unitary
operator $\hat{U}(\theta)$ introduces the correlations between qubits. Since
current quantum gates can only operate on a few qubits at a time,
$\hat{U}(\theta)$ is decomposed as a tensor product of many unitary operators,
denoted as entanglers, with each entangler acting on a few qubits, i.e.,
$\hat{U}(\theta)=\prod^{\mathcal{D}}_{j}\hat{U}(j)$ where
$\hat{U}_{j}=e^{i\theta_{j}\hat{A}_{j}}$ are the entanglers and $\mathcal{D}$
is the number of entanglers in the ansatz. In the widely used UCCSD ansatz,
the $\\{\hat{A}_{j}\\}$ can be readily obtained by mapping the single/double
Fermionic excitation operators into qubit-space as shown in Appendix II. With
the ansatz, the ground-state (GS) energy can be obtained by imposing the
variational principle,

$E_{GS}=\min_{\theta}\bra{\Phi_{r}}\hat{U}^{\dagger}(\theta)\hat{H}\hat{U}(\theta)\ket{\Phi_{r}},$
(2)

where the expectation value is measured on a quantum circuit. However, the
contribution of each entangler to the GS energy is different. Thus, the UCCSD
ansatz can be made even more compact by the adaptive construction Grimsley
_et al._ (2019); Tang  _et al._ (2021a) in the qubit-ADAPT-VQE, at the cost of
more VQE calculations. This algorithm ends up with a compact sequence of
entanglers $\hat{U}(\theta)=\prod^{D}_{j=1}\hat{U}_{j}(\theta_{j})$ with
$D<\mathcal{D}$, where $\hat{U}_{j}(\theta_{j})=e^{i\theta_{j}\hat{P}_{j}}$
and the corresponding $\\{\hat{P}_{j}\\}$ is chosen from the operator pool
$\\{\hat{A}_{j}\\}$ based on their contributions to the energy weighted by the
gradients Grimsley  _et al._ (2019).

Alternatively, according to the qubit clustering introduced in Sec. II.1, the
entanglers used in the ansatz can be grouped into different sets,

$\hat{U}(\theta)=\\{\hat{U}_{c_{1}},\cdots,\hat{U}_{c_{M}},\hat{U}_{c_{12}},\cdots,\hat{U}_{c_{M-1,M}}\\}$
(3)

where $M$ is the number of clusters. In this way, the ordering of entanglers
in $\hat{U}(\theta)$ is changed and the corresponding parameters can be
different as well. The entangler set $\\{\hat{U}_{c_{i}}\\}$ only acts on the
corresponding cluster, and {$\hat{U}_{c_{ij}}\\}$ represents the entanglers
that couple different clusters. Consequently, the energy minimization in Eq. 2
can be rewritten as

$\displaystyle E_{GS}=$
$\displaystyle\min_{\theta}\bra{\Phi_{r}}\prod_{i}\hat{U}^{\dagger}_{c_{i}}\left[\prod_{i\neq
j}\hat{U}^{\dagger}_{c_{ij}}\hat{H}\hat{U}_{c_{ij}}\right]\hat{U}_{c_{i}}\ket{\Phi_{r}}$
$\displaystyle=$
$\displaystyle\min_{\theta}\bra{\Phi_{r}}\prod_{i}\hat{U}^{\dagger}_{c_{i}}\hat{H}^{d}\hat{U}_{c_{i}}\ket{\Phi_{r}}$
(4)

where $\hat{H}_{d}\equiv\prod_{i\neq
j}\hat{U}^{\dagger}_{c_{ij}}\hat{H}\hat{U}_{c_{ij}}$ is the dressed
Hamiltonian. Because the $\hat{H}_{d}$ can be represented in terms of Pauli
words and each Pauli word can be trivially rewritten as the tensor product of
different clusters, i.e.,
$\hat{H}_{d}=\sum_{k}\alpha_{k}\hat{P}_{k}=\sum_{k}\alpha_{k}\prod_{i}\hat{P}_{k}(c_{i})$
and the initial state can be decomposed as $\prod_{i}\Phi_{r}(c_{i})$, the
energy in Eq. II.2 can be further rewritten as

$E_{GS}=\min_{\theta}\sum_{k}\alpha_{k}\left[\prod_{i}\bra{\Phi_{r}(c_{i})}\hat{U}^{\dagger}_{c_{i}}\hat{P}_{k}(c_{i})\hat{U}_{c_{i}}\ket{\Phi_{r}(c_{i})}\right].$
(5)

Here each expectation value
$\bra{\Phi_{r}(c_{i})}\hat{U}^{\dagger}_{c_{i}}\hat{P}_{k}(c_{i})\hat{U}_{c_{i}}\ket{\Phi_{r}(c_{i})}$
only involves the qubits in the cluster $i$. Hence, the energy of a larger
problem can be measured on $M$ smaller quantum circuits after Hamiltonian
dressing, as shown by the schematic diagram in Figure 1. In practice, like the
qubit-ADAPT-VQE algorithm, clusterVQE adaptively grows the ansatz in Eq. II.2.
The ansatz used in the $j^{th}$ iteration is denoted as
$\hat{\mathcal{U}}_{j}$. However, in contrast to the growing ansatz in qubit-
ADAPT-VQE Grimsley  _et al._ (2019) and growing Hamiltonian in iQCC Ryabinkin
_et al._ (2020), the ansatz in ClusterVQE is decomposed onto smaller quantum
circuits, and only the entanglers between different clusters are used to dress
the Hamiltonian. Consequently, the iterative construction breaks the original
VQE algorithm with deeper and wider quantum circuits into a series of quantum
circuits with a shallower depth and a smaller number of qubits, as illustrated
by Figure 1.

The entangler $\hat{U}_{k}$ chosen in the ansatz $\hat{\mathcal{U}}_{j}$ for
the $j^{th}$ iteration is determined by its contribution to the energy
weighted by the gradient of energy with respect to the parameter, i.e., the
$\hat{U}_{k}$ terms with largest gradients $\frac{\partial
E}{\partial\theta_{k}}$ are selected Grimsley  _et al._ (2019), where the
gradient is defined as

$\displaystyle\frac{\partial E}{\partial\theta_{k}}=$
$\displaystyle\frac{\partial}{\partial\theta_{k}}\bra{\Phi_{r}}\hat{U}^{\dagger}_{k}(\theta_{k})\hat{H}_{d,j-1}\hat{U}_{k}(\theta_{k})\ket{\Phi_{r}}|_{\theta_{k}=0}$
$\displaystyle=$ $\displaystyle
i\bra{\Phi_{r}}\hat{U}^{\dagger}_{k}(\theta_{k})[\hat{H}_{d,j-1},\hat{P}_{k}]\hat{U}_{k}(\theta_{k})\ket{\Phi_{r}}.$
(6)

Where $\hat{H}_{d,j-1}$ is the dressed Hamiltonian of the $j-1$ iteration. It
should be noted that $\hat{H}_{d,j}$ may be the same as the $\hat{H}_{d,j-1}$
if the selected entangler of the $j^{th}$ iteration is not used to dress the
Hamiltonian. Since the measurement of each gradient term is independent, all
the gradients can be measured in parallel with several uncoupled quantum
computers Grimsley  _et al._ (2019).

Figure 2:  Flowchart of the ClusterVQE algorithm. The fifth step in the red-
dashed box is skipped if MI is provided.

The flowchart of the ClusterVQE scheme is presented in Figure 2. In the first
step, the qubit Hamiltonian and operator pool are prepared: a) Provided
classical mean-field calculation of the molecule, map the second quantized
Hamiltonian into its qubit counterpart by using either Jordan-Wigner (JW)
Jordan and Wigner (1928), Bravyi-Kitaev (BK) Bravyi and Kitaev (2002), or
parity Seeley  _et al._ (2012b); Bravyi  _et al._ (2017) encoding. b) Generate
the UCCSD excitation operators and map them into Pauli words. Here only Pauli
words with the different flip indices (Appendix. II) are chosen and added into
the excitation pool $\\{\hat{A}_{i}\\}$. In addition, clustering of qubits is
performed if MI is provided. Otherwise, the qubits in the first iteration are
simply grouped according to their spin index. If the largest gradient in the
second step is smaller than the convergence criteria $\epsilon$, the energy
minimization criteria are satisfied and the algorithm exits from the iterative
loop. The fifth step is only triggered if clustering on-the-fly is desired.

Since some of the entanglers in $\hat{\mathcal{U}}_{j}$ are used to dress the
Hamiltonian, the circuit depth is smaller than that of qubit-ADAPT-VQE.
Besides, because the clustering minimizes the entanglement between different
clusters, the number of dressing entanglers is much smaller than that of iQCC.
Consequently, ClusterVQE can balance the computational costs between the
Quantum Processing Unit (QPU) and CPU. The circuit depth/width is further
reduced compared to qubit-ADAPT-VQE, and exponential scaling in growing of the
Hamiltonian in iQCC is avoided. As a result, ClusterVQE is more robust to
noise, as shown and discussed later, making the ClusterVQE algorithm more
suitable for NISQ devices than qubit-ADAPT-VQE. Most importantly, compared to
previous approaches, ClusterVQE can solve larger problems on smaller quantum
circuits.

###  II.3 Analytical gradients for VQE optimization

In order to improve the algorithm performance in the presence of noise, an
analytical gradient is essential for classical optimization. It should be
noted that the gradient for VQE optimization is different from the gradients
in Eq. II.2 which is used to chose entanglers for each iteration. Ref. Li  _et
al._ (2019b) suggests calculation of the analytical gradients by performing a
measurement on an ancillary qubit Romero  _et al._ (2018); Dallaire-Demers
_et al._ (2018). A similar approach has been developed in this work for the
analytical gradient measurements but without introducing an ancillary qubit.
The gradient of energy with respect to the variational parameter reads as
follows.

$\frac{\partial
E}{\partial\theta_{i}}=2\text{Im}\left\\{\bra{\Phi_{r}}\hat{U}^{\dagger}\hat{H}\hat{P}^{\text{eff}}_{i}\hat{U}\ket{\Phi_{r}}\right\\}$
(7)

where $\hat{P}^{\text{eff}}_{i}$ is the dressed Pauli word as shown in
Appendix III. Therefore, the energy gradient can be measured in the same way
as the energy measurement by replacing $\hat{H}$ with
$\hat{H}\hat{P}^{\text{eff}}_{j}$ and no ancillary qubit is required. On a
noiseless simulator, VQE calculations with numerical gradients or analytical
gradients give the same results. However, analytical gradients are much more
stable on noisy simulators and quantum computers.

It should be noted that even though the ClusterVQE method can reduce the
number of qubits and ansatz complexity by distributing the measurements on
smaller quantum circuits, the number of measurements increases with the number
of entanglers used in the dressed Hamiltonian. To further reduce the number of
measurements, future work may consider recently developed measurement
reduction techniques Bonet-Monroig  _et al._ (2019); Yen and Izmaylov (2020);
Huang and Kueng (2019); Zhao  _et al._ (2020); Huang  _et al._ (2020); Yen
_et al._ (2020); Izmaylov  _et al._ (2019); Huggins  _et al._ (2019).

##  III Results and discussion

To demonstrate the validity of the ClusterVQE algorithm, we have simulated
several molecules on both quantum simulator and IBM quantum devices. In order
to perform the simulations, an in-house modified Qiskit code Abraham and et al
(2019) has been developed. The Jordan-Wigner operator transformation Jordan
and Wigner (1928) and the Limited-memory BFGS Bound (L-BFGS-B) optimizer
Morales and Nocedal (2011) were used in this work. The obtained results are
compared with that of exact diagonalization. In order to make the comparison
over different methods (VQE, qubit-ADAPT-VQE, iQCC) meaningful, the same
ansatz, qubit-UCCSD (details in Appendix II), is used for all simulations
below. The convergence criteria of $\epsilon=10^{-4}$ and $5\times 10^{-3}$
are used for LiH and N2, respectively, if not specified elsewhere. For the N2
molecules, the maximum number of iterations is 25, unless convergence occurs
earlier.

Figure 3:  GS Energy errors (a,c) and (b,d) number of Pauli words per each
iteration for the LiH molecule with different Li-H bond lengths, 1.547 Å(a,b)
and 2.4 Å(c,d), obtained from the qubit-ADAPT-VQE, iQCC, and ClusterVQE
algorithms. Figure 4:  (a) GS energy, (b) energy error, and (c) number of
Pauli words in the qubit Hamiltonian of LiH for different bond lengths
obtained from different methods. For VQE, one layer of the qUCCSD ansatz is
used. One layer of the qUCCSD ansatz is also used as the operator pool for the
other methods.

###  III.1 LiH molecule

We start with simulating the ground state of the LiH molecule. LiH requires 10
qubits with the minimal STO-3G basis set and frozen core orbitals. Figures
3(a) and (c) compare the energy convergence of the qubit-ADAPT-VQE,
ClusterVQE, and iQCC approaches for different Li-H bond lengths (1.547 Åand
2.4 Å). Because iQCC only implements one entangler on the quantum circuit, the
number of parameters to be optimized is also smaller compared to that of
qubit-ADAPT-VQE and ClusterVQE. Consequently, a smaller number of optimization
cycles is required per iteration. For the iQCC method, once the entanglers in
the previous iteration are used to dress the Hamiltonian, these parameters are
fixed. In contrast, the parameterized ansatz in the qubit-ADAPT-VQE and
ClusterVQE approaches are to be re-optimized at each iteration, making them
more flexible to achieve better convergence. As a result, the convergence of
iQCC is slightly worse than that of qubit-ADAPT-VQE for a given number of
iterations, as shown in Figures 3(a) and (c). For ClusterVQE, because only a
few operators are used to dress the Hamiltonian, the performance of ClusterVQE
is almost the same as that of qubit-ADAPT-VQE. Because the Hamiltonian of the
iQCC and ClusterVQE methods is dressed by the entanglers, the number of Pauli
words in the dressed Hamiltonian is increased. However, compared to iQCC, only
two operators are used to dress the Hamiltonian in ClusterVQE as shown, the
number of Pauli words of the dressed Hamiltonian does not increase
significantly. In contrast, the number of Pauli words in iQCC increases
exponentially with the number of iterations before saturation, see Figures
3(b) and (d). We also compared performance of the qubit-ADAPT-VQE, iQCC, and
ClusterVQE approaches on the ground-state potential energy surface (PES) of
LiH as shown in Figure 4. Overall, ClusterVQE can achieve the same energy
error as that of qubit-ADAPT-VQE but with smaller quantum circuits. And
ClusterVQE uses less Pauli words compared to iQCC, as shown in Figure 4.

We also compared the performance of the different clustering methods used
(Appendix. I). For most geometries of the LiH molecule, these three methods
produce the same clustering. However, at the bond length of 2.4 Å, the Metis
graph partition method results in less optimal clustering. Consequently, more
entanglers are needed to dress the Hamiltonian and ClusterVQE with the Metis
method ends up with 808 Pauli words. This comparison suggests that the QCD and
MI-based clustering methods provide better and more stable performance in
clustering and help improve the ClusterVQE performance.

Figure 5:  GS Energy errors (a,c) and number of Pauli words (b,d) per each
iteration for the N2 molecule with different bond lengths, 1.09 Åand 1.6 Å,
obtained from the qubit-ADAPT-VQE (blue), iQCC (purple) and ClusterVQE (green)
methods.

###  III.2 N2 molecule

We next demonstrated the validity of ClusterVQE in the strongly correlated
regime. The N2 molecule is a prototype system for testing multi-reference
methods in quantum chemistry due to its strong electronic correlations
particularly when the bond is stretched. Here N2 at equilibrium (1.09 Å) and
stretched configurations (1.6 Å) are used as examples. The 14-qubit N2/STO-3G
system is simulated after freezing the lowest two spin-orbitals. As shown in
Figures 5(a) and (c), the energy convergence of ClusterVQE is almost the same
as that of qubit-ADAPT-VQE and is better than that of iQCC. Moreover, the size
of the dressed Hamiltonian in ClusterVQE is much smaller than that of iQCC. It
should be noted that the stretched N2 molecule has stronger correlations,
however, accompanied with a better separation between clusters. Consequently,
less dressing is required for the stretched N2, and ClusterVQE only slightly
increases the size of the dressed Hamiltonian. In contrast, the size of the
Hamiltonian exceeds 15,000 Pauli words in iQCC after 25 iterations. These
results suggest that, even for strongly correlated systems, our ClusterVQE
reduces circuit complexity without a significant increase in Hamiltonian size.

###  III.3 Noisy simulation and clustering on-the-fly

We next examine the dissociation curve of the LiH molecule in the presence of
noise to demonstrate ClusterVQE’s resilience to this common feature of all
NISQ architectures. Figure 6 shows the dissociation curve of the LiH molecule
calculated with different algorithms. The energy values are averaged over 5
runs. As shown by the curves, in the presence of noise, the dissociation
curves obtained from the VQE, qubit-ADAPT-VQE, iQCC, and ClusterVQE methods
are different from the exact solution. Since the VQE algorithm requires more
parameters and the deepest circuit, the noise influence is the largest across
the set, showing the strongest deviation from the reference. The qubit-ADAPT-
VQE method, in contrast, is able to grow a compact ansatz with fewer
parameters and lower circuit depth. The influence of noise is suppressed. On
the other hand, the iQCC algorithm uses much shallower circuit depth. As a
result, the iQCC is most immune to the noise, and the corresponding
dissociation curve is the closest to the reference. Next is ClusterVQE which
partially dresses the Hamiltonian and grows the ansatz. The circuit depth of
ClusterVQE lies between that of iQCC and qubit-ADAPT-VQE. Consequently, the
energy error in the presence of noise is also in between that of the iQCC and
qubit-ADAPT-VQE methods.

Figure 6:  (a) Dissociation curves of LiH and (b) corresponding errors with
exact diagonalization (black), VQE (red), qubit-ADAPT-VQE (blue), iQCC
(purple), and ClusterVQE (green) methods on a noisy simulator. The custom
noise in Qiskit (Abraham and et al, 2019) was used, where the error rates of
single-qubit and two-qubit gates are 0.001 and 0.005, respectively. The means
and error bars are calculated from 5 independent simulations. The bottom panel
shows the corresponding errors.

We also compare the performance of the four algorithms with different error
rates. Figures S1(a) and (b) compare the energy errors with respect to
different single- and two-qubit error rates, respectively. The single-qubit
gates usually have much lower error rates on real quantum devices. For
example, the single- and two-qubit gate error rates of the IBM-Yorktown IBM
device are around $\sim 10^{-3}$ and $\sim 10^{-2}$, respectively.
Consequently, the error of single-qubit gates have a much smaller influence on
the results. In contrast, the results are significantly affected by the error
of two-qubit gates. In general, the error of ClusterVQE lies between that of
the qubit-ADAPT-VQE and iQCC methods on the noisy simulator, as expected.

In the previous examples, the MI is calculated based on the exact wavefunction
(WF) for proof-of-principle. But, the exact WF is not known before the VQE
calculation. However, an approximate calculation of MI provides a sufficient
estimate of the correlation between different qubits, as shown in the previous
DMRG-type selection of active space Stein and Reiher (2016) and our recent
development of the PermVQE algorithm Tkachenko  _et al._ (2021). The
approximate MI can be estimated on classical computers by using either MP2 or
DMRG (with small bond-dimension) methods. Here, we also demonstrate that
ClusterVQE even works with MI calculated on-the-fly. At each iteration, the MI
is updated and used for re-clustering the qubits. As shown by the red dots in
Figures 3 (a) and (c), the performance of ClusterVQE on-the-fly is almost the
same as that using the exact MI. The clustering is suboptimal only for the
first three iterations and the dressed Hamiltonian has a slightly larger size
as a result. The compositions of each cluster are changing before a stable
clustering is found. Nevertheless, these results illustrate that ClusterVQE
also works when clustering on-the-fly.

Figure 7:  Simulation of the LiH molecule (10 qubits) on the IBMQ-Rome quantum
device (5 qubits). The error rate of the CNOT gate is $0.0067\sim 0.026$ and
0.014 in average at runtime. The inset plot shows the VQE optimization for the
first, third and fifth ClusterVQE iterations.

###  III.4 Simulation on a quantum computer

Finally, we implement ClusterVQE on an actual quantum backend. The IBMQ-Rome
device with 5 qubits was chosen to simulate the LiH molecule described by 10
qubits (lowest energy spatial orbital is assumed to be frozen). In order to
run on IBMQ-Rome, we set the size of each cluster to be 5. Figure 7 shows the
energies as a function of ClusterVQE iterations without error mitigation, and
the inset plot shows the VQE optimization for the first, third, and fifth
ClusterVQE iterations, respectively. Due to the large CNOT error rate, the
energies go up with the first a few ClusterVQE iterations without error
mitigation. Namely, the number of entanglers encoded on the quantum circuits
increases with the number of iterations and introduces more errors. However,
the energy error does not increase monotonically with ClusterVQE iterations,
suggesting that the ClusterVQE can compete with noise even without error
mitigation.

The ClusterVQE performance on real quantum devices can be further improved
with recently developed error mitigation techniques Cincio  _et al._ (2021a).
For proof-of-principle, we first perform ClusterVQE on the ideal simulator to
generate the circuits of the final iteration and run the final circuits on
both the noisy simulator and real quantum device (IBMQ-Rome) for error
mitigation. We perform the mitigation using Clifford data regression (CDR)
Czarnik  _et al._ (2020), 111CDR training set is constructed using 4 near-
Clifford circuits with 1 non-Clifford gate in each of 2 clusters. The training
circuits are chosen by replacing non-Clifford gates with the closest Clifford
gates.. For the noisy simulation, we use a noise model obtained by gate set
tomography of the IBMQ Ourense device Cincio  _et al._ (2021b). The results
are listed in Table. 1, which shows that error mitigation reduces the error by
several orders of magnitude. Such simulations demonstrate that the ClusterVQE
method, with error mitigation, can readily achieve chemical accuracy on real
quantum devices.

Table 1:  Energy errors with and without error mitigation. We show results for different numbers of shots per measurement of a Hamiltonian term.  Backend | Shots | Energy error (au)
---|---|---
Not mitigated | Mitigated
Noisy simulator | $204080$ | 0.225 | $5.6\times 10^{-2}$
Noisy simulator | $\infty$ | 0.221 | $1.0\times 10^{-5}$
IBMQ-Rome | 5120 | 0.236 | $1.98\times 10^{-3}$

##  IV Conclusion

To summarize, we developed a novel ClusterVQE algorithm to reduce quantum
circuit complexity, including both the number of qubits and circuit depth. Our
algorithm groups qubits into clusters based on maximal correlation between
them and further takes into account the residual entanglement between clusters
in a new ”dressed” Hamiltonian. Hence, ClusterVQE naturally implements the
benefits of the previously developed PermVQE algorithm, which significantly
reduces the number of controlled-NOT gates to connect strongly correlated
qubits Tkachenko  _et al._ (2021). The entanglers that couple different
clusters are further used to dress the Hamiltonian and remove the entanglement
between different clusters. Consequently, the Hamiltonian can be decomposed as
a tensor product of each cluster. Each cluster can then be measured separately
with a smaller number of qubits and shallower circuit depth. Therefore,
ClusterVQE is a particularly promising algorithmic development for quantum
chemistry applications, where it is able to reduce the quantum hardware
requirements significantly. Moreover, our algorithm has the flexibility of
reducing qubits to any number at the expense of additional computational costs
on classical computers. Since fewer entanglers are encoded on quantum
circuits, ClusterVQE is more robust to noise than qubit-ADAPT-VQE, being one
of the most circuit-efficient VQE algorithms, which is verified by our
numerical experiments on both the noisy simulator and real quantum platforms.
Since the correlation between different clusters is minimized, only a few
operators are used to dress the Hamiltonian. Consequently, the Hamiltonian
size does not increase significantly compared to that of the iQCC method.

ClusterVQE also uses a newly proposed gradient measurement method for the VQE
without using an ancillary qubit. The analytical gradient can significantly
improve the performance of ClusterVQE in noisy environments. The simulation of
the LiH molecule on the noisy quantum simulator found that the energy error of
the ClusterVQE algorithm is smaller than that of the qubit-ADAPT-VQE method
due to the shallower circuit depth. A groundbreaking example of the exact
error-corrected LiH molecule simulation described by 10 qubits on the 5-qubits
IBMQ-Rome quantum computer with ClusterVQE is demonstrated. Obtained results
further demonstrate the validity and advantage of the ClusterVQE algorithm.
Hence, we believe the algorithm developed in this work can pave the way
towards the adaptation of molecular quantum chemistry simulation on NISQ
devices. It should be mentioned that since the size of the qubit Hamiltonian
increases with the number of operators used to dress the Hamiltonian, it
increases both computational costs on classical computers and measurements on
quantum computers. The recently developed measurement reduction techniques
Bonet-Monroig  _et al._ (2019); Yen and Izmaylov (2020); Huang and Kueng
(2019); Zhao  _et al._ (2020); Huang  _et al._ (2020); Yen  _et al._ (2020);
Izmaylov  _et al._ (2019); Huggins  _et al._ (2019) are likely a path-forward
to further reduce the simulation cost on quantum devices.

Finally, it should be noted that the computational scaling of ClusterVQE is
dependent on the number of dressing operators. Even though the scaling of
computing the dressed Hamiltonian is exponential against the number of
dressing operators Ilya G. Ryabinkin , the scaling of ClusterVQE against
iterations is not exponential in general because not every entangler is used
to dress the Hamiltonian. For the best cases, only a limited number of
entanglers is used to dress the Hamiltonian (LiH, for example). Such a
situation can be intuitively expected in molecular multimers, where the
intermolecular coupling is usually smaller than the intramolecular
correlation. For the worse cases, usually in strongly correlated systems in
which qubit clustering is not optimal, the scaling of ClusterVQE may be
exponential if a larger number of entanglers are used to dress the
Hamiltonian. However, even in strongly correlated molecules, the scaling of
different geometries can be different. For example, the ClusterVQE simulation
of 14-qubits $N_{2}$ at equilibrium ends up with a large number of Pauli words
(Figure 5(b)). However, the simulation of $N_{2}$ with stretched bond has a
much better scaling (Figure 5(d)) even though the stretched configuration has
a stronger correlation. That being said, the bottleneck of our approach relies
on whether the inter-cluster MI can be efficiently minimized. Hence,
ClusterVQE does not introduce exponential scaling for general cases.

###### Acknowledgements.

Research presented in this article was supported by the Laboratory Directed
Research and Development (LDRD) program of Los Alamos National Laboratory
(LANL) under project number 20200056DR. LANL is operated by Triad National
Security, LLC, for the National Nuclear Security Administration of U.S.
Department of Energy (contract no. 89233218CNA000001). We thank LANL
Institutional Computing (IC) program for access to HPC resources. PJC and LC
were partially supported by the U.S. Department of Energy (DOE), Office of
Science, Office of Advanced Scientific Computing Research, under the
Accelerated Research in Quantum Computing (ARQC) program. Quantum resources
from IBMQ Hub are acknowledged. This research was also supported by the U.S.
Department of Energy (DOE) National Nuclear Security Administration (NNSA)
Advanced Simulation and Computing (ASC) program at LANL. We acknowledge the
ASC program at LANL for use of their Ising D-Wave 2000Q quantum computing
resource.

Appendix

##  I Clustering Methods Used

Three clustering methods were explored for this effort. The details of each
are described below.

###  I.1 Metis Graph Partitioning on CPU

The MI matrix can be viewed as a graph. This allows for the use of classical
graph partitioning methods from the Metis library
https://metis.readthedocs.io/en/latest . The part$\\_$graph function in the
Metis library was used and the partition weights were also optimized to
minimize the cut.

###  I.2 MI-based clustering on Quantum Annealer QPU

The MI between qubits can be interpreted as a degree of correlation. Note that
when two random variables are fully correlated, their joint entropy is
minimal, and the MI is maximal. Moreover, MI can be mapped to what others have
coined as the “generalized correlation coefficient”, where all higher order
correlations are taking into account Lange and Grubmüller (2008, 2005);
Rivalta  _et al._ (2012). This generalized correlation coefficient is computed
as:

$r^{MI}_{ij}=[1-\mathrm{e}^{-2I_{ij}/d}]^{1/2}$ (S1)

where $d$ is the dimensionality of the problem, and $I_{ij}$ are is the Mutual
Information matrix elements. Note that, when $I_{ij}$ goes to infinity,
$r^{MI}_{ij}$ approaches 1; whereas if $I_{ij}$ tends to zero, $r^{MI}_{ij}$
tends to 0. Taking into account Equation S2, we see that, by maximizing MI for
a constrained set of qubits, we will be maximizing the correlation among that
same set. Moreover, maximizing the MI restricted to a certain number of
qubits, is somehow related to getting the “most central hub,” according to
node centrality ordering. Albeit for some exceptions, we have verified the
latter using eigenvector centrality metric.

Maximizing the MI restricted to a certain number of qubits can be formulated
as a quadratic unconstrained binary optimization (QUBO) problem. This
formulation can, in turn, benefit from available alternative computational
paradigms tailored to solve these kinds of problems, e.g. Quantum, and Digital
annealers. The cost function to be maximized can be written as:

$C=\textbf{x}(I-\lambda R)\textbf{x}^{t}$ (S2)

where the vector x is a bit-string containing the information of the selected
qubits. This is, if the i-th qubit is selected $x_{i}=1$, and 0 otherwise.
Matrix $I$ contains all the MI elements $I_{ij}$ as defined by Eq. 1, and $R$
is the matrix representation of an operator that “restricts” the total number
of selected qubits Rogers and Singleton (2020). A penalty constant $\lambda$
is used to tune the strength of the constraint. Since the annealing methods
and devices commonly minimize cost function, the problem is rather solved by
minimizing $-C$ leading to a QUBO matrix equal to $Q=-I+\lambda R$

Matrix R, is defined as follows.

$R_{ij}=\begin{cases}1&\quad\text{if}\quad i\neq j\\\ 1-\alpha\times f\times
S&\quad\text{if}\quad i=j\end{cases}$

Here, $f$ is the target fraction of qubits to be selected; $S$ is the total
size of the qubits pool; and $\alpha$ is a constant that depends on the
optimizer (1 or 2 depending if only the upper triangle vs the full QUBO matrix
is used to compute the energy). For all the calculations where this clustering
method was used we have set $f=0.5$, and $\lambda=2.0$.

An explanation of why the $R$ matrix restricts the set is described next.
Consider the case where a general vector $\mathbf{x}$ of length $S=N$ is such
that $\sum_{k}x_{k}=K$. Now let us compute the energy term of that vector
given by matrix $R$. This is
$E_{R}(\mathbf{x},\mathrm{such\,that\,}||x||_{1}=K)=\mathbf{x}^{t}R\mathbf{x}$.
Then, the off diagonal elements of $R$ will contribute with $K(K-1)$. If we
want $f=M/N$, then, the diagonal elements will contribute with $K-\alpha MK$.
By summing the two contributions, we have that $E_{R}(K)=K^{2}-\alpha
MK=K(K-\alpha M)$. This means that the function will have two minima, the
trivial one at K=0, and another at $K=\alpha M$ ($\alpha M$ bits equal to 1).

###  I.3 Quantum Community Detection on Quantum Annealer QPU

Classical community detection as graph clustering maximizes modularity to
split the MI matrix into clusters that minimize correlation or connectivity
between clusters Newman (2006). In quantum community detection (QCD), the
problem is represented as a QUBO that can be run on a D-Wave quantum annealer
as an NP-complete combinatorial optimization problem.

The symmetric MI matrix is viewed as a graph. The adjacency matrix $A$ is
calculated as follows.

$A_{ij}=\begin{cases}0&i=j\\\ |{MI}_{ij}|&i\neq j\end{cases}$

The modularity matrix $B$ is calculated from the adjacency matrix $A$ and
degree of each node (number of connected nodes) $d_{i}$.

$B_{ij}=A_{ij}-\frac{d_{i}d_{j}}{2m}=A_{ij}-\frac{d_{i}d_{j}}{\sum_{i}d_{i}}$
(S3)

While QCD can cluster a graph into 2 or more clusters Ushijima-Mwesigwa  _et
al._ (2017); Negre  _et al._ (2020); Mniszewski  _et al._ (2021), in this work
we are interested in splitting only into 2 clusters. The modularity matrix in
QUBO form is run on the D-Wave 2000Q quantum annealer:

$Q(x)=max(x^{T}Bx)$ (S4)

The resulting $x$ is a string of zeros and ones that maps each of the MI nodes
(representing qubits) to one of two clusters that can be used in the
ClusterVQE algorithm.

Figure S1:  The error in energy ($\Delta E$) as a function of the error rates
of single (a) and two-qubit (b) gates, respectively. The two and single-qubit
gate error in a and b are fixed at 0.01 and 0.001, respectively. The single-
qubit gate on real quantum devices is usually one order of magnitude smaller
than that of two-qubit gates. Hence, single qubit error has a smaller impact
on the results.

##  II Introduction to UCCSD ansatz and Qubit-UCCSD

Within the UCCSD ansatz, the correlation beyond the Hartree-Fock reference
state is introduced via the single and double excitation (and de-excitation)
operators:
$\ket{\Phi(\theta)}=e^{\hat{T}(\theta)-\hat{T}^{\dagger}(\theta)}\ket{\Phi_{r}}$,
where $\ket{\Phi_{r}}$ is the reference state being the HF state in this work.
$\hat{T}(\theta)=\hat{T}_{1}(\theta_{1})+\hat{T}_{2}(\theta_{2})$ is the
excitation operator with $\hat{T}_{1}(\theta_{1})$ and
$\hat{T}_{2}(\theta_{2})$ being the parameterized single and double excitation
operators, respectively,

$\displaystyle\hat{T}_{1}(\theta)=$
$\displaystyle\sum_{ia}\theta^{a}_{i}\hat{c}^{\dagger}_{a}\hat{c}_{i}$
$\displaystyle\hat{T}_{2}(\theta)=$
$\displaystyle\sum_{ij,ab}\theta^{ab}_{ij}\hat{c}^{\dagger}_{a}\hat{c}^{\dagger}_{b}\hat{c}_{i}\hat{c}_{j}.$
(S5)

For simplicity, indices $i,j,\cdots$ are used to represent the occupied
molecular orbitals (MOs), $a,b,\cdots$ represent the unoccupied MOs. Because
the NISQ devices only have gates acting on a few qubits at a time, the UCCSD
operator must be broken up into a sequence of single operators (either single
or double excitation) via the Trotter expansion of an exponent operator Suzuki
(1990), i.e.,
$e^{A+B}=\lim_{n\rightarrow\infty}\left(e^{A/n}e^{B/n}\right)^{n}$. Even
though the truncated Trotter expansion is an approximation to the UCCSD
ansatz, the variational procedure is able to absorb the error according to
recent report O’Malley  _et al._ (2016). Consequently, a finite order of
Trotter expansion is able to reproduce the FCI energies due to the fact that
the many-body interaction can be decomposed as the tensor product of one- and
two-body interactions Grimsley  _et al._ (2019). In this way, the exact FCI
state can be written as the tensor product of single and double excitation
operators,
$\ket{\Phi(\theta)}=\prod^{D}_{k=1}e^{\hat{t}_{k}(\theta)}\ket{\Phi_{r}}$,
where $\hat{t}_{k}(\theta)$ represents either single or double excitation
operators and their anti-Hermitian forms
($\hat{t}^{\dagger}_{k}=-\hat{t}_{k}$). Consequently, the GS energy is readily
obtained by minimizing the parameterized ansatz, i.e.,

$\displaystyle E_{GS}=$
$\displaystyle\min_{\theta}\bra{\Phi(\theta)}\hat{H}\ket{\Phi(\theta)}$
$\displaystyle=$
$\displaystyle\min_{\theta}\bra{\Phi_{r}}\prod^{1}_{k=D}e^{-\hat{t}_{k}(\theta_{k})}\hat{H}\prod^{D}_{k=1}e^{\hat{t}_{k}(\theta_{k})}\ket{\Phi_{r}}$
(S6)

In the UCCSD ansatz, the SD sequences are directly generated from the Trotter
expansion of the $\hat{T}(\theta)-\hat{T}^{\dagger}(\theta)$ operator. In the
recently developed ADAPT-VQE algorithm Grimsley  _et al._ (2019), the sequence
of SD operators in the ansatz are chosen adaptively by measuring their
contribution to the GS energy weighted by the energy gradients with respect to
the parameters. With the adaptive method, a more compact sequence of the
unitary operators is generated and a smaller number of parameters is optimized
at the cost of more VQE simulations. Because the UCCSD ansatz is not hardware-
efficient, a qubit version of the UCCSD ansatz was developed recently in which
the ansatz is broken up into a series of Pauli words and a small portion of
the Pauli words is chosen to adaptively grow the ansatz. Eq. II.2 indicates
the product $\hat{P}_{kj}=[\hat{P}_{k},\hat{P}_{j}]$, that only contains the
$\hat{Z}$ operators, has non-vanishing derivatives because
$\bra{\Phi_{r}}\hat{\sigma}_{i}\ket{\Phi_{r}}=0$ if
$\hat{\sigma}_{i}=\hat{X},\hat{Y}$. Therefore, the Pauli word $\hat{P}_{j}$
that contains $\hat{Z}$ has no contribution to the derivative. And two Pauli
words with the the same flip indices (indices of $\hat{X}/\hat{Y}$) have the
same contribution. Hence, after the fermion-qubit mapping of the UCCSD
operators, we only include those Pauli words with unique flip indices. Thus,
only one of these operators is included in the unitary operators
$\\{\hat{U}_{j}\\}$. In this way, the number of operators and circuit depth
can be reduced.

##  III Analytical gradient in VQE

The energy gradient is used in the parameter optimization on the classical
computers. The gradient can be obtained by the finite-difference method or
analytical gradients. According to our numerical experiment on either a noisy
simulator or real quantum devices, the analytical gradient can significantly
improve the performance of VQE in the presence of noise. In this work, the
analytical gradient is obtained by the measurements on quantum circuits. The
gradient of a parameter is,

$\frac{\partial
E}{\partial\theta_{i}}=2\text{Im}\left\\{\bra{\Phi_{r}}\hat{U}^{\dagger}\hat{H}\frac{\partial\hat{U}}{\partial\theta_{i}}\ket{\Phi_{r}}\right\\}.$
(S7)

Suppose there are $M$ operators in the ansatz, i.e,
$\hat{U}=\prod^{M}_{k=1}\hat{U}_{k}$, the derivative of $\hat{U}$ with respect
to $\theta_{i}$ reads

$\displaystyle\frac{\partial\hat{U}}{\partial\theta_{i}}=$ $\displaystyle
i\prod^{M}_{k=i+1}\hat{U}_{k}\hat{P}_{i}\prod^{i}_{k=1}\hat{U}_{k}$
$\displaystyle=$ $\displaystyle
i\left(\prod^{M}_{k=i+1}\hat{U}_{k}\hat{P}_{i}\prod^{i+1}_{k=M}\hat{U}^{\dagger}_{k}\right)\prod^{M}_{k=1}\hat{U}_{k}\equiv
i\hat{P}^{\text{eff}}_{i}\hat{U},$ (S8)

where $\hat{P}^{\text{eff}}_{i}$ is the dressed operator. Consequently, the
energy gradient can be obtained by

$\frac{\partial
E}{\partial\theta_{i}}=2\text{Im}\left\\{\bra{\Phi_{r}}\hat{U}^{\dagger}\hat{H}\hat{P}^{\text{eff}}_{i}\hat{U}\ket{\Phi_{r}}\right\\}.$
(S9)

In this way, calculation of the dressed operator $\hat{P}^{\text{eff}}_{i}$ on
classical computers is required, while no ancillary qubit on the quantum
circuit is needed.

## References

  * Gan and Harrison (2005) Z. Gan and R. J. Harrison, in  _SC ’05: Proceedings of the 2005 ACM/IEEE Conference on Supercomputing_ (2005) pp. 22–22.
  * Cao  _et al._ (2019) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, _et al._ , Chem. Rev. 119, 10856 (2019).
  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Rev. Mod. Phys. 92, 015003 (2020).
  * Kitaev (1995) A. Y. Kitaev, “Quantum measurements and the abelian stabilizer problem,” (1995), arXiv:quant-ph/9511026 [quant-ph] .
  * Aspuru-Guzik  _et al._ (2005) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704 (2005).
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Wei  _et al._ (2020) S. Wei, H. Li, and G. Long, Research 2020, 1486935 (2020).
  * Moll  _et al._ (2018) N. Moll, P. Barkoutsos, L. S. Bishop, J. M. Chow, A. Cross, D. J. Egger, S. Filipp, A. Fuhrer, J. M. Gambetta, M. Ganzhorn, A. Kandala, A. Mezzacapo, P. Müller, W. Riess, G. Salis, J. Smolin, I. Tavernelli, and K. Temme, Quantum Sci. Technol. 3, 030503 (2018).
  * Cerezo  _et al._ (2020) M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, “Variational quantum algorithms,” (2020), arXiv:2012.09265 [quant-ph] .
  * Bharti  _et al._ (2021) K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, “Noisy intermediate-scale quantum (nisq) algorithms,” (2021), arXiv:2101.08448 [quant-ph] .
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M. H. Yung, X. Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nat. Commun. 5, 4213 (2014).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New J. Phys. 18, 023023 (2016).
  * Li  _et al._ (2019a) Y. Li, J. Hu, X.-M. Zhang, Z. Song, and M.-H. Yung, Adv. Theory Simul. 2, 1800182 (2019a).
  * Fedorov  _et al._ (2021) D. A. Fedorov, B. Peng, N. Govind, and Y. Alexeev, “Vqe method: A short survey and recent developments,” (2021), arXiv:2103.08505 [quant-ph] .
  * Shen  _et al._ (2017) Y. Shen, X. Zhang, S. Zhang, J.-N. Zhang, M.-H. Yung, and K. Kim, Phys. Rev. A 95, 020501 (2017).
  * McClean  _et al._ (2017) J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong, Phys. Rev. A 95, 042308 (2017).
  * Nakanishi  _et al._ (2019) K. M. Nakanishi, K. Mitarai, and K. Fujii, Phys. Rev. Research 1, 033062 (2019).
  * Higgott  _et al._ (2019) O. Higgott, D. Wang, and S. Brierley, Quantum 3, 156 (2019).
  * Kawai and Nakagawa (2020) H. Kawai and Y. O. Nakagawa, Machine Learning: Science and Technology 1, 045027 (2020).
  * Greene-Diniz and Muñoz Ramo (2021) G. Greene-Diniz and D. Muñoz Ramo, Int. J. Quantum Chem. 121, e26352 (2021).
  * Kandala  _et al._ (2019) A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Nature 567, 491 (2019).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Tang  _et al._ (2021a) H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, PRX Quantum 2, 020310 (2021a).
  * McClean  _et al._ (2018) J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Nat. Commun. 9, 4812 (2018).
  * Wang  _et al._ (2021) S. Wang, E. Fontana, M. Cerezo, K. Sharma, A. Sone, L. Cincio, and P. J. Coles, “Noise-induced barren plateaus in variational quantum algorithms,” (2021), arXiv:2007.14384 [quant-ph] .
  * Cerezo  _et al._ (2021) M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Nat. Commun. 12, 1791 (2021).
  * Abbas  _et al._ (2020) A. Abbas, D. Sutter, C. Zoufal, A. Lucchi, A. Figalli, and S. Woerner, “The power of quantum neural networks,” (2020), arXiv:2011.00027 [quant-ph] .
  * Bartlett  _et al._ (1989) R. J. Bartlett, S. A. Kucharski, and J. Noga, Chem. Phys. Lett. 155, 133 (1989).
  * Bartlett and Musiał (2007) R. J. Bartlett and M. Musiał, Rev. Mod. Phys. 79, 291 (2007).
  * Seeley  _et al._ (2012a) J. T. Seeley, M. J. Richard, and P. J. Love, J. Chem. Phys. 137, 224109 (2012a).
  * Barkoutsos  _et al._ (2018) P. K. Barkoutsos, J. F. Gonthier, I. Sokolov, N. Moll, G. Salis, A. Fuhrer, M. Ganzhorn, D. J. Egger, M. Troyer, A. Mezzacapo, S. Filipp, and I. Tavernelli, Phys. Rev. A 98, 022322 (2018).
  * Romero  _et al._ (2018) J. Romero, R. Babbush, J. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, Quantum Sci. Technol. 4, 14008 (2018).
  * Kühn  _et al._ (2019) M. Kühn, S. Zanker, P. Deglmann, M. Marthaler, and H. Weiß, J. Chem. Theory Comput. 15, 4764 (2019), pMID: 31403781\.
  * Ryabinkin  _et al._ (2018) I. G. Ryabinkin, T.-C. Yen, S. N. Genin, and A. F. Izmaylov, J. Chem. Theory Comput. 14, 6317 (2018).
  * Lee  _et al._ (2019) J. Lee, W. J. Huggins, M. Head-Gordon, and K. B. Whaley, J. Chem. Theory Comput. 15, 311 (2019).
  * Matsuzawa and Kurashige (2020) Y. Matsuzawa and Y. Kurashige, J. Chem. Theory Comput. 16, 944 (2020).
  * Xia and Kais (2020) R. Xia and S. Kais, Quantum Sci. Tech. 6, 015001 (2020).
  * Tkachenko  _et al._ (2021) N. V. Tkachenko, J. Sud, Y. Zhang, S. Tretiak, P. M. Anisimov, A. T. Arrasmith, P. J. Coles, L. Cincio, and P. A. Dub, PRX Quantum 2, 020337 (2021).
  * Grimsley  _et al._ (2019) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nat. Commun. 10, 3007 (2019).
  * Tang  _et al._ (2021b) H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, PRX Quantum 2, 020310 (2021b).
  * Ryabinkin  _et al._ (2020) I. G. Ryabinkin, R. A. Lang, S. N. Genin, and A. F. Izmaylov, J. Chem. Theory Comput. 16, 1055 (2020).
  * Lang  _et al._ (2021) R. A. Lang, I. G. Ryabinkin, and A. F. Izmaylov, J. Chem. Theory Comput. 17, 66 (2021).
  * Ryabinkin  _et al._ (2021) I. G. Ryabinkin, A. F. Izmaylov, and S. N. Genin, Quantum Sci. Tech. 6, 024012 (2021).
  * Yordanov  _et al._ (2020) Y. S. Yordanov, V. Armaos, C. H. W. Barnes, and D. R. M. Arvidsson-Shukur, “Iterative qubit-excitation based variational quantum eigensolver,” (2020), arXiv:2011.10540 [quant-ph] .
  * Bravyi  _et al._ (2017) S. Bravyi, J. M. Gambetta, A. Mezzacapo, and K. Temme, “Tapering off qubits to simulate fermionic hamiltonians,” (2017), arXiv:1701.08213 [quant-ph] .
  * Setia  _et al._ (2020) K. Setia, R. Chen, J. E. Rice, A. Mezzacapo, M. Pistoia, and J. D. Whitfield, J. Chem. Theory Comput. 16, 6091 (2020).
  * Yen  _et al._ (2019) T.-C. Yen, R. A. Lang, and A. F. Izmaylov, J. Chem. Phys. 151, 164111 (2019), https://doi.org/10.1063/1.5110682 .
  * Zhang  _et al._ (2021a) F. Zhang, N. Gomes, N. F. Berthusen, P. P. Orth, C.-Z. Wang, K.-M. Ho, and Y.-X. Yao, Phys. Rev. Research 3, 013039 (2021a).
  * Rissler  _et al._ (2006) J. Rissler, R. M. Noack, and S. R. White, Chem. Phys 323, 519 (2006).
  * Huang and Kais (2005) Z. Huang and S. Kais, Chem. Phys. Lett. 413, 1 (2005).
  * Zhang  _et al._ (2021b) Z.-J. Zhang, T. H. Kyaw, J. Kottmann, M. Degroote, and A. Aspuru-Guzik, Quantum Sci.Tech. (2021b).
  * (52) https://metis.readthedocs.io/en/latest, “Metis for python,” .
  * Mniszewski  _et al._ (2021) S. M. Mniszewski, P. A. Dub, S. Tretiak, P. M. Anisimov, Y. Zhang, and C. F. A. Negre, Sci. Rep. 11, 4099 (2021).
  * Booth  _et al._ (2017) M. Booth, S. P. Reinhardt, and A. Roy,  _Partitioning Optimization Problems for Hybrid Classical / Quantum Execution_ , Tech. Rep. (2017).
  * Jordan and Wigner (1928) P. Jordan and E. Wigner, Zeitschrift für Physik 47, 631 (1928).
  * Bravyi and Kitaev (2002) S. B. Bravyi and A. Y. Kitaev, Ann. Phys. 298, 210 (2002).
  * Seeley  _et al._ (2012b) J. T. Seeley, M. J. Richard, and P. J. Love, J. Chem. Phys. 137, 224109 (2012b).
  * Li  _et al._ (2019b) Y. Li, J. Hu, X.-M. Zhang, Z. Song, and M.-H. Yung, Adv. Theory and Simul. 2, 1800182 (2019b).
  * Dallaire-Demers  _et al._ (2018) P.-L. Dallaire-Demers, J. Romero, L. Veis, S. Sim, and A. Aspuru-Guzik, arXiv:1801.01053 (2018).
  * Bonet-Monroig  _et al._ (2019) X. Bonet-Monroig, R. Babbush, and T. E. O’Brien, arXiv:1908.05628 (2019).
  * Yen and Izmaylov (2020) T.-C. Yen and A. F. Izmaylov, arXiv:2007.01234 (2020).
  * Huang and Kueng (2019) R.-Y. Huang and i. Kueng, arXiv:1908.08909 (2019).
  * Zhao  _et al._ (2020) A. Zhao, A. Tranter, W. M. Kirby, S. F. Ung, A. Miyake, and P. J. Love, Phys. Rev. A 101, 062322 (2020).
  * Huang  _et al._ (2020) H.-Y. Huang, R. Kueng, and J. Preskill, arXiv:2002.08953 (2020).
  * Yen  _et al._ (2020) T.-C. Yen, V. Verteletskyi, and A. F. Izmaylov, J. Chem. Theory Comput. 16, 2400 (2020).
  * Izmaylov  _et al._ (2019) A. F. Izmaylov, T.-C. Yen, R. A. Lang, and V. Verteletskyi, J. Chem. Theory Comput. (2019).
  * Huggins  _et al._ (2019) W. J. Huggins, J. McClean, N. Rubin, Z. Jiang, N. Wiebe, K. B. Whaley, and R. Babbush, arXiv:1907.13117 (2019).
  * Abraham and et al (2019) H. Abraham and et al, “Qiskit: An open-source framework for quantum computing,”  (2019).
  * Morales and Nocedal (2011) J. Morales and J. Nocedal, ACM Trans. Math. Softw. 38 (2011).
  * (70) https://quantum-computing.ibm.com .
  * Stein and Reiher (2016) C. J. Stein and M. Reiher, J. Chem. Theory Comput. 12, 1760 (2016).
  * Cincio  _et al._ (2021a) L. Cincio, K. Rudinger, M. Sarovar, and P. J. Coles, PRX Quantum 2, 010324 (2021a).
  * Czarnik  _et al._ (2020) P. Czarnik, A. Arrasmith, P. J. Coles, and L. Cincio, arXiv:2005.10189 (2020).
  * Note (1) CDR training set is constructed using 4 near-Clifford circuits with 1 non-Clifford gate in each of 2 clusters. The training circuits are chosen by replacing non-Clifford gates with the closest Clifford gates.
  * Cincio  _et al._ (2021b) L. Cincio, K. Rudinger, M. Sarovar, and P. J. Coles, PRX Quantum 2, 010324 (2021b).
  * (76) S. N. G. Ilya G. Ryabinkin, Artur F. Izmaylov, .
  * Lange and Grubmüller (2008) O. F. Lange and H. Grubmüller, Proteins 70, 1294 (2008).
  * Lange and Grubmüller (2005) O. F. Lange and H. Grubmüller, “Generalized correlation for biomolecular dynamics,” (2005).
  * Rivalta  _et al._ (2012) I. Rivalta, M. M. Sultan, N.-S. Lee, G. A. Manley, J. P. Loria, and V. S. Batista, PNAS 109, E1428 (2012).
  * Rogers and Singleton (2020) M. L. Rogers and R. L. Singleton, “Floating-Point calculations on a quantum annealer: Division and matrix inversion,” (2020).
  * Newman (2006) M. E. J. Newman, PNAS 103, 8577 (2006).
  * Ushijima-Mwesigwa  _et al._ (2017) H. Ushijima-Mwesigwa, C. F. A. Negre, and S. M. Mniszewski, in  _Proceedings of the Second International Workshop on Post Moores Era Supercomputing (PMES’17)_ (ACM New York, NY, USA, 2017) pp. 22–29.
  * Negre  _et al._ (2020) C. F. A. Negre, H. Ushijima-Mwesigwa, and S. M. Mniszewski, PLOS ONE 15, e0227538 (2020).
  * Suzuki (1990) M. Suzuki, Phys. Lett. A 146, 319 (1990).
  * O’Malley  _et al._ (2016) P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis, Phys. Rev. X 6, 031007 (2016).


[1]organization=Centro de Laseres Pulsados, addressline=C/ Adaja 8,
postcode=37185, city=Villamayor, country=Spain

# Ablation Holes in Tape Targets Induced by Ultra-Intense Laser Pulses

Michael Ehret<EMAIL_ADDRESS>Jon Imanol Apiñaniz  Jose Luis Henares  Roberto
Lera
Diego de Luis  José Antonio Pérez-Hernández  Luca Volpe  Giancarlo Gatti

###### Abstract

We present a theoretical modelling able to predict the dimensions of
$\mathrm{mm}$-sized through-holes produced in interactions of the
$1\text{\,}\mathrm{PW}$ high-power $30\text{\,}\mathrm{fs}$ Ti:Sa laser VEGA-3
with tape targets. We find that sizes of through-holes can be calculated by
assuming that the full heat transferred from the laser-heated electron
population to the target by electron-electron collisions drives the
evaporation of target material. We demonstrate the good reproduction of
experimental results.

###### keywords:

laser-plasma accelerator , hot electron transport , target ablation , high-
repetition-rate targetry , ultra-intense laser , laser-driven evaporation

##  1 Introduction

Advances in laser technology [1, 2, 3, 4] make accessible the relativistic
interaction regime at PW-class Ti:Sa laser systems opening up to laser driven
sources of charged particle beams [5] and X-rays [6]. Bright ion beams can be
generated from solid-density targets by well known mechanisms in the vast
field of laser-plasma acceleration [7] such as Target Normal Sheath
Acceleration (TNSA) and Radiation Pressure Acceleration. Many applications in
science and technology can benefit from laser-driven ion beams such as hadron-
therapy [8, 9], isotope production [10], positron emission tomography [11],
ion beam microscopy [12], particle induced X-ray emission [13] as well as
inertial confinement fusion [14].

A major requirement for several of these applications is a high repetition
rate (HRR) delivery $>1\text{\,}\mathrm{Hz}$ of the laser-driven source.
Difficulty arises as well studied solid density targets are destroyed by the
interaction with ultra-intense high-power laser pulses. The frequent
replacement of the target is a hot topic in HRR targetry, e.g. by generation
of cryogenic ribbons [15], use of liquid jets [16] or by unwinding of a thin
tape under the interaction position [17, 18].

This study regards the characterization of laser-induced ablation of material
from solid tape targets – for cases in which the extent is so severe that
$\mathrm{mm}$-sized through-holes are drilled into the target. So far and to
our best knowledge, dimensions of tape targets are deduced in heuristic manner
by trial and error. We present this modelling of the creation of ablation
holes in large solid targets to enable a-priori studies on tape systems.

##  2 Material and methods

Experiments for this work are conducted in the VEGA-3 laser facility at the
Centro de Laseres Pulsados (CLPU). Here the high power laser pulse is focused
via a $f/10$ off-axis parabola to spots of
$d_{\mathrm{L}}~{}=~{}$12\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ full-width
at half-maximum (FWHM). The pulse is amplified to an energy $E_{\mathrm{L}}$
up to $30\text{\,}\mathrm{J}$, and compressed to a duration
$\tau_{\mathrm{L}}$ down to $28\text{\,}\mathrm{fs}$. Only laser pulse energy
and duration are varied for this study. After compression, the short laser
pulse is transported in high vacuum of
$1\text{\times}{10}^{-6}\text{\,}\mathrm{mbar}$ and the laser focal spot is
aligned onto a target for the generation of secondary sources. Routinely
reached intensities $I_{\mathrm{L}}=E_{\mathrm{L}}/(0.25\cdot\pi
d_{\mathrm{L}}^{2})$ range from
$1\text{\times}{10}^{18}\text{\,}\mathrm{W}\text{\,}{\mathrm{cm}}^{-2}2\text{\times}{10}^{21}\text{\,}\mathrm{W}\text{\,}{\mathrm{cm}}^{-2}$.
The energy on target is extrapolated from throughput calibrations recorded at
low-energy, the focal spot at high energy is estimated to be the same as for
low-energy measurements.

This study uses a tape target, which is a
$9\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick solid density Al stripe of
$10\text{\,}\mathrm{mm}$ width. The tape normal is vertically inclined by
$12.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ with respect to the laser axis. A
slider frame with $3\text{\,}\mathrm{mm}$ aperture ensures the tape to stay in
the laser-focus plane.

In order to understand the laser-driven heating of the target, numerical
simulations are performed for the evolution of the hot electron distribution
generated by laser. We follow a detailed model of target charging in short
laser pulse interactions [19] that intends to predict the expected discharge
current due to laser-heated relativistic electrons. The numerical code
ChoCoLaT2 [20] simulates the distribution of electrons generated by laser
heating on a thin disk target. The successive electron escape is mitigated by
the target potential, based on the driver laser parameters and the interaction
geometry. The model takes into account the collisional cooling of electrons
within cold solid density targets. The energy and time depending hot electron
distribution function $f(E,t)$ describes electrons inside the target and
evolves according to

$\displaystyle\partial_{t}f(E,t)$
$\displaystyle=\frac{h_{\mathrm{L}}(E)\Theta(\tau_{\mathrm{L}}-t)}{\tau_{\mathrm{L}}}-\frac{f(E,t)}{\tau_{\mathrm{ee}}(E)}-g(E,t)$
(1) $\displaystyle h_{\mathrm{L}}(E)$
$\displaystyle\overset{!}{=}\frac{N_{0}}{T_{0}}\exp{\left[-E/T_{0}\right]}$
(2) $\displaystyle N_{0}$ $\displaystyle\overset{!}{=}\int f(E,0)\text{~{}d}E$
(3)

where $h_{\mathrm{L}}(E)$ is a constant exponential source of hot electrons,
$\Theta(t)$ the Heaviside function limiting electron heating to the laser
duration, $\tau_{\mathrm{ee}}(E)$ the energy dependent cooling time and
$g(E,t)$ the rate of electron ejection from the target. The initial hot
electron temperature $T_{0}$ depends on laser wavelength and pulse intensity
[23, 22, 21]; and $N_{0}$ is re-normalized to the energy balance
$N_{0}T_{0}=\eta E_{\mathrm{L}}$ between the total energy of hot electrons in
the target and the absorbed laser energy. Simulations require the conversion
efficiency $\eta$ of laser energy to energy in the hot electron distribution:
for this work a module is added to the code which calculates them
automatically according to [24, 25]

$\eta=1.2\times
10^{15}\cdot\left({I_{\mathrm{L}}\left[{$\mathrm{W}\text{\,}{\mathrm{cm}}^{-2}$}\right]}\right)^{0.75}\quad,$
(4)

a scaling which is successfully applied to the interpretation of experimental
data in the intensity regime from
$1\text{\times}{10}^{18}\text{\,}\mathrm{W}\text{\,}{\mathrm{cm}}^{-2}1\text{\times}{10}^{21}\text{\,}\mathrm{W}\text{\,}{\mathrm{cm}}^{-2}$
[26, 27].

The hot electron cooling time depends on target material properties such as
mass density $\rho_{\mathrm{t}}$, mass number $A_{\mathrm{t}}$, atomic number
$Z_{\mathrm{t}}$, and the hot electron energy distribution that allows to
calculate average speed $\left<v\right>_{\mathrm{e}}$ and energy
$\left<E\right>_{\mathrm{e}}$. Its meticulous calculation is demonstrated in
[28] with an emphasis on cases relevant for this work. For studies
hereinafter, Fe and Ti are added to the library of ChoCoLaT2.

With further slight modifications to the source code, we extract the branching
ratio between energy that leaves the target (transported by escaping
electrons) and energy that remains in the target (by collisional damping). A
section is added to integrate the the second term of Eq. (1) over the full
energy distribution. The cumulative sum over all previous evaluation steps
$\Delta E_{\mathrm{ee}}(t)$ is then added to the output generated by the
simulation for every recorded time-step.

##  3 Theory and calculation

The energy of the focused laser pulse is partially absorbed by target
electrons, which in turn transfer a fraction of their energy to the target
bulk material in collisions with cold bulk electrons. If the energy needed to
evaporate the volume under the heated area is smaller than the energy
transferred, a through-hole can be produced.

The deposited heat is assumed equal to the energy deposited by laser-heated
electrons in collisions $\Delta E_{\mathrm{ee}}$, with efficiency

$\zeta=\frac{\Delta E_{\mathrm{ee}}}{\eta E_{\mathrm{L}}}\quad.$ (5)

The $\mathrm{ps}$-scale duration of the energy transfer from laser to bulk
material is much shorter than the timescale of heat diffusion in metals.
Therefore we assume the heating process to be instantaneous. Further it is
assumed that all the deposited energy is consumed for the evaporation of a
cylindrical volume of target material which includes heating, melting, and
further heating processes. For a cylindrical volume with radius
$R_{\mathrm{t}}$ and height $h_{\mathrm{t}}$, the energy balance equation
reads

$\zeta\eta E_{\mathrm{L}}=\left(h_{\mathrm{t}}\cdot\pi
R_{\mathrm{t}}^{2}\cdot\rho_{\mathrm{t}}\right)\cdot\\\
\left(\int_{T_{0}}^{T_{\mathrm{f}}}C_{\mathrm{s}}\mathrm{\,d}T+S_{\mathrm{t}}+\int_{T_{\mathrm{f}}}^{T_{\mathrm{b}}}C_{\mathrm{s}}\mathrm{\,d}T+E_{\mathrm{t}}\right)\quad,$
(6)

where $\rho_{\mathrm{t}}$ is the mass density of the tape material.
$C_{\mathrm{s}}$ denotes the specific heat of the target material,
$S_{\mathrm{t}}$ is the latent heat of fusion at melting temperature
$T_{\mathrm{f}}$ and $E_{\mathrm{t}}$ is the enthalpy of evaporation at the
boiling temperature $T_{\mathrm{b}}$.

##  4 Results

Tape sections which are shot by laser show clear circular holes opened up by
evaporation of target material and rims of ablation craters which can be
produced by melting and re-condensation, see fig 1. The tape drive moves at
speeds of the order of $\mathrm{cm}\text{\,}{\mathrm{s}}^{-1}$ and does not
influence the shapes of holes, there is no systematic ellipticity.

Figure 1: Tape target with ablation holes (allow view to black background) and
rims of craters (visible through brighter reflectivity of ambient light).
Nominal pulse energies range from $4\text{\,}\mathrm{J}9\text{\,}\mathrm{J}$,
of which $87\text{\,}\mathrm{\char 37\relax}$ arrive on-target.

The hole diameters scales with the on-target energy of laser pulses, as shown
in Fig. 2. Compared are shots of different pulse duration but the same focal
spot size. With only $\mathrm{\SIUnitSymbolMicro m}$ thick targets, material
from large $\mathrm{mm}$-sized areas is evaporated. The variation of results
at different pulse duration is below the standard deviation of hole diameters.
The later is deduced for each shot from the measurement uncertainty and the
difference of the exact hole shapes from perfect circles. Shallow ablation
craters are systematically larger than holes, but their sizes approach each
other towards higher on-target energy.

Figure 2: The diameter of the evaporation hole under variation of the laser
pulse energy at best focus.

Our modelling with Eq. (6) allows to calculate the diameters of holes by
solving the equation for $R_{\mathrm{t}}$. As all shots are conducted at best
focus, we compare experimental results and theoretical predictions in terms of
laser pulse power versus hole diameters in Fig. 3. The theoretical prediction
is in the margin of uncertainty for $37.5\text{\,}\mathrm{\char 37\relax}$ of
the data. With shots at a moderate power of $\approx 150\text{\,}\mathrm{TW}$
each shot evaporates $\approx 0.4\text{\,}\mathrm{mg}$ of Al. HRR operation at
$10\text{\,}\mathrm{Hz}$ would cause debris of $1\text{\,}\mathrm{kg}$ of Al
within $72\text{\,}\mathrm{h}$.

Figure 3: The diameter of the evaporation hole under variation of the laser
pulse power at best focus, compared are experimental results (EXP) and
theoretical predictions (THEO).

##  5 Discussion

Our model allows to deduce theoretical predictions for other laser
characteristics, materials and target dimensions. A set of common target
materials for laser-accelerators based on TNSA is compared in Fig. 4. Note
that the deposition of energy depends strongly on the laser absorption into
hot electrons. The absorption modelling that we apply has a discontinuity at
intensities of
$3.1\text{\times}{10}^{19}\text{\,}\mathrm{W}\text{\,}{\mathrm{cm}}^{-2}$ [24,
25], remaining as discontinuity in our model. This is the reason for a small
kink in the graphs.

Figure 4: The diameter of the evaporation hole under variation of the laser
pulse energy at different FWHM focal diameters $d_{\mathrm{L}}$ with
$30\text{\,}\mathrm{fs}$ pulses, and the target characteristics for
$10\text{\,}\mathrm{mm}$ wide tapes of differing thickness.

Thinner tapes deliver less areal density to the interaction region and suffer
from larger hole diameters. A $6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Al
tape shows much larger holes compared to a
$9\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Al tape. Despite the larger mass
density yielding a higher energy absorption, Fe tapes show smaller hole
diameters when compared to Al tapes of similar characteristics. This is due to
the larger melting and evaporation temperatures respectively.

Hole diameters can be reduced by choosing a tape material with high melting
and evaporation temperature as well with high latent heats, e.g. W. A
$6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ W tape shows similar holes to a
$9\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Al tape. The disadvantage of W with
respect to Al is its high mass density, leading to high rates of hot-electron
cooling by collisions, thus a lower ejection of relativistic electrons and
therefore to a source of lower ion energies at lower number density.

The minimum speed with which the tape should be unwound calculates to
$v_{\mathrm{tape}}=m_{0}\cdot 2R_{\mathrm{t}}$ where $m_{0}$ denotes a factor
which takes into account the mechanical stability of the punctured tape. It is
left to future studies to evaluate this factor in detail for different tape
materials. Similar relations will be found for the tape width
$w_{\mathrm{tape}}=m_{1}\cdot 2R_{\mathrm{t}}$ and the stabilising aperture
which guides the tape $a_{\mathrm{tape}}=m_{2}\cdot 2R_{\mathrm{t}}$.

##  6 Conclusions

Our modelling allows to study future tape target materials based on their
design parameters without the need for costly experimental campaigns. New
materials should have high melting and evaporation temperatures, heat of
fusion and evaporation enthalpy respectively. In parallel, a small electron
electron collision frequency is beneficial for a small heat supply to the
target material. An important parameter for the later is a small mass density
of the material.

For stability of the tape, a priori (a) the unwinding speed can now be
adjusted in order to be larger than expected hole diameters, (b) the tape
width can be chosen much larger than the hole diameter to avoid rupture, and
(c) the guiding aperture can be designed slightly larger than holes in order
to avoid welding of the tape to the mounting – thus rupture. Precisely
estimated targets make it possible to attach secondary functional structures,
e.g. grounding pins.

Further research effort is required but this model can help to understand the
laser induced matter ejection mechanisms in thin foils and high intensities.
For the present data sets, the assumption that all the heat stored in hot
electron population is spent in evaporation in thermal equilibrium is leading
to good agreement with experimental results. Other out-of-equilibrium matter
ejection mechanisms like phase explosion [29] are not taken into account and
the agreement is acceptable. New studies with materials of very different
thermodynamic parameters (specific heat, latent heat of evaporation, heat
conductivity, etc..) will help to further define the model validity.

## Acknowledgements

This work would not have been possible without the help of the laser- and the
engineering team at CLPU. This work received funding from the European Union’s
Horizon 2020 research and innovation program through the European IMPULSE
project under grant agreement No 871161 and from LASERLAB-EUROPE V under grant
agreement No 871124. It benefited from funding from the Ministerio de Ciencia,
Innovación y Universidades in Spain through ICTS Equipment grant No
EQC2018-005230-P, further from grant PID2021-125389OA-I00 funded by MCIN / AEI
/ 10.13039/501100011033 / FEDER, UE and by “ERDF A way of making Europe”, by
the “European Union” and in addition from grants of the Junta de Castilla y
León with No CLP263P20 and No CLP087U16.

.

## Data statement

The raw data and numerical methods that support the findings of this study are
available from the corresponding author upon reasonable request.

## References

  * [1] T. H. Maiman et al., Stimulated Optical Radiation in Ruby, Nature 187 (1966), 493-494. https://doi.org/10.1038/187493a0
  * [2] M. DiDomenico et al., Generation Of Ultrashort Optical Pulses by Mode Locking the YAIG:Nd Laser, Applied Physics Letters 8, 7 (1966), 180-183. https://doi.org/10.1063/1.1754544
  * [3] P. Maine et al., Generation of ultrahigh peak power pulses by chirped pulse amplification, IEEE Journal of Quantum Electronics 24, 2 (1988) 398-403. https://doi.org/10.1109/3.137
  * [4] M. Aoyama et al., 0.85 PW, 33 fs Ti:Sa laser, Optics Letters 28, 17 (2003), 1594-1596. https://doi.org/10.1364/OL.28.001594
  * [5] T. Tajima and V. Malka, Laser plasma accelerators, Plasma Physics and Controlled Fusion 62, 034004 (2012). https://doi.org/10.1088/1361-6587/ab6da4
  * [6] J. A. Chakera et al., A continuous kilohertz Cu Ka source produced by submillijoule femtosecond laser pulses for phase contrast imaging, Applied Physics Letters 93, 261501 (2008). https://doi.org/10.1063/1.3046727
  * [7] M. Borghesi, Ion Acceleration: TNSA and Beyond, Springer Proceedings in Physics 231 (2019). https://doi.org/10.1364/10.1007/978-3-030-25850-4_7
  * [8] M. Dosanjh et al.,Status of hadron therapy in Europe and the role of ENLIGHT, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 571, 1 (2007) 191-194. https://doi.org/10.1016/j.nima.2006.10.060
  * [9] S. V. Bulanov et al., Laser ion acceleration for hadron therapy, Physics-Uspekhi 57, 12 (2014) 1149-1179. https://doi.org/10.3367/ufne.0184.201412a.1265
  * [10] K. Nemoto et al., Laser-triggered ion acceleration and table top isotope production, Applied Physics Letters 78, 5 (2001) 595-597. https://doi.org/10.1063/1.1343845
  * [11] M. I. K. Santala et al., Production of radioactive nuclides by energetic protons generated from intense laser-plasma interactions, Applied Physics Letters 78, 1 (2001) 19-21. https://doi.org/10.1063/1.1335849
  * [12] F. E. Merrill et al., Proton Microscopy at FAIR, AIP Conference Proceedings 1195, 1 (2009) 667-670. https://doi.org/10.1063/1.3295228
  * [13] F. Mirani et al., Integrated quantitative PIXE analysis and EDX spectroscopy using a laser-driven particle source, Science Advances 7, 3 (2021). https://doi.org/10.1126/sciadv.abc8660
  * [14] M. Roth et al., Fast Ignition by Intense Laser-Accelerated Proton Beams, Physical Review Letters 86, 3 (2001) 436-439. https://doi.org/10.1103/PhysRevLett.86.436
  * [15] S. Garcia et al., Continuous production of a thin ribbon of solid hydrogen, Laser and Particle Beams 32, 4 (2014), 569-575. https://doi.org/10.1017/S0263034614000524
  * [16] J. T. Morrison et al., New Journal of Physics 20, 022001 (2018). https://doi.org/10.1088/1367-2630/aaa8d1
  * [17] T. Nayuki et al., Review of Scientific Instruments 74, 3293 (2003). https://doi.org/10.1063/1.1578156
  * [18] M. Noaman-ul-Haq et al., Physical Review Accelerators and Beams 20, 041301 (2017). https://doi.org/10.1103/PhysRevAccelBeams.20.041301
  * [19] A. Poyé et al., Thin target charging in short laser pulse interactions, Physical Review E 98, 033201 (2018). https://doi.org/10.1103/PhysRevE.98.033201
  * [20] A. Poyé et al., ChoCoLaT, CELIA Program Library (2015). https://researchgate.net/publication/284609502_ChoCoLaT
  * [21] S. C. Wilks et al., Absorption of ultra-intense laser pulses, Physical Review Letters 69, 9 (1992) 1383-1386. https://doi.org/10.1103/PhysRevLett.69.1383
  * [22] F. N. Beg et al., A study of picosecond laser-solid interactions up to e19Wcm-2, Physics of Plasmas 4, 2 (1997) 447-457. https://doi.org/10.1063/1.872103
  * [23] R. Fabbro et al., Planar laser-driven ablation: Effect of inhibited electron thermal conduction, The Physics of Fluids 28, 5 (1985) 1463-1481. https://doi.org/10.1063/1.864982
  * [24] J. Yu et al., Physics of Plasmas 6, 1318 (1999).
  * [25] M. H. Key et al., Physics of Plasmas 5, 1966 (1998).
  * [26] J. Schreiber et al., Analytical Model for Ion Acceleration by High-Intensity Laser Pulses, Physical Review Letters 97, 045005 (2006). http://dx.doi.org/10.1103/PhysRevLett.97.045005
  * [27] N. P. Dover et al., Effect of Small Focus on Electron Heating and Proton Acceleration in Ultrarelativistic Laser-Solid Interactions, Physical Review Letters 124, 084802 (2020). https://doi.org/10.1103/PhysRevLett.124.084802
  * [28] A. Poyé et al., Dynamic model of target charging by short laser pulse interactions, Physical Review E 92, 043107 (2015). http://dx.doi.org/10.1103/PhysRevE.92.043107
  * [29] N. Bulgakova et al., Pulsed laser ablation of solids: transition from normal vaporization to phase explosion, Applied Physics Letters A 73 (2001). https://doi.org/10.1007/s003390000686