# Confining continuous manipulations of accelerator beamline optics

Ph. Amstutz<EMAIL_ADDRESS>C. Lechner  T. Plath  University of
Hamburg, Department of Physics, Luruper Chaussee 149, 22761 Hamburg, Germany
S. Ackermann  J. Bödewadt  M. Vogt  Deutsches Elektronen-Synchrotron DESY,
Notkestraße 85, 22603 Hamburg, Germany

###### Abstract

Altering the optics in one section of a linear accelerator beamline will in
general cause an alteration of the optics in all downstream sections. In
circular accelerators, changing the optical properties of any beamline element
will have an impact on the optical functions throughout the whole machine. In
many cases, however, it is desirable to change the optics in a certain
beamline section without disturbing any other parts of the machine. Such a
local optics manipulation can be achieved by adjusting a number of additional
corrector magnets that restore the initial optics after the manipulated
section. In that case, the effect of the manipulation is confined in the
region between the manipulated and the correcting beamline elements.
Introducing a manipulation continuously, while the machine is operating,
therefore requires continuous correction functions to be applied to the
correcting quadrupole magnets. In this paper we present an analytic approach
to calculate such continuous correction functions for six quadrupole magnets
by means of a homotopy method. Besides a detailed derivation of the method, we
present its application to an algebraic example, as well as its implementation
at the seeding experiment sFLASH at the free-electron laser FLASH located at
DESY in Hamburg.

###### pacs:

41.85.-p

##  I Introduction

Many applications of particle accelerators require the dynamic manipulation of
the optical functions in certain regions of the beamline during operation as,
for example, minimizing the $\beta$-function at the interaction point of a
collider experiment, or closing a variable-gap undulator in a synchrotron
light source or a free-electron laser (FEL). Changing the optical properties
of a beamline element, however, causes not only local changes in the optics
but has an impact on the optical functions in all sections downstream of the
adjusted element. Hence, in many cases a manipulation will result in unmatched
optics in these sections and a correction is required to rematch the optics.

Our approach is to completely confine the influence of the manipulation to a
region around the adjusted elements, making the manipulation transparent for
all downstream sections. This can be achieved by appropriately adjusting
adjacent quadrupole magnets. Accelerator simulation codes such as ELEGANT
Borland (2000) or MAD-X Grote  _et al._ (2016) can perform this matching and
determine a suitable correction by numerically solving minimization problems.
Such an approach, however, is ill-suited for calculating continuous correction
functions, which are required to compensate a manipulation that is being
introduced continuously, as we will discuss later. In contrast to this, we
present an analytic approach that is based on the implicit function theorem
(compare for instance Krantz and Parks (2002)) and allows for the
determination of correction parameters as a continuous function of the
introduced disturbance.

Moreover, the implementation of this correction method at the FEL user
facility FLASH is presented. As depicted in Figure 1, FLASH features two
parallel undulator beamlines FLASH1 and FLASH2. Upstream of the FLASH1 main
undulator the sFLASH experiment is located, which, since its installation in
2010, utilizes a variable-gap undulator system to generate seeded FEL
radiation using different seeding techniques Bödewadt  _et al._ (2015). By
closing this undulator in order to start seeded operation, the optics in the
following FLASH1 main undulator is altered, which results in a deteriorated
FEL performance. The presented method was developed to provide the means to
correct the influence of the variable-gap undulator and thus allow
simultaneous operation of FLASH1 and sFLASH. Recently, the precise restoration
of the initial optics could be verified by beam size measurements along the
FLASH1 main undulator, which are presented in the final section of this paper.

Figure 1: Schematic layout of the FLASH facility. The superconducting linear
accelerator delivers electron bunches for the FLASH1 and FLASH2 undulator
beamlines and for the seeding experiment sFLASH.

##  II The method

In linear approximation the effect of a magnetic structure on the transversal
phase space vector of a particle can be described by a $4\times 4$ matrix
$\mathcal{M}$, called the transport matrix Courant and Snyder (1958). It
becomes apparent that any transport matrix is bound to be symplectic,
$\mathcal{M}\in\mathrm{Sp}(4,\mathbb{R})$. Furthermore it can be shown that
the group of symplectic $2N\times 2N$ matrices is a $N(2N+1)$-dimensional
manifold Gilmore (2005),

$\mathrm{dim}\,\mathrm{Sp}(2N,\mathbb{R})=N(2N+1),$ (1)

with $N\in\mathbb{N}$. As a result of this, any symplectic matrix can be
uniquely defined by $N(2N+1)$ of their elements. A general transversal
transfer matrix features therefore $10$ independent elements. However, in many
practical cases the two transversal directions are decoupled. Then
$\mathcal{M}$ reduces to two symplectic block matrices,

$\mathcal{M}\in\mathrm{Sp}(2,\mathbb{R})\oplus\mathrm{Sp}(2,\mathbb{R}),$

each of which consists of 3 independent elements. Decoupled transfer matrices
hence can be fully defined by a total of $6$ elements.

Every beamline can be understood as a sequence of $P$ elements, each of which
is represented by a transport matrix $M_{\textrm{n}}$. In general, the optical
properties of an element can be altered during operation, e.g. by changing the
current of a quadrupole magnet. Formally, we will therefore treat every
transport matrix as a function of a real experimental parameter $\rho_{n}$

$\mathcal{M}_{\textrm{n}}=\mathcal{M}_{\textrm{n}}(\rho_{n}).$

The total optical properties of a beamline are therefore given by the product
of its constituent’s transport matrices

$\mathcal{M}=\prod_{n=0}^{P-1}\mathcal{M}_{P-n}(\rho_{P-n}).$ (2)

A variation of a subset of $S$ machine parameters

$\bm{\sigma}=\bm{\sigma}(t)=\left(\sigma_{1}(t),\dots,\sigma_{S}(t)\right)^{T},$
(3)

will therefore result in deviations from the initial optical properties of the
beamline. Here, $\bm{\sigma}(t)$ denotes the course of the disturbance
parameters as a function of a single real parameter $t$. In many experimental
situations such a disturbance is undesirable and the ability to freely change
certain machine parameters without influencing the overall optical properties
of a beamline section would be of great benefit. To achieve this, the
disturbance has to be corrected by another subset of machine parameters
$\bm{\tau}$. These parameters will have to be controlled such that they
fulfill the condition

$\mathcal{M}(\bm{\sigma}(t),\bm{\tau}(t))=\mathcal{M}(\bm{\sigma}^{0},\bm{\tau}^{0}),$
(4)

where $\bm{\sigma}^{0}=\bm{\sigma}(0)$ and $\bm{\tau}^{0}=\bm{\tau}(0)$ denote
the constant initial machine parameters. If this condition is fulfilled, the
transport matrix of the beamline is constant and in particular is not affected
by the disturbances.

As we have seen, $\mathcal{M}$ is a symplectic matrix and therefore has a
reduced number of independent elements. By defining a vector-valued function
$\bm{B}\colon\,\mathbb{R}^{P}\to\,\mathbb{R}^{N(2N+1)}$ with the components

$B_{i}(\bm{\sigma},\bm{\tau})=\mathcal{M}_{i}(\bm{\sigma},\bm{\tau})-\mathcal{M}_{i}(\bm{\sigma^{0}},\bm{\tau^{0}}),$

where the $\mathcal{M}_{i}$ are a set of independent elements of
$\mathcal{M}$, equation (4) can be written as

$\bm{B}(\bm{\sigma}(t),\bm{\tau}(t))=0.$ (5)

Even for simple beamlines sections, consisting only of a few elements,
$\bm{B}(\bm{\sigma},\bm{\tau})$ will be highly intricate function. In general,
it is therefore not possible to solve equation (5) for $\bm{\tau}(t)$ by
algebraic means, given an arbitrary disturbance function $\bm{\sigma}(t)$. For
a given set of numerical values for the disturbance parameters it is, however,
possible to find a suitable approximation for correction parameters by
utilizing numerical minimization methods. This approach is readily realized
with the help of existing accelerator simulation codes. To find a numerical
approximation to the correction function using this method, one would have to
compute the correction parameters for a number of disturbance values along
$\bm{\sigma}(t)$ and interpolate. However, two successive interpolation points
are uncorrelated, as they were found by two independent numerical processes.
Hence, even with a high number of interpolation points one can not expect that
this method will produce a smooth, well-behaved correction function. In
contrast to this, we will derive a method that yields continuous correction
functions. Loosely speaking, our approach is to track how a known root of
$\bm{B}$ evolves along $\bm{\sigma}(t)$, rather than to recalculate the root
for each step.

The mathematical foundation of this method is the implicit function theorem
Ref. Krantz and Parks (2002). It states, in short, that equation (5)
implicitly defines a unique correction function $\bm{\tau}(t)$ in a
neighborhood around those machine states
$(\bm{\sigma}^{\star},\bm{\tau}^{\star})$, at which the Jacobian matrix is
non-singular,

$\mathrm{det}\,\dfrac{\partial\bm{B}}{\partial\bm{\tau}}(\bm{\sigma}^{\star},\bm{\tau}^{\star})\neq
0,$ (6)

and equation (5) is fulfilled

$\bm{B}(\bm{\sigma}^{\star},\bm{\tau}^{\star})=0.$ (7)

For the implicit function theorem to be applicable $\bm{B}$ and $\bm{\tau}$
have to be of the same dimension. Only then the Jacobian matrix in equation
(6) is square and its determinant can be formed. The number of disturbance
parameters $\bm{\sigma}$, in contrast, is not restricted. In the case of
uncoupled matrices, this means that a disturbance in any number of parameters
may be corrected by adjusting only six correction parameters.

Let us for the moment assume that for a certain $\bm{\sigma}(t)$ the
corresponding correction function $\bm{\tau}(t)$ exists. By differentiating
equation (5) with respect to $t$

$\dfrac{\textrm{d}B}{\textrm{d}t}=0=\dfrac{\partial\bm{B}}{\partial\bm{\sigma}}\dfrac{\textrm{d}\bm{\sigma}}{\textrm{d}t}+\dfrac{\partial\bm{B}}{\partial\bm{\tau}}\dfrac{\textrm{d}\bm{\tau}}{\textrm{d}t},$
(8)

an expression for the derivative of the correction function can be obtained

$\dfrac{\textrm{d}\bm{\tau}}{\textrm{d}t}=-\left(\dfrac{\partial\bm{B}}{\partial\bm{\tau}}\right)^{-1}\dfrac{\partial\bm{B}}{\partial\bm{\sigma}}\,\dfrac{\textrm{d}\bm{\sigma}}{\textrm{d}t}.$
(9)

This defines a system of coupled non-linear ordinary differential equations
for the components of $\bm{\tau}(t)$. We have therefore shifted the problem of
algebraically solving the intricate system (5) for $\bm{\tau}(t)$ to finding
the solution of the even more intricate differential equations (9). While it
still will be unfeasible to find an analytic solution, this problem now lends
itself readily to well-known numerical approximation methods. This approach
can be considered a numerical homotopy method, which are extensively treated
in Ref. Allgower and Georg (2003).

However, any iterative solution method will inevitably fail at those points
where the Jacobian (6) is singular. If the correction function for a certain
set of correction parameters contains such a singularity, these parameters are
not suitable to fully compensate the disturbance. In this case, another set of
correction parameters has to be chosen.

##  III An example

(a) Exemplary initial configuration.

(b) Corrected state. Despite the change of sign of $f_{2}$, the overall
optical properties are identical to the initial state.

Figure 2: Schematic representation of the four-lens system in different
configurations.

Symbolic transfer matrices of real-world lattices typically consist of long
analytical expressions, which renders obtaining solutions to (9) very
challenging. To promote better understanding of the proposed method we shall
nevertheless provide the reader with an algebraic, worked out example of
application. For this purpose, we will apply our method to a beamline
consisting of four quadrupole magnets, each separated by a distance $l$. To
simplify the calculations in this example, we restrict the problem to
preserving the optical properties in just one plane of motion and treat the
quadrupole magnets as thin lenses. The symbolic $2\times 2$ transfer matrix
$\mathcal{M}$ in this plane of motion is given by

$\mathcal{M}=\mathcal{L}(f_{4})\cdot\mathcal{D}(l)\cdot\mathcal{L}(f_{3})\cdot\mathcal{D}(l)\cdot\mathcal{L}(f_{2})\cdot\mathcal{D}(l)\cdot\mathcal{L}(f_{1}),$

where $\mathcal{L}(f)$ is the known transport matrix of a thin lens with focal
strength $f$, while $\mathcal{D}(l)$ represents a drift space,

$\mathcal{L}(f)=\begin{pmatrix}1&0\\\
f&1\end{pmatrix}\qquad\mathcal{D}(l)=\begin{pmatrix}1&l\\\ 0&1\end{pmatrix}.$

Suppose the system is in an arbitrary initial configuration $f_{i}=f_{i}^{*}$,
as for example depicted in Figure 2(a). If now the focal strength of one lens
is altered, the resulting optics will deviate from the initial case. The goal
of this example is to calculate functions for the focal strengths $f_{1}$,
$f_{3}$ and $f_{4}$ that correct any influence variation of $f_{2}$ has on the
transfer matrix of the system.

According to equation (1), $\mathcal{M}$ has three independent elements.
Therefore, we can choose three arbitrary elements of $\mathcal{M}$ for the
definition of $\bm{B}$. We might choose

$\displaystyle\bm{B}=(\mathcal{M}_{11},\mathcal{M}_{12},\mathcal{M}_{22})^{T}=$
$\displaystyle\scalebox{0.95}{$\left(\begin{array}[]{c}l\left(f_{2}\left(f_{3}l+2\right)+f_{1}\left(2f_{3}l+f_{2}l\left(f_{3}l+2\right)+3\right)+f_{3}\right)+1\\\
l\left(2f_{3}l+f_{2}l\left(f_{3}l+2\right)+3\right)\\\
3f_{4}l+2f_{3}l\left(f_{4}l+1\right)+f_{2}l\left(f_{3}l+f_{4}l\left(f_{3}l+2\right)+1\right)+1\\\
\end{array}\right)$}.$

We define the correction parameters to be the focal strengths of the remaining
lenses

$\tau=(f_{1},f_{3},f_{4})^{T}.$

In this case, equation (9) yields the following differential equations

$\begin{pmatrix}\dfrac{\textrm{d}f_{1}}{\textrm{d}f_{2}}\\\
\dfrac{\textrm{d}f_{3}}{\textrm{d}f_{2}}\\\
\dfrac{\textrm{d}f_{4}}{\textrm{d}f_{2}}\end{pmatrix}=\left(\begin{array}[]{c}-\dfrac{f_{3}l+2}{\left(f_{2}l+2\right)\left(2f_{3}l+f_{2}l\left(f_{3}l+2\right)+3\right)}\\\
-\dfrac{f_{3}l+2}{f_{2}l+2}\\\
\dfrac{1}{2f_{3}l+f_{2}l\left(f_{3}l+2\right)+3}\\\ \end{array}\right).$

The solution to these equations for the initial conditions
$f_{i}(f_{2}^{*})=f_{i}^{*}$ are

$\displaystyle f_{1}$
$\displaystyle=\dfrac{f^{*}_{1}\left(f_{2}l+2\right)\left(2f^{*}_{3}l+f^{*}_{2}l\left(f^{*}_{3}l+2\right)+3\right)-\left(f_{2}-f^{*}_{2}\right)\left(f^{*}_{3}l+2\right)}{\left(f_{2}l+2\right)\left(2f^{*}_{3}l+f^{*}_{2}l\left(f^{*}_{3}l+2\right)+3\right)}$
$\displaystyle f_{3}$
$\displaystyle=\dfrac{f^{*}_{2}\left(f^{*}_{3}l+2\right)+2f^{*}_{3}-2f_{2}}{f_{2}l+2}$
$\displaystyle f_{4}$
$\displaystyle=\dfrac{f_{2}-f^{*}_{2}}{2f^{*}_{3}l+f^{*}_{2}l\left(f^{*}_{3}l+2\right)+3}+f^{*}_{4}.$

For every value of $f_{2}$ these functions yield values for the remaing focal
strengths, so that the optical properties of the system are equivalent to
those of the initial setting $f_{i}^{*}$, see for example Figure 2(b). It can
be seen that the first two functions have a pole at $f_{2}=-2/l$, independent
on the initial condition. At this point the Jacobian determinant in equation
(6) vanishes and no continuous correction beyond this point is possible.

##  IV Implementation at sFLASH

Figure 3: Positions of the quadrupole magnets close to the SFLASH undulator.
Note that the magnets in the SFUND section can not be used as correctors.
Magnets used in the example shown in Figure 4 are highlighted in red. Drawing
not to scale.

The super-conducting linear accelerator FLASH at DESY in Hamburg delivers
trains of electron bunches into two parallel undulator beamlines, FLASH1 and
FLASH2. While FLASH2 is currently being prepared for user operation, FLASH1
has been delivering soft x-ray photons for user experiments since 2005
Honkavaara  _et al._ (2014). Upstream of the FLASH1 main undulator an
experimental setup for FEL seeding (sFLASH) is located, see Figure 1. It
features four variable-gap undulator modules and thus is the ideal test bench
for the study of several different seeding schemes Bödewadt  _et al._ (2015).

Closing the seeding undulator intensifies its natural focusing effect and will
therefore impact the beam optics further downstream of this undulator and in
particular at the entrance of the main undulator. At lower beam energies the
change in the $\beta$-function can easily be in the order of several tens of
percent so that the matching into the FODO structure of the main undulator is
corrupted. The result is a significantly deteriorated FEL performance.
Therefore, the disturbance introduced by the seeding undulator has to be
compensated by adjusting the field strength of appropriate quadrupole magnets.

Initially, the presented method was developed to correct for the influence of
the variable-gap undulator used by the seeding experiment at FLASH. Due to
their magnetic structure, planar undulators exhibit a systematic focusing
effect in the vertical plane (here, denoted by $y$), perpendicular to their
pole faces. For our purpose, the additional weak defocusing in the horizontal
direction ($x$), stemming from the finite width of real-world undulators, can
be neglected. Employing the matrix formalism introduced above, the effect of
an undulator of length $l$ and undulator period $\lambda_{\mathrm{u}}$ on a
particle’s trajectory vector $u=\left(x,x^{\prime},y,y^{\prime}\right)^{T}$ is
represented by the transfer matrix

$\mathcal{M}_{\mathrm{und}}=\begin{pmatrix}1&l&0&0\\\ 0&1&0&0\\\
0&0&\mathrm{cos}\sqrt{\kappa}l&\frac{\mathrm{sin}\sqrt{\kappa}l}{\sqrt{\kappa}}\\\
0&0&-\sqrt{\kappa}\,\mathrm{sin}\sqrt{\kappa}l&\mathrm{cos}\sqrt{\kappa}l\end{pmatrix},$
(10)

where $\kappa$ is the focusing strength

$\kappa\approx 2\left(\dfrac{K\,\pi}{\gamma\lambda_{\mathrm{u}}}\right)^{2},$
(11)

defined by the undulator parameter $K$ and $\gamma$, the particle’s Lorentz
factor Quattromini  _et al._ (2012). By varying the value of the gap, the $K$
parameter of the undulator system used by the sFLASH experiment can be set to
a value between $0$ and $2.72$ Tischer  _et al._ (2010), introducing the
undesirable focusing effect. We therefore identify $K$ as the disturbance
parameter.

Our implementation utilizes the computer algebra system Mathematica Wolfram
Research, Inc. (2014), which allows for easily calculating the required
symbolic transfer matrix according to equation (2), as well as setting up the
differential equations (9), the numerical solution of which is then found by
means of a basic Runge-Kutta algorithm. There are $14$ quadrupole magnets in
the vicinity of the sFLASH undulator (see Figure 3) that are included in the
Mathematica model and are available as correctors. As our method produces
correction functions for exactly $6$ parameters there are

$\binom{14}{6}=3003$

combinations of corrector magnets to choose from. Investigation shows that a
lot of these combinations are prone to produce solutions featuring
impracticably high changes in the quadrupole strengths or solutions containing
a singularity, as mentioned above. However, we were able to identify a set of
known good corrector combinations for optics similar to the standard optics
used in nominal SASE operation. An example of the correction functions that
are obtained by this method is shown in Figure 4.

Surely, for any given manipulation, correction functions exist which involve
more than six corrector magnets. In that case, the solution of equation (5) is
no longer unique because, as we have seen, solutions can be constructed using
any combination of six of the magnets and leave the others constant, which
results in our method being not directly applicable. However, any number of
additional constraints, other than those given by the constancy of the
transfer matrix, can easily be incorporated into our method by appending them
to $\bm{B}$ in equation (5) and selecting the same number of additional
correction parameters.

Figure 4: Example of a correction function calculated using the presented
method. Here, the strengths $k$ of the six named quadrupole magnets are used
as correction parameters to compensate for the optics disturbance caused by
the sFLASH undulator being closed from $K=0$ to $K=K_{\textrm{max}}=2.63$ at a
beam energy of $693\text{\,}\mathrm{MeV}$.

While our method assures the constancy of the transfer matrix between the
first and the last corrector magnet, it however makes no statement about the
behavior of the optical functions within this correction section. Using the
particle tracking code ELEGANT, the optics can be checked for unwanted
features and the efficacy of the correction can be verified, see Figure 5.

Figure 5: Development of the $\beta$-functions in dependence on the
longitudinal coordinate $s$. The undisturbed state (dotted) is compared to the
disturbed state with (solid line) and without (dashed) application of the
correction. Grey areas in the upper plot mark the position of the variable-gap
undulator segments. In this example, a different set of correction magnets is
used. Their positions are indicated in the lower plot. Note that at the end of
the beamline the corrected optical functions agree with the open case.

##  V Measurements

Diagnostics in the FLASH1 main undulator area include wire scanners at the
ends of the undulator modules. With the help of these devices the beam
diameter can be determined along the undulator. The relationship between beam
size $\sigma$ and $\beta$-function in the vertical and horizontal direction is
given by the well known equation

$\sigma_{u}(z)=\sqrt{\varepsilon_{u}\,\beta_{u}(z)},\quad u\in\\{x,y\\},$ (12)

where $z$ is the longitudinal position and $\varepsilon_{u}$ denotes the beam
emittance. A successful correction can therefore be verified experimentally,
by comparing the beam sizes along the undulator.

The measurements presented in Figure 6 have been conducted with a 1.0-GeV-
electron beam, while the disturbance was introduced by closing all four
seeding undulator segments to their lowest gap value.

Figure 6: Transverse beam sizes measured using the wire scanners between the
six FLASH main undulator segments (grey boxes in the middle). The dotted curve
shows the beam sizes at theory optics without any disturbance. Measurements of
the disturbed case, where the seeding undulator is closed, are depicted by the
dashed curve. The solid curve shows beam sizes with corrections applied. At
the time of measurement, the wire scanner between the last two undulator
segments was unavailable [Figures\ref{fig:correction};
\ref{fig:betaplot}and\ref{fig:wirescans}werecreatedusingmatplotlib.]Hunter2007.

As expected, closing the seeding undulator primarily affects the beam size in
the vertical direction. After the optics correction is applied, the beam sizes
correspond to those measured with an open undulator within the range of
measuring accuracy. Thereby, it is shown experimentally that our method is
able to restore the initial course of the optical functions in the main
undulator with satisfactory precision.

As mentioned above, the motivation leading to the development of the presented
method was to allow simultaneous operation of both the sFLASH and FLASH1 main
undulator. In recent efforts to establish procedures for three beamline lasing
(FLASH1, FLASH2, sFLASH) the described method is in frequent use Plath  _et
al._ (2016). Once FLASH1 SASE operation is established, the sFLASH variable-
gap undulator segments are closed to their target gap values, causing the SASE
signal of the main undulator to vanish. Subsequently, the calculated
correction is applied to the chosen corrector quadrupole magnets. This on-line
procedure was able to recover 70% of the initial FEL photon pulse energy with
only a few orbit corrections applied.

##  VI Conclusion and Outlook

We developed a novel approach to compensate for variations in particle beam
optics that would otherwise disturb the optics in other parts of the machine,
if let uncorrected. It produces continuous, numerical functions for six
machine parameters, compensating the unwanted effect of an arbitrary number of
disturbance parameters. Because of the continuity of the correction function,
it is suitable to be applied simultaneously to the disturbance. Effectively,
this method allows to freely modify the optics locally within a beamline
section, without altering the optics in the rest of the machine.

The method has been implemented at FLASH, with the objective of compensating
the influence of the variable-gap undulator used by the sFLASH experiment and
is in frequent use. We showed measurements confirming the precise restoration
of the initial optics.

A next experimental step will be to apply the correction steadily as the
disturbance arises, which in the end enables us to close the seeding undulator
and start seeding experiments without any notable changes of electron beam
quality in the FLASH1 main undulator. This compensation constitutes an
important milestone for the planned simultaneous operation of a seeding
experiment at FLASH1. The method will generally be useful for two FELs
operated independently of each other in cascaded mode. Ultimately, it will be
beneficial for all optics manipulations that need to stay restricted to a
local section of a long beamline.

###### Acknowledgements.

We would like to thank all colleagues participating in FLASH operation. This
work has been supported by Federal Ministry of Education and Research of
Germany under contract No. 05K1GU4 and the German Research Foundation
programme graduate school 1355.

## References

  * Borland (2000) M. Borland, “ELEGANT: A flexible SDDS-compliant code for accelerator simulation,” (2000).
  * Grote  _et al._ (2016) H. Grote _et al._ ,  _The MAD-X Program (Methodical Accelerator Design)_ , CERN, Geneva (2016).
  * Krantz and Parks (2002) S. Krantz and H. Parks,  _The Implicit Function Theorem: History, Theory, and Applications_ (Birkhäuser, 2002).
  * Bödewadt  _et al._ (2015) J. Bödewadt _et al._ , in  _Proceedings of IPAC2015_ , TUPC3 (2015).
  * Courant and Snyder (1958) E. Courant and H. Snyder, Annals of Physics 3, 1 (1958).
  * Gilmore (2005) R. Gilmore,  _Lie Groups, Lie Algebras, and Some of Their Applications_ , Dover Books on Mathematics (Dover Publications, 2005).
  * Allgower and Georg (2003) E. Allgower and K. Georg,  _Introduction to Numerical Continuation Methods_ , Classics in Applied Mathematics (Society for Industrial and Applied Mathematics, 2003).
  * Honkavaara  _et al._ (2014) K. Honkavaara _et al._ , in  _Proceedings of FEL2014_ , WEB05 (2014).
  * Quattromini  _et al._ (2012) M. Quattromini _et al._ , Phys. Rev. ST Accel. Beams 15, 080704 (2012).
  * Tischer  _et al._ (2010) M. Tischer _et al._ , in  _Proceedings IPAC2010_ , WEPD014 (2010).
  * Wolfram Research, Inc. (2014) Wolfram Research, Inc., “Mathematica,” (2014).
  * Hunter (2007) J. D. Hunter, Computing In Science & Engineering 9, 90 (2007).
  * Plath  _et al._ (2016) T. Plath _et al._ , TBP (2016).


# Shortcut to Chemically Accurate Quantum Computing via Density-based Basis-
set Correction

Diata Traore  Sorbonne Université, LCT, UMR 7616 CNRS, 75005 Paris, France
Qubit Pharmaceuticals, Advanced Research Department, 75014 Paris, France
Olivier Adjoua  Sorbonne Université, LCT, UMR 7616 CNRS, 75005 Paris, France
César Feniou  Sorbonne Université, LCT, UMR 7616 CNRS, 75005 Paris, France
Qubit Pharmaceuticals, Advanced Research Department, 75014 Paris, France
Ioanna-Maria Lygatsika  Present address: CEA, DAM, DIF, F-91297 Arpajon,
France; and Université Paris-Saclay, LMCE, 91680 Bruyères-le-Châtel, France
Sorbonne Université, LJLL, UMR 7598 CNRS, 75005 Paris, France  Sorbonne
Université, LCT, UMR 7616 CNRS, 75005 Paris, France  Yvon Maday  Sorbonne
Université, LJLL, UMR 7598 CNRS, 75005 Paris, France  Institut Universitaire
de France, 75005 Paris, France  Evgeny Posenitskiy  Qubit Pharmaceuticals,
Advanced Research Department, 75014 Paris, France  Kerstin Hammernik  NVIDIA
Corporation, Santa Clara, CA, USA  Alberto Peruzzo  Qubit Pharmaceuticals,
Advanced Research Department, 75014 Paris, France  Julien Toulouse  Sorbonne
Université, LCT, UMR 7616 CNRS, 75005 Paris, France  Institut Universitaire de
France, 75005 Paris, France  Emmanuel Giner  Sorbonne Université, LCT, UMR
7616 CNRS, 75005 Paris, France  Jean-Philip Piquemal  jean-
<EMAIL_ADDRESS>Sorbonne Université, LCT, UMR 7616
CNRS, 75005 Paris, France  Qubit Pharmaceuticals, Advanced Research
Department, 75014 Paris, France

###### Abstract

Quantum computing promises a computational advantage over classical methods in
electronic-structure calculations, with expected applications in drug design
and materials science. Accessing a quantitative description of chemical
systems while minimizing quantum resources, such as the number of qubits, is
an essential challenge given the limited capabilities of current quantum
processors. We provide a shortcut towards quantum computations at chemical
accuracy by approaching the complete-basis-set (CBS) limit through integrating
density-functional theory into quantum algorithms via density-based basis-set
corrections coupled to basis-sets crafted on-the-fly and specifically adapted
to a given system/user-defined qubit budget. The approach self-consistently
accelerates the basis-set convergence, improving electronic densities, ground-
state energies, and first-order properties such as dipole moments. It can also
serve as a classical, a posteriori, energy correction to quantum hardware
calculations. The strategy is assessed using GPU-accelerated state-vector
emulation up to 32 qubits. We converge the ground-state energies of four
atomic and molecular systems (He, Be, H2, LiH) within chemical accuracy of the
CBS full-configuration-interaction reference, while offering a systematic
increase of accuracy beyond a double-zeta basis-set quality for various
molecules up to the H8 hydrogen chain. We also obtain dissociation curves for
H2 and LiH that reach the CBS limit whereas for the challenging simulation of
the N2 triple-bond breaking we achieve a near-triple-zeta basis-set quality at
the cost of a minimal basis set. This hybrid strategy allows us to obtain
quantitative results that would otherwise require brute-force quantum
calculations using far more than 100 logical qubits, thereby opening up
opportunities to explore real-world chemistry with reasonable computational
resources.

Quantum computing (QC) offers a promising approach to solving electronic-
structure problems, with algorithms like Quantum Phase Estimation (QPE) and
Variational Quantum Eigensolver (VQE) proving effective in performing ground-
state quantum-chemistry wave-function calculations [1, 2, 3, 4, 5, 6]. The
electronic Hamiltonian is expressed in second quantization, employing an
encoding that maps one spin-orbital to one qubit. This mapping allows one to
represent an exponentially large Hilbert space using only a linear number of
qubits. However, to achieve accurate and practically valuable predictions for
chemical systems in real-world applications, the molecular Hamiltonian should
be expressed and solved using extensive basis sets of one-electron orbital
functions. The number of qubits required for such calculations quickly exceeds
the available capacities on current Noisy Intermediate Scale Quantum (NISQ)
devices, upcoming early Fault-Tolerant Quantum Computing (FTQC) devices, and
high-performance classical emulators. Therefore, while quantum chemistry has
long been stated as a promising application for quantum computing, the
endeavors have so far been limited to small molecular systems and minimal
basis sets. However, in practice, minimal basis sets fail to be predictive for
ground-state energies and chemically useful calculations require at least
significantly larger double-zeta basis sets. Requirements for computing
molecular properties are even more drastic as they tend to converge more
slowly with the size of the basis set than energies.

Figure 1: Architecture of the hybrid quantum-classical scheme introducing the
DBBSC method to wave-function QC calculations. The QPU workload can be either
performed by quantum hardware or replaced by a GPU-accelerated quantum
emulator. The input are the definition of the system and the basis set
(standard basis sets or our SABS can be used). In Strategy 1, the QPU/GPU
computes the expectation value of the standard Hamiltonian
$\langle\hat{H}\rangle$ and the non-self-consistent correlation basis-set
correction and the HF basis-set correction are shifted to the CPU of the
classical computer. In Strategy 2, the model includes a self-consistent basis-
set correction potential which is iteratively optimized between QPU/GPU and
CPU. The insets depict the effect of the basis-set corrections for the two
strategies.

In electronic-structure theory, the exact solution is defined by the full-
configuration interaction (FCI) method in the (infinite) complete-basis-set
(CBS) limit. However, in practice, the solution derived from a finite basis
set is used, which inherently suffers from truncation errors. This errors can
be substantial for small basis sets, but for large enough basis sets it is
possible to reach the target of chemical accuracy on energy differences, i.e.
1 kcal/mol (1.6 mHa). Unfortunately, employing sufficiently large basis sets
becomes very rapidly impractical for large systems, particularly those with
more than a few dozen of atoms. As a result, traditional quantum chemistry has
developed a variety of computational methods to approach chemical accuracy at
a reasonable cost.

In particular, basis-set correction techniques such as the density-based
basis-set correction (DBBSC) method [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19], relying on density-functional theory (DFT), have proven effective for
calculating ground- and excited-state energies, and also dipole moments, for a
variety of systems including atoms, small organic molecules, and strongly
correlated systems. To date, the DBBSC method has only been applied to still
relatively large basis sets, more specifically the family of Dunning basis
sets [20] for which chemical accuracy with respect to the CBS limit can be
reached starting from a triple-zeta basis set. We note that the approximations
developed within the DBBSC method target the CBS limit within a given wave-
function ansatz. This means that they do not address the intrinsic errors of
the wave-function method itself, such as the effects of neglected higher
excitations in truncated configuration-interaction (CI) or coupled-cluster
methods.

In this work, we propose the integration of the DBBSC method with quantum
algorithms to expedite reaching the CBS limit and achieve chemical accuracy on
complex molecular systems. This strategy, which natively limits the required
qubit counts, can be applied to quantum algorithms that tackle the ground-
state quantum chemistry problem, such as the QPE or VQE algorithms. While QPE
can guarantee ground-state energy with arbitrary high precision given a
carefully chosen initial state, it demands large quantum circuits and will
only be viable in the FTQC era. Conversely, VQE lacks convergence guarantees
but employs smaller circuits, aligning better with this study’s aim of
advancing short-term quantum computers toward practical quantum-chemistry
applications.

The paper is organized as follows. In Section I, we present two variants of
the application of the DBBSC method to wave-function QC calculations, denoted
as Strategy 1 and Strategy 2: (1) a basis-set correction a posteriori added to
the solution of the quantum algorithm, integrating two contributions, namely a
basis-set correlation density-functional correction and a basis-set Hartree-
Fock (HF) correction, and (2) a self-consistent scheme integrating the DBBSC
method to the quantum algorithm that dynamically modifies the one-electron
density used in the basis-set correlation correction. Strategy 1 offers the
possibility to correct any wave-function QC energy calculation through a
simple additive correction. Strategy 2 enables one to self-consistently access
to an improved electronic density, offering both improved energies and first-
order molecular properties. Also, we introduce a new type of system-adapted
basis sets (SABS) for Gaussian-type orbitals (GTOs), with size comparable to
minimal basis sets. These new methodologies enable us to perform the first
investigation of DBBSC corrections in the minimal basis-set regime.

In Section II, we provide relevant tests of our approach for different atomic
and molecular systems using several families of basis sets. We carry out
numerical simulations using graphics-processing-unit (GPU) accelerated QC
sparse emulation on up to 32 qubits, exploring the applicability of this
method to converge ground-state energies, dissociation curves, and dipole
moments. We consistently observe significant improvements over typical
quantum-algorithm approaches, reaching an accuracy level that would have
otherwise required hundreds of qubits. We thus expect the present approach to
become a standard part of quantum-enhanced wave-function calculations,
particularly the approach of Strategy 1 which can immediately provide large
improvements to existing results with relatively simple efforts.

Finally, Section III contains our conclusions. Additional details and results
are provided in the Supplementary Information (SI).

##  I Methodology

The overall methodological procedure is illustrated in Figure 1. The procedure
starts with the definition of the system, and a standard basis set or our
SABS. After defining the second-quantized Hamiltonian, the quantum state is
prepared on a quantum processing unit (QPU) or a GPU-accelerated quantum
emulator. The basis-set corrections are calculated on the classical CPU using
one of the two DBBSC variants labeled as Strategy 1 and Strategy 2.

###  I.1 Density-based basis-set correction method

In the infinite-dimensional (antisymmetric) $N$-electron Hilbert space, ${\cal
H}=\bigwedge^{N}L^{2}(\mathbb{R}^{3},\mathbb{C})$, we consider an atomic or
molecular system with Hamiltonian

$\hat{H}=\hat{T}+\hat{W}_{\text{ee}}+\hat{V}_{\text{ne}},$ (1)

where $\hat{T}$ is the kinetic-energy operator, $\hat{W}_{\text{ee}}$ is the
Coulomb electron-electron operator, and $\hat{V}_{\text{ne}}$ is the nuclei-
electron potential operator. The exact ground-state energy is defined as

$E_{0}=\min_{\Psi\in{\cal W}}\langle\Psi,\hat{H}\Psi\rangle,$ (2)

where ${\cal
W}=\\{\Psi\in\bigwedge^{N}H^{1}(\mathbb{R}^{3},\mathbb{C})\;|\;\langle\Psi,\Psi\rangle=1\\}$
is the space of admissible wave functions, $H^{1}(\mathbb{R}^{3},\mathbb{C})$
is the first-order Sobolev space, and $\langle\cdot,\cdot\rangle$ designates
the standard inner product of ${\cal H}$. In quantum-chemistry calculations,
we normally introduce a one-electron basis set ${\cal B}\subset
H^{1}(\mathbb{R}^{3},\mathbb{C})$ and we work in the finite-dimensional
$N$-electron Hilbert space generated by this basis set, i.e. ${\cal H}^{\cal
B}=\bigwedge^{N}\text{span}({\cal B})$. The FCI ground-state energy is then
defined as

$E_{\text{FCI}}^{\cal B}=\min_{\Psi\in{\cal W}^{\cal
B}}\langle\Psi,\hat{H}\Psi\rangle,$ (3)

where ${\cal W}^{\cal B}=\\{\Psi\in{\cal H}^{\cal
B}\;|\;\langle\Psi,\Psi\rangle=1\\}$ is the space of ${\cal B}$-representable
wave functions. In the CBS limit, the FCI ground-state energy tends to the
exact ground-state energy, i.e. $E_{\text{FCI}}^{{\cal
B}\to\text{CBS}}=E_{0}$, but the convergence is infamously slow due to short-
range electron correlation [21].

In the DBBSC method [7, 13], for the given basis set ${\cal B}$, we introduce
the following approximation to the ground-state energy

$E_{0}^{\cal B}=\min_{\Psi\in{\cal W}^{\cal
B}}\left(\langle\Psi,\hat{H}\Psi\rangle+\bar{E}^{\cal B}[n_{\Psi}]\right),$
(4)

where $\bar{E}^{\cal B}[n_{\Psi}]$ is a basis-set correction density
functional evaluated at the one-electron density of $\Psi$. This basis-set
correction density functional must vanish in the CBS limit so that
$E_{0}^{\cal B}$ properly converges to the exact ground-state energy, i.e.
$E_{0}^{{\cal B}\to\text{CBS}}=E_{0}$, but it must be such that it accelerates
the convergence to the CBS limit. In Ref. [8], based on the Perdew-Burke-
Ernzerhof (PBE) correlation density functional [22], such a basis-set
correction density functional was constructed in a semilocal form

$\bar{E}^{\cal B}[n]=\int_{\mathbb{R}^{3}}\bar{e}^{\cal
B}(n(\textbf{r}),\nabla n(\textbf{r}))\text{d}\textbf{r},$ (5)

where $n(\textbf{r})$ and $\nabla n(\textbf{r})$ are the density and density
gradient at point r, and the function $\bar{e}^{\cal B}(n,\nabla n)$,
including how the dependence on the basis set ${\cal B}$ is included, can be
found in Ref. [8]. We will now discuss how the DBBSC method can be adapted to
wave-function QC calculations.

####  I.1.1 Strategy 1: a posteriori basis-set correction

In Strategy 1, we use a non-self-consistent approximation to Eq. (4). The idea
is to calculate the FCI energy $E_{\text{FCI}}^{\cal B}$, or a good
approximation to it, on a quantum computer, and then add a posteriori the
basis-set correlation correction of Eq. (5) and a HF basis-set correction,
both calculated on a classical computer. This leads to the following
approximation to the ground-state energy

$E_{1}^{\cal B}=E_{\text{FCI}}^{\cal B}+\bar{E}^{\cal B}[n_{\text{HF}}^{\cal
B}]+\Delta E_{\text{HF}}^{\cal B}.$ (6)

As in previous works [7, 8], the basis-set correlation correction is evaluated
at the HF density in the basis set ${\cal B}$, leading to a calculation of
$\bar{E}^{\cal B}[n_{\text{HF}}^{\cal B}]$ with a marginal computational cost
with respect to the cost of calculating $E_{\text{FCI}}^{\cal B}$. The PBE-
based basis-set correlation correction only corrects for basis-set
incompleteness errors due to short-range electron correlation. However, when
using small basis sets, a significant part of the basis-set error also comes
from the fact the HF part of the energy is not converged to the CBS limit.
Similarly to Refs. [18, 19], we thus add a HF basis-set correction. In the
present work, we choose it simply as the difference between the HF energy in
the CBS limit, that we estimate as the value $E_{\text{HF}}^{\text{5Z}}$
obtained with the cc-pV5Z basis set, and the HF energy
$E_{\text{HF}}^{\mathcal{B}}$ in the basis set $\mathcal{B}$

$\Delta
E^{\mathcal{B}}_{\text{HF}}=E_{\text{HF}}^{\text{5Z}}-E_{\text{HF}}^{\mathcal{B}}.$
(7)

We note that it is also possible to avoid performing a HF calculation with the
cc-pV5Z basis set for estimating the CBS limit of the HF energy by using a
complementary auxiliary basis set, as in Refs. [18, 19], but we have not found
it necessary for the present work.

####  I.1.2 Strategy 2: Self-consistent basis-set correction

In Strategy 2, we use the self-consistent version of the DBBSC method in Eq.
(4), using the basis-set correlation correction of Eq. (5) in the self-
consistent part of the calculation, and we only add a posteriori the fixed HF
basis-set correction. This leads to the following approximation to the ground-
state energy

$E_{2}^{\cal B}=\min_{\Psi\in{\cal W}^{\cal
B}}\left(\langle\Psi,\hat{H}\Psi\rangle+\bar{E}^{\cal
B}[n_{\Psi}]\right)+\Delta E^{\mathcal{B}}_{\text{HF}}.$ (8)

The minimization in Eq. (8) leads to the following self-consistent Schrödinger
equation [13]

$\hat{P}^{\cal B}\hat{\tilde{H}}^{\cal B}[n_{\Psi^{\cal B}}]\Psi^{\cal
B}={\cal E}^{\cal B}\Psi^{\cal B},$ (9)

where $\hat{P}^{\cal B}$ is the projector on the $N$-electron Hilbert space
generated by the basis set ${\cal B}$, i.e. ${\cal H}^{\cal B}$, and
$\hat{\tilde{H}}^{\cal B}[n]$ is the effective Hamiltonian

$\hat{\tilde{H}}^{\cal B}[n]=\hat{H}+\hat{\bar{V}}^{\cal B}[n].$ (10)

Here, $\hat{\bar{V}}^{\cal B}[n]$ is the basis-set correction one-electron
potential operator

$\hat{\bar{V}}^{\cal
B}[n]=\int_{\mathbb{R}^{3}}\bar{v}^{\mathcal{B}}[n](\textbf{r})\;\hat{n}(\textbf{r})\text{d}\textbf{r},$
(11)

where
$\bar{v}^{\mathcal{B}}[n](\textbf{r})=\delta\bar{E}^{\mathcal{B}}[n]/\delta
n(\textbf{r})$ is the derivative of the basis-set correlation correction with
respect to the density, and $\hat{n}(\textbf{r})$ is the density operator.

The idea is now to solve iteratively Eq. (9) on a quantum computer. The
potential $\bar{v}^{\mathcal{B}}[n](\textbf{r})$, and therefore the
Hamiltonian $\hat{\tilde{H}}^{\cal B}[n]$, is iteratively updated with the
density $n_{i}^{\mathcal{B}}$ of the wave-function ansatz solution
$\Psi^{{\cal B}}_{i}$ of the $i^{\text{th}}$ iteration. The convergence
criterion is reached when the difference between the last two iterated energy
eigenvalues ${\cal E}_{i}^{\cal B}$ is less than 0.1 mHa. In practice, this is
fulfilled after two iterations.

Finally, the self-consistent basis-set corrected energy $E_{2}^{\cal B}$ can
be directly computed from the last energy eigenvalue ${\cal E}^{\cal B}$ and
density $n^{\cal B}=n_{\Psi^{\cal B}}$ coming out from the quantum solver as
follows

$E_{2}^{\mathcal{B}}={\cal E}^{\cal B}+\bar{E}^{\mathcal{B}}[n^{\cal
B}]-\int_{\mathbb{R}^{3}}\bar{v}^{\mathcal{B}}[n^{\cal B}](\textbf{r})n^{\cal
B}(\textbf{r})\text{d}\textbf{r}+\Delta E^{\mathcal{B}}_{\text{HF}},$ (12)

where the HF basis-set correction $\Delta E^{\mathcal{B}}_{\text{HF}}$ is
identical to the one used in Strategy 1.

####  I.1.3 GPU-accelerated QPU emulation and shift of the basis-set
correction to classical resources

The DBBSC method does not increase the qubit count as the basis-set correction
is performed classically using a density functional and a HF correction. In
practice, such computations present no advantage to be performed at the QC
level since they would consume a large number of qubits without accuracy
benefit over its classical counterpart. Therefore a hybrid QC wave-
function/classical DBBSC approach is highly preferable. In this work,
VQE/DBBSC computations are performed classically on GPU-accelerated computing
nodes, the GPUs replacing the quantum processing unit (QPU) through the use of
a classical quantum emulator. The basis-set correction contributions can be
then shifted to the unused CPUs to maximize the computational efficiency while
the computationally challenging VQE algorithm is fully offloaded on GPUs for
optimal time-to-solution performances. Of course, the proposed hybrid schemes
are also fully tractable on a real quantum computer, as the basis-set
correction contributions can be entirely shifted to classical resources
harnessing GPUs/CPUs to preserve the qubit count while the core wave-function
ansatz is maintained on the QPU.

###  I.2 System-adapted basis-set generation

Within the DBBSC method, the natural choice for basis sets is the Dunning
correlation-consistent basis-set family that offers a path to a reliable
convergence to the CBS limit, as demonstrated in initial DBBSC classical
quantum-chemistry studies [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19].
For QC calculations, such basis sets are unfortunately usually out of reach of
quantum hardware/emulators as they are associated with very large qubit
counts. In the present section, we address the task of generating atomic-
orbital (AO) basis sets under a basis-size budget, with controllable accuracy
[23]. We propose a mathematical formulation of this goal in terms of a
constraint optimization problem, as follows.

Given a molecule composed of $N_{\text{atm}}$ atoms, with fixed nuclear
coordinates $\\{\textbf{r}_{a}\\}_{1\leq a\leq N_{\text{atm}}}$, let

$\mathcal{B}=\\{\chi_{\mu}\\}_{1\leq\mu\leq
N_{\text{bas}}}=\bigcup_{a=1}^{N_{\text{atm}}}\\{\chi^{a}_{\mu}(\textbf{r}-\textbf{r}_{a})\\}_{1\leq\mu\leq
N^{a}_{\text{bas}}}$

denote a basis set of atom-centered GTOs. It is assumed that a one-electron
density, denoted by $n_{0}^{\mathcal{B}}$, is available after a converged
ground-state energy minimization procedure (e.g., HF or CI). Our main purpose
is to a posteriori extract subsets of the given “large” basis set
$\mathcal{B}$, denoted by $\mathcal{B}_{I}=\\{\chi_{\mu}\\}_{\mu\in
I}\subseteq\mathcal{B}$ for any index subset
$I\subseteq\\{1,\ldots,N_{\text{bas}}\\}$, achieving  (i) a target size, and
(ii) minimal accuracy loss on the given density $n_{0}^{\mathcal{B}}$ used as
a reference.  In the present computing scheme, since the DBBSC method uses the
cc-pV5Z basis set for the HF basis-set correction, we choose $\mathcal{B}$ =
cc-pV5Z. In practice, this choice is motivated by the fact that such a
quintuple-zeta basis set fairly approaches chemical accuracy compared to
quasi-exact all-electron numeric atom-centered orbital computations [24].
Hence, given a target basis size $M$ smaller than $N_{\text{bas}}$, we seek
the optimal index subset, denoted by $I_{M}$, that minimizes the best
approximation error of the density $n_{0}^{\mathcal{B}}$ over the set of
$\mathcal{B}_{I}$-representable one-electron densities, for any
$I\subseteq\\{1,\ldots,N_{\text{bas}}\\}$ with $|I|=M$, where $|\cdot|$
denotes the cardinal of a set, or, in other words, we solve an optimization
problem under constraints:

$I_{M}:=\arg\min_{\begin{subarray}{c}I\subseteq\\{1,\ldots,N_{\text{bas}}\\}\\\
|I|=M\end{subarray}}\,\min_{n^{\mathcal{B}_{I}}}\,\|n_{0}^{\mathcal{B}}-n^{\mathcal{B}_{I}}\|,$
(13)

where $\|\cdot\|$ is a given norm for functions over $\mathbb{R}^{3}$. In the
present work, we use the norm induced by the Coulomb inner product $\langle
u,v\rangle_{\text{c}}:=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\text{d}\textbf{r}\text{d}\textbf{r}^{\prime}\,u(\textbf{r})v(\textbf{r}^{\prime})|\textbf{r}-\textbf{r}^{\prime}|^{-1}$.
Let us emphasize that Eq. (13) requires knowledge of the density
$n_{0}^{\mathcal{B}}$, which is assumed to be pre-computed on a classical
computer. Such quantity is available in the framework of the DBBSC scheme
since HF computations are systematically performed at the cc-pV5Z level to
estimate the CBS limit of the HF energy.

The problem in Eq. (13) can be solved as follows. Our approach is to identify
a greedy procedure for discarding elements of the full AO-product set that
spans the space containing the reference density $n_{0}^{\mathcal{B}}$, which
admits the expansion

$n_{0}^{\mathcal{B}}=\sum_{\mu,\nu=1}^{N_{\text{bas}}}D_{\mu\nu}\chi_{\mu}\chi_{\nu},$

where $D=CC^{\top}$ is the density matrix and $C$ is the $N_{\text{bas}}\times
N_{\text{occ}}$ matrix of the coefficients in the AO basis of the
$N_{\text{occ}}$ occupied molecular orbitals. We plan to achieve this by
eliminating linear dependencies present in the AO-product set [25]. The
pivoted Cholesky decomposition (PCD) of a matrix [26] is an algebraic tool for
eliminating linear dependencies occurring between matrix rows (resp. columns),
which may be interpreted as an iterative greedy procedure for discarding
elements that do not contribute to the full row (resp. column) space, up to an
orthogonal projection error tolerance. PCD has been previously applied to the
auxiliary basis-set generation for density fitting [27, 28, 25]. In the
present work, we employ PCD within a new scheme, named system-adapted basis-
set (SABS) generation, for solving the problem in Eq. (13).

Table 1: Walltime (in minutes) required to grow the wave-function ansatz of
size $N_{\text{adapt}}$ for a given molecular system using Hyperion-1 state-
vector emulator on $N_{\text{gpus}}$ NVIDIA GPUs. The results have been
obtained on A100 (80 GB) and H100 (80 GB) GPUs using CUDA Toolkit 12.0 and
NVIDIA HPC SDK 23.3.

Molecule / Basis set | $N_{\text{adapt}}$ | $N_{\text{qubits}}$ | $N_{\text{gpus}}$ | A100 walltime [min] | H100 walltime [min]
---|---|---|---|---|---
H2O/6-31G | 500 | 24 | 1 | 644 | 503
H12/STO-3G | 500 | 24 | 1 | 174 | 134
H14/STO-3G | 300 | 28 | 8 | 283 | 184
H16/STO-3G | 100 | 32 | 128 | 450 | 147

Let us formulate our scheme in detail. Prior to the AO-product selection and
in order to ensure orbital symmetry of the resulting SABS, we pre-process the
initial basis $\mathcal{B}$ and first contract all angular components (e.g.
all three p-type components $\text{p}_{x},\text{p}_{y},\text{p}_{z}$) of GTOs.
To this end, we consider the partition $\\{B_{i}\\}_{i=1}^{N_{{\text{orb}}}}$
of $\\{1,\ldots,N_{\text{bas}}\\}$, each $B_{i}$ being an index block
containing all angular components of a single GTO in $\mathcal{B}$, and the
$N_{\text{orb}}\times N_{\text{bas}}$ contraction matrix $P$, defined for any
$1\leq i\leq N_{\text{orb}}$ as $P_{ij}=1$ if $j\in B_{i}$ and zero otherwise.
Next, we define the four-index tensor $T$ with entries

$T_{pqrs}=\sum_{\mu,\nu,\kappa,\lambda=1}^{N_{\text{bas}}}P_{p\mu}P_{q\nu}D_{\mu\nu}\,\langle\chi_{\mu}\chi_{\nu},\chi_{\kappa}\chi_{\lambda}\rangle_{\text{c}}\,D_{\kappa\lambda}P_{r\kappa}P_{s\lambda},$

and fold pairwise its first two and its last two dimensions, in order to form
the $N_{\text{prod}}\times N_{\text{prod}}$ Gram matrix of weighted AO
products, denoted by $G$, with $N_{\text{prod}}=(N_{\text{orb}})^{2}$. As a
last pre-processing step, we discard rows and columns of $G$ corresponding to
products made of components not centered on the same atom and denote the
resulting submatrix $A$. Now, PCD is applied to $A$, using the machine epsilon
as a tolerance threshold, yielding an index set sorted in pivot-descending
order. Given a target $M$, we define the selected AO-product index set,
denoted by $S^{\text{Chol}}_{M}$, as the $M$-first pivot indices. Lastly, we
recover the underlying AO index set

$J_{M}=\bigcup_{(p_{1},p_{2})\in S^{\text{Chol}}_{M}}\\{p_{1},p_{2}\\},$

and post-process it in order to ensure that the same GTO-parameter set is
assigned to all atoms of the same chemical type. For this purpose, the new
basis set associated to a chemical type is the union of parameters of those
GTOs in $\mathcal{B}_{J_{M}}$ that are centered on any atom of that type. This
yields a solution $I_{M}\supseteq J_{M}$ to our problem in Eq. (13) and the
resulting SABS is $\mathcal{B}_{I_{M}}$.

Note that our generation scheme directly fixes the size $M$ of the selected
products, i.e. $|S_{M}^{\text{Chol}}|=M$. The actual size of the AO basis set,
equal to $|I_{M}|$, is only implicitly controlled during our procedure. In
practice, as numerical results show, $|I_{M}|$ is very close to $M$ for s- and
p-type basis sets, while it remains the same order as $M$ for higher angular-
momentum orbital types. Overall, the SABS generation approach is extremely
fast and offers access to compact basis sets, specifically adapted to a given
system and user-defined qubit budget. Examples of SABS generation can be found
in the SI.

We note that the adaptation and generation of basis sets to the molecular
geometry has been the subject of several publications exploring other
strategies. Among them, we can cite some recent works [29, 30], but also some
other research in the context of quantum computing focusing on the need to
limit the qubit requirements through basis-set reoptimization such as Refs.
[31, 32]. Such reoptimizations cannot be performed on the fly and obtaining
such improvements can take a large amount of time when dealing with large
basis sets. Alternatively, Kottman and Aspuru [33] proposed a basis-set-free
approach through an adaptive representation using pair-natural orbitals which
was tested up to 22 qubits.

###  I.3 Computational Details

In the present study, we perform ADAPT-VQE computations [34] using the Qubit-
Excitation-Based pool of operators, which is considered a standard [35].
Additional details about the ADAPT-VQE methodology can be found in the SI. The
communication of the Hamiltonian from the CPU to the QPU/GPU is done by using
standard FCIDUMP files to communicate the one- and two-electron integrals to
the QPU/GPU software in order to construct the fermion operators. TREXIO files
[36] are also useful to communicate a wider range of relevant information.

ADAPT-VQE computations were performed using the Hyperion-1 GPU-accelerated
state-vector sparse emulator [37] up to 32 qubits. Hyperion-1 uses classical
computing systems and is grounded on an efficient multi-GPU ensemble of fast
custom sparse linear-algebra libraries accelerating Hyperion-1’s
exact/noiseless simulations. In this paper, computations were performed on
NVIDIA DGX A100 nodes (8x 80GB A100 GPUs per node) and NVIDIA DGX H100 nodes
(8x 80GB A100 GPUs per node). QC calculations being strongly memory-dependent,
a single GPU can carry out a 20-qubit ADAPT-VQE simulation depending on the
nature (i.e. Hamiltonian sparsity) of the system whereas a single node (8
GPUs) can handle up to 28 qubits. Multi-node computations are required beyond
such a qubit count. Further details about Hyperion-1 and its full capabilities
will be given in a forthcoming publication. Convergence for all ADAPT-VQE
computations were set to $10^{-6}$ Ha. Most computations started from a HF
initial state. For selected ones (indicated in the text and Tables), we
started the ADAPT-VQE procedure from a rough configuration-interaction
perturbatively-selected-iteratively (CIPSI) [38] initial state (converged to
only $10^{-2}$ Ha) to save some computational time within Hyperion-1. This
reflects a commonly adopted strategy where a multi-determinant initial state
is employed instead of a single HF determinant to increase the ground-state
support in the initial state. The Quantum State Preparation of such
classically-derived CIPSI wave functions has been studied in the context of
VQE [39] and QPE [40]. Also, we report in Table 1, the walltime required for
several ground-state energy calculations. As one can see, the walltime is not
only a function of the number of qubits, sparsity is also important. For
example, the Hamiltonian of the H12 molecule is sparser than the one of the
water molecule resulting in a faster convergence. Beside the increased
computing power of H100 GPUs leading to improved time-to-solution, our results
also highlight the importance of fast node-to-node interconnects when
performing large-scale quantum emulation. Indeed, the benefit of H100 over
A100 is striking for the largest 32 qubits simulations on 16 nodes where an
improvement of a factor 3 was observed on the DGX H100 systems. Such speedup
is therefore also partially related to higher node-to-node bandwidth observed
on DGX H100 versus A100 systems.

Figure 2: Ground-state energy errors with respect to the extrapolated CBS
limit. Outside the green box: ADAPT-VQE calculations. Inside the green box:
FCI calculations from a classical computer. Dot markers correspond to energies
without basis-set corrections. Cross markers correspond to energies corrected
with Strategy 1 [Eq. (6)]. On the x-axis, the labels $\mathcal{B}+$ symbolize
basis-set corrected values. The blue box corresponds to a range of 1.6 mHa
around the extrapolated CBS limit. The errors for each ADAPT-VQE computations
with respect to the FCI energy for a given basis set are reported in the SI.

##  II Results

We computed the ground-state energies, dissociation curves, and dipole moments
for the N2, H2O, LiH, and H2 molecules using both basis-set correction
strategies. For all systems except H2, the 1s molecular orbitals were frozen,
i.e. we use the frozen-core approximation. Correspondingly, we use the frozen-
core version of the DBBSC method as defined in Ref. [8], in which in
particular the contribution to the density coming from the core electrons is
neglected in the basis-set correlation density functional. The error in each
calculation is quantified as the deviation from the CBS limit, which itself is
determined using a two-point extrapolation scheme based on the cc-pVQZ and cc-
pV5Z basis sets [21]. For the dipole-moment calculations we employed only
Strategy 2, using an expectation value over the dipole-moment operator,
avoiding a finite-difference approach prone to numerical errors [41].

The classical computations, including basis-set corrections and the
calculation of reference energies and dipole moments, were carried out using
Quantum Package 2.0 [42]. In this software, the FCI energies are approximated
by the energy from a CIPSI wave function to which a second-order perturbation-
theory (PT2) correction is added. Given the demonstrated nearly-FCI quality of
this approximation for the systems we studied, we simplify our terminology by
referring to this approach as simply FCI rather than CIPSI+PT2.

The SABS are labeled VXZ-Y where X is respectively D, T, Q, 5, or 6 for cc-
pVDZ, cc-pVTZ, up to cc-pV6Z basis sets, and Y is the target size. We report
in the SI the classically computed energies for the N2, H2O, LiH, and H2
molecules.

###  II.1 Ground-state Energies

Ground-state energies of the H2, LiH, H2O, and N2 molecules close to their
respective equilibrium geometries are presented in Figure 2 and Table 2.
Details regarding the ADAPT-VQE iterations and the associated “ADAPT” values
are provided in the SI, as well as additional tests on atoms and hydrogen
chains.

In Figure 2, a general trend is observed: the basis-set corrected ground-state
energies with the small basis sets, i.e. STO-3G, 6-31G, or pcseg-0 [43], align
with values between the cc-pVDZ and cc-pVTZ basis-set levels, while requiring
much less qubits than these latter basis sets. As we can see from Table 2,
both strategies provide the same quantitative improvements. These results are
consistent with the conclusion from Ref. [13] where the non-self-consistent
approximation was found sufficient for calculating energies.

For the best cases, the basis-set corrected ground-state energies have errors
in the order of tens of mHa: 40 mHa for H2O with 24 qubits, 60 mHa for N2 with
16 qubits, and less than 10 mHa for LiH and H2 with less than 10 qubits. The
results stay consistent for the other systems available in the SI (H4, H6,
H8). For H2 using the cc-pVDZ basis set, the basis-set corrected ground-state
energy reaches chemical accuracy relative to the CBS limit, while requiring 20
logical qubits. Similar convergence to the CBS limit is observed in the SI for
He and Be. Without the basis-set correction, reaching a similar accuracy would
have required typically more than a hundred of qubits.

In the initial work on the DBBSC method [7], it was found that basis-set
corrected ground-state energies of atoms reach chemical accuracy with the cc-
pVTZ basis set. Similarly, for the molecules studied here, if one had access
to a few hundred logical qubits required for the cc-pVTZ basis set, achieving
chemical accuracy would be theoretically feasible for all cases. To support
this claim, we report in the SI classically computed basis-set corrected FCI
ground-state energies using Dunning basis sets. The error relative to the CBS
limit for the LiH molecule is 0.2 mHa with the cc-pVTZ basis set. For H2O and
N2, we achieve errors of 3 mHa and 1 mHa, respectively, with the cc-pVTZ basis
set.

Table 2: Ground-state energies (in Ha) for H2O, N2, LiH, H2, and H8 calculated
by FCI, ADAPT-VQE (denoted as ADAPT), and basis-set corrected ADAPT-VQE
according to Strategy 1, denoted as A+PBE+$\Delta$HF, and to Strategy 2,
denoted as SC(A+PBE) and SC(A+PBE)+$\Delta$HF (i.e., without and with the HF
basis-set correction, respectively). Here, PBE refers to the PBE-based
correlation basis-set correction and $\Delta$HF refers to the HF basis-set
correction, and SC stands for “self-consistent”. The frozen-core approximation
has been used for H2O, N2, and LiH. The CBS limits are estimated by two-point
extrapolations from cc-pVQZ and cc-pV5Z calculations.

H2O | $N_{\text{qubits}}$ | FCI | ADAPT | A+PBE+$\Delta$HF | SC(A+PBE) | SC(A+PBE)+$\Delta$HF
---|---|---|---|---|---|---
STO-3G | 12 | -75.01250 | -75.01250 | -76.30232 | -75.197880 | -76.30191
pcseg-0 | 24 | -75.90855 | -75.90843111500 iterations with a 14668-determinant CIPSI initial state. | -76.33681 | -76.03999 | -76.33279
6-31G | 24 | -76.11995 | -76.119892221000 iterations with a 10879-determinant CIPSI initial state. | -76.28035 | -76.23717 | -76.32025
V5Z-10 | 24 | -76.12626 | -76.12409 | -76.32705 | -76.27418333The value is -76.27636 when using a 11016-determinant CIPSI initial state. | -76.32505
V5Z-11 | 30 | -76.15902 | -76.15165 | -76.33704 | -76.28622 | -76.33570
cc-pVDZ | 46 | -76.24165 | - | - | - | -
VTZ | 114 | -76.33250 | - | - | - | -
cc-pVQZ | 228 | -76.35985 | - | - | - | -
cc-pV5Z | 400 | -76.36877 | - | - | - | -
CBS | - | -76.37812 | - | - | - | -
N2 | $N_{\text{qubits}}$ | FCI | ADAPT | A+PBE+$\Delta$HF | SC(A+PBE) | SC(A+PBE)+$\Delta$HF
STO-3G | 16 | -107.65251 | -107.65251 | -109.36630 | -107.86974 | -109.36661
V5Z-6 | 16 | -108.88869 | -108.88869 | -109.34552 | -109.09850 | -109.34608
cc-pVDZ | 52 | -109.27698 | - | - | - | -
cc-pVTZ | 116 | -109.37527 | - | - | - | -
cc-pVQZ | 216 | -109.40558 | - | - | - | -
cc-pV5Z | 360 | -109.41505 | - | - | - | -
CBS | - | -109.42498 | - | - | - | -
LiH | $N_{\text{qubits}}$ | FCI | ADAPT | A+PBE+$\Delta$HF | SC(A+PBE) | SC(A+PBE)+$\Delta$HF
STO-3G | 10 | -7.88218 | -7.88218 | -8.02160 | -7.89590 | -8.02119
pcseg-0 | 14 | -7.98139 | -7.98139 | -8.0160 | -7.99166 | -8.01561
6-31G | 20 | -7.99800 | -7.99800 | -8.01668 | -8.00806 | -8.01611
V5Z-4 | 10 | -7.99287 | -7.99287 | -8.01758 | -8.00562 | -8.01690
V5Z-7 | 16 | -7.99793 | -7.99793 | -8.01710 | -8.01109 | -8.01643
V5Z-10 | 28 | -8.01302 | -8.01302 | -8.02575 | -8.02134 | -8.02540
cc-pVDZ | 26 | -8.01438 | - | - | - | -
cc-pVTZ | 86 | -8.02234 | - | - | - | -
cc-pVQZ | 190 | -8.02386 | - | - | - | -
cc-pV5Z | 290 | -8.02433 | - | - | - | -
CBS | - | -8.02482 | - | - | - | -
H2 | $N_{\text{qubits}}$ | FCI | ADAPT | A+PBE+$\Delta$HF | SC(A+PBE) | SC(A+PBE)+$\Delta$HF
STO-3G | 4 | -1.13415 | -1.13415 | -1.17606 | -1.15590 | -1.17594
6-31G | 8 | -1.15003 | -1.15003 | -1.16911 | -1.16196 | -1.16913
cc-pVDZ | 20 | -1.16275 | -1.16275 | -1.17239 | -1.16858 | -1.17246
V5Z-8 | 24 | -1.16613 | -1.16613 | -1.17315 | -1.17170 | -1.17320
cc-pVTZ | 56 | -1.17041 | - | - | - | -
cc-pVQZ | 120 | -1.17182 | - | - | - | -
cc-pV5Z | 220 | -1.17223 | - | - | - | -
CBS | - | -1.17265 | - | - | - | -
H8 | $N_{\text{qubits}}$ | FCI | ADAPT | A+PBE+$\Delta$HF | SC(A+PBE) | SC(A+PBE)+$\Delta$HF
STO-3G | 16 | -4.24339 | -4.24320 | -4.45483 | -4.32764 | -4.45354
6-31G | 32 | -4.37032 | -4.35752 | -4.43488 | -4.41275 | -4.43502
cc-pVDZ | 80 | -4.42756 | - | - | - | -
cc-pVTZ | 222 | -4.47121 | - | - | - | -
cc-pVQZ | 474 | -4.47702 | - | - | - | -
cc-pV5Z | 864 | - | - | - | - | -
CBS | - | - | - | - | - | -

Let us now discuss the results obtained with our SABS for LiH, N2, and H2O.
They were chosen to match the sizes of the small basis sets previously
discussed. Specifically for H2O, employing the V5Z-10 basis set (24 qubits)
achieves results slightly superior to those obtained with the 6-31G basis set.
For H2O we observed a reduction in the HF basis-set correction by 40 mHa when
moving from the 6-31G to the V5Z-10 basis set, by over 100 mHa when comparing
the STO-3G to the V5Z-4 basis set for the LiH molecule, and more than an
Hartree for the N2 molecule moving from the STO-3G to the V5Z-6 basis set.
Based on these findings, we anticipate that further exploration of this
strategy could lead to a basis-set correction scheme requiring only the
correlation basis-set correction term. Of course, the interest of SABS is to
be systematically improvable toward the parent Dunning basis set. It is
important to note that this SABS systematic improvement is not limited to
quantum algorithms such as ADAPT-VQE and is also present for the HF and FCI
calculations. Since the SABS approach strongly reduces the computational cost
while maintaining accuracy, it should offer further applications of the CIPSI
approach in classical quantum chemistry. In practice, a cc-pVTZ-like quality
is reached for N2 with a V5Z-11 basis set (32 qubits) and for H2O with a
V5Z-11 basis set (30 qubits). A V5Z-10 basis set (28 qubits) achieves a cc-
pV5Z-like ground-state energy accuracy for LiH. Overall, the SABS always
provide the best ADAPT energies with the self-consistent correlation basis-set
correction (before addition of the HF basis-set correction).

Overall, concerning the basis sets, it is important to point out a key
anomaly, i.e. the remarkable performance of the minimal STO-3G basis set. This
was expected as Pople already discussed the outperformance of this basis set
in the seventies [44]. Indeed, Davidson and Feller detailed in their 1986
review [45] the existence of error compensations and stated that “the smaller
the basis the more ab initio calculations assume an empirical flavor”. In
practice, among the available minimal basis-set possibilities, STO-3G proved
to be the most robust and cost-efficient choice for the basis-set corrections.
Also, in the DBBSC schemes, the STO-3G basis set highly benefits from the HF
basis-set correction which reduces as the basis set and number of qubits
increases. When compared to SABS, for example, STO-3G always display smaller
self-consistent basis-set corrected ADAPT values and larger HF basis-set
corrections. In practice, Pople also noted that STO-3G is especially good for
energies around the equilibrium geometry, and larger basis sets are clearly
required to describe accurately the notoriously more difficult dissociation
curves and dipole moments.

Finally, one last aspect of the analysis of such large basis-set simulations
is related to the convergence of ADAPT-VQE. Indeed, in VQE-type computations
[5], there is no formal guarantee of convergence and ADAPT-VQE belongs to such
heuristic family of methods. If, for systems like H2, LiH or N2 (see SI,
Figures 3-14), convergence is obtained in a few dozen of iterations, for more
complex systems such as H2O, the number of required iterations strongly
increases when a double-zeta basis set is used. Clearly, Figure 2 in the SI
shows that more than a thousand iterations is required to achieve full
convergence starting from the HF initial state. This heuristic aspect makes
that no anticipation of the exact required number of iterations can be made.
This may be manageable on a real quantum processor, but the use of classical
emulation makes each iteration more costly in term of time-to-solution than
the previous one, limiting the overall convergence capabilities. It is
possible to force convergence by replacing the HF starting point by a CIPSI
one. Figure 2 in the SI shows that chemical accuracy can be easily reached
with a well-converged CIPSI solution. However, ADAPT-VQE struggles to improve
the solution which is expected due to the enormous size of the parameter
space. Furthermore, the CIPSI-based initial state already contains significant
amount of correlation. This represents a hard challenge for the ADAPT-VQE
procedure, which now needs to more carefully pick the next ansatz operator to
improve the existing quantum state. In any case, starting from an unconverged
CIPSI wave function is a robust solution to strongly reduce the overall
computational time. One example is given for H2O and the V5Z-10 basis set
which initially led to a not fully converged result. With a better CIPSI
starting point, it is possible to recover a few milli-Hartrees (see footnotes
in Table 2). It is also possible to change the operator pool but the point
here is that convergence becomes challenging when tackling complex electronic
structures and such computations would not be possible without GPU-accelerated
emulation.

Figure 3:  Dissociation curves of H2, LiH and N2 molecules. Blue and orange
lines correspond to Dunning basis sets using the FCI method. The notation VXZ
stands for the basis set cc-pVXZ. Green and purple plain lines correspond to
ADAPT-VQE calculations without basis-set corrections. Green and purple dash-
dotted lines correspond to basis-set corrected dissociation curves. Raw data
are available in the SI.

###  II.2 Dissociation Curves

We pursue with the dissociation energies reported in Figure 3. Clearly, the
simple ADAPT-VQE/STO-3G level of theory appears quite far from an accurate
description of the dissociation (compared to the calculations in the largest
basis sets), which is achieved with FCI/triple-zeta (and beyond) levels. The
non-corrected ADAPT-VQE values are represented with bold lines whereas the
corrected values are reported with dashed curves. These dissociation curves
are extremely challenging as they involve several regimes of correlation going
from weak correlation happening typically at short distances to strong
correlation effects at large distances. First, we notice that for all the
cases, the basis-set corrected values around the equilibrium are always
substantially closer to the large cc-pV5Z reference values than the
uncorrected ones and improve with the basis-set size and the use of SABS. For
larger distances, the basis-set requirements appear even more stringent. In
practice, the three dissociation curves manage to be converged to high
accuracy thanks to the use of SABS: we obtain nearly a triple-zeta quality for
N2 using the basis-set correction with the V5Z-6 basis set, nearly a cc-pV5Z
quality for H2 using the basis-set correction with the V5Z-8 basis set, and
nearly a cc-pV5Z quality for LiH using the basis-set correction with a V5Z-7
basis set. For H2, the basis-set corrected cc-pVDZ curve matches the cc-pV5Z
reference up to a distance of 2.5 Å and a slight discrepancy appears at long
distances. This is consistent with classical computations leading to a
convergence to the CBS limit only when the basis-set correction is applied to
the cc-pVTZ basis set (114 qubits). It is possible to fix this issue more
affordably by using a V5Z-8 basis set requiring only 24 qubits. For LiH, it is
also possible to fix the small discrepancies observed on the V5Z-7 basis-set
dissociation curve using a larger V5Z-10 basis set (28 qubits) which reaches
the CBS limit.

Table 3: Dipole moments (in atomic units) of LiH and H2O calculated by FCI and
self-consistently basis-set corrected ADAPT-VQE, denoted as SC(A+PBE) and
SC(A+PBE)+$\Delta$HF (without and with and the HF basis-set correction,
respectively). Here, $\Delta$HF corresponds to the HF basis-set correction to
the dipole moment, calculated as the difference between the HF dipole moment
in the aug-cc-pV5Z basis set and the HF dipole moment in the considered basis
set. More data are available in the SI.

LiH | $N_{\text{qubits}}$ | FCI | SC(A+PBE) | SC(A+PBE)+$\Delta$HF
---|---|---|---|---
STO-3G | 10 | -1.81835 | -1.86299 | -2.31321
pcseg-0 | 14 | -2.33313 | -2.37289 | -2.25650
6-31G | 20 | -2.16646 | -2.20768 | -2.23674
V5Z-4 | 10 | -2.37818 | -2.44145 | -2.22886
V5Z-7 | 16 | -2.23095 | -2.27789 | -2.25618
V5Z-10 | 28 | -2.24997 | -2.27458 | -2.31438
cc-pVDZ | 26 | -2.25566 | - | -
cc-pVTZ | 86 | -2.29998 | - | -
cc-pVQZ | 190 | -2.30361 | - | -
cc-pV5Z | 290 | -2.30647 | - | -
H2O | $N_{\text{qubits}}$ | FCI | SC(A+PBE) | SC(A+PBE)+$\Delta$HF
STO-3G | 12 | -0.63584 | -0.67084 | -0.77162
pcseg-0 | 24 | -0.95822 | -0.99450 | -0.77065
6-31G | 24 | -0.99020 | -1.01898 | -0.76342
V5Z-10 | 24 | -0.99305 | -1.01887 | -0.77634
V5Z-11 | 30 | -0.99185 | -1.02170 | -0.77737
cc-pVDZ | 46 | -0.76073 | - | -
cc-pVTZ | 114 | -0.75013 | - | -
cc-pVQZ | 228 | -0.74994 | - | -
cc-pV5Z | 400 | -0.74241 | - | -

Finally, it is important to highlight the good performance of the DBBSC method
for the prediction of the dissociation curve of the triple-bonded N2 molecule.
Indeed, such computation is well-documented in the literature and known as
extremely difficult as it requires both a multi-reference treatment and
various weak correlation effects going from short-range to charge polarization
effects [46, 47, 48, 49]. At the cost of a minimal basis set (i.e. 16 qubits),
our N2 DBBSC computations achieved nearly a cc-pVTZ quality, providing an
accuracy which would have required around 100 logical qubits in the context of
a brute-force simulation. To the best of our knowledge, these results are the
most accurate using a quantum algorithm when compared to the recent results
from the literature. Indeed, a group of researchers managed to predict this
dissociation curve using a Local Unitary Cluster Jastrow (LUCJ) ansatz coupled
to double-zeta basis sets (6-31G and cc-pVDZ) [50]. The computation required
the use of massive computational resources, namely hundreds of compute nodes
of the Fugaku classical supercomputer coupled to a QPU. Alternatively, in the
context of the fermionic quantum emulator, N2 simulations used ADAPT-VQE
coupled to the 6-31G and def2-SVP (56 qubits) basis sets to perform
computation via approximate Matrix Product States (MPS) [51]. In the present
work, we achieve a better accuracy with far less qubits and, again, these
results can be systematically improved using larger SABS. Indeed, going beyond
a cc-pVTZ basis set is possible using the next (larger) SABS in term of size,
i.e. V5Z-11 (32 qubits).

###  II.3 Dipole Moments

We also explore the idea to extend the basis-set correction scheme in QC
calculations of molecular properties such as the dipole moment. As pointed out
by Halkier et al. [52], the dipole moment also suffers from a slow basis-set
convergence, which can be thought of as an indirect impact of the missing
short-range correlation effects in a finite basis set. Therefore, it is
relevant to apply the basis-set correction to the dipole moment, as already
shown in Ref. [13]. We report in Table 3 the dipole moments of H2O and LiH. We
calculate the basis-set corrected dipole moment as expectation value over the
dipole-moment operator over the last ADAPT-VQE wave function coming out of the
self-consistent basis-set correction method. We also add a posteriori the HF
basis-set correction to the dipole moment, calculated as the difference
between the HF dipole moment in the aug-cc-pV5Z basis set and the HF dipole
moment in the considered basis set. From Table 3, we see that the correlation
basis-set correction is not sufficient to converge the dipole moments.
However, the addition of the HF basis-set correction significantly improves
the basis-set convergence. For both molecules, strong improvements are
observed with respect to the ADAPT-VQE/STO-3G level. Further improvements are
observed as the size of the basis set increases. Finally, we remark the ADAPT-
VQE slow convergence in term of iterations is also the source of small errors
in the dipole moments.

##  III Conclusion

We demonstrated the applicability of the DBBSC method to QC algorithms for
quantum chemistry. Using ADAPT-VQE, these approaches were shown to be able to
systematically improve minimal basis-set results by predicting ground-state
energies that are intermediate between double-zeta and triple-zeta FCI
qualities. Overall, the presented self-consistent basis-set corrected ground-
state energies are very close to their non-self-consistent counterparts. As
the latter a posteriori basis-set correction approach can be very easily
applied to any wave-function QC calculation performed on real quantum
hardware, it provides an affordable improvement strategy for current QC
chemistry computations. The self-consistent basis-set correction scheme is
still useful since it permits the calculation of properties such as dipole
moments thanks to the availability of improved QC densities.

An additional reduction of the required basis size is provided by our SABS
approach. Such a “black-box” pivoted-Cholesky strategy for the on-the fly
generation of basis sets has been shown to be competitive and often superior
to available standard choices. It offers systematic computational savings for
large basis sets, reducing significantly the qubit requirements. Thus the
SABS’ usefulness is not limited to QC algorithms as they also offer
systematically improvable solutions for performing reference classical
computations. Indeed, since SABS are user-defined truncated versions of the
Dunning basis sets that can match a qubit budget, they also offer a reduced
cost access to improved accuracy for any quantum-chemistry method.

Our methodology allows us to compute chemically meaningful energies and
properties on systems that would have required far more than 100 logical
qubits. For example, the computation of the H2 total energy at the FCI/cc-pV5Z
level that would have required more than 220 logical qubits, can be achieved
here with only 24 qubits using our basis-set correction scheme and our SABS
technique. Overall, we were able to converge four systems to the FCI/CBS limit
including He, Be, H2, and LiH. We were also able to provide accurate
dissociation curves for H2, LiH, and N2. Such computations required the use of
only a single GPU. Since most quantum-chemistry studies are presently out of
reach of quantum computers, this research opens the path to more affordable
quantitative quantum-chemistry simulations of small molecules using QC
algorithms. In particular, the present a posteriori basis-set corrections can
be easily added to any type of STO-3G VQE fermionic computations on real
hardware [5, 53, 54], allowing one to improve significantly their accuracy at
very little computational cost. Since adaptive simulations on real hardware
are making some progress [55] while hardware itself improves, basis-set
corrected simulations should be progressively possible on future quantum
computers providing a route to FCI/CBS quality computations. This strategy is
particularly suited for resources demanding computations that converge slowly
and require large basis sets associated to large qubits counts. Finally, the
DBBSC method is not limited to ground-state computations and can be extended
to excited states via a linear-response formalism[17]. This will be the
subject of future QC research.

In this context, beside the state-of-the-art ADAPT-VQE hybrid quantum-
classical algorithm, it would be interesting to revisit accurate fixed wave-
function ansätze [56], such as UCCSDT [57] and others, to analyze their
convergence when coupled to the DBBSC method. In practice, the presented
basis-set correction framework is not restricted to a given QC wave-function
ansatz and can leverage any future QC algorithmic improvements. Our present
DBBSC/SABS methodology still requires to use qubits either through quantum
hardware or through the use of a classical quantum emulator. The qubit count
is therefore our main limitation. The computations of this paper used up to 32
qubits and represent a proof-of-concept study of what one can do with quantum
emulation. However, state-vector simulations have limitations due to memory.
Indeed, if they are theoretically possible up to 40-50 qubits on very large
exascale supercomputers, they start to become relatively unpractical for
quantum-chemistry simulations when reaching 36 qubits due to the high
computational resources and time-to-solution requirements. To explore further
the chemical electronic space with quantum algorithms, we are currently
upgrading our Hyperion-1 framework to increase our emulated qubit counts by
going beyond the state-vector formalism thanks to various elements from the
CUDA-Q SDK [58, 59] developed by NVIDIA. We are also presently testing various
implementations of the DBBSC algorithms on available quantum hardware. To
conclude, by reducing the number of qubits required to reach the CBS limit, we
expect to tackle predictive real-world quantum-chemistry applications with
strategies applicable to both NISQ and FTQC algorithms.

#### Acknowledgments

We thank the PEPR EPIQ - Quantum Software (ANR-22-PETQ-0007, J.-P.P.) for
funding D. Traore’s postdoctoral position. I.-M. Lygatsika’s PhD position has
been funded by the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program (grant No 810367, project EMC2,
J.-P. P. and Y.M.). Computations have been performed at IDRIS (Jean Zay) on
GENCI Grant no A0150712052 (J.-P.P.), at Scaleway (Jero) and on the SuperPOD
reference cluster EOS at NVIDIA. The authors thank D. Ruau and C. Hundt
(NVIDIA) for continuous support, as well as the NVIDIA Quantum team for
technical discussions.

## References

  * [1] A. Y. Kitaev, “Quantum measurements and the abelian stabilizer problem,” Electron. Colloquium Comput. Complex., vol. TR96, 1995.
  * [2] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, “Simulated quantum computation of molecular energies,” Science, vol. 309, no. 5741, pp. 1704–1707, 2005.
  * [3] M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information.  Cambridge University Press, 2010.
  * [4] J. D. Whitfield, J. Biamonte, and A. Aspuru-Guzik, “Simulation of electronic structure hamiltonians using quantum computers,” Molecular Physics, vol. 109, no. 5, pp. 735–750, 2011.
  * [5] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, “A variational eigenvalue solver on a photonic quantum processor,” Nature Communications, vol. 5, no. 1, pp. 1–7, 2014.
  * [6] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, et al., “The variational quantum eigensolver: a review of methods and best practices,” Physics Reports, vol. 986, pp. 1–128, 2022.
  * [7] E. Giner, B. Pradines, A. Ferté, R. Assaraf, A. Savin, and J. Toulouse, “Curing basis-set convergence of wave-function theory using density-functional theory: A systematically improvable approach,” Journal of Chemical Physics, vol. 149, p. 194301, 2018.
  * [8] P.-F. Loos, B. Pradines, A. Scemama, J. Toulouse, and E. Giner, “A density-based basis-set correction for wave function theory,” Journal of Physical Chemistry Letters, vol. 10, no. 11, pp. 2931–2937, 2019.
  * [9] E. Giner, A. Scemama, J. Toulouse, and P.-F. Loos, “Chemically accurate excitation energies with small basis sets,” Journal of Chemical Physics, vol. 151, no. 14, p. 144118, 2019.
  * [10] P.-F. Loos, B. Pradines, A. Scemama, E. Giner, and J. Toulouse, “Density-Based Basis-Set Incompleteness Correction for GW Methods,” Journal of Chemical Theory and Computation, vol. 16, pp. 1018–1028, Feb. 2020.
  * [11] E. Giner, A. Scemama, P.-F. Loos, and J. Toulouse, “A basis-set error correction based on density-functional theory for strongly correlated molecular systems,” Journal of Chemical Physics, vol. 152, p. 174104, May 2020.
  * [12] Y. Yao, E. Giner, J. Li, J. Toulouse, and C. J. Umrigar, “Almost exact energies for the gaussian-2 set with the semistochastic heat-bath configuration interaction method,” Journal of Chemical Physics, vol. 153, p. 124117, 2020.
  * [13] E. Giner, D. Traore, B. Pradines, and J. Toulouse, “Self-consistent density-based basis-set correction: How much do we lower total energies and improve dipole moments?,” Journal of Chemical Physics, vol. 155, p. 044109, Jul 2021.
  * [14] Y. Yao, E. Giner, T. A. Anderson, J. Toulouse, and C. J. Umrigar, “Accurate energies of transition metal atoms, ions, and monoxides using selected configuration interaction and density-based basis-set corrections,” Journal of Chemical Physics, vol. 155, p. 204104, 2021.
  * [15] D. Traore, E. Giner, and J. Toulouse, “Basis-set correction based on density-functional theory: Rigorous framework for a one-dimensional model,” Journal of Chemical Physics, vol. 156, no. 4, 2022.
  * [16] D. Traore, J. Toulouse, and E. Giner, “Basis-set correction for coupled-cluster estimation of dipole moments,” Journal of Chemical Physics, vol. 156, no. 17, 2022.
  * [17] D. Traore, E. Giner, and J. Toulouse, “Basis-set correction based on density-functional theory: Linear-response formalism for excited-state energies,” Journal of Chemical Physics, vol. 158, no. 23, 2023.
  * [18] D. Mester and M. Kállay, “Basis set limit of ccsd(t) energies: Explicit correlation versus density-based basis-set correction,” Journal of Chemical Theory and Computation, vol. 19, no. 22, pp. 8210–8222, 2023.
  * [19] A. Heßelmann, E. Giner, P. Reinhardt, P. J. Knowles, H.-J. Werner, and J. Toulouse, “A density-fitting implementation of the density-based basis-set correction method,” Journal of Computational Chemistry, vol. 45, no. 15, pp. 1247–1253, 2024.
  * [20] T. H. Dunning, “Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen,” Journal of Chemical Physics, vol. 90, p. 1007, 1989.
  * [21] T. Helgaker, W. Klopper, H. Koch, and J. Noga, “Basis-set convergence of correlated calculations on water,” Journal of Chemical Physics, vol. 106, no. 23, pp. 9639–9646, 1997.
  * [22] J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Physical Review Letters, vol. 77, p. 3865, 1996.
  * [23] I.-M. Lygatsika, Méthodes numériques pour les discrétisations Gaussiennes des problèmes en structure électronique.  Phd thesis, Sorbonne Université, Paris, France, October 2023.
  * [24] S. R. Jensen, S. Saha, J. A. Flores-Livas, W. Huhn, V. Blum, S. Goedecker, and L. Frediani, “The Elephant in the Room of Density Functional Theory Calculations,” Journal of Physical Chemistry Letters, vol. 8, no. 7, pp. 1449–1457, 2017.  PMID: 28291362.
  * [25] S. Lehtola, “Curing basis set overcompleteness with pivoted Cholesky decompositions,” Journal of Chemical Physics, vol. 151, no. 24, 2019.
  * [26] C. Lucas, “LAPACK-style codes for level 2 and 3 pivoted Cholesky factorizations,” LAPACK Working, vol. 4, no. 5, 2004.
  * [27] F. Aquilante, L. Boman, J. Boström, H. Koch, R. Lindh, A. S. de Merás, and T. B. Pedersen, “Cholesky decomposition techniques in electronic structure theory,” Linear-Scaling Techniques in Computational Chemistry and Physics: Methods and Applications, pp. 301–343, 2011.
  * [28] T. B. Pedersen, S. Lehtola, I. Fdez. Galván, and R. Lindh, “The versatility of the Cholesky decomposition in electronic structure theory,” Wiley Interdisciplinary Reviews: Computational Molecular Science, vol. 14, no. 1, p. e1692, 2024.
  * [29] O. Schutt and J. VandeVondele, “Machine learning adaptive basis sets for efficient large scale density functional theory simulation,” Journal of Chemical Theory and Computation, vol. 14, no. 8, pp. 4168–4175, 2018.
  * [30] Y. Mao, P. R. Horn, N. Mardirossian, T. Head-Gordon, C.-K. Skylaris, and M. Head-Gordon, “Approaching the basis set limit for DFT calculations using an environment-adapted minimal basis with perturbation theory: Formulation, proof of concept, and a pilot implementation,” Journal of Chemical Physics, vol. 145, no. 4, 2016.
  * [31] W. Wang and J. D. Whitfield, “Basis set generation and optimization in the nisq era with quiqbox.jl,” Journal of Chemical Theory and Computation, vol. 19, no. 22, pp. 8032–8052, 2023.  PMID: 37924295.
  * [32] H.-Y. Kwon, G. M. Curtin, Z. Morrow, C. Kelley, and E. Jakubikova, “Adaptive basis sets for practical quantum computing,” International Journal of Quantum Chemistry, vol. 123, no. 14, p. e27123, 2023.
  * [33] J. S. Kottmann, P. Schleich, T. Tamayo-Mendoza, and A. Aspuru-Guzik, “Reducing qubit requirements while maintaining numerical precision for the variational quantum eigensolver: A basis-set-free approach,” Journal of Physical Chemistry Letters, vol. 12, no. 1, pp. 663–673, 2021.  PMID: 33393305.
  * [34] 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, p. 3007, 2019.
  * [35] Y. S. Yordanov, V. Armaos, C. H. Barnes, and D. R. Arvidsson-Shukur, “Qubit-excitation-based adaptive variational quantum eigensolver,” Communications Physics, vol. 4, no. 1, pp. 1–11, 2021.
  * [36] E. Posenitskiy, V. G. Chilkuri, A. Ammar, M. Hapka, K. Pernal, R. Shinde, E. J. Landinez Borda, C. Filippi, K. Nakano, O. Kohulák, et al., “TREXIO: A file format and library for quantum chemistry,” Journal of Chemical Physics, vol. 158, no. 17, 2023.
  * [37] O. Adjoua and C. Feniou, “Sorbonne Université, CNRS and Qubit Pharmaceuticals,” 2023.
  * [38] B. Huron, J. Malrieu, and P. Rancurel, “Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wavefunctions,” Journal of Chemical Physics, vol. 58, no. 12, pp. 5745–5759, 1973.
  * [39] C. Feniou, M. Hassan, D. Traore, E. Giner, Y. Maday, and J.-P. Piquemal, “Overlap-ADAPT-VQE: Practical Quantum Chemistry on Quantum Computers via Overlap-Guided Compact Ansätze,” Communications Physics, vol. 6, 2023.
  * [40] C. Feniou, O. Adjoua, B. Claudon, J. Zylberman, E. Giner, and J.-P. Piquemal, “Sparse quantum state preparation for strongly correlated systems,” Journal of Physical Chemistry Letters, vol. 15, no. 11, pp. 3197–3205, 2024.
  * [41] I. Kassal and A. Aspuru-Guzik, “Quantum algorithm for molecular properties and geometry optimization,” Journal of Chemical Physics, vol. 131, no. 22, 2009.
  * [42] Y. Garniron, T. Applencourt, K. Gasperich, A. Benali, A. Ferté, J. Paquier, B. Pradines, R. Assaraf, P. Reinhardt, J. Toulouse, P. Barbaresco, N. Renon, G. David, J.-P. Malrieu, M. Véril, M. Caffarel, P.-F. Loos, E. Giner, and A. Scemama, “Quantum Package 2.0: An Open-Source Determinant-Driven Suite of Programs,” Journal of Chemical Theory and Computation, vol. 15, no. 6, pp. 3591–3609, 2019.  PMID: 31082265.
  * [43] F. Jensen, “Unifying general and segmented contracted basis sets. Segmented polarization consistent basis sets,” Journal of chemical theory and computation, vol. 10, no. 3, pp. 1074–1085, 2014.
  * [44] J. A. Pople, Modern theoretical chemistry.  Plenum, New York, 1976.
  * [45] E. R. Davidson and D. Feller, “Basis set selection for molecular calculations,” Chemical Reviews, vol. 86, no. 4, pp. 681–696, 1986.
  * [46] N. Tubman, C. Mejuto Zaera, J. Epstein, D. Hait, D. Levine, W. Huggins, Z. Jiang, J. McClean, R. Babbush, M. Head-Gordon, et al., “Postponing the orthogonality catastrophe: efficient state preparation for electronic structure simulations on quantum devices,” Bulletin of the American Physical Society, vol. 64, 2019.
  * [47] J. Almlöf, B. J. Deleeuw, P. R. Taylor, C. W. Bauschlicher Jr, and P. Siegbahn, “The dissociation energy of N2,” International Journal of Quantum Chemistry, vol. 36, no. S23, pp. 345–354, 1989.
  * [48] K. A. Peterson and T. H. Dunning Jr, “Intrinsic errors in several ab initio methods: the dissociation energy of N2,” Journal of Physical Chemistry, vol. 99, no. 12, pp. 3898–3901, 1995.
  * [49] J.-P. Malrieu, N. Guihéry, C. J. Calzado, and C. Angeli, “Bond electron pair: Its relevance and analysis from the quantum chemistry point of view,” Journal of Computional Chemistry, vol. 28, pp. 35–50, Jan. 2007.
  * [50] J. Robledo-Moreno, M. Motta, H. Haas, A. Javadi-Abhari, P. Jurcevic, W. Kirby, S. Martiel, K. Sharma, S. Sharma, T. Shirakawa, et al., “Chemistry beyond exact solutions on a quantum-centric supercomputer,” arXiv preprint arXiv:2405.05068, 2024.
  * [51] J. Provazza, K. Gunst, H. Zhai, G. K.-L. Chan, T. Shiozaki, N. C. Rubin, and A. F. White, “Fast emulation of fermionic circuits with matrix product states,” Journal of Chemical Theory and Computation, 2024.
  * [52] A. Halkier, W. Klopper, T. Helgaker, and P. Jørgensen, “Basis-set convergence of the molecular electric dipole moment,” Journal of Chemical Physics, vol. 111, no. 10, pp. 4424–4430, 1999.
  * [53] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” Nature, vol. 549, no. 7671, pp. 242–246, 2017.
  * [54] W. Sennane, J.-P. Piquemal, and M. J. Rančić, “Calculating the ground-state energy of benzene under spatial deformations with noisy quantum computing,” Physical Review A, vol. 107, no. 1, p. 012416, 2023.
  * [55] C. Feniou, B. Claudon, M. Hassan, A. Courtat, O. Adjoua, Y. Maday, and J.-P. Piquemal, “Greedy gradient-free adaptive variational quantum algorithms on a noisy intermediate scale quantum computer,” arXiv preprint arXiv:2306.17159, 2023.
  * [56] M. Haidar, M. J. Rančić, T. Ayral, Y. Maday, and J.-P. Piquemal, “Open source variational quantum eigensolver extension of the quantum learning machine for quantum chemistry,” WIREs Computational Molecular Science, vol. 13, no. 5, p. e1664, 2023.
  * [57] M. Haidar, M. J. Rančić, Y. Maday, and J.-P. Piquemal, “Extension of the trotterized unitary coupled cluster to triple excitations,” Journal of Physical Chemistry A, vol. 127, no. 15, pp. 3543–3550, 2023.
  * [58] H. Bayraktar, A. Charara, D. Clark, S. Cohen, T. Costa, Y.-L. L. Fang, Y. Gao, J. Guan, J. Gunnels, A. Haidar, et al., “cuquantum sdk: A high-performance library for accelerating quantum science,” in 2023 IEEE International Conference on Quantum Computing and Engineering (QCE), vol. 1, pp. 1050–1061, IEEE, 2023.
  * [59] J.-S. Kim, A. McCaskey, B. Heim, M. Modani, S. Stanwyck, and T. Costa, “Cuda quantum: The platform for integrated quantum-classical computing,” in 2023 60th ACM/IEEE Design Automation Conference (DAC), pp. 1–4, IEEE, 2023.